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