Public Member Functions | Public Attributes | Private Member Functions | Static Private Member Functions | List of all members
bpfit::Path Class Reference

Track parameters (s, X0, KE, beta, ...) along a particular path in the detector. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-12-03/BreakPointFitter/func/Path.h"

Public Member Functions

 Path (const geo::GeometryBase &geo, const double *p0, const double *p1, double s0=0)
 Construct a simple straight-line path from p0 to p1. More...
 
 Path (const geo::GeometryBase &geo, unsigned int n, double s0=0)
 Construct a path which will have n points. More...
 
 Path (const geo::GeometryBase &geo, std::vector< TVector3 > const &trajectory, double s0=0)
 Construct a path using a vector of TVector3. More...
 
void SetPoint (unsigned int i, const double *xyz)
 Add a space point to the path. More...
 
unsigned int Scrub (double eps=0.01)
 Remove points from the trajectory that have negligible distances between them. More...
 
void AdjustEndPoint (double f)
 Adjust the track end point to account for dead material beyond the end point. More...
 
void AdjustStartPoint (double f)
 
double Length () const
 Return the total path length. More...
 
double X0 (double s) const
 Return the number of radiation lengths crossed by the path. More...
 
TVector3 Dir (double s=0.0)
 Return the direction of the track after it has travelled a distance s. More...
 
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. More...
 
void PionParams (double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
 
void ProtonParams (double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
 

Public Attributes

const geo::GeometryBasefGeo
 Pointer to the detector geometry. More...
 
std::vector< TVector3 > fXYZ
 x/y/z points along the path More...
 
double fS0
 s value at path start More...
 
bool fTablesOK
 Are the tables OK? More...
 
std::vector< double > fS
 Accumulated path length. More...
 
std::vector< double > fX0
 Path length / radiation length. More...
 
std::vector< double > fMuonKE
 Muon kinetic energy. More...
 
std::vector< double > fPionKE
 Pion kinetic energy. More...
 
std::vector< double > fProtonKE
 Proton kinetic energy. More...
 

Private Member Functions

unsigned int ScrubOne (double eps)
 
double IntegrateT (double ds, const TGeoMaterial *mat, double m, double Tin) const
 
void IntegrateLeg (const double *p1, const double *p2, double *stot, double *x0, double *Tmuon, double *Tpion, double *Tproton) const
 
void Integrate () const
 
void FindLayers (const double *x0, const double *dx, int c, std::vector< double > &s, std::vector< const TGeoMaterial * > &m) const
 

Static Private Member Functions

static void TPBG (double t, double m, double *T, double *p, double *beta, double *gamma)
 
static double Interp (double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
 

Detailed Description

Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.

Definition at line 18 of file Path.h.

Constructor & Destructor Documentation

Path::Path ( const geo::GeometryBase geo,
const double *  p0,
const double *  p1,
double  s0 = 0 
)

Construct a simple straight-line path from p0 to p1.

Parameters
geo- pointer to the detector geometry
p0- x/y/z coordinates of start point
p1- x/y/z coordinates of end point
s0- value of s at location p0

Definition at line 22 of file Path.cxx.

References SetPoint().

25  :
26  fGeo(geo),
27  fXYZ(2),
28  fS0(s0),
29  fTablesOK(false)
30 {
31  this->SetPoint(0,p0);
32  this->SetPoint(1,p1);
33 }
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
double fS0
s value at path start
Definition: Path.h:156
void SetPoint(unsigned int i, const double *xyz)
Add a space point to the path.
Definition: Path.cxx:60
Path::Path ( const geo::GeometryBase geo,
unsigned int  n,
double  s0 = 0 
)

Construct a path which will have n points.

Parameters
geo- pointer to the detector geometry
n- number of (xyz) points along path

Definition at line 37 of file Path.cxx.

39  :
40  fGeo(geo),
41  fXYZ(n),
42  fS0(s0),
43  fTablesOK(false)
44 { }
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
double fS0
s value at path start
Definition: Path.h:156
Path::Path ( const geo::GeometryBase geo,
std::vector< TVector3 > const &  trajectory,
double  s0 = 0 
)

Construct a path using a vector of TVector3.

Parameters
geo- pointer to the detector geometry
n- number of (xyz) points along path

Definition at line 48 of file Path.cxx.

50  :
51  fGeo(geo),
52  fXYZ(trajectory),
53  fS0(s0),
54  fTablesOK(false)
55 {
56 }
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
double fS0
s value at path start
Definition: Path.h:156

Member Function Documentation

void Path::AdjustEndPoint ( double  f)

Adjust the track end point to account for dead material beyond the end point.

Parameters
f- fractionally, how far to move the end point into the adjacent material

Definition at line 208 of file Path.cxx.

References febshutoff_auto::end, MakeMiniprodValidationCuts::f, FindLayers(), fTablesOK, fXYZ, MECModelEnuComparisons::i, Interp(), push_back(), site_stats_from_log::reverse, and std::sqrt().

Referenced by bpfit::BreakPoint::FitTracks().

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 }
void FindLayers(const double *x0, const double *dx, int c, std::vector< double > &s, std::vector< const TGeoMaterial * > &m) const
Definition: Path.cxx:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
static double Interp(double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
Definition: Path.cxx:653
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
base_types push_back(int_type())
void Path::AdjustStartPoint ( double  f)

Definition at line 304 of file Path.cxx.

References febshutoff_auto::end, MakeMiniprodValidationCuts::f, FindLayers(), fTablesOK, fXYZ, MECModelEnuComparisons::i, Interp(), push_back(), site_stats_from_log::reverse, and std::sqrt().

Referenced by bpfit::BreakPoint::FitTracks().

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 }
void FindLayers(const double *x0, const double *dx, int c, std::vector< double > &s, std::vector< const TGeoMaterial * > &m) const
Definition: Path.cxx:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
static double Interp(double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
Definition: Path.cxx:653
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
base_types push_back(int_type())
TVector3 Path::Dir ( double  s = 0.0)

Return the direction of the track after it has travelled a distance s.

Parameters
s- distance (cm) along the path

Definition at line 120 of file Path.cxx.

References d, fS, fS0, fTablesOK, fXYZ, MECModelEnuComparisons::i, Integrate(), bpfit::kBAD_TRACK, X, Y, and Z.

Referenced by bpfit::BreakPoint::FitTracks().

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 }
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
const XML_Char * s
Definition: expat.h:262
Float_t d
Definition: plot.C:236
double fS0
s value at path start
Definition: Path.h:156
void Integrate() const
Definition: Path.cxx:587
Float_t X
Definition: plot.C:38
const unsigned int kBAD_TRACK
a generic catch-all category for tracking failures
Definition: BPException.h:29
void Path::FindLayers ( const double *  x0,
const double *  dx,
int  c,
std::vector< double > &  s,
std::vector< const TGeoMaterial * > &  m 
) const
private

Definition at line 146 of file Path.cxx.

References geo::GeometryBase::DEdge(), fGeo, MECModelEnuComparisons::i, geo::GeometryBase::IsActive(), calib::j, geo::GeometryBase::MaterialsBetweenPoints(), r(), and x1.

Referenced by AdjustEndPoint(), and AdjustStartPoint().

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 }
Float_t x1[n_points_granero]
Definition: compare.C:5
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
const XML_Char * s
Definition: expat.h:262
double DEdge(const double *x0, const double *dx, double *exitp=0) const
double dx[NP][NC]
const double j
Definition: BetheBloch.cxx:29
A very simple service to remember what detector we&#39;re working in.
TRandom3 r(0)
bool IsActive(const TGeoMaterial *m) const
void Path::Integrate ( ) const
private

Update all the look-up tables

Definition at line 587 of file Path.cxx.

References clear, febshutoff_auto::end, fMuonKE, fPionKE, fProtonKE, fS, fTablesOK, fX0, fXYZ, gsMmuon, gsMpion, gsMproton, MECModelEnuComparisons::i, IntegrateLeg(), plot_validation_datamc::p1, plot_validation_datamc::p2, push_back(), and site_stats_from_log::reverse.

Referenced by Dir(), Length(), MuonParams(), PionParams(), ProtonParams(), and X0().

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 }
std::vector< double > fMuonKE
Muon kinetic energy.
Definition: Path.h:167
vector< vector< double > > clear
static const double gsMmuon
Definition: Path.cxx:16
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
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())
static const double gsMpion
Definition: Path.cxx:17
std::vector< double > fX0
Path length / radiation length.
Definition: Path.h:166
static const double gsMproton
Definition: Path.cxx:18
std::vector< double > fPionKE
Pion kinetic energy.
Definition: Path.h:168
std::vector< double > fProtonKE
Proton kinetic energy.
Definition: Path.h:169
void Path::IntegrateLeg ( const double *  p1,
const double *  p2,
double *  stot,
double *  x0,
double *  Tmuon,
double *  Tpion,
double *  Tproton 
) const
private

Integrate along one straight-line segment of a path

Parameters
p1- start x/y/z point [cm]
p2- end x/y/z point [cm]
stot- path length, updated on return [cm]
x0- path length in units of rad. len., updated on return.
Tmuon- muon KE, updated on return [GeV]
Tpion- pion KE, updated on return [GeV]
Tproton- proton KE, updated on return [GeV]

Definition at line 541 of file Path.cxx.

References fGeo, fMuonKE, fPionKE, fProtonKE, fS, fX0, GetName(), gsMmuon, gsMpion, gsMproton, MECModelEnuComparisons::i, IntegrateT(), bpfit::kBAD_RADL, geo::GeometryBase::MaterialsBetweenPoints(), push_back(), and string.

Referenced by Integrate().

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 }
std::vector< double > fMuonKE
Muon kinetic energy.
Definition: Path.h:167
static const double gsMmuon
Definition: Path.cxx:16
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
Definition: Path.h:154
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
std::string GetName(int i)
base_types push_back(int_type())
A very simple service to remember what detector we&#39;re working in.
static const double gsMpion
Definition: Path.cxx:17
std::vector< double > fX0
Path length / radiation length.
Definition: Path.h:166
static const double gsMproton
Definition: Path.cxx:18
const unsigned int kBAD_RADL
non-positive radiation length
Definition: BPException.h:19
std::vector< double > fPionKE
Pion kinetic energy.
Definition: Path.h:168
double IntegrateT(double ds, const TGeoMaterial *mat, double m, double Tin) const
Definition: Path.cxx:482
Eigen::MatrixXd mat
std::vector< double > fProtonKE
Proton kinetic energy.
Definition: Path.h:169
enum BeamMode string
double Path::IntegrateT ( double  ds,
const TGeoMaterial *  mat,
double  m,
double  Tin 
) const
private

Integrate the dT = int[ dEdx(T), dx=0...ds ]

Parameters
ds- step size [cm]
mat- material name
m- mass of particle [GeV]
Tin- kinetic energy at start of step [GeV]

Definition at line 482 of file Path.cxx.

References APDHVSetting::dummy, dx, check_time_usage::float, MECModelEnuComparisons::i, m, and T.

Referenced by IntegrateLeg().

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 }
double dx[NP][NC]
Tables of mean energy loss in NOvA materials.
Definition: dEdxTable.h:13
A very simple service to remember what detector we&#39;re working in.
double T
Definition: Xdiff_gwt.C:5
Eigen::MatrixXd mat
double Path::Interp ( double  x,
const std::vector< double > &  X,
const std::vector< double > &  Y,
bool  check_limits = true 
)
staticprivate

Definition at line 653 of file Path.cxx.

References bpfit::kBAD_KE, and submit_syst::x.

Referenced by AdjustEndPoint(), AdjustStartPoint(), MuonParams(), PionParams(), ProtonParams(), and X0().

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 }
const unsigned int kBAD_KE
negative KE value < tol
Definition: BPException.h:18
Float_t Y
Definition: plot.C:38
Float_t X
Definition: plot.C:38
double Path::Length ( ) const

Return the total path length.

Definition at line 104 of file Path.cxx.

References fS, fTablesOK, and Integrate().

Referenced by bpfit::BreakPoint::FillTrackNt().

105 {
106  if (!fTablesOK) this->Integrate();
107  return fS[fS.size()-1];
108 }
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
void Integrate() const
Definition: Path.cxx:587
void Path::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.

Parameters
s- position on path in cm
T- on return, kinetic energy in GeV
p- on return, momentum in GeV/c
beta- on return, beta=v/c=p/E
gamma- on return, gamma=1/sqrt(1-b^2)=E/m

Definition at line 436 of file Path.cxx.

References fMuonKE, fS, fS0, fTablesOK, gsMmuon, Integrate(), Interp(), confusionMatrixTree::t, and TPBG().

Referenced by bpfit::BreakPoint::FillTrackNt(), bpfit::BreakPoint::FitTracks(), and bpfit::ScatteringSurfaces::GetStepData().

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 }
std::vector< double > fMuonKE
Muon kinetic energy.
Definition: Path.h:167
const char * p
Definition: xmltok.h:285
static const double gsMmuon
Definition: Path.cxx:16
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
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
const XML_Char * s
Definition: expat.h:262
double fS0
s value at path start
Definition: Path.h:156
void Integrate() const
Definition: Path.cxx:587
double T
Definition: Xdiff_gwt.C:5
void Path::PionParams ( double  s,
double *  T = 0,
double *  p = 0,
double *  beta = 0,
double *  gamma = 0 
) const

Definition at line 449 of file Path.cxx.

References fPionKE, fS, fS0, fTablesOK, gsMpion, Integrate(), Interp(), confusionMatrixTree::t, and TPBG().

Referenced by bpfit::BreakPoint::FillTrackNt(), bpfit::BreakPoint::FitTracks(), and bpfit::ScatteringSurfaces::GetStepData().

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 }
const char * p
Definition: xmltok.h:285
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
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
const XML_Char * s
Definition: expat.h:262
double fS0
s value at path start
Definition: Path.h:156
static const double gsMpion
Definition: Path.cxx:17
void Integrate() const
Definition: Path.cxx:587
double T
Definition: Xdiff_gwt.C:5
std::vector< double > fPionKE
Pion kinetic energy.
Definition: Path.h:168
void Path::ProtonParams ( double  s,
double *  T = 0,
double *  p = 0,
double *  beta = 0,
double *  gamma = 0 
) const

Definition at line 462 of file Path.cxx.

References fProtonKE, fS, fS0, fTablesOK, gsMproton, Integrate(), Interp(), confusionMatrixTree::t, and TPBG().

Referenced by bpfit::BreakPoint::FillTrackNt(), bpfit::BreakPoint::FitTracks(), and bpfit::ScatteringSurfaces::GetStepData().

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 }
const char * p
Definition: xmltok.h:285
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
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
const XML_Char * s
Definition: expat.h:262
double fS0
s value at path start
Definition: Path.h:156
void Integrate() const
Definition: Path.cxx:587
static const double gsMproton
Definition: Path.cxx:18
double T
Definition: Xdiff_gwt.C:5
std::vector< double > fProtonKE
Proton kinetic energy.
Definition: Path.h:169
unsigned int Path::Scrub ( double  eps = 0.01)

Remove points from the trajectory that have negligible distances between them.

Parameters
eps- tolerance in cm^2 for deciding if points are close
Returns
The number of points removed from the lists.

Scrub is guranteed to always leave the end points in the list.

Definition at line 90 of file Path.cxx.

References getGoodRuns4SAM::n, and ScrubOne().

Referenced by bpfit::BreakPoint::FitTracks().

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 }
unsigned int ScrubOne(double eps)
Definition: Path.cxx:68
unsigned int Path::ScrubOne ( double  eps)
private

Definition at line 68 of file Path.cxx.

References d, fTablesOK, fXYZ, PandAna.Demos.pi0_spectra::p0, and plot_validation_datamc::p1.

Referenced by Scrub().

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 }
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
Float_t d
Definition: plot.C:236
void Path::SetPoint ( unsigned int  i,
const double *  xyz 
)

Add a space point to the path.

Parameters
i- point index [0...n)
xyz- x/y/z location triplet

Definition at line 60 of file Path.cxx.

References fTablesOK, fXYZ, and MECModelEnuComparisons::i.

Referenced by bpfit::BreakPoint::FitTracks(), and Path().

61 {
62  fXYZ[i].SetXYZ(xyz[0],xyz[1],xyz[2]);
63  fTablesOK = false;
64 }
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
void Path::TPBG ( double  t,
double  m,
double *  T,
double *  p,
double *  beta,
double *  gamma 
)
staticprivate

Given kinetic energy (GeV) and mass (GeV) compute other parameters relavent for relativistic kinematics

Definition at line 408 of file Path.cxx.

References E, MECModelEnuComparisons::g, bpfit::kBAD_KE, m, P, std::sqrt(), and confusionMatrixTree::t.

Referenced by MuonParams(), PionParams(), and ProtonParams().

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 }
const unsigned int kBAD_KE
negative KE value < tol
Definition: BPException.h:18
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Double_t beta
#define P(a, b, c, d, e, x)
Float_t E
Definition: plot.C:20
double T
Definition: Xdiff_gwt.C:5
double Path::X0 ( double  s) const

Return the number of radiation lengths crossed by the path.

Parameters
s- position along the path (cm)

Definition at line 112 of file Path.cxx.

References fS, fS0, fTablesOK, fX0, Integrate(), and Interp().

Referenced by bpfit::ScatteringSurfaces::GetStepData().

113 {
114  if (!fTablesOK) this->Integrate();
115  return this->Interp(s-fS0, fS, fX0);
116 }
static double Interp(double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
Definition: Path.cxx:653
bool fTablesOK
Are the tables OK?
Definition: Path.h:164
std::vector< double > fS
Accumulated path length.
Definition: Path.h:165
const XML_Char * s
Definition: expat.h:262
double fS0
s value at path start
Definition: Path.h:156
std::vector< double > fX0
Path length / radiation length.
Definition: Path.h:166
void Integrate() const
Definition: Path.cxx:587

Member Data Documentation

const geo::GeometryBase& bpfit::Path::fGeo

Pointer to the detector geometry.

Definition at line 154 of file Path.h.

Referenced by FindLayers(), and IntegrateLeg().

std::vector<double> bpfit::Path::fMuonKE
mutable

Muon kinetic energy.

Definition at line 167 of file Path.h.

Referenced by Integrate(), IntegrateLeg(), and MuonParams().

std::vector<double> bpfit::Path::fPionKE
mutable

Pion kinetic energy.

Definition at line 168 of file Path.h.

Referenced by Integrate(), IntegrateLeg(), and PionParams().

std::vector<double> bpfit::Path::fProtonKE
mutable

Proton kinetic energy.

Definition at line 169 of file Path.h.

Referenced by Integrate(), IntegrateLeg(), and ProtonParams().

std::vector<double> bpfit::Path::fS
mutable

Accumulated path length.

Definition at line 165 of file Path.h.

Referenced by Dir(), Integrate(), IntegrateLeg(), Length(), MuonParams(), PionParams(), ProtonParams(), and X0().

double bpfit::Path::fS0

s value at path start

Definition at line 156 of file Path.h.

Referenced by Dir(), MuonParams(), PionParams(), ProtonParams(), and X0().

bool bpfit::Path::fTablesOK
mutable

Are the tables OK?

Definition at line 164 of file Path.h.

Referenced by AdjustEndPoint(), AdjustStartPoint(), Dir(), Integrate(), Length(), MuonParams(), PionParams(), ProtonParams(), ScrubOne(), SetPoint(), and X0().

std::vector<double> bpfit::Path::fX0
mutable

Path length / radiation length.

Definition at line 166 of file Path.h.

Referenced by Integrate(), IntegrateLeg(), and X0().

std::vector<TVector3> bpfit::Path::fXYZ

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