ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcPadPlaneReadout.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcPadPlaneReadout.cc
2 
3 #include <g4detectors/PHG4Cell.h> // for PHG4Cell
5 #include <g4detectors/PHG4CellDefs.h> // for genkey, keytype
9 
10 #include <g4main/PHG4Hit.h> // for PHG4Hit
12 
13 #include <phool/PHRandomSeed.h>
14 
15 // Move to new storage containers
16 #include <trackbase/TrkrDefs.h> // for hitkey, hitse...
17 #include <trackbase/TrkrHit.h> // for TrkrHit
18 #include <trackbase/TrkrHitSet.h>
20 #include <trackbase/TrkrHitv2.h> // for TrkrHit
21 
22 #include <tpc/TpcDefs.h>
23 
24 #include <TSystem.h>
25 
26 #include <gsl/gsl_randist.h>
27 #include <gsl/gsl_rng.h> // for gsl_rng_alloc
28 
29 #include <cmath>
30 #include <iostream>
31 #include <map> // for _Rb_tree_cons...
32 #include <utility> // for pair
33 
34 class PHCompositeNode;
35 class TrkrHitTruthAssoc;
36 
37 namespace
38 {
40  template <class T>
41  inline constexpr T square(const T &x)
42  {
43  return x * x;
44  }
45 
47  template <class T>
48  inline T gaus(const T &x, const T &sigma)
49  {
50  return std::exp(-square(x / sigma) / 2) / (sigma * std::sqrt(2 * M_PI));
51  }
52 
53 } // namespace
54 
56  : PHG4TpcPadPlane(name)
57 {
59 
60  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
61  gsl_rng_set(RandomGenerator, PHRandomSeed()); // fixed seed is handled in this funtcion
62 
63  return;
64 }
65 
67 {
68  gsl_rng_free(RandomGenerator);
69 }
70 
72 {
73  if (Verbosity()) std::cout << "PHG4TpcPadPlaneReadout: CreateReadoutGeometry: " << std::endl;
74 
75  for (int iregion = 0; iregion < 3; ++iregion)
76  {
77  for (int layer = MinLayer[iregion]; layer < MinLayer[iregion] + NTpcLayers[iregion]; ++layer)
78  {
79  if (Verbosity())
80  {
81  std::cout << " layer " << layer << " MinLayer " << MinLayer[iregion] << " region " << iregion
82  << " radius " << MinRadius[iregion] + ((double) (layer - MinLayer[iregion]) + 0.5) * Thickness[iregion]
83  << " thickness " << Thickness[iregion]
84  << " NZbins " << NZBins << " zmin " << MinZ << " zstep " << ZBinWidth
85  << " phibins " << NPhiBins[iregion] << " phistep " << PhiBinWidth[iregion] << std::endl;
86  }
87 
88  PHG4CylinderCellGeom *layerseggeo = new PHG4CylinderCellGeom();
89  layerseggeo->set_layer(layer);
90  layerseggeo->set_radius(MinRadius[iregion] + ((double) (layer - MinLayer[iregion]) + 0.5) * Thickness[iregion]);
91  layerseggeo->set_thickness(Thickness[iregion]);
93  layerseggeo->set_zbins(NZBins);
94  layerseggeo->set_zmin(MinZ);
95  layerseggeo->set_zstep(ZBinWidth);
96  layerseggeo->set_phibins(NPhiBins[iregion]);
97  layerseggeo->set_phistep(PhiBinWidth[iregion]);
98  // Chris Pinkenburg: greater causes huge memory growth which causes problems
99  // on our farm. If you need to increase this - TALK TO ME first
100  if (NPhiBins[iregion] * NZBins > 5100000)
101  {
102  std::cout << "increase Tpc cellsize, number of cells "
103  << NPhiBins[iregion] * NZBins << " for layer " << layer
104  << " exceed 5.1M limit" << std::endl;
105  gSystem->Exit(1);
106  }
107  seggeo->AddLayerCellGeom(layerseggeo);
108  }
109  }
110 
111  GeomContainer = seggeo;
112 
113  return 0;
114 }
115 
116 // This is obsolete, it uses the old PHG4Cell containers
117 void PHG4TpcPadPlaneReadout::MapToPadPlane(PHG4CellContainer *g4cells, const double x_gem, const double y_gem, const double z_gem, PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/)
118 {
119  // One electron per call of this method
120  // The x_gem and y_gem values have already been randomized within the transverse drift diffusion width
121  // The z_gem value already reflects the drift time of the primary electron from the production point, and is randomized within the longitudinal diffusion witdth
122 
123  double phi = atan2(y_gem, x_gem);
124  if (phi > +M_PI) phi -= 2 * M_PI;
125  if (phi < -M_PI) phi += 2 * M_PI;
126 
127  rad_gem = sqrt(x_gem * x_gem + y_gem * y_gem);
128  //std::cout << "Enter old MapToPadPlane with rad_gem " << rad_gem << std::endl;
129 
130  unsigned int layernum = 0;
131 
132  // Find which readout layer this electron ends up in
133 
135  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter = layerrange.first;
136  layeriter != layerrange.second;
137  ++layeriter)
138  {
139  double rad_low = layeriter->second->get_radius() - layeriter->second->get_thickness() / 2.0;
140  double rad_high = layeriter->second->get_radius() + layeriter->second->get_thickness() / 2.0;
141 
142  if (rad_gem > rad_low && rad_gem < rad_high)
143  {
144  // capture the layer where this electron hits sthe gem stack
145  LayerGeom = layeriter->second;
146  layernum = LayerGeom->get_layer();
147  if (Verbosity() > 1000 && layernum == print_layer)
148  std::cout << " g4hit id " << hiter->first << " rad_gem " << rad_gem << " rad_low " << rad_low << " rad_high " << rad_high
149  << " layer " << hiter->second->get_layer() << " want to change to " << layernum << std::endl;
150  hiter->second->set_layer(layernum); // have to set here, since the stepping action knows nothing about layers
151  }
152  }
153 
154  if (layernum == 0)
155  {
156  return;
157  }
158 
159  // store phi bins and zbins upfront to avoid repetitive checks on the phi methods
160  const auto phibins = LayerGeom->get_phibins();
161  const auto zbins = LayerGeom->get_zbins();
162 
163  // Create the distribution function of charge on the pad plane around the electron position
164 
165  // The resolution due to pad readout includes the charge spread during GEM multiplication.
166  // this now defaults to 400 microns during construction from Tom (see 8/11 email).
167  // Use the setSigmaT(const double) method to update...
168  // We use a double gaussian to represent the smearing due to the SAMPA chip shaping time - default values of fShapingLead and fShapingTail are for 80 ns SAMPA
169 
170  // amplify the single electron in the gem stack
171  //===============================
172 
173  double nelec = getSingleEGEMAmplification();
174 
175  // Distribute the charge between the pads in phi
176  //====================================
177 
178  if (Verbosity() > 200)
179  std::cout << " populate phi bins for "
180  << " layernum " << layernum
181  << " phi " << phi
182  << " sigmaT " << sigmaT
183  << " zigzag_pads " << zigzag_pads
184  << std::endl;
185 
186  pad_phibin.clear();
187  pad_phibin_share.clear();
188  if (zigzag_pads)
190  else
192 
193  // Normalize the shares so they add up to 1
194  double norm1 = 0.0;
195  for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad)
196  {
197  double pad_share = pad_phibin_share[ipad];
198  norm1 += pad_share;
199  }
200  for (unsigned int iphi = 0; iphi < pad_phibin.size(); ++iphi)
201  pad_phibin_share[iphi] /= norm1;
202 
203  // Distribute the charge between the pads in z
204  //====================================
205  if (Verbosity() > 100 && layernum == print_layer)
206  std::cout << " populate z bins for layernum " << layernum
207  << " with z_gem " << z_gem << " sigmaL[0] " << sigmaL[0] << " sigmaL[1] " << sigmaL[1] << std::endl;
208 
209  adc_zbin.clear();
210  adc_zbin_share.clear();
212 
213  // Normalize the shares so that they add up to 1
214  double znorm = 0.0;
215  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
216  {
217  double bin_share = adc_zbin_share[iz];
218  znorm += bin_share;
219  }
220  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
221  adc_zbin_share[iz] /= znorm;
222 
223  // Fill cells
224  //========
225  // These are used to do a quick clustering for checking
226  double phi_integral = 0.0;
227  double z_integral = 0.0;
228  double weight = 0.0;
229 
230  for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad)
231  {
232  int pad_num = pad_phibin[ipad];
233  double pad_share = pad_phibin_share[ipad];
234 
235  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
236  {
237  int zbin_num = adc_zbin[iz];
238  double adc_bin_share = adc_zbin_share[iz];
239 
240  // Divide electrons from avalanche between bins
241  float neffelectrons = nelec * (pad_share) * (adc_bin_share);
242  if (neffelectrons < neffelectrons_threshold) continue; // skip signals that will be below the noise suppression threshold
243 
244  if (zbin_num >= zbins) std::cout << " Error making key: adc_zbin " << zbin_num << " nzbins " << zbins << std::endl;
245  if (pad_num >= phibins) std::cout << " Error making key: pad_phibin " << pad_num << " nphibins " << phibins << std::endl;
246 
247  // collect information to do simple clustering. Checks operation of PHG4CylinderCellTpcReco, and
248  // is also useful for comparison with PHG4TpcClusterizer result when running single track events.
249  // The only information written to the cell other than neffelectrons is zbin and pad number, so get those from geometry
250  double zcenter = LayerGeom->get_zcenter(zbin_num);
251  double phicenter = LayerGeom->get_phicenter(pad_num);
252  phi_integral += phicenter * neffelectrons;
253  z_integral += zcenter * neffelectrons;
254  weight += neffelectrons;
255  if (Verbosity() > 1000 && layernum == print_layer)
256  std::cout << " zbin_num " << zbin_num << " zcenter " << zcenter << " pad_num " << pad_num << " phicenter " << phicenter
257  << " neffelectrons " << neffelectrons << " neffelectrons_threshold " << neffelectrons_threshold << std::endl;
258 
259  // Add edep from this electron for this zbin and phi bin combination to the appropriate cell
260  PHG4CellDefs::keytype key = PHG4CellDefs::SizeBinning::genkey(layernum, zbin_num, pad_num);
261  PHG4Cell *cell = g4cells->findCell(key);
262  if (!cell)
263  {
264  cell = new PHG4Cellv1(key);
265  g4cells->AddCell(cell);
266  }
267  cell->add_edep(neffelectrons);
268  cell->add_edep(hiter->first, neffelectrons); // associates g4hit with this edep
269  if (Verbosity() > 100 && layernum == 50) cell->identify();
270  } // end of loop over adc Z bins
271  } // end of loop over zigzag pads
272 
273  /*
274  // Capture the input values at the gem stack and the quick clustering results, elecron-by-electron
275  if (Verbosity() > 0)
276  {
277  assert(ntpad);
278  ntpad->Fill(layernum, phi, phi_integral / weight, z_gem, z_integral / weight);
279  }
280  */
281 
282  if (Verbosity() > 100)
283  if (layernum == print_layer)
284  {
285  std::cout << " hit " << hit << " quick centroid for this electron " << std::endl;
286  std::cout << " phi centroid = " << phi_integral / weight << " phi in " << phi << " phi diff " << phi_integral / weight - phi << std::endl;
287  std::cout << " z centroid = " << z_integral / weight << " z in " << z_gem << " z diff " << z_integral / weight - z_gem << std::endl;
288  // For a single track event, this captures the distribution of single electron centroids on the pad plane for layer print_layer.
289  // The centroid of that should match the cluster centroid found by PHG4TpcClusterizer for layer print_layer, if everything is working
290  // - matches to < .01 cm for a few cases that I checked
291  /*
292  assert(nthit);
293  nthit->Fill(hit, layernum, phi, phi_integral / weight, z_gem, z_integral / weight, weight);
294  */
295  }
296 
297  hit++;
298 
299  return;
300 }
301 
303 {
304  // Jin H.: For the GEM gain in sPHENIX TPC,
305  // Bob pointed out the PHENIX HBD measured it as the Polya function with theta parameter = 0.8.
306  // Just talked with Tom too, he suggest us to start the TPC modeling with simpler exponential function
307  // with lambda parameter of 1/2000, (i.e. Polya function with theta parameter = 0, q_bar = 2000). Please note, this gain variation need to be sampled for each initial electron individually.
308  // Summing over ~30 initial electrons, the distribution is pushed towards more Gauss like.
309  // Bob A.: I like Tom's suggestion to use the exponential distribution as a first approximation
310  // for the single electron gain distribution -
311  // and yes, the parameter you're looking for is of course the slope, which is the inverse gain.
312  double nelec = gsl_ran_exponential(RandomGenerator, averageGEMGain);
313 
314  return nelec;
315 }
316 
317 void PHG4TpcPadPlaneReadout::MapToPadPlane(TrkrHitSetContainer *single_hitsetcontainer, TrkrHitSetContainer *hitsetcontainer, TrkrHitTruthAssoc * /*hittruthassoc*/, const double x_gem, const double y_gem, const double z_gem, PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/)
318 {
319  // One electron per call of this method
320  // The x_gem and y_gem values have already been randomized within the transverse drift diffusion width
321  // The z_gem value already reflects the drift time of the primary electron from the production point, and is randomized within the longitudinal diffusion witdth
322 
323  double phi = atan2(y_gem, x_gem);
324  if (phi > +M_PI) phi -= 2 * M_PI;
325  if (phi < -M_PI) phi += 2 * M_PI;
326 
327  rad_gem = sqrt(x_gem * x_gem + y_gem * y_gem);
328  //std::cout << "Enter new MapToPadPlane with rad_gem " << rad_gem << std::endl;
329 
330  unsigned int layernum = 0;
331 
332  // Find which readout layer this electron ends up in
333 
335  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter = layerrange.first;
336  layeriter != layerrange.second;
337  ++layeriter)
338  {
339  double rad_low = layeriter->second->get_radius() - layeriter->second->get_thickness() / 2.0;
340  double rad_high = layeriter->second->get_radius() + layeriter->second->get_thickness() / 2.0;
341 
342  if (rad_gem > rad_low && rad_gem < rad_high)
343  {
344  // capture the layer where this electron hits sthe gem stack
345  LayerGeom = layeriter->second;
346  layernum = LayerGeom->get_layer();
347  if (Verbosity() > 1000)
348  std::cout << " g4hit id " << hiter->first << " rad_gem " << rad_gem << " rad_low " << rad_low << " rad_high " << rad_high
349  << " layer " << hiter->second->get_layer() << " want to change to " << layernum << std::endl;
350  hiter->second->set_layer(layernum); // have to set here, since the stepping action knows nothing about layers
351  }
352  }
353 
354  if (layernum == 0)
355  {
356  return;
357  }
358 
359  // store phi bins and zbins upfront to avoid repetitive checks on the phi methods
360  const auto phibins = LayerGeom->get_phibins();
361  const auto zbins = LayerGeom->get_zbins();
362 
363  // Create the distribution function of charge on the pad plane around the electron position
364 
365  // The resolution due to pad readout includes the charge spread during GEM multiplication.
366  // this now defaults to 400 microns during construction from Tom (see 8/11 email).
367  // Use the setSigmaT(const double) method to update...
368  // We use a double gaussian to represent the smearing due to the SAMPA chip shaping time - default values of fShapingLead and fShapingTail are for 80 ns SAMPA
369 
370  // amplify the single electron in the gem stack
371  //===============================
372 
373  double nelec = getSingleEGEMAmplification();
374 
375  // Distribute the charge between the pads in phi
376  //====================================
377 
378  if (Verbosity() > 200)
379  std::cout << " populate phi bins for "
380  << " layernum " << layernum
381  << " phi " << phi
382  << " sigmaT " << sigmaT
383  << " zigzag_pads " << zigzag_pads
384  << std::endl;
385 
386  pad_phibin.clear();
387  pad_phibin_share.clear();
388  if (zigzag_pads)
390  else
392 
393  // Normalize the shares so they add up to 1
394  double norm1 = 0.0;
395  for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad)
396  {
397  double pad_share = pad_phibin_share[ipad];
398  norm1 += pad_share;
399  }
400  for (unsigned int iphi = 0; iphi < pad_phibin.size(); ++iphi)
401  pad_phibin_share[iphi] /= norm1;
402 
403  // Distribute the charge between the pads in z
404  //====================================
405  if (Verbosity() > 100 && layernum == print_layer)
406  std::cout << " populate z bins for layernum " << layernum
407  << " with z_gem " << z_gem << " sigmaL[0] " << sigmaL[0] << " sigmaL[1] " << sigmaL[1] << std::endl;
408 
409  adc_zbin.clear();
410  adc_zbin_share.clear();
412 
413  // Normalize the shares so that they add up to 1
414  double znorm = 0.0;
415  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
416  {
417  double bin_share = adc_zbin_share[iz];
418  znorm += bin_share;
419  }
420  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
421  adc_zbin_share[iz] /= znorm;
422 
423  // Fill HitSetContainer
424  //===============
425  // These are used to do a quick clustering for checking
426  double phi_integral = 0.0;
427  double z_integral = 0.0;
428  double weight = 0.0;
429 
430  for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad)
431  {
432  int pad_num = pad_phibin[ipad];
433  double pad_share = pad_phibin_share[ipad];
434 
435  for (unsigned int iz = 0; iz < adc_zbin.size(); ++iz)
436  {
437  int zbin_num = adc_zbin[iz];
438  double adc_bin_share = adc_zbin_share[iz];
439 
440  // Divide electrons from avalanche between bins
441  float neffelectrons = nelec * (pad_share) * (adc_bin_share);
442  if (neffelectrons < neffelectrons_threshold) continue; // skip signals that will be below the noise suppression threshold
443 
444  if (zbin_num >= zbins) std::cout << " Error making key: adc_zbin " << zbin_num << " nzbins " << zbins << std::endl;
445  if (pad_num >= phibins) std::cout << " Error making key: pad_phibin " << pad_num << " nphibins " << phibins << std::endl;
446 
447  // collect information to do simple clustering. Checks operation of PHG4CylinderCellTpcReco, and
448  // is also useful for comparison with PHG4TpcClusterizer result when running single track events.
449  // The only information written to the cell other than neffelectrons is zbin and pad number, so get those from geometry
450  double zcenter = LayerGeom->get_zcenter(zbin_num);
451  double phicenter = LayerGeom->get_phicenter(pad_num);
452  phi_integral += phicenter * neffelectrons;
453  z_integral += zcenter * neffelectrons;
454  weight += neffelectrons;
455  if (Verbosity() > 1 && layernum == print_layer)
456  std::cout << " zbin_num " << zbin_num << " zcenter " << zcenter << " pad_num " << pad_num << " phicenter " << phicenter
457  << " neffelectrons " << neffelectrons << " neffelectrons_threshold " << neffelectrons_threshold << std::endl;
458 
459  // new containers
460  //============
461  // We add the Tpc TrkrHitsets directly to the node using hitsetcontainer
462  // We need to create the TrkrHitSet if not already made - each TrkrHitSet should correspond to a Tpc readout module
463  // The hitset key includes the layer, sector, side
464 
465  // Get the side - 0 for negative z, 1 for positive z
466  unsigned int side = 0;
467  if (zcenter > 0) side = 1;
468  // get the Tpc readout sector - there are 12 sectors with how many pads each?
469  unsigned int pads_per_sector = phibins / 12;
470  unsigned int sector = pad_num / pads_per_sector;
471  TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(layernum, sector, side);
472  // Use existing hitset or add new one if needed
473  TrkrHitSetContainer::Iterator hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey);
474  TrkrHitSetContainer::Iterator single_hitsetit = single_hitsetcontainer->findOrAddHitSet(hitsetkey);
475 
476  // generate the key for this hit, requires zbin and phibin
477  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey((unsigned int) pad_num, (unsigned int) zbin_num);
478  // See if this hit already exists
479  TrkrHit *hit = nullptr;
480  hit = hitsetit->second->getHit(hitkey);
481  if (!hit)
482  {
483  // create a new one
484  hit = new TrkrHitv2();
485  hitsetit->second->addHitSpecificKey(hitkey, hit);
486  }
487  // Either way, add the energy to it -- adc values will be added at digitization
488  //std::cout << " PadPlaneReadout: adding energy " << neffelectrons << " for layer " << layernum << " pad_num " << pad_num << " zbin_num " << zbin_num << " hitkey " << hitkey << std::endl;
489  hit->addEnergy(neffelectrons);
490 
491  // repeat for the single_hitsetcontainer
492  // See if this hit already exists
493  TrkrHit *single_hit = nullptr;
494  single_hit = single_hitsetit->second->getHit(hitkey);
495  if (!single_hit)
496  {
497  // create a new one
498  single_hit = new TrkrHitv2();
499  single_hitsetit->second->addHitSpecificKey(hitkey, single_hit);
500  }
501  // Either way, add the energy to it -- adc values will be added at digitization
502  single_hit->addEnergy(neffelectrons);
503 
504  /*
505  if (Verbosity() > 0)
506  {
507  assert(nthit);
508  nthit->Fill(layernum, pad_num, zbin_num, neffelectrons);
509  }
510  */
511 
512  } // end of loop over adc Z bins
513  } // end of loop over zigzag pads
514 
515  /*
516  // Capture the input values at the gem stack and the quick clustering results, elecron-by-electron
517  if (Verbosity() > 0)
518  {
519  assert(ntpad);
520  ntpad->Fill(layernum, phi, phi_integral / weight, z_gem, z_integral / weight);
521  }
522  */
523 
524  if (Verbosity() > 100)
525  if (layernum == print_layer)
526  {
527  std::cout << " hit " << hit << " quick centroid for this electron " << std::endl;
528  std::cout << " phi centroid = " << phi_integral / weight << " phi in " << phi << " phi diff " << phi_integral / weight - phi << std::endl;
529  std::cout << " z centroid = " << z_integral / weight << " z in " << z_gem << " z diff " << z_integral / weight - z_gem << std::endl;
530  // For a single track event, this captures the distribution of single electron centroids on the pad plane for layer print_layer.
531  // The centroid of that should match the cluster centroid found by PHG4TpcClusterizer for layer print_layer, if everything is working
532  // - matches to < .01 cm for a few cases that I checked
533 
534  /*
535  assert(nthit);
536  nthit->Fill(hit, layernum, phi, phi_integral / weight, z_gem, z_integral / weight, weight);
537  */
538  }
539 
540  hit++;
541 
542  return;
543 }
544 
545 void PHG4TpcPadPlaneReadout::populate_rectangular_phibins(const unsigned int /*layernum*/, const double phi, const double cloud_sig_rp, std::vector<int> &pad_phibin, std::vector<double> &pad_phibin_share)
546 {
547  double cloud_sig_rp_inv = 1. / cloud_sig_rp;
548 
549  const int phibin = LayerGeom->get_phibin(phi);
550  const int nphibins = LayerGeom->get_phibins();
551 
552  double radius = LayerGeom->get_radius();
553  double phidisp = phi - LayerGeom->get_phicenter(phibin);
554  double phistepsize = LayerGeom->get_phistep();
555 
556  // bin the charge in phi - consider phi bins up and down 3 sigma in r-phi
557  int n_rp = int(3 * cloud_sig_rp / (radius * phistepsize) + 1);
558  for (int iphi = -n_rp; iphi != n_rp + 1; ++iphi)
559  {
560  int cur_phi_bin = phibin + iphi;
561  // correcting for continuity in phi
562  if (cur_phi_bin < 0)
563  cur_phi_bin += nphibins;
564  else if (cur_phi_bin >= nphibins)
565  cur_phi_bin -= nphibins;
566  if ((cur_phi_bin < 0) || (cur_phi_bin >= nphibins))
567  {
568  std::cout << "PHG4CylinderCellTpcReco => error in phi continuity. Skipping" << std::endl;
569  continue;
570  }
571  // Get the integral of the charge probability distribution in phi inside the current phi step
572  double phiLim1 = 0.5 * M_SQRT2 * ((iphi + 0.5) * phistepsize * radius - phidisp * radius) * cloud_sig_rp_inv;
573  double phiLim2 = 0.5 * M_SQRT2 * ((iphi - 0.5) * phistepsize * radius - phidisp * radius) * cloud_sig_rp_inv;
574  double phi_integral = 0.5 * (erf(phiLim1) - erf(phiLim2));
575 
576  pad_phibin.push_back(cur_phi_bin);
577  pad_phibin_share.push_back(phi_integral);
578  }
579 
580  return;
581 }
582 
583 void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int layernum, const double phi, const double cloud_sig_rp, std::vector<int> &pad_phibin, std::vector<double> &pad_phibin_share)
584 {
585  const double radius = LayerGeom->get_radius();
586  const double phistepsize = LayerGeom->get_phistep();
587  const auto phibins = LayerGeom->get_phibins();
588 
589  // make the charge distribution gaussian
590  double rphi = phi * radius;
591  if (Verbosity() > 100)
592  if (LayerGeom->get_layer() == print_layer)
593  {
594  std::cout << " populate_zigzag_phibins for layer " << layernum << " with radius " << radius << " phi " << phi
595  << " rphi " << rphi << " phistepsize " << phistepsize << std::endl;
596  std::cout << " fcharge created: radius " << radius << " rphi " << rphi << " cloud_sig_rp " << cloud_sig_rp << std::endl;
597  }
598 
599  // Get the range of phi values that completely contains all pads that touch the charge distribution - (nsigmas + 1/2 pad width) in each direction
600  const double philim_low = phi - (_nsigmas * cloud_sig_rp / radius) - phistepsize;
601  const double philim_high = phi + (_nsigmas * cloud_sig_rp / radius) + phistepsize;
602 
603  // Find the pad range that covers this phi range
604  int phibin_low = LayerGeom->get_phibin(philim_low);
605  int phibin_high = LayerGeom->get_phibin(philim_high);
606  int npads = phibin_high - phibin_low;
607 
608  if (Verbosity() > 1000)
609  if (layernum == print_layer)
610  std::cout << " zigzags: phi " << phi << " philim_low " << philim_low << " phibin_low " << phibin_low
611  << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads << std::endl;
612 
613  if (npads < 0 || npads > 9) npads = 9; // can happen if phibin_high wraps around. If so, limit to 10 pads and fix below
614 
615  // Calculate the maximum extent in r-phi of pads in this layer. Pads are assumed to touch the center of the next phi bin on both sides.
616  const double pad_rphi = 2.0 * LayerGeom->get_phistep() * radius;
617 
618  // Make a TF1 for each pad in the phi range
619  using PadParameterSet = std::array<double, 2>;
620  std::array<PadParameterSet, 10> pad_parameters;
621  std::array<int, 10> pad_keep;
622  for (int ipad = 0; ipad <= npads; ipad++)
623  {
624  int pad_now = phibin_low + ipad;
625  // check that we do not exceed the maximum number of pads, wrap if necessary
626  if (pad_now >= phibins) pad_now -= phibins;
627 
628  pad_keep[ipad] = pad_now;
629  const double rphi_pad_now = LayerGeom->get_phicenter(pad_now) * radius;
630  pad_parameters[ipad] = {{pad_rphi / 2.0, rphi_pad_now}};
631 
632  if (Verbosity() > 1000)
633  if (layernum == print_layer)
634  std::cout << " zigzags: make fpad for ipad " << ipad << " pad_now " << pad_now << " pad_rphi/2 " << pad_rphi / 2.0
635  << " rphi_pad_now " << rphi_pad_now << std::endl;
636  }
637 
638  // Now make a loop that steps through the charge distribution and evaluates the response at that point on each pad
639  std::array<double, 10> overlap = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
640 
641  // use analytic integral
642  for (int ipad = 0; ipad <= npads; ipad++)
643  {
644  const double pitch = pad_parameters[ipad][0];
645  const double x_loc = pad_parameters[ipad][1] - rphi;
646  const double sigma = cloud_sig_rp;
647 
648  // calculate fraction of the total charge on this strip
649  /*
650  this corresponds to integrating the charge distribution Gaussian function (centered on rphi and of width cloud_sig_rp),
651  convoluted with a strip response function, which is triangular from -pitch to +pitch, with a maximum of 1. at stript center
652  */
653  overlap[ipad] =
654  (pitch - x_loc) * (std::erf(x_loc / (M_SQRT2 * sigma)) - std::erf((x_loc - pitch) / (M_SQRT2 * sigma))) / (pitch * 2) + (pitch + x_loc) * (std::erf((x_loc + pitch) / (M_SQRT2 * sigma)) - std::erf(x_loc / (M_SQRT2 * sigma))) / (pitch * 2) + (gaus(x_loc - pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch + (gaus(x_loc + pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch;
655  }
656 
657  // now we have the overlap for each pad
658  for (int ipad = 0; ipad <= npads; ipad++)
659  {
660  pad_phibin.push_back(pad_keep[ipad]);
661  pad_phibin_share.push_back(overlap[ipad]);
662  if (rad_gem < output_radius) std::cout << " zigzags: for pad " << ipad << " integral is " << overlap[ipad] << std::endl;
663  }
664 
665  return;
666 }
667 
668 void PHG4TpcPadPlaneReadout::populate_zbins(const double z, const std::array<double, 2> &cloud_sig_zz, std::vector<int> &adc_zbin, std::vector<double> &adc_zbin_share)
669 {
670  int zbin = LayerGeom->get_zbin(z);
671  if (zbin < 0 || zbin > LayerGeom->get_zbins())
672  {
673  //std::cout << " z bin is outside range, return" << std::endl;
674  return;
675  }
676 
677  double zstepsize = LayerGeom->get_zstep();
678  double zdisp = z - LayerGeom->get_zcenter(zbin);
679 
680  if (Verbosity() > 1000)
681  std::cout << " input: z " << z << " zbin " << zbin << " zstepsize " << zstepsize << " z center " << LayerGeom->get_zcenter(zbin) << " zdisp " << zdisp << std::endl;
682 
683  // Because of diffusion, hits can be shared across the membrane, so we allow all z bins
684  int min_cell_zbin = 0;
685  int max_cell_zbin = NZBins - 1;
686 
687  double cloud_sig_zz_inv[2];
688  cloud_sig_zz_inv[0] = 1. / cloud_sig_zz[0];
689  cloud_sig_zz_inv[1] = 1. / cloud_sig_zz[1];
690 
691  int zsect = 0;
692  if (z < 0)
693  zsect = -1;
694  else
695  zsect = 1;
696 
697  int n_zz = int(3 * (cloud_sig_zz[0] + cloud_sig_zz[1]) / (2.0 * zstepsize) + 1);
698  if (Verbosity() > 1000) std::cout << " n_zz " << n_zz << " cloud_sigzz[0] " << cloud_sig_zz[0] << " cloud_sig_zz[1] " << cloud_sig_zz[1] << std::endl;
699  for (int iz = -n_zz; iz != n_zz + 1; ++iz)
700  {
701  int cur_z_bin = zbin + iz;
702  if ((cur_z_bin < min_cell_zbin) || (cur_z_bin > max_cell_zbin)) continue;
703 
704  if (Verbosity() > 1000)
705  std::cout << " iz " << iz << " cur_z_bin " << cur_z_bin << " min_cell_zbin " << min_cell_zbin << " max_cell_zbin " << max_cell_zbin << std::endl;
706 
707  double z_integral = 0.0;
708  if (iz == 0)
709  {
710  // the crossover between lead and tail shaping occurs in this bin
711  int index1 = -1;
712  int index2 = -1;
713  if (zsect == -1)
714  {
715  index1 = 0;
716  index2 = 1;
717  }
718  else
719  {
720  index1 = 1;
721  index2 = 0;
722  }
723 
724  double zLim1 = 0.0;
725  double zLim2 = 0.5 * M_SQRT2 * (-0.5 * zstepsize - zdisp) * cloud_sig_zz_inv[index1];
726  // 1/2 * the erf is the integral probability from the argument Z value to zero, so this is the integral probability between the Z limits
727  double z_integral1 = 0.5 * (erf(zLim1) - erf(zLim2));
728 
729  if (Verbosity() > 1000)
730  if (LayerGeom->get_layer() == print_layer)
731  std::cout << " populate_zbins: cur_z_bin " << cur_z_bin << " center z " << LayerGeom->get_zcenter(cur_z_bin)
732  << " index1 " << index1 << " zLim1 " << zLim1 << " zLim2 " << zLim2 << " z_integral1 " << z_integral1 << std::endl;
733 
734  zLim2 = 0.0;
735  zLim1 = 0.5 * M_SQRT2 * (0.5 * zstepsize - zdisp) * cloud_sig_zz_inv[index2];
736  double z_integral2 = 0.5 * (erf(zLim1) - erf(zLim2));
737 
738  if (Verbosity() > 1000)
739  if (LayerGeom->get_layer() == print_layer)
740  std::cout << " populate_zbins: cur_z_bin " << cur_z_bin << " center z " << LayerGeom->get_zcenter(cur_z_bin)
741  << " index2 " << index2 << " zLim1 " << zLim1 << " zLim2 " << zLim2 << " z_integral2 " << z_integral2 << std::endl;
742 
743  z_integral = z_integral1 + z_integral2;
744  }
745  else
746  {
747  // The non zero bins are entirely in the lead or tail region
748  // lead or tail depends on which side of the membrane
749  int index = 0;
750  if (iz < 0)
751  {
752  if (zsect == -1)
753  index = 0;
754  else
755  index = 1;
756  }
757  else
758  {
759  if (zsect == -1)
760  index = 1;
761  else
762  index = 0;
763  }
764  double zLim1 = 0.5 * M_SQRT2 * ((iz + 0.5) * zstepsize - zdisp) * cloud_sig_zz_inv[index];
765  double zLim2 = 0.5 * M_SQRT2 * ((iz - 0.5) * zstepsize - zdisp) * cloud_sig_zz_inv[index];
766  z_integral = 0.5 * (erf(zLim1) - erf(zLim2));
767 
768  if (Verbosity() > 1000)
769  if (LayerGeom->get_layer() == print_layer)
770  std::cout << " populate_zbins: z_bin " << cur_z_bin << " center z " << LayerGeom->get_zcenter(cur_z_bin)
771  << " index " << index << " zLim1 " << zLim1 << " zLim2 " << zLim2 << " z_integral " << z_integral << std::endl;
772  }
773 
774  adc_zbin.push_back(cur_z_bin);
775  adc_zbin_share.push_back(z_integral);
776  }
777 
778  return;
779 }
780 
782 {
783  set_default_int_param("ntpc_layers_inner", 16);
784  set_default_int_param("ntpc_layers_mid", 16);
785  set_default_int_param("ntpc_layers_outer", 16);
786 
787  set_default_int_param("tpc_minlayer_inner", 7);
788 
789  set_default_double_param("tpc_minradius_inner", 30.0); // cm
790  set_default_double_param("tpc_minradius_mid", 40.0);
791  set_default_double_param("tpc_minradius_outer", 60.0);
792 
793  set_default_double_param("tpc_maxradius_inner", 40.0); // cm
794  set_default_double_param("tpc_maxradius_mid", 60.0);
795  set_default_double_param("tpc_maxradius_outer", 77.0); // from Tom
796 
797  set_default_double_param("neffelectrons_threshold", 1.0);
798  set_default_double_param("maxdriftlength", 105.5); // cm
799  set_default_double_param("drift_velocity", 8.0 / 1000.0); // cm/ns
800  set_default_double_param("tpc_adc_clock", 53.0); // ns, for 18.8 MHz clock
801 
802  set_default_double_param("gem_cloud_sigma", 0.04); // cm = 400 microns
803  set_default_double_param("sampa_shaping_lead", 32.0); // ns, for 80 ns SAMPA
804  set_default_double_param("sampa_shaping_tail", 48.0); // ns, for 80 ns SAMPA
805 
806  set_default_int_param("ntpc_phibins_inner", 1152);
807  set_default_int_param("ntpc_phibins_mid", 1536);
808  set_default_int_param("ntpc_phibins_outer", 2304);
809 
810  set_default_int_param("zigzag_pads", 1);
811 
812  // GEM Gain
813  /*
814  hp (2020/09/04): gain changed from 2000 to 1400, to accomodate gas mixture change
815  from Ne/CF4 90/10 to Ne/CF4 50/50, and keep the average charge per particle per pad constant
816  */
817  set_default_double_param("gem_amplification", 1400);
818  return;
819 }
820 
822 {
823  NTpcLayers[0] = get_int_param("ntpc_layers_inner");
824  NTpcLayers[1] = get_int_param("ntpc_layers_mid");
825  NTpcLayers[2] = get_int_param("ntpc_layers_outer");
826 
827  MinLayer[0] = get_int_param("tpc_minlayer_inner");
828  MinLayer[1] = MinLayer[0] + NTpcLayers[0];
829  MinLayer[2] = MinLayer[1] + NTpcLayers[1];
830 
831  neffelectrons_threshold = get_double_param("neffelectrons_threshold");
832 
833  MinRadius[0] = get_double_param("tpc_minradius_inner");
834  MinRadius[1] = get_double_param("tpc_minradius_mid");
835  MinRadius[2] = get_double_param("tpc_minradius_outer");
836 
837  MaxRadius[0] = get_double_param("tpc_maxradius_inner");
838  MaxRadius[1] = get_double_param("tpc_maxradius_mid");
839  MaxRadius[2] = get_double_param("tpc_maxradius_outer");
840 
841  Thickness[0] = NTpcLayers[0] <= 0 ? 0 : (MaxRadius[0] - MinRadius[0]) / NTpcLayers[0];
842  Thickness[1] = NTpcLayers[1] <= 0 ? 0 : (MaxRadius[1] - MinRadius[1]) / NTpcLayers[1];
843  Thickness[2] = NTpcLayers[2] <= 0 ? 0 : (MaxRadius[2] - MinRadius[2]) / NTpcLayers[2];
844 
845  MaxRadius[0] = get_double_param("tpc_maxradius_inner");
846  MaxRadius[1] = get_double_param("tpc_maxradius_mid");
847  MaxRadius[2] = get_double_param("tpc_maxradius_outer");
848 
849  sigmaT = get_double_param("gem_cloud_sigma");
850  sigmaL[0] = get_double_param("sampa_shaping_lead") * get_double_param("drift_velocity");
851  sigmaL[1] = get_double_param("sampa_shaping_tail") * get_double_param("drift_velocity");
852 
853  tpc_drift_velocity = get_double_param("drift_velocity");
854  tpc_adc_clock = get_double_param("tpc_adc_clock");
856  MaxZ = get_double_param("maxdriftlength");
857  MinZ = -MaxZ;
858  NZBins = (int) ((MaxZ - MinZ) / ZBinWidth) + 1;
859 
860  NPhiBins[0] = get_int_param("ntpc_phibins_inner");
861  NPhiBins[1] = get_int_param("ntpc_phibins_mid");
862  NPhiBins[2] = get_int_param("ntpc_phibins_outer");
863 
864  PhiBinWidth[0] = 2.0 * M_PI / (double) NPhiBins[0];
865  PhiBinWidth[1] = 2.0 * M_PI / (double) NPhiBins[1];
866  PhiBinWidth[2] = 2.0 * M_PI / (double) NPhiBins[2];
867 
868  zigzag_pads = get_int_param("zigzag_pads");
869 
870  averageGEMGain = get_double_param("gem_amplification");
871 }