Path.cxx

Go to the documentation of this file.

///
/// \file Path.cxx 
/// \brief Information about materials along a path in the detector
///
#include "BreakPointFitter/func/Path.h"
#include <algorithm>
#include <cmath>
// ROOT includes
#include "TGeoMaterial.h"
// NOvA includes
#include "GeometryObjects/GeometryBase.h"
#include "BreakPointFitter/func/dEdxTable.h"
#include "BreakPointFitter/func/BPException.h"
using namespace bpfit;

static const double gsMmuon   = 105.658E-3;
static const double gsMpion   = 139.570E-3;
static const double gsMproton = 938.272E-3;

//......................................................................

Path::Path(const geo::GeometryBase& geo,
           const double*            p0,
           const double*            p1,
           double                   s0) :
  fGeo(geo),
  fXYZ(2),
  fS0(s0),
  fTablesOK(false)
{
  this->SetPoint(0,p0);
  this->SetPoint(1,p1);
}

//......................................................................

Path::Path(const geo::GeometryBase& geo, 
           unsigned int               n,
           double                    s0) :
  fGeo(geo),
  fXYZ(n),
  fS0(s0),
  fTablesOK(false)
{ }

//......................................................................

Path::Path(const geo::GeometryBase& geo, 
           std::vector<TVector3> const& trajectory,
           double                    s0) :
  fGeo(geo),
  fXYZ(trajectory),
  fS0(s0),
  fTablesOK(false)
{ 
}

//......................................................................    

void Path::SetPoint(unsigned int i, const double* xyz) 
{
  fXYZ[i].SetXYZ(xyz[0],xyz[1],xyz[2]);
  fTablesOK = false;
}

//......................................................................

unsigned int Path::ScrubOne(double eps) 
{
  if (fXYZ.size()<3) return 0;

  std::vector<TVector3>::iterator p0(fXYZ.begin());
  std::vector<TVector3>::iterator p1(p0+1);
  std::vector<TVector3>::iterator pEnd(fXYZ.end());
  
  TVector3 d;
  for (; p1!=pEnd; ++p0, ++p1) {
    d.SetXYZ(p0->X()-p1->X(), p0->Y()-p1->Y(), p0->Z()-p1->Z());
    if (d.Mag2()<eps) {
      fXYZ.erase(p1);
      fTablesOK = false;
      return 1;
    }
  }
  return 0;
}

//......................................................................

unsigned int Path::Scrub(double eps) 
{
  unsigned int n;
  unsigned int ntot=0;
  for (;;) {
    n = this->ScrubOne(eps);
    if (n==0) return ntot;
    ntot += n;
  }
  return ntot;
}

//......................................................................

double Path::Length() const 
{
  if (!fTablesOK) this->Integrate();
  return fS[fS.size()-1];
}

//......................................................................

double Path::X0(double s) const 
{
  if (!fTablesOK) this->Integrate();
  return this->Interp(s-fS0, fS, fX0);
}

//......................................................................

TVector3 Path::Dir(double s)
{
  if (!fTablesOK) this->Integrate();
  
  s -= fS0;
  TVector3 d(0,0,1);
  for (unsigned int i=0; i+1<fS.size(); ++i) {
    if (fS[i]<=s && s<=fS[i+1]) {
      d.SetXYZ(fXYZ[i+1].X()-fXYZ[i].X(),
               fXYZ[i+1].Y()-fXYZ[i].Y(),
               fXYZ[i+1].Z()-fXYZ[i].Z());
      if(d.Mag() == 0.0) {
        // Something terrible has happened if the direction vector is
        // zero. Throw an exception to prevent EVERYTHING from crashing
        // down...
        throw BPException(__FILE__, __LINE__, kBAD_TRACK, 0);
      }
      d.SetMag(1.0);
      return d;
    }
  }
  return d;
}

//......................................................................

void Path::FindLayers(const double*                    x0,
                      const double*                    dx,
                      int                               c,
                      std::vector<double>&              s,
                      std::vector<const TGeoMaterial*>& m) const
{
  if (!(c==0||c==1)) abort();

  unsigned int i=0;
  unsigned int j=0;
  double       x1[3];

  std::vector<double>              ds;
  std::vector<const TGeoMaterial*> mt;
  bool   done  = false;
  double rseed = 15.0;   // Boundaries should be about this far away [cm]
  double rmax  = fGeo.DEdge(x0, dx); // ...and can't be further than this
  double r     = rseed;
  for (done=false; !done && r<rmax; r*=2) {
    x1[0] = x0[0] + r*dx[0];
    x1[1] = x0[1] + r*dx[1];
    x1[2] = x0[2] + r*dx[2];

    ds.clear();
    mt.clear();
    fGeo.MaterialsBetweenPoints(x0, x1, ds, mt);
    if (mt.size()<2) continue;
    for (i=1; i<mt.size(); ++i) {
      //
      // Case 0 searches for the first transition from either dead to
      // active or active to dead. Appropriate for searching in the
      // backwards track direction.
      //
      // Case 1 searches for the first dead to active
      // transition. Appropriate for searching in the forward track
      // direction
      //
      if (c==0 && (fGeo.IsActive(mt[i])!=fGeo.IsActive(mt[i-1]))) {
        done = true;
        break;
      }
      else if (c==1 && (fGeo.IsActive(mt[i])&&!fGeo.IsActive(mt[i-1]))) {
        done = true;
        break;
      }
    }
    //
    // A rare special case: If we stop in the middle of the last
    // layer, we may not have counted its entire thickness. In that
    // case, we aren't actually done and should expand and try again.
    //
    if ((i+1)==mt.size()) done=false;
  }
  for (j=0; j<mt.size(); ++j) {
    if (j==i) break;
    s.push_back(ds[j]);
    m.push_back(mt[j]);
  }
}

//......................................................................

void Path::AdjustEndPoint(double f)
{
  //
  // Nothing I can do with a path this short
  //
  if (fXYZ.size()<2) return;
  
  unsigned int i;
  //
  // Find the end point and directions from the track end point
  //
  unsigned int iend1 = fXYZ.size()-1;
  unsigned int iend0 = iend1-1;
  double x0[3];
  double dxF[3];
  double dxB[3];
  double mag;
  x0[0]  = fXYZ[iend1].X();
  x0[1]  = fXYZ[iend1].Y();
  x0[2]  = fXYZ[iend1].Z();
  dxF[0] = fXYZ[iend1].X()-fXYZ[iend0].X();
  dxF[1] = fXYZ[iend1].Y()-fXYZ[iend0].Y();
  dxF[2] = fXYZ[iend1].Z()-fXYZ[iend0].Z();
  mag = sqrt(dxF[0]*dxF[0] + dxF[1]*dxF[1] + dxF[2]*dxF[2]);
  if (mag>0.0) {
    dxF[0] /= mag;
    dxF[1] /= mag;
    dxF[2] /= mag;
  }
  else {
    return; // Something has gone wrong, bail on trying to make this adjustment.
  }
  dxB[0] = -dxF[0];
  dxB[1] = -dxF[1];
  dxB[2] = -dxF[2];
  //
  // Find the material layers between track end point and the most
  // upstream and most down stream locations where the track could
  // have ended.
  //
  int upstream = 0;
  int dnstream = 1;
  std::vector<double>              sup, sdn;
  std::vector<const TGeoMaterial*> mup, mdn;
  this->FindLayers(x0, dxB, upstream, sup, mup);
  this->FindLayers(x0, dxF, dnstream, sdn, mdn);
  if (sup.size()==0 && sdn.size()==0) {
    //
    // Sometimes we don't find any adjacent layers inside the detector
    // - for example, the track is calculated to have ended outside
    // the detector. In this case, we cannot/should not make any
    // adjustments to the track end point.
    //
    return;
  }
  //
  // Build the table of stopping power vs. thickness from the list of
  // materials
  //
  double s_sum   = 0;        // Integrated path length [cm]
  double rho_sum = 0;        // Integrated column density [g/cm^2]
  std::vector<double> rho_i; // Integrated c.d. at layer i
  std::vector<double> s_i;   // Integrated path length at layer i
  rho_i.push_back(rho_sum);
  s_i.  push_back(s_sum);
  //
  // First leg of the integral is backwards
  //
  for (i=0; i<sup.size(); ++i) {
    s_sum   -= sup[i];
    rho_sum -= sup[i]*mup[i]->GetDensity()*mup[i]->GetZ()/mup[i]->GetA();
    s_i.  push_back(s_sum);
    rho_i.push_back(rho_sum);
  }
  reverse(s_i.  begin(), s_i.  end());
  reverse(rho_i.begin(), rho_i.end());
  s_sum   = 0;
  rho_sum = 0;
  for (i=0; i<sdn.size(); ++i) {
    s_sum   += sdn[i];
    rho_sum += sdn[i]*mdn[i]->GetDensity()*mdn[i]->GetZ()/mdn[i]->GetA();
    s_i.  push_back(s_sum);
    rho_i.push_back(rho_sum);
  }
  double rho_mid = (1-f)*rho_i[0] + f*rho_i[rho_i.size()-1];
  double sadd    = Path::Interp(rho_mid, rho_i, s_i, false);
  if (sadd>0.0) {
    fXYZ.push_back(TVector3(x0[0]+dxF[0]*sadd,
                            x0[1]+dxF[1]*sadd,
                            x0[2]+dxF[2]*sadd));
  }
  fTablesOK = false;
}

//......................................................................

void Path::AdjustStartPoint(double f)
{
  //
  // Nothing I can do with a path this short
  //
  if (fXYZ.size()<2) return;

  unsigned int i;
  //
  // Find the end point and directions from the track end point
  //
  unsigned int iend0 = 0;
  unsigned int iend1 = 1;
  double x0[3];
  double dxF[3];
  double dxB[3];
  double mag;
  x0[0]  = fXYZ[iend0].X();
  x0[1]  = fXYZ[iend0].Y();
  x0[2]  = fXYZ[iend0].Z();
  dxF[0] = fXYZ[iend1].X()-fXYZ[iend0].X();
  dxF[1] = fXYZ[iend1].Y()-fXYZ[iend0].Y();
  dxF[2] = fXYZ[iend1].Z()-fXYZ[iend0].Z();
  mag = sqrt(dxF[0]*dxF[0] + dxF[1]*dxF[1] + dxF[2]*dxF[2]);
  if (mag>0.0) {
    dxF[0] /= mag;
    dxF[1] /= mag;
    dxF[2] /= mag;
  }
  else {
    return; // Something has gone wrong, bail on trying to make this adjustment.
  }
  dxB[0] = -dxF[0];
  dxB[1] = -dxF[1];
  dxB[2] = -dxF[2];
  //
  // Find the material layers between track start point and the most
  // upstream and most down stream locations where the track could
  // have started.
  //
  int upstream = 1;
  int dnstream = 0;
  std::vector<double>              sup, sdn;
  std::vector<const TGeoMaterial*> mup, mdn;
  this->FindLayers(x0, dxB, upstream, sup, mup);
  this->FindLayers(x0, dxF, dnstream, sdn, mdn);
  if (sup.size()==0 && sdn.size()==0) {
    //
    // Sometimes we don't find any adjacent layers inside the detector
    // - for example, the track is calculated to have started outside
    // the detector. In this case, we cannot/should not make any
    // adjustments to the track end point.
    //
    return;
  }
  //
  // Build the table of stopping power vs. thickness from the list of
  // materials
  //
  double s_sum   = 0;        // Integrated path length [cm]
  double rho_sum = 0;        // Integrated column density [g/cm^2]
  std::vector<double> rho_i; // Integrated c.d. at layer i
  std::vector<double> s_i;   // Integrated path length at layer i
  rho_i.push_back(rho_sum);
  s_i.  push_back(s_sum);
  //
  // First leg of the integral is backwards
  //
  for (i=0; i<sup.size(); ++i) {
    s_sum   -= sup[i];
    rho_sum -= sup[i]*mup[i]->GetDensity();
    s_i.  push_back(s_sum);
    rho_i.push_back(rho_sum);
  }
  reverse(s_i.  begin(), s_i.  end());
  reverse(rho_i.begin(), rho_i.end());
  s_sum   = 0;
  rho_sum = 0;
  for (i=0; i<sdn.size(); ++i) {
    s_sum   += sdn[i];
    rho_sum += sdn[i]*mdn[i]->GetDensity();
    s_i.  push_back(s_sum);
    rho_i.push_back(rho_sum);
  }
  double rho_mid = (1-f)*rho_i[0] + f*rho_i[rho_i.size()-1];
  double sadd    = Path::Interp(rho_mid, rho_i, s_i, false);
  if (sadd<0.0) {
    fXYZ.insert(fXYZ.begin(),TVector3(x0[0]+dxF[0]*sadd,
                                      x0[1]+dxF[1]*sadd,
                                      x0[2]+dxF[2]*sadd));
  }
  else {
    fXYZ[0] = TVector3(x0[0]+dxF[0]*sadd,
                       x0[1]+dxF[1]*sadd,
                       x0[2]+dxF[2]*sadd);
  }
  fTablesOK = false;
}

//......................................................................
///
/// Given kinetic energy (GeV) and mass (GeV) compute other parameters
/// relavent for relativistic kinematics
///
void Path::TPBG(double  t,
                double  m,
                double* T,
                double* p,
                double* beta,
                double* gamma)
{
  //
  // Integration sometimes overshoots 0 -- clamp it here to 0.0
  //
  static double eps = 1.0E-5; // eps is allowable error
  if (t<0.0) {
    if (t>-eps) t = 0.0;
    else
      throw BPException(__FILE__, __LINE__, kBAD_KE, t);
  }
  double E = t+m;
  double g = m/E;
  double P = E*sqrt(1.0-g*g);
  
  if (T)     *T = t;
  if (p)     *p = P;
  if (beta)  *beta  = P/E;
  if (gamma) *gamma = E/m;
}

//......................................................................

void Path::MuonParams(double s,
                      double* T, 
                      double* p,
                      double* beta,
                      double* gamma) const
{
  if (!fTablesOK) this->Integrate();
  double t = this->Interp(s-fS0, fS, fMuonKE);
  this->TPBG(t, gsMmuon, T, p, beta, gamma);
}

//......................................................................

void Path::PionParams(double s, 
                      double* T, 
                      double* p,
                      double* beta,
                      double* gamma) const
{
  if (!fTablesOK) this->Integrate();
  double t = this->Interp(s-fS0, fS, fPionKE);
  this->TPBG(t, gsMpion, T, p, beta, gamma);
}

//......................................................................

void Path::ProtonParams(double s, 
                        double* T, 
                        double* p,
                        double* beta,
                        double* gamma) const 
{
  if (!fTablesOK) this->Integrate();
  double t = this->Interp(s-fS0, fS, fProtonKE);
  this->TPBG(t, gsMproton, T, p, beta, gamma);
}

//......................................................................
///
/// Integrate the dT = int[ dEdx(T), dx=0...ds ]
///
/// \param ds  - step size [cm]
/// \param mat - material name
/// \param m   - mass of particle [GeV]
/// \param Tin - kinetic energy at start of step [GeV]
///
double Path::IntegrateT(double              ds, 
                        const TGeoMaterial* mat, 
                        double              m,
                        double              Tin) const
{
  unsigned int i;
  //
  // Control the step size. Keep the steps no bigger than a fraction
  // of a g/cm^2 in column density to avoid integration errors.
  //
  double          mx = 0.1; // g/cm^2
  unsigned int nstep = rint(ds*mat->GetDensity()/mx)+1;
  double          dx = ds/(float)nstep;

  // Cache these because otherwise we spend 40% of our time computing
  // the logs of the same kinetic energies and dEdxs for every track.
  static dEdxTable dedx_scintillator("Scintillator");
  static dEdxTable dedx_pvc("PVC");
  static dEdxTable dedx_glue("Glue");
  static dEdxTable dedx_steel("Steel");
  static dEdxTable dedx_vacuum("Vacuum");

  dEdxTable * dedx = &dedx_scintillator;
  if     (!strcmp(mat->GetName(), "Scintillator")) dedx = &dedx_scintillator;
  else if(!strcmp(mat->GetName(), "PVC"         )) dedx = &dedx_pvc;
  else if(!strcmp(mat->GetName(), "Glue"        )) dedx = &dedx_glue;
  else if(!strcmp(mat->GetName(), "Steel"       )) dedx = &dedx_steel;
  else if(!strcmp(mat->GetName(), "Vacuum"      )) dedx = &dedx_vacuum;
  else{
    // Create a dummy table to generate the warning, then throw it away
    // and use vacuum just like the warnings says.
    dEdxTable dummy(mat->GetName());
    dedx = &dedx_vacuum;
  }

  //
  // Take the steps leap-frog style
  //
  double T    = Tin;
  double Tmid = 0.0;
  for (i=0; i<nstep; ++i) {
    Tmid = T + 0.5*dx*(*dedx)(T,m);
    T += dx*(*dedx)(Tmid,m);
  }
  return T;
}

//......................................................................
///
/// Integrate along one straight-line segment of a path
///
/// \param p1      - start x/y/z point [cm]
/// \param p2      - end   x/y/z point [cm]
/// \param stot    - path length, updated on return [cm]
/// \param x0      - path length in units of rad. len., updated on return.
/// \param Tmuon   - muon KE, updated on return [GeV]
/// \param Tpion   - pion KE, updated on return [GeV]
/// \param Tproton - proton KE, updated on return [GeV]
///
void Path::IntegrateLeg(const double* p1,
                        const double* p2,
                        double* stot,
                        double* x0,
                        double* Tmuon,
                        double* Tpion,
                        double* Tproton) const
{
  unsigned int i;
  std::vector<double>              ds;
  std::vector<const TGeoMaterial*> mat;
  fGeo.MaterialsBetweenPoints(p1, p2, ds, mat);

  for (i=0; i<ds.size(); ++i) {
    *stot += ds[i];

    double rlen = mat[i]->GetRadLen();
    //
    // Something changed somewhere in ??? between 14-09-29 and
    // 14-10-08 and now TGeoMaterial believes that the radiation
    // length in Vacuum is zero.  So for now I'm going to intercept
    // and correct this number.  I hope that a commit will follow this
    // sometime later removing this hard coded case.
    //
    if(mat[i]->GetName() == std::string("Vacuum")) rlen = 1.0e30;
    if (rlen <= 0.0) {
      throw BPException(__FILE__, __LINE__, kBAD_RADL, rlen);
    }

    *x0 += ds[i]/rlen;
    *Tmuon   = this->IntegrateT(ds[i], mat[i], gsMmuon,   *Tmuon);
    *Tpion   = this->IntegrateT(ds[i], mat[i], gsMpion,   *Tpion);
    *Tproton = this->IntegrateT(ds[i], mat[i], gsMproton, *Tproton);

    fS.       push_back(*stot);
    fX0.      push_back(*x0);
    fMuonKE.  push_back(*Tmuon);
    fPionKE.  push_back(*Tpion);
    fProtonKE.push_back(*Tproton);
  }
}

//......................................................................
///
/// Update all the look-up tables
///
void Path::Integrate() const
{
  fS.       clear();
  fX0.      clear();
  fMuonKE.  clear();
  fPionKE.  clear();
  fProtonKE.clear();
  //
  // Where to start the integration of the kinetic energies is a
  // somewhat interesting question. A few relevant energy scales are:
  //
  // The dEdx tables are tabulated down to muon kinetic energies of
  // 1e-5 GeV (beta*gamma = 0.014) so we really can't expect to start
  // below that. That's about 0.01% of the particle mass.
  //
  // We expect ~2% resolution for muons of ~2 GeV ~= 40 MeV so it
  // would make sense to start at a negligibly small fraction of that
  // 40 MeV. That suggests something like T0 ~= (40/106)*1% ~= 0.4% of
  // the particle mass.
  //
  // Splitting the difference gives 0.05% of particle mass.
  //
  double fT   = 0.05E-2;
  double Tmu  = fT*gsMmuon;
  double Tpi  = fT*gsMpion;
  double Tp   = fT*gsMproton;
  double stot = 0.0;
  double x0   = 0.0;
  fS.       push_back(stot);
  fX0.      push_back(x0);
  fMuonKE.  push_back(Tmu);
  fPionKE.  push_back(Tpi);
  fProtonKE.push_back(Tp);
  //
  // Walk the path backwards since we are integrating up the energy
  // loss curve. Reverse the tables when done.
  //
  unsigned int i1 = fXYZ.size()-1;
  unsigned int i2 = i1-1;
  double p1[3], p2[3];

  for (; i1!=0; --i1, --i2) {
    p1[0] = fXYZ[i1].X(); p1[1] = fXYZ[i1].Y(); p1[2] = fXYZ[i1].Z();
    p2[0] = fXYZ[i2].X(); p2[1] = fXYZ[i2].Y(); p2[2] = fXYZ[i2].Z();
 
    this->IntegrateLeg(p1, p2, &stot, &x0, &Tmu, &Tpi, &Tp);
  }
  //
  // Put all the tables in time order
  //
  unsigned int i;
  std::reverse(fS.begin(), fS.end());
  for (i=0; i<fS.size(); ++i) fS[i] = stot-fS[i];

  std::reverse(fX0.begin(), fX0.end());
  for (i=0; i<fX0.size(); ++i) fX0[i] = x0-fX0[i];

  std::reverse(fMuonKE.  begin(), fMuonKE.  end());
  std::reverse(fPionKE.  begin(), fPionKE.  end());
  std::reverse(fProtonKE.begin(), fProtonKE.end());

  fTablesOK = true;
}

//......................................................................

double Path::Interp(double x, 
                    const std::vector<double>& X,
                    const std::vector<double>& Y,
                    bool check_limits)
{
  unsigned int i1=0, i2=1;
  if      (x< X[0])          { i1 = 0;          i2 = 1; }
  else if (x>=X[X.size()-1]) { i1 = X.size()-2; i2 = X.size()-1; }
  else {
    i2 = std::upper_bound(X.begin(), X.end(), x) - X.begin();
    i1 = i2-1;
  }
  double rtn = Y[i1]+(Y[i2]-Y[i1])*(x-X[i1])/(X[i2]-X[i1]);
  if (check_limits) {
    //
    // Trap illegal values here. Tables are all kinetic energies so
    // nan's and negative values are illegal. Note that tracking errors
    // of ~1 cell will result is uncertainties of up to about 12 MeV, so
    // I've set the tolerance (eps) to 50 MeV to be safe.
    //
    static const double eps = 50.0E-3;
    if (rtn<0.0) {
      if (rtn>-eps) rtn = 0.0;
      else throw BPException(__FILE__, __LINE__, kBAD_KE, rtn);
    }
  }
  return rtn;
}

////////////////////////////////////////////////////////////////////////