ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SvtxTruthEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SvtxTruthEval.cc
1 #include "SvtxTruthEval.h"
2 
3 #include "BaseTruthEval.h"
4 
5 #include <g4main/PHG4Hit.h>
8 
9 #include <g4main/PHG4Particle.h>
10 
12 #include <trackbase/TrkrDefs.h>
13 
14 #include <tpc/TpcDefs.h>
15 #include <intt/InttDefs.h>
16 #include <mvtx/MvtxDefs.h>
18 
21 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
23 
24 #include <mvtx/CylinderGeom_Mvtx.h>
25 
26 #include <intt/CylinderGeomIntt.h>
27 
28 
29 #include <phool/getClass.h>
30 #include <phool/phool.h> // for PHWHERE
31 
32 #include <TVector3.h>
33 
34 #include <cassert>
35 #include <cmath> // for sqrt, NAN, fabs
36 #include <cstdlib> // for abs
37 #include <iostream>
38 #include <map>
39 #include <set>
40 #include <utility>
41 
42 using namespace std;
43 
45  : _basetrutheval(topNode)
46  , _truthinfo(nullptr)
47  , _strict(false)
48  , _verbosity(0)
49  , _errors(0)
50  , _do_cache(true)
51  , _cache_all_truth_hits()
52  , _cache_all_truth_hits_g4particle()
53  , _cache_all_truth_clusters_g4particle()
54  , _cache_get_innermost_truth_hit()
55  , _cache_get_outermost_truth_hit()
56  , _cache_get_primary_particle_g4hit()
57 {
58  get_node_pointers(topNode);
59  iclus = 0;
60 }
61 
63 {
64  if (_verbosity > 1)
65  {
66  if ((_errors > 0) || (_verbosity > 1))
67  {
68  cout << "SvtxTruthEval::~SvtxTruthEval() - Error Count: " << _errors << endl;
69  }
70  }
71 }
72 
74 {
75  _cache_all_truth_hits.clear();
81 
82  _basetrutheval.next_event(topNode);
83 
84  get_node_pointers(topNode);
85 }
86 
88 std::set<PHG4Hit*> SvtxTruthEval::all_truth_hits()
89 {
90  if (!has_node_pointers())
91  {
92  ++_errors;
93  return std::set<PHG4Hit*>();
94  }
95 
96  if (_do_cache)
97  {
98  if (!_cache_all_truth_hits.empty())
99  {
100  return _cache_all_truth_hits;
101  }
102  }
103 
104  // since the SVTX can be composed of two different trackers this is a
105  // handy function to spill out all the g4hits from both "detectors"
106 
107  std::set<PHG4Hit*> truth_hits;
108 
109  // loop over all the g4hits in the cylinder layers
110  if (_g4hits_svtx)
111  {
113  g4iter != _g4hits_svtx->getHits().second;
114  ++g4iter)
115  {
116  PHG4Hit* g4hit = g4iter->second;
117  truth_hits.insert(g4hit);
118  }
119  }
120 
121  // loop over all the g4hits in the ladder layers
122  if (_g4hits_tracker)
123  {
125  g4iter != _g4hits_tracker->getHits().second;
126  ++g4iter)
127  {
128  PHG4Hit* g4hit = g4iter->second;
129  truth_hits.insert(g4hit);
130  }
131  }
132 
133  // loop over all the g4hits in the maps layers
134  if (_g4hits_maps)
135  {
137  g4iter != _g4hits_maps->getHits().second;
138  ++g4iter)
139  {
140  PHG4Hit* g4hit = g4iter->second;
141  truth_hits.insert(g4hit);
142  }
143  }
144 
145  // loop over all the g4hits in the maicromegas layers
146  if (_g4hits_mms)
147  {
148  for (PHG4HitContainer::ConstIterator g4iter = _g4hits_mms->getHits().first;
149  g4iter != _g4hits_mms->getHits().second;
150  ++g4iter)
151  {
152  PHG4Hit* g4hit = g4iter->second;
153  truth_hits.insert(g4hit);
154  }
155  }
156 
157  if (_do_cache) _cache_all_truth_hits = truth_hits;
158 
159  return truth_hits;
160 }
161 
163 {
164  if (!has_node_pointers())
165  {
166  ++_errors;
167  return std::set<PHG4Hit*>();
168  }
169 
170  if (_strict)
171  {
172  assert(particle);
173  }
174  else if (!particle)
175  {
176  ++_errors;
177  return std::set<PHG4Hit*>();
178  }
179  // if( _cache_all_truth_hits_g4particle.count(particle)==0){
182  }
183 
184  if (_do_cache)
185  {
186  std::map<PHG4Particle*, std::set<PHG4Hit*> >::iterator iter =
187  _cache_all_truth_hits_g4particle.find(particle);
188  if (iter != _cache_all_truth_hits_g4particle.end())
189  {
190  return iter->second;
191  }
192  }
193  //return empty if we dont find anything int he cache
194 
195  std::set<PHG4Hit*> truth_hits;
196  return truth_hits;
197 }
198 
200 
201  std::multimap<const int, PHG4Hit*> temp_clusters_from_particles;
202  // loop over all the g4hits in the cylinder layers
203  if (_g4hits_svtx)
204  {
206  g4iter != _g4hits_svtx->getHits().second;
207  ++g4iter){
208  PHG4Hit* g4hit = g4iter->second;
209  temp_clusters_from_particles.insert(make_pair(g4hit->get_trkid(),g4hit));
210  }
211  }
212 
213  // loop over all the g4hits in the ladder layers
214  if (_g4hits_tracker)
215  {
217  g4iter != _g4hits_tracker->getHits().second;
218  ++g4iter)
219  {
220  PHG4Hit* g4hit = g4iter->second;
221  temp_clusters_from_particles.insert(make_pair(g4hit->get_trkid(),g4hit));
222  }
223  }
224 
225  // loop over all the g4hits in the maps ladder layers
226  if (_g4hits_maps)
227  {
229  g4iter != _g4hits_maps->getHits().second;
230  ++g4iter)
231  {
232  PHG4Hit* g4hit = g4iter->second;
233  temp_clusters_from_particles.insert(make_pair(g4hit->get_trkid(),g4hit));
234  }
235  }
236 
237  // loop over all the g4hits in the micromegas layers
238  if (_g4hits_mms)
239  {
240  for (PHG4HitContainer::ConstIterator g4iter = _g4hits_mms->getHits().first;
241  g4iter != _g4hits_mms->getHits().second;
242  ++g4iter)
243  {
244  PHG4Hit* g4hit = g4iter->second;
245  temp_clusters_from_particles.insert(make_pair(g4hit->get_trkid(),g4hit));
246  }
247  }
248 
249 
250  PHG4TruthInfoContainer::ConstRange range = _truthinfo->GetParticleRange();
251  for(PHG4TruthInfoContainer::ConstIterator iter = range.first;
252  iter != range.second; ++iter){
253  PHG4Particle* g4particle = iter->second;
254  std::set<PHG4Hit*> truth_hits;
255  std::multimap<const int, PHG4Hit*>::const_iterator lower_bound = temp_clusters_from_particles.lower_bound(g4particle->get_track_id());
256  std::multimap<const int, PHG4Hit*>::const_iterator upper_bound = temp_clusters_from_particles.upper_bound(g4particle->get_track_id());
257  std::multimap<const int, PHG4Hit*>::const_iterator cfp_iter;
258  for(cfp_iter = lower_bound;cfp_iter != upper_bound;++cfp_iter){
259  PHG4Hit* g4hit = cfp_iter->second;
260  truth_hits.insert(g4hit);
261  }
262  _cache_all_truth_hits_g4particle.insert(make_pair(g4particle, truth_hits));
263  }
264 
265 }
266 
267 std::map<unsigned int, std::shared_ptr<TrkrCluster> > SvtxTruthEval::all_truth_clusters(PHG4Particle* particle)
268 {
269  if (!has_node_pointers())
270  {
271  ++_errors;
272  return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
273  }
274 
275  if (_strict)
276  {
277  assert(particle);
278  }
279  else if (!particle)
280  {
281  ++_errors;
282  return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
283  }
284 
285  if (_do_cache)
286  {
287  std::map<PHG4Particle*, std::map<unsigned int, std::shared_ptr<TrkrCluster> > >::iterator iter =
289  if (iter != _cache_all_truth_clusters_g4particle.end())
290  {
291  return iter->second;
292  }
293  }
294 
295  if(_verbosity > 1)
296  cout << PHWHERE << " Truth clustering for particle " << particle->get_track_id() << endl;;
297 
298  // get all g4hits for this particle
299  std::set<PHG4Hit*> g4hits = all_truth_hits(particle);
300 
301  float ng4hits = g4hits.size();
302  if(ng4hits == 0)
303  return std::map<unsigned int, std::shared_ptr<TrkrCluster> >();
304 
305  // container for storing truth clusters
306  //std::map<unsigned int, TrkrCluster*> truth_clusters;
307  std::map<unsigned int, std::shared_ptr<TrkrCluster>> truth_clusters;
308 
309  // convert truth hits for this particle to truth clusters in each layer
310  // loop over layers
311 
312  unsigned int layer;
313  for(layer = 0; layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc +_nlayers_mms; ++layer)
314  {
315  float gx = NAN;
316  float gy = NAN;
317  float gz = NAN;
318  float gt = NAN;
319  float gedep = NAN;
320 
321  std::vector<PHG4Hit*> contributing_hits;
322  std::vector<double> contributing_hits_energy;
323  std::vector<std::vector<double>> contributing_hits_entry;
324  std::vector<std::vector<double>> contributing_hits_exit;
325  LayerClusterG4Hits(g4hits, contributing_hits, contributing_hits_energy, contributing_hits_entry, contributing_hits_exit, layer, gx, gy, gz, gt, gedep);
326  if(!(gedep > 0)) continue;
327 
328  // we have the cluster in this layer from this truth particle
329  // add the cluster to a TrkrCluster object
330  TrkrDefs::cluskey ckey;
331  if(layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc) // in TPC
332  {
333  unsigned int side = 0;
334  if(gz > 0) side = 1;
335  // need dummy sector here
336  unsigned int sector = 0;
337  ckey = TpcDefs::genClusKey(layer, sector, side, iclus);
338  }
339  else if(layer < _nlayers_maps) // in MVTX
340  {
341  unsigned int stave = 0;
342  unsigned int chip = 0;
343  ckey = MvtxDefs::genClusKey(layer, stave, chip, iclus);
344  }
345  else if(layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) // in INTT
346  {
347  // dummy ladder and phi ID
348  unsigned int ladderzid = 0;
349  unsigned int ladderphiid = 0;
350  ckey = InttDefs::genClusKey(layer, ladderzid, ladderphiid,iclus);
351  }
352  else if(layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc) // in MICROMEGAS
353  {
354  unsigned int tile = 0;
357  TrkrDefs::hitsetkey hkey = MicromegasDefs::genHitSetKey(layer, segtype, tile);
358  ckey = MicromegasDefs::genClusterKey(hkey, iclus);
359  }
360  else
361  {
362  std::cout << PHWHERE << "Bad layer number: " << layer << std::endl;
363  continue;
364  }
365 
366  std::shared_ptr<TrkrClusterv1> clus(new TrkrClusterv1());
367  clus->setClusKey(ckey);
368  iclus++;
369 
370  // estimate cluster ADC value
371  unsigned int adc_value = getAdcValue(gedep);
372  //std::cout << " gedep " << gedep << " adc_value " << adc_value << std::endl;
373  clus->setAdc(adc_value);
374  clus->setPosition(0, gx);
375  clus->setPosition(1, gy);
376  clus->setPosition(2, gz);
377  clus->setGlobal();
378 
379  // record which g4hits contribute to this truth cluster
380  for(unsigned int i=0; i< contributing_hits.size(); ++i)
381  {
382  _truth_cluster_truth_hit_map.insert(make_pair(ckey,contributing_hits[i]));
383  }
384 
385  // Estimate the size of the truth cluster
386  float g4phisize = NAN;
387  float g4zsize = NAN;
388  G4ClusterSize(layer, contributing_hits_entry, contributing_hits_exit, g4phisize, g4zsize);
389 
390  for(int i1=0;i1<3;++i1)
391  for(int i2=0;i2<3;++i2)
392  {
393  clus->setSize(i1, i2, 0.0);
394  clus->setError(i1, i2, 0.0);
395  }
396  clus->setError(0,0,gedep); // stores truth energy
397  clus->setSize(1, 1, g4phisize);
398  clus->setSize(2, 2, g4zsize);
399  clus->setError(1, 1, g4phisize/sqrt(12));
400  clus->setError(2, 2, g4zsize/sqrt(12.0));
401 
402  truth_clusters.insert(make_pair(layer, clus));
403 
404  } // end loop over layers for this particle
405 
406  if (_do_cache) _cache_all_truth_clusters_g4particle.insert(make_pair(particle, truth_clusters));
407 
408  return truth_clusters;
409 }
410 
411 void SvtxTruthEval::LayerClusterG4Hits(std::set<PHG4Hit*> truth_hits, std::vector<PHG4Hit*> &contributing_hits, std::vector<double> &contributing_hits_energy, std::vector<std::vector<double>> &contributing_hits_entry, std::vector<std::vector<double>> &contributing_hits_exit, float layer, float &x, float &y, float &z, float &t, float &e)
412 {
413  bool use_geo = true;
414 
415  // Given a set of g4hits, cluster them within a given layer of the TPC
416 
417  float gx = 0.0;
418  float gy = 0.0;
419  float gz = 0.0;
420  float gr = 0.0;
421  float gt = 0.0;
422  float gwt = 0.0;
423 
424  if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc) // in TPC
425  {
426  //cout << "layer = " << layer << " _nlayers_maps " << _nlayers_maps << " _nlayers_intt " << _nlayers_intt << endl;
427 
428  // This calculates the truth cluster position for the TPC from all of the contributing g4hits from a g4particle, typically 2-4 for the TPC
429  // Complicated, since only the part of the energy that is collected within a layer contributes to the position
430  //===============================================================================
431 
433  // get layer boundaries here for later use
434  // radii of layer boundaries
435  float rbin = GeoLayer->get_radius() - GeoLayer->get_thickness() / 2.0;
436  float rbout = GeoLayer->get_radius() + GeoLayer->get_thickness() / 2.0;
437 
438  if(_verbosity > 1)
439  cout << " TruthEval::LayerCluster hits for layer " << layer << " with rbin " << rbin << " rbout " << rbout << endl;
440 
441  // we do not assume that the truth hits know what layer they are in
442  for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
443  iter != truth_hits.end();
444  ++iter)
445  {
446 
447  PHG4Hit* this_g4hit = *iter;
448  float rbegin = sqrt(this_g4hit->get_x(0) * this_g4hit->get_x(0) + this_g4hit->get_y(0) * this_g4hit->get_y(0));
449  float rend = sqrt(this_g4hit->get_x(1) * this_g4hit->get_x(1) + this_g4hit->get_y(1) * this_g4hit->get_y(1));
450  //cout << " Eval: g4hit " << this_g4hit->get_hit_id() << " layer " << layer << " rbegin " << rbegin << " rend " << rend << endl;
451 
452  // make sure the entry point is at lower radius
453  float xl[2];
454  float yl[2];
455  float zl[2];
456 
457  if (rbegin < rend)
458  {
459  xl[0] = this_g4hit->get_x(0);
460  yl[0] = this_g4hit->get_y(0);
461  zl[0] = this_g4hit->get_z(0);
462  xl[1] = this_g4hit->get_x(1);
463  yl[1] = this_g4hit->get_y(1);
464  zl[1] = this_g4hit->get_z(1);
465  }
466  else
467  {
468  xl[0] = this_g4hit->get_x(1);
469  yl[0] = this_g4hit->get_y(1);
470  zl[0] = this_g4hit->get_z(1);
471  xl[1] = this_g4hit->get_x(0);
472  yl[1] = this_g4hit->get_y(0);
473  zl[1] = this_g4hit->get_z(0);
474  swap(rbegin, rend);
475  //cout << "swapped in and out " << endl;
476  }
477 
478  // check that the g4hit is not completely outside the cluster layer. Just skip this g4hit if it is
479  if (rbegin < rbin && rend < rbin)
480  continue;
481  if (rbegin > rbout && rend > rbout)
482  continue;
483 
484 
485  if(_verbosity > 1)
486  {
487  cout << " keep g4hit with rbegin " << rbegin << " rend " << rend
488  << " xbegin " << xl[0] << " xend " << xl[1]
489  << " ybegin " << yl[0] << " yend " << yl[1]
490  << " zbegin " << zl[0] << " zend " << zl[1]
491  << endl;
492  }
493 
494 
495  float xin = xl[0];
496  float yin = yl[0];
497  float zin = zl[0];
498  float xout = xl[1];
499  float yout = yl[1];
500  float zout = zl[1];
501 
502  float t = NAN;
503 
504  if (rbegin < rbin)
505  {
506  // line segment begins before boundary, find where it crosses
507  t = line_circle_intersection(xl, yl, zl, rbin);
508  if (t > 0)
509  {
510  xin = xl[0] + t * (xl[1] - xl[0]);
511  yin = yl[0] + t * (yl[1] - yl[0]);
512  zin = zl[0] + t * (zl[1] - zl[0]);
513  }
514  }
515 
516  if (rend > rbout)
517  {
518  // line segment ends after boundary, find where it crosses
519  t = line_circle_intersection(xl, yl, zl, rbout);
520  if (t > 0)
521  {
522  xout = xl[0] + t * (xl[1] - xl[0]);
523  yout = yl[0] + t * (yl[1] - yl[0]);
524  zout = zl[0] + t * (zl[1] - zl[0]);
525  }
526  }
527 
528  double rin = sqrt(xin*xin + yin*yin);
529  double rout = sqrt(xout*xout + yout*yout);
530 
531  // we want only the fraction of edep inside the layer
532  double efrac = this_g4hit->get_edep() * (rout - rin) / (rend - rbegin);
533  gx += (xin + xout) * 0.5 * efrac;
534  gy += (yin + yout) * 0.5 * efrac;
535  gz += (zin + zout) * 0.5 * efrac;
536  gt += this_g4hit->get_avg_t() * efrac;
537  gr += (rin + rout) * 0.5 * efrac;
538  gwt += efrac;
539 
540  if(_verbosity > 1)
541  {
542  cout << " rin " << rin << " rout " << rout
543  << " xin " << xin << " xout " << xout << " yin " << yin << " yout " << yout << " zin " << zin << " zout " << zout
544  << " edep " << this_g4hit->get_edep()
545  << " this_edep " << efrac << endl;
546  cout << " xavge " << (xin+xout) * 0.5 << " yavge " << (yin+yout) * 0.5 << " zavge " << (zin+zout) * 0.5 << " ravge " << (rin+rout) * 0.5
547  << endl;
548  }
549 
550  // Capture entry and exit points
551  std::vector<double> entry_loc;
552  entry_loc.push_back(xin);
553  entry_loc.push_back(yin);
554  entry_loc.push_back(zin);
555  std::vector<double> exit_loc;
556  exit_loc.push_back(xout);
557  exit_loc.push_back(yout);
558  exit_loc.push_back(zout);
559 
560  // this_g4hit is inside the layer, add it to the vectors
561  contributing_hits.push_back(this_g4hit);
562  contributing_hits_energy.push_back( this_g4hit->get_edep() * (zout - zin) / (zl[1] - zl[0]) );
563  contributing_hits_entry.push_back(entry_loc);
564  contributing_hits_exit.push_back(exit_loc);
565 
566  } // loop over contributing hits
567 
568  if(gwt == 0)
569  {
570  e = gwt;
571  return; // will be discarded
572  }
573 
574  gx /= gwt;
575  gy /= gwt;
576  gz /= gwt;
577  gr = (rbin + rbout) * 0.5;
578  gt /= gwt;
579 
580  if(_verbosity > 1)
581  {
582  cout << " weighted means: gx " << gx << " gy " << gy << " gz " << gz << " gr " << gr << " e " << gwt << endl;
583  }
584 
585  if(use_geo)
586  {
587  // The energy weighted values above have significant scatter due to fluctuations in the energy deposit from Geant
588  // Calculate the geometric mean positions instead
589  float rentry = 999.0;
590  float xentry = 999.0;
591  float yentry = 999.0;
592  float zentry = 999.0;
593  float rexit = - 999.0;
594  float xexit = -999.0;
595  float yexit = -999.0;
596  float zexit = -999.0;
597 
598  for(unsigned int ientry = 0; ientry < contributing_hits_entry.size(); ++ientry)
599  {
600  float tmpx = contributing_hits_entry[ientry][0];
601  float tmpy = contributing_hits_entry[ientry][1];
602  float tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
603 
604  if(tmpr < rentry)
605  {
606  rentry = tmpr;
607  xentry = contributing_hits_entry[ientry][0];
608  yentry = contributing_hits_entry[ientry][1];
609  zentry = contributing_hits_entry[ientry][2];
610  }
611 
612  tmpx = contributing_hits_exit[ientry][0];
613  tmpy = contributing_hits_exit[ientry][1];
614  tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);
615 
616  if(tmpr > rexit)
617  {
618  rexit = tmpr;
619  xexit = contributing_hits_exit[ientry][0];
620  yexit = contributing_hits_exit[ientry][1];
621  zexit = contributing_hits_exit[ientry][2];
622  }
623  }
624 
625  float geo_r = (rentry+rexit)*0.5;
626  float geo_x = (xentry+xexit)*0.5;
627  float geo_y = (yentry+yexit)*0.5;
628  float geo_z = (zentry+zexit)*0.5;
629 
630  if(rexit > 0)
631  {
632  gx = geo_x;
633  gy = geo_y;
634  gz = geo_z;
635  gr = geo_r;
636  }
637 
638 
639  if(_verbosity > 1)
640  {
641  cout << " rentry " << rentry << " rexit " << rexit
642  << " xentry " << xentry << " xexit " << xexit << " yentry " << yentry << " yexit " << yexit << " zentry " << zentry << " zexit " << zexit << endl;
643 
644  cout << " geometric means: geo_x " << geo_x << " geo_y " << geo_y << " geo_z " << geo_z << " geo r " << geo_r << " e " << gwt << endl << endl;
645  }
646  }
647 
648  } // if TPC
649  else
650  {
651  // not TPC, one g4hit per cluster
652  for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
653  iter != truth_hits.end();
654  ++iter)
655  {
656 
657  PHG4Hit* this_g4hit = *iter;
658 
659  if(this_g4hit->get_layer() != (unsigned int) layer) continue;
660 
661  gx = this_g4hit->get_avg_x();
662  gy = this_g4hit->get_avg_y();
663  gz = this_g4hit->get_avg_z();
664  gt = this_g4hit->get_avg_t();
665  gwt += this_g4hit->get_edep();
666 
667  // Capture entry and exit points
668  std::vector<double> entry_loc;
669  entry_loc.push_back(this_g4hit->get_x(0));
670  entry_loc.push_back(this_g4hit->get_y(0));
671  entry_loc.push_back(this_g4hit->get_z(0));
672  std::vector<double> exit_loc;
673  exit_loc.push_back(this_g4hit->get_x(1));
674  exit_loc.push_back(this_g4hit->get_y(1));
675  exit_loc.push_back(this_g4hit->get_z(1));
676 
677  // this_g4hit is inside the layer, add it to the vectors
678  contributing_hits.push_back(this_g4hit);
679  contributing_hits_energy.push_back( this_g4hit->get_edep() );
680  contributing_hits_entry.push_back(entry_loc);
681  contributing_hits_exit.push_back(exit_loc);
682  }
683  } // not TPC
684 
685  x = gx;
686  y = gy;
687  z = gz;
688  t = gt;
689  e = gwt;
690 
691  return;
692 }
693 
694 void SvtxTruthEval::G4ClusterSize(unsigned int layer, std::vector<std::vector<double>> contributing_hits_entry,std::vector<std::vector<double>> contributing_hits_exit, float &g4phisize, float &g4zsize)
695 {
696 
697  // sort the contributing g4hits in radius
698  double inner_radius = 100.;
699  double inner_x = NAN;
700  double inner_y = NAN;
701  double inner_z = NAN;;
702 
703  double outer_radius = 0.;
704  double outer_x = NAN;
705  double outer_y = NAN;
706  double outer_z = NAN;
707 
708  for(unsigned int ihit=0;ihit<contributing_hits_entry.size(); ++ihit)
709  {
710  double rad1 = sqrt(pow(contributing_hits_entry[ihit][0], 2) + pow(contributing_hits_entry[ihit][1], 2));
711  if(rad1 < inner_radius)
712  {
713  inner_radius = rad1;
714  inner_x = contributing_hits_entry[ihit][0];
715  inner_y = contributing_hits_entry[ihit][1];
716  inner_z = contributing_hits_entry[ihit][2];
717  }
718 
719  double rad2 = sqrt(pow(contributing_hits_exit[ihit][0], 2) + pow(contributing_hits_exit[ihit][1], 2));
720  if(rad2 > outer_radius)
721  {
722  outer_radius = rad2;
723  outer_x = contributing_hits_exit[ihit][0];
724  outer_y = contributing_hits_exit[ihit][1];
725  outer_z = contributing_hits_exit[ihit][2];
726  }
727  }
728 
729  double inner_phi = atan2(inner_y, inner_x);
730  double outer_phi = atan2(outer_y, outer_x);
731  double avge_z = (outer_z + inner_z) / 2.0;
732 
733  // Now fold these with the expected diffusion and shaping widths
734  // assume spread is +/- equals this many sigmas times diffusion and shaping when extending the size
735  double sigmas = 2.0;
736 
737  double radius = (inner_radius + outer_radius)/2.;
738  if(radius > 28 && radius < 80) // TPC
739  {
741 
742  double tpc_length = 211.0; // cm
743  double drift_velocity = 8.0 / 1000.0; // cm/ns
744 
745  // Phi size
746  //======
747  double diffusion_trans = 0.006; // cm/SQRT(cm)
748  double phidiffusion = diffusion_trans * sqrt(tpc_length / 2. - fabs(avge_z));
749 
750  double added_smear_trans = 0.085; // cm
751  double gem_spread = 0.04; // 400 microns
752 
753  if(outer_phi < inner_phi) swap(outer_phi, inner_phi);
754 
755  // convert diffusion from cm to radians
756  double g4max_phi = outer_phi + sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) / radius;
757  double g4min_phi = inner_phi - sigmas * sqrt( pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2) ) / radius;
758 
759  // find the bins containing these max and min z edges
760  unsigned int phibinmin = layergeom->get_phibin(g4min_phi);
761  unsigned int phibinmax = layergeom->get_phibin(g4max_phi);
762  unsigned int phibinwidth = phibinmax - phibinmin + 1;
763  g4phisize = (double) phibinwidth * layergeom->get_phistep() * layergeom->get_radius();
764 
765  // Z size
766  //=====
767  double g4max_z = 0;
768  double g4min_z = 0;
769 
770  outer_z = fabs(outer_z);
771  inner_z = fabs(inner_z);
772 
773  double diffusion_long = 0.015; // cm/SQRT(cm)
774  double zdiffusion = diffusion_long * sqrt(tpc_length / 2. - fabs(avge_z)) ;
775  double zshaping_lead = 32.0 * drift_velocity; // ns * cm/ns = cm
776  double zshaping_tail = 48.0 * drift_velocity;
777  double added_smear_long = 0.105; // cm
778 
779  // largest z reaches gems first, make that the outer z
780  if(outer_z < inner_z) swap(outer_z, inner_z);
781  g4max_z = outer_z + sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_lead, 2));
782  g4min_z = inner_z - sigmas*sqrt(pow(zdiffusion,2) + pow(added_smear_long,2) + pow(zshaping_tail, 2));
783 
784  // find the bins containing these max and min z edges
785  unsigned int binmin = layergeom->get_zbin(g4min_z);
786  unsigned int binmax = layergeom->get_zbin(g4max_z);
787  if(binmax < binmin) swap(binmax, binmin);
788  unsigned int binwidth = binmax - binmin + 1;
789 
790  // multiply total number of bins that include the edges by the bin size
791  g4zsize = (double) binwidth * layergeom->get_zstep();
792  }
793  else if(radius > 5 && radius < 20) // INTT
794  {
795  // All we have is the position and layer number
796 
797  CylinderGeomIntt *layergeom = dynamic_cast<CylinderGeomIntt *>(_intt_geom_container->GetLayerGeom(layer));
798 
799  // inner location
800  double world_inner[3] = {inner_x, inner_y, inner_z};
801  TVector3 world_inner_vec = {inner_x, inner_y, inner_z};
802 
803  int segment_z_bin, segment_phi_bin;
804  layergeom->find_indices_from_world_location(segment_z_bin, segment_phi_bin, world_inner);
805 
806  TVector3 local_inner_vec = layergeom->get_local_from_world_coords(segment_z_bin, segment_phi_bin, world_inner_vec);
807  double yin = local_inner_vec[1];
808  double zin = local_inner_vec[2];
809  int strip_y_index, strip_z_index;
810  layergeom->find_strip_index_values(segment_z_bin, yin, zin, strip_y_index, strip_z_index);
811 
812  // outer location
813  double world_outer[3] = {outer_x, outer_y, outer_z};
814  TVector3 world_outer_vec = {outer_x, outer_y, outer_z};
815 
816  layergeom->find_indices_from_world_location(segment_z_bin, segment_phi_bin, world_outer);
817 
818  TVector3 local_outer_vec = layergeom->get_local_from_world_coords(segment_z_bin, segment_phi_bin, world_outer_vec);
819  double yout = local_outer_vec[1];
820  double zout = local_outer_vec[2];
821  int strip_y_index_out, strip_z_index_out;
822  layergeom->find_strip_index_values(segment_z_bin, yout, zout, strip_y_index_out, strip_z_index_out);
823 
824  int strips = abs(strip_y_index_out - strip_y_index) + 1;
825  int cols = abs(strip_z_index_out - strip_z_index) + 1;
826 
827 
828  double strip_width = (double) strips * layergeom->get_strip_y_spacing(); // cm
829  double strip_length = (double) cols * layergeom->get_strip_z_spacing(); // cm
830 
831  g4phisize = strip_width;
832  g4zsize = strip_length;
833 
834  /*
835  if(Verbosity() > 1)
836  cout << " INTT: layer " << layer << " strips " << strips << " strip pitch " << layergeom->get_strip_y_spacing() << " g4phisize "<< g4phisize
837  << " columns " << cols << " strip_z_spacing " << layergeom->get_strip_z_spacing() << " g4zsize " << g4zsize << endl;
838  */
839  }
840  else if(radius > 80) // MICROMEGAS
841  {
842  // made up for now
843  g4phisize = 300e-04;
844  g4zsize = 300e-04;
845  }
846  else // MVTX
847  {
848  unsigned int stave, stave_outer;
849  unsigned int chip, chip_outer;
850  int row, row_outer;
851  int column, column_outer;
852 
853  // add diffusion to entry and exit locations
854  double max_diffusion_radius = 25.0e-4; // 25 microns
855  double min_diffusion_radius = 8.0e-4; // 8 microns
856 
857  CylinderGeom_Mvtx *layergeom = dynamic_cast<CylinderGeom_Mvtx *>(_mvtx_geom_container->GetLayerGeom(layer));
858 
859  TVector3 world_inner = {inner_x, inner_y, inner_z};
860  std::vector<double> world_inner_vec = { world_inner[0], world_inner[1], world_inner[2] };
861  layergeom->get_sensor_indices_from_world_coords(world_inner_vec, stave, chip);
862  TVector3 local_inner = layergeom->get_local_from_world_coords(stave, chip, world_inner);
863 
864  TVector3 world_outer = {outer_x, outer_y, outer_z};
865  std::vector<double> world_outer_vec = { world_outer[0], world_outer[1], world_outer[2] };
866  layergeom->get_sensor_indices_from_world_coords(world_outer_vec, stave_outer, chip_outer);
867  TVector3 local_outer = layergeom->get_local_from_world_coords(stave_outer, chip_outer, world_outer);
868 
869  double diff = max_diffusion_radius * 0.6; // factor of 0.6 gives decent agreement with low occupancy reco clusters
870  if(local_outer[0] < local_inner[0])
871  diff = -diff;
872  local_outer[0] += diff;
873  local_inner[0] -= diff;
874 
875  double diff_outer = min_diffusion_radius * 0.6;
876  if(local_outer[2] < local_inner[2])
877  diff_outer = -diff_outer;
878  local_outer[2] += diff_outer;
879  local_inner[2] -= diff_outer;
880 
881  layergeom->get_pixel_from_local_coords(local_inner, row, column);
882  layergeom->get_pixel_from_local_coords(local_outer, row_outer, column_outer);
883 
884  if(row_outer < row) swap(row_outer, row);
885  unsigned int rows = row_outer - row + 1;
886  g4phisize = (double) rows * layergeom->get_pixel_x();
887 
888  if(column_outer < column) swap(column_outer, column);
889  unsigned int columns = column_outer - column + 1;
890  g4zsize = (double) columns * layergeom->get_pixel_z();
891 
892  /*
893  if(Verbosity() > 1)
894  cout << " MVTX: layer " << layer << " rows " << rows << " pixel x " << layergeom->get_pixel_x() << " g4phisize "<< g4phisize
895  << " columns " << columns << " pixel_z " << layergeom->get_pixel_z() << " g4zsize " << g4zsize << endl;
896  */
897 
898  }
899 
900 }
901 
903 {
904  std::set<PHG4Hit *> g4hit_set;
905 
906 
907  std::pair<std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator,
908  std::multimap<TrkrDefs::cluskey,PHG4Hit*>::iterator> ret = _truth_cluster_truth_hit_map.equal_range(ckey);
909 
910  for(std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator jter = ret.first; jter != ret.second; ++jter)
911  {
912  g4hit_set.insert(jter->second);
913  }
914 
915  return g4hit_set;
916 }
917 
919 {
920  if (!has_node_pointers())
921  {
922  ++_errors;
923  return nullptr;
924  }
925 
926  if (_strict)
927  {
928  assert(particle);
929  }
930  else if (!particle)
931  {
932  ++_errors;
933  return nullptr;
934  }
935 
936  PHG4Hit* innermost_hit = nullptr;
937  float innermost_radius = FLT_MAX;
938 
939  std::set<PHG4Hit*> truth_hits = all_truth_hits(particle);
940  for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
941  iter != truth_hits.end();
942  ++iter)
943  {
944  PHG4Hit* candidate = *iter;
945  float x = candidate->get_x(0); // use entry points
946  float y = candidate->get_y(0); // use entry points
947  float r = sqrt(x * x + y * y);
948  if (r < innermost_radius)
949  {
950  innermost_radius = r;
951  innermost_hit = candidate;
952  }
953  }
954 
955  return innermost_hit;
956 }
957 
959 {
960  if (!has_node_pointers())
961  {
962  ++_errors;
963  return nullptr;
964  }
965 
966  if (_strict)
967  {
968  assert(particle);
969  }
970  else if (!particle)
971  {
972  ++_errors;
973  return nullptr;
974  }
975 
976  PHG4Hit* outermost_hit = nullptr;
977  float outermost_radius = FLT_MAX * -1.0;
978 
979  if (_do_cache)
980  {
981  std::map<PHG4Particle*, PHG4Hit*>::iterator iter =
982  _cache_get_outermost_truth_hit.find(particle);
983  if (iter != _cache_get_outermost_truth_hit.end())
984  {
985  return iter->second;
986  }
987  }
988 
989  std::set<PHG4Hit*> truth_hits = all_truth_hits(particle);
990  for (std::set<PHG4Hit*>::iterator iter = truth_hits.begin();
991  iter != truth_hits.end();
992  ++iter)
993  {
994  PHG4Hit* candidate = *iter;
995  float x = candidate->get_x(1); // use exit points
996  float y = candidate->get_y(1); // use exit points
997  float r = sqrt(x * x + y * y);
998  if (r > outermost_radius)
999  {
1000  outermost_radius = r;
1001  outermost_hit = candidate;
1002  }
1003  }
1004  if (_do_cache) _cache_get_outermost_truth_hit.insert(make_pair(particle, outermost_hit));
1005 
1006  return outermost_hit;
1007 }
1008 
1010 {
1011  return _basetrutheval.get_particle(g4hit);
1012 }
1013 
1015 {
1016  return _basetrutheval.get_embed(particle);
1017 }
1018 
1020 {
1021  return _basetrutheval.get_vertex(particle);
1022 }
1023 
1025 {
1026  return _basetrutheval.is_primary(particle);
1027 }
1028 
1030 {
1031  if (!has_node_pointers())
1032  {
1033  ++_errors;
1034  return nullptr;
1035  }
1036 
1037  if (_strict)
1038  {
1039  assert(g4hit);
1040  }
1041  else if (!g4hit)
1042  {
1043  ++_errors;
1044  return nullptr;
1045  }
1046 
1047  if (_do_cache)
1048  {
1049  std::map<PHG4Hit*, PHG4Particle*>::iterator iter =
1051  if (iter != _cache_get_primary_particle_g4hit.end())
1052  {
1053  return iter->second;
1054  }
1055  }
1056 
1058 
1059  if (_do_cache) _cache_get_primary_particle_g4hit.insert(make_pair(g4hit, primary));
1060 
1061  if (_strict)
1062  {
1063  assert(primary);
1064  }
1065  else if (!primary)
1066  {
1067  ++_errors;
1068  }
1069 
1070  return primary;
1071 }
1072 
1074 {
1075  return _basetrutheval.get_primary_particle(particle);
1076 }
1077 
1079 {
1080  return _basetrutheval.get_particle(trackid);
1081 }
1082 
1084 {
1085  return _basetrutheval.is_g4hit_from_particle(g4hit, particle);
1086 }
1087 
1089 {
1090  return _basetrutheval.are_same_particle(p1, p2);
1091 }
1092 
1094 {
1095  return _basetrutheval.are_same_vertex(vtx1, vtx2);
1096 }
1097 
1099 {
1100  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1101 
1102  _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
1103  _g4hits_svtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
1104  _g4hits_tracker = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
1105  _g4hits_maps = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
1106 
1107  _mms_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
1108  _tpc_geom_container = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1109  _intt_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1110  _mvtx_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1111 
1112  return;
1113 }
1114 
1116 {
1117  if (_strict)
1118  assert(_truthinfo);
1119  else if (!_truthinfo)
1120  return false;
1121 
1122  if (_strict)
1124  else if (!_g4hits_mms && !_g4hits_svtx && !_g4hits_tracker && !_g4hits_maps)
1125  return false;
1126 
1127  return true;
1128 }
1129 
1130 float SvtxTruthEval::line_circle_intersection(float x[], float y[], float z[], float radius)
1131 {
1132  // parameterize the line in terms of t (distance along the line segment, from 0-1) as
1133  // x = x0 + t * (x1-x0); y=y0 + t * (y1-y0); z = z0 + t * (z1-z0)
1134  // parameterize the cylinder (centered at x,y = 0,0) as x^2 + y^2 = radius^2, then
1135  // (x0 + t*(x1-z0))^2 + (y0+t*(y1-y0))^2 = radius^2
1136  // (x0^2 + y0^2 - radius^2) + (2x0*(x1-x0) + 2y0*(y1-y0))*t + ((x1-x0)^2 + (y1-y0)^2)*t^2 = 0 = C + B*t + A*t^2
1137  // quadratic with: A = (x1-x0)^2+(y1-y0)^2 ; B = 2x0*(x1-x0) + 2y0*(y1-y0); C = x0^2 + y0^2 - radius^2
1138  // solution: t = (-B +/- sqrt(B^2 - 4*A*C)) / (2*A)
1139 
1140  float A = (x[1] - x[0]) * (x[1] - x[0]) + (y[1] - y[0]) * (y[1] - y[0]);
1141  float B = 2.0 * x[0] * (x[1] - x[0]) + 2.0 * y[0] * (y[1] - y[0]);
1142  float C = x[0] * x[0] + y[0] * y[0] - radius * radius;
1143  float tup = (-B + sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1144  float tdn = (-B - sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1145 
1146  // The limits are 0 and 1, but we allow a little for floating point precision
1147  float t;
1148  if (tdn >= -0.0e-4 && tdn <= 1.0004)
1149  t = tdn;
1150  else if (tup >= -0.0e-4 && tup <= 1.0004)
1151  t = tup;
1152  else
1153  {
1154  cout << PHWHERE << " **** Oops! No valid solution for tup or tdn, tdn = " << tdn << " tup = " << tup << endl;
1155  cout << " radius " << radius << " rbegin " << sqrt(x[0] * x[0] + y[0] * y[0]) << " rend " << sqrt(x[1] * x[1] + y[1] * y[1]) << endl;
1156  cout << " x0 " << x[0] << " x1 " << x[1] << endl;
1157  cout << " y0 " << y[0] << " y1 " << y[1] << endl;
1158  cout << " z0 " << z[0] << " z1 " << z[1] << endl;
1159  cout << " A " << A << " B " << B << " C " << C << endl;
1160 
1161  t = -1;
1162  }
1163 
1164  return t;
1165 }
1166 
1167 unsigned int SvtxTruthEval::getAdcValue(double gedep)
1168 {
1169  // see TPC digitizer for algorithm
1170 
1171  // drift electrons per GeV of energy deposited in the TPC
1172  double Ne_dEdx = 1.56; // keV/cm
1173  double CF4_dEdx = 7.00; // keV/cm
1174  double Ne_NTotal = 43; // Number/cm
1175  double CF4_NTotal = 100; // Number/cm
1176  double Tpc_NTot = 0.5*Ne_NTotal + 0.5*CF4_NTotal;
1177  double Tpc_dEdx = 0.5*Ne_dEdx + 0.5*CF4_dEdx;
1178  double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
1179  double electrons_per_gev = Tpc_ElectronsPerKeV * 1e6;
1180 
1181  double gem_amplification = 1400; // GEM output electrons per drifted electron
1182  double input_electrons = gedep * electrons_per_gev * gem_amplification;
1183 
1184  // convert electrons after GEM to ADC output
1185  double ChargeToPeakVolts = 20;
1186  double ADCSignalConversionGain = ChargeToPeakVolts * 1.60e-04 * 2.4; // 20 (or 30) mV/fC * fC/electron * scaleup factor
1187  double adc_input_voltage = input_electrons * ADCSignalConversionGain; // mV, see comments above
1188  unsigned int adc_output = (unsigned int) (adc_input_voltage * 1024.0 / 2200.0); // input voltage x 1024 channels over 2200 mV max range
1189  if (adc_output > 1023) adc_output = 1023;
1190 
1191  return adc_output;
1192 }