TreeFiller/src/FillerMCParticles.cc

// $Id: FillerMCParticles.cc,v 1.7 2009/02/26 17:04:03 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/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","source")),
  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.

  tws.AddBranch(mitName_.c_str(),&mcParticles_);

  if (genActive_)
    OS()->add(genMap_,genMapName_.c_str());
  if (simActive_)
    OS()->add(simMap_,simMapName_.c_str());
  if (trackingActive_)
    OS()->add(trackingMap_,trackingMapName_.c_str());
}

//--------------------------------------------------------------------------------------------------
void FillerMCParticles::FillDataBlock(const edm::Event      &event, 
                                      const edm::EventSetup &setup)
{
  // Loop over HepMC particle 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);
    }
    
  }
  
  
  //loop over aod GenParticle candidates and fill their information
  
  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());
        //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 OAK 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() ) {
        //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 HepMC particle 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);
          }
        }

      }

    }

  }

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