TreeFiller/src/FillerMCParticles.cc

// $Id: FillerMCParticles.cc,v 1.6 2008/10/06 13:35:56 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerMCParticles.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.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)),
  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::GenParticleMap:0),
  simMap_(simActive_?new mithep::SimTrackMap:0),
  trackingMap_(trackingActive_?new mithep::TrackingParticleMap:0)
{
  // Constructor.
}

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

  delete mcParticles_;
  delete genMap_;
  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();
  
    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 (simActive_) {
    simTidMap_.clear();
    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) {

      SimTrackRef theRef(hSimTrackProduct, iM-simTracks.begin());
      simTidMap_[iM->trackId()] = theRef; 
      mithep::MCParticle *outSimParticle;
      
      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(theRef, 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());
        SimTrackRef theSimTrackRef = simTidMap_[theSimTrack.trackId()];
        mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrackRef);
        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_) {
  
    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 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) {
        SimTrackRef theRef(hSimTrackProduct, iM-simTracks.begin());
        mithep::MCParticle *simDaughter = simMap_->GetMit(theRef);
        const SimVertex &theVertex = simVertexes.at(iM->vertIndex());
        if (theVertex.parentIndex()>=0) {
          mithep::MCParticle *simParent = simMap_->GetMit(simTidMap_[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.7
Committed:	Thu Feb 26 17:04:03 2009 UTC (16 years, 2 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_008pre1
Changes since 1.6:	+2 -2 lines
Log Message:	Switch from Reset to Delete calls on arrays, since we now use some heap
#	User	Rev	Content
1	bendavid	1.7	// $Id: FillerMCParticles.cc,v 1.6 2008/10/06 13:35:56 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			#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
7	bendavid	1.4	#include "SimDataFormats/Track/interface/SimTrack.h"
8			#include "SimDataFormats/Track/interface/SimTrackContainer.h"
9			#include "SimDataFormats/Vertex/interface/SimVertex.h"
10			#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
11	bendavid	1.1	#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
12			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
13			#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
14			#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
15			#include "MitAna/DataTree/interface/Names.h"
16
17			using namespace std;
18			using namespace edm;
19			using namespace mithep;
20
21			//--------------------------------------------------------------------------------------------------
22			FillerMCParticles::FillerMCParticles(const ParameterSet &cfg, const char *name, bool active) :
23			BaseFiller(cfg, name, active),
24			genActive_(Conf().getUntrackedParameter<bool>("genActive",true)),
25			simActive_(Conf().getUntrackedParameter<bool>("simActive",true)),
26	bendavid	1.4	trackingActive_(Conf().getUntrackedParameter<bool>("trackingActive",false)),
27	bendavid	1.1	genEdmName_(Conf().getUntrackedParameter<string>("genEdmName","source")),
28	bendavid	1.4	simEdmName_(Conf().getUntrackedParameter<string>("simEdmName","g4SimHits")),
29			trackingEdmName_(Conf().getUntrackedParameter<string>("trackingEdmName","mergedtruth:MergedTrackTruth")),
30	loizides	1.3	genMapName_(Conf().getUntrackedParameter<string>("genMapName","GenMap")),
31			simMapName_(Conf().getUntrackedParameter<string>("simMapName","SimMap")),
32	bendavid	1.4	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
33	bendavid	1.1	mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkMCPartBrn)),
34			mcParticles_(new mithep::MCParticleArr(250)),
35	loizides	1.3	genMap_(genActive_?new mithep::GenParticleMap:0),
36	bendavid	1.4	simMap_(simActive_?new mithep::SimTrackMap:0),
37			trackingMap_(trackingActive_?new mithep::TrackingParticleMap:0)
38	bendavid	1.1	{
39			// Constructor.
40			}
41
42			//--------------------------------------------------------------------------------------------------
43			FillerMCParticles::~FillerMCParticles()
44			{
45			// Destructor.
46
47			delete mcParticles_;
48			delete genMap_;
49			delete simMap_;
50	bendavid	1.4	delete trackingMap_;
51	bendavid	1.1	}
52
53			//--------------------------------------------------------------------------------------------------
54			void FillerMCParticles::BookDataBlock(TreeWriter &tws)
55			{
56	loizides	1.3	// Add branch to tree and publish our maps.
57	bendavid	1.1
58			tws.AddBranch(mitName_.c_str(),&mcParticles_);
59	loizides	1.2
60	loizides	1.3	if (genActive_)
61			OS()->add(genMap_,genMapName_.c_str());
62			if (simActive_)
63			OS()->add(simMap_,simMapName_.c_str());
64	bendavid	1.4	if (trackingActive_)
65			OS()->add(trackingMap_,trackingMapName_.c_str());
66	bendavid	1.1	}
67
68			//--------------------------------------------------------------------------------------------------
69			void FillerMCParticles::FillDataBlock(const edm::Event &event,
70	loizides	1.2	const edm::EventSetup &setup)
71	bendavid	1.1	{
72			// Loop over HepMC particle and fill their information.
73
74	bendavid	1.7	mcParticles_->Delete();
75	bendavid	1.1
76			if (genActive_) {
77	loizides	1.3	genMap_->Reset();
78	bendavid	1.1
79			Handle<edm::HepMCProduct> hHepMCProduct;
80			GetProduct(genEdmName_, hHepMCProduct, event);
81
82			const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
83
84			//loop over all hepmc particles and copy their information
85			for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
86			pgen != GenEvent.particles_end(); ++pgen) {
87
88			HepMC::GenParticle mcPart = (pgen);
89			if(!mcPart) continue;
90
91			mithep::MCParticle *genParticle = mcParticles_->Allocate();
92			new (genParticle) mithep::MCParticle(mcPart->momentum().x(),mcPart->momentum().y(),
93			mcPart->momentum().z(),mcPart->momentum().e(),
94			mcPart->pdg_id(),
95			mcPart->status());
96
97			genParticle->SetIsGenerated();
98
99			genMap_->Add(mcPart->barcode(), genParticle);
100			}
101
102			}
103
104			if (simActive_) {
105	bendavid	1.4	simTidMap_.clear();
106	loizides	1.3	simMap_->Reset();
107	bendavid	1.4
108			Handle<edm::SimTrackContainer> hSimTrackProduct;
109			GetProduct(simEdmName_, hSimTrackProduct, event);
110
111			const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
112
113	bendavid	1.1	// loop through all simParticles
114	bendavid	1.4	for (SimTrackContainer::const_iterator iM = simTracks.begin();
115			iM != simTracks.end(); ++iM) {
116
117			SimTrackRef theRef(hSimTrackProduct, iM-simTracks.begin());
118			simTidMap_[iM->trackId()] = theRef;
119	bendavid	1.1	mithep::MCParticle *outSimParticle;
120
121	bendavid	1.4	if (genActive_ && (iM->genpartIndex() >= 0) ) {
122			outSimParticle = genMap_->GetMit(iM->genpartIndex());
123			//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());
124	bendavid	1.1	}
125			else {
126			outSimParticle = mcParticles_->Allocate();
127	bendavid	1.4	new (outSimParticle) mithep::MCParticle(iM->momentum().px(),
128			iM->momentum().py(),
129			iM->momentum().pz(),
130			iM->momentum().e(),
131			iM->type(),
132			-99);
133	bendavid	1.1	}
134
135			outSimParticle->SetIsSimulated();
136			simMap_->Add(theRef, outSimParticle);
137	bendavid	1.4
138			}
139			}
140
141			//loop through TrackingParticles and build association map to OAK MCParticles
142			if (trackingActive_) {
143			trackingMap_->Reset();
144
145			Handle<TrackingParticleCollection> hTrackingParticleProduct;
146			GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
147
148	bendavid	1.5	const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
149	bendavid	1.4
150			// loop through all TrackingParticles
151	bendavid	1.5	for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
152	bendavid	1.4	iM != trackingParticles.end(); ++iM) {
153
154			if ( simActive_ && iM->g4Tracks().size() ) {
155			//printf("Tracking particle has %i SimTracks\n",iM->g4Tracks().size());
156			// if (iM->g4Tracks().size()>1) {
157			// for (std::vector<SimTrack>::const_iterator iST = iM->g4Tracks().begin();
158			// iST != iM->g4Tracks().end(); ++iST) {
159			// printf("g4 trackid = %i\n",iST->trackId());
160			// }
161			// }
162			TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
163			const SimTrack &theSimTrack = iM->g4Tracks().at(iM->g4Tracks().size()-1);
164			//printf("trackId = %i\n",theSimTrack.trackId());
165			SimTrackRef theSimTrackRef = simTidMap_[theSimTrack.trackId()];
166			mithep::MCParticle *outSimParticle = simMap_->GetMit(theSimTrackRef);
167			trackingMap_->Add(theRef, outSimParticle);
168			}
169
170	bendavid	1.1	}
171			}
172
173			mcParticles_->Trim();
174			}
175
176			//--------------------------------------------------------------------------------------------------
177			void FillerMCParticles::ResolveLinks(const edm::Event &event,
178	loizides	1.2	const edm::EventSetup &setup)
179	bendavid	1.1	{
180			// Loop over HepMC particle and resolve their links.
181
182			if (genActive_) {
183
184			Handle<edm::HepMCProduct> hHepMCProduct;
185			GetProduct(genEdmName_, hHepMCProduct, event);
186
187			const HepMC::GenEvent &GenEvent = hHepMCProduct->getHepMCData();
188
189			for (HepMC::GenEvent::particle_const_iterator pgen = GenEvent.particles_begin();
190			pgen != GenEvent.particles_end(); ++pgen) {
191
192			HepMC::GenParticle mcPart = (pgen);
193			if(!mcPart) continue;
194
195			//check if genpart has a decay vertex
196			HepMC::GenVertex *dVertex = mcPart->end_vertex();
197			if (!dVertex) continue;
198
199			//find corresponding mithep genparticle parent in association table
200			mithep::MCParticle *genParent = genMap_->GetMit(mcPart->barcode());
201
202			//set decay vertex
203			//division by 10.0 is needed due to HepMC use of mm instead of cm for distance units
204			genParent->SetVertex(dVertex->point3d().x()/10.0,
205			dVertex->point3d().y()/10.0,
206			dVertex->point3d().z()/10.0);
207
208			//loop through daugthers
209	loizides	1.3	for (HepMC::GenVertex::particles_out_const_iterator pgenD =
210			dVertex->particles_out_const_begin();
211			pgenD != dVertex->particles_out_const_end(); ++pgenD) {
212	bendavid	1.1	HepMC::GenParticle mcDaughter = (pgenD);
213			mithep::MCParticle *genDaughter = genMap_->GetMit(mcDaughter->barcode());
214			genParent->AddDaughter(genDaughter);
215			if (!(genDaughter->HasMother()))
216			genDaughter->SetMother(genParent);
217			}
218			}
219			}
220
221	bendavid	1.4	//loop over SimTracks and resolve links
222	bendavid	1.1	if (simActive_) {
223	bendavid	1.4	Handle<edm::SimTrackContainer> hSimTrackProduct;
224			GetProduct(simEdmName_, hSimTrackProduct, event);
225
226			const edm::SimTrackContainer simTracks = *(hSimTrackProduct.product());
227
228			Handle<std::vector<SimVertex> > hSimVertexProduct;
229			GetProduct(simEdmName_, hSimVertexProduct, event);
230
231			const edm::SimVertexContainer simVertexes = *(hSimVertexProduct.product());
232
233
234
235	loizides	1.2
236	bendavid	1.1	// loop through all simParticles
237	bendavid	1.4	for (SimTrackContainer::const_iterator iM = simTracks.begin();
238			iM != simTracks.end(); ++iM) {
239
240			if (iM->vertIndex()>=0) {
241			SimTrackRef theRef(hSimTrackProduct, iM-simTracks.begin());
242			mithep::MCParticle *simDaughter = simMap_->GetMit(theRef);
243			const SimVertex &theVertex = simVertexes.at(iM->vertIndex());
244			if (theVertex.parentIndex()>=0) {
245			mithep::MCParticle *simParent = simMap_->GetMit(simTidMap_[theVertex.parentIndex()]);
246			simParent->SetVertex(theVertex.position().x(),theVertex.position().y(),theVertex.position().z());
247	bendavid	1.6	//make sure we don't double count the decay tree
248			if ( !simParent->HasDaughter(simDaughter) ) {
249			simParent->AddDaughter(simDaughter);
250			}
251			if ( !simDaughter->HasMother() ) {
252			simDaughter->SetMother(simParent);
253			}
254	bendavid	1.1	}
255	bendavid	1.4
256	bendavid	1.1	}
257	bendavid	1.4
258	bendavid	1.1	}
259	bendavid	1.4
260	bendavid	1.1	}
261	bendavid	1.4
262			//Debug output, will be uncommented and switchable later
263			// if (trackingActive_) {
264			//
265			// Handle<TrackingParticleCollection> hTrackingParticleProduct;
266			// GetProduct(trackingEdmName_, hTrackingParticleProduct, event);
267			//
268			// const TrackingParticleCollection trackingParticles = *(hTrackingParticleProduct.product());
269			//
270			// // loop through all simParticles
271			// for (TrackingParticleCollection::const_iterator iM = trackingParticles.begin();
272			// iM != trackingParticles.end(); ++iM) {
273			//
274			// TrackingParticleRef theRef(hTrackingParticleProduct, iM-trackingParticles.begin());
275			//
276			// mithep::MCParticle *simParent = trackingMap_->GetMit(theRef);
277			// printf("MCParticle pdg = %i, Daughter Pdgs: ", simParent->PdgId());
278			// for (UInt_t i=0; i<simParent->NDaughters(); ++i) {
279			// printf("%i, ", simParent->Daughter(i)->PdgId());
280			// }
281			// printf("\n");
282			//
283			// printf("TrackingParticle pdg = %i, Daughter Pdgs: ", iM->pdgId());
284			// for (TrackingVertexRefVector::iterator v= iM->decayVertices().begin();
285			// v != iM->decayVertices().end(); ++v) {
286			// for (TrackingParticleRefVector::iterator pd = v->get()->daughterTracks().begin();
287			// pd != v->get()->daughterTracks().end(); ++pd) {
288			// ;//printf("%i, ", pd->get()->pdgId());
289			// }
290			// }
291			// printf("\n");
292			// }
293			// }
294	bendavid	1.1	}