20 std::vector<showerreco::ShowerCluster_t>
const& clusters)
26 std::vector<util::PxPoint> fStartPoint;
27 std::vector<util::PxPoint> fEndPoint;
28 std::vector<double> fOmega2D;
31 std::vector<double> fEnergy(gser.Nplanes(), -1);
32 std::vector<double> fMIPEnergy(gser.Nplanes(), -1);
33 std::vector<double> fdEdx(gser.Nplanes(), -1);
34 std::vector<unsigned char> fPlaneID;
37 for (
auto const& cl : clusters) {
38 fStartPoint.push_back(cl.start_point);
39 fEndPoint.push_back(cl.end_point);
40 fOmega2D.push_back(cl.angle_2d);
41 fPlaneID.push_back(cl.plane_id);
43 std::cout <<
" planes : " << cl.plane_id <<
" " << cl.start_point.t <<
" " 44 << cl.start_point.w <<
" " << cl.end_point.t <<
" " << cl.end_point.w
45 <<
" angle2d " << cl.angle_2d << std::endl;
50 int index_to_use[2] = {0, 1};
52 double best_length = 0;
53 double good_length = 0;
54 for (
size_t ip = 0; ip < fPlaneID.size(); ip++) {
55 double dist = fabs(fEndPoint[ip].
w - fStartPoint[ip].
w);
56 if (dist > best_length) {
57 good_length = best_length;
58 index_to_use[1] = index_to_use[0];
61 best_plane = fPlaneID.at(ip);
64 else if (dist > good_length) {
71 double xphi = 0, xtheta = 0;
73 gser.Get3DaxisN(fPlaneID[index_to_use[0]],
74 fPlaneID[index_to_use[1]],
75 fOmega2D[index_to_use[0]] * TMath::Pi() / 180.,
76 fOmega2D[index_to_use[1]] * TMath::Pi() / 180.,
80 if (
fVerbosity) std::cout <<
" new angles: " << xphi <<
" " << xtheta << std::endl;
84 gser.GetXYZ(&fStartPoint[index_to_use[0]], &fStartPoint[index_to_use[1]], xyz);
86 if (
fVerbosity) std::cout <<
" XYZ: " << xyz[0] <<
" " << xyz[1] <<
" " << xyz[2] << std::endl;
89 for (
size_t cl_index = 0; cl_index < fPlaneID.size(); ++cl_index) {
90 int plane = fPlaneID.at(cl_index);
91 double newpitch = gser.PitchInView(plane, xphi, xtheta);
92 if (plane == best_plane) best_length *= newpitch / gser.WireToCm();
94 if (
fVerbosity) std::cout << std::endl <<
" Plane: " << plane << std::endl;
97 double totLowEnergy = 0;
98 double totHighEnergy = 0;
99 double totMIPEnergy = 0;
101 double dEdx_av = 0, dedx_final = 0;
102 int npoints_first = 0, npoints_sec = 0;
106 if (fabs(fOmega2D.at(cl_index)) < 90) direction = 1;
108 auto const& hitlist = clusters.at(cl_index).hit_vector;
109 std::vector<util::PxHit> local_hitlist;
110 local_hitlist.reserve(hitlist.size());
112 for (
const auto& theHit : hitlist) {
114 double hitElectrons = 0;
117 clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
122 clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
132 if (plane < 2) multiplier = 2;
133 if (!
fUseArea) { totMIPEnergy += theHit.peak * 0.0061 * multiplier; }
135 totMIPEnergy += theHit.charge * 0.00336 * multiplier;
138 if (
fVerbosity && dEdx_new > 1.9 && dEdx_new < 2.1)
139 std::cout <<
"dEdx_new " << dEdx_new <<
" " << dEdx_new / theHit.charge * newpitch <<
" " 140 << theHit.charge * 0.0033 * multiplier / newpitch << std::endl;
145 fOmega2D.at(cl_index), &(fStartPoint.at(cl_index)), &theHit, OnlinePoint);
147 double ortdist = gser.Get2DDistance(&OnlinePoint, &theHit);
148 double linedist = gser.Get2DDistance(&OnlinePoint, &(fStartPoint.at(cl_index)));
151 double wdist = ((theHit.w - fStartPoint.at(cl_index).w) * newpitch) *
155 std::cout <<
" CALORIMETRY:" 156 <<
" Pitch " << newpitch <<
" dist: " << wdist <<
" dE/dx: " << dEdx_new
158 <<
" average: " << totEnergy <<
"hit: wire, time " << theHit.
w / gser.WireToCm()
159 <<
" " << theHit.t / gser.TimeToCm() <<
"total energy" << totEnergy
167 && ((direction == 1 && theHit.w > fStartPoint.at(cl_index).w) ||
168 (direction == -1 && theHit.w < fStartPoint.at(cl_index).w)) &&
172 local_hitlist.push_back(theHit);
176 std::cout <<
" CALORIMETRY:" 177 <<
" Pitch " << newpitch <<
" dist: " << wdist <<
" dE/dx: " << dEdx_new
179 <<
" average: " << dEdx_av / npoints_first <<
" hit: wire, time " 180 << theHit.w <<
" " << theHit.t <<
" line,ort " << linedist <<
" " << ortdist
181 <<
" direction " << direction << std::endl;
189 double fRMS_corr = 0;
192 if (npoints_first > 0) { mevav2cm = dEdx_av / npoints_first; }
195 for (
auto const& theHit : local_hitlist) {
199 clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
204 clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
207 fRMS_corr += (dEdx - mevav2cm) * (dEdx - mevav2cm);
210 if (npoints_first > 0) { fRMS_corr = std::sqrt(fRMS_corr / npoints_first); }
215 for (
auto const& theHit : local_hitlist) {
219 clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
224 clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
227 if (((dEdx > (mevav2cm - fRMS_corr)) && (dEdx < (mevav2cm + fRMS_corr))) ||
234 if (npoints_sec > 0) { dedx_final /= npoints_sec; }
237 std::cout <<
" total ENERGY, birks: " << totEnergy <<
" MeV " 238 <<
" |average: " << dedx_final << std::endl
239 <<
" Energy: lo:" << totLowEnergy <<
" hi: " << totHighEnergy
240 <<
" MIP corrected " << totMIPEnergy << std::endl;
245 fEnergy[plane] = totEnergy;
246 fMIPEnergy[plane] = totMIPEnergy;
247 fdEdx[plane] = dedx_final;
257 double dirs[3] = {0};
258 gser.GetDirectionCosines(xphi, xtheta, dirs);
259 TVector3 vdirs(dirs[0], dirs[1], dirs[2]);
261 TVector3 vxyz(xyz[0], xyz[1], xyz[2]);
static const double DEFAULT_ECorr
void set_total_energy(const std::vector< double > &q)
double ElectronsFromADCArea(double area, unsigned short plane) const
double ElectronsFromADCPeak(double adc, unsigned short plane) const
constexpr double kGeVToElectrons
23.6eV per ion pair, 1e9 eV/GeV
void set_direction(const TVector3 &dir)
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
void set_length(const double &l)
recob::Shower RecoOneShower(geo::GeometryCore const &geom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< showerreco::ShowerCluster_t > const &)
Function to reconstruct a shower.
Description of geometry of one entire detector.
void set_total_best_plane(const int q)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
void set_total_MIPenergy(const std::vector< double > &q)
Contains all timing reference information for the detector.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Class def header for a class ShowerRecoAlgBase.
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Class def header for a class ShowerCalo.
calo::CalorimetryAlg * fCaloAlg
void set_start_point(const TVector3 &xyz)
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Collection of Physical constants used in LArSoft.
void set_dedx(const std::vector< double > &q)