ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTruthSiliconAssociation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTruthSiliconAssociation.cc
2 
4 
6 #include <phool/getClass.h>
7 #include <phool/phool.h>
8 
12 #include <trackbase/TrkrHitSet.h>
14 #include <trackbase/TrkrDefs.h>
21 
22 #include <g4main/PHG4Hit.h> // for PHG4Hit
23 #include <g4main/PHG4Particle.h> // for PHG4Particle
25 #include <g4main/PHG4HitDefs.h> // for keytype
27 #include <g4main/PHG4VtxPoint.h>
28 
29 #include "AssocInfoContainer.h"
30 
31 using namespace std;
32 
33 //____________________________________________________________________________..
35  SubsysReco(name)
36 {
37  //cout << "PHTruthSiliconAssociation::PHTruthSiliconAssociation(const std::string &name) Calling ctor" << endl;
38 }
39 
40 //____________________________________________________________________________..
42 {
43 
44 }
45 
46 //____________________________________________________________________________..
48 {
49 
51 }
52 
53 //____________________________________________________________________________..
55 {
56  int ret = GetNodes(topNode);
57 
58  return ret;
59 }
60 
61 //____________________________________________________________________________..
63 {
64  if (Verbosity() >= 1)
65  cout << "PHTruthSiliconAssociation::process_event(PHCompositeNode *topNode) Processing Event" << endl;
66 
67  // Loop over all SvtxTracks from the CA seeder
68  // These should contain all TPC clusters already
69 
70  const unsigned int original_track_map_lastkey = _track_map->empty() ? 0:std::prev(_track_map->end())->first;
71  if(Verbosity() > 0) std::cout << " track map last key = " << original_track_map_lastkey << std::endl;
72 
73  for (auto phtrk_iter = _track_map->begin();
74  phtrk_iter != _track_map->end();
75  ++phtrk_iter)
76  {
77  // we may add tracks to the map, so we stop at the last original track
78  if(phtrk_iter->first > original_track_map_lastkey) break;
79 
80  _tracklet = phtrk_iter->second;
81  if (Verbosity() >= 1)
82  {
83  std::cout
84  << __LINE__
85  << ": Processing seed itrack: " << phtrk_iter->first
86  << ": nhits: " << _tracklet-> size_cluster_keys()
87  << ": Total tracks: " << _track_map->size()
88  << ": phi: " << _tracklet->get_phi()
89  << endl;
90  }
91 
92  // this one is always there
94 
95  // identify the best truth track match(es) for this seed track
96  std::vector<PHG4Particle*> g4particle_vec = getG4PrimaryParticle(_tracklet);
97  if(Verbosity() > 0) std::cout << " g4particle_vec.size() " << g4particle_vec.size() << std::endl;
98 
99  if(g4particle_vec.size() < 1) continue;
100 
101  bool test_phi_matching = false; // normally false
102  if(test_phi_matching)
103  {
104  // for getting the pT dependence of dphi to eliminate the bias in phi from PHTpcTracker
105  if(g4particle_vec.size() == 1)
106  {
107  // print out the eta and phi values for analysis in case where match is unique
108 
109  double si_phi = atan2(g4particle_vec[0]->get_py(), g4particle_vec[0]->get_px());
110  double si_eta = asinh(g4particle_vec[0]->get_pz() / sqrt( pow(g4particle_vec[0]->get_px(),2)+pow( g4particle_vec[0]->get_py(),2)));
111  double si_pt = sqrt(pow(g4particle_vec[0]->get_px(),2)+pow( g4particle_vec[0]->get_py(),2));
112 
113  double tpc_phi = atan2(_tracklet->get_py(), _tracklet->get_px());
114  double tpc_eta = _tracklet->get_eta();
115 
116  cout << " Try: " << " pt " << si_pt << " tpc_phi " << tpc_phi << " si_phi " << si_phi
117  << " tpc_eta " << tpc_eta << " si_eta " << si_eta << endl;
118 
119  }
120  }
121 
122  // make copies of the original track for later use
123  std::vector<SvtxTrack*> extraTrack;
124  for(unsigned int ig4=0;ig4 < g4particle_vec.size()-1; ++ig4)
125  {
126  SvtxTrack *newTrack = new SvtxTrack_v2();
127  // Not the first g4particle in the list, we need to add a new copy of the track to the track map and add the silicon clusters to that
128  const unsigned int lastTrackKey = ( _track_map->empty() ? 0:std::prev(_track_map->end())->first ) + ig4;
129  //std::cout << " extra track key " << lastTrackKey + 1 << std::endl;
130 
131  newTrack->set_id(lastTrackKey + 1);
132  _seed_track_map->addAssoc(_tracklet->get_id(), newTrack->get_id()) ;
133 
134  newTrack->set_charge(_tracklet->get_charge());
135  newTrack->set_px(_tracklet->get_px());
136  newTrack->set_py(_tracklet->get_py());
137  newTrack->set_pz(_tracklet->get_pz());
138  newTrack->set_x(_tracklet->get_x());
139  newTrack->set_y(_tracklet->get_y());
140  newTrack->set_z(_tracklet->get_z());
141  for(int i = 0; i < 6; ++i)
142  {
143  for(int j = 0; j < 6; ++j)
144  {
145  newTrack->set_error(i,j, _tracklet->get_error(i,j));
146  }
147  }
148 
149  // loop over associated clusters to get hits for original TPC track, copy to new track
151  iter != _tracklet->end_cluster_keys();
152  ++iter)
153  {
154  TrkrDefs::cluskey cluster_key = *iter;
155  newTrack->insert_cluster_key(cluster_key);
156  }
157 
158  extraTrack.push_back(newTrack);
159  }
160 
161  // we just add the silicon clusters to the original track for the first returned g4particle
162 
163  if (Verbosity() >= 1)
164  {
165  std::cout << " original track:" << endl;
166  _tracklet->identify();
167  }
168 
169  // identify the clusters that are associated with this g4particle
170  PHG4Particle* g4particle = g4particle_vec[0];
171  std::set<TrkrDefs::cluskey> clusters = getSiliconClustersFromParticle(g4particle);
172 
173  for (std::set<TrkrDefs::cluskey>::iterator jter = clusters.begin();
174  jter != clusters.end();
175  ++jter)
176  {
177  TrkrDefs::cluskey cluster_key = *jter;
178  unsigned int layer = TrkrDefs::getLayer(cluster_key);
179  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
180  if (Verbosity() >= 1) std::cout << " found silicon cluster with key: " << cluster_key << " in layer " << layer << std::endl;
181 
182  // Identify the MVTX and INTT clusters and add them to the SvtxTrack cluster key list
183  if(trkrid == TrkrDefs::mvtxId || trkrid == TrkrDefs::inttId)
184  {
185  if (Verbosity() >= 1) std::cout << " cluster belongs to MVTX or INTT, add to track " << std::endl;
186  _tracklet->insert_cluster_key(cluster_key);
188  }
189  }
190 
191  if (Verbosity() >= 1)
192  {
193  std::cout << " updated original track:" << std::endl;
194  _tracklet->identify();
195  std::cout << " new cluster keys size " << _tracklet->size_cluster_keys() << endl;
196  std::cout << "Done with original track " << phtrk_iter->first << std::endl;
197  }
198 
199  // now add any extra copies of the track
200  if(g4particle_vec.size() > 1)
201  {
202  for(unsigned int ig4=0;ig4 < g4particle_vec.size()-1; ++ ig4)
203  {
204  PHG4Particle* g4particle = g4particle_vec[ig4];
205 
206  const int vertexID = g4particle->get_vtx_id();
207 
208  if (Verbosity() >= 1)
209  std::cout << " ig4 " << ig4 << " g4particleID " << g4particle->get_track_id()
210  << " px " << g4particle->get_px()
211  << " py " << g4particle->get_py()
212  << " phi " << atan2(g4particle->get_py(), g4particle->get_px() )
213  << std::endl;
214 
215  // identify the clusters that are associated with this g4particle
216  std::set<TrkrDefs::cluskey> clusters = getSiliconClustersFromParticle(g4particle);
217 
218  // Not the first g4particle in the list, we need to add a new copy of the track to the track map and add the silicon clusters to that
219  const unsigned int lastTrackKey = _track_map->empty() ? 0:std::prev(_track_map->end())->first;
220  //std::cout << " lastTrackKey " << lastTrackKey + 1 << std::endl;
221 
222  extraTrack[ig4]->set_id(lastTrackKey + 1);
223 
224  if (Verbosity() >= 1)
225  {
226  std::cout << " original (copy) track:" << endl;
227  extraTrack[ig4]->identify();
228  }
229 
230  _track_map->insert(extraTrack[ig4]);
231 
232  for (std::set<TrkrDefs::cluskey>::iterator jter = clusters.begin();
233  jter != clusters.end();
234  ++jter)
235  {
236  TrkrDefs::cluskey cluster_key = *jter;
237  unsigned int layer = TrkrDefs::getLayer(cluster_key);
238  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
239  if (Verbosity() >= 1) std::cout << " found cluster with key: " << cluster_key << " in layer " << layer << std::endl;
240 
241  // Identify the MVTX and INTT clusters and add them to the SvtxTrack cluster key list
242  if(trkrid == TrkrDefs::mvtxId || trkrid == TrkrDefs::inttId)
243  {
244  if (Verbosity() >= 1) std::cout << " cluster belongs to MVTX or INTT, add to track " << std::endl;
245  extraTrack[ig4]->insert_cluster_key(cluster_key);
246  _assoc_container->SetClusterTrackAssoc(cluster_key, extraTrack[ig4]->get_id());
247  }
248  }
249 
250  // set the combined track vertex to the truth vertex for the silicon truth track
251  double random = 0;//((double) rand() / (RAND_MAX)) * 0.05;
252  // make it negative sometimes
253  if(rand() % 2)
254  random *= -1;
255  // assign the track position using the truth vertex for this track
256  auto g4vertex = _g4truth_container->GetVtx(vertexID);
257  extraTrack[ig4]->set_x(g4vertex->get_x() * (1 + random));
258  extraTrack[ig4]->set_y(g4vertex->get_y() * (1 + random));
259  extraTrack[ig4]->set_z(g4vertex->get_z() * (1 + random));
260 
261  if (Verbosity() >= 1)
262  {
263  std::cout << " updated additional track:" << std::endl;
264  extraTrack[ig4]->identify();
265  std::cout << " new cluster keys size " << extraTrack[ig4]->size_cluster_keys() << endl;
266  std::cout << "Done with extra track " << extraTrack[ig4]->get_id() << std::endl;
267  }
268  }
269  }
270  }
271 
272  // loop over all tracks and copy the silicon clusters to the corrected cluster map
275 
276  if (Verbosity() >= 1)
277  cout << "PHTruthSiliconAssociation::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
278 
280 }
281 
282 //____________________________________________________________________________..
284 {
285  //cout << "PHTruthSiliconAssociation::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << endl;
287 }
288 
289 //____________________________________________________________________________..
290 int PHTruthSiliconAssociation::EndRun(const int /*runnumber*/)
291 {
293 }
294 
295 //____________________________________________________________________________..
297 {
299 }
300 
301 //____________________________________________________________________________..
303 {
305 }
306 
307 //____________________________________________________________________________..
308 void PHTruthSiliconAssociation::Print(const std::string &/*what*/) const
309 {
310  //cout << "PHTruthSiliconAssociation::Print(const std::string &what) const Printing info for " << what << endl;
311 }
312 
314 {
315  //---------------------------------
316  // Get Objects off of the Node Tree
317  //---------------------------------
318 
319  _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
320  if (!_g4truth_container)
321  {
322  cerr << PHWHERE << " ERROR: Can't find node G4TruthInfo" << endl;
324  }
325 
326  _corrected_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"CORRECTED_TRKR_CLUSTER");
328  {
329  std::cout << " Found CORRECTED_TRKR_CLUSTER node " << std::endl;
330  }
331 
332  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
333  if (!_cluster_map)
334  {
335  cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << endl;
337  }
338  _hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
339  if(!_hitsets)
340  {
341  std::cout << PHWHERE << "No hitset container on node tree. Bailing."
342  << std::endl;
344  }
345 
346  _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
347  if (!_track_map)
348  {
349  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << endl;
351  }
352 
353  _assoc_container = findNode::getClass<AssocInfoContainer>(topNode, "AssocInfoContainer");
354  if (!_assoc_container)
355  {
356  cerr << PHWHERE << " ERROR: Can't find AssocInfoContainer." << endl;
358  }
359 
360  _hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
361  if (!_hit_truth_map)
362  {
363  cerr << PHWHERE << " ERROR: Can't find TrkrHitTruthAssoc." << endl;
365  }
366 
367  _cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
368  if (!_cluster_hit_map)
369  {
370  cerr << PHWHERE << " ERROR: Can't find TrkrClusterHitAssoc." << endl;
372  }
373 
374 
376  PHNodeIterator iter(topNode);
377  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
378 
380  if (!dstNode)
381  {
382  std::cerr << "DST Node missing, quitting" << std::endl;
383  throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes");
384  }
385 
387  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
388 
390  if (!svtxNode)
391  {
392  svtxNode = new PHCompositeNode("SVTX");
393  dstNode->addNode(svtxNode);
394  }
395 
399  svtxNode->addNode(node);
400 
401  _g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
402  _g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
403  _g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
404 
406 }
407 
409 {
410  // Find the best g4particle match to the clusters associated with this reco track
411 
412  std::vector<PHG4Particle*> g4part_vec;
413  std::set<int> pid;
414  std::multimap<int, int> pid_count;
415  int minfound = 200; // require at least this many hits to be put on the list of contributing particles
416 
417  // loop over associated clusters to get hits
419  iter != track->end_cluster_keys();
420  ++iter)
421  {
422  TrkrDefs::cluskey cluster_key = *iter;
423 
424  // get all reco hits for this cluster
425  //TrkrClusterHitAssoc::ConstRange
426  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
427  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
428  //for(TrkrClusterHitAssoc::ConstIterator clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
429  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
430  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
431  {
432  TrkrDefs::hitkey hitkey = clushititer->second;
433  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
435 
436  // get all of the g4hits for this hitkey
437  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
438  _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
439  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
440  {
441  // extract the g4 hit key
442  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
443  PHG4Hit * g4hit = nullptr;
444  unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
445  switch( trkrid )
446  {
447  case TrkrDefs::tpcId: g4hit = _g4hits_tpc->findHit(g4hitkey); break;
448  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(g4hitkey); break;
449  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(g4hitkey); break;
450  default: break;
451  }
452  if( g4hit )
453  {
454  // get the g4particle ID for this g4hit
455  pid.insert(g4hit->get_trkid());
456  // this is for counting the number of times this g4particle was associated
457  int cnt = 1;
458  pid_count.insert(std::make_pair(g4hit->get_trkid(), cnt));
459  }
460  } // end loop over g4hits associated with hitsetkey and hitkey
461  } // end loop over hits associated with cluskey
462  } // end loop over clusters associated with this track
463 
464 
465  // loop over the particle id's, and count the number for each one
466  for( auto it = pid.begin(); it != pid.end(); ++it)
467  {
468  if(*it < 0) continue; // ignore secondary particles
469 
470  int nfound = 0;
471  std::pair<std::multimap<int, int>::iterator, std::multimap<int,int>::iterator> this_pid = pid_count.equal_range(*it);
472  for(auto cnt_it = this_pid.first; cnt_it != this_pid.second; ++cnt_it)
473  {
474  nfound++;
475  }
476 
477  if(Verbosity() >= 1) std::cout << " pid: " << *it << " nfound " << nfound << std::endl;
478  if(nfound > minfound)
479  {
480  g4part_vec.push_back(_g4truth_container->GetParticle(*it));
481  }
482  }
483 
484  return g4part_vec;
485 
486 }
487 
489 {
490  // Find the reco clusters in the silicon layers from this g4particle
491 
492  std::set<TrkrDefs::cluskey> clusters;
493 
494  // loop over all the clusters
495  auto hitsetrange = _hitsets->getHitSets();
496  for (auto hitsetitr = hitsetrange.first;
497  hitsetitr != hitsetrange.second;
498  ++hitsetitr){
499  auto range = _cluster_map->getClusters(hitsetitr->first);
500  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
501  TrkrDefs::cluskey cluster_key = clusIter->first;
502  unsigned int layer = TrkrDefs::getLayer(cluster_key);
503 
504  if(layer > 6) continue; // we need the silicon layers only
505 
506  // get all truth hits for this cluster
507  //TrkrClusterHitAssoc::ConstRange
508  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
509  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
510  //for(TrkrClusterHitAssoc::ConstIterator
511  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
512  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
513  {
514  TrkrDefs::hitkey hitkey = clushititer->second;
515  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
517 
518  // get all of the g4hits for this hitkey
519  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
520  _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
521  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
522  {
523  // extract the g4 hit key
524  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
525  PHG4Hit * g4hit = nullptr;
526  unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
527  switch( trkrid )
528  {
529  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(g4hitkey); break;
530  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(g4hitkey); break;
531  default: break;
532  }
533  if( g4hit )
534  {
535  // get the g4particle for this g4hit
536  int trkid = g4hit->get_trkid();
537  if(trkid == g4particle->get_track_id())
538  {
539  clusters.insert(cluster_key);
540  }
541  }
542  } // end loop over g4hits associated with hitsetkey and hitkey
543  } // end loop over hits associated with cluskey
544  } // end loop over all cluster keys in silicon layers
545  }//end loop over hitsets
546  return clusters;
547 
548 }
549 
551 {
552  // loop over final track map, copy silicon clusters to corrected cluster map
553  for (auto phtrk_iter = _track_map->begin();
554  phtrk_iter != _track_map->end();
555  ++phtrk_iter)
556  {
557  SvtxTrack *track = phtrk_iter->second;
558 
559  // loop over associated clusters to get keys for silicon cluster
561  iter != track->end_cluster_keys();
562  ++iter)
563  {
564  TrkrDefs::cluskey cluster_key = *iter;
565  unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
566  if(trkrid == TrkrDefs::mvtxId || trkrid == TrkrDefs::inttId)
567  {
568  TrkrCluster *cluster = _cluster_map->findCluster(cluster_key);
569  TrkrCluster *newclus = _corrected_cluster_map->findOrAddCluster(cluster_key)->second;
570 
571  newclus->setSubSurfKey(cluster->getSubSurfKey());
572  newclus->setAdc(cluster->getAdc());
573 
574  newclus->setActsLocalError(0,0,cluster->getActsLocalError(0,0));
575  newclus->setActsLocalError(1,0,cluster->getActsLocalError(1,0));
576  newclus->setActsLocalError(0,1,cluster->getActsLocalError(0,1));
577  newclus->setActsLocalError(1,1,cluster->getActsLocalError(1,1));
578 
579  newclus->setLocalX(cluster->getLocalX());
580  newclus->setLocalY(cluster->getLocalY());
581  }
582  }
583  }
584 }