LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Protected Member Functions | Protected Attributes | Private Attributes | List of all members
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, 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::Shower > Reconstruct (geo::GeometryCore const &geom, 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, 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::CalorimetryAlg * fCaloAlg
 
bool _Ecorrection {true}
 
bool fVerbosity {true}
 
double fcalodEdxlength {1000}
 
double fdEdxlength {2.4}
 
bool fUseArea {true}
 

Detailed Description

Definition at line 20 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  }
showerreco::ShowerCluster_t
Definition: ShowerRecoAlgBase.h:37
showerreco::ShowerRecoAlgBase::fInputClusters
std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
Input clusters.
Definition: ShowerRecoAlgBase.h:78
void showerreco::ShowerRecoAlg::CaloAlgo ( calo::CalorimetryAlg *  alg)
inline

Definition at line 27 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

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

Function to reorganize input cluster information.

Definition at line 67 of file ShowerRecoAlgBase.h.

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

67 {}
std::vector<::recob::Shower > showerreco::ShowerRecoAlgBase::Reconstruct ( geo::GeometryCore const &  geom,
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().

38  {
39  ProcessInputClusters();
40 
41  std::vector<::recob::Shower> output;
42  output.reserve(fInputClusters.size());
43 
44  for (auto const& clusters : fInputClusters)
45  output.push_back(RecoOneShower(geom, clockData, detProp, clusters));
46 
47  return output;
48  }
showerreco::ShowerRecoAlgBase::ProcessInputClusters
virtual void ProcessInputClusters()
Function to reorganize input cluster information.
Definition: ShowerRecoAlgBase.h:67
showerreco::ShowerRecoAlgBase::fInputClusters
std::vector< std::vector< showerreco::ShowerCluster_t > > fInputClusters
Input clusters.
Definition: ShowerRecoAlgBase.h:78
showerreco::ShowerRecoAlgBase::RecoOneShower
virtual ::recob::Shower RecoOneShower(geo::GeometryCore const &geom, 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,
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.

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

Function to set whether to use E correction.

Definition at line 30 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

30 { _Ecorrection = on; }
void showerreco::ShowerRecoAlg::SetUseArea ( bool  on)
inline

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

Definition at line 23 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

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

Definition at line 25 of file ShowerRecoAlg.h.

Referenced by ShowerReco3D::ShowerReco3D().

25 { fVerbosity = verbose; }

Member Data Documentation

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

Definition at line 40 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

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

Definition at line 39 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

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

Definition at line 42 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

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

Definition at line 43 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

std::vector<std::vector<showerreco::ShowerCluster_t> > showerreco::ShowerRecoAlgBase::fInputClusters
protectedinherited

Input clusters.

Definition at line 78 of file ShowerRecoAlgBase.h.

Referenced by showerreco::ShowerRecoAlgBase::AppendInputClusters(), showerreco::ShowerRecoAlgBase::Reconstruct(), and showerreco::ShowerRecoAlgBase::Reset().

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

Definition at line 44 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

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

Definition at line 41 of file ShowerRecoAlg.h.

Referenced by RecoOneShower().

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