ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DSTEmulator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DSTEmulator.cc
1 
6 #include "DSTEmulator.h"
8 #include "DSTCompressor.h"
9 
11 #include <g4main/PHG4Hit.h>
13 #include <g4main/PHG4Particle.h>
15 #include <intt/InttDefs.h>
17 #include <mvtx/MvtxDefs.h>
18 #include <phool/getClass.h>
19 #include <phool/PHCompositeNode.h>
20 #include <phool/PHNodeIterator.h>
21 #include <phool/recoConsts.h>
22 #include <tpc/TpcDefs.h>
23 #include <trackbase/TrkrDefs.h>
27 #include <trackbase/TrkrHit.h>
28 #include <trackbase/TrkrHitSet.h>
33 
34 #include <Acts/Utilities/Units.hpp>
36 
37 #include <TFile.h>
38 #include <TNtuple.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( 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 
76  /* //! radius
77  float get_r( PHG4Hit* hit, int i )
78  { return get_r( hit->get_x(i), hit->get_y(i) ); }
79  */
80  /*
82  template< float (PHG4Hit::*accessor)(int) const>
83  float interpolate( std::set<PHG4Hit*> hits, float rextrap )
84  {
85  // calculate all terms needed for the interpolation
86  // need to use double everywhere here due to numerical divergences
87  double sw = 0;
88  double swr = 0;
89  double swr2 = 0;
90  double swx = 0;
91  double swrx = 0;
92 
93  bool valid( false );
94  for( const auto& hit:hits )
95  {
96 
97  const double x0 = (hit->*accessor)(0);
98  const double x1 = (hit->*accessor)(1);
99  if( std::isnan( x0 ) || std::isnan( x1 ) ) continue;
100 
101  const double w = hit->get_edep();
102  if( w < 0 ) continue;
103 
104  valid = true;
105  const double r0 = get_r( hit, 0 );
106  const double r1 = get_r( hit, 1 );
107 
108  sw += w*2;
109  swr += w*(r0 + r1);
110  swr2 += w*(square(r0) + square(r1));
111  swx += w*(x0 + x1);
112  swrx += w*(r0*x0 + r1*x1);
113  }
114 
115  if( !valid ) return NAN;
116 
117  const auto alpha = (sw*swrx - swr*swx);
118  const auto beta = (swr2*swx - swr*swrx);
119  const auto denom = (sw*swr2 - square(swr));
120 
121  return ( alpha*rextrap + beta )/denom;
122  }
123  */
124  /*
126  inline int is_primary( PHG4Particle* particle )
127  { return particle->get_parent_id() == 0; }
128  */
130  int64_t get_mask( SvtxTrack* track )
131  { return std::accumulate( track->begin_cluster_keys(), track->end_cluster_keys(), int64_t(0),
132  []( int64_t value, const TrkrDefs::cluskey& key ) {
133  return TrkrDefs::getLayer(key)<64 ? value|(1LL<<TrkrDefs::getLayer(key)) : value;
134  } );
135  }
136 
138  template<int type>
139  int get_clusters( SvtxTrack* track )
140  {
141  return std::count_if( track->begin_cluster_keys(), track->end_cluster_keys(),
142  []( const TrkrDefs::cluskey& key ) { return TrkrDefs::getTrkrId(key) == type; } );
143  }
144 
147  {
149 
150  trackStruct.charge = track->get_charge();
151  trackStruct.nclusters = track->size_cluster_keys();
152  trackStruct.mask = get_mask( track );
153  trackStruct.nclusters_mvtx = get_clusters<TrkrDefs::mvtxId>( track );
154  trackStruct.nclusters_intt = get_clusters<TrkrDefs::inttId>( track );
155  trackStruct.nclusters_tpc = get_clusters<TrkrDefs::tpcId>( track );
156  trackStruct.nclusters_micromegas = get_clusters<TrkrDefs::micromegasId>( track );
157 
158  trackStruct.chisquare = track->get_chisq();
159  trackStruct.ndf = track->get_ndf();
160 
161  trackStruct.x = track->get_x();
162  trackStruct.y = track->get_y();
163  trackStruct.z = track->get_z();
164  trackStruct.r = get_r( trackStruct.x, trackStruct.y );
165  trackStruct.phi = std::atan2( trackStruct.y, trackStruct.x );
166 
167  trackStruct.px = track->get_px();
168  trackStruct.py = track->get_py();
169  trackStruct.pz = track->get_pz();
170  trackStruct.pt = get_pt( trackStruct.px, trackStruct.py );
171  trackStruct.p = get_p( trackStruct.px, trackStruct.py, trackStruct.pz );
172  trackStruct.eta = get_eta( trackStruct.p, trackStruct.pz );
173 
174  return trackStruct;
175  }
176 #ifdef JUNK
177 
179  {
181  cluster_struct.layer = TrkrDefs::getLayer(key);
182  cluster_struct.x = cluster->getX();
183  cluster_struct.y = cluster->getY();
184  cluster_struct.z = cluster->getZ();
185  cluster_struct.r = get_r( cluster_struct.x, cluster_struct.y );
186  cluster_struct.phi = std::atan2( cluster_struct.y, cluster_struct.x );
187  cluster_struct.phi_error = cluster->getPhiError();
188  cluster_struct.z_error = cluster->getZError();
189  std::cout << " (x|y|z|r|l) "
190  << cluster_struct.x << " | "
191  << cluster_struct.y << " | "
192  << cluster_struct.z << " | "
193  << cluster_struct.r << " | "
194  << cluster_struct.layer << " | "
195  << std::endl;
196  std::cout << " (xl|yl) "
197  << cluster->getLocalX() << " | "
198  << cluster->getLocalY()
199  << std::endl;
200  return cluster_struct;
201  }
202 
204  void add_trk_information( TrackEvaluationContainerv1::ClusterStruct& cluster, SvtxTrackState* state )
205  {
206  // need to extrapolate the track state to the right cluster radius to get proper residuals
207  const auto trk_r = get_r( state->get_x(), state->get_y() );
208  const auto dr = cluster.r - trk_r;
209  const auto trk_drdt = (state->get_x()*state->get_px() + state->get_y()*state->get_py())/trk_r;
210  const auto trk_dxdr = state->get_px()/trk_drdt;
211  const auto trk_dydr = state->get_py()/trk_drdt;
212  const auto trk_dzdr = state->get_pz()/trk_drdt;
213 
214  // store state position
215  cluster.trk_x = state->get_x() + dr*trk_dxdr;
216  cluster.trk_y = state->get_y() + dr*trk_dydr;
217  cluster.trk_z = state->get_z() + dr*trk_dzdr;
218  cluster.trk_r = get_r( cluster.trk_x, cluster.trk_y );
219  cluster.trk_phi = std::atan2( cluster.trk_y, cluster.trk_x );
220 
221  /* store local momentum information */
222  cluster.trk_px = state->get_px();
223  cluster.trk_py = state->get_py();
224  cluster.trk_pz = state->get_pz();
225 
226  /*
227  store state angles in (r,phi) and (r,z) plans
228  they are needed to study space charge distortions
229  */
230  const auto cosphi( std::cos( cluster.trk_phi ) );
231  const auto sinphi( std::sin( cluster.trk_phi ) );
232  const auto trk_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
233  const auto trk_pr = state->get_px()*cosphi + state->get_py()*sinphi;
234  const auto trk_pz = state->get_pz();
235  cluster.trk_alpha = std::atan2( trk_pphi, trk_pr );
236  cluster.trk_beta = std::atan2( trk_pz, trk_pr );
237  cluster.trk_phi_error = state->get_phi_error();
238  cluster.trk_z_error = state->get_z_error();
239 
240  }
241 
243  void add_cluster_size( TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key, TrkrClusterHitAssoc* cluster_hit_map )
244  {
245  if( !cluster_hit_map ) return;
246  const auto range = cluster_hit_map->getHits(clus_key);
247 
248  // store full size
249  cluster.size = std::distance( range.first, range.second );
250 
251  const auto detId = TrkrDefs::getTrkrId(clus_key);
252  if(detId == TrkrDefs::micromegasId)
253  {
254 
255  // for micromegas the directional cluster size depends on segmentation type
256  auto segmentation_type = MicromegasDefs::getSegmentationType(clus_key);
257  if( segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_Z ) cluster.z_size = cluster.size;
258  else cluster.phi_size = cluster.size;
259 
260  } else {
261 
262  // for other detectors, one must loop over the constituting hits
263  std::set<int> phibins;
264  std::set<int> zbins;
265  for(const auto& [first, hit_key]:range_adaptor(range))
266  {
267  switch( detId )
268  {
269  default: break;
270  case TrkrDefs::mvtxId:
271  {
272  phibins.insert( MvtxDefs::getRow( hit_key ) );
273  zbins.insert( MvtxDefs::getCol( hit_key ) );
274  break;
275  }
276  case TrkrDefs::inttId:
277  {
278  phibins.insert( InttDefs::getRow( hit_key ) );
279  zbins.insert( InttDefs::getCol( hit_key ) );
280  break;
281  }
282  case TrkrDefs::tpcId:
283  {
284  phibins.insert( TpcDefs::getPad( hit_key ) );
285  zbins.insert( TpcDefs::getTBin( hit_key ) );
286  break;
287  }
288  }
289  }
290  cluster.phi_size = phibins.size();
291  cluster.z_size = zbins.size();
292  }
293  }
294 
295 
297  void add_cluster_energy( TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key,
298  TrkrClusterHitAssoc* cluster_hit_map,
299  TrkrHitSetContainer* hitsetcontainer )
300  {
301 
302  // check container
303  if(!(cluster_hit_map && hitsetcontainer)) return;
304 
305  // for now this is only filled for micromegas
306  const auto detId = TrkrDefs::getTrkrId(clus_key);
307  if(detId != TrkrDefs::micromegasId) return;
308 
309  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(clus_key);
310  const auto hitset = hitsetcontainer->findHitSet( hitset_key );
311  if( !hitset ) return;
312 
313  const auto range = cluster_hit_map->getHits(clus_key);
314  cluster.energy_max = 0;
315  cluster.energy_sum = 0;
316 
317  for( const auto& pair:range_adaptor(range))
318  {
319  const auto hit = hitset->getHit( pair.second );
320  if( hit )
321  {
322  const auto energy = hit->getEnergy();
323  cluster.energy_sum += energy;
324  if( energy > cluster.energy_max ) cluster.energy_max = energy;
325  }
326  }
327 
328  }
329 
330  // add truth information
331  void add_truth_information( TrackEvaluationContainerv1::ClusterStruct& cluster, std::set<PHG4Hit*> hits )
332  {
333  // need to extrapolate g4hits position to the cluster r
334  /* this is done by using a linear extrapolation, with a straight line going through all provided G4Hits in and out positions */
335  const auto rextrap = cluster.r;
336  cluster.truth_size = hits.size();
337  std::cout << "inter" << "\n";
338 
339  cluster.truth_x = interpolate<&PHG4Hit::get_x>( hits, rextrap );
340  cluster.truth_y = interpolate<&PHG4Hit::get_y>( hits, rextrap );
341  cluster.truth_z = interpolate<&PHG4Hit::get_z>( hits, rextrap );
342  cluster.truth_r = get_r( cluster.truth_x, cluster.truth_y );
343  cluster.truth_phi = std::atan2( cluster.truth_y, cluster.truth_x );
344 
345  /* add truth momentum information */
346  cluster.truth_px = interpolate<&PHG4Hit::get_px>( hits, rextrap );
347  cluster.truth_py = interpolate<&PHG4Hit::get_py>( hits, rextrap );
348  cluster.truth_pz = interpolate<&PHG4Hit::get_pz>( hits, rextrap );
349 
350  std::cout << "inter2" << "\n";
351 
352  /*
353  store state angles in (r,phi) and (r,z) plans
354  they are needed to study space charge distortions
355  */
356  const auto cosphi( std::cos( cluster.truth_phi ) );
357  const auto sinphi( std::sin( cluster.truth_phi ) );
358  const auto truth_pphi = -cluster.truth_px*sinphi + cluster.truth_py*cosphi;
359  const auto truth_pr = cluster.truth_px*cosphi + cluster.truth_py*sinphi;
360 
361  cluster.truth_alpha = std::atan2( truth_pphi, truth_pr );
362  cluster.truth_beta = std::atan2( cluster.truth_pz, truth_pr );
363  if(std::isnan(cluster.truth_alpha) || std::isnan(cluster.truth_beta))
364  {
365  // recalculate
366  double truth_alpha = 0;
367  double truth_beta = 0;
368  double sum_w = 0;
369  for( const auto& hit:hits )
370  {
371  const auto px = hit->get_x(1) - hit->get_x(0);
372  const auto py = hit->get_y(1) - hit->get_y(0);
373  const auto pz = hit->get_z(1) - hit->get_z(0);
374  const auto pphi = -px*sinphi + py*cosphi;
375  const auto pr = px*cosphi + py*sinphi;
376 
377  const auto w = hit->get_edep();
378  if( w < 0 ) continue;
379 
380  sum_w += w;
381  truth_alpha += w*std::atan2( pphi, pr );
382  truth_beta += w*std::atan2( pz, pr );
383  }
384  truth_alpha /= sum_w;
385  truth_beta /= sum_w;
386  cluster.truth_alpha = truth_alpha;
387  cluster.truth_beta = truth_beta;
388  }
389 
390  }
391 
392  // ad}d truth information
393  void add_truth_information( TrackEvaluationContainerv1::TrackStruct& track, PHG4Particle* particle )
394  {
395  if( particle )
396  {
397  track.is_primary = is_primary( particle );
398  track.pid = particle->get_pid();
399  track.truth_px = particle->get_px();
400  track.truth_py = particle->get_py();
401  track.truth_pz = particle->get_pz();
402  track.truth_pt = get_pt( track.truth_px, track.truth_py );
403  track.truth_p = get_p( track.truth_px, track.truth_py, track.truth_pz );
404  track.truth_eta = get_eta( track.truth_p, track.truth_pz );
405  }
406  }
407 #endif
408 }
409 
410 //_____________________________________________________________________
411 DSTEmulator::DSTEmulator( const std::string& name, const std::string &filename, int inBits,
412  int inSabotage, bool compress):
413  SubsysReco( name)
414  , _filename(filename)
415  , _tfile(nullptr)
416  , nBits(inBits)
417  , sabotage(inSabotage)
418  , apply_compression(compress)
419 {}
420 
421 //_____________________________________________________________________
423 {
424  if (_tfile) delete _tfile;
425  _tfile = new TFile(_filename.c_str(), "RECREATE");
426 
427  _dst_data = new TNtuple("dst_data", "dst data","event:seed:"
428  "c3x:c3y:c3z:c3p:t3x:t3y:t3z:"
429  "c2x:c2y:c2r:c2l:t2x:t2y:t2r:t2l:"
430  "d2x:d2y:dr:"
431  "cmp_d2x:cmp_d2y:"
432  "pt:eta:phi:charge");
433 
434  // find DST node
435  PHNodeIterator iter(topNode);
436  auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
437  if (!dstNode)
438  {
439  std::cout << "DSTEmulator::Init - DST Node missing" << std::endl;
441  }
442 
443  // get EVAL node
444  iter = PHNodeIterator(dstNode);
445  auto evalNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "EVAL"));
446  if( !evalNode )
447  {
448  // create
449  std::cout << "DSTEmulator::Init - EVAL node missing - creating" << std::endl;
450  evalNode = new PHCompositeNode( "EVAL" );
451  dstNode->addNode(evalNode);
452  }
453 
454  auto newNode = new PHIODataNode<PHObject>( new TrackEvaluationContainerv1, "TrackEvaluationContainer","PHObject");
455  evalNode->addNode(newNode);
456 
457  // m_compressor = new DSTCompressor(4.08407e-02,
458  // 7.46530e-02,
459  // 5.14381e-05,
460  // 2.06291e-01,
461  // with new distortion setting
462  // m_compressor = new DSTCompressor(-2.70072e-02,
463  // 2.49574e-02,
464  // 1.12803e-03,
465  // 5.91965e-02,
466  // with mininum layer 7
467  m_compressor = new DSTCompressor(6.96257e-04,
468  3.16806e-02,
469  7.32860e-05,
470  5.93230e-02,
471  nBits,
472  nBits);
474 }
475 
476 //_____________________________________________________________________
479 
480 //_____________________________________________________________________
482 {
483  // load nodes
484  auto res = load_nodes(topNode);
485  if( res != Fun4AllReturnCodes::EVENT_OK ) return res;
486 
487  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
488  "ActsTrackingGeometry");
489  if(!m_tGeometry)
490  {
491  std::cout << PHWHERE
492  << "ActsTrackingGeometry not found on node tree. Exiting"
493  << std::endl;
495  }
496 
497  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
498  "ActsSurfaceMaps");
499  if(!m_surfMaps)
500  {
501  std::cout << PHWHERE
502  << "ActsSurfaceMaps not found on node tree. Exiting"
503  << std::endl;
505  }
506 
507  // cleanup output container
508  if( m_container ) m_container->Reset();
509 
510  evaluate_tracks();
511 
512  // clear maps
513  m_g4hit_map.clear();
515 }
516 
517 //_____________________________________________________________________
519 {
520  _tfile->cd();
521  _dst_data->Write();
522  _tfile->Close();
523  delete _tfile;
525 }
526 
528 {
529  // get global position from Acts transform
530  auto globalpos = m_transform.getGlobalPosition(cluster,
531  m_surfMaps,
532  m_tGeometry);
533 
534  // check if TPC distortion correction are in place and apply
535  // if( m_dcc ) { globalpos = m_distortionCorrection.get_corrected_position( globalpos, m_dcc ); }
536 
537  return globalpos;
538 }
539 
540 //_____________________________________________________________________
542 {
543 
544  // get necessary nodes
545  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMapTPCOnly");
546  if(m_track_map)
547  {
548  std::cout << " DSTEmulator: Using TPC Only Track Map node " << std::endl;
549  }
550  else
551  {
552  std::cout << " DSTEmulator: TPC Only Track Map node not found, using default" << std::endl;
553  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
554  }
555  // cluster map
556 
557  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"CORRECTED_TRKR_CLUSTER");
558  if(m_cluster_map){
559  if(m_cluster_map->size()>0)
560  {
561  std::cout << " DSTEmulator: Using CORRECTED_TRKR_CLUSTER node " << std::endl;
562  }
563  else
564  {
565  std::cout << " DSTEmulator: CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER" << std::endl;
566  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
567  }
568  }else{
569 
570  std::cout << " DSTEmulator: CORRECTED_TRKR_CLUSTER node not found at all, using TRKR_CLUSTER" << std::endl;
571  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
572  }
573 
574 
575  // cluster hit association map
576  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
577 
578  // cluster hit association map
579  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,"TRKR_HITTRUTHASSOC");
580 
581  // local container
582  m_container = findNode::getClass<TrackEvaluationContainerv1>(topNode, "TrackEvaluationContainer");
583 
584  // hitset container
585  m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
586 
587  // g4hits
588  m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
589  m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
590  m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
591  m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
592 
593  // g4 truth info
594  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
595 
597 
598 }
599 
600 //_____________________________________________________________________
602 {
603 
604  if( !( m_track_map && m_cluster_map && m_container ) ) return;
605 
606  int iseed = 0;
608  if (rc->FlagExist("RANDOMSEED")){
609  iseed = rc->get_IntFlag("RANDOMSEED");
610  }
611 
612  // clear array
614 
615  for( const auto& trackpair:*m_track_map )
616  {
617 
618  const auto track = trackpair.second;
619  auto track_struct = create_track( track );
620 
621  // truth information
622  const auto [id,contributors] = get_max_contributor( track );
623  track_struct.contributors = contributors;
624 
625  // get particle
626  auto particle = m_g4truthinfo->GetParticle(id);
627  track_struct.embed = get_embed(particle);
628  // add_truth_information(track_struct, particle);
629 
630  // running iterator over track states, used to match a given cluster to a track state
631  auto state_iter = track->begin_states();
632 
633  // loop over clusters
634  for( auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
635  {
636 
637  const auto& cluster_key = *key_iter;
638  auto cluster = m_cluster_map->findCluster( cluster_key );
640 
641  if( !cluster )
642  {
643  std::cout << "DSTEmulator::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl;
644  continue;
645  }
646 
647  if(TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::tpcId) continue;
648 
649  // create new cluster struct
650  // auto cluster_struct = create_cluster( cluster_key, cluster );
651 
652  // find track state that is the closest to cluster
653  /* this assumes that both clusters and states are sorted along r */
654  // const auto radius( cluster_struct.r );
655  const Acts::Vector3D globalpos_d = getGlobalPosition(cluster);
656  float radius = get_r( globalpos_d[0], globalpos_d[1] );
657  float clu_phi = std::atan2( globalpos_d[0], globalpos_d[1] );
658  std::cout << "radius " << radius << std::endl;
659  float dr_min = -1;
660  for( auto iter = state_iter; iter != track->end_states(); ++iter )
661  {
662  const auto dr = std::abs( radius - get_r( iter->second->get_x(), iter->second->get_y() ) );
663  if( dr_min < 0 || dr < dr_min )
664  {
665  state_iter = iter;
666  dr_min = dr;
667  } else break;
668  }
669  // Got cluster and got state
670  /*
671  */
672 
673  std::cout << "NEW (xg|yg) "
674  << globalpos_d[0] << " | "
675  << globalpos_d[1]
676  << std::endl;
677  std::cout << "NEW (xl|yl) "
678  << cluster->getLocalX() << " | "
679  << cluster->getLocalY()
680  << std::endl;
681 
682  //state to local
683  // store track state in cluster struct
684  // add_trk_information( cluster_struct, state_iter->second );
685  // void add_trk_information( TrackEvaluationContainerv1::ClusterStruct& cluster, SvtxTrackState* state )
686  // need to extrapolate the track state to the right cluster radius to get proper residuals
687  const auto trk_r = get_r( state_iter->second->get_x(), state_iter->second->get_y() );
688  std::cout << " trk_r " << trk_r << std::endl;
689  const auto dr = get_r( globalpos_d[0], globalpos_d[1] ) - trk_r;
690  std::cout << " dr " << dr << std::endl;
691  const auto trk_drdt = (state_iter->second->get_x()*state_iter->second->get_px() + state_iter->second->get_y()*state_iter->second->get_py())/trk_r;
692  std::cout << " trk_drdt " << trk_drdt << std::endl;
693  const auto trk_dxdr = state_iter->second->get_px()/trk_drdt;
694  std::cout << " trk_dxdr " << trk_dxdr << std::endl;
695  const auto trk_dydr = state_iter->second->get_py()/trk_drdt;
696  std::cout << " trk_dydr " << trk_dydr << std::endl;
697  const auto trk_dzdr = state_iter->second->get_pz()/trk_drdt;
698  std::cout << " trk_dzdr " << trk_dzdr << std::endl;
699 
700  // store state position
701  /* */
702  float trk_x = state_iter->second->get_x() + dr*trk_dxdr;
703  float trk_y = state_iter->second->get_y() + dr*trk_dydr;
704  float trk_z = state_iter->second->get_z() + dr*trk_dzdr;
705  std::cout << "trk_x " << state_iter->second->get_x() << "trk_y" << state_iter->second->get_y() << "trk_z " << state_iter->second->get_z() << std::endl;
706  // float trk_r = get_r( trk_x, trk_y );
707  //cluster.trk_phi = std::atan2( cluster.trk_y, cluster.trk_x );
708 
709 
710  /* store local momentum information
711  cluster.trk_px = state->get_px();
712  cluster.trk_py = state->get_py();
713  cluster.trk_pz = state->get_pz();
714  */
715  auto layer = TrkrDefs::getLayer(cluster_key);
716  /*
717  std::cout << " track 3D (x|y|z|r|l) "
718  << trk_x << " | "
719  << trk_y << " | "
720  << trk_z << " | "
721  << trk_r << " | "
722  << layer << " | "
723  << std::endl;
724 
725  std::cout << " cluster 3D (x|y|z|r|l) "
726  << cluster->getX() << " | "
727  << cluster->getY() << " | "
728  << cluster->getZ() << " | "
729  << get_r(cluster->getX() ,cluster->getY() ) << " | "
730  << layer << " | "
731  << std::endl;
732  */
733  //Get Surface
734  Acts::Vector3D global(trk_x, trk_y, trk_z);
735  // TrkrDefs::subsurfkey subsurfkey = cluster->getSubSurfKey();
736 
737  // std::cout << " subsurfkey: " << subsurfkey << std::endl;
738  std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
739  mapIter = m_surfMaps->tpcSurfaceMap.find(layer);
740 
741  if(mapIter == m_surfMaps->tpcSurfaceMap.end()){
742  std::cout << PHWHERE
743  << "Error: hitsetkey not found in clusterSurfaceMap, layer = " <<
744  trk_r//layer
745  << " hitsetkey = "
746  << hitsetkey << std::endl;
747  continue;// nullptr;
748  }
749  std::cout << " g0: " << global[0] << " g1: " << global[1] << " g2:" << global[2] << std::endl;
750  // double global_phi = atan2(global[1], global[0]);
751  // double global_z = global[2];
752 
753 
754  // Predict which surface index this phi and z will correspond to
755  // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
756  std::vector<Surface> surf_vec = mapIter->second;
757 
758  Acts::Vector3D world(globalpos_d[0], globalpos_d[1],globalpos_d[2]);
759  double world_phi = atan2(world[1], world[0]);
760  double world_z = world[2];
761 
762  double fraction = (world_phi + M_PI) / (2.0 * M_PI);
763  double rounded_nsurf = round( (double) (surf_vec.size()/2) * fraction - 0.5);
764  unsigned int nsurf = (unsigned int) rounded_nsurf;
765  if(world_z < 0)
766  nsurf += surf_vec.size()/2;
767 
768  Surface surface = surf_vec[nsurf];
769 
770  Acts::Vector3D center = surface->center(m_tGeometry->geoContext)
772 
773  // no conversion needed, only used in acts
774  // Acts::Vector3D normal = surface->normal(m_tGeometry->geoContext);
775  double TrkRadius = sqrt(trk_x * trk_x + trk_y * trk_y);
776  double rTrkPhi = TrkRadius * atan2(trk_y, trk_x);//trkphi;
777  double surfRadius = sqrt(center(0)*center(0) + center(1)*center(1));
778  double surfPhiCenter = atan2(center[1], center[0]);
779  double surfRphiCenter = surfPhiCenter * surfRadius;
780  double surfZCenter = center[2];
781 
782  double trklocX = 0;
783  double trklocY = 0;
784  float delta_r = 0;
785 
786  delta_r = TrkRadius - surfRadius;
787  trklocX = rTrkPhi - surfRphiCenter;
788  trklocY = trk_z - surfZCenter;
789  /*
790  std::cout << " clus 2D: (xl|yl) "
791  << cluster->getLocalX() << " | "
792  << cluster->getLocalY()
793  << std::endl;
794 
795  std::cout << " trk 2D: (xl|yl) "
796  << trklocX << " | "
797  << trklocY
798  << std::endl;
799  */
800  float delta_x = trklocX - cluster->getLocalX();
801  float delta_y = trklocY - cluster->getLocalY();
802 
803  float comp_dx = compress_dx(delta_x);
804  float comp_dy = compress_dy(delta_y);
805 
806  // Sabotage the tracking by multiplying the residuals with a random number
807  if (sabotage == -1) {
808  comp_dx = rnd.Uniform(-1, 1) * delta_x;
809  comp_dy = rnd.Uniform(-1, 1) * delta_y;
810  } else if (sabotage == -2) {
811  comp_dx = rnd.Uniform(-10, 10);
812  comp_dy = rnd.Uniform(-10, 10);
813  }
814 
815  /* "dst_data", "dst data","event:seed:"
816  "c3x:c3y:c3z:t3x:t3y:t3z:"
817  "c2x:c2y:c2r:c2l:t2x:t2y:t2r:t2l:"
818  "d2x:d2y:"
819  "cmp_d2x:cmp_d2y:"
820  "pt:eta:phi"
821  */
822  float data[] = {1,
823  (float)iseed,
824  (float)globalpos_d[0],
825  (float)globalpos_d[1],
826  (float)globalpos_d[2],
827  clu_phi,
828  trk_x,
829  trk_y,
830  trk_z,
831  cluster->getLocalX(),
832  cluster->getLocalY(),
833  get_r((float)globalpos_d[0],(float)globalpos_d[1]),
834  (float)layer,
835  (float)trklocX,
836  (float)trklocY,
837  get_r(trk_x,trk_y),
838  (float)layer,
839  delta_x,
840  delta_y,
841  delta_r,
842  comp_dx,
843  comp_dy,
844  track_struct.pt,
845  track_struct.eta,
846  track_struct.phi,
847  (float) track_struct.charge
848  };
849  _dst_data->Fill(data);
850  std::cout << "filled" << "\n";
851 
852  if (apply_compression) {
853  cluster->setLocalX(trklocX - comp_dx);
854  cluster->setLocalY(trklocY - comp_dy);
855  }
856 
857 
858 #ifdef JUNK
859  /*
860  store state angles in (r,phi) and (r,z) plans
861  they are needed to study space charge distortions
862  */
863  /*
864  const auto cosphi( std::cos( cluster.trk_phi ) );
865  const auto sinphi( std::sin( cluster.trk_phi ) );
866  const auto trk_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
867  const auto trk_pr = state->get_px()*cosphi + state->get_py()*sinphi;
868  const auto trk_pz = state->get_pz();
869  cluster.trk_alpha = std::atan2( trk_pphi, trk_pr );
870  cluster.trk_beta = std::atan2( trk_pz, trk_pr );
871  cluster.trk_phi_error = state->get_phi_error();
872  cluster.trk_z_error = state->get_z_error();
873  */
874 #endif
875  }
876  }
877 
878  std::cout << "DSTEmulator::evaluate_tracks - tracks: " << m_container->tracks().size() << std::endl;
879 
880 
881 }
882 
883 //_____________________________________________________________________
884 float DSTEmulator::compress_dx( float in_val )
885 {
886  unsigned short key = m_compressor->compressPhi(in_val);
887  float out_val = m_compressor->decompressPhi(key);
888  return out_val;
889 }
890 
891 //_____________________________________________________________________
892 float DSTEmulator::compress_dy( float in_val )
893 {
894  unsigned short key = m_compressor->compressZ(in_val);
895  float out_val = m_compressor->decompressZ(key);
896  return out_val;
897 }
898 
899 //_____________________________________________________________________
901 {
902 
903  // check maps
904  if( !( m_cluster_hit_map && m_hit_truth_map ) ) return G4HitSet();
905 
906  // check if in map
907  auto map_iter = m_g4hit_map.lower_bound( cluster_key );
908  if( map_iter != m_g4hit_map.end() && cluster_key == map_iter->first )
909  { return map_iter->second; }
910 
911  // find hitset associated to cluster
912  G4HitSet out;
913  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
914  for(const auto& [first,hit_key]:range_adaptor( m_cluster_hit_map->getHits(cluster_key)))
915  {
916 
917  // store hits to g4hit associations
918  TrkrHitTruthAssoc::MMap g4hit_map;
919  m_hit_truth_map->getG4Hits( hitset_key, hit_key, g4hit_map );
920 
921  // find corresponding g4 hist
922  for( auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter )
923  {
924 
925  const auto g4hit_key = truth_iter->second.second;
926  PHG4Hit* g4hit = nullptr;
927 
928  switch( TrkrDefs::getTrkrId( hitset_key ) )
929  {
930  case TrkrDefs::mvtxId:
931  if( m_g4hits_mvtx ) g4hit = m_g4hits_mvtx->findHit( g4hit_key );
932  break;
933 
934  case TrkrDefs::inttId:
935  if( m_g4hits_intt ) g4hit = m_g4hits_intt->findHit( g4hit_key );
936  break;
937 
938  case TrkrDefs::tpcId:
939  if( m_g4hits_tpc ) g4hit = m_g4hits_tpc->findHit( g4hit_key );
940  break;
941 
943  if( m_g4hits_micromegas ) g4hit = m_g4hits_micromegas->findHit( g4hit_key );
944  break;
945 
946  default: break;
947  }
948 
949  if( g4hit ) out.insert( g4hit );
950  else std::cout << "DSTEmulator::find_g4hits - g4hit not found " << g4hit_key << std::endl;
951 
952  }
953  }
954 
955  // insert in map and return
956  return m_g4hit_map.insert( map_iter, std::make_pair( cluster_key, std::move( out ) ) )->second;
957 
958 }
959 
960 //_____________________________________________________________________
961 std::pair<int,int> DSTEmulator::get_max_contributor( SvtxTrack* track ) const
962 {
963  if(!(m_track_map && m_cluster_map && m_g4truthinfo)) return {0,0};
964 
965  // maps MC track id and number of matching g4hits
966  using IdMap = std::map<int,int>;
967  IdMap contributor_map;
968 
969  // loop over clusters
970  for( auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
971  {
972  const auto& cluster_key = *key_iter;
973  for( const auto& hit:find_g4hits( cluster_key ) )
974  {
975  const int trkid = hit->get_trkid();
976  auto iter = contributor_map.lower_bound( trkid );
977  if( iter == contributor_map.end() || iter->first != trkid )
978  {
979  contributor_map.insert(iter, std::make_pair(trkid,1));
980  } else ++iter->second;
981  }
982  }
983 
984  if( contributor_map.empty() ) return {0,0};
985  else return *std::max_element(
986  contributor_map.cbegin(), contributor_map.cend(),
987  []( const IdMap::value_type& first, const IdMap::value_type& second )
988  { return first.second < second.second; } );
989 
990 }
991 
992 //_____________________________________________________________________
993 int DSTEmulator::get_embed( PHG4Particle* particle ) const
994 { return (m_g4truthinfo && particle) ? m_g4truthinfo->isEmbeded( particle->get_primary_id() ):0; }