25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_rng.h>
39 , ChargeToPeakVolts(20)
40 , ADCSignalConversionGain(std::numeric_limits<float>::signaling_NaN())
41 , ADCNoiseConversionGain(std::numeric_limits<float>::signaling_NaN())
44 std::cout <<
Name() <<
" random seed: " << seed << std::endl;
50 std::cout <<
"Creating PHG4TpcDigitizer with name = " << name << std::endl;
79 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
92 std::cout <<
"====================== PHG4TpcDigitizer::InitRun() =====================" << std::endl;
93 for (std::map<int, unsigned int>::iterator iter =
_max_adc.begin();
97 std::cout <<
" Max ADC in Layer #" << iter->first <<
" = " << iter->second << std::endl;
99 for (std::map<int, float>::iterator iter =
_energy_scale.begin();
103 std::cout <<
" Energy per ADC in Layer #" << iter->first <<
" = " << 1.0e6 * iter->second <<
" keV" << std::endl;
105 std::cout <<
"===========================================================================" << std::endl;
124 if (!geom_container)
return;
128 layeriter != layerrange.second;
131 int layer = layeriter->second->get_layer();
133 float pitch = (layeriter->second)->get_phistep() * (layeriter->second)->get_radius();
134 float length = (layeriter->second)->get_zstep();
136 float minpath = pitch;
137 if (length < minpath) minpath =
length;
138 if (thickness < minpath) minpath =
thickness;
139 float mip_e = 0.003876 * minpath;
153 unsigned int print_layer = 18;
215 findNode::getClass<PHG4CylinderCellGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
218 std::cout <<
PHWHERE <<
"ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
222 TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
223 if (!trkrhitsetcontainer)
225 std::cout <<
"Could not locate TRKR_HITSET node, quit! " << std::endl;
229 TrkrHitTruthAssoc *hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
232 std::cout <<
PHWHERE <<
" Could not locate TRKR_HITTRUTHASSOC node, quit! " << std::endl;
243 hitset_iter != hitset_range.second;
252 if (layer == print_layer)
253 std::cout <<
"new: PHG4TpcDigitizer: processing hits for layer " << layer <<
" hitsetkey " << hitsetkey << std::endl;
265 for (
int iphi = 0; iphi < nphibins; iphi++)
274 hit_iter != hit_range.second;
296 for (
int iz = 0; iz < nzbins; iz++)
297 z_sorted_hits.push_back(std::vector<TrkrHitSet::ConstIterator>());
300 if (layer == print_layer) std::cout << std::endl;
310 if (layer == print_layer)
313 std::cout <<
"iphi " << iphi <<
" adding existing hit to z vector for layer " << layer <<
" zbin " << zbin <<
" hitkey " << hitkey <<
" pad " <<
TpcDefs::getPad(hitkey)
326 for (
int iz = 0; iz < nzbins; iz++)
335 adc_input.push_back(adc_input_voltage + noise_voltage);
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;
354 if (layer == print_layer)
355 std::cout <<
"noise hit: layer " << layer <<
" iphi " << iphi <<
" iz " << iz
357 <<
" adc_input " <<
adc_input[iz] << std::endl;
362 std::cout <<
"Impossible value of is_populated, iz = " << iz <<
" is_populated = " <<
is_populated[iz] << std::endl;
372 for (
int iz = 0; iz < nzbins / 2; iz++)
374 if (iz < binpointer)
continue;
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;
385 for (
int izup = 0; izup < 5; izup++)
387 if (iz + izup < nzbins / 2 && iz + izup >= 0)
389 unsigned int adc_output = (
unsigned int) (
adc_input[iz + izup] * 1024.0 / 2200.0);
390 if (
adc_input[iz + izup] < 0) adc_output = 0;
391 if (adc_output > 1023) adc_output = 1023;
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;
401 hit = hitset_iter->second->getHit(hitkey);
407 hitset_iter->second->addHitSpecificKey(hitkey, hit);
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;
429 hit = hitset_iter->second->getHit(hitkey);
442 binpointer = nzbins - 1;
443 for (
int iz = nzbins - 1; iz >= nzbins / 2; iz--)
445 if (iz > binpointer)
continue;
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;
456 for (
int izup = 0; izup < 5; izup++)
458 if (iz - izup < nzbins && iz - izup >= nzbins / 2)
460 unsigned int adc_output = (
unsigned int) (
adc_input[iz - izup] * 1024.0 / 2200.0);
461 if (
adc_input[iz - izup] < 0) adc_output = 0;
462 if (adc_output > 1023) adc_output = 1023;
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;
472 hit = hitset_iter->second->getHit(hitkey);
478 hitset_iter->second->addHitSpecificKey(hitkey, hit);
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;
499 hit = hitset_iter->second->getHit(hitkey);
517 std::cout <<
"From PHG4TpcDigitizer: hitsetcontainer dump at end before cleaning:" << std::endl;
519 std::vector<std::pair<TrkrDefs::hitsetkey, TrkrDefs::hitkey>> delete_hitkey_list;
524 hitset_iter != hitset_range_now.second;
533 std::cout <<
"PHG4TpcDigitizer: hitset with key: " << hitsetkey <<
" in layer " << layer <<
" with sector " << sector <<
" side " << side << std::endl;
539 hit_iter != hit_range.second;
543 TrkrHit *tpchit = hit_iter->second;
547 <<
" adc " << tpchit->
getAdc() << std::endl;
549 if (tpchit->
getAdc() == 0)
553 std::cout <<
" -- this hit not digitized - delete it" << std::endl;
556 delete_hitkey_list.push_back(std::make_pair(hitsetkey, hitkey));
562 for (
unsigned int i = 0; i < delete_hitkey_list.size(); i++)
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;
572 hittruthassoc->
removeAssoc(delete_hitkey_list[i].first, delete_hitkey_list[i].second);
577 std::cout <<
"From PHG4TpcDigitizer: hitsetcontainer dump at end after cleaning:" << std::endl;
581 hitset_iter != hitset_range_now.second;
587 if (layer != print_layer)
continue;
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;
597 hit_iter != hit_range.second;
601 TrkrHit *tpchit = hit_iter->second;
604 <<
" adc " << tpchit->
getAdc() << std::endl;
606 if (tpchit->
getAdc() == 0)
608 std::cout <<
" Oops! -- this hit not digitized and not deleted!" << std::endl;