UpMuRecoAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //// Class: UpMuRecoAna
3 //// Module Type: analyzer
4 //// File: UpMuRecoAna_module.cc
5 ////
6 //////////////////////////////////////////////////////////////////////////
7 
8 
13 
21 #include "fhiclcpp/ParameterSet.h"
22 
23 #include <iostream>
24 #include <set>
25 #include <cmath>
26 #include <algorithm>
27 #include <boost/math/special_functions.hpp>
28 #include <boost/math/special_functions/gamma.hpp>
29 #include "TNtuple.h"
30 
31 #include "Geometry/Geometry.h"
35 #include "MEFinder/MEClusters.h"
36 #include "MCCheater/BackTracker.h"
37 #include "RecoBase/CellHit.h"
38 #include "RecoBase/Track.h"
39 #include "RecoBase/Prong.h"
40 #include "RecoBase/Cluster.h"
41 #include "RecoBase/RecoHit.h"
42 #include "RecoBase/Vertex.h"
43 #include "Calibrator/Calibrator.h"
44 #include "Utilities/func/MathUtil.h"
45 #include "Simulation/FLSHit.h"
46 #include "Simulation/FLSHitList.h"
47 #include "Simulation/Particle.h"
51 
52 namespace upmuana {
53  class UpMuRecoAna;
54 
55  inline double getErr(double PE) {
56  return 231594.0/(2267.16+pow(PE,2.01011)) + 9.55689;
57  }
58 
59  inline double getDist(double x1, double y1, double z1,
60  double x2, double y2, double z2) {
61  return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
62  }
63 
64  bool recoHitComp (rb::RecoHit lhs, rb::RecoHit rhs) {
65  if (lhs.Y() == rhs.Y()) {
66  if (lhs.Z() == rhs.Z()) {
67  return lhs.T() < rhs.T();
68  }
69  else return lhs.Z() < rhs.Z();
70  }
71  else return lhs.Y() < rhs.Y();
72  }
73 
74  bool pairHitComp (std::pair<rb::CellHit, rb::RecoHit> lhp,
75  std::pair<rb::CellHit, rb::RecoHit> rhp) {
76  return recoHitComp(lhp.second, rhp.second);
77  }
78 
79  inline int getDiblock(int plane) { return (plane/64)+1; } // start at 1
80 }
81 
83 public:
84  explicit UpMuRecoAna(fhicl::ParameterSet const & p);
85  virtual ~UpMuRecoAna();
86 
87  void analyze(art::Event const & e) override;
88 
89  void beginJob() override;
90 
91 private:
101 
103 
104  double _TrackLen;
105  unsigned _TrackHitsXY;
106  unsigned _TrackHitsX;
107  unsigned _TrackHitsY;
108  int _dX;
109  int _dY;
110  int _dZ;
111  double _Chi2;
112  double _LLR;
113  double _R2X;
114  double _R2Y;
115  double _MinSlope;
116  double _MaxSlope;
117  bool fIsSim;
120 
122 
123  double cw, pw;
124 
125  double getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
126  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
127  std::pair<rb::CellHit,rb::RecoHit>)> ,
128  std::vector<rb::RecoHit> &outliers,
129  double &P_up, double &P_dn, double &chi2, double &slope);
130  void LLR(std::vector<double>& eT,
131  std::vector<double>& mT,
132  std::vector<double>& mTErr, double& slope, double& chi2,
133  double& P_up, double& P_dn, std::vector<rb::RecoHit> &in,
134  std::vector<rb::RecoHit> &outliers);
135 
136  void LinFit(const std::vector<double>& x,
137  const std::vector<double>& y, double *fitpar);
138  void LinFit(const std::vector<double>& x,
139  const std::vector<double>& y,
140  const std::vector<double>& ye, double *fitpar);
141 
142 
143  struct hit_info {
144  float Run, SubRun, Event, TrackID, HitID, X, Y, Z, T,
147  };
148 
149  struct time_window {
150  float startT;
151  float endT;
152  };
153 
154  float containmentType(rb::Track);
156 
157  TNtuple* fTrackNtuple;
158  TNtuple* fHitNtuple;
159  TNtuple* fVertexNtuple;
160  TNtuple* fEventNtuple;
161  TNtuple* fParticleNtuple;
162  TNtuple* fSliceNtuple;
163  TNtuple* fFLSNtuple;
164  unsigned int fNTracks = 0;
165 };
166 
168  : EDAnalyzer (p)
169  , fHitsLabel (p.get<std::string>("HitModuleLabel" ))
170  , fHitsInstance (p.get<std::string>("HitInstanceLabel"))
171  , fSliceModuleLabel (p.get<std::string>("SliceModuleLabel"))
172  , fSliceInstance (p.get<std::string>("SliceInstanceLabel"))
173  , fTrackModuleLabel (p.get<std::string>("TrackModuleLabel"))
174  , fTrackToSliceInstance(p.get<std::string>("TrackToSliceInstance"))
175  , fFuzzyKVertexLabel (p.get<std::string>("FuzzyKVertexLabel"))
176  , fFuzzyKVertexInstance(p.get<std::string>("FuzzyKVertexInstance"))
177  , fMichelELabel (p.get<std::string>("MichelELabel"))
178  , fMinY (p.get<double>("MinY"))
179  , fMinX (p.get<double>("MinX"))
180  , fMinZ (p.get<double>("MinZ"))
181  , fMaxY (p.get<double>("MaxY"))
182  , fMaxX (p.get<double>("MaxX"))
183  , fMaxZ (p.get<double>("MaxZ"))
184  , _TrackLen (p.get<double>("TrackLen"))
185  , _TrackHitsXY (p.get<unsigned>("TrackHitsXY"))
186  , _TrackHitsX (p.get<unsigned>("TrackHitsX"))
187  , _TrackHitsY (p.get<unsigned>("TrackHitsY"))
188  , _dX (p.get<int>("dX"))
189  , _dY (p.get<int>("dY"))
190  , _dZ (p.get<int>("dZ"))
191  , _Chi2 (p.get<double>("Chi2"))
192  , _LLR (p.get<double>("LLR"))
193  , _R2X (p.get<double>("R2X"))
194  , _R2Y (p.get<double>("R2Y"))
195  , _MinSlope (p.get<double>("MinSlope"))
196  , _MaxSlope (p.get<double>("MaxSlope"))
197  , fIsSim (p.get<bool>("IsSim"))
198  , fParticleModuleLabel (p.get<std::string>("ParticleModuleLabel"))
199  , fParticleInstance (p.get<std::string>("ParticleInstance"))
200  , fMCTruthLabel (p.get<std::string>("MCTruthLabel"))
201 {
202  cw = 3.97; //cell width in cm
203  pw = 6.654; //plane width in cm
204 }
205 
207 {}
208 
210 {
212 
213  fTrackNtuple = tfs->make<TNtuple>("track_ntuple", "Track Ntuple",
214  "Run:SubRun:Event:SliceID:TrackID:Nhits:NRecohits:NOutliers:ProbUp:ProbDn:LLR:Chi2:Slope:LLRX:Chi2X:SlopeX:LLRY:Chi2Y:SlopeY:R2X:R2Y:StartX:StartY:StartZ:StartT:EndX:EndY:EndZ:EndT:TrackHitsX:TrackHitsY:Length:dirX:dirY:dirZ:eleAngle:totalE:containment:avgTX:avgTY:meRecoGeV:meRecoNHits:meRecoX:meRecoY:meRecoZ:meRecoDist:meRecoDeltaT:meRecoAtTrackEnd");
215  fHitNtuple = tfs->make<TNtuple>("hit_ntuple", "Hit Ntuple",
216  "Run:SubRun:Event:TrackID:HitID:X:Y:Z:T:deltaT:eT:PE:PECorr:ADC:Plane:Cell:View:GoodTiming:Reco:Outlier:dirX:dirY:dirZ:eleAngle:GeV:Diblock");
217 
218  fVertexNtuple = tfs->make<TNtuple>("vertex_ntuple", "Vertex Ntuple",
219  "Run:SubRun:Event:Ntracks:totalHits:totalE:containment:X:Y:Z:T");
220 
221  fEventNtuple = tfs->make<TNtuple>("event_ntuple", "Event Ntuple",
222  "Run:SubRun:Event:Ntracks:NcontainedT:Nvertices:NcontainedV:avgNhits:avgNRecohits:avgLen:hasUpMu:hasFCUpMu:hasTGUpMu:nDiblocks");
223 
224  fSliceNtuple = tfs->make<TNtuple>("slice_ntuple", "Slice Ntuple",
225  "Run:SubRun:Event:SliceID:Ntracks:containment");
226 
227  if (fIsSim){
228  fParticleNtuple = tfs->make<TNtuple>("part_ntuple", "Particle Ntuple",
229  "Run:SubRun:Event:ParticleID:PDG:Length:StartX:StartY:StartZ:StartT:StartE:EndX:EndY:EndZ:EndT:EndE:StartPX:StartPY:StartPZ:StartEleAngle:containment:MotherID:TrackID:EDep:nFLS");
230  /*
231  fFLSNtuple = tfs->make<TNtuple>("fls_ntuple", "FLS Ntuple",
232  "Run:SubRun:Event:TrackID:PDG:Cell:Plane:EntryX:EntryY:EntryZ:EntryT:ExitX:ExitY:ExitZ:ExitT:EDep:XAverage:YAverage:ZAverage");*/
233  }
234 }
235 
237 {
238  // getting all hits for backtracker
239 
241  e.getByLabel( fSliceModuleLabel, slicecol);
243 
244  for(int i = 0; i < (int)slicecol->size(); i++){
245  art::Ptr<rb::Cluster> sli(slicecol, i);
246  slces.push_back(sli);
247  }
248 
249 
250  // A vector of all hits in the event is needed by backtracker
251 
252  int nSlices = slces.size();
253 
255  for(int s = 0; s < nSlices; s++){
256  for(int c = 0; c<(int)slces[s]->NCell(); c++)
257  allHits.push_back( slces[s]->Cell(c) );
258  }
259  //seems redundent to get all hits twice in two different forms,
260  //need to improve this
261 
262 
265  //art::PtrVector< rb::CellHit> allHits;
266  e.getByLabel(fHitsLabel, fHitsInstance, hits);
267 
268  //art::Handle<std::vector<rb::Cluster> > slices;
269  //e.getByLabel(fSliceModuleLabel, fSliceInstance, slices);
270 
272  e.getByLabel(fTrackModuleLabel, "", tracks);
273  art::FindOneP<rb::Cluster> fmSlice(tracks, e,
274  art::InputTag(fTrackModuleLabel,""));
275 
276  // art::FindOneP<rb::Cluster> fohl(tracks, e, fTrackModuleLabel);
277  // assert(fohl.size() == tracks->size());
278 
279  art::Handle<std::vector<rb::Prong> > fuzzyKprongs;
280  e.getByLabel(fFuzzyKVertexLabel, fFuzzyKVertexInstance, fuzzyKprongs);
281  art::FindManyP<rb::Vertex> fmVert(fuzzyKprongs, e,
282  art::InputTag(fFuzzyKVertexLabel,fFuzzyKVertexInstance));
283  std::vector<art::Ptr<rb::Vertex> > vertices;
284 
285 
286  float hasUpMu=0, hasFCUpMu=0, hasTGUpMu=0;
287 
288  size_t nContainedT = 0, nVertices = 0, nContainedV = 0;
289  float avgNhits = 0, avgNRecohits = 0, avgLen = 0;
290 
291  for (size_t i_prong=0; i_prong<fuzzyKprongs->size(); ++i_prong) {
292  auto theVertexVec = fmVert.at(i_prong);
293 
294  for(size_t i_vecvert=0; i_vecvert < theVertexVec.size(); ++i_vecvert) {
295 
296  art::Ptr<rb::Vertex> theVertex = theVertexVec.at(i_vecvert);
297 
298  bool seen = false;
299  // check if this vertex has been seen before
300  for (size_t i_vert=0; i_vert<vertices.size(); ++i_vert) {
301  if (&(*vertices.at(i_vert)) == &(*theVertex)) {
302  seen = true; break;
303  }
304  else
305  vertices.push_back(theVertex);
306  }
307  if (seen) continue;
308 
309  ++nVertices;
310  /* float nTracks = theVertex->NTracks(),*/
311  float totalHits=0,
312  vertTotalE=0, containment=4; // assume fully contained
313  /*
314  //std::cout << "Clusters: " << theVertex->NClusters() << std::endl;
315 
316  for (size_t i_vprong=0; i_vprong<nTracks; ++i_vprong) {
317  rb::Track theProng = *theVertex->Track(i_vprong);
318  //if (containmentType(theProng) != 4)
319  // containment = 3; // in-produced
320 
321  totalHits += theProng.NCell();
322 
323  // get total E
324  for (size_t i_hit=0; i_hit<theProng.NCell(); ++i_hit) {
325  rb::RecoHit theRecoHit = theProng.RecoHit(i_hit);
326 
327  if (!theRecoHit.IsCalibrated()) continue;
328 
329  vertTotalE += theRecoHit.GeV();
330  } // i_hit
331  } // i_vprong
332  */
333 
334  // simplified containment criteria: is the vertex inside the detector?
335  if (theVertex->GetX() < fMaxX && theVertex->GetX() > fMinX &&
336  theVertex->GetY() < fMaxY && theVertex->GetY() > fMinY &&
337  theVertex->GetZ() < fMaxZ && theVertex->GetZ() > fMinZ)
338  containment = 4;
339  else containment = 1;
340 
341  if (containment == 4) ++nContainedV;
342 
343  float vertex_entries[11] = {(float)e.id().run(), (float)e.id().subRun(),
344  (float)e.id().event(), 0/*nTracks*/, totalHits, vertTotalE, containment,
345  (float)theVertex->GetX(), (float)theVertex->GetY(),
346  (float)theVertex->GetZ(), (float)theVertex->GetT()};
347  fVertexNtuple->Fill(vertex_entries);
348  } // i_vecvert
349  } // i_prong
350 
351  unsigned int nGoodTiming = 0;
352  for (size_t i_hit=0; i_hit < hits->size(); ++i_hit) {
353  if (hits->at(i_hit).GoodTiming()) ++nGoodTiming;
354  }
355 
356  //std::cout << "Number of hits: " << hits->size() << std::endl;
357  //std::cout << "Number of hits with fine timing: " << nGoodTiming << std::endl;
358  //std::cout << "Number of tracks: " << tracks->size() << std::endl;
359  //std::cout << "Number of vertices: " << nVertices << std::endl;
360 
361  std::vector<const rb::Cluster*> slices;
362  std::vector<unsigned> tracks_per_slice;
363  std::vector<float> slice_containment;
364 
365  std::vector<hit_info> allHitInfo; //ntuple info for all hits
366  std::vector<time_window> goodTimeWindows;
367 
368  // in case we want all hit info
369  //time_window full_event_window;
370  //full_event_window.startT = -5e6;
371  //full_event_window.endT = 5e6;
372  //goodTimeWindows.push_back(full_event_window);
373  art::FindManyP<me::TrkME> fo(tracks, e, fMichelELabel);
374  std::vector<float> activeDiblocks; //keeps track of diblocks with a hit
375  for (size_t i_track=0; i_track < tracks->size(); ++i_track) {
376  bool keepHits = true;
377  rb::Track theTrack = tracks->at(i_track);
378  float containment = containmentType(theTrack);
379  // find the slice number, if any
380  float slice=-1;
381  bool found_slice = false;
382  const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
383  for (size_t i_slice=0; i_slice<slices.size(); ++i_slice) {
384  //std::cout << "i_slice " << i_slice << std::endl;
385  if (theSlice == slices.at(i_slice))
386  { slice = (float)i_slice; found_slice = true;
387  ++tracks_per_slice.at(i_slice);
388  if(containment != 4) slice_containment.at(i_slice) = 1;
389  break; }
390  }
391  if (!found_slice) {
392  slice = (float)slices.size();
393  slices.push_back(theSlice);
394  tracks_per_slice.push_back(1);
395  if(containment!=4)
396  slice_containment.push_back(1);
397  else
398  slice_containment.push_back(4);
399  }
400 
401  if (!theTrack.Is3D()) continue;
402 
403  //std::cout << "3D track id: " << theTrack.ID() << std::endl;
404 
405  // get hits in track
406  art::PtrVector<rb::CellHit> trackHits = theTrack.AllCells();
407  std::set< std::pair<rb::CellHit,rb::RecoHit>,
408  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
409  std::pair<rb::CellHit,rb::RecoHit>)>
410  sortedTrackHits (pairHitComp);
411  std::set< std::pair<rb::CellHit,rb::RecoHit>,
412  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
413  std::pair<rb::CellHit,rb::RecoHit>)>
414  sortedTrackHitsX (pairHitComp);
415  std::set< std::pair<rb::CellHit,rb::RecoHit>,
416  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
417  std::pair<rb::CellHit,rb::RecoHit>)>
418  sortedTrackHitsY (pairHitComp);
419 
420  std::vector<double> x_hit; //maybe should change to short?
421  std::vector<double> y_hit;
422  std::vector<double> zy_hit;
423  std::vector<double> zx_hit;
424  //Get Michel Cluster associated with track
425  std::vector<art::Ptr<me::TrkME> > michelE = fo.at(i_track);
426  //initialize to garbage values
427  //will change if ME is found
428  float meRecoGeV = -5;
429  float meRecoHits = -5;
430  float meRecoPos[3] = {-5, -5, -5};
431  float meRecoDist = -5;
432  float meRecoDeltaT = -5;
433  float meAtEnd = 1;
434  if(!michelE.empty()){
435  meRecoGeV = michelE[0]->TotalGeV();
436  meRecoHits = michelE[0]->NCell();
437  meRecoPos[0] = michelE[0]->MeanX();
438  meRecoPos[1] = michelE[0]->MeanY();
439  meRecoPos[2] = michelE[0]->MeanZ();
440  meRecoDist = michelE[0]->DistToTrk();
441  meRecoDeltaT = michelE[0]->DeltaT();
442  } //end Michel if
443 
444 
445  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
446  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
447  if (!theHit->GoodTiming()) continue;
448  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
449  if (!theRecoHit.IsCalibrated()) continue;
450  sortedTrackHits.insert(std::make_pair(*theHit,theRecoHit));
451  unsigned short iC = theHit->Cell();
452  unsigned short iP = theHit->Plane();
453  geo::View_t view = theHit->View();
454  if(view == geo::kY){
455  x_hit.push_back(iC);
456  zx_hit.push_back(iP);
457  sortedTrackHitsY.insert(std::make_pair(*theHit,theRecoHit));
458  }
459  else if(view == geo::kX){
460  y_hit.push_back(iC);
461  zy_hit.push_back(iP);
462  sortedTrackHitsX.insert(std::make_pair(*theHit,theRecoHit));
463  }
464  } // i_hit loop
465 
466  // std::cout << "***Sorted RecoHits***" << std::endl;
467  // for (auto i_hit=sortedTrackHits.begin();
468  // i_hit != sortedTrackHits.end(); ++i_hit)
469  // std::cout << "Y = " << i_hit->Y() << std::endl;
470  // std::cout << "*********************" << std::endl;
471 
472  double fitpar[3];
473  LinFit(x_hit, zx_hit, fitpar);
474  double r2x = fitpar[2];
475  LinFit(y_hit, zy_hit, fitpar);
476  double r2y = fitpar[2];
477  unsigned nxhit = x_hit.size();
478  unsigned nyhit = y_hit.size();
479 
480  std::vector<rb::RecoHit> outliers, outliersX, outliersY;
481  double P_up, P_dn, chi2, slope;
482 
483  double llr = getLLR(sortedTrackHits, outliers, P_up, P_dn, chi2,
484  slope);
485 
486  double chi2X, slopeX, P_upX, P_dnX, chi2Y, slopeY, P_upY, P_dnY,
487  llrX = getLLR(sortedTrackHitsX, outliersX, P_upX, P_dnX,
488  chi2X, slopeX),
489  llrY = getLLR(sortedTrackHitsY, outliersY, P_upY, P_dnY,
490  chi2Y, slopeY);
491 
494  if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
495  else { start = rb::RecoHit(); end = rb::RecoHit(); }
496  if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->second;
497  else end = start;
498  float dist = theTrack.TotalLength();
499  //getDist(start.X(), start.Y(), start.Z(),
500  //end.X(), end.Y(), end.Z());
501  //std::cout << "P_up: " << P_up << ", P_dn: " << P_dn << ", LLR: " << llr
502  // << std::endl;
503 
504  // averages of all hit directions
505  float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
506 
507  //Checks to see if hits should be kept
508  if(chi2 > _Chi2) keepHits = false;
509  if(dist < _TrackLen) keepHits = false;
510  if(TMath::Abs(start.X() - end.X()) < _dX*cw) keepHits = false;
511  if(TMath::Abs(start.Y() - end.Y()) < _dY*cw) keepHits = false;
512  if(TMath::Abs(start.Z() - end.Z()) < _dZ*pw) keepHits = false;
513  if(r2x < _R2X) keepHits = false;
514  if(r2y < _R2Y) keepHits = false;
515  if(nxhit < _TrackHitsX) keepHits = false;
516  if(nyhit < _TrackHitsY) keepHits = false;
517  if((nxhit+nyhit) < _TrackHitsXY) keepHits = false;
518 
519  if(keepHits){
521  window.startT = start.T() - 2000;
522  window.endT = end.T() + 2000;
523  goodTimeWindows.push_back(window);
524  }
525 
526  float avgTX=0, avgTY=0;
527 
528  // now fill hit ntuple
529  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
530  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
531  float goodtime = theHit->GoodTiming() ? 1 : 0;
532  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
533  float PE = theHit->PE(),
534  reco = theRecoHit.IsCalibrated() ? 1 : 0,
535  PECorr = theRecoHit.IsCalibrated() ? theRecoHit.PECorr() :
536  theHit->PE(),
537  outlier = 0;
538  TVector3 theDir = theTrack.InterpolateDir(theRecoHit.Z());
539  float eleAngle = atan(theDir.y()/sqrt(theDir.x()*theDir.x() +
540  theDir.z()*theDir.z()));
541 
542  tDirX += theDir.x()/trackHits.size();
543  tDirY += theDir.y()/trackHits.size();
544  tDirZ += theDir.z()/trackHits.size();
545 
546  float GeV = 0;
547 
548  // only calibrated hits are used for the fit
549  if (reco) {
550  GeV = theRecoHit.GeV();
551  for (size_t i_out=0; i_out < outliers.size(); ++i_out) {
552  if (!outliers.at(i_out).IsCalibrated()) continue;
553  if (outliers.at(i_out).PECorr() == theRecoHit.PECorr() &&
554  outliers.at(i_out).X() == theRecoHit.X() &&
555  outliers.at(i_out).Y() == theRecoHit.Y() &&
556  outliers.at(i_out).Z() == theRecoHit.Z() &&
557  outliers.at(i_out).T() == theRecoHit.T())
558  { outlier = 1; break; }
559  }
560  }
561 
562  trackTotalE += GeV;
563  if (theHit->View() == geo::kX)
564  avgTX += theRecoHit.T()/sortedTrackHitsX.size();
565  else if (theHit->View() == geo::kY)
566  avgTY += theRecoHit.T()/sortedTrackHitsY.size();
567 
568  float deltaT = (float) getErr(PE);
569  float eT = (float)(getDist(theRecoHit.X(), theRecoHit.Y(), theRecoHit.Z(),
570  start.X(), start.Y(), start.Z())/29.97),
571  ADC = theHit->ADC();
572  //store hit info for later
573  hit_info info;
574  info.Run = (float)e.id().run();
575  info.SubRun = (float)e.id().subRun();
576  info.Event = (float)e.id().event();
577  info.TrackID = (float)i_track;
578  info.HitID = (float)i_hit;
579  info.X = (float)theRecoHit.X();
580  info.Y = (float)theRecoHit.Y();
581  info.Z = (float)theRecoHit.Z();
582  info.T = (float)theRecoHit.T() - sortedTrackHits.begin()->second.T();
583  info.deltaT = deltaT;
584  info.eT = eT;
585  info.PE = PE;
586  info.PECorr = PECorr;
587  info.ADC = ADC;
588  info.Plane = (float)theHit->Plane();
589  info.Cell = (float)theHit->Cell();
590  info.View = (float)theHit->View();
591  info.GoodTiming = goodtime;
592  info.Reco = reco;
593  info.Outlier = outlier;
594  info.dirX = (float)theDir.x();
595  info.dirY = (float)theDir.y();
596  info.dirZ = (float)theDir.z();
597  info.eleAngle = eleAngle;
598  info.GeV = GeV;
599  info.absT = theRecoHit.T();
600  info.Diblock = getDiblock(theHit->Plane());
601  allHitInfo.push_back(info);
602 
603  if(std::find(activeDiblocks.begin(),activeDiblocks.end(),
604  info.Diblock)==activeDiblocks.end())
605  activeDiblocks.push_back(info.Diblock);
606  } // i_hit loop
607 
608 
609  tEleAngle = atan(tDirY/sqrt(tDirX*tDirX + tDirZ*tDirZ));
610 
611  // update event-wide variables
612  if (containment == 4) ++nContainedT;
613  avgNhits += trackHits.size()/tracks->size();
614  avgNRecohits += sortedTrackHits.size()/tracks->size();
615  avgLen += dist/tracks->size();
616 
617  float track_entries[48]=
618  { (float)e.id().run(), (float)e.id().subRun(),
619  (float)e.id().event(), slice, (float)i_track,
620  (float)trackHits.size(), (float)sortedTrackHits.size(),
621  (float)outliers.size(), (float)P_up, (float)P_dn, (float)llr,
622  (float)chi2, (float)slope, (float)llrX, (float)chi2X, (float)slopeX,
623  (float)llrY, (float)chi2Y, (float)slopeY, (float)r2x, (float)r2y,
624  start.X(), start.Y(), start.Z(), start.T(), end.X(), end.Y(),
625  end.Z(), end.T(), (float)nxhit,(float)nyhit, dist, tDirX, tDirY,
626  tDirZ, tEleAngle, trackTotalE, containment, avgTX, avgTY,
627  meRecoGeV, meRecoHits, meRecoPos[0], meRecoPos[1],
628  meRecoPos[2], meRecoDist, meRecoDeltaT, meAtEnd
629  };
630  fTrackNtuple->Fill(track_entries);
631 
632  } // i_track loop
633 
634  // fill slice ntuple
635  for (size_t i_slice = 0; i_slice < slices.size(); ++i_slice) {
636  float slice_entries[6] = { (float)e.id().run(), (float)e.id().subRun(),
637  (float)e.id().event(), (float)i_slice,
638  (float)tracks_per_slice.at(i_slice),
639  (float)slice_containment.at(i_slice)};
640  fSliceNtuple->Fill(slice_entries);
641  }
642 
643  for(size_t i_hitInfo = 0; i_hitInfo < allHitInfo.size(); ++i_hitInfo){
644  hit_info info = allHitInfo.at(i_hitInfo);
645  bool added = false;
646  for(size_t i_window = 0; i_window < goodTimeWindows.size(); i_window++){
647  if (added) break;
648  time_window window = goodTimeWindows.at(i_window);
649  if(info.absT > window.startT && info.absT < window.endT){
650  float hit_entries[26] = {info.Run, info.SubRun, info.Event,
651  info.TrackID, info.HitID, info.X,
652  info.Y, info.Z, info.T, info.deltaT,
653  info.eT, info.PE, info.PECorr, info.ADC,
654  info.Plane, info.Cell,
655  info.View, info.GoodTiming, info.Reco,
656  info.Outlier, info.dirX, info.dirY,
657  info.dirZ, info.eleAngle, info.GeV, info.Diblock};
658  fHitNtuple->Fill(hit_entries);
659  added = true;
660  }
661  }// i_window
662  }// i_hitInfo
663 
665  if (fIsSim) {
666  e.getByLabel(fParticleModuleLabel, fParticleInstance, particles);
667  std::vector<sim::Particle> michelElectrons;
668  for (size_t i_part=0; i_part<particles->size(); ++i_part) {
669  sim::Particle thePart = particles->at(i_part);
670  float startEleAngle = atan(thePart.Py() / sqrt(thePart.Px()*thePart.Px() +
671  thePart.Pz()*thePart.Pz())),
672  containment = containmentType(thePart);
673 
674  if ((thePart.PdgCode()==13 || thePart.PdgCode()==-13) &&
675  thePart.Py() > 0 && fabs(startEleAngle)*180/M_PI > 10.0 &&
676  thePart.E() > 0.5) {
677  hasUpMu = 1;
678  if (containment == 4) hasFCUpMu=1;
679  if (containment == 1) hasTGUpMu=1;
680  }
681 
682  float length =
683  sqrt((thePart.Vx()-thePart.EndX())*(thePart.Vx()-thePart.EndX())
684  +(thePart.Vy()-thePart.EndY())*(thePart.Vy()-thePart.EndY())
685  +(thePart.Vz()-thePart.EndZ())*(thePart.Vz()-thePart.EndZ()));
686 
687  //std::vector<sim::FLSHit> flsHits= bt->ParticleToFLSHit(thePart.TrackId(), thePart.PdgCode());
688  std::vector<sim::FLSHit> flsHits= bt->ParticleToFLSHit(thePart.TrackId());
689 
690  //float detLength = -1;
691  //float entryT = 9000000;
692  //float exitT = 0;
693  //sim::FLSHit first;
694  //sim::FLSHit last;
695  float eDep = 0;
696  for(size_t i_fls = 0; i_fls < flsHits.size(); i_fls++){
697  sim::FLSHit h = flsHits[i_fls];
698  eDep += h.GetEdep();
699  /*
700  float T = (h.GetEntryT() + h.GetExitT()) / 2; //avg T of flshit
701 
702  if(T < entryT){
703  entryT = T;
704  first = h;
705  }
706  if(T > exitT){
707  exitT = T;
708  last = h;
709  }
710 
711  float fls_entries[19] =
712  {(float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
713  (float)h.GetTrackID(), (float)h.GetPDG(), (float)h.GetCellID(),
714  (float)h.GetPlaneID(), (float)h.GetEntryX(), (float)h.GetEntryY(),
715  (float)h.GetEntryZ(), (float)h.GetEntryT(), (float)h.GetExitX(),
716  (float)h.GetExitY(), (float)h.GetExitZ(), (float)h.GetExitT(),
717  (float)h.GetEdep(), (float)h.GetXAverage(), (float)h.GetYAverage(),
718  (float)h.GetZAverage()
719  };
720  */
721  // fFLSNtuple->Fill(fls_entries);
722 
723  } //end loop on FLS hits
724 
725  /*
726  if(flsHits.size() > 1){ //particles often have no hits
727  detLength = sqrt((entryX - exitX)*(entryX - exitX) +
728  (entryY - exitY)*(entryY - exitY) +
729  (entryZ - exitZ)*(entryZ - exitZ));
730  }
731 
732  if(flsHits.size() > 5){
733  UInt_t firstPlane = first.GetPlaneID();
734  UInt_t lastPlane = last.GetPlaneID();
735  UInt_t firstCell = first.GetCellID();
736  UInt_t lastCell = last.GetCellID();
737  art::ServiceHandle<geo::Geometry> geom;
738  Double_t firstPos[4];
739  //double firstW;
740  Double_t lastPos[4];
741  //double lastW;
742  const geo::PlaneGeo* firstP = geom->Plane(firstPlane);
743  const geo::PlaneGeo* lastP = geom->Plane(lastPlane);
744  const geo::CellGeo* firstC = firstP->Cell(firstCell);
745  const geo::CellGeo* lastC = lastP->Cell(lastCell);
746 
747  firstC->GetCenter(firstPos);
748  lastC->GetCenter(lastPos);
749 
750  TVector3 traj = {thePart.Vx()-thePart.EndX(),
751  thePart.Vy()-thePart.EndY(),
752  thePart.Vz()-thePart.EndZ()};
753 
754  double scale = (firstPos[2]-lastPos[2])/traj.Z();
755  TVector3 path = traj*scale;
756  detLength = path.Mag();
757  }
758  }*/
759  float particle_entries[25] =
760  { (float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
761  (float)i_part, (float)thePart.PdgCode(), length,
762  (float)thePart.Vx(), (float)thePart.Vy(), (float)thePart.Vz(),
763  (float)thePart.T(), (float)thePart.E(), (float)thePart.EndX(),
764  (float)thePart.EndY(), (float)thePart.EndZ(),
765  (float)thePart.EndT(),(float)thePart.EndE(),
766  (float)thePart.Px(), (float)thePart.Py(), (float)thePart.Pz(),
767  startEleAngle, containment,(float)thePart.Mother(),
768  (float)thePart.TrackId(), eDep, (float) flsHits.size()};
769  fParticleNtuple->Fill(particle_entries);
770  }
771  }
772 
773  float event_entries[14] = { (float)e.id().run(), (float)e.id().subRun(),
774  (float)e.id().event(), (float)tracks->size(), (float)nContainedT,
775  (float)nVertices, (float)nContainedV, avgNhits, avgNRecohits,
776  avgLen, hasUpMu, hasFCUpMu, hasTGUpMu,
777  (float) activeDiblocks.size() };
778  fEventNtuple->Fill(event_entries);
779 
780 
781 } // analyze() function
782 
783 
784 double upmuana::UpMuRecoAna::getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
785  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
786  std::pair<rb::CellHit,rb::RecoHit>)>
787  sortedHits,
788  std::vector<rb::RecoHit> &outliers,
789  double &P_up, double &P_dn, double &chi2,
790  double &slope)
791 {
792  if (sortedHits.size() < 2) {
793  P_up = 1e-30;
794  P_dn = 1e-30;
795  chi2 = 1e6;
796  slope = -1e6;
797  return 0;
798  }
799 
800  std::vector<double> measuredTimes;
801  std::vector<double> expectedTimes;
802  std::vector<double> wgts;
803 
804  // first hit
805  measuredTimes.push_back(0.0);
806  expectedTimes.push_back(0.0);
807  double err = getErr(sortedHits.begin()->first.PE());
808  wgts.push_back(1.0/err/err);
809 
810  double startX = sortedHits.begin()->second.X(),
811  startY = sortedHits.begin()->second.Y(),
812  startZ = sortedHits.begin()->second.Z(),
813  startT = sortedHits.begin()->second.T();
814 
815  std::vector<rb::RecoHit> in;
816 
817  for (auto i_hit = ++sortedHits.begin();
818  i_hit != sortedHits.end();
819  ++i_hit) {
820  in.push_back(i_hit->second);
821  measuredTimes.push_back(i_hit->second.T() - startT);
822  double dist = getDist(i_hit->second.X(), i_hit->second.Y(),
823  i_hit->second.Z(), startX, startY, startZ);
824  expectedTimes.push_back(dist/29.97);
825  err = getErr(i_hit->first.PE());
826  wgts.push_back(1.0/err/err);
827  }
828 
829  LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
830  in, outliers);
831 
832  return log(P_up/P_dn);
833 }
834 
835 void upmuana::UpMuRecoAna::LLR(std::vector<double>& eT,
836  std::vector<double>& mT,
837  std::vector<double>& wgts, double& slope,
838  double& chi2, double& P_up, double& P_dn,
839  std::vector<rb::RecoHit>& in,
840  std::vector<rb::RecoHit>& outliers)
841 {
842  // eT - x, mT - y
843  size_t n = mT.size();
844 
845  // y = a + bx
846  double a, b;
847  if (eT.size() >= 2)
848  chi2 = util::LinFit(eT, mT, wgts, b, a);
849  else {
850  chi2 = 1e6;
851  slope = -1e6;
852  P_up = 1e-30;
853  P_dn = 1e-30;
854  return;
855  }
856  //a = fitpar[0];
857  //b = fitpar[1];
858 
859  size_t totAllowOutlier = n/10;
860  size_t nOutliers = 0;
861  std::vector<double> x_filt, y_filt, w_filt;
862  for (size_t i=0; i < n; i++) {
863  double y_fit = a + b*eT.at(i);
864  if ( fabs(mT.at(i) - y_fit) < 5*(1.0/sqrt(wgts.at(i)))
865  || nOutliers>totAllowOutlier)
866  {
867  x_filt.push_back(eT.at(i));
868  y_filt.push_back(mT.at(i));
869  w_filt.push_back(wgts.at(i));
870  }
871  else {
872  ++nOutliers;
873  outliers.push_back(in[i]);
874  }
875  }
876 
877  if (x_filt.size() >= 2)
878  chi2 = util::LinFit(x_filt, y_filt, w_filt, b, a);
879  else {
880  chi2 = 1e6;
881  slope = -1e6;
882  P_up = 1e-30;
883  P_dn = 1e-30;
884  }
885  //a = fitpar[0];
886  //b = fitpar[1];
887  n = x_filt.size();
888  if (n < 5) {
889  slope = 0;
890  chi2 = 999;
891  P_up = 1e-30;
892  P_dn = 1;
893  return;
894  }
895 
896  slope = b;
897  //chi2 = 0;
898  //for (size_t i=0; i<n; i++) {
899  // double y_exp = a + b*x_filt.at(i);
900  // chi2 += pow((y_filt.at(i)-y_exp) / ye_filt.at(i), 2.0);
901  //}
902  chi2 /= (double)(n-2);
903 
904  double one_over_ee = 0.0,
905  x_over_ee = 0.0,
906  y_over_ee = 0.0;
907 
908  for (size_t i=0; i<n; ++i) {
909  //double e = ye_filt.at(i);
910  one_over_ee += w_filt.at(i);
911  x_over_ee += x_filt.at(i)*w_filt.at(i);
912  y_over_ee += y_filt.at(i)*w_filt.at(i);
913  }
914 
915  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
916  double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;
917 
918  double up_chi2 = 0.0, dn_chi2 = 0.0;
919  for (size_t i=0; i<n; ++i) {
920  double e = 1.0/sqrt(w_filt.at(i));
921  double up_expected = up_inter + x_filt.at(i);
922  double dn_expected = dn_inter - x_filt.at(i);
923  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
924  dn_chi2 += pow((y_filt.at(i)-dn_expected) / e, 2.0);
925  }
926 
927  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0),
928  prob_dn = boost::math::gamma_q((double)(n-2)/2.0,dn_chi2/2.0);
929 
930  if (prob_up < 1e-30) prob_up = 1e-30;
931  if (prob_dn < 1e-30) prob_dn = 1e-30;
932 
933  P_up = prob_up;
934  P_dn = prob_dn;
935 
936 }
937 
938 /*void upmuana::UpMuRecoAna::LinFit(std::vector<double>& x,
939  std::vector<double>& y,
940  std::vector<double>& y_err,
941  double *fitpar) {
942  const size_t n = x.size();
943  double xsum=0, ysum=0, xxsum=0, yysum=0, xysum=0, eesum=0;
944 
945  for (size_t i=0; i<n; ++i) {
946  xsum += x[i]/y_err[i]/y_err[i];
947  ysum += y[i]/y_err[i]/y_err[i];
948  xxsum += x[i]*x[i]/y_err[i]/y_err[i];
949  yysum += y[i]*y[i]/y_err[i]/y_err[i];
950  xysum += x[i]*y[i]/y_err[i]/y_err[i];
951  eesum += 1./y_err[i]/y_err[i];
952  }
953 
954  const double b = (xysum*eesum - xsum*ysum)/(xxsum*eesum - xsum*xsum);
955  const double a = (xysum - b*xxsum)/xsum;
956  const double r = (xysum - xsum*ysum)/sqrt(xxsum - xsum*xsum)
957  *(yysum - ysum*ysum);
958 
959  fitpar[0] = a;
960  fitpar[1] = b;
961  fitpar[2] = r*r;
962  }*/
963 
965 {
966  TLorentzVector start = TLorentzVector(thePart.Position());
967  TLorentzVector stop = TLorentzVector(thePart.EndPosition());
968 
969  if (start.Y() > stop.Y()) // assume upward-going
970  { TLorentzVector hold = start; start = stop; stop = hold; }
971 
972  if (start.Y() < fMinY || start.X() < fMinX || start.X() > fMaxX ||
973  start.Z() < fMinZ || start.Z() > fMaxZ) { // entered from outside
974  if (stop.Y() > fMaxY ||
975  stop.X() < fMinX ||
976  stop.X() > fMaxX ||
977  stop.Z() < fMinZ ||
978  stop.Z() > fMaxZ) { return 1; }//through-going
979  return 2; // stopped
980  } // started in detector
981  if (stop.Y() > fMaxY || stop.X() < fMinX || stop.X() > fMaxX ||
982  stop.Z() < fMinZ || stop.Z() > fMaxZ)
983  { return 3; }// in-produced
984 
985  return 4; // fully contained
986 
987 }
988 
990 {
991  TVector3 start = theTrack.Start();
992  TVector3 stop = theTrack.Stop();
993 
994  //std::cout << "Start: (" << start.X() <<", " << start.Y() << ", " << start.Z() << ")" << std::endl;
995  //std::cout << "Stop: (" << stop.X() << ", " << stop.Y() << ", " << stop.Z() << ")" << std::endl;
996 
997  if (start.y() > stop.y()) // assume upward-going, so switch
998  { TVector3 hold = start; start = stop; stop = hold; }
999 
1000  //std::cout << "Start: (" << start.X() <<", " << start.Y() << ", " << start.Z() << ")" << std::endl;
1001  //std::cout << "Stop: (" << stop.X() << ", " << stop.Y() << ", " << stop.Z() << ")" << std::endl;
1002 
1003  if (start.y() < fMinY || start.x() < fMinX || start.x() > fMaxX ||
1004  start.z() < fMinZ || start.z() > fMaxZ) { // entered from outside
1005  if (stop.y() > fMaxY ||
1006  stop.x() < fMinX ||
1007  stop.x() > fMaxX ||
1008  stop.z() < fMinZ ||
1009  stop.z() > fMaxZ) { return 1; }//through-going
1010  return 2; // stopped
1011  } // started in detector
1012  if (stop.y() > fMaxY || stop.x() < fMinX || stop.x() > fMaxX ||
1013  stop.z() < fMinZ || stop.z() > fMaxZ)
1014  { return 3; }// in-produced
1015 
1016  return 4; // fully contained
1017 
1018 }
1019 
1020 
1021 void upmuana::UpMuRecoAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
1022  // http://en.wikipedia.org/wiki/Simple_linear_regression
1023  const auto n = x_val.size();
1024  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
1025  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
1026  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n; // <xx>
1027  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n; // <yy>
1028  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n; // <xy>
1029 
1030  const auto b = (xy - x*y)/(xx - x*x); // slope
1031  const auto a = y - b*x; // intercept
1032  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
1033  fitpar[0] = a;
1034  fitpar[1] = b;
1035  fitpar[2] = r*r;
1036 
1037 }
1038 
1039 void upmuana::UpMuRecoAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, const std::vector<double>& y_err, double *fitpar){
1040  // http://en.wikipedia.org/wiki/Simple_linear_regression
1041  int n = x_val.size();
1042  double x = 0;
1043  double y = 0;
1044  double xx = 0;
1045  double yy = 0;
1046  double xy = 0;
1047  double ee = 0;
1048 
1049  for ( int i=0; i<n; i++ ){
1050 
1051  x = x + x_val[i]/y_err[i]/y_err[i];
1052  y = y + y_val[i]/y_err[i]/y_err[i];
1053 
1054  xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
1055  yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
1056  xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
1057  ee = ee + 1./y_err[i]/y_err[i];
1058  }
1059 
1060  const auto b = (xy*ee - x*y)/(xx*ee - x*x); // slope
1061  const auto a = (xy - b*xx)/x; // intercept
1062  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
1063  fitpar[0] = a;
1064  fitpar[1] = b;
1065  fitpar[2] = r*r;
1066 
1067 
1068 }
1069 
float T() const
Definition: RecoHit.h:39
double E(const int i=0) const
Definition: MCParticle.h:232
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
const XML_Char XML_Encoding * info
Definition: expat.h:530
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
Double_t xx
Definition: macro.C:12
double Py(const int i=0) const
Definition: MCParticle.h:230
double EndE() const
Definition: MCParticle.h:243
float containmentType(rb::Track)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
double EndZ() const
Definition: MCParticle.h:227
int getDiblock(int plane)
Float_t y1[n_points_granero]
Definition: compare.C:5
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double GetX() const
Definition: Vertex.h:23
unsigned short Plane() const
Definition: CellHit.h:39
Float_t x1[n_points_granero]
Definition: compare.C:5
int Mother() const
Definition: MCParticle.h:212
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double Px(const int i=0) const
Definition: MCParticle.h:229
float Z() const
Definition: RecoHit.h:38
UpMuRecoAna(fhicl::ParameterSet const &p)
A collection of associated CellHits.
Definition: Cluster.h:47
void LLR(std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
double GetY() const
Definition: Vertex.h:24
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
RunNumber_t run() const
Definition: EventID.h:98
double EndY() const
Definition: MCParticle.h:226
double chi2()
double getErr(double PE)
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double GetZ() const
Definition: Vertex.h:25
int TrackId() const
Definition: MCParticle.h:209
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double dist
Definition: runWimpSim.h:113
T atan(T number)
Definition: d0nt_math.hpp:66
std::string fTrackToSliceInstance
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
std::string fFuzzyKVertexInstance
void hits()
Definition: readHits.C:15
length
Definition: demo0.py:21
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double getLLR(std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>, std::vector< rb::RecoHit > &outliers, double &P_up, double &P_dn, double &chi2, double &slope)
const double a
bool recoHitComp(rb::RecoHit lhs, rb::RecoHit rhs)
double T(const int i=0) const
Definition: MCParticle.h:223
float PE() const
Definition: CellHit.h:42
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
static constexpr Double_t GeV
Definition: Munits.h:291
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
reference at(size_type n)
Definition: PtrVector.h:365
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double EndT() const
Definition: MCParticle.h:228
Vertex location in position and time.
size_type size() const
Definition: PtrVector.h:308
EventNumber_t event() const
Definition: EventID.h:116
void analyze(art::Event const &e) override
bool pairHitComp(std::pair< rb::CellHit, rb::RecoHit > lhp, std::pair< rb::CellHit, rb::RecoHit > rhp)
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
SubRunNumber_t subRun() const
Definition: EventID.h:110
T * make(ARGS...args) const
ifstream in
Definition: comparison.C:7
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float X() const
Definition: RecoHit.h:36
double Pz(const int i=0) const
Definition: MCParticle.h:231
const hit & b
Definition: hits.cxx:21
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
double Vz(const int i=0) const
Definition: MCParticle.h:222
virtual TVector3 InterpolateDir(double z) const
Definition: Track.cxx:348
TRandom3 r(0)
float Y() const
Definition: RecoHit.h:37
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
double GetT() const
Definition: Vertex.h:26
float PECorr() const
Definition: RecoHit.cxx:47
bool GoodTiming() const
Definition: CellHit.h:48
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
double EndX() const
Definition: MCParticle.h:225
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
Float_t e
Definition: plot.C:35
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12
size_t size() const
EventID id() const
Definition: Event.h:56
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71