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