26 #include <gsl/gsl_randist.h>
27 #include <gsl/gsl_rng.h>
48 inline T gaus(
const T &
x,
const T &sigma)
50 return std::exp(-
square(x / sigma) / 2) / (sigma * std::sqrt(2 *
M_PI));
73 if (
Verbosity()) std::cout <<
"PHG4TpcPadPlaneReadout: CreateReadoutGeometry: " << std::endl;
75 for (
int iregion = 0; iregion < 3; ++iregion)
81 std::cout <<
" layer " <<
layer <<
" MinLayer " <<
MinLayer[iregion] <<
" region " << iregion
102 std::cout <<
"increase Tpc cellsize, number of cells "
104 <<
" exceed 5.1M limit" << std::endl;
123 double phi = atan2(y_gem, x_gem);
127 rad_gem = sqrt(x_gem * x_gem + y_gem * y_gem);
130 unsigned int layernum = 0;
136 layeriter != layerrange.second;
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;
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);
179 std::cout <<
" populate phi bins for "
180 <<
" layernum " << layernum
195 for (
unsigned int ipad = 0; ipad <
pad_phibin.size(); ++ipad)
200 for (
unsigned int iphi = 0; iphi <
pad_phibin.size(); ++iphi)
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;
215 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
220 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
226 double phi_integral = 0.0;
227 double z_integral = 0.0;
230 for (
unsigned int ipad = 0; ipad <
pad_phibin.size(); ++ipad)
235 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
241 float neffelectrons = nelec * (pad_share) * (adc_bin_share);
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;
252 phi_integral += phicenter * neffelectrons;
253 z_integral += zcenter * neffelectrons;
254 weight += neffelectrons;
256 std::cout <<
" zbin_num " << zbin_num <<
" zcenter " << zcenter <<
" pad_num " << pad_num <<
" phicenter " << phicenter
268 cell->
add_edep(hiter->first, neffelectrons);
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;
323 double phi = atan2(y_gem, x_gem);
327 rad_gem = sqrt(x_gem * x_gem + y_gem * y_gem);
330 unsigned int layernum = 0;
336 layeriter != layerrange.second;
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;
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);
379 std::cout <<
" populate phi bins for "
380 <<
" layernum " << layernum
395 for (
unsigned int ipad = 0; ipad <
pad_phibin.size(); ++ipad)
400 for (
unsigned int iphi = 0; iphi <
pad_phibin.size(); ++iphi)
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;
415 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
420 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
426 double phi_integral = 0.0;
427 double z_integral = 0.0;
430 for (
unsigned int ipad = 0; ipad <
pad_phibin.size(); ++ipad)
435 for (
unsigned int iz = 0; iz <
adc_zbin.size(); ++iz)
441 float neffelectrons = nelec * (pad_share) * (adc_bin_share);
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;
452 phi_integral += phicenter * neffelectrons;
453 z_integral += zcenter * neffelectrons;
454 weight += neffelectrons;
456 std::cout <<
" zbin_num " << zbin_num <<
" zcenter " << zcenter <<
" pad_num " << pad_num <<
" phicenter " << phicenter
466 unsigned int side = 0;
467 if (zcenter > 0) side = 1;
469 unsigned int pads_per_sector =
phibins / 12;
470 unsigned int sector = pad_num / pads_per_sector;
480 hit = hitsetit->second->getHit(hitkey);
485 hitsetit->second->addHitSpecificKey(hitkey, hit);
494 single_hit = single_hitsetit->second->getHit(hitkey);
499 single_hitsetit->second->addHitSpecificKey(hitkey, single_hit);
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;
547 double cloud_sig_rp_inv = 1. / cloud_sig_rp;
557 int n_rp =
int(3 * cloud_sig_rp / (radius * phistepsize) + 1);
558 for (
int iphi = -n_rp; iphi != n_rp + 1; ++iphi)
560 int cur_phi_bin = phibin + iphi;
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))
568 std::cout <<
"PHG4CylinderCellTpcReco => error in phi continuity. Skipping" << std::endl;
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));
576 pad_phibin.push_back(cur_phi_bin);
577 pad_phibin_share.push_back(phi_integral);
590 double rphi = phi *
radius;
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;
600 const double philim_low = phi - (
_nsigmas * cloud_sig_rp /
radius) - phistepsize;
601 const double philim_high = phi + (
_nsigmas * cloud_sig_rp /
radius) + phistepsize;
606 int npads = phibin_high - phibin_low;
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;
613 if (npads < 0 || npads > 9) npads = 9;
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++)
624 int pad_now = phibin_low + ipad;
628 pad_keep[ipad] = pad_now;
630 pad_parameters[ipad] = {{pad_rphi / 2.0, rphi_pad_now}};
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;
639 std::array<double, 10> overlap = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
642 for (
int ipad = 0; ipad <= npads; ipad++)
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;
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;
658 for (
int ipad = 0; ipad <= npads; ipad++)
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;
681 std::cout <<
" input: z " << z <<
" zbin " << zbin <<
" zstepsize " << zstepsize <<
" z center " <<
LayerGeom->
get_zcenter(zbin) <<
" zdisp " << zdisp << std::endl;
684 int min_cell_zbin = 0;
685 int max_cell_zbin =
NZBins - 1;
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];
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)
701 int cur_z_bin = zbin + iz;
702 if ((cur_z_bin < min_cell_zbin) || (cur_z_bin > max_cell_zbin))
continue;
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;
707 double z_integral = 0.0;
725 double zLim2 = 0.5 * M_SQRT2 * (-0.5 * zstepsize - zdisp) * cloud_sig_zz_inv[index1];
727 double z_integral1 = 0.5 * (erf(zLim1) - erf(zLim2));
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;
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));
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;
743 z_integral = z_integral1 + z_integral2;
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));
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;
774 adc_zbin.push_back(cur_z_bin);
775 adc_zbin_share.push_back(z_integral);