ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InttHitReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InttHitReco.cc
1 #include "PHG4InttHitReco.h"
2 
4 
5 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
7 
8 #include <trackbase/TrkrDefs.h>
9 #include <trackbase/TrkrHit.h> // for TrkrHit
10 #include <trackbase/TrkrHitv2.h> // for TrkrHit
11 #include <trackbase/TrkrHitSet.h>
16 
17 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
18 
19 #include <intt/InttDefs.h>
20 
21 #include <g4main/PHG4Hit.h>
23 #include <g4main/PHG4Utils.h>
24 
26 #include <fun4all/SubsysReco.h> // for SubsysReco
27 
28 #include <phool/PHCompositeNode.h>
29 #include <phool/PHIODataNode.h>
30 #include <phool/PHNode.h> // for PHNode
31 #include <phool/PHNodeIterator.h>
32 #include <phool/PHObject.h> // for PHObject
33 #include <phool/getClass.h>
34 #include <phool/phool.h> // for PHWHERE
35 
36 #include <TSystem.h>
37 
38 #include <cmath>
39 #include <cstdlib>
40 #include <iostream>
41 #include <map> // for _Rb_tree_const_it...
42 #include <memory> // for allocator_traits<...
43 #include <utility> // for pair, swap, make_...
44 #include <vector> // for vector
45 
46 using namespace std;
47 
49  : SubsysReco(name)
50  , PHParameterInterface(name)
51  , m_Detector("INTT")
52  , m_Tmin(NAN)
53  , m_Tmax(NAN)
54 {
56 
57  m_HitNodeName = "G4HIT_" + m_Detector;
58  m_CellNodeName = "G4CELL_" + m_Detector;
59  m_GeoNodeName = "CYLINDERGEOM_" + m_Detector;
60  m_LocalOutVec = gsl_vector_alloc(3);
61  m_PathVec = gsl_vector_alloc(3);
62  m_SegmentVec = gsl_vector_alloc(3);
63 }
64 
66 {
67  gsl_vector_free(m_LocalOutVec);
68  gsl_vector_free(m_PathVec);
69  gsl_vector_free(m_SegmentVec);
70 }
71 
73 {
74  PHNodeIterator iter(topNode);
75 
76  // Looking for the DST node
77  PHCompositeNode *dstNode;
78  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
79  if (!dstNode)
80  {
81  std::cout << PHWHERE << "DST Node missing, exiting." << std::endl;
82  gSystem->Exit(1);
83  exit(1);
84  }
85 
86  PHCompositeNode *runNode;
87  runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
88  if (!runNode)
89  {
90  cout << Name() << "RUN Node missing, exiting." << endl;
91  gSystem->Exit(1);
92  exit(1);
93  }
94  PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
95  if (!parNode)
96  {
97  cout << Name() << "PAR Node missing, exiting." << endl;
98  gSystem->Exit(1);
99  exit(1);
100  }
101  string paramnodename = "G4CELLPARAM_" + m_Detector;
102 
103  PHNodeIterator runiter(runNode);
104  PHCompositeNode *RunDetNode =
105  dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode",
106  m_Detector));
107  if (!RunDetNode)
108  {
109  RunDetNode = new PHCompositeNode(m_Detector);
110  runNode->addNode(RunDetNode);
111  }
112 
113  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
114  if (!g4hit)
115  {
116  std::cout << "Could not locate g4 hit node " << m_HitNodeName << std::endl;
117  exit(1);
118  }
119 
120  auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
121  if (!hitsetcontainer)
122  {
123  PHNodeIterator dstiter(dstNode);
124  PHCompositeNode *DetNode =
125  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
126  if (!DetNode)
127  {
128  DetNode = new PHCompositeNode("TRKR");
129  dstNode->addNode(DetNode);
130  }
131 
132  hitsetcontainer = new TrkrHitSetContainerv1;
133  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
134  DetNode->addNode(newNode);
135  }
136 
137  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
138  if (!hittruthassoc)
139  {
140  PHNodeIterator dstiter(dstNode);
141  PHCompositeNode *DetNode =
142  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
143  if (!DetNode)
144  {
145  DetNode = new PHCompositeNode("TRKR");
146  dstNode->addNode(DetNode);
147  }
148 
149  hittruthassoc = new TrkrHitTruthAssocv1;
150  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
151  DetNode->addNode(newNode);
152  }
153 
154  PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, m_GeoNodeName);
155  if (!geo)
156  {
157  std::cout << "Could not locate geometry node " << m_GeoNodeName << std::endl;
158  exit(1);
159  }
160 
161  if (Verbosity() > 0)
162  {
163  geo->identify();
164  }
165 
167  SaveToNodeTree(RunDetNode, paramnodename);
168  // save this to the parNode for use
169  PHNodeIterator parIter(parNode);
170  PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", m_Detector));
171  if (!ParDetNode)
172  {
173  ParDetNode = new PHCompositeNode(m_Detector);
174  parNode->addNode(ParDetNode);
175  }
176  PutOnParNode(ParDetNode, m_GeoNodeName);
177  m_Tmin = get_double_param("tmin");
178  m_Tmax = get_double_param("tmax");
179 
181 }
182 
184 {
185  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
186  if (!g4hit)
187  {
188  std::cout << "Could not locate g4 hit node " << m_HitNodeName << std::endl;
189  exit(1);
190  }
191 
192  // Get the TrkrHitSetContainer node
193  auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
194  if (!hitsetcontainer)
195  {
196  cout << "Could not locate TRKR_HITSET node, quit! " << endl;
197  exit(1);
198  }
199 
200  // Get the TrkrHitTruthAssoc node
201  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
202  if (!hittruthassoc)
203  {
204  cout << "Could not locate TRKR_HITTRUTHASSOC node, quit! " << endl;
205  exit(1);
206  }
207 
208  PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, m_GeoNodeName);
209  if (!geo)
210  {
211  std::cout << "Could not locate geometry node " << m_GeoNodeName << std::endl;
212  exit(1);
213  }
214  // loop over all of the layers in the hit container
215  // we need the geometry object for this layer
216  if (Verbosity() > 2) cout << " PHG4InttHitReco: Loop over hits" << endl;
217  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
218  for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
219  {
220  const int sphxlayer = hiter->second->get_detid();
221  CylinderGeomIntt *layergeom = dynamic_cast<CylinderGeomIntt *>(geo->GetLayerGeom(sphxlayer));
222 
223  // checking ADC timing integration window cut
224  // uses default values for now
225  // these should depend on layer radius
226  if (hiter->second->get_t(0) > m_Tmax)
227  continue;
228  if (hiter->second->get_t(1) < m_Tmin)
229  continue;
230 
231  // I made this (small) diffusion up for now, we will get actual values for the Intt later
232  double diffusion_width = 5.0e-04; // diffusion radius 5 microns, in cm
233 
234  const int ladder_z_index = hiter->second->get_ladder_z_index();
235  const int ladder_phi_index = hiter->second->get_ladder_phi_index();
236 
237  // What we have is a hit in the sensor. We have not yet assigned the strip(s) that were hit, we do that here
238  //========================================================================
239 
240  int strip_y_index_in, strip_z_index_in, strip_y_index_out, strip_z_index_out;
241  layergeom->find_strip_index_values(ladder_z_index, hiter->second->get_local_y(0), hiter->second->get_local_z(0), strip_y_index_in, strip_z_index_in);
242  layergeom->find_strip_index_values(ladder_z_index, hiter->second->get_local_y(1), hiter->second->get_local_z(1), strip_y_index_out, strip_z_index_out);
243 
244  if (Verbosity() > 5)
245  {
246  // check to see if we get back the positions from these strip index values
247  double check_location[3] = {-1, -1, -1};
248  layergeom->find_strip_center_localcoords(ladder_z_index, strip_y_index_in, strip_z_index_in, check_location);
249  cout << " G4 entry location = " << hiter->second->get_local_x(0) << " " << hiter->second->get_local_y(0) << " " << hiter->second->get_local_z(0) << endl;
250  cout << " Check entry location = " << check_location[0] << " " << check_location[1] << " " << check_location[2] << endl;
251  layergeom->find_strip_center_localcoords(ladder_z_index, strip_y_index_out, strip_z_index_out, check_location);
252  cout << " G4 exit location = " << hiter->second->get_local_x(1) << " " << hiter->second->get_local_y(1) << " " << hiter->second->get_local_z(1) << endl;
253  cout << " Check exit location = " << check_location[0] << " " << check_location[1] << " " << check_location[2] << endl;
254  }
255 
256  // Now we find how many strips were crossed by this track, and divide the energy between them
257  int minstrip_z = strip_z_index_in;
258  int maxstrip_z = strip_z_index_out;
259  if (minstrip_z > maxstrip_z) swap(minstrip_z, maxstrip_z);
260 
261  int minstrip_y = strip_y_index_in;
262  int maxstrip_y = strip_y_index_out;
263  if (minstrip_y > maxstrip_y) swap(minstrip_y, maxstrip_y);
264 
265  // Use an algorithm similar to the one for the MVTX pixels, since it facilitates adding charge diffusion
266  // for now we assume small charge diffusion
267  vector<int> vybin;
268  vector<int> vzbin;
269  //vector<double> vlen;
270  vector<pair<double, double> > venergy;
271 
272  //====================================================
273  // Beginning of charge sharing implementation
274  // Find tracklet line inside sensor
275  // Divide tracklet line into n segments (vary n until answer stabilizes)
276  // Find centroid of each segment
277  // Diffuse charge at each centroid
278  // Apportion charge between neighboring pixels
279  // Add the pixel energy contributions from different track segments together
280  //====================================================
281 
282  // skip this hit if it involves an unreasonable number of pixels
283  // this skips it if either the xbin or ybin range traversed is greater than 8 (for 8 adding two pixels at each end makes the range 12)
284  if (maxstrip_y - minstrip_y > 12 || maxstrip_z - minstrip_z > 12)
285  {
286  continue;
287  }
288  // this hit is skipped above if this dimensioning would be exceeded
289  double stripenergy[13][13] = {}; // init to 0
290  double stripeion[13][13] = {}; // init to 0
291 
292  int nsegments = 10;
293  // Loop over track segments and diffuse charge at each segment location, collect energy in pixels
294  // Get the entry point of the hit in sensor local coordinates
295  gsl_vector_set(m_PathVec, 0, hiter->second->get_local_x(0));
296  gsl_vector_set(m_PathVec, 1, hiter->second->get_local_y(0));
297  gsl_vector_set(m_PathVec, 2, hiter->second->get_local_z(0));
298  gsl_vector_set(m_LocalOutVec, 0, hiter->second->get_local_x(1));
299  gsl_vector_set(m_LocalOutVec, 1, hiter->second->get_local_y(1));
300  gsl_vector_set(m_LocalOutVec, 2, hiter->second->get_local_z(1));
301  gsl_vector_sub(m_PathVec, m_LocalOutVec);
302  for (int i = 0; i < nsegments; i++)
303  {
304  // Find the tracklet segment location
305  // If there are n segments of equal length, we want 2*n intervals
306  // The 1st segment is centered at interval 1, the 2nd at interval 3, the nth at interval 2n -1
307  double interval = 2 * (double) i + 1;
308  double frac = interval / (double) (2 * nsegments);
309  gsl_vector_memcpy(m_SegmentVec, m_PathVec);
310  gsl_vector_scale(m_SegmentVec, frac);
311  gsl_vector_add(m_SegmentVec, m_LocalOutVec);
312  // Caculate the charge diffusion over this drift distance
313  // increases from diffusion width_min to diffusion_width_max
314  double diffusion_radius = diffusion_width;
315 
316  if (Verbosity() > 5)
317  cout << " segment " << i
318  << " interval " << interval
319  << " frac " << frac
320  << " local_in.X " << hiter->second->get_local_x(0)
321  << " local_in.Z " << hiter->second->get_local_z(0)
322  << " local_in.Y " << hiter->second->get_local_y(0)
323  << " pathvec.X " << gsl_vector_get(m_PathVec, 0)
324  << " pathvec.Z " << gsl_vector_get(m_PathVec, 2)
325  << " pathvec.Y " << gsl_vector_get(m_PathVec, 1)
326  << " segvec.X " << gsl_vector_get(m_SegmentVec, 0)
327  << " segvec.Z " << gsl_vector_get(m_SegmentVec, 2)
328  << " segvec.Y " << gsl_vector_get(m_SegmentVec, 1) << endl
329  << " diffusion_radius " << diffusion_radius
330  << endl;
331 
332  // Now find the area of overlap of the diffusion circle with each pixel and apportion the energy
333  for (int iz = minstrip_z; iz <= maxstrip_z; iz++)
334  {
335  for (int iy = minstrip_y; iy <= maxstrip_y; iy++)
336  {
337  // Find the pixel corners for this pixel number
338  double location[3] = {-1, -1, -1};
339  layergeom->find_strip_center_localcoords(ladder_z_index, iy, iz, location);
340  // note that (y1,z1) is the top left corner, (y2,z2) is the bottom right corner of the pixel - circle_rectangle_intersection expects this ordering
341  double y1 = location[1] - layergeom->get_strip_y_spacing() / 2.0;
342  double y2 = location[1] + layergeom->get_strip_y_spacing() / 2.0;
343  double z1 = location[2] + layergeom->get_strip_z_spacing() / 2.0;
344  double z2 = location[2] - layergeom->get_strip_z_spacing() / 2.0;
345 
346  // here m_SegmentVec.1 (Y) and m_SegmentVec.2 (Z) are the center of the circle, and diffusion_radius is the circle radius
347  // circle_rectangle_intersection returns the overlap area of the circle and the pixel. It is very fast if there is no overlap.
348  double striparea_frac = PHG4Utils::circle_rectangle_intersection(y1, z1, y2, z2, gsl_vector_get(m_SegmentVec, 1), gsl_vector_get(m_SegmentVec, 2), diffusion_radius) / (M_PI * (diffusion_radius * diffusion_radius));
349  // assume that the energy is deposited uniformly along the tracklet length, so that this segment gets the fraction 1/nsegments of the energy
350  stripenergy[iy - minstrip_y][iz - minstrip_z] += striparea_frac * hiter->second->get_edep() / (float) nsegments;
351  if (hiter->second->has_property(PHG4Hit::prop_eion))
352  {
353  stripeion[iy - minstrip_y][iz - minstrip_z] += striparea_frac * hiter->second->get_eion() / (float) nsegments;
354  }
355  if (Verbosity() > 5)
356  {
357  cout << " strip y index " << iy << " strip z index " << iz
358  << " strip area fraction of circle " << striparea_frac << " accumulated pixel energy " << stripenergy[iy - minstrip_y][iz - minstrip_z]
359  << endl;
360  }
361  }
362  }
363  } // end loop over segments
364  // now we have the energy deposited in each pixel, summed over all tracklet segments. We make a vector of all pixels with non-zero energy deposited
365 
366  for (int iz = minstrip_z; iz <= maxstrip_z; iz++)
367  {
368  for (int iy = minstrip_y; iy <= maxstrip_y; iy++)
369  {
370  if (stripenergy[iy - minstrip_y][iz - minstrip_z] > 0.0)
371  {
372  vybin.push_back(iy);
373  vzbin.push_back(iz);
374  pair<double, double> tmppair = make_pair(stripenergy[iy - minstrip_y][iz - minstrip_z], stripeion[iy - minstrip_y][iz - minstrip_z]);
375  venergy.push_back(tmppair);
376  if (Verbosity() > 1)
377  cout << " Added ybin " << iy << " zbin " << iz << " to vectors with energy " << stripenergy[iy - minstrip_y][iz - minstrip_z] << endl;
378  }
379  }
380  }
381 
382  //===================================
383  // End of charge sharing implementation
384  //===================================
385 
386  for (unsigned int i1 = 0; i1 < vybin.size(); i1++) // loop over all fired cells
387  {
388  // We add the Intt TrkrHitsets directly to the node using hitsetcontainer
389 
390  // We need to create the TrkrHitSet if not already made - each TrkrHitSet should correspond to a sensor for the Intt ?
391  // The hitset key includes the layer, the ladder_z_index (sensors numbered 0-3) and ladder_phi_index (azimuthal location of ladder) for this hit
392  TrkrDefs::hitsetkey hitsetkey = InttDefs::genHitSetKey(sphxlayer, ladder_z_index, ladder_phi_index);
393  // Use existing hitset or add new one if needed
394  TrkrHitSetContainer::Iterator hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey);
395 
396  // generate the key for this hit
397  TrkrDefs::hitkey hitkey = InttDefs::genHitKey(vzbin[i1], vybin[i1]);
398  // See if this hit already exists
399  TrkrHit *hit = hitsetit->second->getHit(hitkey);
400  if (!hit)
401  {
402  // Otherwise, create a new one
403  //hit = new InttHit();
404  hit = new TrkrHitv2();
405  hitsetit->second->addHitSpecificKey(hitkey, hit);
406  }
407 
408  // Either way, add the energy to it
409  if (Verbosity() > 2)
410  cout << "add energy " << venergy[i1].first << " to intthit " << endl;
411  hit->addEnergy(venergy[i1].first * TrkrDefs::InttEnergyScaleup);
412 
413  // Add this hit to the association map
414  hittruthassoc->addAssoc(hitsetkey, hitkey, hiter->first);
415 
416  if (Verbosity() > 2)
417  cout << "PHG4InttHitReco: added hit wirh hitsetkey " << hitsetkey << " hitkey " << hitkey << " g4hitkey " << hiter->first << " energy " << hit->getEnergy() << endl;
418  }
419 
420  } // end loop over g4hits
421 
422  // print the list of entries in the association table
423  if (Verbosity() > 0)
424  {
425  cout << "From PHG4InttHitReco: " << endl;
426  hitsetcontainer->identify();
427  hittruthassoc->identify();
428  }
430 }
431 
433 {
434  // if we ever need separate timing windows, don't patch around here!
435  // use PHParameterContainerInterface which
436  // provides for multiple layers/detector types
437  set_default_double_param("tmax", 80.0); // FVTX NIM paper Fig 32
438  set_default_double_param("tmin", -20.0); // FVTX NIM paper Fig 32, collision has a timing spread around the triggered event. Accepting negative time too.
439  return;
440 }