TreeFiller/src/FillerMCParticles.cc

// $Id: FillerMCParticles.cc,v 1.12 2009/03/19 16:14:29 loizides 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/HepMCProduct/interface/HepMCProduct.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.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 "MitAna/DataTree/interface/Names.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)),
  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)),
  mcParticles_(new mithep::MCParticleArr(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 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 (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 (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() ) {
        if (verbose_>1) {
          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());
            }
          }
          TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
          const SimTrack &theSimTrack = iM->g4Tracks().at(iM->g4Tracks().size()-1);
          printf("trackId = %i\n",theSimTrack.trackId());
          mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrack.trackId());
          trackingMap_->Add(theRef, outSimParticle);
        }
      }
    }
  }

  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.13
Committed:	Sun Mar 22 11:35:45 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009, Mit_008
Changes since 1.12:	+6 -2 lines
Log Message:	Fix branchsize to meaningful minimum.
#	User	Rev	Content
1	loizides	1.13	// $Id: FillerMCParticles.cc,v 1.12 2009/03/19 16:14:29 loizides Exp $
2	bendavid	1.1
3			#include "MitProd/TreeFiller/interface/FillerMCParticles.h"
4			#include "FWCore/MessageLogger/interface/MessageLogger.h"
5			#include "DataFormats/Common/interface/Handle.h"
6	bendavid	1.8	#include "DataFormats/Common/interface/RefToPtr.h"
7	bendavid	1.1	#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
8	bendavid	1.8	#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
9			#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
10	bendavid	1.4	#include "SimDataFormats/Track/interface/SimTrack.h"
11			#include "SimDataFormats/Track/interface/SimTrackContainer.h"
12			#include "SimDataFormats/Vertex/interface/SimVertex.h"
13			#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
14	bendavid	1.1	#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
15			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
16			#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
17			#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
18			#include "MitAna/DataTree/interface/Names.h"
19
20			using namespace std;
21			using namespace edm;
22			using namespace mithep;
23
24			//--------------------------------------------------------------------------------------------------
25			FillerMCParticles::FillerMCParticles(const ParameterSet &cfg, const char *name, bool active) :
26			BaseFiller(cfg, name, active),
27			genActive_(Conf().getUntrackedParameter<bool>("genActive",true)),
28	bendavid	1.8	useAodGen_(Conf().getUntrackedParameter<bool>("useAodGen",true)),
29	bendavid	1.1	simActive_(Conf().getUntrackedParameter<bool>("simActive",true)),
30	bendavid	1.4	trackingActive_(Conf().getUntrackedParameter<bool>("trackingActive",false)),
31	loizides	1.12	genEdmName_(Conf().getUntrackedParameter<string>("genEdmName","genParticles")),
32	bendavid	1.4	simEdmName_(Conf().getUntrackedParameter<string>("simEdmName","g4SimHits")),
33	loizides	1.9	trackingEdmName_(Conf().getUntrackedParameter<string>("trackingEdmName",
34			"mergedtruth:MergedTrackTruth")),
35	loizides	1.3	genMapName_(Conf().getUntrackedParameter<string>("genMapName","GenMap")),
36			simMapName_(Conf().getUntrackedParameter<string>("simMapName","SimMap")),
37	bendavid	1.4	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
38	bendavid	1.1	mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkMCPartBrn)),
39			mcParticles_(new mithep::MCParticleArr(250)),
40	bendavid	1.8	genMap_(genActive_?new mithep::GenParticleBarcodeMap:0),
41			aodGenMap_((genActive_&&useAodGen_)?new mithep::AODGenParticleMap:0),
42			simMap_(simActive_?new mithep::SimTrackTidMap:0),
43	bendavid	1.4	trackingMap_(trackingActive_?new mithep::TrackingParticleMap:0)
44	bendavid	1.1	{
45			// Constructor.
46			}
47
48			//--------------------------------------------------------------------------------------------------
49			FillerMCParticles::~FillerMCParticles()
50			{
51			// Destructor.
52
53			delete mcParticles_;
54			delete genMap_;
55	bendavid	1.8	delete aodGenMap_;
56	bendavid	1.1	delete simMap_;
57	bendavid	1.4	delete trackingMap_;
58	bendavid	1.1	}
59
60			//--------------------------------------------------------------------------------------------------
61			void FillerMCParticles::BookDataBlock(TreeWriter &tws)
62			{
63	loizides	1.3	// Add branch to tree and publish our maps.
64	bendavid	1.1
65	loizides	1.13	Int_t brsize = tws.GetDefaultBrSize();
66			if (brsize<128*1024)
67			brsize=128*1024;
68
69			tws.AddBranch(mitName_,&mcParticles_,brsize);
70	loizides	1.9	OS()->add<mithep::MCParticleArr>(mcParticles_,mitName_);
71	loizides	1.2
72	loizides	1.9	if (genActive_ && !genMapName_.empty()) {
73			genMap_->SetBrName(mitName_);
74			if (genMap_)
75			OS()->add(genMap_,genMapName_);
76			}
77			if (simActive_ && !simMapName_.empty()) {
78			simMap_->SetBrName(mitName_);
79			if (simMap_)
80			OS()->add(simMap_,simMapName_);
81			}
82			if (trackingActive_ && !trackingMapName_.empty()) {
83			trackingMap_->SetBrName(mitName_);
84			if (trackingMap_)
85			OS()->add(trackingMap_,trackingMapName_);
86			}
87	bendavid	1.1	}
88
89			//--------------------------------------------------------------------------------------------------
90			void FillerMCParticles::FillDataBlock(const edm::Event &event,
91	loizides	1.2	const edm::EventSetup &setup)
92	bendavid	1.1	{
93	loizides	1.9	// Loop over generated or simulated particles and fill their information.
94	bendavid	1.1
95	bendavid	1.7	mcParticles_->Delete();
96	bendavid	1.1
97	bendavid	1.8	if (genActive_)
98	loizides	1.3	genMap_->Reset();
99	bendavid	1.1
100	bendavid	1.8	if (genActive_ && !useAodGen_) {
101	bendavid	1.1	Handle<edm::HepMCProduct> hHepMCProduct;
102			GetProduct(genEdmName_, hHepMCProduct, event);
103
104			const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
105
106	loizides	1.9	// loop over all hepmc particles and copy their information
107	bendavid	1.1	for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
108			pgen != GenEvent.particles_end(); ++pgen) {
109
110			HepMC::GenParticle mcPart = (pgen);
111			if(!mcPart) continue;
112
113			mithep::MCParticle *genParticle = mcParticles_->Allocate();
114			new (genParticle) mithep::MCParticle(mcPart->momentum().x(),mcPart->momentum().y(),
115	loizides	1.11	mcPart->momentum().z(),mcPart->momentum().e(),
116			mcPart->pdg_id(),mcPart->status());
117	bendavid	1.1	genParticle->SetIsGenerated();
118			genMap_->Add(mcPart->barcode(), genParticle);
119			}
120			}
121
122	bendavid	1.8	if (genActive_ && useAodGen_) {
123	loizides	1.9
124	bendavid	1.8	aodGenMap_->Reset();
125
126			Handle<reco::GenParticleCollection> hGenPProduct;
127			GetProduct(genEdmName_, hGenPProduct, event);
128
129	loizides	1.9	// element by element aligned collection of hepmc barcodes associated to
130			// the genparticles
131	bendavid	1.8	Handle<std::vector<int> > genBarcodes;
132			if (simActive_)
133			GetProduct(genEdmName_, genBarcodes, event);
134
135			const reco::GenParticleCollection genParticles = *(hGenPProduct.product());
136
137	loizides	1.11	// loop over all genparticles and copy their information
138	bendavid	1.8	for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
139			pgen != genParticles.end(); ++pgen) {
140
141			mithep::MCParticle *mcPart = mcParticles_->Allocate();
142			new (mcPart) mithep::MCParticle(pgen->px(),pgen->py(),
143			pgen->pz(),pgen->energy(),
144			pgen->pdgId(),
145			pgen->status());
146
147			mcPart->SetIsGenerated();
148
149	loizides	1.9	// add hepmc barcode association, needed to merge in sim particles
150	bendavid	1.8	if (simActive_) {
151			int genBarcode = genBarcodes->at(pgen - genParticles.begin());
152			genMap_->Add(genBarcode, mcPart);
153			}
154
155			edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
156			aodGenMap_->Add(thePtr, mcPart);
157			}
158			}
159
160	bendavid	1.1	if (simActive_) {
161	loizides	1.3	simMap_->Reset();
162	bendavid	1.4
163			Handle<edm::SimTrackContainer> hSimTrackProduct;
164			GetProduct(simEdmName_, hSimTrackProduct, event);
165
166			const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
167
168	bendavid	1.1	// loop through all simParticles
169	bendavid	1.4	for (SimTrackContainer::const_iterator iM = simTracks.begin();
170			iM != simTracks.end(); ++iM) {
171
172	bendavid	1.8	mithep::MCParticle *outSimParticle = 0;
173	bendavid	1.1
174	bendavid	1.4	if (genActive_ && (iM->genpartIndex() >= 0) ) {
175			outSimParticle = genMap_->GetMit(iM->genpartIndex());
176	loizides	1.9	if (verbose_>1) {
177			printf("Sim px,py,pz,e = %f,%f,%f,%f\nGen px,py,pz,e = %f,%f,%f,%f\n",
178			iM->momentum().px(),iM->momentum().py(),
179			iM->momentum().pz(),iM->momentum().e(),
180			outSimParticle->Px(),outSimParticle->Py(),
181			outSimParticle->Pz(),outSimParticle->E());
182			}
183			} else {
184	bendavid	1.1	outSimParticle = mcParticles_->Allocate();
185	bendavid	1.4	new (outSimParticle) mithep::MCParticle(iM->momentum().px(),
186			iM->momentum().py(),
187			iM->momentum().pz(),
188			iM->momentum().e(),
189	loizides	1.9	iM->type(), -99);
190	bendavid	1.1	}
191
192			outSimParticle->SetIsSimulated();
193	bendavid	1.8	simMap_->Add(iM->trackId(), outSimParticle);
194	bendavid	1.4	}
195			}
196
197	loizides	1.10	// loop through TrackingParticles and build association map to MCParticles
198	bendavid	1.4	if (trackingActive_) {
199			trackingMap_->Reset();
200
201			Handle<TrackingParticleCollection> hTrackingParticleProduct;
202			GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
203
204	bendavid	1.5	const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
205	bendavid	1.4
206			// loop through all TrackingParticles
207	bendavid	1.5	for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
208	bendavid	1.4	iM != trackingParticles.end(); ++iM) {
209
210			if ( simActive_ && iM->g4Tracks().size() ) {
211	loizides	1.9	if (verbose_>1) {
212			printf("Tracking particle has %i SimTracks\n",iM->g4Tracks().size());
213			if (iM->g4Tracks().size()>1) {
214			for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
215			iST != iM->g4Tracks().end(); ++iST) {
216			printf("g4 trackid = %i\n",iST->trackId());
217			}
218			}
219			TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
220			const SimTrack &theSimTrack = iM->g4Tracks().at(iM->g4Tracks().size()-1);
221			printf("trackId = %i\n",theSimTrack.trackId());
222			mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrack.trackId());
223			trackingMap_->Add(theRef, outSimParticle);
224			}
225	bendavid	1.4	}
226	bendavid	1.1	}
227			}
228
229			mcParticles_->Trim();
230			}
231
232			//--------------------------------------------------------------------------------------------------
233			void FillerMCParticles::ResolveLinks(const edm::Event &event,
234	loizides	1.2	const edm::EventSetup &setup)
235	bendavid	1.1	{
236	loizides	1.9	// Loop over generated and simulated particles and resolve their links.
237	bendavid	1.1
238	bendavid	1.8	if (genActive_ && !useAodGen_) {
239	bendavid	1.1
240			Handle<edm::HepMCProduct> hHepMCProduct;
241			GetProduct(genEdmName_, hHepMCProduct, event);
242
243			const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
244
245			for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
246			pgen != GenEvent.particles_end(); ++pgen) {
247
248			HepMC::GenParticle mcPart = (pgen);
249			if(!mcPart) continue;
250
251	loizides	1.11	// check if genpart has a decay vertex
252	bendavid	1.1	HepMC::GenVertex *dVertex = mcPart->end_vertex();
253			if (!dVertex) continue;
254
255	loizides	1.11	// find corresponding mithep genparticle parent in association table
256	bendavid	1.1	mithep::MCParticle *genParent = genMap_->GetMit(mcPart->barcode());
257
258	loizides	1.11	// set decay vertex
259			// division by 10.0 is needed due to HepMC use of mm instead of cm for distance units
260	bendavid	1.1	genParent->SetVertex(dVertex->point3d().x()/10.0,
261			dVertex->point3d().y()/10.0,
262			dVertex->point3d().z()/10.0);
263
264	loizides	1.11	// loop through daugthers
265	loizides	1.3	for (HepMC::GenVertex::particles_out_const_iterator pgenD =
266			dVertex->particles_out_const_begin();
267			pgenD != dVertex->particles_out_const_end(); ++pgenD) {
268	bendavid	1.1	HepMC::GenParticle mcDaughter = (pgenD);
269			mithep::MCParticle *genDaughter = genMap_->GetMit(mcDaughter->barcode());
270			genParent->AddDaughter(genDaughter);
271			if (!(genDaughter->HasMother()))
272			genDaughter->SetMother(genParent);
273			}
274			}
275			}
276
277	loizides	1.11	// loop over aod GenParticle candidates and resolve their links
278	bendavid	1.8	if (genActive_ && useAodGen_) {
279
280			Handle<reco::GenParticleCollection> hGenPProduct;
281			GetProduct(genEdmName_, hGenPProduct, event);
282
283			const reco::GenParticleCollection genParticles = *(hGenPProduct.product());
284
285	loizides	1.11	// loop over all genparticles and copy their information
286	bendavid	1.8	for (reco::GenParticleCollection::const_iterator pgen = genParticles.begin();
287			pgen != genParticles.end(); ++pgen) {
288
289			int nDaughters = pgen->numberOfDaughters();
290
291			if (nDaughters>0) {
292
293			edm::Ptr<reco::GenParticle> thePtr(hGenPProduct, pgen - genParticles.begin());
294			MCParticle *mcMother = aodGenMap_->GetMit(thePtr);
295			for (int i=0; i<nDaughters; ++i) {
296			const reco::Candidate *genDaughter = pgen->daughter(i);
297			MCParticle *mcDaughter = aodGenMap_->GetMit(refToPtr(pgen->daughterRef(i)));
298	loizides	1.11	// set mother decay vertex
299	bendavid	1.8	if (i==0)
300			mcMother->SetVertex(genDaughter->vx(),genDaughter->vy(),genDaughter->vz());
301
302	loizides	1.11	// set mother-daughter links
303	bendavid	1.8	mcMother->AddDaughter(mcDaughter);
304			if (!mcDaughter->HasMother())
305			mcDaughter->SetMother(mcMother);
306
307			}
308			}
309			}
310
311			}
312
313	loizides	1.11	// loop over SimTracks and resolve links
314	bendavid	1.1	if (simActive_) {
315	bendavid	1.4	Handle<edm::SimTrackContainer> hSimTrackProduct;
316			GetProduct(simEdmName_, hSimTrackProduct, event);
317
318			const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
319
320			Handle<std::vector<SimVertex> > hSimVertexProduct;
321			GetProduct(simEdmName_, hSimVertexProduct, event);
322
323			const edm::SimVertexContainer simVertexes = *(hSimVertexProduct.product());
324
325	bendavid	1.1	// loop through all simParticles
326	bendavid	1.4	for (SimTrackContainer::const_iterator iM = simTracks.begin();
327			iM != simTracks.end(); ++iM) {
328
329			if (iM->vertIndex()>=0) {
330	bendavid	1.8	mithep::MCParticle *simDaughter = simMap_->GetMit(iM->trackId());
331	bendavid	1.4	const SimVertex &theVertex = simVertexes.at(iM->vertIndex());
332			if (theVertex.parentIndex()>=0) {
333	bendavid	1.8	mithep::MCParticle *simParent = simMap_->GetMit(theVertex.parentIndex());
334	loizides	1.9	simParent->SetVertex(theVertex.position().x(),
335			theVertex.position().y(),
336			theVertex.position().z());
337	bendavid	1.6	//make sure we don't double count the decay tree
338			if ( !simParent->HasDaughter(simDaughter) ) {
339			simParent->AddDaughter(simDaughter);
340			}
341			if ( !simDaughter->HasMother() ) {
342			simDaughter->SetMother(simParent);
343			}
344	bendavid	1.1	}
345	loizides	1.9	}
346			}
347			}
348	bendavid	1.4
349	loizides	1.9	if (verbose_>1 && trackingActive_) {
350			Handle<TrackingParticleCollection> hTrackingParticleProduct;
351			GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
352
353			const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
354
355			// loop through all simParticles
356			for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
357			iM != trackingParticles.end(); ++iM) {
358
359			TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
360
361			mithep::MCParticle *simParent = trackingMap_->GetMit(theRef);
362			printf("MCParticle pdg = %i, Daughter Pdgs: ", simParent->PdgId());
363			for (UInt_t i=0; i<simParent->NDaughters(); ++i) {
364			printf("%i, ", simParent->Daughter(i)->PdgId());
365	bendavid	1.1	}
366	loizides	1.9	printf("\n");
367
368			printf("TrackingParticle pdg = %i, Daughter Pdgs: ", iM->pdgId());
369			for (TrackingVertexRefVector::iterator v= iM->decayVertices().begin();
370			v != iM->decayVertices().end(); ++v) {
371			for (TrackingParticleRefVector::iterator pd = v->get()->daughterTracks().begin();
372			pd != v->get()->daughterTracks().end(); ++pd) {
373			printf("%i, ", pd->get()->pdgId());
374			}
375			}
376			printf("\n");
377	bendavid	1.1	}
378			}
379			}