Reference Manual

Go to the documentation of this file.
00001 /* \class GenParticleProducer
00002  *
00003  * \author Luca Lista, INFN
00004  *
00005  * Convert HepMC GenEvent format into a collection of type
00006  * CandidateCollection containing objects of type GenParticle
00007  *
00008  * \version $Id: GenParticleProducer.cc,v 1.17 2010/10/20 23:23:45 wmtan Exp $
00009  *
00010  */
00011 #include "FWCore/Framework/interface/EDProducer.h"
00012 #include "FWCore/Utilities/interface/InputTag.h"
00013 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00014 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00015 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00016 
00017 #include <vector>
00018 #include <string>
00019 #include <map>
00020 #include <set>
00021 
00022 namespace edm { class ParameterSet; }
00023 namespace HepMC { class GenParticle; class GenEvent; }
00024 
00025 class GenParticleProducer : public edm::EDProducer {
00026  public:
00028   GenParticleProducer( const edm::ParameterSet & );
00030   ~GenParticleProducer();
00031 
00033   virtual void produce( edm::Event& e, const edm::EventSetup& );
00034   int chargeTimesThree( int ) const;
00035   bool convertParticle(reco::GenParticle& cand, const HepMC::GenParticle * part);
00036   bool fillDaughters(reco::GenParticleCollection& cand, const HepMC::GenParticle * part, size_t index);
00037   bool fillIndices(const HepMC::GenEvent * mc, std::vector<const HepMC::GenParticle*>& particles, std::vector<int>& barCodeVector, int offset);
00038   std::map<int, size_t> barcodes_;
00039   reco::GenParticleRefProd ref_;
00040 
00041  private:
00043   edm::InputTag src_;
00044   std::vector<std::string> vectorSrc_;
00045 
00047   bool firstEvent_; 
00049   bool abortOnUnknownPDGCode_;
00051   bool saveBarCodes_;
00053   std::vector<int> chargeP_, chargeM_;
00054   std::map<int, int> chargeMap_;
00055 
00057   bool doSubEvent_;
00058   bool useCF_;
00059 
00060 };
00061 
00062 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00063 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00064 
00065 //#include "SimDataFormats/HiGenData/interface/SubEventMap.h"
00066 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
00067 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
00068 
00069 #include "DataFormats/Common/interface/Handle.h"
00070 #include "FWCore/Framework/interface/ESHandle.h"
00071 #include "FWCore/Framework/interface/Event.h"
00072 #include "FWCore/Framework/interface/EventSetup.h"
00073 #include "FWCore/Utilities/interface/EDMException.h"
00074 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00075 #include <fstream>
00076 #include <algorithm>
00077 using namespace edm;
00078 using namespace reco;
00079 using namespace std;
00080 using namespace HepMC;
00081 
00082 static const int PDGCacheMax = 32768;
00083 static const double mmToCm = 0.1;
00084 
00085 GenParticleProducer::GenParticleProducer( const ParameterSet & cfg ) :
00086   firstEvent_(true), 
00087   abortOnUnknownPDGCode_( cfg.getUntrackedParameter<bool>( "abortOnUnknownPDGCode", true ) ),
00088   saveBarCodes_( cfg.getUntrackedParameter<bool>( "saveBarCodes", false ) ),
00089   chargeP_( PDGCacheMax, 0 ), chargeM_( PDGCacheMax, 0 ),
00090   doSubEvent_(cfg.getUntrackedParameter<bool>( "doSubEvent", false )),
00091   useCF_(cfg.getUntrackedParameter<bool>( "useCrossingFrame", false ))
00092 {
00093   produces<GenParticleCollection>();
00094   if( saveBarCodes_ ) {
00095     std::string alias( cfg.getParameter<std::string>( "@module_label" ) );
00096     produces<vector<int> >().setBranchAlias( alias + "BarCodes" );
00097   }                               
00098 
00099   if(doSubEvent_){
00100      vectorSrc_ = cfg.getParameter<std::vector<std::string> >( "srcVector" );
00101      //     produces<SubEventMap>();
00102   }else if(useCF_) src_ = cfg.getUntrackedParameter<InputTag>( "src" , InputTag("mix","generator"));
00103   else src_ = cfg.getParameter<InputTag>( "src" );
00104   
00105 }
00106 
00107 GenParticleProducer::~GenParticleProducer() { 
00108 }
00109 
00110 int GenParticleProducer::chargeTimesThree( int id ) const {
00111   if( std::abs( id ) < PDGCacheMax ) 
00112     return id > 0 ? chargeP_[ id ] : chargeM_[ - id ];
00113   map<int, int>::const_iterator f = chargeMap_.find( id );
00114   if ( f == chargeMap_.end() )  {
00115     if ( abortOnUnknownPDGCode_ )
00116       throw edm::Exception( edm::errors::LogicError ) 
00117         << "invalid PDG id: " << id << endl;
00118     else
00119       return HepPDT::ParticleID(id).threeCharge();
00120   }
00121   return f->second;
00122 }
00123 
00124 void GenParticleProducer::produce( Event& evt, const EventSetup& es ) {
00125 
00126   if (firstEvent_) {
00127      ESHandle<HepPDT::ParticleDataTable> pdt;
00128      es.getData( pdt );
00129      for( HepPDT::ParticleDataTable::const_iterator p = pdt->begin(); p != pdt->end(); ++ p ) {
00130        const HepPDT::ParticleID & id = p->first;
00131        int pdgId = id.pid(), apdgId = std::abs( pdgId );
00132        int q3 = id.threeCharge();
00133        if ( apdgId < PDGCacheMax && pdgId > 0 ) {
00134          chargeP_[ apdgId ] = q3;
00135          chargeM_[ apdgId ] = -q3;
00136        } else if ( apdgId < PDGCacheMax ) {
00137          chargeP_[ apdgId ] = -q3;
00138          chargeM_[ apdgId ] = q3;
00139        } else {
00140          chargeMap_[ pdgId ] = q3;
00141          chargeMap_[ -pdgId ] = -q3;
00142        }
00143      }
00144      firstEvent_ = false; 
00145    }
00146       
00147    barcodes_.clear();
00148    
00149    size_t totalSize = 0;
00150    const GenEvent * mc = 0;   
00151    std::vector<Handle<HepMCProduct> > heps;
00152    MixCollection<HepMCProduct>* cfhepmcprod = 0;
00153    size_t npiles = vectorSrc_.size();
00154 
00155    if(useCF_){
00156       Handle<CrossingFrame<HepMCProduct> > cf;
00157       evt.getByLabel(InputTag("mix","generator"),cf);
00158       cfhepmcprod = new MixCollection<HepMCProduct>(cf.product());
00159       npiles = cfhepmcprod->size();
00160       for(unsigned int icf = 0; icf < npiles; ++icf){
00161          totalSize += cfhepmcprod->getObject(icf).GetEvent()->particles_size();
00162       }
00163    }else if (doSubEvent_){
00164       for(size_t i = 0; i < npiles; ++i){
00165         //       cout<<"Tag "<<vectorSrc_[i]<<endl;
00166          Handle<HepMCProduct> handle;
00167          heps.push_back(handle);
00168          evt.getByLabel( vectorSrc_[i], heps[i] );
00169          totalSize += heps[i]->GetEvent()->particles_size();
00170       }
00171    }else{
00172       Handle<HepMCProduct> mcp;
00173       evt.getByLabel( src_, mcp );
00174       mc = mcp->GetEvent();
00175       if( mc == 0 ) 
00176          throw edm::Exception( edm::errors::InvalidReference ) 
00177             << "HepMC has null pointer to GenEvent" << endl;
00178       totalSize  = mc->particles_size();
00179    }
00180       
00181    // initialise containers 
00182    const size_t size = totalSize;
00183   vector<const HepMC::GenParticle *> particles( size );
00184   auto_ptr<GenParticleCollection> candsPtr( new GenParticleCollection( size ) );
00185   //  auto_ptr<SubEventMap> subsPtr( new SubEventMap() );
00186   auto_ptr<vector<int> > barCodeVector( new vector<int>( size ) );
00187   ref_ = evt.getRefBeforePut<GenParticleCollection>();
00188   GenParticleCollection & cands = * candsPtr;
00189   //  SubEventMap & subs = *subsPtr;
00190   size_t offset = 0;
00191   size_t suboffset = 0;
00192 
00194   if(doSubEvent_ || useCF_){
00195      for(size_t i = 0; i < npiles; ++i){
00196         barcodes_.clear();
00197         if(useCF_) mc = cfhepmcprod->getObject(i).GetEvent();
00198         else mc = heps[i]->GetEvent();
00199 
00200         //Look whether heavy ion/signal event
00201         bool isHI = false;
00202         const HepMC::HeavyIon * hi = mc->heavy_ion();
00203         if(hi && hi->Ncoll_hard() > 1) isHI = true;
00204         size_t num_particles = mc->particles_size();
00205         fillIndices(mc, particles, *barCodeVector, offset);
00206         // fill output collection and save association 
00207         for( size_t i = offset; i < offset + num_particles; ++ i ) {
00208 
00209            const HepMC::GenParticle * part = particles[ i ];
00210            reco::GenParticle & cand = cands[ i ];
00211            // convert HepMC::GenParticle to new reco::GenParticle
00212            convertParticle(cand, part);
00213            cand.resetDaughters( ref_.id() );
00214         }
00215 
00216         for( size_t d = offset; d < offset + num_particles; ++ d ) {
00217            const HepMC::GenParticle * part = particles[ d ];
00218            const GenVertex * productionVertex = part->production_vertex();
00219            int sub_id = 0;
00220            if ( productionVertex != 0 ) {
00221               sub_id = productionVertex->id();
00222               if(!isHI) sub_id = 0;
00223               // search barcode map and attach daughters 
00224               fillDaughters(cands,part,d);
00225            }else{
00226               const GenVertex * endVertex = part->end_vertex();
00227               if(endVertex != 0) sub_id = endVertex->id();
00228               else throw cms::Exception( "SubEventID" )<<"SubEvent not determined. Particle has no production and no end vertex!"<<endl;
00229            }
00230            if(sub_id < 0) sub_id = 0;
00231            int new_id = sub_id + suboffset;
00232            GenParticleRef dref( ref_, d );
00233            //      subs.insert(dref,new_id);   // For SubEventMap
00234            cands[d].setCollisionId(new_id); // For new GenParticle
00235            LogDebug("VertexId")<<"SubEvent offset 3 : "<<suboffset;
00236         }
00237         int nsub = -2;
00238         if(isHI){
00239            nsub = hi->Ncoll_hard()+1;
00240            suboffset += nsub;
00241         }else{
00242            suboffset += 1;
00243         }
00244         offset += num_particles;
00245      }
00246   }else{
00247      fillIndices(mc, particles, *barCodeVector, 0);
00248      
00249      // fill output collection and save association
00250      for( size_t i = 0; i < particles.size(); ++ i ) {
00251         const HepMC::GenParticle * part = particles[ i ];
00252         reco::GenParticle & cand = cands[ i ];
00253         // convert HepMC::GenParticle to new reco::GenParticle
00254         convertParticle(cand, part);
00255         cand.resetDaughters( ref_.id() );
00256      }
00257      
00258      // fill references to daughters
00259      for( size_t d = 0; d < cands.size(); ++ d ) {
00260         const HepMC::GenParticle * part = particles[ d ];
00261         const GenVertex * productionVertex = part->production_vertex();
00262         // search barcode map and attach daughters
00263         if ( productionVertex != 0 ) fillDaughters(cands,part,d);
00264         cands[d].setCollisionId(0);
00265      }
00266   }
00267   
00268   evt.put( candsPtr );
00269   if(saveBarCodes_) evt.put( barCodeVector );
00270   //  if(doSubEvent_) evt.put(subsPtr); // For SubEventMap
00271   if(cfhepmcprod) delete cfhepmcprod;
00272 
00273 }
00274 
00275 bool GenParticleProducer::convertParticle(reco::GenParticle& cand, const HepMC::GenParticle * part){
00276    Candidate::LorentzVector p4( part->momentum() );
00277    int pdgId = part->pdg_id();
00278    cand.setThreeCharge( chargeTimesThree( pdgId ) );
00279    cand.setPdgId( pdgId );
00280    cand.setStatus( part->status() );
00281    cand.setP4( p4 );
00282    cand.setCollisionId(0);
00283    const GenVertex * v = part->production_vertex();
00284    if ( v != 0 ) {
00285       ThreeVector vtx = v->point3d();
00286       Candidate::Point vertex( vtx.x() * mmToCm, vtx.y() * mmToCm, vtx.z() * mmToCm );
00287       cand.setVertex( vertex );
00288    } else {
00289       cand.setVertex( Candidate::Point( 0, 0, 0 ) );
00290    }
00291    return true;
00292 }
00293 
00294 bool GenParticleProducer::fillDaughters(reco::GenParticleCollection& cands, const HepMC::GenParticle * part, size_t index){
00295 
00296    const GenVertex * productionVertex = part->production_vertex();
00297    size_t numberOfMothers = productionVertex->particles_in_size();
00298    if ( numberOfMothers > 0 ) {
00299       GenVertex::particles_in_const_iterator motherIt = productionVertex->particles_in_const_begin();
00300       for( ; motherIt != productionVertex->particles_in_const_end(); motherIt++) {
00301          const HepMC::GenParticle * mother = * motherIt;
00302          size_t m = barcodes_.find( mother->barcode() )->second;
00303          cands[ m ].addDaughter( GenParticleRef( ref_, index ) );
00304          cands[ index ].addMother( GenParticleRef( ref_, m ) );
00305       }
00306    }
00307 
00308    return true;
00309 }
00310 
00311 bool GenParticleProducer::fillIndices(const HepMC::GenEvent * mc, vector<const HepMC::GenParticle*>& particles, vector<int>& barCodeVector, int offset){
00312    size_t idx = offset;
00313    HepMC::GenEvent::particle_const_iterator begin = mc->particles_begin(), end = mc->particles_end();
00314    for( HepMC::GenEvent::particle_const_iterator p = begin; p != end; ++ p ) {
00315       const HepMC::GenParticle * particle = * p;
00316       size_t barCode_this_event = particle->barcode();
00317       size_t barCode = barCode_this_event + offset;
00318       if( barcodes_.find(barCode) != barcodes_.end() )
00319          throw cms::Exception( "WrongReference" )
00320             << "barcodes are duplicated! " << endl;
00321       particles[idx] = particle;
00322       barCodeVector[idx] = barCode;
00323       barcodes_.insert( make_pair(barCode_this_event, idx ++) );
00324    }
00325       return true;
00326 }
00327 
00328 
00329 #include "FWCore/Framework/interface/MakerMacros.h"
00330 
00331 DEFINE_FWK_MODULE( GenParticleProducer );
00332
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/PhysicsTools/HepMCCandAlgos/plugins/GenParticleProducer.cc
