ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MvtxHitReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MvtxHitReco.cc
1 // this is the new trackbase version
2 
3 #include "PHG4MvtxHitReco.h"
4 
6 #include <mvtx/MvtxDefs.h>
7 
8 #include <trackbase/TrkrDefs.h>
9 #include <trackbase/TrkrHitv2.h> // for TrkrHit
10 #include <trackbase/TrkrHitSet.h>
13 
14 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
16 
17 #include <phparameter/PHParameterContainerInterface.h>
18 
19 #include <g4main/PHG4Hit.h>
21 #include <g4main/PHG4Utils.h>
22 
24 #include <fun4all/SubsysReco.h> // for SubsysReco
25 
26 #include <phool/PHCompositeNode.h>
27 #include <phool/PHIODataNode.h>
28 #include <phool/PHNode.h> // for PHNode
29 #include <phool/PHNodeIterator.h>
30 #include <phool/PHObject.h> // for PHObject
31 #include <phool/getClass.h>
32 #include <phool/phool.h> // for PHWHERE
33 
34 #include <TVector3.h> // for TVector3, ope...
35 
36 #include <cmath>
37 #include <cstdlib>
38 #include <iostream>
39 #include <memory> // for allocator_tra...
40 #include <vector> // for vector
41 
42 using namespace std;
43 
45  : SubsysReco(name)
47  , detector(name)
48 {
49  SetDefaultParameters(); // sets default timing window
50  if (Verbosity() > 0)
51  cout << "Creating PHG4MvtxHitReco for name = " << name << endl;
52 }
53 
55 {
56  PHNodeIterator iter(topNode);
57 
58  // Looking for the DST node
59  PHCompositeNode *dstNode;
60  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
61  if (!dstNode)
62  {
63  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
64  exit(1);
65  }
66  hitnodename = "G4HIT_" + detector;
67  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
68  if (!g4hit)
69  {
70  cout << "Could not locate g4 hit node " << hitnodename << endl;
71  exit(1);
72  }
73 
74  geonodename = "CYLINDERGEOM_" + detector;
75  PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
76  if (!geo)
77  {
78  cout << "Could not locate geometry node " << geonodename << endl;
79  exit(1);
80  }
81  if (Verbosity() > 0)
82  {
83  geo->identify();
84  }
85 
86  auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
87  if (!hitsetcontainer)
88  {
89  PHNodeIterator dstiter(dstNode);
90  PHCompositeNode *DetNode =
91  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
92  if (!DetNode)
93  {
94  DetNode = new PHCompositeNode("TRKR");
95  dstNode->addNode(DetNode);
96  }
97 
98  hitsetcontainer = new TrkrHitSetContainerv1;
99  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
100  DetNode->addNode(newNode);
101  }
102 
103  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
104  if (!hittruthassoc)
105  {
106  PHNodeIterator dstiter(dstNode);
107  PHCompositeNode *DetNode =
108  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
109  if (!DetNode)
110  {
111  DetNode = new PHCompositeNode("TRKR");
112  dstNode->addNode(DetNode);
113  }
114 
115  hittruthassoc = new TrkrHitTruthAssocv1;
116  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
117  DetNode->addNode(newNode);
118  }
119 
121 }
122 
124 {
125  //cout << PHWHERE << "Entering process_event for PHG4MvtxHitReco" << endl;
126 
127  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
128  if (!g4hit)
129  {
130  cout << "Could not locate g4 hit node " << hitnodename << endl;
131  exit(1);
132  }
133 
134  PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
135  if (!geo)
136  {
137  cout << "Could not locate geometry node " << geonodename << endl;
138  exit(1);
139  }
140 
141  // Get the TrkrHitSetContainer node
142  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
143  if (!trkrhitsetcontainer)
144  {
145  cout << "Could not locate TRKR_HITSET node, quit! " << endl;
146  exit(1);
147  }
148 
149  // Get the TrkrHitTruthAssoc node
150  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
151  if (!hittruthassoc)
152  {
153  cout << "Could not locate TRKR_HITTRUTHASSOC node, quit! " << endl;
154  exit(1);
155  }
156 
157  // loop over all of the layers in the g4hit container
159  pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
160  for (layer = layer_begin_end.first; layer != layer_begin_end.second; ++layer)
161  {
162  //cout << "---------- PHG4MvtxHitReco: Looping over layers " << endl;
163 
164  // loop over the hits in this layer
166  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
167 
168  // we need the geometry object for this layer
169  CylinderGeom_Mvtx *layergeom = dynamic_cast<CylinderGeom_Mvtx *>(geo->GetLayerGeom(*layer));
170  if (!layergeom)
171  exit(1);
172 
173  if (Verbosity() > 2)
174  layergeom->identify();
175 
176  // Get some layer parameters for later use
177  double xpixw = layergeom->get_pixel_x();
178  double xpixw_half = xpixw / 2.0;
179  double zpixw = layergeom->get_pixel_z();
180  double zpixw_half = zpixw / 2.0;
181  int maxNX = layergeom->get_NX();
182  int maxNZ = layergeom->get_NZ();
183 
184  // Now loop over all g4 hits for this layer
185  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
186  {
187  //cout << "From PHG4MvtxHitReco: Call hit print method: " << endl;
188  if (Verbosity() > 4)
189  hiter->second->print();
190 
191  // checking ADC timing integration window cut
192  if (*layer > 2)
193  {
194  cout << PHWHERE << "Mvtx layers only go up to three! Quit." << endl;
195  exit(1);
196  }
197  if (Verbosity() > 1)
198  cout << " layer " << *layer << " t0 " << hiter->second->get_t(0) << " t1 " << hiter->second->get_t(1)
199  << " tmin " << tmin_max[*layer].first << " tmax " << tmin_max[*layer].second
200  << endl;
201  if (hiter->second->get_t(0) > tmin_max[*layer].second) continue;
202  if (hiter->second->get_t(1) < tmin_max[*layer].first) continue;
203 
204  // get_property_int(const PROPERTY prop_id) const {return INT_MIN;}
205  int stave_number = hiter->second->get_property_int(PHG4Hit::prop_stave_index);
206  int half_stave_number = hiter->second->get_property_int(PHG4Hit::prop_half_stave_index);
207  int module_number = hiter->second->get_property_int(PHG4Hit::prop_module_index);
208  int chip_number = hiter->second->get_property_int(PHG4Hit::prop_chip_index);
209 
210  TVector3 local_in(hiter->second->get_local_x(0), hiter->second->get_local_y(0), hiter->second->get_local_z(0));
211  TVector3 local_out(hiter->second->get_local_x(1), hiter->second->get_local_y(1), hiter->second->get_local_z(1));
212  TVector3 midpoint((local_in.X() + local_out.X()) / 2.0, (local_in.Y() + local_out.Y()) / 2.0, (local_in.Z() + local_out.Z()) / 2.0);
213 
214  if (Verbosity() > 4)
215  {
216  cout << endl
217  << " world entry point position: " << hiter->second->get_x(0) << " " << hiter->second->get_y(0) << " " << hiter->second->get_z(0) << endl;
218  cout << " world exit point position: " << hiter->second->get_x(1) << " " << hiter->second->get_y(1) << " " << hiter->second->get_z(1) << endl;
219  cout << " local coords of entry point from G4 " << hiter->second->get_local_x(0) << " " << hiter->second->get_local_y(0) << " " << hiter->second->get_local_z(0) << endl;
220  TVector3 world_in(hiter->second->get_x(0), hiter->second->get_y(0), hiter->second->get_z(0));
221  TVector3 local_in_check = layergeom->get_local_from_world_coords(stave_number, half_stave_number, module_number, chip_number, world_in);
222  cout << " local coords of entry point from geom (a check) " << local_in_check.X() << " " << local_in_check.Y() << " " << local_in_check.Z() << endl;
223  cout << " local coords of exit point from G4 " << hiter->second->get_local_x(1) << " " << hiter->second->get_local_y(1) << " " << hiter->second->get_local_z(1) << endl;
224  cout << " local coords of exit point from geom (a check) " << local_out.X() << " " << local_out.Y() << " " << local_out.Z() << endl;
225  cout << endl;
226  }
227 
228  if (Verbosity() > 2)
229  {
230  // As a check, get the positions of the hit pixels in world coordinates from the geo object
231  TVector3 location_in = layergeom->get_world_from_local_coords(stave_number, half_stave_number, module_number, chip_number, local_in);
232  TVector3 location_out = layergeom->get_world_from_local_coords(stave_number, half_stave_number, module_number, chip_number, local_out);
233 
234  cout << endl
235  << " PHG4MvtxHitReco: Found world entry location from geometry for "
236  << " stave number " << stave_number
237  << " half stave number " << half_stave_number
238  << " module number " << module_number
239  << " chip number " << chip_number
240  << endl
241  << " x = " << location_in.X()
242  << " y = " << location_in.Y()
243  << " z = " << location_in.Z()
244  << " radius " << sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
245  << " angle " << atan(location_in.Y() / location_in.X())
246  << endl;
247  cout << " PHG4MvtxHitReco: The world entry location from G4 was "
248  << endl
249  << " x = " << hiter->second->get_x(0)
250  << " y = " << hiter->second->get_y(0)
251  << " z = " << hiter->second->get_z(0)
252  << " radius " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2))
253  << " angle " << atan(hiter->second->get_y(0) / hiter->second->get_x(0))
254  << endl;
255  cout << " difference in x = " << hiter->second->get_x(0) - location_in.X()
256  << " difference in y = " << hiter->second->get_y(0) - location_in.Y()
257  << " difference in z = " << hiter->second->get_z(0) - location_in.Z()
258  << " difference in radius = " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2)) - sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
259  << " in angle = " << atan(hiter->second->get_y(0) / hiter->second->get_x(0)) - atan(location_in.Y() / location_in.X())
260  << endl
261  << endl;
262 
263  cout << " PHG4MvtxHitReco: Found world exit location from geometry for "
264  << " stave number " << stave_number
265  << " half stave number " << half_stave_number
266  << " module number" << module_number
267  << endl
268  << " x = " << location_out.X()
269  << " y = " << location_out.Y()
270  << " z = " << location_out.Z()
271  << " radius " << sqrt(pow(location_out.X(), 2) + pow(location_out.Y(), 2))
272  << " angle " << atan(location_out.Y() / location_out.X())
273  << endl;
274  cout << " PHG4MvtxHitReco: The world exit location from G4 was "
275  << endl
276  << " x = " << hiter->second->get_x(1)
277  << " y = " << hiter->second->get_y(1)
278  << " z = " << hiter->second->get_z(1)
279  << " radius " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2))
280  << " angle " << atan(hiter->second->get_y(1) / hiter->second->get_x(1))
281  << endl;
282  cout << " difference in radius = " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2)) - sqrt(pow(location_out.X(), 2) + pow(location_out.Y(), 2))
283  << " in angle = " << atan(hiter->second->get_y(1) / hiter->second->get_x(1)) - atan(location_out.Y() / location_out.X())
284  << endl
285  << endl;
286  }
287 
288  // Get the pixel number of the entry location
289  int pixel_number_in = layergeom->get_pixel_from_local_coords(local_in);
290  // Get the pixel number of the exit location
291  int pixel_number_out = layergeom->get_pixel_from_local_coords(local_out);
292 
293  if (pixel_number_in < 0 || pixel_number_out < 0)
294  {
295  cout << "Oops! got negative pixel number in layer " << layergeom->get_layer()
296  << " pixel_number_in " << pixel_number_in
297  << " pixel_number_out " << pixel_number_out
298  << " local_in = " << local_in.X() << " " << local_in.Y() << " " << local_in.Z()
299  << " local_out = " << local_out.X() << " " << local_out.Y() << " " << local_out.Z()
300  << endl;
301  }
302 
303  if (Verbosity() > 0)
304  cout << "entry pixel number " << pixel_number_in << " exit pixel number " << pixel_number_out << endl;
305 
306  vector<int> vpixel;
307  vector<int> vxbin;
308  vector<int> vzbin;
309  vector<pair<double, double> > venergy;
310  //double trklen = 0.0;
311 
312  //===================================================
313  // OK, now we have found which sensor the hit is in, extracted the hit
314  // position in local sensor coordinates, and found the pixel numbers of the
315  // entry point and exit point
316 
317  //====================================================
318  // Beginning of charge sharing implementation
319  // Find tracklet line inside sensor
320  // Divide tracklet line into n segments (vary n until answer stabilizes)
321  // Find centroid of each segment
322  // Diffuse charge at each centroid
323  // Apportion charge between neighboring pixels
324  // Add the pixel energy contributions from different track segments together
325  //====================================================
326 
327  TVector3 pathvec = local_in - local_out;
328 
329  // See figure 7.3 of the thesis by Lucasz Maczewski (arXiv:10053.3710) for diffusion simulations in a MAPS epitaxial layer
330  // The diffusion widths below were inspired by those plots, corresponding to where the probability drops off to 1/3 of the peak value
331  // However note that we make the simplifying assumption that the probability distribution is flat within this diffusion width,
332  // while in the simulation it is not
333  //double diffusion_width_max = 35.0e-04; // maximum diffusion radius 35 microns, in cm
334  //double diffusion_width_min = 12.0e-04; // minimum diffusion radius 12 microns, in cm
335  double diffusion_width_max = 25.0e-04; // maximum diffusion radius 35 microns, in cm
336  double diffusion_width_min = 8.0e-04; // minimum diffusion radius 12 microns, in cm
337 
338  double ydrift_max = pathvec.Y();
339  int nsegments = 4;
340 
341  // we want to make a list of all pixels possibly affected by this hit
342  // we take the entry and exit locations in local coordinates, and build
343  // a rectangular array of pixels that encompasses both, with "nadd" pixels added all around
344 
345  int xbin_in = layergeom->get_pixel_X_from_pixel_number(pixel_number_in);
346  int zbin_in = layergeom->get_pixel_Z_from_pixel_number(pixel_number_in);
347  int xbin_out = layergeom->get_pixel_X_from_pixel_number(pixel_number_out);
348  int zbin_out = layergeom->get_pixel_Z_from_pixel_number(pixel_number_out);
349 
350  int xbin_max, xbin_min;
351  int nadd = 2;
352  if (xbin_in > xbin_out)
353  {
354  xbin_max = xbin_in + nadd;
355  xbin_min = xbin_out - nadd;
356  }
357  else
358  {
359  xbin_max = xbin_out + nadd;
360  xbin_min = xbin_in - nadd;
361  }
362 
363  int zbin_max, zbin_min;
364  if (zbin_in > zbin_out)
365  {
366  zbin_max = zbin_in + nadd;
367  zbin_min = zbin_out - nadd;
368  }
369  else
370  {
371  zbin_max = zbin_out + nadd;
372  zbin_min = zbin_in - nadd;
373  }
374 
375  // need to check that values of xbin and zbin are within the valid range
376  // YCM: Fix pixel range: Xbin (row) 0 to 511, Zbin (col) 0 to 1023
377  if (xbin_min < 0) xbin_min = 0;
378  if (zbin_min < 0) zbin_min = 0;
379  if (xbin_max >= maxNX) xbin_max = maxNX-1;
380  if (zbin_max >= maxNZ) xbin_max = maxNZ-1;
381 
382  if (Verbosity() > 1)
383  {
384  cout << " xbin_in " << xbin_in << " xbin_out " << xbin_out << " xbin_min " << xbin_min << " xbin_max " << xbin_max << endl;
385  cout << " zbin_in " << zbin_in << " zbin_out " << zbin_out << " zbin_min " << zbin_min << " zbin_max " << zbin_max << endl;
386  }
387 
388  // skip this hit if it involves an unreasonable number of pixels
389  // 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)
390  if (xbin_max - xbin_min > 12 || zbin_max - zbin_min > 12)
391  continue;
392 
393  // this hit is skipped earlier if this dimensioning would be exceeded
394  double pixenergy[12][12] = {}; // init to 0
395  double pixeion[12][12] = {}; // init to 0
396 
397  // Loop over track segments and diffuse charge at each segment location, collect energy in pixels
398  for (int i = 0; i < nsegments; i++)
399  {
400  // Find the tracklet segment location
401  // If there are n segments of equal length, we want 2*n intervals
402  // The 1st segment is centered at interval 1, the 2nd at interval 3, the nth at interval 2n -1
403  double interval = 2 * (double) i + 1;
404  double frac = interval / (double) (2 * nsegments);
405  TVector3 segvec(pathvec.X() * frac, pathvec.Y() * frac, pathvec.Z() * frac);
406  segvec = segvec + local_out;
407 
408  // Find the distance to the back of the sensor from the segment location
409  // That projection changes only the value of y
410  double ydrift = segvec.Y() - local_out.Y();
411 
412  // Caculate the charge diffusion over this drift distance
413  // increases from diffusion width_min to diffusion_width_max
414  double ydiffusion_radius = diffusion_width_min + (ydrift / ydrift_max) * (diffusion_width_max - diffusion_width_min);
415 
416  if (Verbosity() > 5)
417  cout << " segment " << i
418  << " interval " << interval
419  << " frac " << frac
420  << " local_in.X " << local_in.X()
421  << " local_in.Z " << local_in.Z()
422  << " local_in.Y " << local_in.Y()
423  << " pathvec.X " << pathvec.X()
424  << " pathvec.Z " << pathvec.Z()
425  << " pathvec.Y " << pathvec.Y()
426  << " segvec.X " << segvec.X()
427  << " segvec.Z " << segvec.Z()
428  << " segvec.Y " << segvec.Y()
429  << " ydrift " << ydrift
430  << " ydrift_max " << ydrift_max
431  << " ydiffusion_radius " << ydiffusion_radius
432  << endl;
433 
434  // Now find the area of overlap of the diffusion circle with each pixel and apportion the energy
435  for (int ix = xbin_min; ix <= xbin_max; ix++)
436  {
437  for (int iz = zbin_min; iz <= zbin_max; iz++)
438  {
439  // Find the pixel corners for this pixel number
440  int pixnum = layergeom->get_pixel_number_from_xbin_zbin(ix, iz);
441 
442  if (pixnum < 0)
443  {
444  cout << " pixnum < 0 , pixnum = " << pixnum << endl;
445  cout << " ix " << ix << " iz " << iz << endl;
446  cout << " xbin_min " << xbin_min << " zbin_min " << zbin_min
447  << " xbin_max " << xbin_max << " zbin_max " << zbin_max
448  << endl;
449  cout << " maxNX " << maxNX << " maxNZ " << maxNZ
450  << endl;
451  }
452 
453  TVector3 tmp = layergeom->get_local_coords_from_pixel(pixnum);
454  // note that (x1,z1) is the top left corner, (x2,z2) is the bottom right corner of the pixel - circle_rectangle_intersection expects this ordering
455  double x1 = tmp.X() - xpixw_half;
456  double z1 = tmp.Z() + zpixw_half;
457  double x2 = tmp.X() + xpixw_half;
458  double z2 = tmp.Z() - zpixw_half;
459 
460  // here segvec.X and segvec.Z are the center of the circle, and diffusion_radius is the circle radius
461  // circle_rectangle_intersection returns the overlap area of the circle and the pixel. It is very fast if there is no overlap.
462  double pixarea_frac = PHG4Utils::circle_rectangle_intersection(x1, z1, x2, z2, segvec.X(), segvec.Z(), ydiffusion_radius) / (M_PI * pow(ydiffusion_radius, 2));
463  // assume that the energy is deposited uniformly along the tracklet length, so that this segment gets the fraction 1/nsegments of the energy
464  pixenergy[ix - xbin_min][iz - zbin_min] += pixarea_frac * hiter->second->get_edep() / (float) nsegments;
465  if (hiter->second->has_property(PHG4Hit::prop_eion))
466  {
467  pixeion[ix - xbin_min][iz - zbin_min] += pixarea_frac * hiter->second->get_eion() / (float) nsegments;
468  }
469  if (Verbosity() > 5)
470  {
471  cout << " pixnum " << pixnum << " xbin " << ix << " zbin " << iz
472  << " pixel_area fraction of circle " << pixarea_frac << " accumulated pixel energy " << pixenergy[ix - xbin_min][iz - zbin_min]
473  << endl;
474  }
475  }
476  }
477  } // end loop over segments
478 
479  // 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
480  for (int ix = xbin_min; ix <= xbin_max; ix++)
481  {
482  for (int iz = zbin_min; iz <= zbin_max; iz++)
483  {
484  if (pixenergy[ix - xbin_min][iz - zbin_min] > 0.0)
485  {
486  int pixnum = layergeom->get_pixel_number_from_xbin_zbin(ix, iz);
487  vpixel.push_back(pixnum);
488  vxbin.push_back(ix);
489  vzbin.push_back(iz);
490  pair<double, double> tmppair = make_pair(pixenergy[ix - xbin_min][iz - zbin_min], pixeion[ix - xbin_min][iz - zbin_min]);
491  venergy.push_back(tmppair);
492  if (Verbosity() > 1)
493  cout << " Added pixel number " << pixnum << " xbin " << ix << " zbin " << iz << " to vectors with energy " << pixenergy[ix - xbin_min][iz - zbin_min] << endl;
494  }
495  }
496  }
497 
498  //===================================
499  // End of charge sharing implementation
500  //===================================
501 
502  // loop over all fired cells for this g4hit and add them to the TrkrHitSet
503  for (unsigned int i1 = 0; i1 < vpixel.size(); i1++) // loop over all fired cells
504  {
505  // This is the new storage object version
506  //====================================
507 
508  // We need to create the TrkrHitSet if not already made - each TrkrHitSet should correspond to a chip for the Mvtx
509  TrkrDefs::hitsetkey hitsetkey = MvtxDefs::genHitSetKey(*layer, stave_number, chip_number);
510  // Use existing hitset or add new one if needed
511  TrkrHitSetContainer::Iterator hitsetit = trkrhitsetcontainer->findOrAddHitSet(hitsetkey);
512 
513  // generate the key for this hit
514  TrkrDefs::hitkey hitkey = MvtxDefs::genHitKey(vzbin[i1], vxbin[i1]);
515  // See if this hit already exists
516  TrkrHit *hit = nullptr;
517  hit = hitsetit->second->getHit(hitkey);
518  if (!hit)
519  {
520  // Otherwise, create a new one
521  hit = new TrkrHitv2();
522  hitsetit->second->addHitSpecificKey(hitkey, hit);
523  }
524 
525  // Either way, add the energy to it
526  hit->addEnergy(venergy[i1].first * TrkrDefs::MvtxEnergyScaleup);
527 
528  // now we update the TrkrHitTruthAssoc map - the map contains <hitsetkey, std::pair <hitkey, g4hitkey> >
529  // There is only one TrkrHit per pixel, but there may be multiple g4hits
530  // How do we know how much energy from PHG4Hit went into TrkrHit? We don't, have to sort it out in evaluator to save memory
531 
532  // How do we check if this association already exists?
533  //cout << "PHG4MvtxHitReco: adding association entry for hitkey " << hitkey << " and g4hitkey " << hiter->first << endl;
534  hittruthassoc->addAssoc(hitsetkey, hitkey, hiter->first);
535 
536  } // end loop over hit cells
537  } // end loop over g4hits for this layer
538 
539  } // end loop over layers
540 
541  // print the list of entries in the association table
542  if (Verbosity() > 2)
543  {
544  cout << "From PHG4MvtxHitReco: " << endl;
545  hittruthassoc->identify();
546  }
547 
549 }
550 
551 void PHG4MvtxHitReco::set_timing_window(const int detid, const double tmin, const double tmax)
552 {
553  // first have to erase the default value
554  std::map<int, std::pair<double, double> >::iterator it = tmin_max.find(detid);
555  tmin_max.erase(it);
556  // now replace it with the new value
557  tmin_max.insert(std::make_pair(detid, std::make_pair(tmin, tmax)));
558 
559  cout << "PHG4MvtxHitReco: Set Mvtx timing window parameters from macro for layer = " << detid << " to tmin = " << tmin_max[detid].first << " tmax = " << tmin_max[detid].second << endl;
560 
561  return;
562 }
563 
565 {
566  cout << "PHG4MvtxHitReco: Setting Mvtx timing window defaults to tmin = -5000 and tmax = 5000 ns" << endl;
567  for (int ilayer = 0; ilayer < 3; ilayer++)
568  {
569  tmin_max.insert(std::make_pair(ilayer, std::make_pair(-5000, 5000)));
570  }
571  return;
572 }