All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
RecoMCTruthLinker.cc
Go to the documentation of this file.
1 #include "RecoMCTruthLinker.h"
2 #include <iostream>
3 
4 #include <cstdlib>
5 
6 // LCIO
7 #include <EVENT/LCCollection.h>
8 #include <EVENT/Track.h>
9 #include <EVENT/Cluster.h>
10 #include <EVENT/SimCalorimeterHit.h>
11 #include <EVENT/CalorimeterHit.h>
12 #include <EVENT/SimTrackerHit.h>
13 #include <EVENT/TrackerHit.h>
14 #include <EVENT/ReconstructedParticle.h>
15 #include <EVENT/LCRelation.h>
16 #include <UTIL/LCRelationNavigator.h>
17 
18 #include "UTIL/LCTOOLS.h"
19 #include <cstdio>
20 
21 #include "gearimpl/Vector3D.h"
22 
23 #include <math.h>
24 #include <map>
25 #include <algorithm>
26 
27 #ifdef MARLIN_USE_AIDA
28 #include <marlin/AIDAProcessor.h>
29 #include <AIDA/AIDA.h>
30 #endif
31 
32 
33 using namespace lcio ;
34 using namespace marlin ;
35 
36 
37 
38 RecoMCTruthLinker aRecoMCTruthLinker;
39 
40 
41 struct MCPKeep : public LCIntExtension<MCPKeep> {} ;
42 
43 typedef std::map< MCParticle* , int > MCPMap ;
44 
45 typedef std::map< Track* , int > TrackMap ;
46 
47 typedef std::map< MCParticle* , double > MCPMapDouble ;
48 
49 typedef std::map< MCParticle* , MCParticle* > Remap_as_you_go;
50 
51 RecoMCTruthLinker::RecoMCTruthLinker() : Processor("RecoMCTruthLinker") {
52 
53  // modify processor description
54  _description = "links RecontructedParticles to the MCParticle based on number of hits used" ;
55 
56 
57  IntVec pdgVecDef ;
58 
59  pdgVecDef.push_back( 22 ) ; // gamma
60  pdgVecDef.push_back( 111 ) ; // pi0
61  pdgVecDef.push_back( 310 ) ; // K0s
62 
63 
64  // >>>>> reconstructed or simulated thingis
65 
66 
67  registerInputCollection( LCIO::MCPARTICLE,
68  "MCParticleCollection" ,
69  "Name of the MCParticle input collection" ,
70  _mcParticleCollectionName ,
71  std::string("MCParticle") ) ;
72 
73  registerInputCollection( LCIO::TRACK,
74  "TrackCollection" ,
75  "Name of the Tracks input collection" ,
76  _trackCollectionName ,
77  std::string("MarlinTrkTracks") ) ;
78 
79  registerInputCollection( LCIO::CLUSTER,
80  "ClusterCollection" ,
81  "Name of the Clusters input collection" ,
82  _clusterCollectionName ,
83  std::string("PandoraClusters") ) ;
84 
85  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
86  "RecoParticleCollection" ,
87  "Name of the ReconstructedParticles input collection" ,
88  _recoParticleCollectionName ,
89  std::string("PandoraPFOs") ) ;
90 
91  // <<<<<<
92 
93 // >>>>> tracker hits and sim<->seen relations
94 
95  StringVec trackerHitsInputColNamesDefault ;
96  trackerHitsInputColNamesDefault.push_back("VXDCollection ");
97  trackerHitsInputColNamesDefault.push_back("SITCollection");
98  trackerHitsInputColNamesDefault.push_back("FTD_PIXELCollection");
99  trackerHitsInputColNamesDefault.push_back("FTD_STRIPCollection");
100  trackerHitsInputColNamesDefault.push_back("TPCCollection");
101  trackerHitsInputColNamesDefault.push_back("SETCollection");
102 
103  registerInputCollections( LCIO::SIMTRACKERHIT,
104  "SimTrackerHitCollections" ,
105  "Names of the SimTrackerHits input collection" ,
106  _simTrkHitCollectionNames ,
107  trackerHitsInputColNamesDefault ) ;
108 
109  registerProcessorParameter("UseTrackerHitRelations",
110  "true: use relations for TrackerHits, false : use getRawHits ",
111  _use_tracker_hit_relations,
112  bool(true));
113 
114 
115  StringVec trackerHitsRelInputColNamesDefault;
116  trackerHitsRelInputColNamesDefault.push_back( "VXDTrackerHitRelations" );
117  trackerHitsRelInputColNamesDefault.push_back( "SITSpacePointRelations" );
118  trackerHitsRelInputColNamesDefault.push_back( "FTDPixelTrackerHitRelations" );
119  trackerHitsRelInputColNamesDefault.push_back( "FTDSpacePointRelations" );
120  trackerHitsRelInputColNamesDefault.push_back( "TPCTrackerHitRelations" );
121  trackerHitsRelInputColNamesDefault.push_back( "SETSpacePointRelations" );
122 
123 
124  registerInputCollections("LCRelation",
125  "TrackerHitsRelInputCollections",
126  "Name of the lcrelation collections, that link the TrackerHits to their SimTrackerHits.",
127  _colNamesTrackerHitRelations,
128  trackerHitsRelInputColNamesDefault );
129 
130  // <<<<<<
131 
132  // >>>>> Calorimeter hits and sim<->seen relations
133 
134  StringVec caloHitsInputColNamesDefault ;
135  caloHitsInputColNamesDefault.push_back("BeamCalCollection");
136  caloHitsInputColNamesDefault.push_back("EcalBarrelSiliconCollection");
137  caloHitsInputColNamesDefault.push_back("EcalBarrelSiliconPreShowerCollection");
138  caloHitsInputColNamesDefault.push_back("EcalEndcapRingCollection");
139  caloHitsInputColNamesDefault.push_back("EcalEndcapRingPreShowerCollection");
140  caloHitsInputColNamesDefault.push_back("EcalEndcapSiliconCollection");
141  caloHitsInputColNamesDefault.push_back("EcalEndcapSiliconPreShowerCollection");
142  caloHitsInputColNamesDefault.push_back("HcalBarrelRegCollection");
143  caloHitsInputColNamesDefault.push_back("HcalEndCapRingsCollection");
144  caloHitsInputColNamesDefault.push_back("HcalEndCapsCollection");
145  caloHitsInputColNamesDefault.push_back("LHcalCollection");
146  caloHitsInputColNamesDefault.push_back("LumiCalCollection");
147  caloHitsInputColNamesDefault.push_back("MuonBarrelCollection");
148  caloHitsInputColNamesDefault.push_back("MuonEndCapCollection");
149 
150  registerInputCollections( LCIO::SIMCALORIMETERHIT,
151  "SimCaloHitCollections" ,
152  "Names of the SimCaloHits input collections" ,
153  _simCaloHitCollectionNames ,
154  caloHitsInputColNamesDefault ) ;
155 
156  StringVec caloHitsRelInputColNamesDefault;
157  caloHitsRelInputColNamesDefault.push_back( "RelationCaloHit" );
158  caloHitsRelInputColNamesDefault.push_back( "RelationLcalHit" );
159  caloHitsRelInputColNamesDefault.push_back( "RelationLHcalHit" );
160  caloHitsRelInputColNamesDefault.push_back( "RelationBCalHit" ) ;
161  caloHitsRelInputColNamesDefault.push_back( "RelationMuonHit" );
162 
163  registerInputCollections( "LCRelation",
164  "SimCalorimeterHitRelationNames" ,
165  "Name of the lcrelation collections, that link the SimCalorimeterHit to CalorimeterHit" ,
166  _caloHitRelationNames ,
167  caloHitsRelInputColNamesDefault );
168 
169  // <<<<<<
170 
171  // >>>> Output true<->seen realations
172 
173  registerOutputCollection( LCIO::LCRELATION,
174  "MCTruthTrackLinkName" ,
175  "Name of the trackMCTruthLink output collection" ,
176  _mCTruthTrackLinkName ,
177  std::string("") ) ;
178 
179 
180  registerOutputCollection( LCIO::LCRELATION,
181  "TrackMCTruthLinkName" ,
182  "Name of the trackMCTruthLink output collection - not created if empty()" ,
183  _trackMCTruthLinkName ,
184  std::string("") ) ;
185 
186 
187  registerOutputCollection( LCIO::LCRELATION,
188  "ClusterMCTruthLinkName" ,
189  "Name of the clusterMCTruthLink output collection - not created if empty()" ,
190  _clusterMCTruthLinkName ,
191  std::string("") ) ;
192 
193 
194  registerOutputCollection( LCIO::LCRELATION,
195  "MCTruthClusterLinkName" ,
196  "Name of the MCTruthClusterLink output collection" ,
197  _mCTruthClusterLinkName ,
198  std::string("") ) ;
199 
200 
201  registerProcessorParameter("FullRecoRelation",
202  "true: All reco <-> true relations are given, with weight = 10000*calo weight+"
203  "track weight (weights in permill). false: Only highest contributor linked,"
204  "and only to tracks, not clusters if there are any tracks",
205  _FullRecoRelation,
206  bool(false)
207  );
208 
209  registerOutputCollection( LCIO::LCRELATION,
210  "RecoMCTruthLinkName" ,
211  "Name of the RecoMCTruthLink output collection - not created if empty()" ,
212  _recoMCTruthLinkName ,
213  std::string("RecoMCTruthLink") ) ;
214 
215  registerOutputCollection( LCIO::LCRELATION,
216  "MCTruthRecoLinkName" ,
217  "Name of the MCTruthRecoLink output collection" ,
218  _mCTruthRecoLinkName ,
219  std::string("") ) ;
220 
221 
222  registerOutputCollection( LCIO::LCRELATION,
223  "CalohitMCTruthLinkName" ,
224  "Name of the updated calo-hit MCTruthLink output collection - not created if empty()" ,
225  _calohitMCTruthLinkName ,
226  std::string("") ) ;
227 
228  // <<<<<<<
229  // Output skimmed MCParticles
230 
231  registerProcessorParameter( "KeepDaughtersPDG" ,
232  "PDG codes of particles of which the daughters will be "
233  "kept in the skimmmed MCParticle collection" ,
234  _pdgVec ,
235  pdgVecDef ) ;
236 
237  registerOutputCollection( LCIO::MCPARTICLE,
238  "MCParticlesSkimmedName" ,
239  "Name of the skimmed MCParticle output collection - not created if empty()" ,
240  _mcParticlesSkimmedName ,
241  std::string("") ) ;
242  // <<<<<<<
243 
244  // various steering parameters
245 
246  registerProcessorParameter( "daughtersECutMeV" ,
247  "energy cut for daughters that are kept from KeepDaughtersPDG" ,
248  _eCutMeV,
249  float( 10. )
250  ) ;
251 
252 
253  registerProcessorParameter( "SaveBremsstrahlungPhotons" ,
254  "save photons from Brems" ,
255  _saveBremsstrahlungPhotons,
256  bool(false)
257  ) ;
258 
259  registerProcessorParameter( "UsingParticleGun" ,
260  "If Using Particle Gun Ignore Gen Stat" ,
261  _using_particle_gun,
262  bool(false)
263  ) ;
264 
265 
266 
267  registerProcessorParameter( "BremsstrahlungEnergyCut" ,
268  "energy cut for Brems that are kept" ,
269  _bremsstrahlungEnergyCut,
270  float( 1. )
271  ) ;
272 
273  registerProcessorParameter( "InvertedNonDestructiveInteractionLogic" ,
274  "Work-around Mokka bug in vertex-is-not-endpoint-of-parent flag (logic inverted)" ,
275  _invertedNonDestructiveInteractionLogic,
276  bool(false)
277  ) ;
278 
279 
280 
281 }
282 
283 
284 void RecoMCTruthLinker::init() {
285 
286  // usually a good idea to
287  printParameters<MESSAGE>() ;
288 
289  for( IntVec::iterator it = _pdgVec.begin() ; it != _pdgVec.end() ; ++it ){
290 
291  if( *it < 0 ) {
292 
293  streamlog_out( WARNING ) << " init: negative PDG given - only abs value is used : " << *it << std::endl ;
294  }
295 
296  _pdgSet.insert( abs( *it ) ) ;
297  }
298 
299  // don't write outpur collections that have an empty name
300  _OutputClusterTruthRelation = ! _clusterMCTruthLinkName.empty() ;
301  _OutputTruthClusterRelation = ! _mCTruthClusterLinkName.empty() ;
302  _OutputCalohitRelation = ! _calohitMCTruthLinkName.empty() ;
303  _OutputTrackTruthRelation = ! _trackMCTruthLinkName.empty() ;
304  _OutputTruthTrackRelation = ! _mCTruthTrackLinkName.empty() ;
305  _OutputTruthRecoRelation = ! _mCTruthRecoLinkName.empty() ;
306 
307 
308  if( ! _use_tracker_hit_relations ){
309 
310  streamlog_out( WARNING ) << " ====== UseTrackerHitRelations=false => not using TrackerHit-SimTrackerHit-Relations but getRawHits() instead - \n"
311  << " this is probably not what you want (only for backward compatibility with very old files...)" << std::endl ;
312  }
313 
314  _nRun = 0 ;
315  _nEvt = 0 ;
316 }
317 
318 
319 void RecoMCTruthLinker::processRunHeader( LCRunHeader* /*run*/) {
320 
321  _nRun++ ;
322 }
323 
324 
325 void RecoMCTruthLinker::processEvent( LCEvent * evt ) {
326 
327  _nEvt ++ ;
328 
329 
330  streamlog_out( MESSAGE0 ) << " processEvent " << evt->getEventNumber() << " - " << evt->getRunNumber()
331  << std::endl ;
332 
333  LCCollection* mcpCol = evt->getCollection( _mcParticleCollectionName ) ;
334 
335 
336  // find track to MCParticle relations
337 
338  LCCollection* trackCol = 0 ;
339  LCCollection* ttrlcol = 0 ;
340  LCCollection* trtlcol = 0 ;
341 
342  bool haveTracks = true ;
343 
344  try{ trackCol = evt->getCollection( _trackCollectionName ) ; } catch(DataNotAvailableException&){ haveTracks=false ; }
345 
346  if( ! haveTracks ) {
347  streamlog_out( DEBUG9 ) << " Track collection : " << _trackCollectionName
348  << " not found - cannot create relation " << std::endl ;
349 
350 
351  } else { //if ( haveTracks ) {
352 
353  trackLinker( evt, mcpCol , trackCol , &ttrlcol , &trtlcol);
354 
355  if (_OutputTrackTruthRelation )
356  evt->addCollection( ttrlcol , _trackMCTruthLinkName ) ;
357  if (_OutputTruthTrackRelation )
358  evt->addCollection( trtlcol , _mCTruthTrackLinkName ) ;
359  }
360 
361  // find cluster to MCParticle relations and the updated calohit to MCParticle relations.
362 
363  LCCollection* clusterCol = 0 ;
364 
365  bool haveClusters = true ;
366  bool haveCaloHitRel = true ;
367 
368  LCCollection* ctrlcol = 0;
369  LCCollection* trclcol = 0;
370  LCCollection* chittrlcol = 0;
371 
372  try{ clusterCol = evt->getCollection( _clusterCollectionName ) ; } catch(DataNotAvailableException&){ haveClusters= false ; }
373 
374  if( ! haveClusters ) {
375  streamlog_out( DEBUG9 ) << " Cluster collection : " << _clusterCollectionName
376  << " not found - cannot create relation " << std::endl ;
377  }
378 
379 
380  if( haveTracks && haveClusters ) {
381 
382  clusterLinker( evt, mcpCol , clusterCol,
383  &ctrlcol , &trclcol , &chittrlcol );
384 
385  if (_OutputClusterTruthRelation ) evt->addCollection( ctrlcol , _clusterMCTruthLinkName ) ;
386  if (_OutputTruthClusterRelation ) evt->addCollection( trclcol , _mCTruthClusterLinkName );
387  if (_OutputCalohitRelation ) evt->addCollection( chittrlcol , _calohitMCTruthLinkName ) ;
388  }
389 
390 
391 
392 
393  // combine track and cluster to MCParticle relations to the reconstructed particle
394  // to MCParticle relation
395  LCCollection* particleCol = 0 ;
396  bool haveRecoParticles = true ;
397 
398  try { particleCol = evt->getCollection( _recoParticleCollectionName ); } catch(DataNotAvailableException&){ haveRecoParticles = false ; }
399 
400  if( ! haveRecoParticles ) {
401  streamlog_out( DEBUG9 ) << " ReconstructedParticle collection : " << _recoParticleCollectionName
402  << " not found - cannot create relation " << std::endl ;
403  }
404 
405  LCCollection* ptrlcol = 0;
406  LCCollection* trplcol = 0;
407  if( haveRecoParticles && haveTracks && haveClusters && haveCaloHitRel && !_recoMCTruthLinkName.empty() ) {
408 
409  particleLinker( mcpCol , particleCol, ttrlcol, ctrlcol, trtlcol, trclcol, &ptrlcol, &trplcol);
410 
411  evt->addCollection( ptrlcol , _recoMCTruthLinkName ) ;
412  if (_OutputTruthRecoRelation ) evt->addCollection( trplcol , _mCTruthRecoLinkName ) ;
413 
414  }
415 
416  if( haveTracks && haveClusters && haveCaloHitRel && !_mcParticlesSkimmedName.empty() ) {
417 
418  LCCollectionVec* skimVec = new LCCollectionVec( LCIO::MCPARTICLE ) ;
419 
420  makeSkim( mcpCol , ttrlcol, ctrlcol , &skimVec );
421  evt->addCollection( skimVec , _mcParticlesSkimmedName ) ;
422  if( streamlog_level(DEBUG5) ){ linkPrinter ( skimVec , particleCol, ptrlcol, trplcol ); }
423  }
424 
425  //If either collection has not been added to the event, we have to delete it now!
426  //Don't delete them before, because they are used
427  if(!_OutputClusterTruthRelation) { delete ctrlcol; }
428  if(!_OutputTruthClusterRelation) { delete trclcol; }
429  if(!_OutputTrackTruthRelation) { delete ttrlcol; }
430  if(!_OutputTruthTrackRelation) { delete trtlcol; }
431  if(!_OutputCalohitRelation) { delete chittrlcol; }
432  if(!_OutputTruthRecoRelation ) { delete trplcol ; }
433 
434  //cleanup relationNavigators and transient collections
435  delete _navMergedTrackerHitRel; _navMergedTrackerHitRel = nullptr;
436  delete _navMergedCaloHitRel; _navMergedCaloHitRel = nullptr;
437  delete _mergedCaloHitRelCol; _mergedCaloHitRelCol = nullptr;
438  delete _mergedTrackerHitRelCol; _mergedTrackerHitRelCol = nullptr;
439 
440 
441 }
442 void RecoMCTruthLinker::trackLinker( LCEvent * evt, LCCollection* mcpCol , LCCollection* trackCol,
443  LCCollection** ttrlcol, LCCollection** trtlcol) {
444 
445 
446 
447  // merge all the SimTrackerHit - TrackerHit relations into on collection and set up the combined navigator
448  mergeTrackerHitRelations(evt);
449 
450 
451  LCRelationNavigator trackTruthRelNav(LCIO::TRACK , LCIO::MCPARTICLE ) ;
452 
453  // the inverse relation from MCTruth particles to tracks
454  // weight is the realtive number of hits from a given MCParticle on the track
455  LCRelationNavigator truthTrackRelNav(LCIO::MCPARTICLE , LCIO::TRACK ) ;
456 
457  //========== fill a map with #SimTrackerHits per MCParticle ==================
458  MCPMap simHitMap ; // counts total simhits for every MCParticle
459  for( unsigned i=0,iN=_simTrkHitCollectionNames.size() ; i<iN ; ++i){
460 
461  const LCCollection* col = 0 ;
462  try{ col = evt->getCollection( _simTrkHitCollectionNames[i] ) ; } catch(DataNotAvailableException&) {}
463  if( col )
464  for( int j=0, jN= col->getNumberOfElements() ; j<jN ; ++j ) {
465 
466  SimTrackerHit* simHit = (SimTrackerHit*) col->getElementAt( j ) ;
467  MCParticle* mcp = simHit->getMCParticle() ;
468  simHitMap[ mcp ] ++ ;
469  }
470  }
471  //===========================================================================
472 
473  // loop over reconstructed tracks
474  int nTrack = trackCol->getNumberOfElements() ;
475 
476  int ifoundch =0;
477 
478  streamlog_out( DEBUG6 ) << " *** Sorting out Track<->MCParticle using simHit<->MCParticle." << std::endl;
479 
480  for(int i=0;i<nTrack;++i){
481 
482  Track* trk = dynamic_cast<Track*> ( trackCol->getElementAt(i) ) ;
483 
484  // charged particle or V0 -> analyse track. We need to find all seen hits
485  // this track is made of, wich sim hits each of the seen hits came from,
486  // and finally which true particle actually created each sim hit
487 
488  MCPMap mcpMap ; // mcpMap is a map seen <-> true particle
489 
490  int nSimHit = 0 ;
491 
492  const TrackerHitVec& trkHits = trk->getTrackerHits() ;
493 
494  for( TrackerHitVec::const_iterator hitIt = trkHits.begin() ; hitIt != trkHits.end() ; ++hitIt ) {
495 
496  TrackerHit* hit = * hitIt ; // ... and a seen hit ...
497 
498  const LCObjectVec& simHits = _use_tracker_hit_relations ? *(this->getSimHits(hit)) : hit->getRawHits() ;
499  MCParticle* mcp2 = 0;
500  for( LCObjectVec::const_iterator objIt = simHits.begin() ; objIt != simHits.end() ; ++objIt ) {
501 
502 
503  SimTrackerHit* simHit = dynamic_cast<SimTrackerHit*>( *objIt ) ; // ...and a sim hit ...
504 
505  MCParticle* mcp = simHit->getMCParticle() ; // ... and a true particle !
506  if ( simHits.size() > 1 ) {
507  if ( mcp2 != 0 && mcp2 != mcp ) {
508  // In this case, the mcp:s will count double !!
509  streamlog_out( DEBUG3 ) << " ghost/double " << mcp << " " << mcp2 << " " << hit->getCellID0() <<std::endl;
510  }
511  mcp2=mcp;
512  }
513  if ( mcp != 0 ) {
514  mcpMap[ mcp ]++ ; // count the hit caused by this true particle
515  } else {
516  streamlog_out( WARNING ) << " tracker SimHit without MCParticle ?! " << std::endl ;
517  }
518 
519  ++nSimHit ; // total hit count
520 
521  }
522  }
523 
524 
525  if( nSimHit == 0 ){
526 
527  streamlog_out( WARNING ) << " No simulated tracker hits found. Set UseTrackerHitRelations to true in steering file to enable using TrackerHit relations if they are available." << std::endl ;
528  streamlog_out( WARNING ) << trk->id() <<" " << i << " " << trk->getTrackerHits().size() << std::endl ;
529 
530  for( TrackerHitVec::const_iterator hitIt = trkHits.begin() ; hitIt != trkHits.end() ; ++hitIt ) {
531 
532  TrackerHit* hit = * hitIt ; // ... and a seen hit ...
533  streamlog_out( WARNING ) << hit->getPosition()[0] << " " << hit->getPosition()[1] << " " << hit->getPosition()[2] << " " << std::endl ;
534  const LCObjectVec& simHits = *(this->getSimHits(hit)) ;
535  streamlog_out( WARNING ) << simHits.size() << std::endl ;
536  }
537  continue ; // won't find a particle
538  }
539 
540 
541  // find the mc particle with the largest #hits
542  // also store all genstat=1 particles and
543  // all particles of type to be saved
544 
545  MCParticle* mother = 0;
546  MCParticleVec theMCPs ; // vector that will contain all true particles contributing
547  theMCPs.reserve( 1000 ) ;
548  std::vector<int> MCPhits;
549  MCPhits.reserve( 1000 ) ;
550  int ifound = 0;
551 
552  for( MCPMap::iterator it = mcpMap.begin() ; // iterate trough the map, map->first is the
553  it != mcpMap.end() ; ++it ){ // true particle, map->second is the number of
554  // times this true particle got mapped, ie. the
555  // number of hits it produced.
556 
557 
558  mother = ( it->first->getParents().size()!=0 ? dynamic_cast<MCParticle*>(it->first->getParents()[0]) : 0 ) ; // mother of the true particle.
559 
560  if ( _using_particle_gun || it->first->getGeneratorStatus() == 1 ) { // genstat 1 particle, ie. it is a bona fide
561  // creating true particle: enter it into the list,
562  // and note how many hits it produced.
563  theMCPs.push_back( it->first ) ;
564  MCPhits.push_back( it->second ) ;
565  ifound++;
566 
567  } else { // not genstat 1. Wat should we do with it ?
568 
569  if ( mother != 0 ) { // if it has a parent, save it
570 
571  theMCPs.push_back(it->first);
572  MCPhits.push_back(it->second);
573  ifound++;
574 
575  } else {
576 
577  streamlog_out( WARNING ) << " track has hit(s) from a non-generator particle with no parents ?!! " << std::endl;
578  }
579  }
580  } // end of loop over map
581 
582  // hits in track
583  unsigned nHit = 0;
584 
585  for (unsigned ihit=0; ihit<trkHits.size(); ++ihit) {
586  if ( UTIL::BitSet32( trkHits[ihit]->getType() )[ UTIL::ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT ]) {
587  nHit += 2;
588  } else {
589  nHit += 1;
590  }
591  }
592 
593  // finally calculate the weight of each true particle to the seen
594  // (= hits_from_this_true/ all_hits),
595  // and add the weighted reltion
596  for (int iii=0 ; iii<ifound ; iii++ ) {
597 
598 
599  float weight = float(MCPhits[iii] )/float(nHit) ;
600 
601  trackTruthRelNav.addRelation( trk , theMCPs[iii] , weight ) ;
602 
603  int Total_SimHits_forMCP = simHitMap[ theMCPs[iii] ];
604 
605  float inv_weight = float(MCPhits[iii] ) / Total_SimHits_forMCP ;
606 
607  truthTrackRelNav.addRelation( theMCPs[iii] , trk , inv_weight ) ;
608 
609 
610  streamlog_out( DEBUG4 ) << " track " << trk->id() << " has " << MCPhits[iii] << " hits of "
611  << nSimHit << " SimHits ("
612  << nHit << " TrackerHits) "
613  << " weight = " << weight << " , "
614  << " inv rel weight = " << inv_weight << " "
615  << " of MCPart with pdg : " << theMCPs[iii]->getPDG()
616  << " , SimHits " << Total_SimHits_forMCP
617  << " and genstat : " << theMCPs[iii]->getGeneratorStatus()
618  << " id: " << theMCPs[iii]->id()
619  << std::endl ;
620 
621 
622 
623  }
624 
625  ifoundch=ifound;
626  }
627  // seen-true relation complete. add the collection
628 
629  streamlog_out( DEBUG6 ) << " *** Sorting out Track<->MCParticle : DONE " << std::endl;
630  streamlog_out( DEBUG6 ) << " *** track linking complete, create collection " << std::endl;
631 
632  *trtlcol = truthTrackRelNav.createLCCollection() ;
633  *ttrlcol = trackTruthRelNav.createLCCollection() ;
634 
635 }
636 
637 
638 
639 void RecoMCTruthLinker::clusterLinker( LCEvent * evt, LCCollection* mcpCol , LCCollection* clusterCol,
640  LCCollection** ctrlcol, LCCollection** trclcol,LCCollection** chittrlcol) {
641 
642 
643 
644 
645 
646  mergeCaloHitRelations(evt);
647  LCRelationNavigator clusterTruthRelNav(LCIO::CLUSTER , LCIO::MCPARTICLE ) ;
648  LCRelationNavigator truthClusterRelNav( LCIO::MCPARTICLE , LCIO::CLUSTER ) ;
649  LCRelationNavigator chitTruthRelNav(LCIO::CALORIMETERHIT , LCIO::MCPARTICLE ) ;
650 
651 
652  Remap_as_you_go remap_as_you_go ; // map from true tracks linked to hits to
653  // those that really should have been linked
654  MCPMapDouble simHitMapEnergy ; // counts total simhits for every MCParticle
655 
656 
657  streamlog_out( DEBUG6 ) << " *** Sorting out simHit<->MCParticle connections, and find corresponding calo hits " << std::endl;
658 
659  // This step is needed for cluster (contrary to the tracks above), because the first few branchings in the
660  // shower is kept in the MCParticle collection. We therefor need to back-track to the particle actually
661  // entering the calorimeter. The originator of a hit sometimes IS a generator particle, in which case
662  // there is no problem. Otherwise, the criteria to keep a particle as a hit originator is:
663  // The particle is an ancestor of the one initial linked to the hit.
664  // AND
665  // [
666  // The particle is a generator particle
667  // OR
668  // The particle starts in the tracker, ends in calorimeter. That it starts in the tracker
669  // is determined by the fact that it's mother ends there. However, note the posibility of
670  // "non-destuctive interactions", where a particles mother does *not* end at the point where
671  // the particle is created !
672  // OR
673  // The particle, or one of it's ancestors is a back-scatterer. This case is further treated in the
674  // loop over clusters, to catch the rather common case that a back-scattered particle re-enters the
675  // calorimeter close to it's start-start point (think 15 MeV charged particle!) and the hits it
676  // produces is in the same cluster as that of the initiator of the shower the backscatter come from.
677  // (This can't be detected in this first loop, since we don't know about clusters here)
678  // ]
679  // A map is set up, so obviously the first check is wether the MCP of the hit has already been
680  // treated when looking at some previous hit, in which case one just thake that association.
681 
682 
683  for( unsigned i=0,iN=_simCaloHitCollectionNames.size() ; i<iN ; ++i){
684  const LCCollection* col = 0 ;
685  try{ col = evt->getCollection( _simCaloHitCollectionNames[i] ) ; } catch(DataNotAvailableException&) {}
686  if( col ) {
687 
688  streamlog_out( DEBUG6 ) << std::endl;
689  streamlog_out( DEBUG6 ) << " ================= " << std::endl;
690  streamlog_out( DEBUG6 ) << " Treating sim-hits in " << _simCaloHitCollectionNames[i] << ". It has "
691  << col->getNumberOfElements() << " hits " << std::endl;
692 
693 
694  for( int j=0, jN= col->getNumberOfElements() ; j<jN ; ++j ) {
695 
696  SimCalorimeterHit* simHit = (SimCalorimeterHit*) col->getElementAt( j ) ;
697  LCObjectVec caloHits = _navMergedCaloHitRel->getRelatedFromObjects(simHit);
698  if ( caloHits.size() == 0 ) { continue ;}
699  if ( caloHits.size() != 1 ) { streamlog_out( DEBUG9 ) << " Sim hit with nore than one calo hit ? " << std::endl; }
700 
701  CalorimeterHit* caloHit = dynamic_cast<CalorimeterHit*>(caloHits[0]);
702 
703  // for split hits, need to get calib factor over sum of all rec hits (DJ)
704  double allRecoEn(0);
705  for ( size_t jk=0; jk<caloHits.size(); jk++) {
706  allRecoEn+= ((CalorimeterHit*) caloHits[jk])->getEnergy();
707  }
708  double calib_factor = allRecoEn/simHit->getEnergy();
709 
710  streamlog_out( DEBUG4 ) << std::endl;
711  streamlog_out( DEBUG4 ) << " ================= " << std::endl;
712  streamlog_out( DEBUG4 ) << " Treating hit " << j << " sim hit id " << simHit->id() << " calo hit id "
713  << caloHit->id() << " nb contributions " << simHit->getNMCContributions() << std::endl;
714  streamlog_out( DEBUG3 ) << " ncalo hits : " << caloHits.size() << " calib factor " << calib_factor
715  << " position " << simHit->getPosition()[0] << " " <<
716  simHit->getPosition()[1] << " " <<
717  simHit->getPosition()[2] << " " << std::endl;
718 
719  for(int k=0;k<simHit->getNMCContributions() ;k++){
720  MCParticle* mcp = simHit->getParticleCont( k ) ;
721 
722  streamlog_out( DEBUG2 ) << " a calo hit to treat " << std::endl;
723 
724  if ( mcp == 0 ) {
725 
726 
727  streamlog_out( DEBUG6 ) <<" collection nb " << i << " nhits " << jN << " hit " << j << " contr " << k << std::endl;
728  streamlog_out( DEBUG6 ) <<" Sim hit with no mcp. N contrib is " <<
729  simHit->getNMCContributions() <<std::endl;
730  streamlog_out( DEBUG6 ) <<" Hit Position: " << simHit->getPosition()[0] << " " <<
731  simHit->getPosition()[1] << " " <<
732  simHit->getPosition()[2] << " " << std::endl;
733  if ( simHit->getNMCContributions() > 0 ) {
734  for (int l=0;l<simHit->getNMCContributions() ;l++){
735  streamlog_out( DEBUG6 ) <<" Contrib " << l << " : " <<std::endl;
736  streamlog_out( DEBUG6 ) <<" Energy: " << simHit->getEnergyCont(l) <<std::endl;
737  streamlog_out( DEBUG6 ) <<" PDG: " << simHit->getPDGCont(l) <<std::endl;
738  streamlog_out( DEBUG6 ) <<" MCPart: " << simHit->getParticleCont(l) <<std::endl;
739  }
740  }
741 
742 
743 
744  continue;
745  }
746 
747  double e = simHit->getEnergyCont( k ) * calib_factor;
748 
749  streamlog_out( DEBUG3 ) <<" initial true contributor "<< k << " id " << mcp->id()
750  <<" (with E = " << mcp->getEnergy() << " and pdg " << mcp->getPDG() << " ) e hit: " << e <<std::endl;
751 
752  if ( remap_as_you_go.find(mcp) != remap_as_you_go.end() ) {
753  // very first condition: I already know what to do from some earlier hit created by mcp.
754  mcp=remap_as_you_go.find(mcp)->second;
755  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
756  << caloHit->id() << " attributed to " << mcp->id() <<
757  " because it's origin mcp has already been treated : originator case 0 " <<std::endl;
758  } else {
759  MCParticle* mother = 0;
760  MCParticle* this_Kid = mcp ; // ... and a true particle !
761 
762  //Particle gun particles dont have parents, but are "created in the simulation" and have genStat 0
763  if ( mcp-> getGeneratorStatus() == 0 && ( mcp->getParents().empty() == false ) ) {
764 
765  // not from generator, find which true particle this
766  // hit should really be attributed to, by tracking back the history.
767  // (the other case, ie. if the if above is false is "case 1")
768 
769  // Two cases to treat:
770  // For some reason, the hit is not attributed to the
771  // incomming particle, but to some particle in the shower.
772  // Just track back to the incomming particle, which might (usually)
773  // be a gen-stat 1 particle from the main vertex. It can also
774  // be from a decay or interaction in the tracking made by
775  // Geant : genstat=0, mother isDecayedInTracker, (or
776  // the production vertex is in the track, but the mother
777  // continues on - ticky case, see comment futher down), or a
778  // decyed particle from the generator (genstat 2) that
779  // hits the calo before decaying (lambdas, K^0_S)
780 
781  // Or: for back-scatterers the calorimiter hit is attributed to the
782  // the last particle *even if this particle both
783  // started and ended in the tracker* !!!! Then we back-track
784  // untill we find a particle which at least started inside
785  // the calorimeter, and then go on backtracking as above.
786  // This case is triggered by the particle linked to the
787  // hit being DecayedInTracker, hence the case where a
788  // back-scatter actually ends in the calorimeter is
789  // treated as the first case.
790 
791 
792 
793  streamlog_out( DEBUG2 ) << " simHit " <<simHit->id()<<","<< j <<
794  " not created by generator particle. backtracking ..."
795  << std::endl;
796  streamlog_out( DEBUG2 ) <<" "<<mcp->id()<<" gs "<<mcp->getGeneratorStatus()<<
797  " dint "<<mcp->isDecayedInTracker()<<
798  " bs "<<mcp->isBackscatter()<<
799  " ndi "<<mcp->vertexIsNotEndpointOfParent()<<
800  " npar "<<mcp->getParents().size()<<
801  " pdg "<<mcp->getPDG()<<
802  " "<< mcp->getVertex()[0]<<
803  " "<<mcp->getVertex()[1]<<
804  " "<<mcp->getVertex()[2]<<
805  " "<< mcp->getEndpoint()[0]<<
806  " "<<mcp->getEndpoint()[1]<<
807  " "<<mcp->getEndpoint()[2]<<std::endl;
808 
809 
810 
811  mother= dynamic_cast<MCParticle*>(mcp->getParents()[0]);
812 
813 
814  if ( !this_Kid->isBackscatter() && mother!= 0 &&
815  mother->getParents().size()>0 &&
816  mother->getGeneratorStatus() ==0 &&
817  !mother->isDecayedInTracker() ) {
818 
819  streamlog_out( DEBUG2 ) << " goes into originator loop " << std::endl;
820  }
821  while ( !this_Kid->isBackscatter() && mother!= 0 && mother->getParents().size()>0 &&
822  mother->getGeneratorStatus() ==0 &&
823  !mother->isDecayedInTracker() ) {
824  // back-track as long as there is a non-generator
825  // mother, or the mother decayed in the tracker
826  // (=> the kid is the particle entering the calorimeter.)
827 
828  // case shower-particle
829  streamlog_out( DEBUG1 ) <<" in originator loop " << std::endl;
830 
831  if ( this_Kid->vertexIsNotEndpointOfParent() != _invertedNonDestructiveInteractionLogic) {
832  MCParticle* oma=dynamic_cast<MCParticle*>(mother->getParents()[0]);
833  if ( oma->isDecayedInTracker() ) {
834  streamlog_out( DEBUG1 ) <<" break out : gandmother "<<oma->id()<<
835  " gs "<<oma->getGeneratorStatus()<<
836  " dint "<<oma->isDecayedInTracker()<<
837  " bs "<<oma->isBackscatter()<<
838  " ndi "<<oma->vertexIsNotEndpointOfParent()<<
839  " npar "<<oma->getParents().size()<<
840  " pdg "<<oma->getPDG()<<std::endl;
841  break ;
842  }
843  }
844  this_Kid=mother ;
845  mother= dynamic_cast<MCParticle*>(mother->getParents()[0]); // (assume only one...)
846 
847  streamlog_out( DEBUG1 ) <<" shower-part mother "<<mother->id()<<
848  " gs "<<mother->getGeneratorStatus()<<
849  " dint "<<mother->isDecayedInTracker()<<
850  " bs "<<mother->isBackscatter()<<
851  " ndi "<<mother->vertexIsNotEndpointOfParent()<<
852  " npar "<<mother->getParents().size()<<
853  " pdg "<<mother->getPDG()<<std::endl;
854 
855 
856 
857  }
858 
859  // Further treatment (basically determining if it this_Kid or mother that enetered the
860  // calorimeter) based on why we left the while-loop
861 
862  // here one of the while conditions is false, ie. at least one of
863  // " kid is back-scatter", "no mother" , "mother has no parents", "mother is from
864  // generator", or "mother did decay in tracker" must be true. We know that this_Kid
865  // fulfills all the while-conditions except the first: obvious if at least one iteration of
866  // the loop was done, since it was the mother in the previous iteration,
867  // but also true even if the loop wasn't transversed, due to the conditions
868  // to at all enter this block of code (explicitly must be simulator particle with
869  // mother, implicitly must have ended in the calo, since it did make calo hits.)
870  if (this_Kid->isBackscatter() ) {
871  // case 2: Kid is back-scatterer. It has thus started in a calo, and entered
872  // from there into the tracking volume, and did cause hits after leaving the tracker
873  // volume again -> this_Kid started before the calo, and is the one
874  remap_as_you_go[mcp]=this_Kid;
875  mcp=this_Kid;
876  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
877  << caloHit->id() << " attributed to kid " << mcp->id() <<
878  " because it's origin is a back-scatter : originator case 2 " <<std::endl;
879  streamlog_out( DEBUG2 ) <<" "<<this_Kid->id()<<
880  " gs "<<this_Kid->getGeneratorStatus()<<
881  " dint "<<this_Kid->isDecayedInTracker()<<
882  " bs "<<this_Kid->isBackscatter()<<
883  " npar "<<this_Kid->getParents().size()<<
884  " pdg "<<this_Kid->getPDG()<<std::endl;
885 
886 
887  } else if ( mother->isDecayedInTracker() ) { // the clear-cut case:
888  remap_as_you_go[mcp]=this_Kid;
889  mcp=this_Kid; // this_Kid started before the calo, and is the one
890  // the hit should be attributed to
891 
892 
893  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
894  << caloHit->id() << " attributed to kid " << mcp->id() <<
895  " because it's origin is in tracker : originator case 3 " <<std::endl;
896 
897 
898  } else { // the other three cases, ie. one or sveral of "no mother", "no grand-mother",
899  // "generator particle" + that we know that "mother decayed in calo"
900  if ( mother == 0 ) { // ... which of course implies no grand-mother, and no gen stat
901  // of the mother as well -> should not be possible !
902 
903  streamlog_out( WARNING ) << " MCparticle " << this_Kid->id() <<
904  " is a simulation particle, created in the calorimeter by nothing . "<< std::endl;
905 
906 
907  remap_as_you_go[mcp]=this_Kid;
908  mcp=this_Kid; // can't do better than that.
909 
910  } else { // here we know: "mother exists, but decayed in calo". In addition, two posibilities:
911  // mother is generator particle, or there was no grand-parents. One or both
912  // must be true here. Here it gets complicated, because what we want to know is
913  // whether this_Kid started in the tracker or not. Unluckily, we don't know that
914  // directly, we only know where the mother ended. IF ithe mother ended in the tracker,
915  // there is no problem, and has already been treated, but if it ended in the calo, it is
916  // still possible that this_Kid came from a "non-destructive interaction" with the tracke-detector
917  // material. This we now try to figure out.
918 
919  if ( this_Kid->vertexIsNotEndpointOfParent() == _invertedNonDestructiveInteractionLogic ) {
920  // This bizare condition is due to a bug in LCIO (at least for the DBD samples. this_Kid->vertexIsNotEndpointOfParent()
921  // should be true in the "non-destructive interaction", but it isn't: actually it is "false", but is "true" for the
922  // for particles that *do* originate at the end-vertex of their parent. This is a bug in Mokka.
923  // Hence the above ensures that there was NO "non-destructive interaction", and it is clear this_Kid was created at the end-point
924  // of the mother. The mother is either a generator particle (to be saved), or the "Eve" of the decay-chain (or both).
925  // So mother is the one to save and assign the hits to:
926  remap_as_you_go[mcp]=mother;
927  mcp=mother; // mother started at ip, and reached the calo, and is the one
928  // the hit should be attributed to.
929 
930  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
931  << caloHit->id() << " attributed to mother " << mcp->id() <<
932  " because it is a generator particle or started in tracker : originator case 4 " <<std::endl;
933 
934  } else { // here we DO have a "non-destructive interaction". Unluckily, we cant directly know if this took place in the
935  // tracker (in which case we should keep this_Kid as the mcp to save and assign hits to), or not.
936  // We will play a few clean tricks to find the cases where it either certain that the
937  // "non-destructive interation" was in the tracker, or that it was in the calo. This reduces
938  // the number of uncertain cases to play dirty tricks with.
939 
940  // Clean tricks to play: look at the sisters of this_Kid: with some luck one of them is a promptly decaying
941  // particle, eg. a pi0.
942  // This sister will be flagged as decayed in calo/tracker, and from that we know for certain that the
943  // "non-destructive interaction" was in the calo/tracker.
944  // It can also be that one of the sisters is flagged as a back-scatter, which only happens in the calo.
945  // If the particles grand-mother is decayed in calo, and the mother isn't from a "non-destructive interaction",
946  // the "non-destructive interaction" was in the calo.
947 
948  // Finally, we can check the distance of the end-point of the mother (sure to be in the calo) to
949  // the vertex of this_Kid. If this is small, this_Kid *probably* started in the calo.
950 
951 
952  int starts_in_tracker = 0 ;
953  unsigned lll=0 ;
954  int has_pi0 = 0 ;
955  int oma_in_calo = 0;
956  double rdist =0.;
957  int has_bs = 0;
958 
959  for ( unsigned kkk=0 ; kkk < mother->getDaughters().size() ; kkk++ ) {
960  MCParticle* sister = dynamic_cast<MCParticle*>(mother->getDaughters()[kkk]);
961  if ( sister == this_Kid ) continue;
962  if ( abs(sister->getVertex()[0]-this_Kid->getVertex()[0]) > 0.1 ||
963  abs(sister->getVertex()[1]-this_Kid->getVertex()[1]) > 0.1 ||
964  abs(sister->getVertex()[2]-this_Kid->getVertex()[2]) > 0.1 ) continue; // must check that it is the same vertex:
965  // several "non-destructive interactions" can
966  // take place (think delta-rays !)
967  if ( sister->isBackscatter()) {
968  has_bs = 1 ;
969  lll=kkk;
970  break ;
971  } else if ( sister->isDecayedInTracker() ) {
972  starts_in_tracker = 1 ;
973  lll=kkk;
974  break ;
975  }
976  // any pi0:s at all ? (it doesn't matter that we break at the two cases above,
977  // because if we do, it doesn't matter if there are
978  // pi0 sisters or not !)
979  if ( sister->getPDG() == 111 ) {
980  has_pi0 = 1 ;
981  }
982  }
983  // if not already clear-cut, calculate distance vertext to mother end-point
984  if ( starts_in_tracker != 1 && has_bs != 1 && has_pi0 != 1 && oma_in_calo != 1 ) {
985  rdist=sqrt(pow(mother->getEndpoint()[0]-this_Kid->getVertex()[0],2)+
986  pow(mother->getEndpoint()[1]-this_Kid->getVertex()[1],2)+
987  pow(mother->getEndpoint()[2]-this_Kid->getVertex()[2],2));
988  if ( mother->getParents().size() != 0 ) {
989  MCParticle* oma=dynamic_cast<MCParticle*>(mother->getParents()[0]);
990  if ( oma->isDecayedInCalorimeter() ) {
991  oma_in_calo = 1 ;
992  streamlog_out( DEBUG1 ) << " grandmother in calo " << std::endl;
993  }
994  }
995  }
996  streamlog_out( DEBUG1 ) << " " << starts_in_tracker << " " << has_pi0 << " " << has_bs << " " << oma_in_calo << std::endl;
997  // }
998  if ( starts_in_tracker == 1 ) { // this_Kid is a clear-cut hit-originator
999 
1000  remap_as_you_go[mcp]=this_Kid;
1001  mcp=this_Kid; // this_Kid started before the calo, and is the one
1002  // the hit should be attributed to
1003 
1004  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
1005  << caloHit->id() << " attributed to kid " << mcp->id() <<
1006  " because it's origin could be deduced to be in tracker : originator case 5 " <<std::endl;
1007 
1008  streamlog_out( DEBUG2 ) << " Details of case 5: "
1009  << this_Kid->getVertex()[0] << " "
1010  << this_Kid->getVertex()[1] << " " << this_Kid->getVertex()[2] << " "
1011  << mother->getEndpoint()[0] << " "
1012  << mother->getEndpoint()[1] << " " << mother->getEndpoint()[2] << " "
1013  << mother->getGeneratorStatus() << " "
1014  << this_Kid->vertexIsNotEndpointOfParent() << std::endl;
1015  streamlog_out( DEBUG1 ) << " starts in tracker " << std::endl;
1016 
1017  } else if ( has_pi0 != 0 || has_bs != 0 || oma_in_calo != 0 ) { // clear-cut case of this_Kid starting in the calo.
1018  // We do know that the mother
1019  // is a generator particle and/or the "Eve" of the cascade,
1020  // so we should attribute hits to the
1021  // mother and save it.
1022 
1023  remap_as_you_go[mcp]=mother;
1024  mcp=mother; // mother started at ip, and reached the calo, and is the one
1025  // the hit should be attributed to.
1026 
1027  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
1028  << caloHit->id() << " attributed to mother " << mcp->id() <<
1029  " because it's origin could be deduced to be in tracker : originator case 6 " <<std::endl;
1030  streamlog_out( DEBUG2 ) << " Case 6 details: kid starts in calo " << " "
1031  << this_Kid->getVertex()[1] << " " << this_Kid->getVertex()[2] << " " << rdist << " "
1032  << has_pi0 << " " << has_bs << " " << oma_in_calo << std::endl;
1033 
1034 
1035  } else { // un-clear case: no pi0 nor back-scatteres among the sisters to help to decide.
1036  // Use distance this_Kid-startpoint to mother-endpoint. We know that the latter is
1037  // in the calo, so if this is small, guess that the start point of this_Kid is
1038  // also in the calo. Calos are dense, so typically in the case the "non-destructive interaction"
1039  // is in the calo, one would guess that the distance is small, ie. large distance ->
1040  // unlikely that it was in the calo.
1041 
1042  if ( rdist > 200. ) { // guess "non-destructive interaction" not in calo -> this_Kid is originator
1043 
1044  remap_as_you_go[mcp]=this_Kid;
1045  mcp=this_Kid; // this_Kid started before the calo, and is the one
1046  // the hit should be attributed to
1047  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
1048  << caloHit->id() << " attributed to kid " << mcp->id() <<
1049  " because it's origin is guessed to be in tracker : originator case 7 " <<std::endl;
1050  streamlog_out( DEBUG2 ) << " Case 7 details: guess kid starts in tracker " << " "
1051  << this_Kid->getVertex()[1] << " " << this_Kid->getVertex()[2] << " " << rdist << " "
1052  << has_pi0 << " " << std::endl;
1053 
1054  } else { // guess in calo -> mother is originator
1055  remap_as_you_go[mcp]=mother;
1056  mcp=mother; // mother started at ip, and reached the calo, and is the one
1057  // the hit should be attributed to.
1058  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
1059  << caloHit->id() << " attributed to mother " << mcp->id() <<
1060  " because it's origin is guessed in tracker : originator case 8 " <<std::endl;
1061  streamlog_out( DEBUG2 ) << " Case 8 details: guess kid starts in calo " << " "
1062  << this_Kid->getVertex()[1] << " " << this_Kid->getVertex()[2] << " " << rdist << " "
1063  << has_pi0 << " " << std::endl;
1064  }
1065 
1066  }
1067  }
1068  }
1069  }
1070 
1071  } else {
1072  // first case: the hit generating mcp itself already fulfills the firest
1073  // criterium, ie. it is a generator particle
1074  streamlog_out( DEBUG3 ) << " Hit " << simHit->id() << " / "
1075  << caloHit->id() << " attributed to " << mcp->id() <<
1076  " because it's origin is a generator particle : originator case 1 " <<std::endl;
1077 
1078  remap_as_you_go[mcp]=mcp;
1079  } // genstat if - then - else
1080  }
1081 
1082  streamlog_out( DEBUG4 ) <<" Final assignment for contribution " << k << " to " << simHit->id() << " / "
1083  << caloHit->id() << " : " << mcp->id() << std::endl;
1084  streamlog_out( DEBUG4 ) <<" gs "<<mcp->getGeneratorStatus()<<
1085  " dint "<<mcp->isDecayedInTracker()<<
1086  " bs "<<mcp->isBackscatter()<<
1087  " npar "<<mcp->getParents().size()<<
1088  " pdg "<<mcp->getPDG()<<
1089  " "<<mcp->getEndpoint()[0]<<
1090  " "<<mcp->getEndpoint()[1]<<
1091  " "<<mcp->getEndpoint()[2]<<std::endl;
1092 
1093  // decided which mcp this sim-hit should be associated to :
1094  // already_known:
1095 
1096  simHitMapEnergy[ mcp ] += e ;
1097  // chitTruthRelNav.addRelation( caloHit , mcp , e ) ;
1098  for ( size_t jk=0; jk<caloHits.size(); jk++) {
1099  chitTruthRelNav.addRelation( caloHits[jk] , mcp , e*( (CalorimeterHit*) caloHits[jk])->getEnergy()/allRecoEn ) ;
1100  }
1101 
1102  } // mc-contributon-to-simHit loop
1103  }
1104  }
1105  } // sim-calo-hits loop
1106 
1107  streamlog_out( DEBUG6 ) << " *** Sorting out simHit<->MCParticle connections : DONE " << std::endl;
1108  streamlog_out( DEBUG6 ) << " *** Sorting out Cluster<->MCParticle using simHit<->MCParticle, re-assigning the latter in some rare cases." << std::endl;
1109 
1110 
1111  // loop over reconstructed particles
1112  int nCluster = clusterCol->getNumberOfElements() ;
1113 
1114  streamlog_out( DEBUG6 ) << std::endl;
1115  streamlog_out( DEBUG6 ) << " ================= " << std::endl;
1116  streamlog_out( DEBUG6 ) <<" Treating clusters. There are "<< nCluster <<" of them " << std::endl;
1117 
1118  std::vector<Cluster*> missingMC ;
1119  missingMC.reserve( nCluster ) ;
1120  int ifoundclu =0;
1121  // now for the clusters
1122 
1123 
1124  MCParticle* mother = 0;
1125 
1126  MCPMapDouble mcpEnergyTot ;
1127 
1128 
1129  for(int i=0;i<nCluster;i++){
1130 
1131  MCPMapDouble mcpEnergy ;
1132  double eTot = 0 ;
1133  Cluster* clu = dynamic_cast<Cluster*> ( clusterCol->getElementAt(i) ) ;
1134 
1135 
1136 
1137 
1138  // We need to find all seen hits this clutser is made of, which sim hits each
1139  // of the seen hits came from, and finally which true particles actually created
1140  // each sim hit. Contrary to the sim tracker hits above, a sim-calo hit can be
1141  // made by several true particles. They also have a signal size (energy) value.
1142  // In addition, the true particle creating sometimes needs to be back-tracked
1143  // to the particle actually entering tha calorimeter.
1144 
1145 
1146 
1147 
1148  const CalorimeterHitVec& cluHits = clu->getCalorimeterHits() ;
1149 
1150  streamlog_out( DEBUG5 ) << std::endl;
1151  streamlog_out( DEBUG5 ) << " ================= " << std::endl;
1152  streamlog_out( DEBUG5 ) <<"Cluster clu = "<< clu->id() << " (i = " << i << " ) with " << cluHits.size() << " hits " << std::endl;
1153 
1154  double ecalohitsum=0.;
1155  double ecalohitsum_unknown=0.;
1156  int no_sim_hit = 0;
1157  for( CalorimeterHitVec::const_iterator hitIt = cluHits.begin() ;
1158  hitIt != cluHits.end() ; ++hitIt ) {
1159 
1160 
1161  CalorimeterHit* hit = * hitIt ; // ... a calo seen hit ...
1162  ecalohitsum+= hit->getEnergy();
1163 
1164 
1165  const LCObjectVec& simHits = *(this->getCaloHits(hit)) ;
1166 
1167  streamlog_out( DEBUG4 ) << std::endl;
1168  streamlog_out( DEBUG4 ) <<" Treating hit = "<< hit->id() << " e " << hit->getEnergy()<< " nb sim hits : " <<
1169  simHits.size() << std::endl;
1170 
1171  double ehit = 0.0;
1172  int nsimhit = 0;
1173  for( LCObjectVec::const_iterator objIt = simHits.begin() ;
1174  objIt != simHits.end() ; ++objIt ){
1175 
1176  SimCalorimeterHit* simHit = dynamic_cast<SimCalorimeterHit*>( *objIt ) ; // ... and a sim hit ....
1177 
1178  double calib_factor = hit->getEnergy()/simHit->getEnergy(); // In this place, this defn of calib_factor is still appropriate for split hits (DJ)
1179 
1180  streamlog_out( DEBUG3 ) <<" simhit = "<< simHit->id() << " has " << simHit->getNMCContributions()
1181  << " contributors " << std::endl;
1182 
1183  nsimhit++;
1184  for(int j=0;j<simHit->getNMCContributions() ;j++){
1185 
1186  MCParticle* mcp = simHit->getParticleCont( j ) ;
1187  double e = simHit->getEnergyCont( j ) * calib_factor ;
1188  streamlog_out( DEBUG3 ) <<" true contributor = "<< mcp->id() << " e: " << e
1189  <<" mapped to " <<remap_as_you_go.find(mcp)->second->id() << std::endl;
1190  if ( mcp == 0 ) {
1191  streamlog_out( DEBUG7 ) <<" simhit = "<< simHit << " has no creator " <<std::endl;
1192  // streamlog_out( DEBUG7 ) <<" true contributor = "<< mcp << " e: " << e
1193  // <<" mapped to " <<remap_as_you_go.find(mcp)->second << std::endl;
1194  continue ;
1195  }
1196  mcp=remap_as_you_go.find(mcp)->second;
1197  mcpEnergy[ mcp ] += e ;// count the hit-energy caused by this true particle
1198  eTot += e ; // total energy
1199  ehit+= e;
1200  } // mc-contributon-to-simHit loop
1201  } // simHit loop
1202 
1203  streamlog_out( DEBUG4 )<< " summed contributed e: " << ehit << " ratio : " << ehit/hit->getEnergy()
1204  << " nsimhit " << nsimhit <<std::endl;
1205  if ( nsimhit == 0 ) {
1206 
1207  streamlog_out( DEBUG5 ) << " Warning: no simhits for calohit " << hit <<
1208  ". Will have to guess true particle ... " << std::endl;
1209  no_sim_hit = 1;
1210  ecalohitsum_unknown+= hit->getEnergy();
1211  streamlog_out( DEBUG5 )<< " sim-less calohit E and position " << hit->getEnergy() << " "
1212  << hit->getPosition()[0] << " " << hit->getPosition()[1] << " " << hit->getPosition()[2] << " " << std::endl;
1213  streamlog_out( DEBUG5 )<< " sim-less calohit clust E and position " << clu->getEnergy()
1214  << " "<< clu->getPosition()[0] << " " << clu->getPosition()[1] << " " << clu->getPosition()[2] << " " << std::endl;
1215  }
1216  } // hit loop
1217  if ( no_sim_hit == 1 ) {
1218  streamlog_out( DEBUG6 ) << " Warning, there are sim-less calohits in cluster " << clu->id() << std::endl;
1219  }
1220  streamlog_out( DEBUG5 ) << std::endl;
1221  streamlog_out( DEBUG5 ) << " Sum of calohit E: " << ecalohitsum << " cluster E "
1222  << clu->getEnergy() << " Sum of Simcalohit E: " << eTot
1223  << " no_sim_hit: " << no_sim_hit << ". Energy from unknow source : "
1224  << ecalohitsum_unknown << std::endl;
1225  if( eTot == 0.0 ){ // fixme - this might happen if clusters are from Lcal/Muon only
1226 
1227 
1228  // save reco particle in missingMC
1229  missingMC.push_back( clu ) ;
1230 
1231  streamlog_out( DEBUG8 ) << " no calorimeter hits found for "
1232  << " cluster " << clu->id() << " e:" << clu->getEnergy()
1233  // << " charge: " << rec->getCharge()
1234  // << " px: " << rec->getMomentum()[0]
1235  // << " py: " << rec->getMomentum()[1]
1236  // << " pz: " << rec->getMomentum()[2]
1237  // << " mass: " << rec->getMass()
1238  // << " *pdg " << rec->getParticleIDUsed()
1239  << std::endl ;
1240 
1241 
1242  continue ;
1243  }
1244 
1245  // At this point, we have a list of all true particles contributiong to this cluster.
1246  // We now need to sum up the energies each particle contributes with and do
1247  // further parsing.
1248 
1249 
1250  // find the mc particle with the largest energy contribution
1251  // also store all genstat=1 particles and
1252  // all particles of type to be saved
1253 
1254  double eMax = 0 ; // energy the most contributing true particle added to
1255  // the cluster
1256 
1257  MCParticleVec theMCPs ; // vector that will contain all true particles contributing,
1258  // if they, ex officio, will be in the skimmed collection.
1259  theMCPs.reserve(1000);
1260  std::vector<double> MCPes; // energy contribution of these
1261  MCPes.reserve(1000);
1262  int ifound = 0; // total number of contributors
1263 
1264  MCParticleVec moreMCPs ; // vector that will contain true particles contributing, that
1265  // normally wouldn't be in the skimmed collection
1266  moreMCPs.reserve(1000);
1267  std::vector<double>moreMCPes; // energy contribution of these
1268  moreMCPes.reserve(1000);
1269  int morefound = 0; // number of such cases.
1270 
1271  mother = 0;
1272 
1273  for( MCPMapDouble::iterator it = mcpEnergy.begin() ; // iterate trough the map.
1274  it != mcpEnergy.end() ; ++it ){
1275  if ( it->first == 0 ) { // ( if == 0, this cluster contains some (but not all) sim-hits with unknown origin.
1276  // If *all* sim-hits would have had unknown origin, we would already have "continue":ed above)
1277  streamlog_out( MESSAGE ) << " SimHit with unknown origin in cluster " << clu << " ( " << clu->id() << " ) " << std::endl;
1278  continue;
1279  }
1280 
1281  if (it->first->getGeneratorStatus() == 1 ) { // genstat 1 particle, ie. it is a bona fide
1282  // creating true particle: enter it into
1283  // the list, and note how much energy it
1284  // contributed.
1285  theMCPs.push_back(it->first); MCPes.push_back(it->second); ifound++;
1286  } else { // not genstat 1. What should we do with it ?
1287  if ( it->first->getParents().size() != 0 ) {
1288  mother= dynamic_cast<MCParticle*>(it->first->getParents()[0]);
1289  } else {
1290  mother = 0 ;
1291  }
1292  if ( mother != 0 ) {
1293  if ( mother->isDecayedInTracker() && // ... so the partic apeared in the
1294  // tracker ...
1295  _pdgSet.find( abs (mother->getPDG())) != _pdgSet.end() ) {
1296  // ... and is of a type we want to
1297  // save (it's mother is in _pdgSet)
1298  // -> also a bona fide creator.
1299  theMCPs.push_back(it->first); MCPes.push_back(it->second); ifound++;
1300  } else { // else: if the mother is a BS, add to the moreMCPs-list, further treated below,
1301  // otherwise keep as a bona fide creator.
1302  streamlog_out( DEBUG2 ) << " case 1 for "<< it->first->id() <<
1303  "(morefound=" << morefound << ")" <<
1304  " mother: " << mother->id() <<
1305  " gs " <<mother->getGeneratorStatus()<<
1306  " dint " << mother->isDecayedInTracker() <<
1307  " bs " << mother->isBackscatter() <<
1308  " pdg " <<mother->getPDG() <<std::endl;
1309  if ( mother->isBackscatter() == 1 ) {
1310  moreMCPs.push_back(it->first); moreMCPes.push_back(it->second); morefound++;
1311  } else {
1312  theMCPs.push_back(it->first); MCPes.push_back(it->second); ifound++;
1313  }
1314  }
1315  } else { // not genstat 1, no mother ?! Also add to the moreMCPs-list.
1316  streamlog_out( DEBUG6 ) << " case 2 for "<< it->first->id() <<
1317  "(morefound=" << morefound << ")" << std::endl;
1318  moreMCPs.push_back(it->first); moreMCPes.push_back(it->second); morefound++;
1319  }
1320  }
1321 
1322 
1323  if( it->second > eMax ){
1324 
1325  eMax = it->second ;
1326  }
1327  }
1328 
1329  // We now have two lists of contributing true particles, one containing those that will be
1330  // put to the skimmed list in any case, one with those that will be put there only is they
1331  // produced a visible signal. In addition, we have the total contribution to the shower from
1332  // each true particle in separate lists.
1333 
1334  if ( morefound > 0 ) {
1335  for (int iii=0 ; iii<morefound ; iii++ ) {
1336  streamlog_out( DEBUG2 ) << " iii, moreMCPes[iii], moreMCPs[iii], gs " << iii <<
1337  " "<< moreMCPes[iii] <<
1338  " "<< moreMCPs[iii]->id()<<
1339  " "<<moreMCPs[iii]->getGeneratorStatus() << std::endl;
1340  }
1341  for (int iii=0 ; iii<ifound ; iii++ ) {
1342  streamlog_out( DEBUG2 ) << " iii, MCPes[iii], theMCPs[iii] " << iii <<
1343  " "<< MCPes[iii] <<" "<< theMCPs[iii]->id() << std::endl;
1344  }
1345  streamlog_out( DEBUG3 )<< " morefond: " << morefound <<std::endl;
1346  }
1347 
1348  // figure out what to do with the cases where true particles not automatically in the
1349  // skimmed list contributed to the cluster: Normally, the back-tracking above has
1350  // done the job, but sometimes one has a neutron or a photon playing pin-ball between
1351  // the calorimeter and the tracker, so that a back-scatter scatters back from the
1352  // tracker into the same cluster it came from. In that case, we should ignore the
1353  // back scatter (normally, of course, a back-scatterer is - if it reaches a
1354  // calorimeter - in a different cluster and should be kept as an originator).
1355  // I couldn't figure out a more efficient way of figuring this out than the
1356  // almost-always-do-nothing loop below ;(
1357 
1358  mother = 0;
1359  for (int iii=0 ; iii<morefound ; iii++ ) {
1360  if ( moreMCPs[iii]->getParents().size() != 0 ) {
1361  mother= dynamic_cast<MCParticle*>(moreMCPs[iii]->getParents()[0]);
1362  streamlog_out( DEBUG2 ) << " iii: " << iii << " mother: " << mother->id() <<std::endl;
1363  } else {
1364  mother = 0 ;
1365  }
1366  streamlog_out( DEBUG2 ) <<" partic vert: "<<moreMCPs[iii]->getVertex()[0]<<
1367  " "<<moreMCPs[iii]->getVertex()[1]<<
1368  " "<<moreMCPs[iii]->getVertex()[2]<<std::endl;
1369  streamlog_out( DEBUG2 ) <<" partic endpoint: "<<moreMCPs[iii]->getEndpoint()[0]<<
1370  " "<<moreMCPs[iii]->getEndpoint()[1]<<
1371  " "<<moreMCPs[iii]->getEndpoint()[2]<<std::endl;
1372  streamlog_out( DEBUG2 ) <<" partic status : gs "<<moreMCPs[iii]->getGeneratorStatus()<<
1373  " dint "<<moreMCPs[iii]->isDecayedInTracker ()<<
1374  " bs "<<moreMCPs[iii]->isBackscatter ()<<
1375  " npar "<<moreMCPs[iii]->getParents().size()<<
1376  " pdg " << moreMCPs[iii]->getPDG() <<
1377  " dinc "<<moreMCPs[iii]->isDecayedInCalorimeter ()<<
1378  " bye "<<moreMCPs[iii]->hasLeftDetector () <<
1379  " stop "<<moreMCPs[iii]->isStopped () << std::endl;
1380 
1381  while ( mother!= 0 && mother->getGeneratorStatus() !=2 ) { // back-track to the
1382  // beginning of the chain
1383 
1384  streamlog_out( DEBUG2 ) << " mother "<< mother->id() <<
1385  " gs " << mother->getGeneratorStatus() <<
1386  " dint " << mother->isDecayedInTracker()<< " " <<
1387  " bs " << mother->isBackscatter()<<
1388  " pdg " << mother->getPDG() <<
1389  " dinc " <<mother->isDecayedInCalorimeter()<< std::endl;
1390  streamlog_out( DEBUG2 ) <<" mother vert: "<<mother->getVertex()[0]<<
1391  " "<<mother->getVertex()[1]<<
1392  " "<<mother->getVertex()[2]<<std::endl;
1393  streamlog_out( DEBUG2 ) <<" mother endpoint: "<<mother->getEndpoint()[0]<<
1394  " "<<mother->getEndpoint()[1]<<
1395  " "<<mother->getEndpoint()[2]<<std::endl;
1396 
1397 
1398  // find out if this ancestor is a back-scatterer actually directly giving rise to hits in this cluster.
1399  // If so, we attribute the hits of moreMCPs[iii] to that true particle instead.
1400 
1401 
1402  for (int kkk=0 ; kkk<ifound ; kkk++){
1403  MCParticle* tmcp = theMCPs[kkk];
1404  if ( tmcp == mother ) {
1405  // find hits related to moreMCPs[iii]
1406  streamlog_out( DEBUG3 ) << " found " << moreMCPs[iii] <<
1407  "(iii= "<<iii <<")" << kkk <<
1408  " to be related to "<<mother->id() <<
1409  " add e : " << moreMCPes[iii] << std::endl;
1410  LCObjectVec hitvec = chitTruthRelNav.getRelatedFromObjects(moreMCPs[iii]);
1411  FloatVec evec=chitTruthRelNav.getRelatedFromWeights(moreMCPs[iii]);
1412  int one_mod_done = 0;
1413  for ( unsigned lll=0 ; lll<hitvec.size() ; lll++ ) {
1414  CalorimeterHit* hit = dynamic_cast<CalorimeterHit*>(hitvec[lll]);
1415  for( CalorimeterHitVec::const_iterator hitIt = cluHits.begin() ;
1416  hitIt != cluHits.end() ; ++hitIt ) {
1417  if ( *hitIt == hit ) { // ... a calo seen hit ...
1418  chitTruthRelNav.removeRelation( hit ,moreMCPs[iii] );
1419  chitTruthRelNav.addRelation( hit , mother , evec[lll] ) ;
1420  one_mod_done = 1;
1421  break ;
1422  }
1423  }
1424 
1425  }
1426  if ( one_mod_done == 1 ) {
1427  MCPes[kkk]+= moreMCPes[iii];
1428  goto endwhile;
1429  } else {
1430  streamlog_out( DEBUG3 ) << " However, no hits related to the current cluster found ??" << std::endl;
1431  }
1432 
1433  }
1434  }
1435  if ( mother->getParents().size() != 0 ) {
1436 
1437  mother= dynamic_cast<MCParticle*>(mother->getParents()[0]);
1438  } else { mother = 0; }
1439  }
1440  endwhile:
1441  if ( mother == 0 || mother->getGeneratorStatus() ==2 ) {
1442 
1443  // no other contributing true particle found among the ancestors
1444  // to moreMCPs[iii], so we add it to the
1445  // list of true particles to be saved.
1446 
1447  streamlog_out( DEBUG3 ) << " No relation found for "<< moreMCPs[iii]->id() <<
1448  ". Keep it as separate originator "<< std::endl;
1449  theMCPs.push_back(moreMCPs[iii]); MCPes.push_back(moreMCPes[iii]); ifound++;
1450  }
1451 
1452  }
1453  streamlog_out( DEBUG5 ) << " cluster " << clu->id() << " , E = " << clu->getEnergy() << " , ifound = " << ifound << std::endl;
1454 
1455  // finally calculate the weight of each true partic to the seen
1456  // (= energy_from_this_true/ total ), and add the weighted reltion.
1457 
1458  float totwgt =0.0;
1459  for (int iii=0 ; iii<ifound ; iii++ ) {
1460  float weight = (MCPes[iii]/eTot)*(clu->getEnergy()-ecalohitsum_unknown)/clu->getEnergy();
1461  mcpEnergyTot[ theMCPs[iii] ] += weight*clu->getEnergy();
1462 
1463  totwgt+= weight;
1464  if( theMCPs[iii] == 0 ) {
1465 
1466  streamlog_out( ERROR ) << " cluster " <<clu->id() << " has " << MCPes[iii] << " GeV of " << eTot
1467  << " GeV [ " << weight << " ] true energy "
1468  << " but no MCParticle " << std::endl ;
1469 
1470  continue ;
1471 
1472  }
1473 
1474  streamlog_out( DEBUG5 ) << " cluster " <<clu->id() << " has " << MCPes[iii] << " of " << eTot
1475  << " [ " << weight << " ] "
1476  << " of MCParticle " << theMCPs[iii]->id() << " with pdg : " << theMCPs[iii]->getPDG()
1477  << " and genstat : " << theMCPs[iii]->getGeneratorStatus()
1478  << std::endl ;
1479 
1480  //Particle gun particles dont have parents but have genStat 0
1481  if ( theMCPs[iii]->getGeneratorStatus() == 0 && ( theMCPs[iii]->getParents().empty() == false) ) {
1482  streamlog_out( DEBUG5 ) << " mother id " << theMCPs[iii]->getParents()[0]->id()
1483  << " and genstat "
1484  << theMCPs[iii]->getParents()[0]->getGeneratorStatus() << std::endl ;
1485 
1486  }
1487 
1488 
1489  clusterTruthRelNav.addRelation( clu , theMCPs[iii] , weight ) ;
1490  // The reverse map. Differs in what the weight means: here it means:
1491  // "this cluster got wgt of all seen cluster energy the true produced"
1492  // (in the other one it means:
1493  // "this true contributed wgt to the total seen energy of the cluster")
1494  weight=(MCPes[iii]/simHitMapEnergy[theMCPs[iii]]);
1495  truthClusterRelNav.addRelation( theMCPs[iii] , clu , weight ) ;
1496 
1497  }
1498  ifoundclu=ifound;
1499  } // cluster loop
1500 
1501  streamlog_out( DEBUG6 ) << " *** Sorting out Cluster<->MCParticle : DONE" << std::endl;
1502  streamlog_out( DEBUG6 ) << " *** Guessing Cluster<->MCParticle in cases where MCParticle<->simHit is broken (LCAL/DBD)" << std::endl;
1503 
1504  // recover missing MCParticles for neutrals, typically LCal. Not fixed in DBD ... :
1505  // attach the MCParticle with smallest angle to cluster
1506 
1507  // (This is from the original RecoMCTruthLinker. I didn't revise it/MB)
1508 
1509  for( unsigned i=0 ; i < missingMC.size() ; ++i ) {
1510 
1511  int nMCP = mcpCol->getNumberOfElements() ;
1512 
1513  Cluster* clu = missingMC[i] ;
1514 
1515 
1516  gear::Vector3D recP( clu->getPosition()[0] , clu->getPosition()[1] ,
1517  clu->getPosition()[2] ) ;
1518 
1519  double recTheta = recP.theta() ;
1520 
1521  double maxProd = 0.0 ;
1522  MCParticle* closestMCP = 0 ;
1523 
1524  for (int j=0; j< nMCP ; j++){
1525 
1526  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( j ) ) ;
1527 
1528  if ( fabs( mcp->getCharge() ) > 0.01 ) {
1529  continue ;
1530  }
1531 
1532  gear::Vector3D mcpP( mcp->getMomentum()[0] , mcp->getMomentum()[1] ,
1533  mcp->getMomentum()[2] ) ;
1534 
1535  if ( fabs( recTheta - mcpP.theta() ) > 0.3 ) {// fixme : proc param...
1536  continue ;
1537  }
1538 
1539 
1540  double prod = mcpP.unit().dot( recP.unit() ) ;
1541 
1542  if ( prod > maxProd ) {
1543  maxProd = prod ;
1544  closestMCP = mcp ;
1545  }
1546  }
1547  if ( maxProd > 0. ) {
1548 
1549  streamlog_out( DEBUG5 )
1550  << " neutral cluster particle recovered"
1551  << clu->getEnergy()
1552  << " maxProd: " << maxProd
1553  << " px: " << closestMCP->getMomentum()[0]
1554  << " py: " << closestMCP->getMomentum()[1]
1555  << " pz: " << closestMCP->getMomentum()[2]
1556  << std::endl ;
1557 
1558  clusterTruthRelNav.addRelation( clu , closestMCP , 1.0 ) ;
1559  truthClusterRelNav.addRelation( closestMCP , clu , 1.0 ) ;
1560  // Could try to also add relation calohit <-> MCPart from this
1561  // fixup. However, this makes more confusion that clarification,
1562  // so, for now at least, leave it out
1563  // const CalorimeterHitVec& cluHits = clu->getCalorimeterHits() ;
1564  //
1565  // for( CalorimeterHitVec::const_iterator hitIt = cluHits.begin() ;
1566  // hitIt != cluHits.end() ; ++hitIt ) {
1567  //
1568  // CalorimeterHit* hit = * hitIt ; // ... a calo seen hit ...
1569  // streamlog_out( DEBUG5 )
1570  // <<"recuperated hit = "<< hit << " e " << hit->getEnergy()<< std::endl;
1571  // chitTruthRelNav.addRelation( hit , closestMCP , hit->getEnergy() ) ;
1572  // }
1573  }
1574 
1575 
1576  }
1577 
1578  // seen-true relation complete. add the collection
1579 
1580  streamlog_out( DEBUG6 ) << " *** Cluster linking complete, create collection " << std::endl;
1581  *trclcol = truthClusterRelNav.createLCCollection() ;
1582  *ctrlcol = clusterTruthRelNav.createLCCollection() ;
1583  *chittrlcol = chitTruthRelNav.createLCCollection() ;
1584 }
1585 
1586 
1587 void RecoMCTruthLinker::particleLinker( LCCollection* mcpCol, LCCollection* particleCol,
1588  LCCollection* ttrlcol, LCCollection* ctrlcol,
1589  LCCollection* trtlcol, LCCollection* trclcol,
1590  LCCollection** ptrlcol, LCCollection** trplcol) {
1591 
1592  LCRelationNavigator particleTruthRelNav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::MCPARTICLE ) ;
1593  LCRelationNavigator truthParticleRelNav( LCIO::MCPARTICLE , LCIO::RECONSTRUCTEDPARTICLE ) ;
1594 
1595  LCRelationNavigator trackTruthRelNav = LCRelationNavigator( ttrlcol );
1596  LCRelationNavigator clusterTruthRelNav = LCRelationNavigator( ctrlcol );
1597  LCRelationNavigator truthTrackRelNav = LCRelationNavigator( trtlcol );
1598  LCRelationNavigator truthClusterRelNav = LCRelationNavigator( trclcol );
1599 
1600  LCObjectVec mcvec;
1601  int nPart = particleCol->getNumberOfElements() ;
1602 
1603 
1604  static FloatVec www ;
1605 
1606  for(int i=0;i<nPart;++i){
1607 
1608  std::map< MCParticle* , int > mcmap;
1609 
1610  MCParticle* mcmax =0;
1611  float maxwgt=0. ;
1612 
1613  ReconstructedParticle* part = dynamic_cast<ReconstructedParticle*> ( particleCol->getElementAt(i) ) ;
1614  TrackVec tracks=part->getTracks();
1615  int ntrk = part->getTracks().size();
1616  ClusterVec clusters=part->getClusters();
1617  int nclu = part->getClusters().size();
1618  streamlog_out( DEBUG3 ) << " ======== " << std::endl;
1619  streamlog_out( DEBUG3 ) << " Treating particle " << part->id() << " with index " << i
1620  << " it has " << ntrk << " tracks, and " << nclu << " clusters " <<std::endl;
1621 
1622  int nhit[100] ;
1623  int nhitT = 0;
1624  if ( ntrk > 1 ) {
1625 
1626  for (int j=0 ; j < ntrk ; j++ ) {
1627  nhit[j] = 0;
1628  for ( unsigned kkk=0 ;kkk<tracks[j]->getSubdetectorHitNumbers().size(); kkk++ ) {
1629  nhit[j]+= tracks[j]->getSubdetectorHitNumbers()[kkk];
1630  nhitT+= tracks[j]->getSubdetectorHitNumbers()[kkk];
1631  }
1632  streamlog_out( DEBUG2 ) << " Track " << tracks[j]->id() << " with index " << j <<" has " << nhit[j] << " hits " << std::endl;
1633  }
1634  streamlog_out( DEBUG2 ) << " Total : " << nhitT << " hits " << std::endl;
1635  } else {
1636  nhit[0]=1 ; nhitT=1;
1637  }
1638  for (int j=0 ; j < ntrk ; j++ ) {
1639 
1640  if ( tracks[j] != 0 ) {
1641  mcvec = trackTruthRelNav.getRelatedToObjects( tracks[j]);
1642  www = trackTruthRelNav.getRelatedToWeights( tracks[j]);
1643  int ntp= mcvec.size();
1644  streamlog_out( DEBUG3 ) << " Track " << tracks[j]->id() << " with index " << j
1645  << " has " << ntp << " true particles " << std::endl;
1646  if ( ntp > 0 ) {
1647  if ( mcvec[0] != 0 ) {
1648 
1649  for ( int kkk=0 ; kkk < ntp ; kkk++ ) {
1650  if ( !_FullRecoRelation && www[kkk]*(float(nhit[j])/float(nhitT)) > maxwgt ) {
1651  maxwgt= www[kkk]*(float(nhit[j])/float(nhitT)) ;
1652  mcmax=dynamic_cast<MCParticle*>(mcvec[kkk]);
1653  }
1654  mcmap[dynamic_cast<MCParticle*>(mcvec[kkk])] +=
1655  int(www[kkk]*1000.*(float(nhit[j])/float(nhitT))+0.5);
1656  streamlog_out( DEBUG2 ) << " Individual track weight to " <<mcvec[kkk]<< " is "
1657  << www[kkk] << ", scaled one is "
1658  << www[kkk]*(float(nhit[j])/float(nhitT))
1659  << " ( loop -index : " << kkk << ")"<< std::endl;
1660  }
1661  }
1662  }
1663  }
1664  }
1665  if ( _FullRecoRelation || ntrk == 0 ) {
1666  double eclu[100] ;
1667  double ecluT = 0.;
1668  if ( nclu > 1 ) {
1669 
1670  for (int j=0 ; j < nclu ; j++ ) {
1671  eclu[j] = clusters[j]->getEnergy();
1672  ecluT+=clusters[j]->getEnergy();
1673  streamlog_out( DEBUG2 ) << " Cluster " << clusters[j]->id() << " with index " << j <<" has energy " << eclu[j] << std::endl;
1674  }
1675  streamlog_out( DEBUG2 ) << " Total : " << ecluT << std::endl;
1676  } else {
1677  eclu[0]=1 ; ecluT=1;
1678  }
1679  for (int j=0 ; j < nclu ; j++ ) {
1680  if ( clusters[j] != 0 ) {
1681  mcvec = clusterTruthRelNav.getRelatedToObjects(clusters[j]);
1682  www = clusterTruthRelNav.getRelatedToWeights(clusters[j]);
1683  int ntp= mcvec.size();
1684  streamlog_out( DEBUG3 ) << " Cluster " << clusters[j]->id() << " with index " << j
1685  << " has " << ntp << " true particles " << std::endl;
1686  if ( ntp > 0 ) {
1687  if ( mcvec[0] != 0 ) {
1688  for ( int kkk=0 ; kkk < ntp ; kkk++ ) {
1689  if ( !_FullRecoRelation && www[kkk]*(eclu[j]/ecluT) > maxwgt ) {
1690  maxwgt= www[kkk]*(eclu[j]/ecluT) ;
1691  mcmax=dynamic_cast<MCParticle*>(mcvec[kkk]);
1692  }
1693  mcmap[dynamic_cast<MCParticle*>(mcvec[kkk])] +=
1694  int(www[kkk]*1000.*(eclu[j]/ecluT)+0.5)*10000;
1695  streamlog_out( DEBUG2 ) << " Individual cluster Weight to " <<mcvec[kkk]<< " is "
1696  << www[kkk] << ", scaled one is "
1697  << www[kkk]*(eclu[j]/ecluT)
1698  << " ( loop -index : " << kkk << ")"<< std::endl;
1699  }
1700  }
1701  }
1702  }
1703  }
1704  }
1705  if ( _FullRecoRelation ) {
1706  for ( std::map< MCParticle* , int >::iterator mcit = mcmap.begin() ;
1707  mcit != mcmap.end() ; mcit++ ) {
1708  // loop all MCparticles releted to the particle
1709  // get the true particle
1710  streamlog_out( DEBUG4 ) << " particle " << part->id() <<" has weight "<<mcit->second
1711  << " (Track: " << int(mcit->second)%10000
1712  << " , Cluster: " << int(mcit->second)/10000 << " ) "
1713  << " of MCParticle with pdg : " << mcit->first->getPDG()
1714  << " and genstat : " << mcit->first->getGeneratorStatus()
1715  << " id: " << mcit->first->id()
1716  << std::endl ;
1717 
1718  particleTruthRelNav.addRelation( part , mcit->first , mcit->second ) ;
1719  }
1720  } else {
1721  if( mcmax != NULL ) {
1722  //AS: FixMe: There is still something going wrong with particles from the particle gun.
1723  //Although there are perfect PFOs no link with the MCParticle is established.
1724  particleTruthRelNav.addRelation( part , mcmax , maxwgt ) ;
1725  streamlog_out( DEBUG3 ) << " particle " << part->id() << " has weight "<<maxwgt
1726  << " of MCParticle with pdg : " << mcmax->getPDG()
1727  << " and genstat : " << mcmax->getGeneratorStatus()
1728  << " id: " << mcmax->id()
1729  << ". Particle charge and ntracks : " << part->getCharge()<<" "<<ntrk
1730  << std::endl ;
1731  } else {
1732  streamlog_out( WARNING ) << " particle has weight "<< maxwgt
1733  << ". Particle charge and ntracks : " << part->getCharge()<<" "<<ntrk
1734  << " but no mcparticle found "
1735  << std::endl ;
1736  }
1737  }
1738  }
1739  // The reverse map. Differs in what the weight means: here it means:
1740  // "this cluster got wgt of all seen cluster energy the true produced"
1741  // (in the other one it means:
1742  // "this true contributed wgt to the total seen energy of the cluster")
1743  LCObjectVec partvec;
1744  static FloatVec www_clu ;
1745  static FloatVec www_trk ;
1746  LCObjectVec cluvec_t;
1747  LCObjectVec trkvec_t;
1748  ClusterVec cluvec_p;
1749  TrackVec trkvec_p;
1750  int nMCP = mcpCol->getNumberOfElements() ;
1751  for(int i=0;i< nMCP;i++){
1752  MCParticle* mcp =dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
1753  partvec = particleTruthRelNav.getRelatedFromObjects(mcp);
1754  cluvec_t = truthClusterRelNav.getRelatedToObjects(mcp);
1755  www_clu = truthClusterRelNav.getRelatedToWeights(mcp);
1756  trkvec_t = truthTrackRelNav.getRelatedToObjects(mcp);
1757  www_trk = truthTrackRelNav.getRelatedToWeights(mcp);
1758 
1759  for ( unsigned j=0 ; j<partvec.size() ; j++ ) {
1760  ReconstructedParticle* msp = dynamic_cast<ReconstructedParticle*>(partvec[j]) ;
1761  cluvec_p = msp->getClusters() ;
1762  trkvec_p = msp->getTracks() ;
1763  float c_wgt=0. ;
1764  for ( unsigned k=0 ; k<cluvec_p.size() ; k++ ) {
1765  for ( unsigned l=0 ; l<cluvec_t.size() ; l++ ) {
1766  if ( cluvec_p[k] == cluvec_t[l] ) {
1767  c_wgt+=www_clu[l];
1768  }
1769  }
1770  }
1771  float t_wgt=0. ;
1772  for ( unsigned k=0 ; k<trkvec_p.size() ; k++ ) {
1773  for ( unsigned l=0 ; l<trkvec_t.size() ; l++ ) {
1774  if ( trkvec_p[k] == trkvec_t[l] ) {
1775  t_wgt+=www_trk[l];
1776  }
1777  }
1778  }
1779  float wgt=int(c_wgt*1000)*10000 + int(t_wgt*1000) ;
1780  truthParticleRelNav.addRelation( mcp, msp , wgt ) ;
1781  streamlog_out( DEBUG4 ) << " True Particle " << mcp->id() << " ( pdg " << mcp->getPDG()
1782  << " ) has weight " << c_wgt <<" / " << t_wgt << "( " << int(wgt) << " ) to particle "
1783  << msp->id() << " with " << cluvec_p.size() << " clusters, and "
1784  << trkvec_p.size() << " tracks " << std::endl;
1785  }
1786  }
1787  streamlog_out( DEBUG6 ) << " particle linking complete, create collection " << std::endl;
1788  *ptrlcol = particleTruthRelNav.createLCCollection() ;
1789  *trplcol = truthParticleRelNav.createLCCollection() ;
1790 }
1791 void RecoMCTruthLinker::makeSkim( LCCollection* mcpCol , LCCollection* ttrlcol, LCCollection* ctrlcol , LCCollectionVec** skimVec){
1792 
1793 
1794  LCRelationNavigator trackTruthRelNav = LCRelationNavigator( ttrlcol );
1795  LCRelationNavigator clusterTruthRelNav = LCRelationNavigator( ctrlcol );
1796 
1797 
1798  //-------------- create skimmed MCParticle collection ------------------------
1799 
1800  // *skimVec = new LCCollectionVec( LCIO::MCPARTICLE ) ;
1801  (*skimVec)->setSubset( true) ; // flag as subset
1802 
1803 
1804  int nMCP = mcpCol->getNumberOfElements() ;
1805 
1806  for(int i=0; i< nMCP ; i++){
1807 
1808  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
1809 
1810 
1811  if( mcp->ext<MCPKeep>() == true ){
1812 
1813  continue ; // particle allready in skim
1814  }
1815 
1816 
1817 
1818 
1819  if ( ! mcp->isCreatedInSimulation() || mcp->getGeneratorStatus() != 0 ) {
1820 
1821  //FIXME: this is a workaround for a Mokka bug: the isCreatedInSimulation
1822  // is set also for generated particles....
1823 
1824  // keep all generated particles (complete event)
1825 
1826  mcp->ext<MCPKeep>() = true ;
1827 
1828  continue ;
1829 
1830  } else { // of those created in the simulation we keep those that actually are reconstructed
1831  // including all parents
1832 
1833  // truthRelNav is the one we created above, remember, so here we make sure that all
1834  // true particles related to seen ones (with the logic we used there) really will
1835  // be in the skimmed collection !
1836 
1837  const LCObjectVec& trackObjects = trackTruthRelNav.getRelatedFromObjects( (LCObject*) mcp ) ;
1838  const LCObjectVec& clusterObjects = clusterTruthRelNav.getRelatedFromObjects( (LCObject*) mcp ) ;
1839 
1840  if( trackObjects.size() > 0 || clusterObjects.size()>0){
1841 
1842  streamlog_out( DEBUG5 ) << " keep MCParticle - e :" << mcp->getEnergy()
1843  << " charge: " << mcp->getCharge()
1844  << " px: " << mcp->getMomentum()[0]
1845  << " py: " << mcp->getMomentum()[1]
1846  << " pz: " << mcp->getMomentum()[2]
1847  << " mass: " << mcp->getMass()
1848  << " pdg " << mcp->getPDG()
1849  << std::endl ;
1850 
1851  // keepMCParticles also flags all parents of a kept particle, guaranteeing that
1852  // the history of any contributor to a detected signal will be kept !
1853 
1854  keepMCParticle( mcp ) ;
1855  }
1856 
1857  } // else
1858 
1859  } // end mcp loop
1860  // --- loop again and add daughters of particles that are in the skimmed list and have a pdg in
1861  // the parameter vector 'KeepDaughtersPDG
1862 
1863  streamlog_out( DEBUG5 ) << " First loop done. Now search for KeepDaughtersPDG:s" << std::endl;
1864 
1865 
1866  if(_saveBremsstrahlungPhotons){
1867  for(int i=0; i< nMCP ; i++){
1868  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
1869  if( ( abs(mcp->getPDG()) == 22 ) && ( mcp->getEnergy() > _bremsstrahlungEnergyCut ) ){
1870  if( mcp->getParents().size() ){
1871  MCParticle* parent = mcp->getParents()[0];
1872  if( abs(parent->getPDG()) == 11){
1873  const float x = mcp->getVertex()[0];
1874  const float y = mcp->getVertex()[1];
1875  const float z = mcp->getVertex()[2];
1876  const float xpe = parent->getEndpoint()[0];
1877  const float ype = parent->getEndpoint()[1];
1878  const float zpe = parent->getEndpoint()[2];
1879  const float dx = x - xpe;
1880  const float dy = y - ype;
1881  const float dz = z - zpe;
1882  const float dr = sqrt(dx*dx+dy*dy+dz*dz);
1883  if(dr>100.){
1884  keepMCParticle( mcp ) ;
1885  }
1886  }
1887  }
1888  }
1889  }
1890  }
1891 
1892 
1893  for(int i=0; i< nMCP ; i++){
1894 
1895  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
1896 
1897  // keep the daughters of all decays in flight of particles in the pdg list (default: gamma, pi0, K0s)
1898  if( mcp->ext<MCPKeep>() == true && mcp->isDecayedInTracker() ){ //&& !mcp->isStopped() ){
1899 
1900  unsigned thePDG = abs( mcp->getPDG() ) ;
1901 
1902  if( _pdgSet.find( thePDG ) != _pdgSet.end() ) {
1903 
1904  const MCParticleVec& daughters = mcp->getDaughters() ;
1905 
1906  streamlog_out( DEBUG5 ) << " keeping daughters of particle with pdg : " << mcp->getPDG() << " : "
1907  << " [" << mcp->getGeneratorStatus() << "] :";
1908  // << " e :" << mcp->getEnergy()
1909  // << " isCreatedInSimulation :" << mcp->isCreatedInSimulation() << std::endl
1910  // << " isBackscatter :" << mcp->isBackscatter() << std::endl
1911  // << " vertexIsNotEndpointOfParent :" << mcp->vertexIsNotEndpointOfParent() << std::endl
1912  // << " isDecayedInTracker :" << mcp->isDecayedInTracker() << std::endl
1913  // << " isDecayedInCalorimeter :" << mcp->isDecayedInCalorimeter() << std::endl
1914  // << " hasLeftDetector :" << mcp->hasLeftDetector() << std::endl
1915  // << " isStopped :" << mcp->isStopped() << " : "
1916 
1917  streamlog_message( DEBUG5 ,
1918  if( mcp->getParents().size() ) ,
1919  " parent pdg : " << mcp->getParents()[0]->getPDG() << " : " ;
1920  ) ;
1921 
1922 
1923  // << std::endl ;
1924 
1925  for( MCParticleVec::const_iterator dIt = daughters.begin() ;
1926  dIt != daughters.end() ; ++dIt ){
1927 
1928 
1929  MCParticle* dau = dynamic_cast<MCParticle*>( *dIt ) ;
1930 
1931  if( dau->getEnergy()*1000. > _eCutMeV ) {
1932 
1933  (*dIt)->ext<MCPKeep>() = true ;
1934 
1935  streamlog_out( DEBUG5 ) << (*dIt)->getPDG() << ", " ;
1936  }
1937  }
1938 
1939  streamlog_out( DEBUG5 ) << std::endl ;
1940 
1941  }
1942  }
1943  }
1944 
1945  streamlog_out( DEBUG4 ) << " All found, add to skimmed list " << std::endl;
1946 
1947  for(int i=0; i< nMCP ; i++){
1948 
1949  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
1950 
1951  if( mcp->ext<MCPKeep>() == true ) {
1952 
1953  (*skimVec)->addElement( mcp ) ;
1954  }
1955  }
1956 
1957  // okidoki, the skimmed collection is complete. Add it.
1958 
1959 }
1960 
1961 void RecoMCTruthLinker::keepMCParticle( MCParticle* mcp ){
1962 
1963  mcp->ext<MCPKeep>() = true ;
1964 
1965 
1966  const MCParticleVec& parents = mcp->getParents() ;
1967 
1968  streamlog_out( DEBUG3 ) << " keepMCParticle keep particle with pdg : " << mcp->getPDG()
1969  << std::endl ;
1970 
1971  for( MCParticleVec::const_iterator pIt = parents.begin() ;
1972  pIt != parents.end() ; ++pIt ){
1973 
1974  if( (*pIt )->ext<MCPKeep>() != true ) { // if parent not yet in skim
1975 
1976  // add it
1977  keepMCParticle( *pIt ) ;
1978  }
1979 
1980  }
1981 }
1982 
1983 void RecoMCTruthLinker::check( LCEvent * evt ) {
1984 
1985  // ---- create some checkplots
1986 
1987 
1988 #ifdef MARLIN_USE_AIDA
1989 
1990  streamlog_out(DEBUG) << " check " << std::endl;
1991  // - define some static histo pointers
1992  // FIXME: these need to become class members eventually ...
1993 
1994  // static AIDA::IHistogram1D* hTrack_z0 ;
1995  static AIDA::ICloud1D* hmcp_etot ;
1996  static AIDA::ICloud1D* hmcp_e ;
1997  static AIDA::ICloud1D* hmcp_n ;
1998  static AIDA::ICloud1D* hmcp_ntot ;
1999 
2000  static AIDA::ICloud1D* hmcpsk_etot ;
2001  static AIDA::ICloud1D* hmcpsk_e ;
2002  static AIDA::ICloud1D* hmcpsk_n ;
2003  static AIDA::ICloud1D* hmcpsk_ntot ;
2004 
2005  if( isFirstEvent() ) {
2006 
2007  hmcp_e = AIDAProcessor::histogramFactory(this)->
2008  createCloud1D( "hmcp_e", " energy/GeV - all " , 100 ) ;
2009  // createHistogram1D( "hmcp_e", " energy/GeV ", 100, 0. , 10. ) ;
2010 
2011  hmcp_etot = AIDAProcessor::histogramFactory(this)->
2012  createCloud1D( "hmcp_etot", " total energy/GeV " , 100 ) ;
2013  // createHistogram1D( "hmcp_etot_e", " energy/GeV ", 1000, 0. , 1000. ) ;
2014 
2015  hmcp_n = AIDAProcessor::histogramFactory(this)->
2016  createCloud1D( "hmcp_n", " # generated stable particles " , 100 ) ;
2017 
2018  hmcp_ntot = AIDAProcessor::histogramFactory(this)->
2019  createCloud1D( "hmcp_ntot", " total # particles " , 100 ) ;
2020 
2021 
2022  hmcpsk_e = AIDAProcessor::histogramFactory(this)->
2023  createCloud1D( "hmcpsk_e", " energy/GeV - all " , 100 ) ;
2024  // createHistogram1D( "hmcpsk_e", " energy/GeV ", 100, 0. , 10. ) ;
2025 
2026  hmcpsk_etot = AIDAProcessor::histogramFactory(this)->
2027  createCloud1D( "hmcpsk_etot", " total energy/GeV " , 100 ) ;
2028  // createHistogram1D( "hmcpsk_etot_e", " energy/GeV ", 1000, 0. , 1000. ) ;
2029 
2030  hmcpsk_n = AIDAProcessor::histogramFactory(this)->
2031  createCloud1D( "hmcpsk_n", " # generated stable particles " , 100 ) ;
2032 
2033  hmcpsk_ntot = AIDAProcessor::histogramFactory(this)->
2034  createCloud1D( "hmcpsk_ntot", " total # particles " , 100 ) ;
2035  }
2036 
2037 
2038  LCCollection* mcpCol = NULL;
2039  try{
2040  mcpCol = evt->getCollection( _mcParticleCollectionName ) ;
2041  } catch (DataNotAvailableException& e) {
2042  streamlog_out(DEBUG9) << "RecoMCTructh::Check(): MCParticle collection \"" << _mcParticleCollectionName << "\" does not exist, skipping" << std::endl;
2043  }
2044 
2045  int nMCP = mcpCol->getNumberOfElements() ;
2046 
2047  double etot = 0.0 ;
2048  int nStable = 0 ;
2049 
2050  for(int i=0; i< nMCP ; i++){
2051 
2052  MCParticle* mcp = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
2053 
2054  if( mcp->getGeneratorStatus() == 1 ) {
2055 
2056  hmcp_e->fill( mcp->getEnergy() ) ;
2057 
2058  etot += mcp->getEnergy() ;
2059 
2060  ++nStable ;
2061  }
2062 
2063  }
2064  hmcp_n->fill( nStable ) ;
2065 
2066  hmcp_ntot->fill( nMCP ) ;
2067 
2068  hmcp_etot->fill( etot ) ;
2069 
2070 
2071  // create the same histos now with the skimmed collection
2072  LCCollection* mcpskCol = NULL;
2073  try{
2074  mcpskCol = evt->getCollection( _mcParticlesSkimmedName );
2075  }
2076  catch (DataNotAvailableException& e){
2077  streamlog_out(DEBUG9) << "RecoMCTructh::Check(): MCParticleSkimmed collection \"" << _mcParticlesSkimmedName << "\" does not exist, skipping" << std::endl;
2078  }
2079 
2080  etot = 0.0 ;
2081  nStable = 0 ;
2082  int nMCPSK = 0;
2083 
2084  if (mcpskCol){
2085  nMCPSK = mcpskCol->getNumberOfElements() ;
2086  for(int i=0; i< nMCPSK ; i++){
2087  MCParticle* mcpsk = dynamic_cast<MCParticle*> ( mcpskCol->getElementAt( i ) ) ;
2088  if( mcpsk->getGeneratorStatus() == 1 ) {
2089  hmcpsk_e->fill( mcpsk->getEnergy() ) ;
2090  etot += mcpsk->getEnergy() ;
2091  ++nStable ;
2092  }
2093  }
2094  }//If Skimmed Collection Exists
2095  //Fill this, even if MCParticles Skimmed is empty, all are 0
2096  hmcpsk_n->fill( nStable ) ;
2097 
2098  hmcpsk_ntot->fill( nMCPSK ) ;
2099 
2100  hmcpsk_etot->fill( etot ) ;
2101 
2102 #endif
2103 
2104 }
2105 
2106 const LCObjectVec* RecoMCTruthLinker::getSimHits( TrackerHit* trkhit, const FloatVec* weights ){
2107 
2108  const LCObjectVec* obj = & _navMergedTrackerHitRel->getRelatedToObjects(trkhit);
2109 
2110  if( obj->empty() == false ) {
2111 
2112  if(weights != 0 ) weights = & _navMergedTrackerHitRel->getRelatedToWeights(trkhit);
2113 
2114  }
2115  else {
2116  streamlog_out( WARNING ) << "getSimHits : TrackerHit : " << trkhit << " has no sim hits related. CellID0 = " <<
2117  trkhit->getCellID0() << " pos = " << trkhit->getPosition()[0] << " " << trkhit->getPosition()[1] << " " <<
2118  trkhit->getPosition()[2] << std::endl ;
2119  }
2120 
2121  return obj;
2122 
2123 }
2124 
2125 
2126 void RecoMCTruthLinker::end(){
2127 
2128  streamlog_out(DEBUG6) << " processed " << _nEvt << " events in " << _nRun << " runs "
2129  << std::endl ;
2130 
2131 }
2132 
2133 
2134 /** helper function to get collection safely */
2135 inline lcio::LCCollection* getCollection(lcio::LCEvent* evt, const std::string name ){
2136 
2137  if( name.size() == 0 )
2138  return 0 ;
2139 
2140  try{
2141 
2142  return evt->getCollection( name ) ;
2143 
2144  } catch( lcio::DataNotAvailableException& e ){
2145 
2146  streamlog_out( DEBUG2 ) << "getCollection : DataNotAvailableException : " << name << std::endl ;
2147 
2148  return 0 ;
2149  }
2150 }
2151 
2152 
2153 void RecoMCTruthLinker::mergeTrackerHitRelations(LCEvent * evt){
2154 
2155  unsigned nCol = _colNamesTrackerHitRelations.size() ;
2156 
2157  //--- copy existing collections to a vector first
2158  std::vector<LCCollection*> colVec ;
2159 
2160  for( unsigned i=0; i < nCol ; ++i) {
2161 
2162  LCCollection* col = getCollection ( evt , _colNamesTrackerHitRelations[i] ) ;
2163 
2164  if( col != 0 ){
2165 
2166  colVec.push_back( col ) ;
2167  } else {
2168 
2169  streamlog_out(DEBUG2) << " mergeTrackerHitRelations: input collection missing : " << _colNamesTrackerHitRelations[i] << std::endl ;
2170  }
2171  }
2172 
2173 
2174 
2175  streamlog_out( DEBUG2 ) << " mergeTrackerHitRelations: copied collection parameters ... " << std::endl ;
2176 
2177 
2178 
2179  nCol = colVec.size() ;
2180 
2181  for( unsigned i=0; i < nCol ; ++i) {
2182 
2183  LCCollection* col = colVec[i] ;
2184 
2185  if( i == 0 ){
2186  // copy collection flags and collection parameters from first collections
2187 
2188  _mergedTrackerHitRelCol = new LCCollectionVec( col->getTypeName() ) ;
2189 
2190  _mergedTrackerHitRelCol->setFlag( col->getFlag() ) ;
2191  _mergedTrackerHitRelCol->setSubset( true );
2192 
2193  StringVec stringKeys ;
2194  col->getParameters().getStringKeys( stringKeys ) ;
2195  for(unsigned ii=0; ii< stringKeys.size() ; ii++ ){
2196  StringVec vals ;
2197  col->getParameters().getStringVals( stringKeys[ii] , vals ) ;
2198  _mergedTrackerHitRelCol->parameters().setValues( stringKeys[ii] , vals ) ;
2199  }
2200  StringVec intKeys ;
2201  col->getParameters().getIntKeys( intKeys ) ;
2202  for(unsigned ii=0; ii< intKeys.size() ; ii++ ){
2203  IntVec vals ;
2204  col->getParameters().getIntVals( intKeys[ii] , vals ) ;
2205  _mergedTrackerHitRelCol->parameters().setValues( intKeys[ii] , vals ) ;
2206  }
2207  StringVec floatKeys ;
2208  col->getParameters().getFloatKeys( floatKeys ) ;
2209  for(unsigned ii=0; ii< floatKeys.size() ; ii++ ){
2210  FloatVec vals ;
2211  col->getParameters().getFloatVals( floatKeys[ii] , vals ) ;
2212  _mergedTrackerHitRelCol->parameters().setValues( floatKeys[ii] , vals ) ;
2213  }
2214 
2215 
2216 
2217  }
2218 
2219  int nEle = col->getNumberOfElements() ;
2220 
2221  for(int j=0; j < nEle ; ++j){
2222 
2223  _mergedTrackerHitRelCol->addElement( col->getElementAt(j) ) ;
2224 
2225  }
2226 
2227  }
2228 
2229  if( nCol != 0 ) _navMergedTrackerHitRel = new LCRelationNavigator( _mergedTrackerHitRelCol );
2230 
2231 }
2232 
2233 
2234 
2235 
2236 
2237 
2238 
2239 void RecoMCTruthLinker::mergeCaloHitRelations(LCEvent * evt){
2240 
2241  unsigned nCol = _caloHitRelationNames.size() ;
2242 
2243  //--- copy existing collections to a vector first
2244  std::vector<LCCollection*> colVec ;
2245 
2246  for( unsigned i=0; i < nCol ; ++i) {
2247 
2248  LCCollection* col = getCollection ( evt , _caloHitRelationNames[i] ) ;
2249 
2250  if( col != 0 ){
2251 
2252  colVec.push_back( col ) ;
2253 
2254  } else {
2255 
2256  streamlog_out(DEBUG2) << " mergeCaloHitRelations: input collection missing : " << _caloHitRelationNames[i] << std::endl ;
2257  }
2258  }
2259 
2260 
2261 
2262  streamlog_out( DEBUG2 ) << " mergeCaloHitRelations: copied collection parameters ... " << std::endl ;
2263 
2264 
2265 
2266  nCol = colVec.size() ;
2267 
2268  for( unsigned i=0; i < nCol ; ++i) {
2269 
2270  LCCollection* col = colVec[i] ;
2271 
2272  if( i == 0 ){
2273  // copy collection flags and collection parameters from first collections
2274 
2275  _mergedCaloHitRelCol = new LCCollectionVec( col->getTypeName() ) ;
2276  _mergedCaloHitRelCol->setFlag( col->getFlag() ) ;
2277  _mergedCaloHitRelCol->setSubset( true );
2278 
2279  StringVec stringKeys ;
2280  col->getParameters().getStringKeys( stringKeys ) ;
2281  for(unsigned ii=0; ii< stringKeys.size() ; ii++ ){
2282  StringVec vals ;
2283  col->getParameters().getStringVals( stringKeys[ii] , vals ) ;
2284  _mergedCaloHitRelCol->parameters().setValues( stringKeys[ii] , vals ) ;
2285  }
2286  StringVec intKeys ;
2287  col->getParameters().getIntKeys( intKeys ) ;
2288  for(unsigned ii=0; ii< intKeys.size() ; ii++ ){
2289  IntVec vals ;
2290  col->getParameters().getIntVals( intKeys[ii] , vals ) ;
2291  _mergedCaloHitRelCol->parameters().setValues( intKeys[ii] , vals ) ;
2292  }
2293  StringVec floatKeys ;
2294  col->getParameters().getFloatKeys( floatKeys ) ;
2295  for(unsigned ii=0; ii< floatKeys.size() ; ii++ ){
2296  FloatVec vals ;
2297  col->getParameters().getFloatVals( floatKeys[ii] , vals ) ;
2298  _mergedCaloHitRelCol->parameters().setValues( floatKeys[ii] , vals ) ;
2299  }
2300 
2301 
2302 
2303  }
2304 
2305  int nEle = col->getNumberOfElements() ;
2306 
2307  for(int j=0; j < nEle ; ++j){
2308 
2309  _mergedCaloHitRelCol->addElement( col->getElementAt(j) ) ;
2310 
2311  }
2312 
2313  }
2314 
2315  if( nCol != 0 ) _navMergedCaloHitRel = new LCRelationNavigator( _mergedCaloHitRelCol );
2316 
2317 }
2318 
2319 const LCObjectVec* RecoMCTruthLinker::getCaloHits( CalorimeterHit* calohit, const FloatVec* weights ){
2320 
2321  const LCObjectVec* obj = & _navMergedCaloHitRel->getRelatedToObjects(calohit);
2322 
2323  if( obj->empty() == false ) {
2324 
2325  if(weights != 0 ) weights = & _navMergedCaloHitRel->getRelatedToWeights(calohit);
2326 
2327  }
2328  else {
2329 
2330  streamlog_out( DEBUG5 ) << "getCaloHits : CalorimeterHit : " << calohit << " has no sim hits related. CellID0 = " <<
2331  calohit->getCellID0() << " pos = " << calohit->getPosition()[0] << " " << calohit->getPosition()[1] << " " <<
2332  calohit->getPosition()[2] << std::endl ;
2333  }
2334  return obj;
2335 
2336 }
2337 
2338 
2339 
2340 void RecoMCTruthLinker::linkPrinter ( LCCollection* mcpCol, LCCollection* particleCol, LCCollection* ptrlcol, LCCollection* trplcol) {
2341 
2342  LCRelationNavigator particleTruthRelNav = LCRelationNavigator( ptrlcol );
2343  LCRelationNavigator truthParticleRelNav = LCRelationNavigator( trplcol );
2344 
2345  if ( ! _FullRecoRelation ) {
2346  // Not useful in this case ...
2347  return;
2348  }
2349  static FloatVec www ;
2350  static FloatVec www_other_way ;
2351 
2352  LCObjectVec mcvec;
2353  int nPart = particleCol->getNumberOfElements() ;
2354  LCObjectVec partvec;
2355 
2356  for(int i=0;i<nPart;++i){
2357 
2358  ReconstructedParticle* part = dynamic_cast<ReconstructedParticle*> ( particleCol->getElementAt(i) ) ;
2359  double clu_e=0.0;
2360  if ( part->getClusters().size() > 0 ) {
2361  Cluster* clu = part->getClusters()[0];
2362  if ( clu != 0 ) {
2363  clu_e=clu->getEnergy();
2364  }
2365  }
2366  mcvec = particleTruthRelNav.getRelatedToObjects( part);
2367  www = particleTruthRelNav.getRelatedToWeights( part);
2368  int ntp= mcvec.size();
2369  streamlog_out( DEBUG5 ) << " Particle " << part->id() << " (q: " << int(part->getCharge()) <<
2370  " ) has " << ntp << " true particles. E= "<< part->getEnergy() << " " << clu_e ;
2371  streamlog_out( DEBUG3 )<< " Index, id, PDG and energy of contributors: " << std::endl;
2372  if ( ntp > 0 ) {
2373  if ( mcvec[0] != 0 ) {
2374  double total_trk_weight = 0.0;
2375  double total_clu_weight = 0.0;
2376  double total_e_from_neutrals = 0.0;
2377  double true_charged_E=0.0;
2378  double total_ecalo_in_this_true=0.;
2379  double total_ecalo_neutral_in_charged=0.;
2380  for ( int kkk=0 ; kkk < ntp ; kkk++ ) {
2381  total_trk_weight+=(int(www[kkk])%10000)/1000.0;
2382  total_clu_weight+=(int(www[kkk])/10000)/1000.0;
2383  MCParticle* mcp = dynamic_cast< MCParticle*>( mcvec[kkk]);
2384  streamlog_out( DEBUG3 ) <<" "<< kkk << " " << mcp->id() << " " << mcp->getPDG() << " " << mcp->getEnergy() << std::endl ;
2385  if ( part->getCharge() != 0 ) {
2386  if ( mcp->getCharge() != 0 && (int(www[kkk])%10000)/1000.0 > 0.6 ) {
2387  true_charged_E=mcp->getEnergy();
2388  double ecalo_from_this_true=((int(www[kkk])/10000)/1000.0)*clu_e;
2389  if ( ecalo_from_this_true != 0.0 ) {
2390  partvec = truthParticleRelNav.getRelatedToObjects( mcp);
2391  www_other_way = truthParticleRelNav.getRelatedToWeights( mcp);
2392  for ( unsigned jjj=0 ; jjj < partvec.size() ; jjj++ ) {
2393  if ( partvec[jjj] == part ) {
2394  total_ecalo_in_this_true = ecalo_from_this_true/(int(www_other_way[jjj]/10000)/1000.0) ;
2395  break;
2396  }
2397  }
2398  }
2399  } else {
2400  total_ecalo_neutral_in_charged+=((int(www[kkk])/10000)/1000.0)*clu_e;
2401  }
2402  } else {
2403  if ( mcp->getCharge() == 0 ) {
2404  total_e_from_neutrals+=((int(www[kkk])/10000)/1000.0)*clu_e;
2405  }
2406  }
2407  }
2408  if ( part->getCharge() != 0 ) {
2409  streamlog_out( DEBUG5 ) << " True E_ch " << true_charged_E << " total E_calo from true " << total_ecalo_in_this_true
2410  << " E_calo from neutral " << total_ecalo_neutral_in_charged << " track weight " <<
2411  total_trk_weight << " cluster weight " << total_clu_weight << std::endl ;
2412  } else {
2413  streamlog_out( DEBUG5 ) << " True E_neutral " << total_e_from_neutrals << " track weight " <<
2414  total_trk_weight << " cluster weight " << total_clu_weight << std::endl ;
2415  }
2416  }
2417  }
2418  }
2419 
2420  int nMCPart = mcpCol->getNumberOfElements() ;
2421 
2422  for(int i=0;i<nMCPart;++i){
2423 
2424  MCParticle* mcpart = dynamic_cast<MCParticle*> ( mcpCol->getElementAt(i) ) ;
2425 
2426  partvec = truthParticleRelNav.getRelatedToObjects( mcpart);
2427  www = truthParticleRelNav.getRelatedToWeights( mcpart);
2428  int nsp= partvec.size();
2429  if ( nsp > 0 ) {
2430  streamlog_out( DEBUG5 ) << " TrueParticle " << mcpart->id() << " (q: " << int(mcpart->getCharge()) <<
2431  " ) has " << nsp << " seen particles. E= "<<mcpart->getEnergy() << " pt= " <<
2432  sqrt(mcpart->getMomentum()[0]* mcpart->getMomentum()[0]+ mcpart->getMomentum()[1]* mcpart->getMomentum()[1]);
2433  streamlog_out( DEBUG3 )<< " index, id, and enregy of reco particle: " << std::endl;
2434  if ( partvec[0] != 0 ) {
2435  double total_trk_weight = 0.0;
2436  double total_clu_weight = 0.0;
2437  for ( int kkk=0 ; kkk < nsp ; kkk++ ) {
2438  total_trk_weight+=(int(www[kkk])%10000)/1000.0;
2439  total_clu_weight+=(int(www[kkk])/10000)/1000.0;
2440  ReconstructedParticle* recopart = dynamic_cast< ReconstructedParticle*>(partvec[kkk]);
2441  streamlog_out( DEBUG3 )<<" "<< kkk << " " << recopart->id() << " " << recopart->getEnergy() << std::endl ;
2442  }
2443  streamlog_out( DEBUG5 ) << " Total track weight " << total_trk_weight << " , Total cluster weight " << total_clu_weight << std::endl ;
2444  }
2445  }
2446  }
2447 }
2448 
2449 
2450 
RecoMCTruthLinker::linkPrinter
virtual void linkPrinter(LCCollection *mcpCol, LCCollection *particleCol, LCCollection *ptrlcol, LCCollection *trplcol)
Definition: RecoMCTruthLinker.cc:2340
RecoMCTruthLinker::_trackMCTruthLinkName
std::string _trackMCTruthLinkName
output collection names
Definition: RecoMCTruthLinker.h:190
RecoMCTruthLinker::_calohitMCTruthLinkName
std::string _calohitMCTruthLinkName
Definition: RecoMCTruthLinker.h:197
RecoMCTruthLinker::end
virtual void end()
Called after data processing for clean up.
Definition: RecoMCTruthLinker.cc:2126
RecoMCTruthLinker::_mergedCaloHitRelCol
LCCollectionVec * _mergedCaloHitRelCol
Definition: RecoMCTruthLinker.h:184
RecoMCTruthLinker::clusterLinker
virtual void clusterLinker(LCEvent *evt, LCCollection *mcpCol, LCCollection *clusterCol, LCCollection **ctrcol, LCCollection **trccol, LCCollection **chittrlcol)
Definition: RecoMCTruthLinker.cc:639
RecoMCTruthLinker::_navMergedTrackerHitRel
LCRelationNavigator * _navMergedTrackerHitRel
Definition: RecoMCTruthLinker.h:183
RecoMCTruthLinker::trackLinker
virtual void trackLinker(LCEvent *evt, LCCollection *mcpCol, LCCollection *trackCol, LCCollection **ttrcol, LCCollection **trtcol)
Definition: RecoMCTruthLinker.cc:442
TrackMap
std::map< Track *, int > TrackMap
Definition: RecoMCTruthLinker.cc:45
RecoMCTruthLinker::_colNamesTrackerHitRelations
StringVec _colNamesTrackerHitRelations
Definition: RecoMCTruthLinker.h:179
Vector3D
CLHEP::Hep3Vector Vector3D
Definition: DistilledPFOCreator.cc:11
RecoMCTruthLinker::processEvent
virtual void processEvent(LCEvent *evt)
Called for every event - the working horse.
Definition: RecoMCTruthLinker.cc:325
sistrip::k
static const float k
Definition: PhysicalConstants.h:73
RecoMCTruthLinker::_trackCollectionName
std::string _trackCollectionName
Definition: RecoMCTruthLinker.h:172
RecoMCTruthLinker::mergeTrackerHitRelations
virtual void mergeTrackerHitRelations(LCEvent *evt)
Definition: RecoMCTruthLinker.cc:2153
RecoMCTruthLinker::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Called for every run.
Definition: RecoMCTruthLinker.cc:319
RecoMCTruthLinker::_FullRecoRelation
bool _FullRecoRelation
output collection steering
Definition: RecoMCTruthLinker.h:199
RecoMCTruthLinker::getSimHits
const LCObjectVec * getSimHits(TrackerHit *trkhit, const FloatVec *weights=NULL)
Definition: RecoMCTruthLinker.cc:2106
RecoMCTruthLinker::_navMergedCaloHitRel
LCRelationNavigator * _navMergedCaloHitRel
Definition: RecoMCTruthLinker.h:185
RecoMCTruthLinker::init
virtual void init()
Called at the begin of the job before anything is read.
Definition: RecoMCTruthLinker.cc:284
MCPMap
std::map< MCParticle *, int > MCPMap
Definition: RecoMCTruthLinker.cc:43
RecoMCTruthLinker::_mergedTrackerHitRelCol
LCCollectionVec * _mergedTrackerHitRelCol
Definition: RecoMCTruthLinker.h:182
RecoMCTruthLinker::makeSkim
virtual void makeSkim(LCCollection *mcpCol, LCCollection *ttrcol, LCCollection *ctrcol, LCCollectionVec **skimVec)
Definition: RecoMCTruthLinker.cc:1791
RecoMCTruthLinker::_clusterMCTruthLinkName
std::string _clusterMCTruthLinkName
Definition: RecoMCTruthLinker.h:192
Remap_as_you_go
std::map< MCParticle *, MCParticle * > Remap_as_you_go
Definition: RecoMCTruthLinker.cc:49
RecoMCTruthLinker::_invertedNonDestructiveInteractionLogic
bool _invertedNonDestructiveInteractionLogic
Definition: RecoMCTruthLinker.h:209
RecoMCTruthLinker::_simTrkHitCollectionNames
StringVec _simTrkHitCollectionNames
Definition: RecoMCTruthLinker.h:176
RecoMCTruthLinker::_recoMCTruthLinkName
std::string _recoMCTruthLinkName
Definition: RecoMCTruthLinker.h:194
RecoMCTruthLinker::check
virtual void check(LCEvent *evt)
Definition: RecoMCTruthLinker.cc:1983
RecoMCTruthLinker::_mCTruthRecoLinkName
std::string _mCTruthRecoLinkName
Definition: RecoMCTruthLinker.h:195
RecoMCTruthLinker::_mcParticlesSkimmedName
std::string _mcParticlesSkimmedName
Definition: RecoMCTruthLinker.h:196
RecoMCTruthLinker::keepMCParticle
void keepMCParticle(MCParticle *mcp)
Definition: RecoMCTruthLinker.cc:1961
RecoMCTruthLinker::_recoParticleCollectionName
std::string _recoParticleCollectionName
Definition: RecoMCTruthLinker.h:174
RecoMCTruthLinker::mergeCaloHitRelations
virtual void mergeCaloHitRelations(LCEvent *evt)
Definition: RecoMCTruthLinker.cc:2239
getCollection
lcio::LCCollection * getCollection(lcio::LCEvent *evt, const std::string name)
helper function to get collection safely
Definition: RecoMCTruthLinker.cc:2135
RecoMCTruthLinker::_mcParticleCollectionName
std::string _mcParticleCollectionName
input collection names
Definition: RecoMCTruthLinker.h:171
RecoMCTruthLinker::particleLinker
virtual void particleLinker(LCCollection *mcpCol, LCCollection *particleCol, LCCollection *ttrcol, LCCollection *ctrlcol, LCCollection *trtlcol, LCCollection *trclcol, LCCollection **ptrlcol, LCCollection **trplcol)
Definition: RecoMCTruthLinker.cc:1587
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
RecoMCTruthLinker::getCaloHits
const LCObjectVec * getCaloHits(CalorimeterHit *calohit, const FloatVec *weights=NULL)
Definition: RecoMCTruthLinker.cc:2319
RecoMCTruthLinker::_caloHitRelationNames
StringVec _caloHitRelationNames
Definition: RecoMCTruthLinker.h:180
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
RecoMCTruthLinker
Optionally creates four collections of LCRelations (&quot;recoMCTruthLink&quot;, &quot;trackMCTruthLink&quot;, &quot;mcTruthTrackLink&quot;, &quot;clusterMCTruthLink&quot;, &quot;&quot;clusterMCTruthLink&quot; and &quot;calohitMCTruthLink") with weighetd relations between true particles and reconstructed particles, tracks, clusters, and calorimeter hits, respectively.
Definition: RecoMCTruthLinker.h:106
RecoMCTruthLinker::_mCTruthTrackLinkName
std::string _mCTruthTrackLinkName
Definition: RecoMCTruthLinker.h:191
RecoMCTruthLinker::_mCTruthClusterLinkName
std::string _mCTruthClusterLinkName
Definition: RecoMCTruthLinker.h:193
aRecoMCTruthLinker
RecoMCTruthLinker aRecoMCTruthLinker
Definition: RecoMCTruthLinker.cc:38
RecoMCTruthLinker::_simCaloHitCollectionNames
StringVec _simCaloHitCollectionNames
Definition: RecoMCTruthLinker.h:177
RecoMCTruthLinker::_clusterCollectionName
std::string _clusterCollectionName
Definition: RecoMCTruthLinker.h:173
sistrip::StringVec
std::vector< std::string > StringVec
Definition: SiStripClus.h:56
MCPMapDouble
std::map< MCParticle *, double > MCPMapDouble
Definition: RecoMCTruthLinker.cc:47