ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcDigitizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcDigitizer.cc
1 #include "PHG4TpcDigitizer.h"
2 
3 #include <trackbase/TrkrDefs.h>
4 #include <trackbase/TrkrHit.h>
5 #include <trackbase/TrkrHitSet.h>
8 #include <trackbase/TrkrHitv2.h>
9 
10 #include <tpc/TpcDefs.h>
11 
14 
16 #include <fun4all/SubsysReco.h> // for SubsysReco
17 #
18 #include <phool/PHCompositeNode.h>
19 #include <phool/PHNode.h> // for PHNode
20 #include <phool/PHNodeIterator.h>
21 #include <phool/PHRandomSeed.h>
22 #include <phool/getClass.h>
23 #include <phool/phool.h> // for PHWHERE
24 
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_rng.h> // for gsl_rng_alloc
27 
28 #include <cstdlib> // for exit
29 #include <iostream>
30 #include <limits>
31 #include <memory> // for allocator_tra...
32 
34  : SubsysReco(name)
35  , TpcMinLayer(7)
36  , ADCThreshold(2700) // electrons
37  , TpcEnc(670) // electrons
38  , Pedestal(50000) // electrons
39  , ChargeToPeakVolts(20) // mV/fC
40  , ADCSignalConversionGain(std::numeric_limits<float>::signaling_NaN()) // will be assigned in PHG4TpcDigitizer::InitRun
41  , ADCNoiseConversionGain(std::numeric_limits<float>::signaling_NaN()) // will be assigned in PHG4TpcDigitizer::InitRun
42 {
43  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
44  std::cout << Name() << " random seed: " << seed << std::endl;
45  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
46  gsl_rng_set(RandomGenerator, seed);
47 
48  if (Verbosity() > 0)
49  {
50  std::cout << "Creating PHG4TpcDigitizer with name = " << name << std::endl;
51  }
52 }
53 
55 {
56  gsl_rng_free(RandomGenerator);
57 }
58 
60 {
62 
63  // Factor that convertes the charge in each z bin into a voltage in each z bin
64  // ChargeToPeakVolts relates TOTAL charge collected to peak voltage, while the cell maker actually distributes the signal
65  // GEM output charge in Z bins using the shaper time response. For 80 ns shaping, the scaleup factor of 2.4 gets the peak voltage right.
66  ADCSignalConversionGain = ChargeToPeakVolts * 1.60e-04 * 2.4; // 20 (or 30) mV/fC * fC/electron * scaleup factor
67  // The noise is by definition the RMS noise width voltage divided by ChargeToPeakVolts
68  ADCNoiseConversionGain = ChargeToPeakVolts * 1.60e-04; // 20 (or 30) mV/fC * fC/electron
69 
70  //-------------
71  // Add Hit Node
72  //-------------
73  PHNodeIterator iter(topNode);
74 
75  // Looking for the DST node
76  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
77  if (!dstNode)
78  {
79  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
81  }
82  PHNodeIterator iter_dst(dstNode);
83 
85 
86  //----------------
87  // Report Settings
88  //----------------
89 
90  if (Verbosity() > 0)
91  {
92  std::cout << "====================== PHG4TpcDigitizer::InitRun() =====================" << std::endl;
93  for (std::map<int, unsigned int>::iterator iter = _max_adc.begin();
94  iter != _max_adc.end();
95  ++iter)
96  {
97  std::cout << " Max ADC in Layer #" << iter->first << " = " << iter->second << std::endl;
98  }
99  for (std::map<int, float>::iterator iter = _energy_scale.begin();
100  iter != _energy_scale.end();
101  ++iter)
102  {
103  std::cout << " Energy per ADC in Layer #" << iter->first << " = " << 1.0e6 * iter->second << " keV" << std::endl;
104  }
105  std::cout << "===========================================================================" << std::endl;
106  }
107 
109 }
110 
112 {
113  DigitizeCylinderCells(topNode);
114 
116 }
117 
119 {
120  // defaults to 8-bit ADC, short-axis MIP placed at 1/4 dynamic range
121 
122  PHG4CylinderCellGeomContainer *geom_container = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
123 
124  if (!geom_container) return;
125 
126  PHG4CylinderCellGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
127  for (PHG4CylinderCellGeomContainer::ConstIterator layeriter = layerrange.first;
128  layeriter != layerrange.second;
129  ++layeriter)
130  {
131  int layer = layeriter->second->get_layer();
132  float thickness = (layeriter->second)->get_thickness();
133  float pitch = (layeriter->second)->get_phistep() * (layeriter->second)->get_radius();
134  float length = (layeriter->second)->get_zstep();
135 
136  float minpath = pitch;
137  if (length < minpath) minpath = length;
138  if (thickness < minpath) minpath = thickness;
139  float mip_e = 0.003876 * minpath;
140 
141  if (_max_adc.find(layer) == _max_adc.end())
142  {
143  _max_adc[layer] = 255;
144  _energy_scale[layer] = mip_e / 64;
145  }
146  }
147 
148  return;
149 }
150 
152 {
153  unsigned int print_layer = 18; // to print diagnostic output for layer 47
154 
155  // Digitizes the Tpc cells that were created in PHG4CylinderCellTpcReco
156  // These contain as edep the number of electrons out of the GEM stack, distributed between Z bins by shaper response and ADC clock window
157  // - i.e. all of the phi and Z bins in a cluster have edep values that add up to the primary charge in the layer times 2000
158 
159  // NOTES:
160  // Modified by ADF June 2018 to do the following:
161  // Add noise to cells before digitizing
162  // Digitize the first adc time bin to exceed the threshold, and the 4 bins after that
163  // If the adc value is still above the threshold after 5 bins, repeat for the next 5 bins
164 
165  // Electron production:
166  // A MIP produces 32 electrons in 1.25 cm of Ne:CF4 gas
167  // The nominal GEM gain is 2000 => 64,000 electrons per MIP through 1.25 cm gas
168  // Thus a MIP produces a charge value out of the GEM stack of 64000/6.242x10^18 = 10.2 fC
169 
170  // SAMPA:
171  // See https://indico.cern.ch/event/489996/timetable/#all.detailed "SAMPA Chip: the New ASIC for the ALICE Tpc and MCH Upgrades", M Bregant
172  // The SAMPA has a maximum output voltage of 2200 mV (but the pedestal is about 200 mV)
173  // The SAMPA shaper is set to 80 ns peaking time
174  // The ADC Digitizes the SAMPA shaper output into 1024 channels
175  // Conversion gains of 20 mV/fC or 30 mV/fC are possible - 1 fC charge input produces a peak volatge out of the shaper of 20 or 30 mV
176  // At 30 mV/fC, the input signal saturates at 2.2 V / 30 mV/fC = 73 fC (say 67 with pedestal not at zero)
177  // At 20 mV/fC, the input signal saturates at 2.2 V / 20 mV/fC = 110 fC (say 100 fC with pedestal not at zero) - assume 20 mV/fC
178  // The equivalent noise charge RMS at 20 mV/fC was measured (w/o detector capacitance) at 490 electrons
179  // - note: this appears to be just the pedestal RMS voltage spread divided by the conversion gain, so it is a bit of a funny number (see below)
180  // - it is better to think of noise and signal in terms of voltage at the input of the ADC
181  // Bregant's slides say 670 electrons ENC for the full chip with 18 pf detector, as in ALICE - should use that number
182 
183  // Signal:
184  // To normalize the signal in each cell, take the entire charge on the pad and multiply by 20 mV/fC to get the adc input AT THE PEAK of the shaper
185  // The cell contents should thus be multipied by the normalization given by:
186  // V_peak = Q_pad (electrons) * 1.6e-04 fC/electron * 20 mV/fC
187  // From the sims, for 80 ns and 18.8 MHz, if we take the input charge and spread it a across the shaping time (which is how it has always been done, and is
188  // not the best way to think about it, because the ADC does not see charge it sees voltage out of a charge integrating preamp followed by a shaper), to get
189  // the voltage at the ADC input right, then the values of Q_(pad,z) have to be scaled up by 2.4
190  // V_(pad,z) = 2.4 * Q_(pad,z) (electrons) * 1.6e-04 fC/electron * 20 mV/fC = Q_(pad,z) * 7.68e-03 (in mV)
191  // ADU_(pad,z) = V_(pad,z) * (1024 ADU / 2200 mV) = V_(pad,z) * 0.465
192  // Remember that Q_(pad,z) is the GEM output charge
193 
194  // Noise:
195  // The ENC is defined as the RMS spread of the ADC pedestal distribution coming out from the SAMPA divided by the corresponding conversion gain.
196  // The full range of the ADC input is 2.2V (which will become 1024 adc counts, i.e. 1024 ADU's).
197  // If you see the RMS of the pedestal in adc counts as 1 at the gain of 20mV/fC, the ENC would be defined by:
198  // (2200 [mV]) * (1/1024) / (20[mV/fC]) / (1.6*10^-4 [fC]) = 671 [electrons]
199  // The RMS noise voltage would be: V_RMS = ENC (electrons) *1.6e-04 fC/electron * 20 mV/fC = ENC (electrons) * 3.2e-03 (in mV)
200  // The ADC readout would be: ADU = V_RMS * (1024 ADU / 2200 mV) = V_RMS * 0.465
201 
202  // The cells that we need to digitize here contain as the energy "edep", which is the number of electrons out of the GEM stack
203  // distributed over the pads and ADC time bins according to the output time distribution of the SAMPA shaper - not what really happens, see above
204  // We convert to volts at the input to the ADC and add noise generated with the RMS value of the noise voltage at the ADC input
205  // We assume the pedestal is zero, for simplicity, so the noise fluctuates above and below zero
206 
207  // Note that zbin = 0 corresponds to -105.5 cm, zbin 248 corresponds to 0 cm, and zbin 497 corresponds to +105.5 cm
208  // increasing time should be bins (497 -> 249) and (0 -> 248)
209 
210  //----------
211  // Get Nodes
212  //----------
213 
214  PHG4CylinderCellGeomContainer *geom_container =
215  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
216  if (!geom_container)
217  {
218  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
219  }
220 
221  // Get the TrkrHitSetContainer node
222  TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
223  if (!trkrhitsetcontainer)
224  {
225  std::cout << "Could not locate TRKR_HITSET node, quit! " << std::endl;
226  exit(1);
227  }
228 
229  TrkrHitTruthAssoc *hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
230  if (!hittruthassoc)
231  {
232  std::cout << PHWHERE << " Could not locate TRKR_HITTRUTHASSOC node, quit! " << std::endl;
233  exit(1);
234  }
235 
236  //-------------
237  // Digitization
238  //-------------
239 
240  // Loop over all hitsets for the Tpc
241  TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
242  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
243  hitset_iter != hitset_range.second;
244  ++hitset_iter)
245  {
246  // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer
247  // get the hitset key so we can find the layer
248  TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
249  const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
250 
251  if (Verbosity() > 2)
252  if (layer == print_layer)
253  std::cout << "new: PHG4TpcDigitizer: processing hits for layer " << layer << " hitsetkey " << hitsetkey << std::endl;
254 
255  // we need the geometry object for this layer
256  PHG4CylinderCellGeom *layergeom = geom_container->GetLayerCellGeom(layer);
257  if (!layergeom)
258  exit(1);
259 
260  // for this hitset, make a vector of a vector of cells for each phibin
261  phi_sorted_hits.clear();
262 
263  // start with an empty vector of cells for each phibin
264  int nphibins = layergeom->get_phibins();
265  for (int iphi = 0; iphi < nphibins; iphi++)
266  {
267  phi_sorted_hits.push_back(std::vector<TrkrHitSet::ConstIterator>());
268  }
269 
270  // get all of the hits from this hitset
271  TrkrHitSet *hitset = hitset_iter->second;
272  TrkrHitSet::ConstRange hit_range = hitset->getHits();
273  for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
274  hit_iter != hit_range.second;
275  ++hit_iter)
276  {
277  // Fill the vector of hits for each phibin
278  unsigned int phibin = TpcDefs::getPad(hit_iter->first);
279  phi_sorted_hits[phibin].push_back(hit_iter);
280  }
281 
282  // For this hitset we have the hits sorted into vectors for each phi
283  // process these vectors one phi bin at a time
284  for (unsigned int iphi = 0; iphi < phi_sorted_hits.size(); iphi++)
285  {
286  if (phi_sorted_hits[iphi].size() == 0)
287  continue;
288 
289  // Populate a vector of hits ordered by Z for each phibin
290  int nzbins = layergeom->get_zbins();
291  is_populated.clear();
292  is_populated.assign(nzbins, 2); // mark all as noise only for now
293  z_sorted_hits.clear();
294 
295  // add an empty vector for each z bin
296  for (int iz = 0; iz < nzbins; iz++)
297  z_sorted_hits.push_back(std::vector<TrkrHitSet::ConstIterator>());
298 
299  if (Verbosity() > 2)
300  if (layer == print_layer) std::cout << std::endl;
301 
302  // add a hit for each z bin that has one
303  for (unsigned int iz = 0; iz < phi_sorted_hits[iphi].size(); iz++)
304  {
305  int zbin = TpcDefs::getTBin(phi_sorted_hits[iphi][iz]->first);
306  is_populated[zbin] = 1; // this bin is a associated with a hit
307  z_sorted_hits[zbin].push_back(phi_sorted_hits[iphi][iz]);
308 
309  if (Verbosity() > 2)
310  if (layer == print_layer)
311  {
312  TrkrDefs::hitkey hitkey = phi_sorted_hits[iphi][iz]->first;
313  std::cout << "iphi " << iphi << " adding existing hit to z vector for layer " << layer << " zbin " << zbin << " hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey)
314  << " z bin " << TpcDefs::getTBin(hitkey) << " energy " << (phi_sorted_hits[iphi][iz]->second)->getEnergy()
315  << std::endl;
316  }
317  }
318 
319  adc_input.clear();
320  adc_hitid.clear();
321  // Now for this phibin we process all bins ordered by Z into hits with noise
322  //======================================================
323  // For this step we take the edep value and convert it to mV at the ADC input
324  // See comments above for how to do this for signal and noise
325 
326  for (int iz = 0; iz < nzbins; iz++)
327  {
328  if (is_populated[iz] == 1)
329  {
330  // This zbin has a hit, add noise
331  float noise = added_noise(); // in electrons
332  float noise_voltage = (Pedestal + noise) * ADCNoiseConversionGain; // mV - from definition of noise charge and pedestal charge
333  float adc_input_voltage = (z_sorted_hits[iz][0]->second)->getEnergy() * ADCSignalConversionGain; // mV, see comments above
334 
335  adc_input.push_back(adc_input_voltage + noise_voltage);
336  adc_hitid.push_back(z_sorted_hits[iz][0]->first);
337 
338  if (Verbosity() > 2)
339  if (layer == print_layer)
340  std::cout << "existing hit: layer " << layer << " iphi " << iphi << " iz " << iz << " edep " << (z_sorted_hits[iz][0]->second)->getEnergy()
341  << " adc gain " << ADCSignalConversionGain << " adc_input_voltage " << adc_input_voltage << " noise voltage " << noise_voltage
342  << " adc_input " << adc_input[iz] << std::endl;
343  }
344  else if (is_populated[iz] == 2)
345  {
346  // This z bin does not have a filled cell, add noise
347  float noise = added_noise(); // returns "electrons"
348  float noise_voltage = (Pedestal + noise) * ADCNoiseConversionGain; // mV - noise "electrons" x conversion gain
349 
350  adc_input.push_back(noise_voltage); // mV
351  adc_hitid.push_back(0); // there is no hit, just add a placeholder in the vector for now, replace it later
352 
353  if (Verbosity() > 2)
354  if (layer == print_layer)
355  std::cout << "noise hit: layer " << layer << " iphi " << iphi << " iz " << iz
356  << " adc gain " << ADCSignalConversionGain << " noise voltage " << noise_voltage
357  << " adc_input " << adc_input[iz] << std::endl;
358  }
359  else
360  {
361  // Cannot happen
362  std::cout << "Impossible value of is_populated, iz = " << iz << " is_populated = " << is_populated[iz] << std::endl;
363  exit(-1);
364  }
365  }
366 
367  // Now we can digitize the entire stream of z bins for this phi bin
368 
369  // start with negative z, the first to arrive is bin 0
370  //================================
371  int binpointer = 0;
372  for (int iz = 0; iz < nzbins / 2; iz++) // 0-247
373  {
374  if (iz < binpointer) continue;
375 
376  if (adc_input[iz] > ADCThreshold * ADCNoiseConversionGain) // convert threshold in "equivalent electrons" to mV
377  {
378  // digitize this bin and the following 4 bins
379 
380  if (Verbosity() > 2)
381  if (layer == print_layer) std::cout << std::endl
382  << "new: (neg z) Above threshold of " << ADCThreshold * ADCNoiseConversionGain << " for phibin " << iphi
383  << " iz " << iz << " with adc_input " << adc_input[iz] << " digitize this and 4 following bins: " << std::endl;
384 
385  for (int izup = 0; izup < 5; izup++)
386  {
387  if (iz + izup < nzbins / 2 && iz + izup >= 0) // stay within the bin limits for negative z
388  {
389  unsigned int adc_output = (unsigned int) (adc_input[iz + izup] * 1024.0 / 2200.0); // input voltage x 1024 channels over 2200 mV max range
390  if (adc_input[iz + izup] < 0) adc_output = 0;
391  if (adc_output > 1023) adc_output = 1023;
392 
393  if (Verbosity() > 2)
394  if (layer == print_layer) std::cout << "new: iphi " << iphi << " (neg z) iz+izup " << iz + izup << " adc_hitid " << adc_hitid[iz + izup]
395  << " adc_input " << adc_input[iz + izup] << " ADCThreshold " << ADCThreshold * ADCNoiseConversionGain
396  << " adc_output " << adc_output << std::endl;
397 
398  // Get the hitkey
399  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, iz + izup);
400  TrkrHit *hit = nullptr;
401  hit = hitset_iter->second->getHit(hitkey);
402 
403  // noise bins do not have TrkrHits associated with them, have to make one
404  if (!hit)
405  {
406  hit = new TrkrHitv2();
407  hitset_iter->second->addHitSpecificKey(hitkey, hit);
408  //hit->addEnergy(adc_input[iz+izup]);
409 
410  if (Verbosity() > 2)
411  if (layer == print_layer) std::cout << "new: adding noise hit for iphi " << iphi << " zbin " << iz + izup
412  << " created new hit with hitkey " << hitkey
413  << " energy " << adc_input[iz + izup] << " adc " << adc_output << std::endl;
414  }
415 
416  hit->setAdc(adc_output);
417 
418  } // end boundary check
419  binpointer++; // skip this bin in future
420  } // end izup loop
421 
422  } // adc threshold if
423  else
424  {
425  // set adc value to zero if there is a hit
426  // Get the hitkey
428  TrkrHit *hit = nullptr;
429  hit = hitset_iter->second->getHit(hitkey);
430  if (hit)
431  {
432  hit->setAdc(0);
433  }
434  // below threshold, move on
435  binpointer++;
436  } // end adc threshold if/else
437 
438  } // end iz loop for negative z
439 
440  // now positive z, the first to arrive is bin 497
441  //===============================
442  binpointer = nzbins - 1;
443  for (int iz = nzbins - 1; iz >= nzbins / 2; iz--) // 495 - 248
444  {
445  if (iz > binpointer) continue;
446 
447  if (adc_input[iz] > ADCThreshold * ADCNoiseConversionGain) // convert threshold in electrons to mV
448  {
449  // digitize this bin and the following 4 bins
450 
451  if (Verbosity() > 2)
452  if (layer == print_layer) std::cout << std::endl
453  << "new: (pos z) Above threshold of " << ADCThreshold * ADCNoiseConversionGain << " for iphi " << iphi << " iz " << iz
454  << " with adc_input " << adc_input[iz] << " digitize this and 4 following bins: " << std::endl;
455 
456  for (int izup = 0; izup < 5; izup++)
457  {
458  if (iz - izup < nzbins && iz - izup >= nzbins / 2)
459  {
460  unsigned int adc_output = (unsigned int) (adc_input[iz - izup] * 1024.0 / 2200.0); // input voltage x 1024 channels over 2200 mV max range
461  if (adc_input[iz - izup] < 0) adc_output = 0;
462  if (adc_output > 1023) adc_output = 1023;
463 
464  if (Verbosity() > 2)
465  if (layer == print_layer) std::cout << "new: iphi " << iphi << " (pos z) iz-izup " << iz - izup << " adc_hitid " << adc_hitid[iz - izup]
466  << " adc_input " << adc_input[iz - izup] << " ADCThreshold " << ADCThreshold * ADCNoiseConversionGain
467  << " adc_output " << adc_output << std::endl;
468 
469  // Get the hitkey
470  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, iz - izup);
471  TrkrHit *hit = nullptr;
472  hit = hitset_iter->second->getHit(hitkey);
473 
474  // noise bins do not have TrkrHits associated with them, have to make one
475  if (!hit)
476  {
477  hit = new TrkrHitv2();
478  hitset_iter->second->addHitSpecificKey(hitkey, hit);
479  //hit->addEnergy(adc_input[iz-izup]);
480 
481  if (Verbosity() > 2)
482  if (layer == print_layer) std::cout << "new: adding noise hit for iphi " << iphi << " zbin " << iz - izup
483  << " created new hit with hitkey " << hitkey
484  << " energy " << adc_input[iz - izup] << " adc " << adc_output << std::endl;
485  }
486 
487  hit->setAdc(adc_output);
488  } // end boundary check
489  binpointer--;
490  } // end izup loop
491 
492  } // adc threshold if
493  else
494  {
495  // set adc value to zero if there is a hit
496  // Get the hitkey
498  TrkrHit *hit = nullptr;
499  hit = hitset_iter->second->getHit(hitkey);
500  if (hit)
501  {
502  hit->setAdc(0);
503  }
504  // below threshold, move on
505  binpointer--;
506  } // end adc threshold if/else
507 
508  } // end iz loop for positive z
509 
510  } // end phibins loop
511 
512  } // end loop over hitsets
513 
514  //======================================================
515  if (Verbosity() > 2)
516  {
517  std::cout << "From PHG4TpcDigitizer: hitsetcontainer dump at end before cleaning:" << std::endl;
518  }
519  std::vector<std::pair<TrkrDefs::hitsetkey, TrkrDefs::hitkey>> delete_hitkey_list;
520 
521  // Clean up undigitized hits - we want all hitsets for the Tpc
522  TrkrHitSetContainer::ConstRange hitset_range_now = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
523  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range_now.first;
524  hitset_iter != hitset_range_now.second;
525  ++hitset_iter)
526  {
527  // we have an iterator to one TrkrHitSet for the Tpc from the trkrHitSetContainer
528  TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
529  const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
530  const int sector = TpcDefs::getSectorId(hitsetkey);
531  const int side = TpcDefs::getSide(hitsetkey);
532  if (Verbosity() > 2)
533  std::cout << "PHG4TpcDigitizer: hitset with key: " << hitsetkey << " in layer " << layer << " with sector " << sector << " side " << side << std::endl;
534 
535  // get all of the hits from this hitset
536  TrkrHitSet *hitset = hitset_iter->second;
537  TrkrHitSet::ConstRange hit_range = hitset->getHits();
538  for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
539  hit_iter != hit_range.second;
540  ++hit_iter)
541  {
542  TrkrDefs::hitkey hitkey = hit_iter->first;
543  TrkrHit *tpchit = hit_iter->second;
544 
545  if (Verbosity() > 2)
546  std::cout << " layer " << layer << " hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey) << " z bin " << TpcDefs::getTBin(hitkey)
547  << " adc " << tpchit->getAdc() << std::endl;
548 
549  if (tpchit->getAdc() == 0)
550  {
551  if (Verbosity() > 20)
552  {
553  std::cout << " -- this hit not digitized - delete it" << std::endl;
554  }
555  // screws up the iterator to delete it here, store the hitkey for later deletion
556  delete_hitkey_list.push_back(std::make_pair(hitsetkey, hitkey));
557  }
558  }
559  }
560 
561  // delete all undigitized hits
562  for (unsigned int i = 0; i < delete_hitkey_list.size(); i++)
563  {
564  TrkrHitSet *hitset = trkrhitsetcontainer->findHitSet(delete_hitkey_list[i].first);
565  const unsigned int layer = TrkrDefs::getLayer(delete_hitkey_list[i].first);
566  hitset->removeHit(delete_hitkey_list[i].second);
567  if (Verbosity() > 20)
568  if (layer == print_layer)
569  std::cout << "removed hit with hitsetkey " << delete_hitkey_list[i].first << " and hitkey " << delete_hitkey_list[i].second << std::endl;
570 
571  // should also delete all entries with this hitkey from the TrkrHitTruthAssoc map
572  hittruthassoc->removeAssoc(delete_hitkey_list[i].first, delete_hitkey_list[i].second);
573  }
574 
575  // Final hitset dump
576  if (Verbosity() > 2)
577  std::cout << "From PHG4TpcDigitizer: hitsetcontainer dump at end after cleaning:" << std::endl;
578  // We want all hitsets for the Tpc
579  TrkrHitSetContainer::ConstRange hitset_range_final = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
580  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range_final.first;
581  hitset_iter != hitset_range_now.second;
582  ++hitset_iter)
583  {
584  // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer
585  TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
586  const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
587  if (layer != print_layer) continue;
588  const int sector = TpcDefs::getSectorId(hitsetkey);
589  const int side = TpcDefs::getSide(hitsetkey);
590  if (Verbosity() > 2 && layer == print_layer)
591  std::cout << "PHG4TpcDigitizer: hitset with key: " << hitsetkey << " in layer " << layer << " with sector " << sector << " side " << side << std::endl;
592 
593  // get all of the hits from this hitset
594  TrkrHitSet *hitset = hitset_iter->second;
595  TrkrHitSet::ConstRange hit_range = hitset->getHits();
596  for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
597  hit_iter != hit_range.second;
598  ++hit_iter)
599  {
600  TrkrDefs::hitkey hitkey = hit_iter->first;
601  TrkrHit *tpchit = hit_iter->second;
602  if (Verbosity() > 2)
603  std::cout << " LAYER " << layer << " hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey) << " z bin " << TpcDefs::getTBin(hitkey)
604  << " adc " << tpchit->getAdc() << std::endl;
605 
606  if (tpchit->getAdc() == 0)
607  {
608  std::cout << " Oops! -- this hit not digitized and not deleted!" << std::endl;
609  }
610  }
611  }
612 
613  //hittruthassoc->identify();
614 
615  return;
616 }
617 
619 {
620  float noise = gsl_ran_gaussian(RandomGenerator, TpcEnc);
621 
622  return noise;
623 }