NumuSandFxs.cxx
Go to the documentation of this file.
2 
4 #include "Geometry/Geometry.h"
7 #include "GeometryObjects/Geo.h"
8 #include "NovaDAQConventions/DAQConventions.h"
9 #include "RecoBase/RecoHit.h"
10 #include "MCCheater/BackTracker.h"
11 #include "RecoBase/Cluster.h"
13 
14 #include "TMath.h"
15 #include "TVector3.h"
16 
18 
19 #define PI 3.14159265
20 
21 namespace numusand
22 {
23 
24  //-------------------------------------------------------
25  NumuSandFxs::NumuSandFxs() //constructor
26  {
27  fKGeo = new trk::KalmanGeoHelper(true);
28  }
29 
30  //-------------------------------------------------------
32  {
33  }
34 
35  //-------------------------------------------------------
36  // getAngle: function to get the angle of the lepton wrt the incoming neutrino.
37  // currently the incoming neutrino direction estimate is hardcoded.
38 
39  float NumuSandFxs::getAngle(TVector3 trackdir)
40  {
42 
43  TVector3 beamdir = geom->NuMIBeamDirection();
44  float ret = trackdir*beamdir;
45 
46  return ret;
47  }
48 
49  //---------------------------------------------------------------------------
50 
51  float NumuSandFxs::fitslope(const std::vector< double > &l, const std::vector< double > &t, const std::vector< double > &q, double &chisq, double &inter, bool linflag) {
52 
53  assert(t.size()==l.size());
54  assert(t.size()>=2);
55 
56  double ztStart[2] = {0.};
57  double ztEnd[2] = {0.};
58  if(linflag==false) {
59  chisq = geo::LinFit(t, l, q,
60  &ztStart[0], &ztStart[1],
61  &ztEnd[0], &ztEnd[1]);
62  }
63  else if(linflag==true) {
64  chisq = geo::LinFitMinDperp(t, l, q,
65  &ztStart[0], &ztStart[1],
66  &ztEnd[0], &ztEnd[1]);
67  }
68  const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);
69 
70  inter = ztStart[1]-(ztSlope*ztStart[0]);
71 
72  return ztSlope;
73  }
74 
75  //----------------------------------------------------------------------------------------
76 
77  float NumuSandFxs::fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c) {
78 
79  assert(t.size()==l.size());
80  assert(t.size()>=2);
81 
82  double chi2 = 0.0;
83  double tave = 0;
84  double lave = 0;
85 
86  for (unsigned int i=1; i<t.size(); ++i) {
87  tave += t[i];
88  lave += l[i];
89  }
90 
91  double norm = t.size();
92 
93  tave = tave / norm;
94  lave = lave / norm;
95 
96  double tstart = tave - 500;
97  double tend = tave + 500;
98  double lstart, lend;
99 
100  assert(flag == 0 || flag == 1);
101  if(flag==0) {
102  lstart = lave - (c * 500);
103  lend = lave + (c * 500);
104  }
105  else {
106  lstart = lave + (c * 500);
107  lend = lave - (c * 500);
108  }
109 
110  for (unsigned int i=0; i<t.size(); ++i) {
111  chi2 += geo::DsqrToLine(t[i],l[i],tstart,lstart,tend,lend);
112  }
113  return chi2;
114  }
115 
116  //----------------------------------------------------------------------------------
117 
119  double tscattmax=0, tscatttot=0, tscatt;
120  TVector3 traj1, traj2, dir, olddir;
121 
122  for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) {
123  traj1 = track->TrajectoryPoint(p);
124  traj2 = track->TrajectoryPoint(p+1);
125  dir = (traj2-traj1).Unit();
126  if((dir-olddir).Mag()<1e-7) continue;
127  if(p>0&&dir.Mag()<=1.0&&olddir.Mag()<=1.0) {
128  tscatt = std::acos(dir.Dot(olddir)) * 180 / PI;
129  if(tscatt > tscattmax) tscattmax = tscatt;
130  tscatttot += tscatt;
131  }
132  olddir = dir;
133  }
134  return tscatttot; // calculated both maximum and total scatter, total seems more useful, ignore max for now
135  }
136 
137  //-----------------------------------------------------------------------------------
138 
140  //float NumuSandFxs::getAvedEdx(const art::Ptr<rb::Track> track) {
141 
143 
144  //double dEdx;
145  float SumDeltaEn = 0.0;
146 
147  //Loop over muon track hits
148  for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
149  const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
150  if(bt->IsNoise(cellhit)) continue;
151  //pc.first = cellhit->Plane();
152  //pc.second = cellhit->Cell();
153 
154  //std::cout<<"_____________Hit is not noisy! "<<std::endl;
155 
156  rb::RecoHit rhit = track->RecoHit(cellhit);
157  if(!rhit.IsCalibrated()) continue;
158 
159  //std::cout<<"_____________Hit is calibrated! "<<std::endl;
160 
161  SumDeltaEn += rhit.GeV();
162  } //end loop over trk cell hits
163 
164  //std::cout<<"_____________SumDeltaEn: "<<SumDeltaEn<<std::endl;
165 
166  float DeltaX = track->TotalLength();
167  float dEdx = 0.0;
168  //some caveats
169  if(DeltaX > 30.0){
170  dEdx = 10e3*(SumDeltaEn)/DeltaX; // convert to MeV
171  }
172 
173  avededx = dEdx; // average dEdx in MeV/cm
174  //return dEdx; // average dEdx in GeV/cm -check units!
175  }
176 
177  //-----------------------------------------------------------------------------------
178 
180  //float NumuSandFxs::getAvedEdx(const art::Ptr<rb::Track> track) {
181 
182  if(track->NCell() <= 9){
183  avededxlast8 = 0.0; // average dEdx in MeV/cm
184  }
185  else{
186 
188 
189  std::vector<std::vector<double> > vecmuRecoHitInfo;
190 
191  double xx=0.0;
192  double yy=0.0;
193  double zz=0.0;
194  double ee=0.0;
195  double fpath = 0.0;
196 
197  std::pair<uint32_t, uint32_t> pc;
198  TVector3 point;
199  std::set<double> planeZBounds;
200  calib::FindZBoundaries(planeZBounds);
201  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track->Trajectory());
202  if(boundaryMap.size() < 1) return;
203 
204  //Loop over muon track hits
205  for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
206  const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
207  if(bt->IsNoise(cellhit)) continue;
208  pc.first = cellhit->Plane();
209  pc.second = cellhit->Cell();
210 
211  rb::RecoHit rhit = track->RecoHit(cellhit);
212  if(!rhit.IsCalibrated()) continue;
213 
214  xx=rhit.X();
215  yy=rhit.Y();
216  zz=rhit.Z();
217  ee=rhit.GeV();
218 
219  std::vector<double> rhitInfo(5);
220  rhitInfo[0]=xx;
221  rhitInfo[1]=yy;
222  rhitInfo[2]=zz;
223  rhitInfo[3]=ee;
224  point = TVector3(xx, yy, zz);
225 
226  fpath = calib::PathLengthInCell(boundaryMap,point,pc);
227  rhitInfo[4]=fpath;
228 
229  vecmuRecoHitInfo.push_back(rhitInfo);
230  } //end loop over trk cell hits
231 
232  // Another caveat
233  if(vecmuRecoHitInfo.size() <= 9){
234  avededxlast8 = 0.0; // average dEdx in MeV/cm
235  }
236  else{
237 
238  // using lambda function to sort it
239  std::sort(vecmuRecoHitInfo.begin(), vecmuRecoHitInfo.end(),
240  [](const std::vector<double>& a, const std::vector<double>& b) {
241  return a[2] < b[2];
242  });
243 
244  float sumEnLast8 = 0.0;
245  float sumPathLast8 =0.0;
246  for(size_t a = vecmuRecoHitInfo.size()-1; a > (vecmuRecoHitInfo.size()-9) ; --a) {
247  sumEnLast8 += 10e3*vecmuRecoHitInfo[a][3]; // In MeV
248  sumPathLast8 += vecmuRecoHitInfo[a][4]; // In cm
249  }
250  float dEdxLast8 = 0.0;
251  // some caveats
252  if(sumPathLast8 > 0.0){
253  dEdxLast8 = (sumEnLast8)/(sumPathLast8);
254  }
255 
256  avededxlast8 = dEdxLast8; // average dEdx in MeV/cm
257  }//end else
258  }//end else
259  }
260 
261  //-----------------------------------------------------------------------------------
262 
264 
265  if(track->NCell() <= 7){
266  avededxlast6 = 0.0; // average dEdx in MeV/cm
267  }
268  else{
269 
271 
272  std::vector<std::vector<double> > vecmuRecoHitInfo;
273 
274  double xx=0.0;
275  double yy=0.0;
276  double zz=0.0;
277  double ee=0.0;
278  double fpath = 0.0;
279 
280  std::pair<uint32_t, uint32_t> pc;
281  TVector3 point;
282  std::set<double> planeZBounds;
283  calib::FindZBoundaries(planeZBounds);
284  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track->Trajectory());
285  if(boundaryMap.size() < 1) return;
286 
287  //Loop over muon track hits
288  for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
289  const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
290  if(bt->IsNoise(cellhit)) continue;
291  pc.first = cellhit->Plane();
292  pc.second = cellhit->Cell();
293 
294  rb::RecoHit rhit = track->RecoHit(cellhit);
295  if(!rhit.IsCalibrated()) continue;
296 
297  xx=rhit.X();
298  yy=rhit.Y();
299  zz=rhit.Z();
300  ee=rhit.GeV();
301 
302  std::vector<double> rhitInfo(5);
303  rhitInfo[0]=xx;
304  rhitInfo[1]=yy;
305  rhitInfo[2]=zz;
306  rhitInfo[3]=ee;
307  point = TVector3(xx, yy, zz);
308 
309  fpath = calib::PathLengthInCell(boundaryMap,point,pc);
310  rhitInfo[4]=fpath;
311 
312  vecmuRecoHitInfo.push_back(rhitInfo);
313  } //end loop over trk cell hits
314 
315  // Another caveat
316  if(vecmuRecoHitInfo.size() <= 7){
317  avededxlast6 = 0.0; // average dEdx in MeV/cm
318  }
319  else{
320 
321  // using lambda function to sort it
322  std::sort(vecmuRecoHitInfo.begin(), vecmuRecoHitInfo.end(),
323  [](const std::vector<double>& a, const std::vector<double>& b) {
324  return a[2] < b[2];
325  });
326 
327  float sumEnLast6 = 0.0;
328  float sumPathLast6 =0.0;
329  for(size_t a = vecmuRecoHitInfo.size()-1; a > (vecmuRecoHitInfo.size()-7) ; --a) {
330  sumEnLast6 += 10e3*vecmuRecoHitInfo[a][3]; // In MeV
331  sumPathLast6 += vecmuRecoHitInfo[a][4]; // In cm
332  }
333  float dEdxLast6 = 0.0;
334  //some caveats
335  if(sumPathLast6 > 0.0){
336  dEdxLast6 = (sumEnLast6)/(sumPathLast6);
337  }
338 
339  avededxlast6 = dEdxLast6; // average dEdx in MeV/cm
340  }//end else
341  }//end else
342  }
343 
344  //-----------------------------------------------------------------------------------
345 
347 
348  //std::cout<<"_________track NCell:"<<track->NCell()<<std::endl;
349 
350  if(track->NCell() <= 5){
351  avededxlast4 = 0.0; // average dEdx in MeV/cm
352  }
353  else{
354 
356 
357  std::vector<std::vector<double> > vecmuRecoHitInfo;
358 
359  double xx=0.0;
360  double yy=0.0;
361  double zz=0.0;
362  double ee=0.0;
363  double fpath = 0.0;
364 
365  std::pair<uint32_t, uint32_t> pc;
366  TVector3 point;
367  std::set<double> planeZBounds;
368  calib::FindZBoundaries(planeZBounds);
369  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track->Trajectory());
370  if(boundaryMap.size() < 1) return;
371 
372  //Loop over muon track hits
373  for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
374  const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
375  if(bt->IsNoise(cellhit)) continue;
376  pc.first = cellhit->Plane();
377  pc.second = cellhit->Cell();
378 
379  rb::RecoHit rhit = track->RecoHit(cellhit);
380  if(!rhit.IsCalibrated()) continue;
381 
382  xx=rhit.X();
383  yy=rhit.Y();
384  zz=rhit.Z();
385  ee=rhit.GeV();
386 
387  std::vector<double> rhitInfo(5);
388  rhitInfo[0]=xx;
389  rhitInfo[1]=yy;
390  rhitInfo[2]=zz;
391  rhitInfo[3]=ee;
392  point = TVector3(xx, yy, zz);
393 
394  fpath = calib::PathLengthInCell(boundaryMap,point,pc);
395  rhitInfo[4]=fpath;
396 
397  //std::cout<<"rhitInfo[4]: "<<rhitInfo[4]<<", rhitInfo[3]: "<<rhitInfo[3]<<std::endl;
398 
399  vecmuRecoHitInfo.push_back(rhitInfo);
400  } //end loop over trk cell hits
401 
402  // Another caveat
403  if(vecmuRecoHitInfo.size() <= 5){
404  avededxlast4 = 0.0; // average dEdx in MeV/cm
405  }
406  else{
407 
408  //double vecsizeA = static_cast<double>(vecmuRecoHitInfo.size());
409  //std::cout<<"________vecsizeA: "<<vecmuRecoHitInfo.size()<<std::endl;
410 
411  // using lambda function to sort it
412  std::sort(vecmuRecoHitInfo.begin(), vecmuRecoHitInfo.end(),
413  [](const std::vector<double>& a, const std::vector<double>& b) {
414  return a[2] < b[2];
415  });
416 
417  float sumEnLast4 = 0.0;
418  float sumPathLast4 =0.0;
419 
420  //double vecsize = static_cast<double>(vecmuRecoHitInfo.size());
421  //std::cout<<"________vecsize: "<<vecsize<<std::endl;
422 
423  for(size_t a = vecmuRecoHitInfo.size()-1; a > (vecmuRecoHitInfo.size()-5) ; --a) {
424  //int aa = static_cast<int>(a);
425  //std::cout<<"_________a:"<< aa <<", --a: "<< --aa
426  //std::cout<<", vecmuRecoHitInfo size: "<< vecmuRecoHitInfo.size()-1
427  //<<", vecmuRecoHitInfo size -5: "<< vecmuRecoHitInfo.size()-5 <<std::endl;
428 
429  sumEnLast4 += 10e3*vecmuRecoHitInfo[a][3]; // In MeV
430  sumPathLast4 += vecmuRecoHitInfo[a][4]; // In cm
431  }
432  float dEdxLast4 = 0.0;
433  // some caveats
434  if(sumPathLast4 > 0.0){
435  dEdxLast4 = (sumEnLast4)/(sumPathLast4);
436  }
437 
438  avededxlast4 = dEdxLast4; // average dEdx in MeV/cm
439  }//end else
440  }//end else
441  }
442 
443  //-----------------------------------------------------------------------------------
444  //-----------------------------------------------------------------------------------
445  // Proxy for hadE[all hits not in main track] - calE[2nd reco track]
447  const art::PtrVector< rb::CellHit > track2Hits,
448  const art::PtrVector< rb::CellHit > sliceHits, float &actvtx) {
449 
450  // largely copied from NumuEAlg->MakeVertexCluster
451  rb::Cluster vertexCluster;
452  for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
453  bool match1 = 0;
454  bool match2 = 0;
455  for (unsigned int track1Hit = 0; track1Hit < track1Hits.size(); ++track1Hit){
456  match1 = ((sliceHits[sliceHit] == track1Hits[track1Hit]));
457  if (match1) break;
458  }//End of loop over hits in track 1
459  for (unsigned int track2Hit = 0; track2Hit < track2Hits.size(); ++track2Hit){
460  match2 = ((sliceHits[sliceHit] == track2Hits[track2Hit]));
461  if (match2) break;
462  }//End of loop over hits in track 2
463 
464  if (!(match1) && !(match2)) vertexCluster.Add(sliceHits[sliceHit]);
465  }//End of loop over hits in the slice
466 
467  //std::cout << "New activity clusters: " << trackHits.size() << " " << sliceHits.size() << " " << vertexCluster.NCell();
468  if(vertexCluster.NCell()>0){
469  actvtx = vertexCluster.TotalGeV();
470  //return vertexCluster.TotalGeV();
471 
472  }
473  else{
474  actvtx = 0;//return 0;
475 
476  }
477  }
478 
479  //-----------------------------------------------------------------------------------
480 
481  void NumuSandFxs::getActivity(const art::Ptr<rb::Track> track, double &startact, double &endact) {
482  double xyz[3];
483  TVector3 tstart = track->Start();
484  TVector3 tend = track->Stop();
487 
488  startact = 0;
489  endact = 0;
490 
491  for(unsigned int hitIdx = 0; hitIdx < track->NCell(); ++hitIdx){
492  art::Ptr<rb::CellHit> chit = track->Cell(hitIdx);
493  rb::RecoHit rhit = track->RecoHit(hitIdx);
494  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
495  if(chit->View()==geo::kX) { // ignore y position of cell center
496  if(sqrt((tstart.X()-xyz[0])*(tstart.X()-xyz[0]) + (tstart.Z()-xyz[2])*(tstart.Z()-xyz[2])) < 50) startact+=rhit.GeV();
497  if(sqrt((tend.X()-xyz[0])*(tend.X()-xyz[0]) + (tend.Z()-xyz[2])*(tend.Z()-xyz[2])) < 50) endact+=rhit.GeV();
498  }
499 
500  if(chit->View()==geo::kY) { // ignore x position of cell center
501  if(sqrt((tstart.Y()-xyz[1])*(tstart.Y()-xyz[1]) + (tstart.Z()-xyz[2])*(tstart.Z()-xyz[2])) < 50) startact+=rhit.GeV();
502  if(sqrt((tend.Y()-xyz[1])*(tend.Y()-xyz[1]) + (tend.Z()-xyz[2])*(tend.Z()-xyz[2])) < 50) endact+=rhit.GeV();
503  }
504  }
505  }
506 
507  //--------------------------------------------------------------------------- start get missing energy
508 
509  bool NumuSandFxs::getMissingE(const art::Ptr<rb::Cluster> slice, float &missE)
510  {
511 
513 
514  if(!bt->HaveTruthInfo()) return false; // No Truth Info!
515 
518 
519  art::PtrVector<rb::CellHit> allhits = slice->AllCells();
520 
521  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(allhits);
522  // const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(allhits);
523 
524  float sum=0;
525 
526  if(parts.size()==0) return false;
527 
528  for(unsigned int i = 0; i < parts.size(); i++) { // loop over particles
529  if(!(abs(parts[i]->PdgCode())==11 || abs(parts[i]->PdgCode())==13 || parts[i]->PdgCode()==2212 || abs(parts[i]->PdgCode())==211 || parts[i]->PdgCode()==22)) continue; // not interesting particle
530  if(parts[i]->E() < 0.001) continue; // ignore things < 1 MeV
531  if(geom->DetId() == novadaq::cnv::kFARDET) {
532  if(livegeom->IsPointLive(parts[i]->Position().Vect())==false || livegeom->IsPointLive(parts[i]->EndPosition().Vect())==false) sum+=parts[i]->E(); // not contained
533  }
534  else if(geom->DetId() == novadaq::cnv::kNEARDET) {
535  if( (livegeom->IsPointLive(parts[i]->Position().Vect())==false && livegeom->InMuonCatcher(parts[i]->Position().Vect())==false) ||
536  (livegeom->IsPointLive(parts[i]->EndPosition().Vect())==false && livegeom->InMuonCatcher(parts[i]->EndPosition().Vect())==false) ) sum+=parts[i]->E(); // not contained
537  }
538  }
539 
540  missE = sum;
541  return true;
542 
543  }
544 
545 } //end namespace
size_t NTrajectoryPoints() const
Definition: Track.h:83
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
Double_t xx
Definition: macro.C:12
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: Geo.cxx:373
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
Definition: CalibUtil.cxx:354
A collection of geometry functions used by the KalmanTrack tracking algorithm.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
void getAveTrackdEdx(const art::Ptr< rb::Track > track, float &avededx)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void getAveTrackdEdxLast4Cells(const art::Ptr< rb::Track > track, float &avededxlast4)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c)
Definition: NumuSandFxs.cxx:77
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
float getScatt(const art::Ptr< rb::Track > track)
void abs(TH1 *hist)
Double_t zz
Definition: plot.C:277
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
T acos(T number)
Definition: d0nt_math.hpp:54
double chi2()
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
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
trk::KalmanGeoHelper * fKGeo
Definition: NumuSandFxs.h:33
Far Detector at Ash River, MN.
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
void getActivityVtx(const art::PtrVector< rb::CellHit > track1Hits, const art::PtrVector< rb::CellHit > track2Hits, const art::PtrVector< rb::CellHit > sliceHits, float &actvtx)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
Float_t E
Definition: plot.C:20
olddir
Are we streaming the files via xrootd? txtcmd="cat %s | xargs -n1 samweb2xrootd > xrootd_inFile...
Definition: runNovaSAM.py:800
const double a
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
void getAveTrackdEdxLast8Cells(const art::Ptr< rb::Track > track, float &avededxlast8)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
#define PI
Definition: NumuSandFxs.cxx:19
TVector3 Unit() const
Sandbox to make objects for numu analysis in CAF.
Definition: FillSandboxes.h:7
size_type size() const
Definition: PtrVector.h:308
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float getAngle(TVector3 trackdir)
Definition: NumuSandFxs.cxx:39
float GeV() const
Definition: RecoHit.cxx:69
void getAveTrackdEdxLast6Cells(const art::Ptr< rb::Track > track, float &avededxlast6)
Float_t norm
TDirectory * dir
Definition: macro.C:5
float X() const
Definition: RecoHit.h:36
bool IsPointLive(TVector3 vertex)
Note the muon catcher is considered bad; use in combination with InMuonCatcher if needed...
float fitslope(const std::vector< double > &l, const std::vector< double > &t, const std::vector< double > &q, double &chisq, double &inter, bool linflag)
Definition: NumuSandFxs.cxx:51
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
const hit & b
Definition: hits.cxx:21
void getActivity(const art::Ptr< rb::Track > track, double &startact, double &endact)
assert(nhit_max >=nhit_nbins)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
bool getMissingE(const art::Ptr< rb::Cluster > slice, float &missE1)
float Y() const
Definition: RecoHit.h:37
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
bool InMuonCatcher(TVector3 vertex)
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
Encapsulate the geometry of one entire detector (near, far, ndos)
TVector3 Vect() const
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.
Definition: CalibUtil.cxx:498