Path.cxx
Go to the documentation of this file.
1 ///
2 /// \file Path.cxx
3 /// \brief Information about materials along a path in the detector
4 ///
6 #include <algorithm>
7 #include <cmath>
8 // ROOT includes
9 #include "TGeoMaterial.h"
10 // NOvA includes
14 using namespace bpfit;
15 
16 static const double gsMmuon = 105.658E-3;
17 static const double gsMpion = 139.570E-3;
18 static const double gsMproton = 938.272E-3;
19 
20 //......................................................................
21 
23  const double* p0,
24  const double* p1,
25  double s0) :
26  fGeo(geo),
27  fXYZ(2),
28  fS0(s0),
29  fTablesOK(false)
30 {
31  this->SetPoint(0,p0);
32  this->SetPoint(1,p1);
33 }
34 
35 //......................................................................
36 
38  unsigned int n,
39  double s0) :
40  fGeo(geo),
41  fXYZ(n),
42  fS0(s0),
43  fTablesOK(false)
44 { }
45 
46 //......................................................................
47 
49  std::vector<TVector3> const& trajectory,
50  double s0) :
51  fGeo(geo),
52  fXYZ(trajectory),
53  fS0(s0),
54  fTablesOK(false)
55 {
56 }
57 
58 //......................................................................
59 
60 void Path::SetPoint(unsigned int i, const double* xyz)
61 {
62  fXYZ[i].SetXYZ(xyz[0],xyz[1],xyz[2]);
63  fTablesOK = false;
64 }
65 
66 //......................................................................
67 
68 unsigned int Path::ScrubOne(double eps)
69 {
70  if (fXYZ.size()<3) return 0;
71 
72  std::vector<TVector3>::iterator p0(fXYZ.begin());
73  std::vector<TVector3>::iterator p1(p0+1);
74  std::vector<TVector3>::iterator pEnd(fXYZ.end());
75 
76  TVector3 d;
77  for (; p1!=pEnd; ++p0, ++p1) {
78  d.SetXYZ(p0->X()-p1->X(), p0->Y()-p1->Y(), p0->Z()-p1->Z());
79  if (d.Mag2()<eps) {
80  fXYZ.erase(p1);
81  fTablesOK = false;
82  return 1;
83  }
84  }
85  return 0;
86 }
87 
88 //......................................................................
89 
90 unsigned int Path::Scrub(double eps)
91 {
92  unsigned int n;
93  unsigned int ntot=0;
94  for (;;) {
95  n = this->ScrubOne(eps);
96  if (n==0) return ntot;
97  ntot += n;
98  }
99  return ntot;
100 }
101 
102 //......................................................................
103 
104 double Path::Length() const
105 {
106  if (!fTablesOK) this->Integrate();
107  return fS[fS.size()-1];
108 }
109 
110 //......................................................................
111 
112 double Path::X0(double s) const
113 {
114  if (!fTablesOK) this->Integrate();
115  return this->Interp(s-fS0, fS, fX0);
116 }
117 
118 //......................................................................
119 
120 TVector3 Path::Dir(double s)
121 {
122  if (!fTablesOK) this->Integrate();
123 
124  s -= fS0;
125  TVector3 d(0,0,1);
126  for (unsigned int i=0; i+1<fS.size(); ++i) {
127  if (fS[i]<=s && s<=fS[i+1]) {
128  d.SetXYZ(fXYZ[i+1].X()-fXYZ[i].X(),
129  fXYZ[i+1].Y()-fXYZ[i].Y(),
130  fXYZ[i+1].Z()-fXYZ[i].Z());
131  if(d.Mag() == 0.0) {
132  // Something terrible has happened if the direction vector is
133  // zero. Throw an exception to prevent EVERYTHING from crashing
134  // down...
135  throw BPException(__FILE__, __LINE__, kBAD_TRACK, 0);
136  }
137  d.SetMag(1.0);
138  return d;
139  }
140  }
141  return d;
142 }
143 
144 //......................................................................
145 
146 void Path::FindLayers(const double* x0,
147  const double* dx,
148  int c,
149  std::vector<double>& s,
150  std::vector<const TGeoMaterial*>& m) const
151 {
152  if (!(c==0||c==1)) abort();
153 
154  unsigned int i=0;
155  unsigned int j=0;
156  double x1[3];
157 
158  std::vector<double> ds;
159  std::vector<const TGeoMaterial*> mt;
160  bool done = false;
161  double rseed = 15.0; // Boundaries should be about this far away [cm]
162  double rmax = fGeo.DEdge(x0, dx); // ...and can't be further than this
163  double r = rseed;
164  for (done=false; !done && r<rmax; r*=2) {
165  x1[0] = x0[0] + r*dx[0];
166  x1[1] = x0[1] + r*dx[1];
167  x1[2] = x0[2] + r*dx[2];
168 
169  ds.clear();
170  mt.clear();
171  fGeo.MaterialsBetweenPoints(x0, x1, ds, mt);
172  if (mt.size()<2) continue;
173  for (i=1; i<mt.size(); ++i) {
174  //
175  // Case 0 searches for the first transition from either dead to
176  // active or active to dead. Appropriate for searching in the
177  // backwards track direction.
178  //
179  // Case 1 searches for the first dead to active
180  // transition. Appropriate for searching in the forward track
181  // direction
182  //
183  if (c==0 && (fGeo.IsActive(mt[i])!=fGeo.IsActive(mt[i-1]))) {
184  done = true;
185  break;
186  }
187  else if (c==1 && (fGeo.IsActive(mt[i])&&!fGeo.IsActive(mt[i-1]))) {
188  done = true;
189  break;
190  }
191  }
192  //
193  // A rare special case: If we stop in the middle of the last
194  // layer, we may not have counted its entire thickness. In that
195  // case, we aren't actually done and should expand and try again.
196  //
197  if ((i+1)==mt.size()) done=false;
198  }
199  for (j=0; j<mt.size(); ++j) {
200  if (j==i) break;
201  s.push_back(ds[j]);
202  m.push_back(mt[j]);
203  }
204 }
205 
206 //......................................................................
207 
209 {
210  //
211  // Nothing I can do with a path this short
212  //
213  if (fXYZ.size()<2) return;
214 
215  unsigned int i;
216  //
217  // Find the end point and directions from the track end point
218  //
219  unsigned int iend1 = fXYZ.size()-1;
220  unsigned int iend0 = iend1-1;
221  double x0[3];
222  double dxF[3];
223  double dxB[3];
224  double mag;
225  x0[0] = fXYZ[iend1].X();
226  x0[1] = fXYZ[iend1].Y();
227  x0[2] = fXYZ[iend1].Z();
228  dxF[0] = fXYZ[iend1].X()-fXYZ[iend0].X();
229  dxF[1] = fXYZ[iend1].Y()-fXYZ[iend0].Y();
230  dxF[2] = fXYZ[iend1].Z()-fXYZ[iend0].Z();
231  mag = sqrt(dxF[0]*dxF[0] + dxF[1]*dxF[1] + dxF[2]*dxF[2]);
232  if (mag>0.0) {
233  dxF[0] /= mag;
234  dxF[1] /= mag;
235  dxF[2] /= mag;
236  }
237  else {
238  return; // Something has gone wrong, bail on trying to make this adjustment.
239  }
240  dxB[0] = -dxF[0];
241  dxB[1] = -dxF[1];
242  dxB[2] = -dxF[2];
243  //
244  // Find the material layers between track end point and the most
245  // upstream and most down stream locations where the track could
246  // have ended.
247  //
248  int upstream = 0;
249  int dnstream = 1;
250  std::vector<double> sup, sdn;
251  std::vector<const TGeoMaterial*> mup, mdn;
252  this->FindLayers(x0, dxB, upstream, sup, mup);
253  this->FindLayers(x0, dxF, dnstream, sdn, mdn);
254  if (sup.size()==0 && sdn.size()==0) {
255  //
256  // Sometimes we don't find any adjacent layers inside the detector
257  // - for example, the track is calculated to have ended outside
258  // the detector. In this case, we cannot/should not make any
259  // adjustments to the track end point.
260  //
261  return;
262  }
263  //
264  // Build the table of stopping power vs. thickness from the list of
265  // materials
266  //
267  double s_sum = 0; // Integrated path length [cm]
268  double rho_sum = 0; // Integrated column density [g/cm^2]
269  std::vector<double> rho_i; // Integrated c.d. at layer i
270  std::vector<double> s_i; // Integrated path length at layer i
271  rho_i.push_back(rho_sum);
272  s_i. push_back(s_sum);
273  //
274  // First leg of the integral is backwards
275  //
276  for (i=0; i<sup.size(); ++i) {
277  s_sum -= sup[i];
278  rho_sum -= sup[i]*mup[i]->GetDensity()*mup[i]->GetZ()/mup[i]->GetA();
279  s_i. push_back(s_sum);
280  rho_i.push_back(rho_sum);
281  }
282  reverse(s_i. begin(), s_i. end());
283  reverse(rho_i.begin(), rho_i.end());
284  s_sum = 0;
285  rho_sum = 0;
286  for (i=0; i<sdn.size(); ++i) {
287  s_sum += sdn[i];
288  rho_sum += sdn[i]*mdn[i]->GetDensity()*mdn[i]->GetZ()/mdn[i]->GetA();
289  s_i. push_back(s_sum);
290  rho_i.push_back(rho_sum);
291  }
292  double rho_mid = (1-f)*rho_i[0] + f*rho_i[rho_i.size()-1];
293  double sadd = Path::Interp(rho_mid, rho_i, s_i, false);
294  if (sadd>0.0) {
295  fXYZ.push_back(TVector3(x0[0]+dxF[0]*sadd,
296  x0[1]+dxF[1]*sadd,
297  x0[2]+dxF[2]*sadd));
298  }
299  fTablesOK = false;
300 }
301 
302 //......................................................................
303 
305 {
306  //
307  // Nothing I can do with a path this short
308  //
309  if (fXYZ.size()<2) return;
310 
311  unsigned int i;
312  //
313  // Find the end point and directions from the track end point
314  //
315  unsigned int iend0 = 0;
316  unsigned int iend1 = 1;
317  double x0[3];
318  double dxF[3];
319  double dxB[3];
320  double mag;
321  x0[0] = fXYZ[iend0].X();
322  x0[1] = fXYZ[iend0].Y();
323  x0[2] = fXYZ[iend0].Z();
324  dxF[0] = fXYZ[iend1].X()-fXYZ[iend0].X();
325  dxF[1] = fXYZ[iend1].Y()-fXYZ[iend0].Y();
326  dxF[2] = fXYZ[iend1].Z()-fXYZ[iend0].Z();
327  mag = sqrt(dxF[0]*dxF[0] + dxF[1]*dxF[1] + dxF[2]*dxF[2]);
328  if (mag>0.0) {
329  dxF[0] /= mag;
330  dxF[1] /= mag;
331  dxF[2] /= mag;
332  }
333  else {
334  return; // Something has gone wrong, bail on trying to make this adjustment.
335  }
336  dxB[0] = -dxF[0];
337  dxB[1] = -dxF[1];
338  dxB[2] = -dxF[2];
339  //
340  // Find the material layers between track start point and the most
341  // upstream and most down stream locations where the track could
342  // have started.
343  //
344  int upstream = 1;
345  int dnstream = 0;
346  std::vector<double> sup, sdn;
347  std::vector<const TGeoMaterial*> mup, mdn;
348  this->FindLayers(x0, dxB, upstream, sup, mup);
349  this->FindLayers(x0, dxF, dnstream, sdn, mdn);
350  if (sup.size()==0 && sdn.size()==0) {
351  //
352  // Sometimes we don't find any adjacent layers inside the detector
353  // - for example, the track is calculated to have started outside
354  // the detector. In this case, we cannot/should not make any
355  // adjustments to the track end point.
356  //
357  return;
358  }
359  //
360  // Build the table of stopping power vs. thickness from the list of
361  // materials
362  //
363  double s_sum = 0; // Integrated path length [cm]
364  double rho_sum = 0; // Integrated column density [g/cm^2]
365  std::vector<double> rho_i; // Integrated c.d. at layer i
366  std::vector<double> s_i; // Integrated path length at layer i
367  rho_i.push_back(rho_sum);
368  s_i. push_back(s_sum);
369  //
370  // First leg of the integral is backwards
371  //
372  for (i=0; i<sup.size(); ++i) {
373  s_sum -= sup[i];
374  rho_sum -= sup[i]*mup[i]->GetDensity();
375  s_i. push_back(s_sum);
376  rho_i.push_back(rho_sum);
377  }
378  reverse(s_i. begin(), s_i. end());
379  reverse(rho_i.begin(), rho_i.end());
380  s_sum = 0;
381  rho_sum = 0;
382  for (i=0; i<sdn.size(); ++i) {
383  s_sum += sdn[i];
384  rho_sum += sdn[i]*mdn[i]->GetDensity();
385  s_i. push_back(s_sum);
386  rho_i.push_back(rho_sum);
387  }
388  double rho_mid = (1-f)*rho_i[0] + f*rho_i[rho_i.size()-1];
389  double sadd = Path::Interp(rho_mid, rho_i, s_i, false);
390  if (sadd<0.0) {
391  fXYZ.insert(fXYZ.begin(),TVector3(x0[0]+dxF[0]*sadd,
392  x0[1]+dxF[1]*sadd,
393  x0[2]+dxF[2]*sadd));
394  }
395  else {
396  fXYZ[0] = TVector3(x0[0]+dxF[0]*sadd,
397  x0[1]+dxF[1]*sadd,
398  x0[2]+dxF[2]*sadd);
399  }
400  fTablesOK = false;
401 }
402 
403 //......................................................................
404 ///
405 /// Given kinetic energy (GeV) and mass (GeV) compute other parameters
406 /// relavent for relativistic kinematics
407 ///
408 void Path::TPBG(double t,
409  double m,
410  double* T,
411  double* p,
412  double* beta,
413  double* gamma)
414 {
415  //
416  // Integration sometimes overshoots 0 -- clamp it here to 0.0
417  //
418  static double eps = 1.0E-5; // eps is allowable error
419  if (t<0.0) {
420  if (t>-eps) t = 0.0;
421  else
422  throw BPException(__FILE__, __LINE__, kBAD_KE, t);
423  }
424  double E = t+m;
425  double g = m/E;
426  double P = E*sqrt(1.0-g*g);
427 
428  if (T) *T = t;
429  if (p) *p = P;
430  if (beta) *beta = P/E;
431  if (gamma) *gamma = E/m;
432 }
433 
434 //......................................................................
435 
436 void Path::MuonParams(double s,
437  double* T,
438  double* p,
439  double* beta,
440  double* gamma) const
441 {
442  if (!fTablesOK) this->Integrate();
443  double t = this->Interp(s-fS0, fS, fMuonKE);
444  this->TPBG(t, gsMmuon, T, p, beta, gamma);
445 }
446 
447 //......................................................................
448 
449 void Path::PionParams(double s,
450  double* T,
451  double* p,
452  double* beta,
453  double* gamma) const
454 {
455  if (!fTablesOK) this->Integrate();
456  double t = this->Interp(s-fS0, fS, fPionKE);
457  this->TPBG(t, gsMpion, T, p, beta, gamma);
458 }
459 
460 //......................................................................
461 
462 void Path::ProtonParams(double s,
463  double* T,
464  double* p,
465  double* beta,
466  double* gamma) const
467 {
468  if (!fTablesOK) this->Integrate();
469  double t = this->Interp(s-fS0, fS, fProtonKE);
470  this->TPBG(t, gsMproton, T, p, beta, gamma);
471 }
472 
473 //......................................................................
474 ///
475 /// Integrate the dT = int[ dEdx(T), dx=0...ds ]
476 ///
477 /// \param ds - step size [cm]
478 /// \param mat - material name
479 /// \param m - mass of particle [GeV]
480 /// \param Tin - kinetic energy at start of step [GeV]
481 ///
482 double Path::IntegrateT(double ds,
483  const TGeoMaterial* mat,
484  double m,
485  double Tin) const
486 {
487  unsigned int i;
488  //
489  // Control the step size. Keep the steps no bigger than a fraction
490  // of a g/cm^2 in column density to avoid integration errors.
491  //
492  double mx = 0.1; // g/cm^2
493  unsigned int nstep = rint(ds*mat->GetDensity()/mx)+1;
494  double dx = ds/(float)nstep;
495 
496  // Cache these because otherwise we spend 40% of our time computing
497  // the logs of the same kinetic energies and dEdxs for every track.
498  static dEdxTable dedx_scintillator("Scintillator");
499  static dEdxTable dedx_pvc("PVC");
500  static dEdxTable dedx_glue("Glue");
501  static dEdxTable dedx_steel("Steel");
502  static dEdxTable dedx_vacuum("Vacuum");
503 
504  dEdxTable * dedx = &dedx_scintillator;
505  if (!strcmp(mat->GetName(), "Scintillator")) dedx = &dedx_scintillator;
506  else if(!strcmp(mat->GetName(), "PVC" )) dedx = &dedx_pvc;
507  else if(!strcmp(mat->GetName(), "Glue" )) dedx = &dedx_glue;
508  else if(!strcmp(mat->GetName(), "Steel" )) dedx = &dedx_steel;
509  else if(!strcmp(mat->GetName(), "Vacuum" )) dedx = &dedx_vacuum;
510  else{
511  // Create a dummy table to generate the warning, then throw it away
512  // and use vacuum just like the warnings says.
513  dEdxTable dummy(mat->GetName());
514  dedx = &dedx_vacuum;
515  }
516 
517  //
518  // Take the steps leap-frog style
519  //
520  double T = Tin;
521  double Tmid = 0.0;
522  for (i=0; i<nstep; ++i) {
523  Tmid = T + 0.5*dx*(*dedx)(T,m);
524  T += dx*(*dedx)(Tmid,m);
525  }
526  return T;
527 }
528 
529 //......................................................................
530 ///
531 /// Integrate along one straight-line segment of a path
532 ///
533 /// \param p1 - start x/y/z point [cm]
534 /// \param p2 - end x/y/z point [cm]
535 /// \param stot - path length, updated on return [cm]
536 /// \param x0 - path length in units of rad. len., updated on return.
537 /// \param Tmuon - muon KE, updated on return [GeV]
538 /// \param Tpion - pion KE, updated on return [GeV]
539 /// \param Tproton - proton KE, updated on return [GeV]
540 ///
541 void Path::IntegrateLeg(const double* p1,
542  const double* p2,
543  double* stot,
544  double* x0,
545  double* Tmuon,
546  double* Tpion,
547  double* Tproton) const
548 {
549  unsigned int i;
550  std::vector<double> ds;
551  std::vector<const TGeoMaterial*> mat;
552  fGeo.MaterialsBetweenPoints(p1, p2, ds, mat);
553 
554  for (i=0; i<ds.size(); ++i) {
555  *stot += ds[i];
556 
557  double rlen = mat[i]->GetRadLen();
558  //
559  // Something changed somewhere in ??? between 14-09-29 and
560  // 14-10-08 and now TGeoMaterial believes that the radiation
561  // length in Vacuum is zero. So for now I'm going to intercept
562  // and correct this number. I hope that a commit will follow this
563  // sometime later removing this hard coded case.
564  //
565  if(mat[i]->GetName() == std::string("Vacuum")) rlen = 1.0e30;
566  if (rlen <= 0.0) {
567  throw BPException(__FILE__, __LINE__, kBAD_RADL, rlen);
568  }
569 
570  *x0 += ds[i]/rlen;
571  *Tmuon = this->IntegrateT(ds[i], mat[i], gsMmuon, *Tmuon);
572  *Tpion = this->IntegrateT(ds[i], mat[i], gsMpion, *Tpion);
573  *Tproton = this->IntegrateT(ds[i], mat[i], gsMproton, *Tproton);
574 
575  fS. push_back(*stot);
576  fX0. push_back(*x0);
577  fMuonKE. push_back(*Tmuon);
578  fPionKE. push_back(*Tpion);
579  fProtonKE.push_back(*Tproton);
580  }
581 }
582 
583 //......................................................................
584 ///
585 /// Update all the look-up tables
586 ///
587 void Path::Integrate() const
588 {
589  fS. clear();
590  fX0. clear();
591  fMuonKE. clear();
592  fPionKE. clear();
593  fProtonKE.clear();
594  //
595  // Where to start the integration of the kinetic energies is a
596  // somewhat interesting question. A few relevant energy scales are:
597  //
598  // The dEdx tables are tabulated down to muon kinetic energies of
599  // 1e-5 GeV (beta*gamma = 0.014) so we really can't expect to start
600  // below that. That's about 0.01% of the particle mass.
601  //
602  // We expect ~2% resolution for muons of ~2 GeV ~= 40 MeV so it
603  // would make sense to start at a negligibly small fraction of that
604  // 40 MeV. That suggests something like T0 ~= (40/106)*1% ~= 0.4% of
605  // the particle mass.
606  //
607  // Splitting the difference gives 0.05% of particle mass.
608  //
609  double fT = 0.05E-2;
610  double Tmu = fT*gsMmuon;
611  double Tpi = fT*gsMpion;
612  double Tp = fT*gsMproton;
613  double stot = 0.0;
614  double x0 = 0.0;
615  fS. push_back(stot);
616  fX0. push_back(x0);
617  fMuonKE. push_back(Tmu);
618  fPionKE. push_back(Tpi);
619  fProtonKE.push_back(Tp);
620  //
621  // Walk the path backwards since we are integrating up the energy
622  // loss curve. Reverse the tables when done.
623  //
624  unsigned int i1 = fXYZ.size()-1;
625  unsigned int i2 = i1-1;
626  double p1[3], p2[3];
627 
628  for (; i1!=0; --i1, --i2) {
629  p1[0] = fXYZ[i1].X(); p1[1] = fXYZ[i1].Y(); p1[2] = fXYZ[i1].Z();
630  p2[0] = fXYZ[i2].X(); p2[1] = fXYZ[i2].Y(); p2[2] = fXYZ[i2].Z();
631 
632  this->IntegrateLeg(p1, p2, &stot, &x0, &Tmu, &Tpi, &Tp);
633  }
634  //
635  // Put all the tables in time order
636  //
637  unsigned int i;
638  std::reverse(fS.begin(), fS.end());
639  for (i=0; i<fS.size(); ++i) fS[i] = stot-fS[i];
640 
641  std::reverse(fX0.begin(), fX0.end());
642  for (i=0; i<fX0.size(); ++i) fX0[i] = x0-fX0[i];
643 
644  std::reverse(fMuonKE. begin(), fMuonKE. end());
645  std::reverse(fPionKE. begin(), fPionKE. end());
646  std::reverse(fProtonKE.begin(), fProtonKE.end());
647 
648  fTablesOK = true;
649 }
650 
651 //......................................................................
652 
653 double Path::Interp(double x,
654  const std::vector<double>& X,
655  const std::vector<double>& Y,
656  bool check_limits)
657 {
658  unsigned int i1=0, i2=1;
659  if (x< X[0]) { i1 = 0; i2 = 1; }
660  else if (x>=X[X.size()-1]) { i1 = X.size()-2; i2 = X.size()-1; }
661  else {
662  i2 = std::upper_bound(X.begin(), X.end(), x) - X.begin();
663  i1 = i2-1;
664  }
665  double rtn = Y[i1]+(Y[i2]-Y[i1])*(x-X[i1])/(X[i2]-X[i1]);
666  if (check_limits) {
667  //
668  // Trap illegal values here. Tables are all kinetic energies so
669  // nan's and negative values are illegal. Note that tracking errors
670  // of ~1 cell will result is uncertainties of up to about 12 MeV, so
671  // I've set the tolerance (eps) to 50 MeV to be safe.
672  //
673  static const double eps = 50.0E-3;
674  if (rtn<0.0) {
675  if (rtn>-eps) rtn = 0.0;
676  else throw BPException(__FILE__, __LINE__, kBAD_KE, rtn);
677  }
678  }
679  return rtn;
680 }
681 
682 ////////////////////////////////////////////////////////////////////////
double X0(double s) const
Return the number of radiation lengths crossed by the path.
Definition: Path.cxx:112
std::vector< double > fMuonKE
Muon kinetic energy.
Definition: Path.h:167
const unsigned int kBAD_KE
negative KE value < tol
Definition: BPException.h:18
double Length() const
Return the total path length.
Definition: Path.cxx:104
void FindLayers(const double *x0, const double *dx, int c, std::vector< double > &s, std::vector< const TGeoMaterial * > &m) const
Definition: Path.cxx:146
TVector3 Dir(double s=0.0)
Return the direction of the track after it has travelled a distance s.
Definition: Path.cxx:120
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
vector< vector< double > > clear
void AdjustEndPoint(double f)
Adjust the track end point to account for dead material beyond the end point.
Definition: Path.cxx:208
static const double gsMmuon
Definition: Path.cxx:16
void ProtonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:462
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
static double Interp(double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
Definition: Path.cxx:653
Double_t beta
static void TPBG(double t, double m, double *T, double *p, double *beta, double *gamma)
Definition: Path.cxx:408
Tables of mean energy loss in NOvA materials.
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
Float_t Y
Definition: plot.C:38
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
Float_t Z
Definition: plot.C:38
#define P(a, b, c, d, e, x)
const XML_Char * s
Definition: expat.h:262
double DEdge(const double *x0, const double *dx, double *exitp=0) const
std::string GetName(int i)
double dx[NP][NC]
void AdjustStartPoint(double f)
Definition: Path.cxx:304
void PionParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:449
Float_t E
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)
void IntegrateLeg(const double *p1, const double *p2, double *stot, double *x0, double *Tmuon, double *Tpion, double *Tproton) const
Definition: Path.cxx:541
base_types push_back(int_type())
Tables of mean energy loss in NOvA materials.
Definition: dEdxTable.h:13
Float_t d
Definition: plot.C:236
double fS0
s value at path start
Definition: Path.h:156
const double j
Definition: BetheBloch.cxx:29
A very simple service to remember what detector we&#39;re working in.
void MuonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Return the estimated muon parameters at path length s.
Definition: Path.cxx:436
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double gsMpion
Definition: Path.cxx:17
Path(const geo::GeometryBase &geo, const double *p0, const double *p1, double s0=0)
Construct a simple straight-line path from p0 to p1.
Definition: Path.cxx:22
std::vector< double > fX0
Path length / radiation length.
Definition: Path.h:166
TRandom3 r(0)
void Integrate() const
Definition: Path.cxx:587
void SetPoint(unsigned int i, const double *xyz)
Add a space point to the path.
Definition: Path.cxx:60
static const double gsMproton
Definition: Path.cxx:18
double T
Definition: Xdiff_gwt.C:5
unsigned int Scrub(double eps=0.01)
Remove points from the trajectory that have negligible distances between them.
Definition: Path.cxx:90
The geometry of one entire detector (near, far, ipnd)
Definition: GeometryBase.h:49
const unsigned int kBAD_RADL
non-positive radiation length
Definition: BPException.h:19
Float_t X
Definition: plot.C:38
Helper for AttenCurve.
Definition: Path.h:10
unsigned int ScrubOne(double eps)
Definition: Path.cxx:68
std::vector< double > fPionKE
Pion kinetic energy.
Definition: Path.h:168
const unsigned int kBAD_TRACK
a generic catch-all category for tracking failures
Definition: BPException.h:29
double IntegrateT(double ds, const TGeoMaterial *mat, double m, double Tin) const
Definition: Path.cxx:482
bool IsActive(const TGeoMaterial *m) const
Eigen::MatrixXd mat
std::vector< double > fProtonKE
Proton kinetic energy.
Definition: Path.h:169