12 #include <phparameter/PHParameterContainerInterface.h>
13 #include <phparameter/PHParametersContainer.h>
43 , chkenergyconservation(0)
44 , sum_energy_before_cuts(0.)
45 , sum_energy_g4hit(0.)
67 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
74 cout <<
Name() <<
"DST Node missing, doing nothing." << endl;
91 cout <<
"Could not locate g4 hit node " <<
hitnodename << endl;
116 cout <<
"Could not locate geometry node " <<
geonodename << endl;
136 map<int, PHG4CylinderGeom *>::const_iterator miter;
137 pair<map<int, PHG4CylinderGeom *>::const_iterator, map<int, PHG4CylinderGeom *>::const_iterator> begin_end = geo->
get_begin_end();
138 map<int, std::pair<double, double> >::iterator sizeiter;
139 for (miter = begin_end.first; miter != begin_end.second; ++miter)
145 cout <<
Name() <<
": No parameters for detid/layer " << layer
146 <<
", hits from this detid/layer will not be accumulated into cells" << endl;
157 cout <<
"no cell sizes for layer " << layer << endl;
172 double etastepsize = (sizeiter->second).first;
175 if (etastepsize > etamax - etamin)
184 double fract = modf((etamax - etamin) / etastepsize, &d_etabins);
190 etastepsize = (etamax - etamin) / d_etabins;
191 (sizeiter->second).first = etastepsize;
192 int etabins = d_etabins;
193 double etalow = etamin;
194 double etahi = etalow + etastepsize;
195 for (
int i = 0; i < etabins; i++)
197 if (etahi > (etamax + 1.
e-6))
199 cout <<
"etahi: " << etahi <<
", etamax: " << etamax << endl;
201 etahi += etastepsize;
205 double phimax =
M_PI;
206 double phistepsize = (sizeiter->second).
second;
208 if (phistepsize >= phimax-phimin)
217 double fract = modf((phimax - phimin) / phistepsize, &d_phibins);
223 phistepsize = (phimax -
phimin) / d_phibins;
224 (sizeiter->second).
second = phistepsize;
227 double phihi = philow + phistepsize;
228 for (
int i = 0; i <
phibins; i++)
232 cout <<
"phihi: " << phihi <<
", phimax: " << phimax << endl;
234 phihi += phistepsize;
236 pair<int, int> phi_z_bin = make_pair(phibins, etabins);
252 double size_r = (sizeiter->second).first;
255 if (size_r >= circumference)
264 double fract = modf(circumference / size_r, &bins_r);
271 size_r = circumference / bins_r;
272 (sizeiter->second).first = size_r;
273 double phistepsize = 2 *
M_PI / bins_r;
275 double phimax = phimin + phistepsize;
277 for (
int i = 0; i <
nbins[0]; i++)
279 if (phimax > (
M_PI + 1
e-9))
281 cout <<
"phimax: " << phimax <<
", M_PI: " <<
M_PI
282 <<
"phimax-M_PI: " << phimax -
M_PI << endl;
284 phimax += phistepsize;
287 if (size_z >= length_in_z)
296 double fract = modf(length_in_z / size_z, &bins_r);
303 pair<int, int> phi_z_bin = make_pair(nbins[0], nbins[1]);
306 size_z = length_in_z / bins_r;
307 (sizeiter->second).
second = size_z;
308 double zlow = layergeom->
get_zmin();
309 double zhigh = zlow +
size_z;
311 for (
int i = 0; i < nbins[1]; i++)
313 if (zhigh > (layergeom->
get_zmax() + 1
e-9))
315 cout <<
"zhigh: " << zhigh <<
", zmax "
317 <<
", zhigh-zmax: " << zhigh - layergeom->
get_zmax()
340 cout <<
"===================== PHG4CylinderCellReco::InitRun() =====================" << endl;
341 cout <<
" " <<
outdetector <<
" Segmentation Description: " << endl;
342 for (std::map<int, int>::iterator iter =
binning.begin();
345 int layer = iter->first;
351 cout <<
" Layer #" <<
binning.begin()->first <<
"-" <<
binning.rbegin()->first << endl;
358 cout <<
" Layer #" << layer << endl;
363 cout <<
"===========================================================================" << endl;
375 cout <<
"Could not locate g4 hit node " <<
hitnodename << endl;
381 cout <<
"could not locate cell node " <<
cellnodename << endl;
392 map<int, std::pair<double, double> >::iterator sizeiter;
394 pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->
getLayers();
401 for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
417 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
421 if (hiter->second->get_t(0) >
tmin_max[*
layer].second)
continue;
422 if (hiter->second->get_t(1) <
tmin_max[*
layer].first)
continue;
424 pair<double, double> etaphi[2];
427 for (
int i = 0; i < 2; i++)
429 etaphi[i] =
PHG4Utils::get_etaphi(hiter->second->get_x(i), hiter->second->get_y(i), hiter->second->get_z(i));
434 if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
438 if (etabin[0] < 0 || etabin[0] >= nzbins || etabin[1] < 0 || etabin[1] >= nzbins)
447 hiter->second->identify();
452 int intphibin = phibin[0];
453 int intetabin = etabin[0];
454 int intphibinout = phibin[1];
455 int intetabinout = etabin[1];
460 double ay = (etaphi[0]).first;
461 double bx = (etaphi[1]).
second;
462 double by = (etaphi[1]).first;
463 if (intphibin > intphibinout)
466 intphibin = intphibinout;
469 if (intetabin > intetabinout)
472 intetabin = intetabinout;
476 double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
489 vector<double> vdedx;
491 if (intphibin == intphibinout && intetabin == intetabinout)
493 if (
Verbosity() > 0) cout <<
"SINGLE CELL FIRED: " << intphibin <<
" " << intetabin << endl;
494 vphi.push_back(intphibin);
495 veta.push_back(intetabin);
496 vdedx.push_back(trklen);
502 for (
int ibp = intphibin; ibp <= intphibinout; ibp++)
506 for (
int ibz = intetabin; ibz <= intetabinout; ibz++)
516 if (
Verbosity() > 0) cout <<
"CELL FIRED: " << ibp <<
" " << ibz <<
" " << intersect.second << endl;
519 vdedx.push_back(intersect.second);
524 if (
Verbosity() > 0) cout <<
"NUMBER OF FIRED CELLS = " << vphi.size() << endl;
527 for (
unsigned int ii = 0; ii < vphi.size(); ii++)
530 vdedx[ii] = vdedx[ii] / trklen;
531 if (
Verbosity() > 0) cout <<
" CELL " << ii <<
" dE/dX = " << vdedx[ii] << endl;
533 if (
Verbosity() > 0) cout <<
" TOTAL TRACK LENGTH = " << tmpsum <<
" " << trklen << endl;
535 for (
unsigned int i1 = 0; i1 < vphi.size(); i1++)
537 int iphibin = vphi[i1];
538 int ietabin = veta[i1];
543 unsigned long long tmp = iphibin;
544 unsigned long long key = tmp << 32;
548 cout <<
" iphibin " << iphibin <<
" ietabin " << ietabin <<
" key 0x" << hex << key << dec << endl;
562 if (!
isfinite(hiter->second->get_edep() * vdedx[i1]))
564 cout <<
"hit 0x" << hex << hiter->first << dec <<
" not finite, edep: "
565 << hiter->second->
get_edep() <<
" weight " << vdedx[i1] << endl;
567 cell->
add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
568 cell->
add_edep(hiter->second->get_edep() * vdedx[i1]);
573 cell->
add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
575 if (!
isfinite(hiter->second->get_edep() * vdedx[i1]))
577 cout <<
PHWHERE <<
" invalid energy dep " << hiter->second->get_edep()
578 <<
" or path length: " << vdedx[i1] << endl;
598 <<
", energy dep: " <<
it->second->get_edep()
605 cout <<
Name() <<
": found " << numcells <<
" eta/phi cells with energy deposition" << endl;
614 cout <<
"logical screwup!!! no sizes for layer " << *layer << endl;
617 double zstepsize = (sizeiter->second).
second;
620 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
624 if (hiter->second->get_t(0) >
tmin_max[*
layer].second)
continue;
625 if (hiter->second->get_t(1) <
tmin_max[*
layer].first)
continue;
635 if (
Verbosity() > 0) cout <<
"--------- new hit in layer # " << *layer << endl;
637 for (
int i = 0; i < 2; i++)
639 xinout[i] = hiter->second->get_x(i);
640 yinout[i] = hiter->second->get_y(i);
641 px[i] = hiter->second->get_px(i);
642 py[i] = hiter->second->get_py(i);
643 phi[i] = atan2(hiter->second->get_y(i), hiter->second->get_x(i));
644 z[i] = hiter->second->get_z(i);
646 zbin[i] = geo->
get_zbin(hiter->second->get_z(i));
648 if (
Verbosity() > 0) cout <<
" " << i <<
" phibin: " << phibin[i] <<
", phi: " << phi[i] <<
", stepsize: " << phistepsize << endl;
649 if (
Verbosity() > 0) cout <<
" " << i <<
" zbin: " << zbin[i] <<
", z = " << hiter->second->get_z(i) <<
", stepsize: " << zstepsize <<
" offset: " <<
zmin_max[*
layer].first << endl;
652 if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
656 if (zbin[0] < 0 || zbin[0] >= nzbins || zbin[1] < 0 || zbin[1] >= nzbins)
663 hiter->second->identify();
668 int intphibin = phibin[0];
669 int intzbin = zbin[0];
670 int intphibinout = phibin[1];
671 int intzbinout = zbin[1];
675 cout <<
" phi bin range: " << intphibin <<
" to " << intphibinout <<
" phi: " << phi[0] <<
" to " << phi[1] << endl;
676 cout <<
" Z bin range: " << intzbin <<
" to " << intzbinout <<
" Z: " << z[0] <<
" to " << z[1] << endl;
677 cout <<
" phi difference: " << (phi[1] - phi[0]) * 1000. <<
" milliradians." << endl;
678 cout <<
" phi difference: " << 2.5 * (phi[1] - phi[0]) * 10000. <<
" microns." << endl;
679 cout <<
" path length = " << sqrt((xinout[1] - xinout[0]) * (xinout[1] - xinout[0]) + (yinout[1] - yinout[0]) * (yinout[1] - yinout[0])) << endl;
680 cout <<
" px = " << px[0] <<
" " << px[1] << endl;
681 cout <<
" py = " << py[0] <<
" " << py[1] << endl;
682 cout <<
" x = " << xinout[0] <<
" " << xinout[1] << endl;
683 cout <<
" y = " << yinout[0] <<
" " << yinout[1] << endl;
692 if (intphibin > intphibinout)
695 intphibin = intphibinout;
698 if (intzbin > intzbinout)
701 intzbin = intzbinout;
705 double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
718 vector<double> vdedx;
720 if (intphibin == intphibinout && intzbin == intzbinout)
724 cout <<
"SINGLE CELL FIRED: " << intphibin <<
" " << intzbin << endl;
726 vphi.push_back(intphibin);
727 vz.push_back(intzbin);
728 vdedx.push_back(trklen);
733 double zstep_half = geo->
get_zstep() / 2.;
734 for (
int ibp = intphibin; ibp <= intphibinout; ibp++)
738 for (
int ibz = intzbin; ibz <= intzbinout; ibz++)
747 if (
Verbosity() > 0) cout <<
"CELL FIRED: " << ibp <<
" " << ibz <<
" " << intersect.second << endl;
750 vdedx.push_back(intersect.second);
755 if (
Verbosity() > 0) cout <<
"NUMBER OF FIRED CELLS = " << vz.size() << endl;
758 for (
unsigned int ii = 0; ii < vz.size(); ii++)
761 vdedx[ii] = vdedx[ii] / trklen;
762 if (
Verbosity() > 0) cout <<
" CELL " << ii <<
" dE/dX = " << vdedx[ii] << endl;
764 if (
Verbosity() > 0) cout <<
" TOTAL TRACK LENGTH = " << tmpsum <<
" " << trklen << endl;
766 for (
unsigned int i1 = 0; i1 < vphi.size(); i1++)
768 int iphibin = vphi[i1];
771 unsigned long long tmp = iphibin;
772 unsigned long long key = tmp << 32;
776 cout <<
" iphibin " << iphibin <<
" izbin " << izbin <<
" key 0x" << hex << key << dec << endl;
787 cout <<
" add energy to existing cell for key = " <<
cellptmap.find(key)->first << endl;
792 cout <<
" NAN lighy yield with vdedx[i1] = " << vdedx[i1]
793 <<
" and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << endl;
800 cout <<
" did not find a previous entry for key = 0x" << hex << key << dec <<
" create a new one" << endl;
806 if (!
isfinite(hiter->second->get_edep() * vdedx[i1]))
808 cout <<
"hit 0x" << hex << hiter->first << dec <<
" not finite, edep: "
809 << hiter->second->
get_edep() <<
" weight " << vdedx[i1] << endl;
811 cell->
add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
812 cell->
add_edep(hiter->second->get_edep() * vdedx[i1]);
818 cout <<
" NAN lighy yield with vdedx[i1] = " << vdedx[i1]
819 <<
" and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << endl;
822 cell->
add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
841 <<
", energy dep: " <<
it->second->get_edep()
848 cout <<
"found " << numcells <<
" z/phi cells with energy deposition" << endl;
856 cout <<
"cellptmap for layer " << *layer <<
" has final length " <<
cellptmap.size();
865 cout <<
" reset it to " <<
cellptmap.size() << endl;
880 cout <<
"size for layer " << detid <<
" already set" << endl;
892 cout <<
"size for layer " << detid <<
" already set" << endl;
903 cell_size[i] = std::make_pair(sizeA, sizeB);
917 double sum_energy_cells = 0.;
918 double sum_energy_stored_hits = 0.;
919 double sum_energy_stored_showers = 0.;
922 for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
924 sum_energy_cells += citer->second->get_edep();
928 sum_energy_stored_hits += iter->second;
933 sum_energy_stored_showers += iter->second;
937 if (sum_energy_stored_hits > 0)
939 if (fabs(sum_energy_cells - sum_energy_stored_hits) / sum_energy_cells > 1
e-6)
941 cout <<
"energy mismatch between cell energy " << sum_energy_cells
942 <<
" and stored hit energies " << sum_energy_stored_hits
946 if (sum_energy_stored_showers > 0)
948 if (fabs(sum_energy_cells - sum_energy_stored_showers) / sum_energy_cells > 1
e-6)
950 cout <<
"energy mismatch between cell energy " << sum_energy_cells
951 <<
" and stored shower energies " << sum_energy_stored_showers
957 cout <<
"energy mismatch between cells: " << sum_energy_cells
963 cout <<
Name() <<
": sum cell energy: " << sum_energy_cells <<
" GeV" << endl;
964 cout <<
Name() <<
": sum shower energy: " << sum_energy_stored_showers <<
" GeV" << endl;
965 cout <<
Name() <<
": sum stored hit energy: " << sum_energy_stored_hits <<
" GeV" << endl;
974 cout <<
Name() <<
": sum cell energy: " << sum_energy_cells <<
" GeV" << endl;
975 cout <<
Name() <<
": sum shower energy: " << sum_energy_stored_showers <<
" GeV" << endl;
976 cout <<
Name() <<
": sum stored hit energy: " << sum_energy_stored_hits <<
" GeV" << endl;