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 44 of file TrajectoryMCSFitter.h.

Member Typedef Documentation

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

Definition at line 114 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 116 of file TrajectoryMCSFitter.h.

References trkf::TrajectoryMCSFitter::Config::angResol, angResol_, trkf::TrajectoryMCSFitter::Config::applySCEcorr, applySCEcorr_, trkf::TrajectoryMCSFitter::Config::eLossMode, eLossMode_, trkf::TrajectoryMCSFitter::Config::fineScanWindow, fineScanWindow_, trkf::TrajectoryMCSFitter::Config::hlParams, hlParams_, trkf::TrajectoryMCSFitter::Config::minHitsPerSegment, minHitsPerSegment_, minNSegs(), minNSegs_, trkf::TrajectoryMCSFitter::Config::nElossSteps, nElossSteps_, pIdHyp_, trkf::TrajectoryMCSFitter::Config::pMax, pMax_, trkf::TrajectoryMCSFitter::Config::pMin, pMin_, trkf::TrajectoryMCSFitter::Config::pStep, pStep_, trkf::TrajectoryMCSFitter::Config::pStepCoarse, pStepCoarse_, segLen(), segLen_, trkf::TrajectoryMCSFitter::Config::segLenTolerance, and segLenTolerance_.

    {
      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 148 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 73 of file TrajectoryMCSFitter.cxx.

References recob::TrackTrajectory::DirectionAtPoint(), geo::GeometryCore::FindTPCAtPosition(), recob::TrackTrajectory::FirstValidPoint(), recob::TrackTrajectory::InvalidIndex, geo::CryostatID::isValid, recob::TrackTrajectory::LastValidPoint(), recob::TrackTrajectory::Length(), recob::TrackTrajectory::LocationAtPoint(), recob::TrackTrajectory::NextValidPoint(), and geo::TPCID::TPC.

Referenced by fitMcs().

{
  //
  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_) {
    geo::TPCID tpcid = geom->FindTPCAtPosition(pos0);
    geo::Vector_t pos0_offset(0., 0., 0.);
    if (tpcid.isValid) { pos0_offset = _SCE->GetCalPosOffsets(pos0, 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_) {
      geo::TPCID tpcid = geom->FindTPCAtPosition(pos1);
      geo::Vector_t pos1_offset(0., 0., 0.);
      if (tpcid.isValid) { pos1_offset = _SCE->GetCalPosOffsets(pos1, 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 231 of file TrajectoryMCSFitter.h.

References angResol_.

    {
      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 203 of file TrajectoryMCSFitter.cxx.

References trkf::TrajectoryMCSFitter::ScanResult::p, and trkf::TrajectoryMCSFitter::ScanResult::pUnc.

{

  //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 153 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 381 of file TrajectoryMCSFitter.cxx.

References E.

Referenced by mass().

{
  // 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 360 of file TrajectoryMCSFitter.cxx.

Referenced by mass().

{
  //
  // 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 166 of file TrajectoryMCSFitter.h.

References pIdHyp_.

Referenced by fitMcs(), trkmkr::KalmanFilterFitTrackMaker::getMomentum(), and trkf::MCSFitProducer::produce().

    {
      return fitMcs(traj, pIdHyp_);
    }

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

inline

Definition at line 170 of file TrajectoryMCSFitter.h.

References fitMcs(), and pIdHyp_.

Referenced by fitMcs().

{ return fitMcs(track, pIdHyp_); }

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

inline

Definition at line 171 of file TrajectoryMCSFitter.h.

References fitMcs(), and pIdHyp_.

    {
      return fitMcs(traj, pIdHyp_);
    }

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

Definition at line 16 of file TrajectoryMCSFitter.cxx.

References util::abs(), trkf::TrajectoryMCSFitter::ScanResult::logL, trkf::TrajectoryMCSFitter::ScanResult::p, trkf::TrajectoryMCSFitter::ScanResult::pUnc, and recob::TrackTrajectory::StartDirection().

{
  //
  // 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 177 of file TrajectoryMCSFitter.h.

References fitMcs(), and recob::Track::Trajectory().

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

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

inline

Definition at line 181 of file TrajectoryMCSFitter.h.

References breakTrajInSegments(), fitMcs(), linearRegression(), mcsLikelihood(), and recob::Trajectory::NPoints().

    {
      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 406 of file TrajectoryMCSFitter.cxx.

Referenced by mass().

{
  //
  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 226 of file TrajectoryMCSFitter.h.

References hlParams_.

    {
      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 230 of file TrajectoryMCSFitter.cxx.

References geo::vect::MiddlePointAccumulatorDim< N >::add(), recob::TrackTrajectory::DirectionAtPoint(), geo::GeometryCore::FindTPCAtPosition(), geo::CryostatID::isValid, recob::TrackTrajectory::LocationAtPoint(), geo::vect::MiddlePointAccumulatorDim< N >::middlePoint(), recob::TrackTrajectory::NextValidPoint(), norm, and geo::TPCID::TPC.

Referenced by fitMcs().

{
  //
  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_) {
      geo::TPCID tpcid = geom->FindTPCAtPosition(tempP);
      geo::Vector_t tempP_offset(0., 0., 0.);
      if (tpcid.isValid) { tempP_offset = _SCE->GetCalPosOffsets(tempP, 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_) {
      geo::TPCID tpcid = geom->FindTPCAtPosition(p);
      geo::Vector_t p_offset(0., 0., 0.);
      if (tpcid.isValid) { p_offset = _SCE->GetCalPosOffsets(p, 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 236 of file TrajectoryMCSFitter.h.

References util::abs(), energyLossBetheBloch(), energyLossLandau(), GetE(), util::kBogusD, and x.

    {
      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 309 of file TrajectoryMCSFitter.cxx.

References util::end().

Referenced by fitMcs().

{
  //
  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 249 of file TrajectoryMCSFitter.h.

References minNSegs_.

Referenced by TrajectoryMCSFitter().

{ return minNSegs_; }

double trkf::TrajectoryMCSFitter::segLen ( ) const

inline

Definition at line 250 of file TrajectoryMCSFitter.h.

References segLen_.

Referenced by TrajectoryMCSFitter().

{ return segLen_; }

double trkf::TrajectoryMCSFitter::segLenTolerance ( ) const

inline

Definition at line 251 of file TrajectoryMCSFitter.h.

References segLenTolerance_.

{ return segLenTolerance_; }

Member Data Documentation

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

private

Definition at line 265 of file TrajectoryMCSFitter.h.

Referenced by DetectorAngularResolution(), and TrajectoryMCSFitter().

bool trkf::TrajectoryMCSFitter::applySCEcorr_

private

Definition at line 268 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::eLossMode_

private

Definition at line 259 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::fineScanWindow_

private

Definition at line 264 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

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

private

Definition at line 266 of file TrajectoryMCSFitter.h.

Referenced by HighlandFirstTerm(), and TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::minHitsPerSegment_

private

Definition at line 257 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::minNSegs_

private

Definition at line 255 of file TrajectoryMCSFitter.h.

Referenced by minNSegs(), and TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::nElossSteps_

private

Definition at line 258 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::pIdHyp_

private

Definition at line 254 of file TrajectoryMCSFitter.h.

Referenced by fitMcs(), and TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pMax_

private

Definition at line 261 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pMin_

private

Definition at line 260 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pStep_

private

Definition at line 263 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pStepCoarse_

private

Definition at line 262 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::segLen_

private

Definition at line 256 of file TrajectoryMCSFitter.h.

Referenced by segLen(), and TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::segLenTolerance_

private

Definition at line 267 of file TrajectoryMCSFitter.h.

Referenced by segLenTolerance(), and TrajectoryMCSFitter().

---

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

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