ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackEvaluation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackEvaluation.cc
1 
6 #include "TrackEvaluation.h"
7 
12 #include <g4main/PHG4Hit.h>
13 #include <g4main/PHG4Hitv1.h>
15 #include <g4main/PHG4Particle.h>
17 #include <intt/InttDefs.h>
20 #include <mvtx/MvtxDefs.h>
21 #include <phool/getClass.h>
22 #include <phool/PHCompositeNode.h>
23 #include <phool/PHNodeIterator.h>
24 #include <tpc/TpcDefs.h>
27 #include <trackbase/TrkrDefs.h>
28 #include <trackbase/TrkrCluster.h>
31 #include <trackbase/TrkrHit.h>
32 #include <trackbase/TrkrHitSet.h>
37 
38 #include <TVector3.h>
39 
40 #include <algorithm>
41 #include <bitset>
42 #include <cassert>
43 #include <iostream>
44 #include <numeric>
45 
46 //_____________________________________________________________________
47 namespace
48 {
49 
51  template<class T> class range_adaptor
52  {
53  public:
54  range_adaptor( const T& range ):m_range(range){}
55  inline const typename T::first_type& begin() {return m_range.first;}
56  inline const typename T::second_type& end() {return m_range.second;}
57  private:
58  T m_range;
59  };
60 
62  template<class T> inline constexpr T square( const T& x ) { return x*x; }
63 
65  template<class T> inline constexpr T get_r( T x, T y ) { return std::sqrt( square(x) + square(y) ); }
66 
68  template<class T> T get_pt( T px, T py ) { return std::sqrt( square(px) + square(py) ); }
69 
71  template<class T> T get_p( T px, T py, T pz ) { return std::sqrt( square(px) + square(py) + square(pz) ); }
72 
74  template<class T> T get_eta( T p, T pz ) { return std::log( (p+pz)/(p-pz) )/2; }
75 
77  struct interpolation_data_t
78  {
79  using list = std::vector<interpolation_data_t>;
80  double x() const { return position.x(); }
81  double y() const { return position.y(); }
82  double z() const { return position.z(); }
83 
84  double px() const { return momentum.x(); }
85  double py() const { return momentum.y(); }
86  double pz() const { return momentum.z(); }
87 
88  TVector3 position;
89  TVector3 momentum;
90  double weight = 1;
91  };
92 
94  template<double (interpolation_data_t::*accessor)() const>
95  double average( const interpolation_data_t::list& hits )
96  {
97  // calculate all terms needed for the interpolation
98  // need to use double everywhere here due to numerical divergences
99  double sw = 0;
100  double swx = 0;
101 
102  bool valid( false );
103  for( const auto& hit:hits )
104  {
105 
106  const double x = (hit.*accessor)();
107  if(std::isnan(x)) continue;
108 
109  const double w = hit.weight;
110  if( w <= 0 ) continue;
111 
112  valid = true;
113  sw += w;
114  swx += w*x;
115  }
116 
117  if( !valid ) return NAN;
118  return swx/sw;
119  }
120 
122  inline int is_primary( PHG4Particle* particle )
123  { return particle->get_parent_id() == 0; }
124 
126  int64_t get_mask( SvtxTrack* track )
127  { return std::accumulate( track->begin_cluster_keys(), track->end_cluster_keys(), int64_t(0),
128  []( int64_t value, const TrkrDefs::cluskey& key ) {
129  return TrkrDefs::getLayer(key)<64 ? value|(1LL<<TrkrDefs::getLayer(key)) : value;
130  } );
131  }
132 
134  template<int type>
135  int get_clusters( SvtxTrack* track )
136  {
137  return std::count_if( track->begin_cluster_keys(), track->end_cluster_keys(),
138  []( const TrkrDefs::cluskey& key ) { return TrkrDefs::getTrkrId(key) == type; } );
139  }
140 
143  {
145 
146  trackStruct.charge = track->get_charge();
147  trackStruct.nclusters = track->size_cluster_keys();
148  trackStruct.mask = get_mask( track );
149  trackStruct.nclusters_mvtx = get_clusters<TrkrDefs::mvtxId>( track );
150  trackStruct.nclusters_intt = get_clusters<TrkrDefs::inttId>( track );
151  trackStruct.nclusters_tpc = get_clusters<TrkrDefs::tpcId>( track );
152  trackStruct.nclusters_micromegas = get_clusters<TrkrDefs::micromegasId>( track );
153 
154  trackStruct.chisquare = track->get_chisq();
155  trackStruct.ndf = track->get_ndf();
156 
157  trackStruct.x = track->get_x();
158  trackStruct.y = track->get_y();
159  trackStruct.z = track->get_z();
160  trackStruct.r = get_r( trackStruct.x, trackStruct.y );
161  trackStruct.phi = std::atan2( trackStruct.y, trackStruct.x );
162 
163  trackStruct.px = track->get_px();
164  trackStruct.py = track->get_py();
165  trackStruct.pz = track->get_pz();
166  trackStruct.pt = get_pt( trackStruct.px, trackStruct.py );
167  trackStruct.p = get_p( trackStruct.px, trackStruct.py, trackStruct.pz );
168  trackStruct.eta = get_eta( trackStruct.p, trackStruct.pz );
169 
170  return trackStruct;
171  }
172 
174  void add_cluster_size( TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key, TrkrClusterHitAssoc* cluster_hit_map )
175  {
176  if( !cluster_hit_map ) return;
177  const auto range = cluster_hit_map->getHits(clus_key);
178 
179  // store full size
180  cluster.size = std::distance( range.first, range.second );
181 
182  const auto detId = TrkrDefs::getTrkrId(clus_key);
183  if(detId == TrkrDefs::micromegasId)
184  {
185 
186  // for micromegas the directional cluster size depends on segmentation type
187  auto segmentation_type = MicromegasDefs::getSegmentationType(clus_key);
188  if( segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_Z ) cluster.z_size = cluster.size;
189  else cluster.phi_size = cluster.size;
190 
191  } else {
192 
193  // for other detectors, one must loop over the constituting hits
194  std::set<int> phibins;
195  std::set<int> zbins;
196  for(const auto& [first, hit_key]:range_adaptor(range))
197  {
198  switch( detId )
199  {
200  default: break;
201  case TrkrDefs::mvtxId:
202  {
203  phibins.insert( MvtxDefs::getRow( hit_key ) );
204  zbins.insert( MvtxDefs::getCol( hit_key ) );
205  break;
206  }
207  case TrkrDefs::inttId:
208  {
209  phibins.insert( InttDefs::getRow( hit_key ) );
210  zbins.insert( InttDefs::getCol( hit_key ) );
211  break;
212  }
213  case TrkrDefs::tpcId:
214  {
215  phibins.insert( TpcDefs::getPad( hit_key ) );
216  zbins.insert( TpcDefs::getTBin( hit_key ) );
217  break;
218  }
219  }
220  }
221  cluster.phi_size = phibins.size();
222  cluster.z_size = zbins.size();
223  }
224  }
225 
226 
228  void add_cluster_energy( TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key,
229  TrkrClusterHitAssoc* cluster_hit_map,
230  TrkrHitSetContainer* hitsetcontainer )
231  {
232 
233  // check container
234  if(!(cluster_hit_map && hitsetcontainer)) return;
235 
236  // for now this is only filled for micromegas
237  const auto detId = TrkrDefs::getTrkrId(clus_key);
238  if(detId != TrkrDefs::micromegasId) return;
239 
240  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(clus_key);
241  const auto hitset = hitsetcontainer->findHitSet( hitset_key );
242  if( !hitset ) return;
243 
244  const auto range = cluster_hit_map->getHits(clus_key);
245  cluster.energy_max = 0;
246  cluster.energy_sum = 0;
247 
248  for( const auto& pair:range_adaptor(range))
249  {
250  const auto hit = hitset->getHit( pair.second );
251  if( hit )
252  {
253  const auto energy = hit->getEnergy();
254  cluster.energy_sum += energy;
255  if( energy > cluster.energy_max ) cluster.energy_max = energy;
256  }
257  }
258 
259  }
260 
261  // ad}d truth information
262  void add_truth_information( TrackEvaluationContainerv1::TrackStruct& track, PHG4Particle* particle )
263  {
264  if( particle )
265  {
266  track.is_primary = is_primary( particle );
267  track.pid = particle->get_pid();
268  track.truth_px = particle->get_px();
269  track.truth_py = particle->get_py();
270  track.truth_pz = particle->get_pz();
271  track.truth_pt = get_pt( track.truth_px, track.truth_py );
272  track.truth_p = get_p( track.truth_px, track.truth_py, track.truth_pz );
273  track.truth_eta = get_eta( track.truth_p, track.truth_pz );
274  }
275  }
276 
277  // calculate intersection between line and circle
278  double line_circle_intersection( const TVector3& p0, const TVector3& p1, double radius )
279  {
280  const double A = square(p1.x() - p0.x()) + square(p1.y() - p0.y());
281  const double B = 2*p0.x()*(p1.x()-p0.x()) + 2*p0.y()*(p1.y()-p0.y());
282  const double C = square(p0.x()) + square(p0.y()) - square(radius);
283  const double delta = square(B)-4*A*C;
284  if( delta < 0 ) return -1;
285 
286  // check first intersection
287  const double tup = (-B + std::sqrt(delta))/(2*A);
288  if( tup >= 0 && tup < 1 ) return tup;
289 
290  // check second intersection
291  const double tdn = (-B-sqrt(delta))/(2*A);
292  if( tdn >= 0 && tdn < 1 ) return tdn;
293 
294  // no valid extrapolation
295  return -1;
296  }
297 
298  // print to stream
299  [[maybe_unused]] std::ostream& operator << (std::ostream& out, const TrackEvaluationContainerv1::ClusterStruct& cluster )
300  {
301  out << "ClusterStruct" << std::endl;
302  out << " cluster: (" << cluster.x << "," << cluster.y << "," << cluster.z << ")" << std::endl;
303  out << " track: (" << cluster.trk_x << "," << cluster.trk_y << "," << cluster.trk_z << ")" << std::endl;
304  out << " truth: (" << cluster.truth_x << "," << cluster.truth_y << "," << cluster.truth_z << ")" << std::endl;
305  return out;
306  }
307 
308  [[maybe_unused]] std::ostream& operator << (std::ostream& out, const TVector3& position)
309  {
310  out << "(" << position.x() << ", " << position.y() << ", " << position.z() << ")";
311  return out;
312  }
313 
314 }
315 
316 //_____________________________________________________________________
318  SubsysReco( name)
319 {}
320 
321 //_____________________________________________________________________
323 {
324  // find DST node
325  PHNodeIterator iter(topNode);
326  auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
327  if (!dstNode)
328  {
329  std::cout << "TrackEvaluation::Init - DST Node missing" << std::endl;
331  }
332 
333  // get EVAL node
334  iter = PHNodeIterator(dstNode);
335  auto evalNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "EVAL"));
336  if( !evalNode )
337  {
338  // create
339  std::cout << "TrackEvaluation::Init - EVAL node missing - creating" << std::endl;
340  evalNode = new PHCompositeNode( "EVAL" );
341  dstNode->addNode(evalNode);
342  }
343 
344  auto newNode = new PHIODataNode<PHObject>( new TrackEvaluationContainerv1, "TrackEvaluationContainer","PHObject");
345  evalNode->addNode(newNode);
346 
348 }
349 
350 //_____________________________________________________________________
353 
354 //_____________________________________________________________________
356 {
357  // load nodes
358  auto res = load_nodes(topNode);
359  if( res != Fun4AllReturnCodes::EVENT_OK ) return res;
360 
361  // cleanup output container
362  if( m_container ) m_container->Reset();
363 
367 
368  // clear maps
369  m_g4hit_map.clear();
371 }
372 
373 //_____________________________________________________________________
376 
377 //_____________________________________________________________________
379 {
380 
381  // acts surface map
382  m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
383  assert( m_surfmaps );
384 
385  // acts geometry
386  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
387  assert( m_tGeometry );
388 
389  // get necessary nodes
390  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
391 
392  // cluster map
393  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
394  if( !m_cluster_map )
395  { m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER"); }
396 
397  // cluster hit association map
398  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
399 
400  // cluster hit association map
401  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,"TRKR_HITTRUTHASSOC");
402 
403  // local container
404  m_container = findNode::getClass<TrackEvaluationContainerv1>(topNode, "TrackEvaluationContainer");
405 
406  // hitset container
407  m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
408 
409  // g4hits
410  m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
411  m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
412  m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
413  m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
414 
415  // g4 truth info
416  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
417 
418  // tpc geometry
419  m_tpc_geom_container = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
420  assert( m_tpc_geom_container );
421 
422  // micromegas geometry
423  m_micromegas_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL" );
424 
426 
427 }
428 
429 //_____________________________________________________________________
431 {
432  if(!m_container) return;
433 
434  // create event struct
436  if( m_hitsetcontainer )
437  {
438  // loop over hitsets
439  for(const auto& [hitsetkey,hitset]:range_adaptor(m_hitsetcontainer->getHitSets()))
440  {
441  const auto trkrId = TrkrDefs::getTrkrId(hitsetkey);
442  const auto layer = TrkrDefs::getLayer(hitsetkey);
444 
445  if(m_cluster_map)
446  {
447 
448  // fill cluster related information
449  const auto clusters = m_cluster_map->getClusters(hitsetkey);
450  const int nclusters = std::distance( clusters.first, clusters.second );
451 
452  switch( trkrId )
453  {
454  case TrkrDefs::mvtxId: event.nclusters_mvtx += nclusters; break;
455  case TrkrDefs::inttId: event.nclusters_intt += nclusters; break;
456  case TrkrDefs::tpcId: event.nclusters_tpc += nclusters; break;
457  case TrkrDefs::micromegasId: event.nclusters_micromegas += nclusters; break;
458  }
459 
460  event.nclusters[layer] += nclusters;
461  }
462  }
463  }
464 
465  // store
466  m_container->addEvent(event);
467 }
468 
469 //_____________________________________________________________________
471 {
472 
474 
475  // clear array
477 
478  // first loop over hitsets
479  for( const auto& [hitsetkey,hitset]:range_adaptor(m_hitsetcontainer->getHitSets()))
480  {
481  for( const auto& [key,cluster]:range_adaptor(m_cluster_map->getClusters(hitsetkey)))
482  {
483  // create cluster structure
484  auto cluster_struct = create_cluster( key, cluster );
485  add_cluster_size( cluster_struct, key, m_cluster_hit_map );
486  add_cluster_energy( cluster_struct, key, m_cluster_hit_map, m_hitsetcontainer );
487 
488  // truth information
489  const auto g4hits = find_g4hits( key );
490  const bool is_micromegas( TrkrDefs::getTrkrId(key) == TrkrDefs::micromegasId );
491  if( is_micromegas )
492  {
493  const int tileid = MicromegasDefs::getTileId(key);
494  add_truth_information_micromegas( cluster_struct, tileid, g4hits );
495  } else {
496  add_truth_information( cluster_struct, g4hits );
497  }
498 
499  // add in array
500  m_container->addCluster( cluster_struct );
501  }
502  }
503 
504 }
505 
506 //_____________________________________________________________________
508 {
509  if( !( m_track_map && m_cluster_map && m_container ) ) return;
510 
511  // clear array
513 
514  for( const auto& trackpair:*m_track_map )
515  {
516 
517  const auto track = trackpair.second;
518  auto track_struct = create_track( track );
519 
520  // truth information
521  const auto [id,contributors] = get_max_contributor( track );
522  track_struct.contributors = contributors;
523 
524  // get particle
525  auto particle = m_g4truthinfo->GetParticle(id);
526  track_struct.embed = get_embed(particle);
527  ::add_truth_information(track_struct, particle);
528 
529  // running iterator over track states, used to match a given cluster to a track state
530  auto state_iter = track->begin_states();
531 
532  // loop over clusters
533  for( auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
534  {
535 
536  const auto& cluster_key = *key_iter;
537  auto cluster = m_cluster_map->findCluster( cluster_key );
538  if( !cluster )
539  {
540  std::cout << "TrackEvaluation::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl;
541  continue;
542  }
543 
544  // create new cluster struct
545  auto cluster_struct = create_cluster( cluster_key, cluster );
546  add_cluster_size( cluster_struct, cluster_key, m_cluster_hit_map );
547  add_cluster_energy( cluster_struct, cluster_key, m_cluster_hit_map, m_hitsetcontainer );
548 
549  // truth information
550  const auto g4hits = find_g4hits( cluster_key );
551  const bool is_micromegas( TrkrDefs::getTrkrId(cluster_key) == TrkrDefs::micromegasId );
552  if( is_micromegas )
553  {
554  const int tileid = MicromegasDefs::getTileId(cluster_key);
555  add_truth_information_micromegas( cluster_struct, tileid, g4hits );
556  } else {
557  add_truth_information( cluster_struct, g4hits );
558  }
559 
560  // find track state that is the closest to cluster
561  /* this assumes that both clusters and states are sorted along r */
562  const auto radius( cluster_struct.r );
563  float dr_min = -1;
564  for( auto iter = state_iter; iter != track->end_states(); ++iter )
565  {
566  const auto dr = std::abs( radius - get_r( iter->second->get_x(), iter->second->get_y() ) );
567  if( dr_min < 0 || dr < dr_min )
568  {
569  state_iter = iter;
570  dr_min = dr;
571  } else break;
572  }
573 
574  // store track state in cluster struct
575  if( is_micromegas )
576  {
577  const int tileid = MicromegasDefs::getTileId(cluster_key);
578  add_trk_information_micromegas( cluster_struct, tileid, state_iter->second );
579  } else {
580  add_trk_information( cluster_struct, state_iter->second );
581  }
582 
583  // add to track
584  track_struct.clusters.push_back( cluster_struct );
585  }
586  m_container->addTrack( track_struct );
587  }
588 }
589 
590 //_____________________________________________________________________
592 {
593 
594  // check maps
595  if( !( m_cluster_hit_map && m_hit_truth_map ) ) return G4HitSet();
596 
597  // check if in map
598  auto map_iter = m_g4hit_map.lower_bound( cluster_key );
599  if( map_iter != m_g4hit_map.end() && cluster_key == map_iter->first )
600  { return map_iter->second; }
601 
602  // find hitset associated to cluster
603  G4HitSet out;
604  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
605  for(const auto& [first,hit_key]:range_adaptor( m_cluster_hit_map->getHits(cluster_key)))
606  {
607 
608  // store hits to g4hit associations
609  TrkrHitTruthAssoc::MMap g4hit_map;
610  m_hit_truth_map->getG4Hits( hitset_key, hit_key, g4hit_map );
611 
612  // find corresponding g4 hist
613  for( auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter )
614  {
615 
616  const auto g4hit_key = truth_iter->second.second;
617  PHG4Hit* g4hit = nullptr;
618 
619  switch( TrkrDefs::getTrkrId( hitset_key ) )
620  {
621  case TrkrDefs::mvtxId:
622  if( m_g4hits_mvtx ) g4hit = m_g4hits_mvtx->findHit( g4hit_key );
623  break;
624 
625  case TrkrDefs::inttId:
626  if( m_g4hits_intt ) g4hit = m_g4hits_intt->findHit( g4hit_key );
627  break;
628 
629  case TrkrDefs::tpcId:
630  if( m_g4hits_tpc ) g4hit = m_g4hits_tpc->findHit( g4hit_key );
631  break;
632 
634  if( m_g4hits_micromegas ) g4hit = m_g4hits_micromegas->findHit( g4hit_key );
635  break;
636 
637  default: break;
638  }
639 
640  if( g4hit ) out.insert( g4hit );
641  else std::cout << "TrackEvaluation::find_g4hits - g4hit not found " << g4hit_key << std::endl;
642 
643  }
644  }
645 
646  // insert in map and return
647  return m_g4hit_map.insert( map_iter, std::make_pair( cluster_key, std::move( out ) ) )->second;
648 
649 }
650 
651 //_____________________________________________________________________
652 std::pair<int,int> TrackEvaluation::get_max_contributor( SvtxTrack* track ) const
653 {
654  if(!(m_track_map && m_cluster_map && m_g4truthinfo)) return {0,0};
655 
656  // maps MC track id and number of matching g4hits
657  using IdMap = std::map<int,int>;
658  IdMap contributor_map;
659 
660  // loop over clusters
661  for( auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
662  {
663  const auto& cluster_key = *key_iter;
664  for( const auto& hit:find_g4hits( cluster_key ) )
665  {
666  const int trkid = hit->get_trkid();
667  auto iter = contributor_map.lower_bound( trkid );
668  if( iter == contributor_map.end() || iter->first != trkid )
669  {
670  contributor_map.insert(iter, std::make_pair(trkid,1));
671  } else ++iter->second;
672  }
673  }
674 
675  if( contributor_map.empty() ) return {0,0};
676  else return *std::max_element(
677  contributor_map.cbegin(), contributor_map.cend(),
678  []( const IdMap::value_type& first, const IdMap::value_type& second )
679  { return first.second < second.second; } );
680 
681 }
682 
683 //_____________________________________________________________________
685 { return (m_g4truthinfo && particle) ? m_g4truthinfo->isEmbeded( particle->get_primary_id() ):0; }
686 
687 //_____________________________________________________________________
689 {
690  // get global coordinates
691  const auto global = m_transformer.getGlobalPosition(cluster,m_surfmaps, m_tGeometry);
692 
694  cluster_struct.layer = TrkrDefs::getLayer(key);
695  cluster_struct.x = global.x();
696  cluster_struct.y = global.y();
697  cluster_struct.z = global.z();
698  cluster_struct.r = get_r( cluster_struct.x, cluster_struct.y );
699  cluster_struct.phi = std::atan2( cluster_struct.y, cluster_struct.x );
700  cluster_struct.phi_error = cluster->getRPhiError()/cluster_struct.r;
701  cluster_struct.z_error = cluster->getZError();
702 
703  return cluster_struct;
704 }
705 
706 //_____________________________________________________________________
708 {
709  // need to extrapolate to the right r
710  const auto trk_r = get_r( state->get_x(), state->get_y() );
711  const auto dr = cluster.r - trk_r;
712  const auto trk_drdt = (state->get_x()*state->get_px() + state->get_y()*state->get_py())/trk_r;
713  const auto trk_dxdr = state->get_px()/trk_drdt;
714  const auto trk_dydr = state->get_py()/trk_drdt;
715  const auto trk_dzdr = state->get_pz()/trk_drdt;
716 
717  // store state position
718  cluster.trk_x = state->get_x() + dr*trk_dxdr;
719  cluster.trk_y = state->get_y() + dr*trk_dydr;
720  cluster.trk_z = state->get_z() + dr*trk_dzdr;
721  cluster.trk_r = get_r( cluster.trk_x, cluster.trk_y );
722  cluster.trk_phi = std::atan2( cluster.trk_y, cluster.trk_x );
723 
724  /* store local momentum information */
725  cluster.trk_px = state->get_px();
726  cluster.trk_py = state->get_py();
727  cluster.trk_pz = state->get_pz();
728 
729  /*
730  store state angles in (r,phi) and (r,z) plans
731  they are needed to study space charge distortions
732  */
733  const auto cosphi( std::cos( cluster.trk_phi ) );
734  const auto sinphi( std::sin( cluster.trk_phi ) );
735  const auto trk_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
736  const auto trk_pr = state->get_px()*cosphi + state->get_py()*sinphi;
737  const auto trk_pz = state->get_pz();
738  cluster.trk_alpha = std::atan2( trk_pphi, trk_pr );
739  cluster.trk_beta = std::atan2( trk_pz, trk_pr );
740  cluster.trk_phi_error = state->get_phi_error();
741  cluster.trk_z_error = state->get_z_error();
742 
743 }
744 
745 //_____________________________________________________________________
747 {
748 
749  // get geometry cylinder from layer
750  const auto layer = cluster.layer;
751  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(m_micromegas_geom_container->GetLayerGeom(layer));
752  assert( layergeom );
753 
754  // convert cluster position to local tile coordinates
755  const TVector3 cluster_world( cluster.x, cluster.y, cluster.z );
756  const TVector3 cluster_local = layergeom->get_local_from_world_coords( tileid, cluster_world );
757 
758  // convert track position to local tile coordinates
759  TVector3 track_world( state->get_x(), state->get_y(), state->get_z() );
760  TVector3 track_local = layergeom->get_local_from_world_coords( tileid, track_world );
761 
762  // convert direction to local tile coordinates
763  const TVector3 direction_world( state->get_px(), state->get_py(), state->get_pz() );
764  const TVector3 direction_local = layergeom->get_local_from_world_vect( tileid, direction_world );
765 
766  // extrapolate to same local y (should be zero) as cluster
767  const auto delta_y = cluster_local.y() - track_local.y();
768  track_local += TVector3(
769  delta_y*direction_local.x()/direction_local.y(),
770  delta_y,
771  delta_y*direction_local.z()/direction_local.y() );
772 
773  // convert back to global coordinates
774  track_world = layergeom->get_world_from_local_coords( tileid, track_local );
775 
776  // store state position
777  cluster.trk_x = track_world.x();
778  cluster.trk_y = track_world.y();
779  cluster.trk_z = track_world.z();
780  cluster.trk_r = get_r( cluster.trk_x, cluster.trk_y );
781  cluster.trk_phi = std::atan2( cluster.trk_y, cluster.trk_x );
782 
783  /* store local momentum information */
784  cluster.trk_px = state->get_px();
785  cluster.trk_py = state->get_py();
786  cluster.trk_pz = state->get_pz();
787 
788  /*
789  store state angles in (r,phi) and (r,z) plans
790  they are needed to study space charge distortions
791  */
792  const auto cosphi( std::cos( cluster.trk_phi ) );
793  const auto sinphi( std::sin( cluster.trk_phi ) );
794  const auto trk_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
795  const auto trk_pr = state->get_px()*cosphi + state->get_py()*sinphi;
796  const auto trk_pz = state->get_pz();
797  cluster.trk_alpha = std::atan2( trk_pphi, trk_pr );
798  cluster.trk_beta = std::atan2( trk_pz, trk_pr );
799  cluster.trk_phi_error = state->get_phi_error();
800  cluster.trk_z_error = state->get_z_error();
801 
802 }
803 
804 //_____________________________________________________________________
806 {
807  // store number of contributing g4hits
808  cluster.truth_size = g4hits.size();
809 
810  // get layer, tpc flag and corresponding layer geometry
811  const auto layer = cluster.layer;
812  const bool is_tpc( layer >= 7 && layer < 55 );
813  const PHG4CylinderCellGeom* layergeom = is_tpc ? m_tpc_geom_container->GetLayerCellGeom(layer):nullptr;
814  const auto rin = layergeom ? layergeom->get_radius()-layergeom->get_thickness()/2:0;
815  const auto rout = layergeom ? layergeom->get_radius()+layergeom->get_thickness()/2:0;
816 
817  // convert hits to list of interpolation_data_t
818  interpolation_data_t::list hits;
819  for( const auto& g4hit:g4hits )
820  {
821  interpolation_data_t::list tmp_hits;
822  const auto weight = g4hit->get_edep();
823  for( int i = 0; i < 2; ++i )
824  {
825  const TVector3 g4hit_world(g4hit->get_x(i), g4hit->get_y(i), g4hit->get_z(i));
826  const TVector3 momentum_world(g4hit->get_px(i), g4hit->get_py(i), g4hit->get_pz(i));
827  tmp_hits.push_back( {.position = g4hit_world, .momentum = momentum_world, .weight = weight } );
828  }
829 
830  if( is_tpc )
831  {
832  // add layer boundary checks
833  // ensure first hit has lowest r
834  auto r0 = get_r(tmp_hits[0].x(),tmp_hits[0].y());
835  auto r1 = get_r(tmp_hits[1].x(),tmp_hits[1].y());
836  if( r0 > r1 )
837  {
838  std::swap(tmp_hits[0],tmp_hits[1]);
839  std::swap(r0, r1);
840  }
841 
842  // do nothing if out of bound
843  if( r1 <= rin || r0 >= rout ) continue;
844 
845  // keep track of original deltar
846  const auto dr_old = r1-r0;
847 
848  // clamp r0 to rin
849  if( r0 < rin )
850  {
851  const auto t = line_circle_intersection( tmp_hits[0].position, tmp_hits[1].position, rin );
852  if( t<0 ) continue;
853 
854  tmp_hits[0].position = tmp_hits[0].position*(1.-t) + tmp_hits[1].position*t;
855  tmp_hits[0].momentum = tmp_hits[0].momentum*(1.-t) + tmp_hits[1].momentum*t;
856  r0 = rin;
857  }
858 
859  if( r1 > rout )
860  {
861  const auto t = line_circle_intersection( tmp_hits[0].position, tmp_hits[1].position, rout );
862  if( t<0 ) continue;
863 
864  tmp_hits[1].position = tmp_hits[0].position*(1.-t) + tmp_hits[1].position*t;
865  tmp_hits[1].momentum = tmp_hits[0].momentum*(1.-t) + tmp_hits[1].momentum*t;
866  r1 = rout;
867  }
868 
869  // update weights, only if clamping occured
870  const auto dr_new = r1-r0;
871  tmp_hits[0].weight *= dr_new/dr_old;
872  tmp_hits[1].weight *= dr_new/dr_old;
873  }
874 
875  // store in global list
876  hits.push_back(std::move(tmp_hits[0]));
877  hits.push_back(std::move(tmp_hits[1]));
878  }
879 
880  // add truth position
881  cluster.truth_x = average<&interpolation_data_t::x>( hits );
882  cluster.truth_y = average<&interpolation_data_t::y>( hits );
883  cluster.truth_z = average<&interpolation_data_t::z>( hits );
884 
885  // add truth momentum information
886  cluster.truth_px = average<&interpolation_data_t::px>( hits );
887  cluster.truth_py = average<&interpolation_data_t::py>( hits );
888  cluster.truth_pz = average<&interpolation_data_t::pz>( hits );
889 
890  cluster.truth_r = get_r( cluster.truth_x, cluster.truth_y );
891  cluster.truth_phi = std::atan2( cluster.truth_y, cluster.truth_x );
892 
893  /*
894  store state angles in (r,phi) and (r,z) plans
895  they are needed to study space charge distortions
896  */
897  const auto cosphi( std::cos( cluster.truth_phi ) );
898  const auto sinphi( std::sin( cluster.truth_phi ) );
899  const auto truth_pphi = -cluster.truth_px*sinphi + cluster.truth_py*cosphi;
900  const auto truth_pr = cluster.truth_px*cosphi + cluster.truth_py*sinphi;
901 
902  cluster.truth_alpha = std::atan2( truth_pphi, truth_pr );
903  cluster.truth_beta = std::atan2( cluster.truth_pz, truth_pr );
904 
905 }
906 
907 //_____________________________________________________________________
908 void TrackEvaluation::add_truth_information_micromegas( TrackEvaluationContainerv1::ClusterStruct& cluster, int tileid, std::set<PHG4Hit*> g4hits ) const
909 {
910  // store number of contributing g4hits
911  cluster.truth_size = g4hits.size();
912 
913  const auto layer = cluster.layer;
914  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(m_micromegas_geom_container->GetLayerGeom(layer));
915  assert( layergeom );
916 
917  // convert cluster position to local tile coordinates
918  const TVector3 cluster_world( cluster.x, cluster.y, cluster.z );
919  const TVector3 cluster_local = layergeom->get_local_from_world_coords( tileid, cluster_world );
920 
921  // convert hits to list of interpolation_data_t
922  interpolation_data_t::list hits;
923  for( const auto& g4hit:g4hits )
924  {
925  const auto weight = g4hit->get_edep();
926 
927  // convert g4hit to planar coordinate
928  PHG4Hitv1 g4hit_copy( g4hit );
929  layergeom->convert_to_planar( tileid, &g4hit_copy );
930 
931  for( int i = 0; i < 2; ++i )
932  {
933 
934  // convert position to local
935  TVector3 g4hit_world(g4hit_copy.get_x(i), g4hit_copy.get_y(i), g4hit_copy.get_z(i));
936  TVector3 g4hit_local = layergeom->get_local_from_world_coords( tileid, g4hit_world );
937 
938  // convert momentum to local
939  TVector3 momentum_world(g4hit_copy.get_px(i), g4hit_copy.get_py(i), g4hit_copy.get_pz(i));
940  TVector3 momentum_local = layergeom->get_local_from_world_vect( tileid, momentum_world );
941 
942  hits.push_back( {.position = g4hit_local, .momentum = momentum_local, .weight = weight } );
943  }
944  }
945 
946  // do position interpolation
947  const TVector3 interpolation_local(
948  average<&interpolation_data_t::x>(hits),
949  average<&interpolation_data_t::y>(hits),
950  average<&interpolation_data_t::z>(hits) );
951 
952  const TVector3 interpolation_world = layergeom->get_world_from_local_coords( tileid, interpolation_local );
953 
954  // do momentum interpolation
955  const TVector3 momentum_local(
956  average<&interpolation_data_t::px>(hits),
957  average<&interpolation_data_t::py>(hits),
958  average<&interpolation_data_t::pz>(hits));
959 
960  const TVector3 momentum_world = layergeom->get_world_from_local_vect( tileid, momentum_local );
961 
962  cluster.truth_x = interpolation_world.x();
963  cluster.truth_y = interpolation_world.y();
964  cluster.truth_z = interpolation_world.z();
965  cluster.truth_r = get_r( cluster.truth_x, cluster.truth_y );
966  cluster.truth_phi = std::atan2( cluster.truth_y, cluster.truth_x );
967 
968  /* add truth momentum information */
969  cluster.truth_px = momentum_world.x();
970  cluster.truth_py = momentum_world.y();
971  cluster.truth_pz = momentum_world.z();
972 
973  /*
974  store state angles in (r,phi) and (r,z) plans
975  they are needed to study space charge distortions
976  */
977  const auto cosphi( std::cos( cluster.truth_phi ) );
978  const auto sinphi( std::sin( cluster.truth_phi ) );
979  const auto truth_pphi = -cluster.truth_px*sinphi + cluster.truth_py*cosphi;
980  const auto truth_pr = cluster.truth_px*cosphi + cluster.truth_py*sinphi;
981 
982  cluster.truth_alpha = std::atan2( truth_pphi, truth_pr );
983  cluster.truth_beta = std::atan2( cluster.truth_pz, truth_pr );
984 
985 }