LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
ShowerRecoAlg.cxx
Go to the documentation of this file.
1 #include "ShowerRecoAlg.h"
3 
10 
11 #include "TMath.h"
12 #include "TVector3.h"
13 
14 namespace showerreco {
15 
17  {
18 
20 
21  fcalodEdxlength=1000;
22  fdEdxlength=2.4;
23  fUseArea = true;
24  fVerbosity = true;
25  _Ecorrection = true;
26  }
27 
28 
29  ::recob::Shower ShowerRecoAlg::RecoOneShower(const std::vector< ::showerreco::ShowerCluster_t>& clusters)
30  {
31 
32  ::recob::Shower result;
33  //
34  // Reconstruct and store
35  //
36  std::vector < util::PxPoint > fStartPoint; // for each plane
37  std::vector < util::PxPoint > fEndPoint; // for each plane
38  std::vector < double > fOmega2D; // for each plane
39 
40  std::vector < double > fEnergy (fGSer->Nplanes(),-1); // for each plane
41  std::vector < double > fMIPEnergy(fGSer->Nplanes(),-1); // for each plane
42  std::vector < double > fdEdx (fGSer->Nplanes(),-1);
43  std::vector <unsigned char> fPlaneID;
44 
45  // First Get Start Points
46  for(auto const & cl : clusters)
47  {
48  fStartPoint.push_back(cl.start_point); // for each plane
49  fEndPoint.push_back(cl.end_point); // for each plane
50  fOmega2D.push_back(cl.angle_2d);
51  fPlaneID.push_back(cl.plane_id);
52  if(fVerbosity) {
53  std::cout << " planes : " << cl.plane_id
54  << " " << cl.start_point.t
55  << " " << cl.start_point.w
56  << " " << cl.end_point.t
57  << " " << cl.end_point.w
58  <<" angle2d "<< cl.angle_2d
59  << std::endl;
60  }
61  }
62 
63 
64  //decide best two planes. for now, using length in wires - flatter is better.
65  int index_to_use[2]={0,1};
66  int best_plane=-1;
67  //int good_plane=-1;
68  double best_length=0;
69  double good_length=0;
70  for (size_t ip=0;ip<fPlaneID.size();ip++)
71  {
72  double dist = fabs( fEndPoint[ip].w - fStartPoint[ip].w );
73  if(dist > best_length )
74  {
75  good_length = best_length;
76  //good_plane = best_plane;
77  index_to_use[1] = index_to_use[0];
78 
79  best_length = dist;
80  best_plane = fPlaneID.at(ip);
81  index_to_use[0] = ip;
82  }
83  else if(dist > good_length)
84  {
85  good_length = dist;
86  //good_plane = fPlaneID.at(ip);
87  index_to_use[1] = ip;
88  }
89  }
90 
91  // Second Calculate 3D angle and effective pitch and start point
92  double xphi=0,xtheta=0;
93 
94 
95  fGSer->Get3DaxisN(fPlaneID[index_to_use[0]],
96  fPlaneID[index_to_use[1]],
97  fOmega2D[index_to_use[0]]*TMath::Pi()/180.,
98  fOmega2D[index_to_use[1]]*TMath::Pi()/180.,
99  xphi,
100  xtheta);
101 
102  if(fVerbosity)
103  std::cout << " new angles: " << xphi << " " << xtheta << std::endl;
104 
105 
106 
107  double xyz[3];
108  // calculate start point here?
109  fGSer-> GetXYZ(&(fStartPoint[index_to_use[0]]),&(fStartPoint[index_to_use[1]]),xyz);
110 
111  if(fVerbosity)
112  std::cout << " XYZ: " << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;
113 
114  // Third calculate dE/dx and total energy for all planes, because why not?
115  //for(auto const &clustit : clusters)
116  for(size_t cl_index=0; cl_index < fPlaneID.size(); ++cl_index)
117  {
118  int plane = fPlaneID.at(cl_index);
119  double newpitch=fGSer->PitchInView(plane,xphi,xtheta);
120  if(plane == best_plane) best_length *= newpitch / fGSer->WireToCm();
121 
122  if(fVerbosity)
123  std::cout << std::endl << " Plane: " << plane << std::endl;
124 
125  double totEnergy=0;
126  double totLowEnergy=0;
127  double totHighEnergy=0;
128  double totMIPEnergy=0;
129  int direction=-1;
130  //double RMS_dedx=0;
131  double dEdx_av=0,dedx_final=0;
132  int npoints_first=0, npoints_sec=0;
133 
134  //override direction if phi (XZ angle) is less than 90 degrees
135  // this is shady, and needs to be rethought.
136  if(fabs(fOmega2D.at(cl_index)) < 90)
137  direction=1;
138 
139  auto const& hitlist = clusters.at(cl_index).hit_vector;
140  std::vector<util::PxHit> local_hitlist;
141  local_hitlist.reserve(hitlist.size());
142 
143  for(const auto& theHit : hitlist){
144 
145  double dEdx_new=0;
146  double hitElectrons = 0;
147  //double Bcorr_half;
148  //double dEdx_sub;
149  // double dEdx_MIP;
150 
151  // dEdx_new = fCaloAlg->dEdx_AREA(theHit, newpitch );
152  //Bcorr_half = 2.*fCaloAlg->dEdx_AREA(theHit->Charge()/2.,theHit->PeakTime(), newpitch, plane); ;
153  //dEdx_sub = fCaloAlg->dEdx_AREA(theHit->Charge()-PION_CORR,theHit->PeakTime(), newpitch, plane); ;
154  // dEdx_MIP = fCaloAlg->dEdx_AREA_forceMIP(theHit, newpitch );
155  if(!fUseArea)
156  {
157  dEdx_new = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
158  hitElectrons = fCaloAlg->ElectronsFromADCPeak(theHit.peak, plane);
159  }
160  else
161  {
162  dEdx_new = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
163  hitElectrons = fCaloAlg->ElectronsFromADCArea(theHit.charge, plane);
164  }
165 
166  hitElectrons *= fCaloAlg->LifetimeCorrection(theHit.t / fGSer->TimeToCm());
167 
168  totEnergy += hitElectrons * 1.e3 / (::util::kGeVToElectrons);
169 
170  //totNewCnrg+=dEdx_MIP;
171  // if(dEdx_new < 3.5 && dEdx_new >0 )
172  // {
173  // totLowEnergy +=dEdx_new*newpitch;
174  // totMIPEnergy += dEdx_new*newpitch;
175  // }
176  // else
177  // {
178  // totHighEnergy += dEdx_new*newpitch;
179  int multiplier=1;
180  if(plane < 2) multiplier =2 ;
181  if(!fUseArea){
182  totMIPEnergy += theHit.peak * 0.0061*multiplier;
183  }
184  else{
185  totMIPEnergy += theHit.charge * 0.00336*multiplier;
186  }
187  // }
188 
189  if(fVerbosity && dEdx_new >1.9 && dEdx_new <2.1)
190  std::cout << "dEdx_new " << dEdx_new << " " <<dEdx_new/theHit.charge*newpitch << " "<< theHit.charge *0.0033*multiplier/newpitch << std::endl;
191 
192  util::PxPoint OnlinePoint;
193  // calculate the wire,time coordinates of the hit projection on to the 2D shower axis
194  fGSer->GetPointOnLine(fOmega2D.at(cl_index),
195  &(fStartPoint.at(cl_index)),
196  &theHit,
197  OnlinePoint);
198 
199  double ortdist=fGSer->Get2DDistance(&OnlinePoint,&theHit);
200  double linedist=fGSer->Get2DDistance(&OnlinePoint,&(fStartPoint.at(cl_index)));
201 
202 
203  //calculate the distance from the vertex using the effective pitch metric
204  double wdist=((theHit.w-fStartPoint.at(cl_index).w)*newpitch)*direction; //wdist is always positive
205 
206  if(fVerbosity && dEdx_new < 3.5 )
207  std::cout << " CALORIMETRY:"
208  << " Pitch " <<newpitch
209  << " dist: " << wdist
210  << " dE/dx: " << dEdx_new << "MeV/cm "
211  << " average: " << totEnergy
212  << "hit: wire, time " << theHit.w/fGSer->WireToCm() << " " << theHit.t/fGSer->TimeToCm()
213  << "total energy" << totEnergy
214  << std::endl;
215 
216  // if( (fabs(wdist)<fcalodEdxlength)&&(fabs(wdist)>0.2)){
217  if( (wdist<fcalodEdxlength)&&(wdist>0.2)){
218 
219  //vdEdx.push_back(dEdx_new);
220  //vresRange.push_back(fabs(wdist));
221  //vdQdx.push_back(theHit->Charge(true)/newpitch);
222  //Trk_Length=wdist;
223 
224  //fTrkPitchC=fNPitch[set][plane];
225  //Kin_En+=dEdx_new*newpitch;
226  //npoints_calo++;
227  //sum+=dEdx_new;
228 
229 
230 
231  //fDistribChargeADC[set].push_back(ch_adc); //vector with the first De/Dx points
232  // fDistribChargeMeV[set].push_back(dEdx_new); //vector with the first De/Dx points
233  // fDistribHalfChargeMeV[set].push_back(Bcorr_half);
234  // fDistribChargeposition[set].push_back(wdist); //vector with the first De/Dx points' positions
235  // fDistribChargeMeVMIPsub[set].push_back(dEdx_sub);
236  //
237 
238 
239  //first pass at average dE/dx
240  if(wdist<fdEdxlength
241  //take no hits before vertex (depending on direction)
242  && ((direction==1 && theHit.w > fStartPoint.at(cl_index).w) || (direction==-1 && theHit.w < fStartPoint.at(cl_index).w) )
243  && ortdist<3.5
244  && linedist < fdEdxlength ){
245 
246  dEdx_av+= dEdx_new;
247  local_hitlist.push_back(theHit);
248 
249  npoints_first++;
250  if(fVerbosity)
251  std::cout << " CALORIMETRY:" << " Pitch " <<newpitch
252  << " dist: " << wdist
253  << " dE/dx: " << dEdx_new << "MeV/cm "
254  << " average: " << dEdx_av/npoints_first
255  << " hit: wire, time " << theHit.w << " " << theHit.t
256  << " line,ort " << linedist << " " << ortdist
257  << " direction " << direction
258  << std::endl;
259  }
260 
261  }//end inside range if statement
262 
263  }// end first loop on hits.
264 
265 
266 
267  double mevav2cm=0;
268  double fRMS_corr=0;
269 
270  //calculate average dE/dx
271  if(npoints_first>0) {
272  mevav2cm=dEdx_av/npoints_first;
273 
274  }
275 
276  //loop only on subset of hits
277  for(auto const & theHit : local_hitlist){
278  double dEdx=0;
279  if(fUseArea)
280  {
281  dEdx = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
282  }
283  else //this will hopefully go away, once all of the calibration factors are calculated.
284  {
285  dEdx = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
286  }
287  //fDistribAfterMin[set].push_back(MinBefore);
288  //fDistribBeforeMin[set].push_back(MinAfter);
289  // auto position = hitIter - local_hitlist.begin() ;
290 
291  fRMS_corr+=(dEdx-mevav2cm)*(dEdx-mevav2cm);
292 
293  }
294 
295 
296 
297  if(npoints_first>0)
298  {
299  fRMS_corr=TMath::Sqrt(fRMS_corr/npoints_first);
300  }
301 
303  //loop only on subset of hits
304 
305  for(auto const & theHit : local_hitlist){
306  double dEdx=0;
307  if(fUseArea)
308  {
309  dEdx = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
310  }
311  else //this will hopefully go away, once all of the calibration factors are calculated.
312  {
313  dEdx = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
314  }
315  //fDistribAfterMin[set].push_back(MinBefore);
316  //fDistribBeforeMin[set].push_back(MinAfter);
317 
318  if( ((dEdx > (mevav2cm-fRMS_corr) ) && (dEdx < (mevav2cm+fRMS_corr) ))
319  ||
320  (newpitch > 0.3*fdEdxlength )
321  ) {
322  dedx_final+= dEdx;
323  npoints_sec++;
324  }
325  }
326 
327  if(npoints_sec>0){
328  dedx_final/=npoints_sec;
329  }
330 
331  if(fVerbosity){
333  std::cout << " total ENERGY, birks: " << totEnergy << " MeV "
334  << " |average: " << dedx_final << std::endl
335  << " Energy: lo:" << totLowEnergy << " hi: " << totHighEnergy
336  << " MIP corrected " << totMIPEnergy << std::endl;
337  }
338  // if Energy correction factor to be used
339  // then get it from ShowerCalo.h
340  if (_Ecorrection) { totEnergy *= showerreco::energy::DEFAULT_ECorr; }
341  fEnergy[plane] = totEnergy; // for each plane
342  fMIPEnergy[plane] = totMIPEnergy;
343  fdEdx[plane] = dedx_final;
344 
345  // break;
346  } // end loop on clusters
347 
348  result.set_total_MIPenergy(fMIPEnergy);
349  result.set_total_best_plane(best_plane);
350  result.set_length(best_length);
351  result.set_total_energy(fEnergy);
352  //result.set_total_energy_err (std::vector< Double_t > q) { fSigmaTotalEnergy = q; }
353 
354  double dirs[3]={0};
355  fGSer->GetDirectionCosines(xphi,xtheta,dirs);
356  TVector3 vdirs(dirs[0],dirs[1],dirs[2]);
357 
358  TVector3 vxyz(xyz[0], xyz[1], xyz[2]);
359 
360  result.set_direction(vdirs);
361  //result.set_direction_err (TVector3 dir_e) { fSigmaDCosStart = dir_e; }
362  result.set_start_point(vxyz);
363  //result.set_start_point_err (TVector3 xyz_e) { fSigmaXYZstart = xyz_e; }
364  result.set_dedx (fdEdx);
365  //result.set_dedx_err (std::vector< Double_t > q) { fSigmadEdx = q; }
366 
367  // done
368  return result;
369  }
370 
371 }
virtual ::recob::Shower RecoOneShower(const std::vector< ::showerreco::ShowerCluster_t > &)
Function to reconstruct a shower.
Class def header for a class ShowerRecoAlg.
static const GeometryUtilities * GetME()
Double_t Get2DDistance(Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const
double dEdx_AREA(art::Ptr< recob::Hit > hit, double pitch, double T0=0) const
::calo::CalorimetryAlg * fCaloAlg
Calorimetry algorithm.
bool _Ecorrection
Boolean -> decide if to use energy correction or not.
Definition: ShowerRecoAlg.h:64
static const double DEFAULT_ECorr
Definition: ShowerCalo.h:28
void set_total_energy(const std::vector< double > &q)
Definition: Shower.h:129
Double_t TimeToCm() const
Double_t WireToCm() const
util::GeometryUtilities * fGSer
Definition: ShowerRecoAlg.h:59
Double_t PitchInView(UInt_t plane, Double_t phi, Double_t theta) const
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:135
bool fVerbosity
Verbosity flag.
void set_length(const double &l)
Definition: Shower.h:141
float dEdx(TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2191
void set_total_best_plane(const int q)
Definition: Shower.h:133
void set_total_MIPenergy(const std::vector< double > &q)
Definition: Shower.h:131
Int_t GetPointOnLine(Double_t slope, Double_t intercept, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const
Class def header for a class ShowerRecoAlgBase.
void GetDirectionCosines(Double_t phi, Double_t theta, Double_t *dirs) const
Class def header for a class ShowerCalo.
void set_start_point(const TVector3 &xyz)
Definition: Shower.h:137
double LifetimeCorrection(double time, double T0=0) const
Float_t w
Definition: plot.C:23
ShowerRecoAlg()
Default constructor.
Int_t Get3DaxisN(Int_t iplane0, Int_t iplane1, Double_t omega0, Double_t omega1, Double_t &phi, Double_t &theta) const
double dEdx_AMP(art::Ptr< recob::Hit > hit, double pitch, double T0=0) const
void set_dedx(const std::vector< double > &q)
Definition: Shower.h:139