trkf::TrajectoryMCSFitter Class Reference

Class for Maximum Likelihood fit of Multiple Coulomb Scattering angles between segments within a Track or Trajectory. More...

#include <TrajectoryMCSFitter.h>

Classes
struct	Config
struct	ScanResult

Public Types
using	Parameters = fhicl::Table< Config >

Public Member Functions
	TrajectoryMCSFitter (int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStepCoarse, double pStep, double fineScanWindow, const std::array< double, 5 > &angResol, const std::array< double, 5 > &hlParams, double segLenTolerance, bool applySCEcorr)
	TrajectoryMCSFitter (const Parameters &p)
recob::MCSFitResult	fitMcs (const recob::TrackTrajectory &traj) const
recob::MCSFitResult	fitMcs (const recob::Track &track) const
recob::MCSFitResult	fitMcs (const recob::Trajectory &traj) const
recob::MCSFitResult	fitMcs (const recob::TrackTrajectory &traj, int pid) const
recob::MCSFitResult	fitMcs (const recob::Track &track, int pid) const
recob::MCSFitResult	fitMcs (const recob::Trajectory &traj, int pid) const
void	breakTrajInSegments (const recob::TrackTrajectory &traj, std::vector< size_t > &breakpoints, std::vector< float > &segradlengths, std::vector< float > &cumseglens) const
void	linearRegression (const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) const
double	mcsLikelihood (double p, double theta0x, std::vector< float > &dthetaij, std::vector< float > &seg_nradl, std::vector< float > &cumLen, bool fwd, int pid) const
const ScanResult	doLikelihoodScan (std::vector< float > &dtheta, std::vector< float > &seg_nradlengths, std::vector< float > &cumLen, bool fwdFit, int pid, float detAngResol) const
const ScanResult	doLikelihoodScan (std::vector< float > &dtheta, std::vector< float > &seg_nradlengths, std::vector< float > &cumLen, bool fwdFit, int pid, float pmin, float pmax, float pstep, float detAngResol) const
double	HighlandFirstTerm (const double p) const
double	DetectorAngularResolution (const double uz) const
double	mass (int pid) const
double	energyLossBetheBloch (const double mass, const double e2) const
double	energyLossLandau (const double mass2, const double E2, const double x) const
double	GetE (const double initial_E, const double length_travelled, const double mass) const
int	minNSegs () const
double	segLen () const
double	segLenTolerance () const

Private Attributes
int	pIdHyp_
int	minNSegs_
double	segLen_
int	minHitsPerSegment_
int	nElossSteps_
int	eLossMode_
double	pMin_
double	pMax_
double	pStepCoarse_
double	pStep_
double	fineScanWindow_
std::array< double, 5 >	angResol_
std::array< double, 5 >	hlParams_
double	segLenTolerance_
bool	applySCEcorr_

Detailed Description

Class for Maximum Likelihood fit of Multiple Coulomb Scattering angles between segments within a Track or Trajectory.

Class for Maximum Likelihood fit of Multiple Coulomb Scattering angles between segments within a Track or Trajectory.

Inputs are: a Track or Trajectory, and various fit parameters (pIdHypothesis, minNumSegments, segmentLength, pMin, pMax, pStep, angResol)

Outputs are: a recob::MCSFitResult, containing: resulting momentum, momentum uncertainty, and best likelihood value (both for fwd and bwd fit); vector of comulative segment (radiation) lengths, vector of scattering angles, and PID hypothesis used in the fit. Note that the comulative segment length is what is used to compute the energy loss, but the segment length is actually slightly different, so the output can be used to reproduce the original results but they will not be identical (but very close).

For configuration options see TrajectoryMCSFitter::Config

Author

G. Cerati (FNAL, MicroBooNE)

Date

2017

Version

1.0

Definition at line 36 of file TrajectoryMCSFitter.h.

Member Typedef Documentation

using trkf::TrajectoryMCSFitter::Parameters = fhicl::Table<Config>

Definition at line 119 of file TrajectoryMCSFitter.h.

Constructor & Destructor Documentation

trkf::TrajectoryMCSFitter::TrajectoryMCSFitter ( int  pIdHyp, int  minNSegs, double  segLen, int  minHitsPerSegment, int  nElossSteps, int  eLossMode, double  pMin, double  pMax, double  pStepCoarse, double  pStep, double  fineScanWindow, const std::array< double, 5 > &  angResol, const std::array< double, 5 > &  hlParams, double  segLenTolerance, bool  applySCEcorr  )

inline

Definition at line 121 of file TrajectoryMCSFitter.h.

                                                                                                                                                                                                                                                                                                                               {
      pIdHyp_ = pIdHyp;
      minNSegs_ = minNSegs;
      segLen_ = segLen;
      minHitsPerSegment_ = minHitsPerSegment;
      nElossSteps_ = nElossSteps;
      eLossMode_ = eLossMode;
      pMin_ = pMin;
      pMax_ = pMax;
      pStepCoarse_ = pStepCoarse;
      pStep_ = pStep;
      fineScanWindow_ = fineScanWindow;
      angResol_ = angResol;
      hlParams_ = hlParams;
      segLenTolerance_ = segLenTolerance;
      applySCEcorr_ = applySCEcorr;
    }

trkf::TrajectoryMCSFitter::TrajectoryMCSFitter ( const Parameters &  p )

inlineexplicit

Definition at line 138 of file TrajectoryMCSFitter.h.

      : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStepCoarse(),p().pStep(),p().fineScanWindow(),p().angResol(),p().hlParams(),p().segLenTolerance(),p().applySCEcorr()) {}

Member Function Documentation

void TrajectoryMCSFitter::breakTrajInSegments	(	const recob::TrackTrajectory &	traj,
		std::vector< size_t > &	breakpoints,
		std::vector< float > &	segradlengths,
		std::vector< float > &	cumseglens
	)		const

Definition at line 61 of file TrajectoryMCSFitter.cxx.

                                                                                                                                                                          {
  //
  art::ServiceHandle<geo::Geometry const> geom;
  auto const* _SCE = (applySCEcorr_ ? lar::providerFrom<spacecharge::SpaceChargeService>() : NULL);
  //
  const double trajlen = traj.Length();
  const double thisSegLen = (trajlen>(segLen_*minNSegs_) ? segLen_ : trajlen/double(minNSegs_) );
  // std::cout << "track with length=" << trajlen << " broken in nseg=" << std::max(minNSegs_,int(trajlen/segLen_)) << " of length=" << thisSegLen << " where segLen_=" << segLen_ << std::endl;
  //
  constexpr double lar_radl_inv = 1./14.0;
  cumseglens.push_back(0.);//first segment has zero cumulative length from previous segments
  double thislen = 0.;
  double totlen = 0.;
  auto nextValid=traj.FirstValidPoint();
  breakpoints.push_back(nextValid);
  auto pos0 = traj.LocationAtPoint(nextValid);
  if (applySCEcorr_) {
    const double Position[3] = {pos0.X(), pos0.Y(), pos0.Z()};
    geo::TPCID tpcid = geom->FindTPCAtPosition(Position);
    geo::Vector_t pos0_offset = geo::Vector_t(0., 0., 0.);
    if(tpcid.isValid) {
      geo::Point_t p0(pos0.X(), pos0.Y(), pos0.Z());
      pos0_offset = _SCE->GetCalPosOffsets(p0, tpcid.TPC);
    }
    pos0.SetX(pos0.X() - pos0_offset.X());
    pos0.SetY(pos0.Y() + pos0_offset.Y());
    pos0.SetZ(pos0.Z() + pos0_offset.Z());
  }
  auto dir0 = traj.DirectionAtPoint(nextValid);
  nextValid = traj.NextValidPoint(nextValid+1);
  int npoints = 0;
  while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
    if (npoints==0) dir0 = traj.DirectionAtPoint(nextValid);
    auto pos1 = traj.LocationAtPoint(nextValid);
    if (applySCEcorr_) {
      const double Position[3] = {pos1.X(), pos1.Y(), pos1.Z()};
      geo::TPCID tpcid = geom->FindTPCAtPosition(Position);
      geo::Vector_t pos1_offset = geo::Vector_t(0., 0., 0.);
      if(tpcid.isValid) {
        geo::Point_t p1(pos1.X(), pos1.Y(), pos1.Z());
        pos1_offset = _SCE->GetCalPosOffsets(p1, tpcid.TPC);
      }
      pos1.SetX(pos1.X() - pos1_offset.X());
      pos1.SetY(pos1.Y() + pos1_offset.Y());
      pos1.SetZ(pos1.Z() + pos1_offset.Z());
    }
    //increments along the initial direction of the segment
    auto step = (pos1-pos0).R();
    thislen += dir0.Dot(pos1-pos0);
    totlen += step;
    pos0=pos1;
    // //fixme: testing alternative approaches here
    // //test1: increments following scatters
    // auto step = (pos1-pos0).R();
    // thislen += step;
    // totlen += step;
    // pos0=pos1;
    // //test2: end-start distance along the initial direction of the segment
    // thislen = dir0.Dot(pos1-pos0);
    // totlen = (pos1-pos0).R();
    //
    npoints++;
    if (thislen>=(thisSegLen-segLenTolerance_)) {
      breakpoints.push_back(nextValid);
      if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
      else segradlengths.push_back(-999.);
      cumseglens.push_back(totlen);
      thislen = 0.;
      npoints = 0;
    }
    nextValid = traj.NextValidPoint(nextValid+1);
  }
  //then add last segment
  if (thislen>0.) {
    breakpoints.push_back(traj.LastValidPoint()+1);
    segradlengths.push_back(thislen*lar_radl_inv);
    cumseglens.push_back(cumseglens.back()+thislen);
  }
  return;
}

double trkf::TrajectoryMCSFitter::DetectorAngularResolution ( const double  uz ) const

inline

Definition at line 170 of file TrajectoryMCSFitter.h.

                                                                   {
      return angResol_[0]/(uz*uz) + angResol_[1]/uz + angResol_[2] + angResol_[3]*uz + angResol_[4]*uz*uz;
    }

const TrajectoryMCSFitter::ScanResult TrajectoryMCSFitter::doLikelihoodScan	(	std::vector< float > &	dtheta,
		std::vector< float > &	seg_nradlengths,
		std::vector< float > &	cumLen,
		bool	fwdFit,
		int	pid,
		float	detAngResol
	)		const

Definition at line 183 of file TrajectoryMCSFitter.cxx.

                                                                                                                           {

  //do a first, coarse scan
  const ScanResult& coarseRes = doLikelihoodScan(dtheta, seg_nradlengths, cumLen, fwdFit, pid, pMin_, pMax_, pStepCoarse_, detAngResol);

  float pmax = std::min(coarseRes.p+fineScanWindow_,pMax_);
  float pmin = std::max(coarseRes.p-fineScanWindow_,pMin_);
  if (coarseRes.pUnc < (std::numeric_limits<float>::max()-1.)) {
    pmax = std::min(coarseRes.p+2*coarseRes.pUnc,pMax_);
    pmin = std::max(coarseRes.p-2*coarseRes.pUnc,pMin_);
  }

  //do the fine grained scan in a smaller region
  const ScanResult& refineRes = doLikelihoodScan(dtheta, seg_nradlengths, cumLen, fwdFit, pid, pmin, pmax, pStep_, detAngResol);

  return refineRes;
}

const TrajectoryMCSFitter::ScanResult TrajectoryMCSFitter::doLikelihoodScan	(	std::vector< float > &	dtheta,
		std::vector< float > &	seg_nradlengths,
		std::vector< float > &	cumLen,
		bool	fwdFit,
		int	pid,
		float	pmin,
		float	pmax,
		float	pstep,
		float	detAngResol
	)		const

Definition at line 142 of file TrajectoryMCSFitter.cxx.

                                                                                                                                                                {
  int    best_idx  = -1;
  float best_logL = std::numeric_limits<float>::max();
  float best_p    = -1.0;
  std::vector<float> vlogL;
  for (float p_test = pmin; p_test <= pmax; p_test+=pstep) {
    float logL = mcsLikelihood(p_test, detAngResol, dtheta, seg_nradlengths, cumLen, fwdFit, pid);
    if (logL < best_logL) {
      best_p    = p_test;
      best_logL = logL;
      best_idx  = vlogL.size();
    }
    vlogL.push_back(logL);
  }
  //
  //uncertainty from left side scan
  float lunc = -1.;
  if (best_idx>0) {
    for (int j=best_idx-1;j>=0;j--) {
      float dLL = vlogL[j]-vlogL[best_idx];
      if ( dLL>=0.5 ) {
        lunc = (best_idx-j)*pstep;
        break;
      }
    }
  }
  //uncertainty from right side scan
  float runc = -1.;
  if (best_idx<int(vlogL.size()-1)) {
    for (unsigned int j=best_idx+1;j<vlogL.size();j++) {
      float dLL = vlogL[j]-vlogL[best_idx];
      if ( dLL>=0.5 ) {
        runc = (j-best_idx)*pstep;
        break;
      }
    }
  }
  return ScanResult(best_p, std::max(lunc,runc), best_logL);
}

double TrajectoryMCSFitter::energyLossBetheBloch	(	const double	mass,
		const double	e2
	)		const

Definition at line 346 of file TrajectoryMCSFitter.cxx.

                                                                                        {
  // stolen, mostly, from GFMaterialEffects.
  constexpr double Iinv = 1./188.E-6;
  constexpr double matConst = 1.4*18./40.;//density*Z/A
  constexpr double me = 0.511;
  constexpr double kappa = 0.307075;
  //
  const double beta2 = (e2-mass*mass)/e2;
  const double gamma2 = 1./(1.0 - beta2);
  const double massRatio = me/mass;
  const double argument = (2.*me*gamma2*beta2*Iinv) * std::sqrt(1+2*std::sqrt(gamma2)*massRatio + massRatio*massRatio);
  //
  double dedx = kappa*matConst/beta2;
  //
  if (mass==0.0) return(0.0);
  if (argument <= exp(beta2)) {
      dedx = 0.;
  } else{
    dedx *= (log(argument)-beta2)*1.E-3; // Bethe-Bloch, converted to GeV/cm
    if (dedx<0.) dedx = 0.;
  }
  return dedx;
}

double TrajectoryMCSFitter::energyLossLandau	(	const double	mass2,
		const double	E2,
		const double	x
	)		const

Definition at line 328 of file TrajectoryMCSFitter.cxx.

                                                                                                     {
  //
  // eq. (33.11) in http://pdg.lbl.gov/2016/reviews/rpp2016-rev-passage-particles-matter.pdf (except density correction is ignored)
  //
  if (x<=0.) return 0.;
  constexpr double Iinv2 = 1./(188.E-6*188.E-6);
  constexpr double matConst = 1.4*18./40.;//density*Z/A
  constexpr double me = 0.511;
  constexpr double kappa = 0.307075;
  constexpr double j = 0.200;
  //
  const double beta2 = (e2-mass2)/e2;
  const double gamma2 = 1./(1.0 - beta2);
  const double epsilon = 0.5*kappa*x*matConst/beta2;
  //
  return 0.001*epsilon*( log(2.*me*beta2*gamma2*epsilon*Iinv2) + j - beta2 );
}

recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::TrackTrajectory &  traj ) const

inline

Definition at line 141 of file TrajectoryMCSFitter.h.

{ return fitMcs(traj,pIdHyp_); }

recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Track &  track ) const

inline

Definition at line 142 of file TrajectoryMCSFitter.h.

{ return fitMcs(track,pIdHyp_); }

recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory &  traj ) const

inline

Definition at line 143 of file TrajectoryMCSFitter.h.

{ return fitMcs(traj,pIdHyp_); }

recob::MCSFitResult TrajectoryMCSFitter::fitMcs	(	const recob::TrackTrajectory &	traj,
		int	pid
	)		const

Definition at line 12 of file TrajectoryMCSFitter.cxx.

                                                                                             {
  //
  // Break the trajectory in segments of length approximately equal to segLen_
  //
  vector<size_t> breakpoints;
  vector<float> segradlengths;
  vector<float> cumseglens;
  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
  //
  // Fit segment directions, and get 3D angles between them
  //
  if (segradlengths.size()<2) return recob::MCSFitResult();
  vector<float> dtheta;
  Vector_t pcdir0;
  Vector_t pcdir1;
  for (unsigned int p = 0; p<segradlengths.size(); p++) {
    linearRegression(traj, breakpoints[p], breakpoints[p+1], pcdir1);
    if (p>0) {
      if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
        dtheta.push_back(-999.);
      } else {
        const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
        //assert(std::abs(cosval)<=1);
        //units are mrad
        double dt = 1000.*acos(cosval);//should we try to use expansion for small angles?
        dtheta.push_back(dt);
      }
    }
    pcdir0 = pcdir1;
  }
  //
  // Perform likelihood scan in forward and backward directions
  //
  vector<float> cumLenFwd;
  vector<float> cumLenBwd;
  for (unsigned int i = 0; i<cumseglens.size()-2; i++) {
    cumLenFwd.push_back(cumseglens[i]);
    cumLenBwd.push_back(cumseglens.back()-cumseglens[i+2]);
  }
  double detAngResol = DetectorAngularResolution(std::abs(traj.StartDirection().Z()));
  const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true , pid, detAngResol);
  const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, pid, detAngResol);
  //
  return recob::MCSFitResult(pid,
                             fwdResult.p,fwdResult.pUnc,fwdResult.logL,
                             bwdResult.p,bwdResult.pUnc,bwdResult.logL,
                             segradlengths,dtheta);
}

recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Track &  track, int  pid  ) const

inline

Definition at line 146 of file TrajectoryMCSFitter.h.

{ return fitMcs(track.Trajectory(),pid); }

recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory &  traj, int  pid  ) const

inline

Definition at line 147 of file TrajectoryMCSFitter.h.

                                                                              {
      recob::TrackTrajectory::Flags_t flags(traj.NPoints());
      const recob::TrackTrajectory tt(traj,std::move(flags));
      return fitMcs(tt,pid);
    }

double TrajectoryMCSFitter::GetE	(	const double	initial_E,
		const double	length_travelled,
		const double	mass
	)		const

Definition at line 370 of file TrajectoryMCSFitter.cxx.

                                                                                                            {
  //
  if (eLossMode_==1) {
    // ELoss mode: MIP (constant)
    constexpr double kcal = 0.002105;
    return (initial_E - kcal*length_travelled);//energy at this segment
  }
  //
  // Non constant energy loss distribution
  const double step_size = length_travelled / nElossSteps_;
  //
  double current_E = initial_E;
  const double m2 = m*m;
  //
  for (auto i = 0; i < nElossSteps_; ++i) {
    if (eLossMode_==2) {
      double dedx = energyLossBetheBloch(m,current_E);
      current_E -= (dedx * step_size);
    } else {
      // MPV of Landau energy loss distribution
      current_E -= energyLossLandau(m2,current_E*current_E,step_size);
    }
    if ( current_E <= m ) {
      // std::cout<<"WARNING: current_E less than mu mass. it is "<<current_E<<std::endl;
      return 0.;
    }
  }
  return current_E;
}

double trkf::TrajectoryMCSFitter::HighlandFirstTerm ( const double  p ) const

inline

Definition at line 167 of file TrajectoryMCSFitter.h.

                                                          {
      return hlParams_[0]/(p*p) + hlParams_[1]/p + hlParams_[2] + hlParams_[3]*p + hlParams_[4]*p*p;
    }

void TrajectoryMCSFitter::linearRegression	(	const recob::TrackTrajectory &	traj,
		const size_t	firstPoint,
		const size_t	lastPoint,
		recob::tracking::Vector_t &	pcdir
	)		const

Definition at line 202 of file TrajectoryMCSFitter.cxx.

                                                                                                                                                   {
  //
  art::ServiceHandle<geo::Geometry const> geom;
  auto const* _SCE = (applySCEcorr_ ? lar::providerFrom<spacecharge::SpaceChargeService>() : NULL);
  //
  int npoints = 0;
  geo::vect::MiddlePointAccumulator middlePointCalc;
  size_t nextValid = firstPoint;
  //fixme explore a max number of points to use for linear regression
  //while (nextValid<std::min(firstPoint+10,lastPoint)) {
  while (nextValid<lastPoint) {
    auto tempP = traj.LocationAtPoint(nextValid);
    if (applySCEcorr_) {
      const double Position[3] = {tempP.X(), tempP.Y(), tempP.Z()};
      geo::TPCID tpcid = geom->FindTPCAtPosition(Position);
      geo::Vector_t tempP_offset = geo::Vector_t(0., 0., 0.);
      if(tpcid.isValid) {
        geo::Point_t ptemp(tempP.X(), tempP.Y(), tempP.Z());
        tempP_offset = _SCE->GetCalPosOffsets(ptemp, tpcid.TPC);
      }
      tempP.SetX(tempP.X() - tempP_offset.X());
      tempP.SetY(tempP.Y() + tempP_offset.Y());
      tempP.SetZ(tempP.Z() + tempP_offset.Z());
    }
    middlePointCalc.add(tempP);
    //middlePointCalc.add(traj.LocationAtPoint(nextValid));
    nextValid = traj.NextValidPoint(nextValid+1);
    npoints++;
  }
  const auto avgpos = middlePointCalc.middlePoint();
  const double norm = 1./double(npoints);
  //
  //assert(npoints>0);
  //
  TMatrixDSym m(3);
  nextValid = firstPoint;
  while (nextValid<lastPoint) {
    auto p = traj.LocationAtPoint(nextValid);
    if (applySCEcorr_) {
      const double Position[3] = {p.X(), p.Y(), p.Z()};
      geo::TPCID tpcid = geom->FindTPCAtPosition(Position);
      geo::Vector_t p_offset = geo::Vector_t(0., 0., 0.);
      if(tpcid.isValid) {
        geo::Point_t point(p.X(), p.Y(), p.Z());
        p_offset = _SCE->GetCalPosOffsets(point, tpcid.TPC);
      }
      p.SetX(p.X() - p_offset.X());
      p.SetY(p.Y() + p_offset.Y());
      p.SetZ(p.Z() + p_offset.Z());
    }
    const double xxw0 = p.X()-avgpos.X();
    const double yyw0 = p.Y()-avgpos.Y();
    const double zzw0 = p.Z()-avgpos.Z();
    m(0, 0) += xxw0*xxw0*norm;
    m(0, 1) += xxw0*yyw0*norm;
    m(0, 2) += xxw0*zzw0*norm;
    m(1, 0) += yyw0*xxw0*norm;
    m(1, 1) += yyw0*yyw0*norm;
    m(1, 2) += yyw0*zzw0*norm;
    m(2, 0) += zzw0*xxw0*norm;
    m(2, 1) += zzw0*yyw0*norm;
    m(2, 2) += zzw0*zzw0*norm;
    nextValid = traj.NextValidPoint(nextValid+1);
  }
  //
  const TMatrixDSymEigen me(m);
  const auto& eigenval = me.GetEigenValues();
  const auto& eigenvec = me.GetEigenVectors();
  //
  int maxevalidx = 0;
  double maxeval = eigenval(0);
  for (int i=1; i<3; ++i) {
    if (eigenval(i)>maxeval) {
      maxevalidx = i;
      maxeval = eigenval(i);
    }
  }
  //
  pcdir = Vector_t(eigenvec(0, maxevalidx),eigenvec(1, maxevalidx),eigenvec(2, maxevalidx));
  if (traj.DirectionAtPoint(firstPoint).Dot(pcdir)<0.) pcdir*=-1.;
  //
}

double trkf::TrajectoryMCSFitter::mass ( int  pid ) const

inline

Definition at line 173 of file TrajectoryMCSFitter.h.

                               {
      if (abs(pid)==13)   { return mumass; }
      if (abs(pid)==211)  { return pimass; }
      if (abs(pid)==321)  { return kmass;  }
      if (abs(pid)==2212) { return pmass;  }
      return util::kBogusD;
    }

double TrajectoryMCSFitter::mcsLikelihood	(	double	p,
		double	theta0x,
		std::vector< float > &	dthetaij,
		std::vector< float > &	seg_nradl,
		std::vector< float > &	cumLen,
		bool	fwd,
		int	pid
	)		const

Definition at line 285 of file TrajectoryMCSFitter.cxx.

                                                                                                                                                                              {
  //
  const int beg  = (fwd ? 0 : (dthetaij.size()-1));
  const int end  = (fwd ? dthetaij.size() : -1);
  const int incr = (fwd ? +1 : -1);
  //
  // bool print = false;
  //
  const double m = mass(pid);
  const double m2 = m*m;
  const double Etot = sqrt(p*p + m2);//Initial energy
  double Eij2 = 0.;
  //
  double const fixedterm = 0.5 * std::log( 2.0 * M_PI );
  double result = 0;
  for (int i = beg; i != end; i+=incr ) {
    if (dthetaij[i]<0) {
      //cout << "skip segment with too few points" << endl;
      continue;
    }
    //
      const double Eij = GetE(Etot,cumLen[i],m);
      Eij2 = Eij*Eij;
    //
    if ( Eij2 <= m2 ) {
      result = std::numeric_limits<double>::max();
      break;
    }
    const double pij = sqrt(Eij2 - m2);//momentum at this segment
    const double beta = sqrt( 1. - ((m2)/(pij*pij + m2)) );
    constexpr double HighlandSecondTerm = 0.038;
    const double tH0 = ( HighlandFirstTerm(pij) / (pij*beta) ) * ( 1.0 + HighlandSecondTerm * std::log( seg_nradl[i] ) ) * sqrt( seg_nradl[i] );
    const double rms = sqrt( 2.0*( tH0 * tH0 + theta0x * theta0x ) );
    if (rms==0.0) {
      std::cout << " Error : RMS cannot be zero ! " << std::endl;
      return std::numeric_limits<double>::max();
    }
    const double arg = dthetaij[i]/rms;
    result += ( std::log( rms ) + 0.5 * arg * arg + fixedterm);
  }
  return result;
}

int trkf::TrajectoryMCSFitter::minNSegs ( ) const

inline

Definition at line 185 of file TrajectoryMCSFitter.h.

{ return minNSegs_; }

double trkf::TrajectoryMCSFitter::segLen ( ) const

inline

Definition at line 186 of file TrajectoryMCSFitter.h.

{ return segLen_; }

double trkf::TrajectoryMCSFitter::segLenTolerance ( ) const

inline

Definition at line 187 of file TrajectoryMCSFitter.h.

{ return segLenTolerance_; }

Member Data Documentation

std::array<double, 5> trkf::TrajectoryMCSFitter::angResol_

private

Definition at line 201 of file TrajectoryMCSFitter.h.

bool trkf::TrajectoryMCSFitter::applySCEcorr_

private

Definition at line 204 of file TrajectoryMCSFitter.h.

int trkf::TrajectoryMCSFitter::eLossMode_

private

Definition at line 195 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::fineScanWindow_

private

Definition at line 200 of file TrajectoryMCSFitter.h.

std::array<double, 5> trkf::TrajectoryMCSFitter::hlParams_

private

Definition at line 202 of file TrajectoryMCSFitter.h.

int trkf::TrajectoryMCSFitter::minHitsPerSegment_

private

Definition at line 193 of file TrajectoryMCSFitter.h.

int trkf::TrajectoryMCSFitter::minNSegs_

private

Definition at line 191 of file TrajectoryMCSFitter.h.

int trkf::TrajectoryMCSFitter::nElossSteps_

private

Definition at line 194 of file TrajectoryMCSFitter.h.

int trkf::TrajectoryMCSFitter::pIdHyp_

private

Definition at line 190 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::pMax_

private

Definition at line 197 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::pMin_

private

Definition at line 196 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::pStep_

private

Definition at line 199 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::pStepCoarse_

private

Definition at line 198 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::segLen_

private

Definition at line 192 of file TrajectoryMCSFitter.h.

double trkf::TrajectoryMCSFitter::segLenTolerance_

private

Definition at line 203 of file TrajectoryMCSFitter.h.

---

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

• larreco/larreco/RecoAlg/TrajectoryMCSFitter.h
• larreco/larreco/RecoAlg/TrajectoryMCSFitter.cxx