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