TreeFiller/src/FillerMCParticles.cc

// $Id: FillerMCParticles.cc,v 1.20 2010/01/07 11:04:32 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerMCParticles.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/RefToPtr.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/METReco/interface/METCollection.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataTree/interface/MCParticleCol.h"
#include "MitAna/DataTree/interface/TrackingParticleCol.h"
#include "MitProd/ObjectService/interface/ObjectService.h"

using namespace std;
using namespace edm;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
FillerMCParticles::FillerMCParticles(const ParameterSet &cfg, const char *name, bool active) : 
  BaseFiller(cfg, name, active),
  genActive_(Conf().getUntrackedParameter<bool>("genActive",true)),
  useAodGen_(Conf().getUntrackedParameter<bool>("useAodGen",true)),
  simActive_(Conf().getUntrackedParameter<bool>("simActive",true)),
  trackingActive_(Conf().getUntrackedParameter<bool>("trackingActive",false)),
  fillTracking_(Conf().getUntrackedParameter<bool>("fillTracking",false)),
  genEdmName_(Conf().getUntrackedParameter<string>("genEdmName","genParticles")),
  simEdmName_(Conf().getUntrackedParameter<string>("simEdmName","g4SimHits")),
  trackingEdmName_(Conf().getUntrackedParameter<string>("trackingEdmName",
                                                        "mergedtruth:MergedTrackTruth")),
  genMapName_(Conf().getUntrackedParameter<string>("genMapName","GenMap")),
  simMapName_(Conf().getUntrackedParameter<string>("simMapName","SimMap")),
  trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
  mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkMCPartBrn)),
  mitTrackingName_(Conf().getUntrackedParameter<string>("mitTrackingName",Names::gkTrackingParticleBrn)),
  mcParticles_(new mithep::MCParticleArr(250)),
  trackingParticles_(new mithep::TrackingParticleArr(250)),
  genMap_(genActive_?new mithep::GenParticleBarcodeMap:0),
  aodGenMap_((genActive_&&useAodGen_)?new mithep::AODGenParticleMap:0),
  simMap_(simActive_?new mithep::SimTrackTidMap:0),
  trackingMap_(trackingActive_?new mithep::TrackingParticleMap:0)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
FillerMCParticles::~FillerMCParticles()
{
  // Destructor.

  delete mcParticles_;
  delete trackingParticles_;
  delete genMap_;
  delete aodGenMap_;
  delete simMap_;
  delete trackingMap_;
}

//--------------------------------------------------------------------------------------------------
void FillerMCParticles::BookDataBlock(TreeWriter &tws)
{
  // Add branch to tree and publish our maps.

  Int_t brsize = tws.GetDefaultBrSize();
  if (brsize<128*1024)
    brsize=128*1024;

  tws.AddBranch(mitName_,&mcParticles_,brsize);
  OS()->add<mithep::MCParticleArr>(mcParticles_,mitName_);

  if (trackingActive_ && fillTracking_) {
    tws.AddBranch(mitTrackingName_,&trackingParticles_,brsize);
    OS()->add<mithep::TrackingParticleArr>(trackingParticles_,mitTrackingName_);
  }
  
  if (genActive_ && !genMapName_.empty()) {
    genMap_->SetBrName(mitName_);
    if (genMap_)
      OS()->add(genMap_,genMapName_);
  }
  if (simActive_ && !simMapName_.empty()) {
    simMap_->SetBrName(mitName_);
    if (simMap_)
      OS()->add(simMap_,simMapName_);
  }
  if (trackingActive_ && !trackingMapName_.empty()) {
    trackingMap_->SetBrName(mitName_);
    if (trackingMap_)
      OS()->add(trackingMap_,trackingMapName_);
  }
}

//--------------------------------------------------------------------------------------------------
void FillerMCParticles::FillDataBlock(const edm::Event      &event, 
                                      const edm::EventSetup &setup)
{
  // Loop over generated or simulated particles and fill their information.

  mcParticles_->Delete();
  
  if (trackingActive_ && fillTracking_) {
    trackingParticles_->Delete();
  }

  if (genActive_)
    genMap_->Reset();
  
  if (genActive_ && !useAodGen_) {
    Handle<edm::HepMCProduct> hHepMCProduct;
    GetProduct(genEdmName_, hHepMCProduct, event);  
  
    const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();  

    // loop over all hepmc particles and copy their information
    for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
        pgen != GenEvent.particles_end(); ++pgen) {
  
      HepMC::GenParticle *mcPart = (*pgen);
      if(!mcPart) 
        continue;
  
      mithep::MCParticle *genParticle = mcParticles_->Allocate();
      new (genParticle) mithep::MCParticle(mcPart->momentum().x(),mcPart->momentum().y(),
                                           mcPart->momentum().z(),mcPart->momentum().e(),
                                           mcPart->pdg_id(),mcPart->status());
      genParticle->SetIsGenerated();                                          
      genMap_->Add(mcPart->barcode(), genParticle);
    }
  }
  
  if (genActive_ && useAodGen_) {
 
    aodGenMap_->Reset();
  
    Handle<reco::GenParticleCollection> hGenPProduct;
    GetProduct(genEdmName_, hGenPProduct, event);  
    
    // element by element aligned collection of hepmc barcodes associated to
    // the genparticles
    Handle<std::vector<int> > genBarcodes;
    if (simActive_)
      GetProduct(genEdmName_, genBarcodes, event);  
  
    const reco::GenParticleCollection genParticles = *(hGenPProduct.product());  
  
    // loop over all genparticles and copy their information
    for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
        pgen != genParticles.end(); ++pgen) {
        
      mithep::MCParticle *mcPart = mcParticles_->Allocate();
      new (mcPart) mithep::MCParticle(pgen->px(),pgen->py(),
                                      pgen->pz(),pgen->energy(),
                                      pgen->pdgId(),pgen->status());
                                            
      mcPart->SetIsGenerated();                                          
      
      // add hepmc barcode association, needed to merge in sim particles
      if (simActive_) {
        int genBarcode = genBarcodes->at(pgen - genParticles.begin());
        genMap_->Add(genBarcode, mcPart);
      }
      
      edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
      aodGenMap_->Add(thePtr, mcPart);
    }
  }
  
  if (simActive_) {
    simMap_->Reset();

    Handle<edm::SimTrackContainer> hSimTrackProduct;
    GetProduct(simEdmName_, hSimTrackProduct, event);

    const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());

    // loop through all simParticles
    for (SimTrackContainer::const_iterator iM = simTracks.begin(); 
        iM != simTracks.end(); ++iM) {

      mithep::MCParticle *outSimParticle = 0;
      
      if (genActive_ && (iM->genpartIndex() >= 0) ) {
        outSimParticle = genMap_->GetMit(iM->genpartIndex());
        if (verbose_>1) {
          printf("Sim px,py,pz,e = %f,%f,%f,%f\nGen px,py,pz,e = %f,%f,%f,%f\n",       
                 iM->momentum().px(),iM->momentum().py(),
                 iM->momentum().pz(),iM->momentum().e(),
                 outSimParticle->Px(),outSimParticle->Py(),
                 outSimParticle->Pz(),outSimParticle->E());
        }
      } else {
        outSimParticle = mcParticles_->Allocate();
        new (outSimParticle) mithep::MCParticle(iM->momentum().px(),
                                                iM->momentum().py(),
                                                iM->momentum().pz(),
                                                iM->momentum().e(),
                                                iM->type(), -99);
      }
      
      outSimParticle->SetIsSimulated();
      simMap_->Add(iM->trackId(), outSimParticle);
    }
  }

  // loop through TrackingParticles and build association map to MCParticles
  if (trackingActive_) {
    trackingMap_->Reset();
  
    Handle<TrackingParticleCollection> hTrackingParticleProduct;
    GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
    
    const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
    
    // loop through all TrackingParticles
    for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
        iM != trackingParticles.end(); ++iM) {

      if ( simActive_ && iM->g4Tracks().size() ) {
        TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
        const SimTrack &theSimTrack = iM->g4Tracks().at(iM->g4Tracks().size()-1);
        mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrack.trackId());
        trackingMap_->Add(theRef, outSimParticle);
        
        if (fillTracking_) {
          mithep::TrackingParticle *outTracking = trackingParticles_->AddNew();
          
          //fill associated sim tracks
          for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
                 iST != iM->g4Tracks().end(); ++iST) {
            outTracking->AddMCPart(simMap_->GetMit(theSimTrack.trackId()));
          }
          
          //fill hit information
          reco::HitPattern hitPattern;
          int nHits = 0;
          
          //loop through sim hits and construct intermediate reco::HitPattern
          for (std::vector<PSimHit>::const_iterator iSimHit = iM->pSimHit_begin();
                iSimHit != iM->pSimHit_end(); ++iSimHit) {

            unsigned int detectorIdIndex = iSimHit->detUnitId();
            DetId detectorId = DetId(detectorIdIndex);
            
            InvalidTrackingRecHit tHit(detectorId,TrackingRecHit::valid);
          
            hitPattern.set(tHit,nHits);
            ++nHits;
          
          }

          //construct hit layer bitmask as used in bambu
          BitMask48 hits;
          // search for all good crossed layers (with or without hit)
          for (Int_t hi=0; hi<hitPattern.numberOfHits(); hi++) {
            uint32_t hit = hitPattern.getHitPattern(hi);
            if (hitPattern.getHitType(hit)<=1) {
              if (hitPattern.trackerHitFilter(hit)) {
                hits.SetBit(hitReader_.Layer(hit));
              }
            }
          }
          
          outTracking->SetHits(hits); 
          
          if (verbose_>1) {
            printf("Hit Layer Mask: ");
            for (Int_t biti=47; biti >= 0; --biti) {
                  printf("%i",hits.TestBit(biti));
            }
            printf("\n");
          }
          
        }
        
        if (verbose_>1) {
          printf("trackId = %i\n",theSimTrack.trackId());
          printf("Tracking particle has %i SimTracks\n",iM->g4Tracks().size());
          if (iM->g4Tracks().size()>1) {
            for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
                 iST != iM->g4Tracks().end(); ++iST) {
              printf("g4 trackid = %i\n",iST->trackId());
            }
          }
        }
      }
    }
  }

  mcParticles_->Trim();
}

//--------------------------------------------------------------------------------------------------
void FillerMCParticles::ResolveLinks(const edm::Event      &event, 
                                     const edm::EventSetup &setup)
{
  // Loop over generated and simulated particles and resolve their links.

  if (genActive_ && !useAodGen_) {
  
    Handle<edm::HepMCProduct> hHepMCProduct;
    GetProduct(genEdmName_, hHepMCProduct, event);
    
    const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();  
  
    for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
        pgen != GenEvent.particles_end(); ++pgen) {
  
      HepMC::GenParticle *mcPart = (*pgen);
      if(!mcPart) 
        continue;
      
      // check if genpart has a decay vertex
      HepMC::GenVertex *dVertex = mcPart->end_vertex();
      if (!dVertex) 
        continue;

      // find corresponding mithep genparticle parent in association table
      mithep::MCParticle *genParent = genMap_->GetMit(mcPart->barcode());
  
      // set decay vertex
      // division by 10.0 is needed due to HepMC use of mm instead of cm for distance units
      genParent->SetVertex(dVertex->point3d().x()/10.0,
                           dVertex->point3d().y()/10.0,
                           dVertex->point3d().z()/10.0);
  
      // loop through daugthers
      for (HepMC::GenVertex::particles_out_const_iterator pgenD = 
             dVertex->particles_out_const_begin(); 
           pgenD != dVertex->particles_out_const_end(); ++pgenD) {
        HepMC::GenParticle *mcDaughter = (*pgenD);
        mithep::MCParticle *genDaughter = genMap_->GetMit(mcDaughter->barcode());
        genParent->AddDaughter(genDaughter);
        if (!(genDaughter->HasMother()))
          genDaughter->SetMother(genParent);
      }
    }
  }
  
  // loop over aod GenParticle candidates and resolve their links
  if (genActive_ && useAodGen_) {
  
    Handle<reco::GenParticleCollection> hGenPProduct;
    GetProduct(genEdmName_, hGenPProduct, event);  
  
    const reco::GenParticleCollection genParticles = *(hGenPProduct.product());  
    // loop over all genparticles and copy their information
    for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
        pgen != genParticles.end(); ++pgen) {
  
      int nDaughters = pgen->numberOfDaughters();
      if (nDaughters>0) {
        edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
        MCParticle *mcMother = aodGenMap_->GetMit(thePtr);
        for (int i=0; i<nDaughters; ++i) {
          const reco::Candidate *genDaughter = pgen->daughter(i);
          MCParticle *mcDaughter = aodGenMap_->GetMit(refToPtr(pgen->daughterRef(i)));
          // set mother decay vertex
          if (i==0)
            mcMother->SetVertex(genDaughter->vx(),genDaughter->vy(),genDaughter->vz());

          // set mother-daughter links
          mcMother->AddDaughter(mcDaughter);
          if (!mcDaughter->HasMother())
            mcDaughter->SetMother(mcMother);
        }
      }
    }
  }
  
  // loop over SimTracks and resolve links
  if (simActive_) {
    Handle<edm::SimTrackContainer> hSimTrackProduct;
    GetProduct(simEdmName_, hSimTrackProduct, event);

    const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
   
    Handle<std::vector<SimVertex> > hSimVertexProduct;
    GetProduct(simEdmName_, hSimVertexProduct, event);

    const edm::SimVertexContainer simVertexes = *(hSimVertexProduct.product());

    // loop through all simParticles
    for (SimTrackContainer::const_iterator iM = simTracks.begin(); 
        iM != simTracks.end(); ++iM) {

      if (iM->vertIndex()>=0) {
        mithep::MCParticle *simDaughter = simMap_->GetMit(iM->trackId());
        const SimVertex &theVertex = simVertexes.at(iM->vertIndex());
        if (theVertex.parentIndex()>=0) {
          mithep::MCParticle *simParent = simMap_->GetMit(theVertex.parentIndex());
          simParent->SetVertex(theVertex.position().x(),
                               theVertex.position().y(),
                               theVertex.position().z());
          //make sure we don't double count the decay tree
          if ( !simParent->HasDaughter(simDaughter) ) {
            simParent->AddDaughter(simDaughter);
          }
          if ( !simDaughter->HasMother() ) {
            simDaughter->SetMother(simParent);
          }
        }
      }
    }
  }

  if (verbose_>1 && trackingActive_) {
    Handle<TrackingParticleCollection> hTrackingParticleProduct;
    GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
    
    const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());  
   
    // loop through all simParticles
    for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin(); 
         iM != trackingParticles.end(); ++iM) {
      
      TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
      
      mithep::MCParticle *simParent = trackingMap_->GetMit(theRef);
      printf("MCParticle pdg = %i, Daughter Pdgs: ", simParent->PdgId());
      for (UInt_t i=0; i<simParent->NDaughters(); ++i) {
        printf("%i, ", simParent->Daughter(i)->PdgId());
      }
      printf("\n");
      
      printf("TrackingParticle pdg = %i, Daughter Pdgs: ", iM->pdgId());
      for (TrackingVertexRefVector::iterator v= iM->decayVertices().begin(); 
           v != iM->decayVertices().end(); ++v) {
        for (TrackingParticleRefVector::iterator pd = v->get()->daughterTracks().begin(); 
             pd != v->get()->daughterTracks().end(); ++pd) {
          printf("%i, ", pd->get()->pdgId());
        }
      }
      printf("\n");
    }
  }
}
Revision:	1.21
Committed:	Thu Mar 18 20:21:00 2010 UTC (15 years, 1 month ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013
Changes since 1.20:	+2 -2 lines
Log Message:	Fix beginrun,beginjob mess
#	Content
1	// $Id: FillerMCParticles.cc,v 1.20 2010/01/07 11:04:32 bendavid Exp $
2
3	#include "MitProd/TreeFiller/interface/FillerMCParticles.h"
4	#include "FWCore/MessageLogger/interface/MessageLogger.h"
5	#include "DataFormats/Common/interface/Handle.h"
6	#include "DataFormats/Common/interface/RefToPtr.h"
7	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
8	#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
9	#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
10	#include "DataFormats/METReco/interface/METCollection.h"
11	#include "SimDataFormats/Track/interface/SimTrack.h"
12	#include "SimDataFormats/Track/interface/SimTrackContainer.h"
13	#include "SimDataFormats/Vertex/interface/SimVertex.h"
14	#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
15	#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
16	#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
17	#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
18	#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
19	#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
20	#include "DataFormats/TrackReco/interface/HitPattern.h"
21	#include "MitAna/DataTree/interface/Names.h"
22	#include "MitAna/DataTree/interface/MCParticleCol.h"
23	#include "MitAna/DataTree/interface/TrackingParticleCol.h"
24	#include "MitProd/ObjectService/interface/ObjectService.h"
25
26	using namespace std;
27	using namespace edm;
28	using namespace mithep;
29
30	//--------------------------------------------------------------------------------------------------
31	FillerMCParticles::FillerMCParticles(const ParameterSet &cfg, const char *name, bool active) :
32	BaseFiller(cfg, name, active),
33	genActive_(Conf().getUntrackedParameter<bool>("genActive",true)),
34	useAodGen_(Conf().getUntrackedParameter<bool>("useAodGen",true)),
35	simActive_(Conf().getUntrackedParameter<bool>("simActive",true)),
36	trackingActive_(Conf().getUntrackedParameter<bool>("trackingActive",false)),
37	fillTracking_(Conf().getUntrackedParameter<bool>("fillTracking",false)),
38	genEdmName_(Conf().getUntrackedParameter<string>("genEdmName","genParticles")),
39	simEdmName_(Conf().getUntrackedParameter<string>("simEdmName","g4SimHits")),
40	trackingEdmName_(Conf().getUntrackedParameter<string>("trackingEdmName",
41	"mergedtruth:MergedTrackTruth")),
42	genMapName_(Conf().getUntrackedParameter<string>("genMapName","GenMap")),
43	simMapName_(Conf().getUntrackedParameter<string>("simMapName","SimMap")),
44	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
45	mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkMCPartBrn)),
46	mitTrackingName_(Conf().getUntrackedParameter<string>("mitTrackingName",Names::gkTrackingParticleBrn)),
47	mcParticles_(new mithep::MCParticleArr(250)),
48	trackingParticles_(new mithep::TrackingParticleArr(250)),
49	genMap_(genActive_?new mithep::GenParticleBarcodeMap:0),
50	aodGenMap_((genActive_&&useAodGen_)?new mithep::AODGenParticleMap:0),
51	simMap_(simActive_?new mithep::SimTrackTidMap:0),
52	trackingMap_(trackingActive_?new mithep::TrackingParticleMap:0)
53	{
54	// Constructor.
55	}
56
57	//--------------------------------------------------------------------------------------------------
58	FillerMCParticles::~FillerMCParticles()
59	{
60	// Destructor.
61
62	delete mcParticles_;
63	delete trackingParticles_;
64	delete genMap_;
65	delete aodGenMap_;
66	delete simMap_;
67	delete trackingMap_;
68	}
69
70	//--------------------------------------------------------------------------------------------------
71	void FillerMCParticles::BookDataBlock(TreeWriter &tws)
72	{
73	// Add branch to tree and publish our maps.
74
75	Int_t brsize = tws.GetDefaultBrSize();
76	if (brsize<128*1024)
77	brsize=128*1024;
78
79	tws.AddBranch(mitName_,&mcParticles_,brsize);
80	OS()->add<mithep::MCParticleArr>(mcParticles_,mitName_);
81
82	if (trackingActive_ && fillTracking_) {
83	tws.AddBranch(mitTrackingName_,&trackingParticles_,brsize);
84	OS()->add<mithep::TrackingParticleArr>(trackingParticles_,mitTrackingName_);
85	}
86
87	if (genActive_ && !genMapName_.empty()) {
88	genMap_->SetBrName(mitName_);
89	if (genMap_)
90	OS()->add(genMap_,genMapName_);
91	}
92	if (simActive_ && !simMapName_.empty()) {
93	simMap_->SetBrName(mitName_);
94	if (simMap_)
95	OS()->add(simMap_,simMapName_);
96	}
97	if (trackingActive_ && !trackingMapName_.empty()) {
98	trackingMap_->SetBrName(mitName_);
99	if (trackingMap_)
100	OS()->add(trackingMap_,trackingMapName_);
101	}
102	}
103
104	//--------------------------------------------------------------------------------------------------
105	void FillerMCParticles::FillDataBlock(const edm::Event &event,
106	const edm::EventSetup &setup)
107	{
108	// Loop over generated or simulated particles and fill their information.
109
110	mcParticles_->Delete();
111
112	if (trackingActive_ && fillTracking_) {
113	trackingParticles_->Delete();
114	}
115
116	if (genActive_)
117	genMap_->Reset();
118
119	if (genActive_ && !useAodGen_) {
120	Handle<edm::HepMCProduct> hHepMCProduct;
121	GetProduct(genEdmName_, hHepMCProduct, event);
122
123	const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
124
125	// loop over all hepmc particles and copy their information
126	for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
127	pgen != GenEvent.particles_end(); ++pgen) {
128
129	HepMC::GenParticle mcPart = (pgen);
130	if(!mcPart)
131	continue;
132
133	mithep::MCParticle *genParticle = mcParticles_->Allocate();
134	new (genParticle) mithep::MCParticle(mcPart->momentum().x(),mcPart->momentum().y(),
135	mcPart->momentum().z(),mcPart->momentum().e(),
136	mcPart->pdg_id(),mcPart->status());
137	genParticle->SetIsGenerated();
138	genMap_->Add(mcPart->barcode(), genParticle);
139	}
140	}
141
142	if (genActive_ && useAodGen_) {
143
144	aodGenMap_->Reset();
145
146	Handle<reco::GenParticleCollection> hGenPProduct;
147	GetProduct(genEdmName_, hGenPProduct, event);
148
149	// element by element aligned collection of hepmc barcodes associated to
150	// the genparticles
151	Handle<std::vector<int> > genBarcodes;
152	if (simActive_)
153	GetProduct(genEdmName_, genBarcodes, event);
154
155	const reco::GenParticleCollection genParticles = *(hGenPProduct.product());
156
157	// loop over all genparticles and copy their information
158	for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
159	pgen != genParticles.end(); ++pgen) {
160
161	mithep::MCParticle *mcPart = mcParticles_->Allocate();
162	new (mcPart) mithep::MCParticle(pgen->px(),pgen->py(),
163	pgen->pz(),pgen->energy(),
164	pgen->pdgId(),pgen->status());
165
166	mcPart->SetIsGenerated();
167
168	// add hepmc barcode association, needed to merge in sim particles
169	if (simActive_) {
170	int genBarcode = genBarcodes->at(pgen - genParticles.begin());
171	genMap_->Add(genBarcode, mcPart);
172	}
173
174	edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
175	aodGenMap_->Add(thePtr, mcPart);
176	}
177	}
178
179	if (simActive_) {
180	simMap_->Reset();
181
182	Handle<edm::SimTrackContainer> hSimTrackProduct;
183	GetProduct(simEdmName_, hSimTrackProduct, event);
184
185	const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
186
187	// loop through all simParticles
188	for (SimTrackContainer::const_iterator iM = simTracks.begin();
189	iM != simTracks.end(); ++iM) {
190
191	mithep::MCParticle *outSimParticle = 0;
192
193	if (genActive_ && (iM->genpartIndex() >= 0) ) {
194	outSimParticle = genMap_->GetMit(iM->genpartIndex());
195	if (verbose_>1) {
196	printf("Sim px,py,pz,e = %f,%f,%f,%f\nGen px,py,pz,e = %f,%f,%f,%f\n",
197	iM->momentum().px(),iM->momentum().py(),
198	iM->momentum().pz(),iM->momentum().e(),
199	outSimParticle->Px(),outSimParticle->Py(),
200	outSimParticle->Pz(),outSimParticle->E());
201	}
202	} else {
203	outSimParticle = mcParticles_->Allocate();
204	new (outSimParticle) mithep::MCParticle(iM->momentum().px(),
205	iM->momentum().py(),
206	iM->momentum().pz(),
207	iM->momentum().e(),
208	iM->type(), -99);
209	}
210
211	outSimParticle->SetIsSimulated();
212	simMap_->Add(iM->trackId(), outSimParticle);
213	}
214	}
215
216	// loop through TrackingParticles and build association map to MCParticles
217	if (trackingActive_) {
218	trackingMap_->Reset();
219
220	Handle<TrackingParticleCollection> hTrackingParticleProduct;
221	GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
222
223	const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
224
225	// loop through all TrackingParticles
226	for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
227	iM != trackingParticles.end(); ++iM) {
228
229	if ( simActive_ && iM->g4Tracks().size() ) {
230	TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
231	const SimTrack &theSimTrack = iM->g4Tracks().at(iM->g4Tracks().size()-1);
232	mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrack.trackId());
233	trackingMap_->Add(theRef, outSimParticle);
234
235	if (fillTracking_) {
236	mithep::TrackingParticle *outTracking = trackingParticles_->AddNew();
237
238	//fill associated sim tracks
239	for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
240	iST != iM->g4Tracks().end(); ++iST) {
241	outTracking->AddMCPart(simMap_->GetMit(theSimTrack.trackId()));
242	}
243
244	//fill hit information
245	reco::HitPattern hitPattern;
246	int nHits = 0;
247
248	//loop through sim hits and construct intermediate reco::HitPattern
249	for (std::vector<PSimHit>::const_iterator iSimHit = iM->pSimHit_begin();
250	iSimHit != iM->pSimHit_end(); ++iSimHit) {
251
252	unsigned int detectorIdIndex = iSimHit->detUnitId();
253	DetId detectorId = DetId(detectorIdIndex);
254
255	InvalidTrackingRecHit tHit(detectorId,TrackingRecHit::valid);
256
257	hitPattern.set(tHit,nHits);
258	++nHits;
259
260	}
261
262	//construct hit layer bitmask as used in bambu
263	BitMask48 hits;
264	// search for all good crossed layers (with or without hit)
265	for (Int_t hi=0; hi<hitPattern.numberOfHits(); hi++) {
266	uint32_t hit = hitPattern.getHitPattern(hi);
267	if (hitPattern.getHitType(hit)<=1) {
268	if (hitPattern.trackerHitFilter(hit)) {
269	hits.SetBit(hitReader_.Layer(hit));
270	}
271	}
272	}
273
274	outTracking->SetHits(hits);
275
276	if (verbose_>1) {
277	printf("Hit Layer Mask: ");
278	for (Int_t biti=47; biti >= 0; --biti) {
279	printf("%i",hits.TestBit(biti));
280	}
281	printf("\n");
282	}
283
284	}
285
286	if (verbose_>1) {
287	printf("trackId = %i\n",theSimTrack.trackId());
288	printf("Tracking particle has %i SimTracks\n",iM->g4Tracks().size());
289	if (iM->g4Tracks().size()>1) {
290	for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
291	iST != iM->g4Tracks().end(); ++iST) {
292	printf("g4 trackid = %i\n",iST->trackId());
293	}
294	}
295	}
296	}
297	}
298	}
299
300	mcParticles_->Trim();
301	}
302
303	//--------------------------------------------------------------------------------------------------
304	void FillerMCParticles::ResolveLinks(const edm::Event &event,
305	const edm::EventSetup &setup)
306	{
307	// Loop over generated and simulated particles and resolve their links.
308
309	if (genActive_ && !useAodGen_) {
310
311	Handle<edm::HepMCProduct> hHepMCProduct;
312	GetProduct(genEdmName_, hHepMCProduct, event);
313
314	const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
315
316	for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
317	pgen != GenEvent.particles_end(); ++pgen) {
318
319	HepMC::GenParticle mcPart = (pgen);
320	if(!mcPart)
321	continue;
322
323	// check if genpart has a decay vertex
324	HepMC::GenVertex *dVertex = mcPart->end_vertex();
325	if (!dVertex)
326	continue;
327
328	// find corresponding mithep genparticle parent in association table
329	mithep::MCParticle *genParent = genMap_->GetMit(mcPart->barcode());
330
331	// set decay vertex
332	// division by 10.0 is needed due to HepMC use of mm instead of cm for distance units
333	genParent->SetVertex(dVertex->point3d().x()/10.0,
334	dVertex->point3d().y()/10.0,
335	dVertex->point3d().z()/10.0);
336
337	// loop through daugthers
338	for (HepMC::GenVertex::particles_out_const_iterator pgenD =
339	dVertex->particles_out_const_begin();
340	pgenD != dVertex->particles_out_const_end(); ++pgenD) {
341	HepMC::GenParticle mcDaughter = (pgenD);
342	mithep::MCParticle *genDaughter = genMap_->GetMit(mcDaughter->barcode());
343	genParent->AddDaughter(genDaughter);
344	if (!(genDaughter->HasMother()))
345	genDaughter->SetMother(genParent);
346	}
347	}
348	}
349
350	// loop over aod GenParticle candidates and resolve their links
351	if (genActive_ && useAodGen_) {
352
353	Handle<reco::GenParticleCollection> hGenPProduct;
354	GetProduct(genEdmName_, hGenPProduct, event);
355
356	const reco::GenParticleCollection genParticles = *(hGenPProduct.product());
357	// loop over all genparticles and copy their information
358	for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
359	pgen != genParticles.end(); ++pgen) {
360
361	int nDaughters = pgen->numberOfDaughters();
362	if (nDaughters>0) {
363	edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
364	MCParticle *mcMother = aodGenMap_->GetMit(thePtr);
365	for (int i=0; i<nDaughters; ++i) {
366	const reco::Candidate *genDaughter = pgen->daughter(i);
367	MCParticle *mcDaughter = aodGenMap_->GetMit(refToPtr(pgen->daughterRef(i)));
368	// set mother decay vertex
369	if (i==0)
370	mcMother->SetVertex(genDaughter->vx(),genDaughter->vy(),genDaughter->vz());
371
372	// set mother-daughter links
373	mcMother->AddDaughter(mcDaughter);
374	if (!mcDaughter->HasMother())
375	mcDaughter->SetMother(mcMother);
376	}
377	}
378	}
379	}
380
381	// loop over SimTracks and resolve links
382	if (simActive_) {
383	Handle<edm::SimTrackContainer> hSimTrackProduct;
384	GetProduct(simEdmName_, hSimTrackProduct, event);
385
386	const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
387
388	Handle<std::vector<SimVertex> > hSimVertexProduct;
389	GetProduct(simEdmName_, hSimVertexProduct, event);
390
391	const edm::SimVertexContainer simVertexes = *(hSimVertexProduct.product());
392
393	// loop through all simParticles
394	for (SimTrackContainer::const_iterator iM = simTracks.begin();
395	iM != simTracks.end(); ++iM) {
396
397	if (iM->vertIndex()>=0) {
398	mithep::MCParticle *simDaughter = simMap_->GetMit(iM->trackId());
399	const SimVertex &theVertex = simVertexes.at(iM->vertIndex());
400	if (theVertex.parentIndex()>=0) {
401	mithep::MCParticle *simParent = simMap_->GetMit(theVertex.parentIndex());
402	simParent->SetVertex(theVertex.position().x(),
403	theVertex.position().y(),
404	theVertex.position().z());
405	//make sure we don't double count the decay tree
406	if ( !simParent->HasDaughter(simDaughter) ) {
407	simParent->AddDaughter(simDaughter);
408	}
409	if ( !simDaughter->HasMother() ) {
410	simDaughter->SetMother(simParent);
411	}
412	}
413	}
414	}
415	}
416
417	if (verbose_>1 && trackingActive_) {
418	Handle<TrackingParticleCollection> hTrackingParticleProduct;
419	GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
420
421	const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
422
423	// loop through all simParticles
424	for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
425	iM != trackingParticles.end(); ++iM) {
426
427	TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
428
429	mithep::MCParticle *simParent = trackingMap_->GetMit(theRef);
430	printf("MCParticle pdg = %i, Daughter Pdgs: ", simParent->PdgId());
431	for (UInt_t i=0; i<simParent->NDaughters(); ++i) {
432	printf("%i, ", simParent->Daughter(i)->PdgId());
433	}
434	printf("\n");
435
436	printf("TrackingParticle pdg = %i, Daughter Pdgs: ", iM->pdgId());
437	for (TrackingVertexRefVector::iterator v= iM->decayVertices().begin();
438	v != iM->decayVertices().end(); ++v) {
439	for (TrackingParticleRefVector::iterator pd = v->get()->daughterTracks().begin();
440	pd != v->get()->daughterTracks().end(); ++pd) {
441	printf("%i, ", pd->get()->pdgId());
442	}
443	}
444	printf("\n");
445	}
446	}
447	}