ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHSiliconTruthTrackSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHSiliconTruthTrackSeeding.cc
2 
3 #include "AssocInfoContainer.h"
4 
5 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack, SvtxTra...
6 #include <trackbase_historic/SvtxTrackMap.h> // for SvtxTrackMap, Svtx...
9 
13 #include <trackbase/TrkrDefs.h>
15 #include <trackbase/TrkrHitSet.h>
17 
18 #include <g4eval/SvtxClusterEval.h>
19 #include <g4eval/SvtxEvalStack.h>
20 #include <g4eval/SvtxEvaluator.h>
21 #include <g4eval/SvtxTrackEval.h>
22 
23 #include <g4main/PHG4Hit.h> // for PHG4Hit
24 #include <g4main/PHG4Particle.h> // for PHG4Particle
26 #include <g4main/PHG4HitDefs.h> // for keytype
28 
30 
31 #include <phool/getClass.h>
32 #include <phool/phool.h>
33 
34 #include <cstdlib> // for exit
35 #include <iostream> // for operator<<, endl
36 #include <map> // for multimap, map<>::c...
37 #include <memory>
38 #include <utility> // for pair
39 #include <cassert>
40 #include <set>
41 
42 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
43 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
44 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << std::endl
45 
46 class PHCompositeNode;
47 
48 using namespace std;
49 
51  : PHTrackSeeding(name)
52 {
53  // we use a separate track map node for the silicon track stubs
54  // std::string track_map_node_name = {"SvtxSiliconTrackMap"};
55  // set_track_map_name(track_map_node_name);
56 }
57 
59 {
60  if(Verbosity() > 0)
61  std::cout << "Enter PHSiliconTruthTrackSeeding:: Setup" << std::endl;
62 
63  // we use a separate track map node for the silicon track stubs
64  std::string track_map_node_name = {"SvtxSiliconTrackMap"};
65  set_track_map_name(track_map_node_name);
66 
67  int ret = PHTrackSeeding::Setup(topNode);
68  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
69 
70  // ret = GetNodes(topNode);
71  // if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
72 
74 }
75 
77 {
78  // We start with all reco clusters in the silicon layers
79  // get the associated truth g4hits and from that the g4particle
80  // make a map of <g4particle, clusters>
81  // Make an SvtxTrack in silicon from g4particle
82  // get momentum from truth
83  // get position from vertex
84  // add clusters to it
85  // Add the track to SvtxSiliconTrackMap
86 
87  typedef std::map<int, std::set<TrkrCluster*> > TrkClustersMap;
88  TrkClustersMap m_trackID_clusters;
89 
90  // loop over all clusters
91  auto hitsetrange = _hitsets->getHitSets();
92  for (auto hitsetitr = hitsetrange.first;
93  hitsetitr != hitsetrange.second;
94  ++hitsetitr){
95  auto range = _cluster_map->getClusters(hitsetitr->first);
96  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
97  TrkrCluster* cluster = clusIter->second;
98  TrkrDefs::cluskey cluskey = clusIter->first;
99  unsigned int trkrid = TrkrDefs::getTrkrId(cluskey);
100 
101  if(trkrid != TrkrDefs::mvtxId && trkrid != TrkrDefs::inttId) continue;
102 
103  unsigned int layer = TrkrDefs::getLayer(cluskey);
104  if(layer > 6)
105  std::cout << PHWHERE << " Warning: layer is screwed up - layer = " << layer << std::endl;
106 
107  if (Verbosity() >= 3)
108  {
109  cout << PHWHERE <<" process cluster ";
110  cluster->identify();
111  }
112 
113  // get the hits for this cluster
114  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
115  hitrange = clusterhitassoc->getHits(cluskey); // returns range of pairs {cluster key, hit key} for this cluskey
116  //for (TrkrClusterHitAssoc::ConstIterator clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
117  for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
118  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
119  {
120  TrkrDefs::hitkey hitkey = clushititer->second;
121  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
123 
124  if (Verbosity() >= 3)
125  {
126  cout << PHWHERE <<" --- process hit with hitkey " << hitkey << " hitsetkey " << hitsetkey << std::endl;
127  }
128 
129  // get all of the g4hits for this hitkey
130  std::multimap<TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
131  hittruthassoc->getG4Hits(hitsetkey, hitkey, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
132  for (std::multimap<TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
133  {
134  // extract the g4 hit key here and add the hits to the set
135  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
136 
137  if (Verbosity() >= 3)
138  {
139  cout << PHWHERE <<" --- process g4hit with key " << g4hitkey << std::endl;
140  }
141 
142  PHG4Hit* phg4hit = nullptr;
143  switch( trkrid )
144  {
145  case TrkrDefs::mvtxId:
146  if (phg4hits_mvtx) phg4hit = phg4hits_mvtx->findHit( g4hitkey );
147  break;
148 
149  case TrkrDefs::inttId:
150  if (phg4hits_intt) phg4hit = phg4hits_intt->findHit( g4hitkey );
151  break;
152  }
153 
154  if( !phg4hit )
155  {
156  std::cout<<PHWHERE<<" unable to find g4hit from key " << g4hitkey << std::endl;
157  continue;
158  }
159 
160  int particle_id = phg4hit->get_trkid();
161 
162  // momentum cut-off
163  if (_min_momentum>0)
164  {
166  if (!particle)
167  {
168  cout << PHWHERE <<" - validity check failed: missing truth particle with ID of "<<particle_id<<". Exiting..."<<endl;
169  exit(1);
170  }
171  const double monentum2 =
172  particle->get_px() * particle->get_px()
173  +
174  particle->get_py() * particle->get_py()
175  +
176  particle->get_pz() * particle->get_pz()
177  ;
178 
179  if (Verbosity() >= 10)
180  {
181  cout << PHWHERE <<" --- check momentum for g4particle "<<particle_id<<" -> cluster "<<cluskey
182  <<" = "<<sqrt(monentum2)<<endl;;
183  particle->identify();
184  }
185 
186  if (monentum2 < _min_momentum * _min_momentum)
187  {
188  if (Verbosity() >= 3)
189  {
190  cout << PHWHERE <<" ignore low momentum g4particle "<<particle_id<<" -> cluster "<<cluskey<<endl;;
191  particle->identify();
192  }
193  continue;
194  }
195  }
196 
197  // add or append this particle/cluster combination to the map
198  TrkClustersMap::iterator it = m_trackID_clusters.find(particle_id);
199 
200  if (it != m_trackID_clusters.end())
201  {
202  // the clusters are stored in a set, no need to check if it is in there already
203  it->second.insert(cluster);
204  if (Verbosity() >= 3)
205  {
206  cout << PHWHERE <<" --- appended to g4particle"<<particle_id<<" new cluster "<<cluskey<<endl;;
207  //cluster->identify();
208  }
209  }
210  else
211  {
212  std::set<TrkrCluster*> clusters;
213  clusters.insert(cluster);
214  m_trackID_clusters.insert(std::pair<int, std::set<TrkrCluster*> >(particle_id, clusters));
215 
216 
217  if (Verbosity() >= 3)
218  {
219  cout << PHWHERE <<" --- added new g4particle "<<particle_id<<" and inserted cluster "<<cluskey<<endl;;
220  //cluster->identify();
221  }
222 
223  }
224  } // loop over g4hits associated with hit
225  } // loop over hits associated with cluster
226  } // loop over clusters
227  } // loop over hitsets
228  //==================================
229  // make the tracks
230 
231  if (Verbosity() >= 2)
232  {
233  cout <<__PRETTY_FUNCTION__
234  <<" Beginning _track_map->size = "<<_track_map->size()<<endl;
235 
237  }
238 
239  // Build tracks, looping over assembled list of truth track ID's and associated reco clusters
240  for (TrkClustersMap::const_iterator trk_clusters_itr = m_trackID_clusters.begin();
241  trk_clusters_itr != m_trackID_clusters.end(); ++trk_clusters_itr)
242  {
243  if (trk_clusters_itr->second.size() < _min_clusters_per_track)
244  continue;
245 
246  // check number of layers also pass the _min_clusters_per_track cut to avoid tight loopers
247  set<uint8_t> layers;
248  for (TrkrCluster* cluster : trk_clusters_itr->second)
249  {
250  assert(cluster);
251  const uint8_t layer = TrkrDefs::getLayer(cluster->getClusKey());
252  layers.insert(layer);
253  }
254  if (Verbosity()>2)
255  {
256  cout <<__PRETTY_FUNCTION__<<" particle "<<trk_clusters_itr->first<<" -> "
257  <<trk_clusters_itr->second.size()<<" clusters covering "<<layers.size()<<" layers."
258  <<" Layer/clusters cuts are > "<<_min_clusters_per_track
259  <<endl;
260  }
261 
262  if (layers.size() >= _min_clusters_per_track)
263  {
264 
265  auto svtx_track = std::make_unique<SvtxTrack_FastSim_v2>();
266  svtx_track->set_id(_track_map->size());
267  svtx_track->set_truth_track_id(trk_clusters_itr->first);
268 
269  // assign truth particle vertex ID to this silicon track stub
270  PHG4Particle* particle = _g4truth_container->GetParticle(trk_clusters_itr->first);
271  int vertexId = particle->get_vtx_id() - 1; // Geant likes to count from 1 for both vertex and g4particle ID, _vertex_map counts from 0
272  if(vertexId < 0)
273  {
274  // Secondary particle, will have the same gembed value as the corresponding primary vertex
275  int track_embed = _g4truth_container->isEmbeded(trk_clusters_itr->first);
276  auto vrange = _g4truth_container->GetPrimaryVtxRange();
277  for (auto iter = vrange.first; iter != vrange.second; ++iter) // all primary vertexes
278  {
279  const int point_id = iter->first;
280  int vert_embed = _g4truth_container->isEmbededVtx(point_id);
281  if(vert_embed == track_embed)
282  vertexId = point_id - 1; // Geant starts counting vertices at 1
283  if(Verbosity() > 3)
284  std::cout << " track_embed " << track_embed << " vert_embed " << vert_embed << " G4 point id " << point_id << " SvtxMap vertexId " << vertexId << std::endl;
285  }
286  }
287  if(vertexId < 0) // should not be possible
288  vertexId = 0;
289 
290  svtx_track->set_vertex_id(vertexId);
291 
292  if(Verbosity() > 0)
293  std::cout << " truth track G4 point id " << particle->get_vtx_id() << " becomes SvtxMap id " << vertexId
294  << " gembed is " << _g4truth_container->isEmbeded(trk_clusters_itr->first) << " for truth particle " << trk_clusters_itr->first << std::endl;
295 
296  // set the track position to the vertex position
297  const SvtxVertex *svtxVertex = _vertex_map->get(vertexId);
298  if(!svtxVertex)
299  {
300  std::cout << PHWHERE << "Failed to get vertex with ID " << vertexId << " from _vertex_map, cannot proceed - skipping this silicon track" << std::endl;
301  continue;
302  }
303  svtx_track->set_x(svtxVertex->get_x());
304  svtx_track->set_y(svtxVertex->get_y());
305  svtx_track->set_z(svtxVertex->get_z());
306 
307  if(Verbosity() > 0)
308  {
309  std::cout << " truth track SvtxMap vertexid is " << vertexId << " for truth particle " << trk_clusters_itr->first << std::endl;
310  std::cout << " track position x,y,z = " << svtxVertex->get_x() << ", " << svtxVertex->get_y() << ", " << svtxVertex->get_z() << std::endl;
311  }
312 
313  // set momentum to the truth value
314  svtx_track->set_px(particle->get_px());
315  svtx_track->set_py(particle->get_py());
316  svtx_track->set_pz(particle->get_pz());
317 
318  if(Verbosity() > 0)
319  std::cout << PHWHERE << "Truth track ID is " << svtx_track->get_truth_track_id() << " particle ID is " << particle->get_pid() << std::endl;
320 
321  // add the silicon clusters
322  for (TrkrCluster* cluster : trk_clusters_itr->second)
323  {
324  svtx_track->insert_cluster_key(cluster->getClusKey());
325  _assoc_container->SetClusterTrackAssoc(cluster->getClusKey(), svtx_track->get_id());
326  }
327 
328  _track_map->insert(svtx_track.get());
329 
330  if (Verbosity() >= 2)
331  {
332  cout <<__PRETTY_FUNCTION__<<" particle "<<trk_clusters_itr->first<<" -> "
333  <<trk_clusters_itr->second.size()<<" clusters"
334  <<" _track_map->size = "<< (_track_map->size()) <<": ";
335  _track_map->identify();
336  }
337  }
338  }
339 
340  if (Verbosity() >= 2)
341  {
342  cout << "Loop over SvtxSiliconTrackMap entries " << endl;
343  for (SvtxTrackMap::Iter iter = _track_map->begin();
344  iter != _track_map->end(); ++iter)
345  {
346  SvtxTrack* svtx_track = iter->second;
347 
348  //svtx_track->identify();
349  //continue;
350 
351  cout << "Track ID: " << svtx_track->get_id() << ", Track pT: "
352  << svtx_track->get_pt() << ", Truth Track/Particle ID: "
353  << svtx_track->get_truth_track_id()
354  << " vertex ID: " << svtx_track->get_vertex_id() << endl;
355  cout << " (px, py, pz) = " << "("<< svtx_track->get_px() << ", " << svtx_track->get_py() << ", " << svtx_track->get_pz() << ")" << endl;
356  cout << " (x, y, z) = " << "("<< svtx_track->get_x() << ", " << svtx_track->get_y() << ", " << svtx_track->get_z() << ")" << endl;
357 
358  //Print associated clusters;
360  svtx_track->begin_cluster_keys();
361  iter != svtx_track->end_cluster_keys(); ++iter)
362  {
363  TrkrDefs::cluskey cluster_key = *iter;
364  TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
365  float radius = sqrt(
366  cluster->getX() * cluster->getX() + cluster->getY() * cluster->getY());
367  cout << " cluster ID: "
368  << cluster->getClusKey() << ", cluster radius: " << radius
369  << endl;
370  }
371  }
372  }
373  //==================================
374 
376 }
377 
379 {
380  /*
381  // If the _use_truth_clusters flag is set, the cluster map now points to the truth clusters
382  // in this module we want to use the reco silicon clusters instead, so we move the pointer
383  if(_use_truth_clusters)
384  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
385 
386  if (!_cluster_map)
387  {
388  cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << endl;
389  return Fun4AllReturnCodes::ABORTEVENT;
390  }
391  */
392 
393  _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
394  if (!_g4truth_container)
395  {
396  cerr << PHWHERE << " ERROR: Can't find node G4TruthInfo" << endl;
398  }
399 
400  clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
401  if (!clusterhitassoc)
402  {
403  cout << PHWHERE << "Failed to find TRKR_CLUSTERHITASSOC node, quit!" << endl;
404  exit(1);
405  }
406 
407  hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
408  if (!hittruthassoc)
409  {
410  cout << PHWHERE << "Failed to find TRKR_HITTRUTHASSOC node, quit!" << endl;
411  exit(1);
412  }
413 
414  using nodePair = std::pair<std::string, PHG4HitContainer*&>;
415  std::initializer_list<nodePair> nodes =
416  {
417  { "G4HIT_TPC", phg4hits_tpc },
418  { "G4HIT_INTT", phg4hits_intt },
419  { "G4HIT_MVTX", phg4hits_mvtx },
420  { "G4HIT_MICROMEGAS", phg4hits_micromegas }
421  };
422 
423  for( auto&& node: nodes )
424  {
425  if( !( node.second = findNode::getClass<PHG4HitContainer>( topNode, node.first ) ) )
426  { std::cerr << PHWHERE << " PHG4HitContainer " << node.first << " not found" << std::endl; }
427  }
428 
430 }
431 
433 {
434  return 0;
435 }