9 #include "TGeoMaterial.h" 14 using namespace bpfit;
49 std::vector<TVector3>
const& trajectory,
62 fXYZ[
i].SetXYZ(xyz[0],xyz[1],xyz[2]);
70 if (
fXYZ.size()<3)
return 0;
72 std::vector<TVector3>::iterator
p0(
fXYZ.begin());
73 std::vector<TVector3>::iterator
p1(
p0+1);
74 std::vector<TVector3>::iterator pEnd(
fXYZ.end());
77 for (; p1!=pEnd; ++
p0, ++
p1) {
78 d.SetXYZ(
p0->X()-p1->X(),
p0->Y()-p1->Y(),
p0->Z()-p1->Z());
96 if (n==0)
return ntot;
107 return fS[
fS.size()-1];
126 for (
unsigned int i=0;
i+1<
fS.size(); ++
i) {
127 if (
fS[
i]<=s && s<=
fS[
i+1]) {
149 std::vector<double>&
s,
150 std::vector<const TGeoMaterial*>&
m)
const 152 if (!(c==0||c==1)) abort();
158 std::vector<double>
ds;
159 std::vector<const TGeoMaterial*> mt;
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];
172 if (mt.size()<2)
continue;
173 for (i=1; i<mt.size(); ++
i) {
197 if ((i+1)==mt.size()) done=
false;
199 for (j=0; j<mt.size(); ++
j) {
213 if (
fXYZ.size()<2)
return;
219 unsigned int iend1 =
fXYZ.size()-1;
220 unsigned int iend0 = iend1-1;
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]);
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) {
269 std::vector<double> rho_i;
270 std::vector<double> s_i;
271 rho_i.push_back(rho_sum);
272 s_i. push_back(s_sum);
276 for (i=0; i<sup.size(); ++
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);
283 reverse(rho_i.begin(), rho_i.end());
286 for (i=0; i<sdn.size(); ++
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);
292 double rho_mid = (1-
f)*rho_i[0] + f*rho_i[rho_i.size()-1];
295 fXYZ.push_back(TVector3(x0[0]+dxF[0]*sadd,
309 if (
fXYZ.size()<2)
return;
315 unsigned int iend0 = 0;
316 unsigned int iend1 = 1;
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]);
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) {
365 std::vector<double> rho_i;
366 std::vector<double> s_i;
367 rho_i.push_back(rho_sum);
368 s_i. push_back(s_sum);
372 for (i=0; i<sup.size(); ++
i) {
374 rho_sum -= sup[
i]*mup[
i]->GetDensity();
375 s_i. push_back(s_sum);
376 rho_i.push_back(rho_sum);
379 reverse(rho_i.begin(), rho_i.end());
382 for (i=0; i<sdn.size(); ++
i) {
384 rho_sum += sdn[
i]*mdn[
i]->GetDensity();
385 s_i. push_back(s_sum);
386 rho_i.push_back(rho_sum);
388 double rho_mid = (1-
f)*rho_i[0] + f*rho_i[rho_i.size()-1];
391 fXYZ.insert(
fXYZ.begin(),TVector3(x0[0]+dxF[0]*sadd,
396 fXYZ[0] = TVector3(x0[0]+dxF[0]*sadd,
418 static double eps = 1.0E-5;
426 double P = E*
sqrt(1.0-g*g);
430 if (beta) *beta = P/
E;
431 if (gamma) *gamma = E/
m;
483 const TGeoMaterial*
mat,
493 unsigned int nstep = rint(ds*mat->GetDensity()/mx)+1;
498 static dEdxTable dedx_scintillator(
"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;
522 for (i=0; i<nstep; ++
i) {
523 Tmid = T + 0.5*dx*(*dedx)(
T,
m);
524 T += dx*(*dedx)(Tmid,
m);
547 double* Tproton)
const 550 std::vector<double>
ds;
551 std::vector<const TGeoMaterial*>
mat;
554 for (i=0; i<ds.size(); ++
i) {
557 double rlen = mat[
i]->GetRadLen();
575 fS. push_back(*stot);
624 unsigned int i1 =
fXYZ.size()-1;
625 unsigned int i2 = i1-1;
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();
632 this->
IntegrateLeg(p1, p2, &stot, &x0, &Tmu, &Tpi, &Tp);
639 for (i=0; i<
fS.size(); ++
i)
fS[i] = stot-
fS[i];
642 for (i=0; i<
fX0.size(); ++
i)
fX0[i] = x0-
fX0[i];
654 const std::vector<double>&
X,
655 const std::vector<double>&
Y,
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; }
662 i2 = std::upper_bound(X.begin(), X.end(),
x) - X.begin();
665 double rtn = Y[i1]+(Y[i2]-Y[i1])*(x-X[i1])/(X[i2]-X[i1]);
673 static const double eps = 50.0E-3;
675 if (rtn>-eps) rtn = 0.0;
double X0(double s) const
Return the number of radiation lengths crossed by the path.
std::vector< double > fMuonKE
Muon kinetic energy.
const unsigned int kBAD_KE
negative KE value < tol
double Length() const
Return the total path length.
void FindLayers(const double *x0, const double *dx, int c, std::vector< double > &s, std::vector< const TGeoMaterial * > &m) const
TVector3 Dir(double s=0.0)
Return the direction of the track after it has travelled a distance s.
Float_t x1[n_points_granero]
vector< vector< double > > clear
void AdjustEndPoint(double f)
Adjust the track end point to account for dead material beyond the end point.
static const double gsMmuon
void ProtonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
const geo::GeometryBase & fGeo
Pointer to the detector geometry.
static double Interp(double x, const std::vector< double > &X, const std::vector< double > &Y, bool check_limits=true)
static void TPBG(double t, double m, double *T, double *p, double *beta, double *gamma)
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?
std::vector< double > fS
Accumulated path length.
std::vector< TVector3 > fXYZ
x/y/z points along the path
#define P(a, b, c, d, e, x)
double DEdge(const double *x0, const double *dx, double *exitp=0) const
std::string GetName(int i)
void AdjustStartPoint(double f)
void PionParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
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
Tables of mean energy loss in NOvA materials.
double fS0
s value at path start
A very simple service to remember what detector we'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.
static const double gsMpion
Path(const geo::GeometryBase &geo, const double *p0, const double *p1, double s0=0)
Construct a simple straight-line path from p0 to p1.
std::vector< double > fX0
Path length / radiation length.
void SetPoint(unsigned int i, const double *xyz)
Add a space point to the path.
static const double gsMproton
unsigned int Scrub(double eps=0.01)
Remove points from the trajectory that have negligible distances between them.
The geometry of one entire detector (near, far, ipnd)
const unsigned int kBAD_RADL
non-positive radiation length
unsigned int ScrubOne(double eps)
std::vector< double > fPionKE
Pion kinetic energy.
const unsigned int kBAD_TRACK
a generic catch-all category for tracking failures
double IntegrateT(double ds, const TGeoMaterial *mat, double m, double Tin) const
bool IsActive(const TGeoMaterial *m) const
std::vector< double > fProtonKE
Proton kinetic energy.