LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
showerreco::ShowerRecoAlg Class Referenceabstract

#include "ShowerRecoAlg.h"

Inheritance diagram for showerreco::ShowerRecoAlg:
showerreco::ShowerRecoAlgBase

Public Member Functions

void SetUseArea (bool on)
 Function to decide if to use Area or Pulse Amplitude for calculations. More...
 
void Verbose (bool verbose)
 
void CaloAlgo (calo::CalorimetryAlg *alg)
 
void setEcorrection (bool on)
 Function to set whether to use E correction. More...
 
recob::Shower RecoOneShower (geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< showerreco::ShowerCluster_t > const &)
 Function to reconstruct a shower. More...
 
virtual void Reset ()
 Function to reset algorithm, to be called @ beginning of each event. More...
 
virtual void AppendInputClusters (const std::vector< cluster::ClusterParamsAlg > &cpan_v)
 Setter for a matched combination of clusters. More...
 
std::vector< recob::ShowerReconstruct (geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp)
 Execute reconstruction. More...
 

Protected Member Functions

virtual void ProcessInputClusters ()
 Function to reorganize input cluster information. More...
 
virtual ::recob::Shower RecoOneShower (geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< showerreco::ShowerCluster_t > &clusters)=0
 Function to reconstruct one shower. More...
 

Protected Attributes

std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
 Input clusters. More...
 

Private Attributes

calo::CalorimetryAlgfCaloAlg
 
bool _Ecorrection {true}
 
bool fVerbosity {true}
 
double fcalodEdxlength {1000}
 
double fdEdxlength {2.4}
 
bool fUseArea {true}
 

Detailed Description

Definition at line 17 of file ShowerRecoAlg.h.

Member Function Documentation

void showerreco::ShowerRecoAlgBase::AppendInputClusters ( const std::vector< cluster::ClusterParamsAlg > &  cpan_v)
virtualinherited

Setter for a matched combination of clusters.

Definition at line 14 of file ShowerRecoAlgBase.cxx.

References showerreco::ShowerRecoAlgBase::fInputClusters.

Referenced by showerreco::ShowerRecoManager::Process().

16  {
17  std::vector<::showerreco::ShowerCluster_t> clusters;
18  clusters.reserve(cpan_v.size());
19 
20  for (auto const& cpan : cpan_v) {
21 
22  clusters.push_back(::showerreco::ShowerCluster_t());
23 
24  (*clusters.rbegin()).start_point = cpan.GetParams().start_point;
25  (*clusters.rbegin()).end_point = cpan.GetParams().end_point;
26  (*clusters.rbegin()).angle_2d = cpan.GetParams().angle_2d;
27  (*clusters.rbegin()).plane_id = cpan.Plane();
28  (*clusters.rbegin()).hit_vector = cpan.GetHitVector();
29  }
30 
31  fInputClusters.push_back(clusters);
32  }
std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
Input clusters.
void showerreco::ShowerRecoAlg::CaloAlgo ( calo::CalorimetryAlg alg)
inline

Definition at line 24 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

24 { fCaloAlg = alg; }
calo::CalorimetryAlg * fCaloAlg
Definition: ShowerRecoAlg.h:37
virtual void showerreco::ShowerRecoAlgBase::ProcessInputClusters ( )
inlineprotectedvirtualinherited

Function to reorganize input cluster information.

Definition at line 66 of file ShowerRecoAlgBase.h.

Referenced by showerreco::ShowerRecoAlgBase::Reconstruct().

66 {}
std::vector<::recob::Shower > showerreco::ShowerRecoAlgBase::Reconstruct ( geo::GeometryCore const &  geom,
geo::WireReadoutGeom const &  wireReadoutGeom,
detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp 
)
inherited

Execute reconstruction.

Definition at line 34 of file ShowerRecoAlgBase.cxx.

References showerreco::ShowerRecoAlgBase::fInputClusters, showerreco::ShowerRecoAlgBase::ProcessInputClusters(), and showerreco::ShowerRecoAlgBase::RecoOneShower().

Referenced by showerreco::ShowerRecoManager::Process().

39  {
41 
42  std::vector<::recob::Shower> output;
43  output.reserve(fInputClusters.size());
44 
45  for (auto const& clusters : fInputClusters)
46  output.push_back(RecoOneShower(geom, wireReadoutGeom, clockData, detProp, clusters));
47 
48  return output;
49  }
virtual void ProcessInputClusters()
Function to reorganize input cluster information.
std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
Input clusters.
virtual ::recob::Shower RecoOneShower(geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< showerreco::ShowerCluster_t > &clusters)=0
Function to reconstruct one shower.
recob::Shower showerreco::ShowerRecoAlg::RecoOneShower ( geo::GeometryCore const &  geom,
geo::WireReadoutGeom const &  wireReadoutGeom,
detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< showerreco::ShowerCluster_t > const &  clusters 
)

Function to reconstruct a shower.

third loop to get only points inside of 1RMS of value.

Definition at line 16 of file ShowerRecoAlg.cxx.

References _Ecorrection, tca::dEdx(), calo::CalorimetryAlg::dEdx_AMP(), calo::CalorimetryAlg::dEdx_AREA(), showerreco::energy::DEFAULT_ECorr, larg4::dist(), calo::CalorimetryAlg::ElectronsFromADCArea(), calo::CalorimetryAlg::ElectronsFromADCPeak(), fCaloAlg, fcalodEdxlength, fdEdxlength, fUseArea, fVerbosity, util::kGeVToElectrons, calo::CalorimetryAlg::LifetimeCorrection(), recob::Shower::set_dedx(), recob::Shower::set_direction(), recob::Shower::set_length(), recob::Shower::set_start_point(), recob::Shower::set_total_best_plane(), recob::Shower::set_total_energy(), recob::Shower::set_total_MIPenergy(), util::PxPoint::w, and w.

22  {
23  recob::Shower result;
24  //
25  // Reconstruct and store
26  //
27  std::vector<util::PxPoint> fStartPoint; // for each plane
28  std::vector<util::PxPoint> fEndPoint; // for each plane
29  std::vector<double> fOmega2D; // for each plane
30 
31  util::GeometryUtilities const gser{geom, wireReadoutGeom, clockData, detProp};
32  std::vector<double> fEnergy(gser.Nplanes(), -1); // for each plane
33  std::vector<double> fMIPEnergy(gser.Nplanes(), -1); // for each plane
34  std::vector<double> fdEdx(gser.Nplanes(), -1);
35  std::vector<unsigned char> fPlaneID;
36 
37  // First Get Start Points
38  for (auto const& cl : clusters) {
39  fStartPoint.push_back(cl.start_point); // for each plane
40  fEndPoint.push_back(cl.end_point); // for each plane
41  fOmega2D.push_back(cl.angle_2d);
42  fPlaneID.push_back(cl.plane_id);
43  if (fVerbosity) {
44  std::cout << " planes : " << cl.plane_id << " " << cl.start_point.t << " "
45  << cl.start_point.w << " " << cl.end_point.t << " " << cl.end_point.w
46  << " angle2d " << cl.angle_2d << std::endl;
47  }
48  }
49 
50  //decide best two planes. for now, using length in wires - flatter is better.
51  int index_to_use[2] = {0, 1};
52  int best_plane = -1;
53  double best_length = 0;
54  double good_length = 0;
55  for (size_t ip = 0; ip < fPlaneID.size(); ip++) {
56  double dist = fabs(fEndPoint[ip].w - fStartPoint[ip].w);
57  if (dist > best_length) {
58  good_length = best_length;
59  index_to_use[1] = index_to_use[0];
60 
61  best_length = dist;
62  best_plane = fPlaneID.at(ip);
63  index_to_use[0] = ip;
64  }
65  else if (dist > good_length) {
66  good_length = dist;
67  index_to_use[1] = ip;
68  }
69  }
70 
71  // Second Calculate 3D angle and effective pitch and start point
72  double xphi = 0, xtheta = 0;
73 
74  gser.Get3DaxisN(fPlaneID[index_to_use[0]],
75  fPlaneID[index_to_use[1]],
76  fOmega2D[index_to_use[0]] * TMath::Pi() / 180.,
77  fOmega2D[index_to_use[1]] * TMath::Pi() / 180.,
78  xphi,
79  xtheta);
80 
81  if (fVerbosity) std::cout << " new angles: " << xphi << " " << xtheta << std::endl;
82 
83  double xyz[3];
84  // calculate start point here?
85  gser.GetXYZ(&fStartPoint[index_to_use[0]], &fStartPoint[index_to_use[1]], xyz);
86 
87  if (fVerbosity) std::cout << " XYZ: " << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;
88 
89  // Third calculate dE/dx and total energy for all planes, because why not?
90  for (size_t cl_index = 0; cl_index < fPlaneID.size(); ++cl_index) {
91  int plane = fPlaneID.at(cl_index);
92  double newpitch = gser.PitchInView(plane, xphi, xtheta);
93  if (plane == best_plane) best_length *= newpitch / gser.WireToCm();
94 
95  if (fVerbosity) std::cout << std::endl << " Plane: " << plane << std::endl;
96 
97  double totEnergy = 0;
98  double totLowEnergy = 0;
99  double totHighEnergy = 0;
100  double totMIPEnergy = 0;
101  int direction = -1;
102  double dEdx_av = 0, dedx_final = 0;
103  int npoints_first = 0, npoints_sec = 0;
104 
105  //override direction if phi (XZ angle) is less than 90 degrees
106  // this is shady, and needs to be rethought.
107  if (fabs(fOmega2D.at(cl_index)) < 90) direction = 1;
108 
109  auto const& hitlist = clusters.at(cl_index).hit_vector;
110  std::vector<util::PxHit> local_hitlist;
111  local_hitlist.reserve(hitlist.size());
112 
113  for (const auto& theHit : hitlist) {
114  double dEdx_new = 0;
115  double hitElectrons = 0;
116  if (!fUseArea) {
117  dEdx_new = fCaloAlg->dEdx_AMP(
118  clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
119  hitElectrons = fCaloAlg->ElectronsFromADCPeak(theHit.peak, plane);
120  }
121  else {
122  dEdx_new = fCaloAlg->dEdx_AREA(
123  clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
124  hitElectrons = fCaloAlg->ElectronsFromADCArea(theHit.charge, plane);
125  }
126 
127  hitElectrons *=
128  fCaloAlg->LifetimeCorrection(clockData, detProp, theHit.t / gser.TimeToCm());
129 
130  totEnergy += hitElectrons * 1.e3 / (::util::kGeVToElectrons);
131 
132  int multiplier = 1;
133  if (plane < 2) multiplier = 2;
134  if (!fUseArea) { totMIPEnergy += theHit.peak * 0.0061 * multiplier; }
135  else {
136  totMIPEnergy += theHit.charge * 0.00336 * multiplier;
137  }
138 
139  if (fVerbosity && dEdx_new > 1.9 && dEdx_new < 2.1)
140  std::cout << "dEdx_new " << dEdx_new << " " << dEdx_new / theHit.charge * newpitch << " "
141  << theHit.charge * 0.0033 * multiplier / newpitch << std::endl;
142 
143  util::PxPoint OnlinePoint;
144  // calculate the wire,time coordinates of the hit projection on to the 2D shower axis
145  gser.GetPointOnLine(
146  fOmega2D.at(cl_index), &(fStartPoint.at(cl_index)), &theHit, OnlinePoint);
147 
148  double ortdist = gser.Get2DDistance(&OnlinePoint, &theHit);
149  double linedist = gser.Get2DDistance(&OnlinePoint, &(fStartPoint.at(cl_index)));
150 
151  //calculate the distance from the vertex using the effective pitch metric
152  double wdist = ((theHit.w - fStartPoint.at(cl_index).w) * newpitch) *
153  direction; //wdist is always positive
154 
155  if (fVerbosity && dEdx_new < 3.5)
156  std::cout << " CALORIMETRY:"
157  << " Pitch " << newpitch << " dist: " << wdist << " dE/dx: " << dEdx_new
158  << "MeV/cm "
159  << " average: " << totEnergy << "hit: wire, time " << theHit.w / gser.WireToCm()
160  << " " << theHit.t / gser.TimeToCm() << "total energy" << totEnergy
161  << std::endl;
162 
163  if ((wdist < fcalodEdxlength) && (wdist > 0.2)) {
164 
165  // first pass at average dE/dx
166  if (wdist < fdEdxlength
167  //take no hits before vertex (depending on direction)
168  && ((direction == 1 && theHit.w > fStartPoint.at(cl_index).w) ||
169  (direction == -1 && theHit.w < fStartPoint.at(cl_index).w)) &&
170  ortdist < 3.5 && linedist < fdEdxlength) {
171 
172  dEdx_av += dEdx_new;
173  local_hitlist.push_back(theHit);
174 
175  npoints_first++;
176  if (fVerbosity)
177  std::cout << " CALORIMETRY:"
178  << " Pitch " << newpitch << " dist: " << wdist << " dE/dx: " << dEdx_new
179  << "MeV/cm "
180  << " average: " << dEdx_av / npoints_first << " hit: wire, time "
181  << theHit.w << " " << theHit.t << " line,ort " << linedist << " " << ortdist
182  << " direction " << direction << std::endl;
183  }
184 
185  } //end inside range if statement
186 
187  } // end first loop on hits.
188 
189  double mevav2cm = 0;
190  double fRMS_corr = 0;
191 
192  //calculate average dE/dx
193  if (npoints_first > 0) { mevav2cm = dEdx_av / npoints_first; }
194 
195  //loop only on subset of hits
196  for (auto const& theHit : local_hitlist) {
197  double dEdx = 0;
198  if (fUseArea) {
199  dEdx = fCaloAlg->dEdx_AREA(
200  clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
201  }
202  else //this will hopefully go away, once all of the calibration factors are calculated.
203  {
204  dEdx = fCaloAlg->dEdx_AMP(
205  clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
206  }
207 
208  fRMS_corr += (dEdx - mevav2cm) * (dEdx - mevav2cm);
209  }
210 
211  if (npoints_first > 0) { fRMS_corr = std::sqrt(fRMS_corr / npoints_first); }
212 
214  //loop only on subset of hits
215 
216  for (auto const& theHit : local_hitlist) {
217  double dEdx = 0;
218  if (fUseArea) {
219  dEdx = fCaloAlg->dEdx_AREA(
220  clockData, detProp, theHit.charge / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
221  }
222  else //this will hopefully go away, once all of the calibration factors are calculated.
223  {
224  dEdx = fCaloAlg->dEdx_AMP(
225  clockData, detProp, theHit.peak / newpitch, theHit.t / gser.TimeToCm(), theHit.plane);
226  }
227 
228  if (((dEdx > (mevav2cm - fRMS_corr)) && (dEdx < (mevav2cm + fRMS_corr))) ||
229  (newpitch > 0.3 * fdEdxlength)) {
230  dedx_final += dEdx;
231  npoints_sec++;
232  }
233  }
234 
235  if (npoints_sec > 0) { dedx_final /= npoints_sec; }
236 
237  if (fVerbosity) {
238  std::cout << " total ENERGY, birks: " << totEnergy << " MeV "
239  << " |average: " << dedx_final << std::endl
240  << " Energy: lo:" << totLowEnergy << " hi: " << totHighEnergy
241  << " MIP corrected " << totMIPEnergy << std::endl;
242  }
243  // if Energy correction factor to be used
244  // then get it from ShowerCalo.h
245  if (_Ecorrection) { totEnergy *= showerreco::energy::DEFAULT_ECorr; }
246  fEnergy[plane] = totEnergy; // for each plane
247  fMIPEnergy[plane] = totMIPEnergy;
248  fdEdx[plane] = dedx_final;
249 
250  // break;
251  } // end loop on clusters
252 
253  result.set_total_MIPenergy(fMIPEnergy);
254  result.set_total_best_plane(best_plane);
255  result.set_length(best_length);
256  result.set_total_energy(fEnergy);
257 
258  double dirs[3] = {0};
259  gser.GetDirectionCosines(xphi, xtheta, dirs);
260  TVector3 vdirs(dirs[0], dirs[1], dirs[2]);
261 
262  TVector3 vxyz(xyz[0], xyz[1], xyz[2]);
263 
264  result.set_direction(vdirs);
265  result.set_start_point(vxyz);
266  result.set_dedx(fdEdx);
267 
268  // done
269  return result;
270  }
static const double DEFAULT_ECorr
Definition: ShowerCalo.h:26
void set_total_energy(const std::vector< double > &q)
Definition: Shower.h:128
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)
Definition: Shower.h:134
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)
Definition: Shower.h:140
void set_total_best_plane(const int q)
Definition: Shower.h:132
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2671
void set_total_MIPenergy(const std::vector< double > &q)
Definition: Shower.h:130
double w
Definition: PxUtils.h:9
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
calo::CalorimetryAlg * fCaloAlg
Definition: ShowerRecoAlg.h:37
void set_start_point(const TVector3 &xyz)
Definition: Shower.h:136
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Float_t w
Definition: plot.C:20
void set_dedx(const std::vector< double > &q)
Definition: Shower.h:138
virtual ::recob::Shower showerreco::ShowerRecoAlgBase::RecoOneShower ( geo::GeometryCore const &  geom,
geo::WireReadoutGeom const &  wireReadoutGeom,
detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
const std::vector< showerreco::ShowerCluster_t > &  clusters 
)
protectedpure virtualinherited

Function to reconstruct one shower.

Referenced by showerreco::ShowerRecoAlgBase::Reconstruct().

void showerreco::ShowerRecoAlgBase::Reset ( )
virtualinherited

Function to reset algorithm, to be called @ beginning of each event.

Definition at line 9 of file ShowerRecoAlgBase.cxx.

References showerreco::ShowerRecoAlgBase::fInputClusters.

Referenced by showerreco::ShowerRecoManager::Reset().

10  {
11  fInputClusters.clear();
12  }
std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
Input clusters.
void showerreco::ShowerRecoAlg::setEcorrection ( bool  on)
inline

Function to set whether to use E correction.

Definition at line 27 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

void showerreco::ShowerRecoAlg::SetUseArea ( bool  on)
inline

Function to decide if to use Area or Pulse Amplitude for calculations.

Definition at line 20 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

20 { fUseArea = on; }
void showerreco::ShowerRecoAlg::Verbose ( bool  verbose)
inline

Definition at line 22 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

22 { fVerbosity = verbose; }

Member Data Documentation

bool showerreco::ShowerRecoAlg::_Ecorrection {true}
private

Definition at line 38 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

calo::CalorimetryAlg* showerreco::ShowerRecoAlg::fCaloAlg
private

Definition at line 37 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

double showerreco::ShowerRecoAlg::fcalodEdxlength {1000}
private

Definition at line 40 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

double showerreco::ShowerRecoAlg::fdEdxlength {2.4}
private

Definition at line 41 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

std::vector<std::vector<showerreco::ShowerCluster_t> > showerreco::ShowerRecoAlgBase::fInputClusters
protectedinherited
bool showerreco::ShowerRecoAlg::fUseArea {true}
private

Definition at line 42 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

bool showerreco::ShowerRecoAlg::fVerbosity {true}
private

Definition at line 39 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().


The documentation for this class was generated from the following files: