LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
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 pStep, double angResol)
 
 TrajectoryMCSFitter (const Parameters &p)
 
recob::MCSFitResult fitMcs (const recob::TrackTrajectory &traj, bool momDepConst=true) const
 
recob::MCSFitResult fitMcs (const recob::Track &track, bool momDepConst=true) const
 
recob::MCSFitResult fitMcs (const recob::Trajectory &traj, bool momDepConst=true) const
 
recob::MCSFitResult fitMcs (const recob::TrackTrajectory &traj, int pid, bool momDepConst=true) const
 
recob::MCSFitResult fitMcs (const recob::Track &track, int pid, bool momDepConst=true) const
 
recob::MCSFitResult fitMcs (const recob::Trajectory &traj, int pid, bool momDepConst=true) 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, bool momDepConst, int pid) const
 
const ScanResult doLikelihoodScan (std::vector< float > &dtheta, std::vector< float > &seg_nradlengths, std::vector< float > &cumLen, bool fwdFit, bool momDepConst, int pid) const
 
double MomentumDependentConstant (const double p) 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
 

Private Attributes

int pIdHyp_
 
int minNSegs_
 
double segLen_
 
int minHitsPerSegment_
 
int nElossSteps_
 
int eLossMode_
 
double pMin_
 
double pMax_
 
double pStep_
 
double angResol_
 

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 segment (radiation) lengths, vector of scattering angles, and PID hypothesis used in the fit.

For configuration options see TrajectoryMCSFitter::Config

Author
G. Cerati (FNAL, MicroBooNE)
Date
2017
Version
1.0

Definition at line 33 of file TrajectoryMCSFitter.h.

Member Typedef Documentation

Constructor & Destructor Documentation

trkf::TrajectoryMCSFitter::TrajectoryMCSFitter ( int  pIdHyp,
int  minNSegs,
double  segLen,
int  minHitsPerSegment,
int  nElossSteps,
int  eLossMode,
double  pMin,
double  pMax,
double  pStep,
double  angResol 
)
inline

Definition at line 93 of file TrajectoryMCSFitter.h.

References trkf::TrajectoryMCSFitter::Config::angResol, angResol_, trkf::TrajectoryMCSFitter::Config::eLossMode, eLossMode_, trkf::TrajectoryMCSFitter::Config::minHitsPerSegment, minHitsPerSegment_, minNSegs_, trkf::TrajectoryMCSFitter::Config::nElossSteps, nElossSteps_, pIdHyp_, trkf::TrajectoryMCSFitter::Config::pMax, pMax_, trkf::TrajectoryMCSFitter::Config::pMin, pMin_, trkf::TrajectoryMCSFitter::Config::pStep, pStep_, and segLen_.

93  {
94  pIdHyp_ = pIdHyp;
95  minNSegs_ = minNSegs;
96  segLen_ = segLen;
97  minHitsPerSegment_ = minHitsPerSegment;
98  nElossSteps_ = nElossSteps;
99  eLossMode_ = eLossMode;
100  pMin_ = pMin;
101  pMax_ = pMax;
102  pStep_ = pStep;
103  angResol_ = angResol;
104  }
trkf::TrajectoryMCSFitter::TrajectoryMCSFitter ( const Parameters p)
inlineexplicit

Definition at line 105 of file TrajectoryMCSFitter.h.

106  : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {}
TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol)

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

References recob::TrackTrajectory::FirstValidPoint(), recob::TrackTrajectory::InvalidIndex, recob::TrackTrajectory::LastValidPoint(), recob::TrackTrajectory::Length(), recob::TrackTrajectory::LocationAtPoint(), max, recob::TrackTrajectory::NextValidPoint(), and R.

Referenced by fitMcs().

59  {
60  //
61  const double trajlen = traj.Length();
62  const int nseg = std::max(minNSegs_,int(trajlen/segLen_));
63  const double thisSegLen = trajlen/double(nseg);
64  // std::cout << "track with length=" << trajlen << " broken in nseg=" << nseg << " of length=" << thisSegLen << " where segLen_=" << segLen_ << std::endl;
65  //
66  constexpr double lar_radl_inv = 1./14.0;
67  cumseglens.push_back(0.);//first segment has zero cumulative length from previous segments
68  double thislen = 0.;
69  auto nextValid=traj.FirstValidPoint();
70  breakpoints.push_back(nextValid);
71  auto pos0 = traj.LocationAtPoint(nextValid);
72  nextValid = traj.NextValidPoint(nextValid+1);
73  int npoints = 0;
74  while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
75  auto pos1 = traj.LocationAtPoint(nextValid);
76  thislen += ( (pos1-pos0).R() );
77  pos0=pos1;
78  npoints++;
79  if (thislen>=thisSegLen) {
80  breakpoints.push_back(nextValid);
81  if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
82  else segradlengths.push_back(-999.);
83  cumseglens.push_back(cumseglens.back()+thislen);
84  thislen = 0.;
85  npoints = 0;
86  }
87  nextValid = traj.NextValidPoint(nextValid+1);
88  }
89  //then add last segment
90  if (thislen>0.) {
91  breakpoints.push_back(traj.LastValidPoint()+1);
92  segradlengths.push_back(thislen*lar_radl_inv);
93  cumseglens.push_back(cumseglens.back()+thislen);
94  }
95  return;
96 }
size_t LastValidPoint() const
Returns the index of the last valid point in the trajectory.
Int_t max
Definition: plot.C:27
double Length(size_t startAt=0) const
Returns the approximate length of the trajectory.
T LocationAtPoint(unsigned int p) const
Position at point p. Use e.g. as:
Double_t R
size_t FirstValidPoint() const
Returns the index of the first valid point in the trajectory.
static constexpr size_t InvalidIndex
Value returned on failed index queries.
size_t NextValidPoint(size_t index) const
Returns the index of the next valid point in the trajectory.
const TrajectoryMCSFitter::ScanResult TrajectoryMCSFitter::doLikelihoodScan ( std::vector< float > &  dtheta,
std::vector< float > &  seg_nradlengths,
std::vector< float > &  cumLen,
bool  fwdFit,
bool  momDepConst,
int  pid 
) const

Definition at line 98 of file TrajectoryMCSFitter.cxx.

References max.

98  {
99  int best_idx = -1;
100  double best_logL = std::numeric_limits<double>::max();
101  double best_p = -1.0;
102  std::vector<float> vlogL;
103  for (double p_test = pMin_; p_test <= pMax_; p_test+=pStep_) {
104  double logL = mcsLikelihood(p_test, angResol_, dtheta, seg_nradlengths, cumLen, fwdFit, momDepConst, pid);
105  if (logL < best_logL) {
106  best_p = p_test;
107  best_logL = logL;
108  best_idx = vlogL.size();
109  }
110  vlogL.push_back(logL);
111  }
112  //
113  //uncertainty from left side scan
114  double lunc = -1.0;
115  if (best_idx>0) {
116  for (int j=best_idx-1;j>=0;j--) {
117  double dLL = vlogL[j]-vlogL[best_idx];
118  if ( dLL<0.5 ) {
119  lunc = (best_idx-j)*pStep_;
120  } else break;
121  }
122  }
123  //uncertainty from right side scan
124  double runc = -1.0;
125  if (best_idx<int(vlogL.size()-1)) {
126  for (unsigned int j=best_idx+1;j<vlogL.size();j++) {
127  double dLL = vlogL[j]-vlogL[best_idx];
128  if ( dLL<0.5 ) {
129  runc = (j-best_idx)*pStep_;
130  } else break;
131  }
132  }
133  return ScanResult(best_p, std::max(lunc,runc), best_logL);
134 }
double mcsLikelihood(double p, double theta0x, std::vector< float > &dthetaij, std::vector< float > &seg_nradl, std::vector< float > &cumLen, bool fwd, bool momDepConst, int pid) const
Int_t max
Definition: plot.C:27
double TrajectoryMCSFitter::energyLossBetheBloch ( const double  mass,
const double  e2 
) const

Definition at line 259 of file TrajectoryMCSFitter.cxx.

References E.

Referenced by mass().

259  {
260  // stolen, mostly, from GFMaterialEffects.
261  constexpr double Iinv = 1./188.E-6;
262  constexpr double matConst = 1.4*18./40.;//density*Z/A
263  constexpr double me = 0.511;
264  constexpr double kappa = 0.307075;
265  //
266  const double beta2 = (e2-mass*mass)/e2;
267  const double gamma2 = 1./(1.0 - beta2);
268  const double massRatio = me/mass;
269  const double argument = (2.*me*gamma2*beta2*Iinv) * std::sqrt(1+2*std::sqrt(gamma2)*massRatio + massRatio*massRatio);
270  //
271  double dedx = kappa*matConst/beta2;
272  //
273  if (mass==0.0) return(0.0);
274  if (argument <= exp(beta2)) {
275  dedx = 0.;
276  } else{
277  dedx *= (log(argument)-beta2)*1.E-3; // Bethe-Bloch, converted to GeV/cm
278  if (dedx<0.) dedx = 0.;
279  }
280  return dedx;
281 }
Float_t E
Definition: plot.C:23
double mass(int pid) const
double TrajectoryMCSFitter::energyLossLandau ( const double  mass2,
const double  E2,
const double  x 
) const

Definition at line 241 of file TrajectoryMCSFitter.cxx.

Referenced by mass().

241  {
242  //
243  // eq. (33.11) in http://pdg.lbl.gov/2016/reviews/rpp2016-rev-passage-particles-matter.pdf (except density correction is ignored)
244  //
245  if (x<=0.) return 0.;
246  constexpr double Iinv2 = 1./(188.E-6*188.E-6);
247  constexpr double matConst = 1.4*18./40.;//density*Z/A
248  constexpr double me = 0.511;
249  constexpr double kappa = 0.307075;
250  constexpr double j = 0.200;
251  //
252  const double beta2 = (e2-mass2)/e2;
253  const double gamma2 = 1./(1.0 - beta2);
254  const double epsilon = 0.5*kappa*x*matConst/beta2;
255  //
256  return 0.001*epsilon*( log(2.*me*beta2*gamma2*epsilon*Iinv2) + j - beta2 );
257 }
Float_t x
Definition: compare.C:6
recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::TrackTrajectory traj,
bool  momDepConst = true 
) const
inline

Definition at line 108 of file TrajectoryMCSFitter.h.

References fitMcs(), and pIdHyp_.

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

108 { return fitMcs(traj,pIdHyp_,momDepConst); }
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Track track,
bool  momDepConst = true 
) const
inline

Definition at line 109 of file TrajectoryMCSFitter.h.

References fitMcs(), and pIdHyp_.

Referenced by fitMcs().

109 { return fitMcs(track,pIdHyp_,momDepConst); }
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory traj,
bool  momDepConst = true 
) const
inline

Definition at line 110 of file TrajectoryMCSFitter.h.

References fitMcs(), and pIdHyp_.

Referenced by fitMcs().

110 { return fitMcs(traj,pIdHyp_,momDepConst); }
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
recob::MCSFitResult TrajectoryMCSFitter::fitMcs ( const recob::TrackTrajectory traj,
int  pid,
bool  momDepConst = true 
) const

Definition at line 11 of file TrajectoryMCSFitter.cxx.

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

11  {
12  //
13  // Break the trajectory in segments of length approximately equal to segLen_
14  //
15  vector<size_t> breakpoints;
16  vector<float> segradlengths;
17  vector<float> cumseglens;
18  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
19  //
20  // Fit segment directions, and get 3D angles between them
21  //
22  if (segradlengths.size()<2) return recob::MCSFitResult();
23  vector<float> dtheta;
24  Vector_t pcdir0;
25  Vector_t pcdir1;
26  for (unsigned int p = 0; p<segradlengths.size(); p++) {
27  linearRegression(traj, breakpoints[p], breakpoints[p+1], pcdir1);
28  if (p>0) {
29  if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
30  dtheta.push_back(-999.);
31  } else {
32  const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
33  //assert(std::abs(cosval)<=1);
34  //units are mrad
35  double dt = 1000.*acos(cosval);//should we try to use expansion for small angles?
36  dtheta.push_back(dt);
37  }
38  }
39  pcdir0 = pcdir1;
40  }
41  //
42  // Perform likelihood scan in forward and backward directions
43  //
44  vector<float> cumLenFwd;
45  vector<float> cumLenBwd;
46  for (unsigned int i = 0; i<cumseglens.size()-2; i++) {
47  cumLenFwd.push_back(cumseglens[i]);
48  cumLenBwd.push_back(cumseglens.back()-cumseglens[i+2]);
49  }
50  const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true, momDepConst, pid);
51  const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid);
52  //
53  return recob::MCSFitResult(pid,
54  fwdResult.p,fwdResult.pUnc,fwdResult.logL,
55  bwdResult.p,bwdResult.pUnc,bwdResult.logL,
56  segradlengths,dtheta);
57 }
void breakTrajInSegments(const recob::TrackTrajectory &traj, std::vector< size_t > &breakpoints, std::vector< float > &segradlengths, std::vector< float > &cumseglens) const
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
void linearRegression(const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) const
Class storing the result of the Maximum Likelihood fit of Multiple Coulomb Scattering angles between ...
Definition: MCSFitResult.h:19
const ScanResult doLikelihoodScan(std::vector< float > &dtheta, std::vector< float > &seg_nradlengths, std::vector< float > &cumLen, bool fwdFit, bool momDepConst, int pid) const
recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Track track,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 113 of file TrajectoryMCSFitter.h.

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

Referenced by fitMcs().

113 { return fitMcs(track.Trajectory(),pid,momDepConst); }
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:101
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
recob::MCSFitResult trkf::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory traj,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 114 of file TrajectoryMCSFitter.h.

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

114  {
116  const recob::TrackTrajectory tt(traj,std::move(flags));
117  return fitMcs(tt,pid,momDepConst);
118  }
Definition: type_traits.h:56
size_t NPoints() const
Returns the number of stored trajectory points.
Definition: Trajectory.h:173
A trajectory in space reconstructed from hits.
std::vector< PointFlags_t > Flags_t
Type of point flag list.
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
double TrajectoryMCSFitter::GetE ( const double  initial_E,
const double  length_travelled,
const double  mass 
) const

Definition at line 283 of file TrajectoryMCSFitter.cxx.

Referenced by mass().

283  {
284  //
285  const double step_size = length_travelled / nElossSteps_;
286  //
287  double current_E = initial_E;
288  const double m2 = m*m;
289  //
290  for (auto i = 0; i < nElossSteps_; ++i) {
291  if (eLossMode_==2) {
292  double dedx = energyLossBetheBloch(m,current_E);
293  current_E -= (dedx * step_size);
294  } else {
295  // MPV of Landau energy loss distribution
296  current_E -= energyLossLandau(m2,current_E*current_E,step_size);
297  }
298  if ( current_E <= m ) {
299  // std::cout<<"WARNING: current_E less than mu mass. it is "<<current_E<<std::endl;
300  return 0.;
301  }
302  }
303  return current_E;
304 }
double energyLossBetheBloch(const double mass, const double e2) const
double energyLossLandau(const double mass2, const double E2, const double x) const
void TrajectoryMCSFitter::linearRegression ( const recob::TrackTrajectory traj,
const size_t  firstPoint,
const size_t  lastPoint,
recob::tracking::Vector_t pcdir 
) const

Definition at line 136 of file TrajectoryMCSFitter.cxx.

References geo::vect::MiddlePointAccumulatorDim< N >::add(), recob::TrackTrajectory::DirectionAtPoint(), recob::TrackTrajectory::LocationAtPoint(), geo::vect::MiddlePointAccumulatorDim< N >::middlePoint(), recob::TrackTrajectory::NextValidPoint(), and norm.

Referenced by fitMcs().

136  {
137  //
138  int npoints = 0;
139  geo::vect::MiddlePointAccumulator middlePointCalc;
140  size_t nextValid = firstPoint;
141  while (nextValid<lastPoint) {
142  middlePointCalc.add(traj.LocationAtPoint(nextValid));
143  nextValid = traj.NextValidPoint(nextValid+1);
144  npoints++;
145  }
146  const auto avgpos = middlePointCalc.middlePoint();
147  const double norm = 1./double(npoints);
148  //
149  //assert(npoints>0);
150  //
151  TMatrixDSym m(3);
152  nextValid = firstPoint;
153  while (nextValid<lastPoint) {
154  const auto p = traj.LocationAtPoint(nextValid);
155  const double xxw0 = p.X()-avgpos.X();
156  const double yyw0 = p.Y()-avgpos.Y();
157  const double zzw0 = p.Z()-avgpos.Z();
158  m(0, 0) += xxw0*xxw0*norm;
159  m(0, 1) += xxw0*yyw0*norm;
160  m(0, 2) += xxw0*zzw0*norm;
161  m(1, 0) += yyw0*xxw0*norm;
162  m(1, 1) += yyw0*yyw0*norm;
163  m(1, 2) += yyw0*zzw0*norm;
164  m(2, 0) += zzw0*xxw0*norm;
165  m(2, 1) += zzw0*yyw0*norm;
166  m(2, 2) += zzw0*zzw0*norm;
167  nextValid = traj.NextValidPoint(nextValid+1);
168  }
169  //
170  const TMatrixDSymEigen me(m);
171  const auto& eigenval = me.GetEigenValues();
172  const auto& eigenvec = me.GetEigenVectors();
173  //
174  int maxevalidx = 0;
175  double maxeval = eigenval(0);
176  for (int i=1; i<3; ++i) {
177  if (eigenval(i)>maxeval) {
178  maxevalidx = i;
179  maxeval = eigenval(i);
180  }
181  }
182  //
183  pcdir = Vector_t(eigenvec(0, maxevalidx),eigenvec(1, maxevalidx),eigenvec(2, maxevalidx));
184  if (traj.DirectionAtPoint(firstPoint).Dot(pcdir)<0.) pcdir*=-1.;
185  //
186 }
T DirectionAtPoint(unsigned int p) const
Direction at point p. Use e.g. as:
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
Helper class to compute the middle point in a point set.
void add(Point const &p)
Accumulates a point.
T LocationAtPoint(unsigned int p) const
Position at point p. Use e.g. as:
Float_t norm
size_t NextValidPoint(size_t index) const
Returns the index of the next valid point in the trajectory.
double trkf::TrajectoryMCSFitter::mass ( int  pid) const
inline

Definition at line 138 of file TrajectoryMCSFitter.h.

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

138  {
139  if (abs(pid)==13) { return mumass; }
140  if (abs(pid)==211) { return pimass; }
141  if (abs(pid)==321) { return kmass; }
142  if (abs(pid)==2212) { return pmass; }
143  return util::kBogusD;
144  }
constexpr double kBogusD
obviously bogus double value
double TrajectoryMCSFitter::mcsLikelihood ( double  p,
double  theta0x,
std::vector< float > &  dthetaij,
std::vector< float > &  seg_nradl,
std::vector< float > &  cumLen,
bool  fwd,
bool  momDepConst,
int  pid 
) const

Definition at line 188 of file TrajectoryMCSFitter.cxx.

References beta, evd::details::end(), and max.

Referenced by fitMcs().

188  {
189  //
190  const int beg = (fwd ? 0 : (dthetaij.size()-1));
191  const int end = (fwd ? dthetaij.size() : -1);
192  const int incr = (fwd ? +1 : -1);
193  //
194  // bool print = false;//(p>1.999 && p<2.001);
195  //
196  const double m = mass(pid);
197  const double m2 = m*m;
198  const double Etot = sqrt(p*p + m2);//Initial energy
199  double Eij2 = 0.;
200  //
201  double const fixedterm = 0.5 * std::log( 2.0 * M_PI );
202  double result = 0;
203  for (int i = beg; i != end; i+=incr ) {
204  if (dthetaij[i]<0) {
205  //cout << "skip segment with too few points" << endl;
206  continue;
207  }
208  //
209  if (eLossMode_==1) {
210  // ELoss mode: MIP (constant)
211  constexpr double kcal = 0.002105;
212  const double Eij = Etot - kcal*cumLen[i];//energy at this segment
213  Eij2 = Eij*Eij;
214  } else {
215  // Non constant energy loss distribution
216  const double Eij = GetE(Etot,cumLen[i],m);
217  Eij2 = Eij*Eij;
218  }
219  //
220  if ( Eij2 <= m2 ) {
222  break;
223  }
224  const double pij = sqrt(Eij2 - m2);//momentum at this segment
225  const double beta = sqrt( 1. - ((m2)/(pij*pij + m2)) );
226  constexpr double tuned_HL_term1 = 11.0038; // https://arxiv.org/abs/1703.06187
227  constexpr double HL_term2 = 0.038;
228  const double tH0 = ( (momDepConst ? MomentumDependentConstant(pij) : tuned_HL_term1) / (pij*beta) ) * ( 1.0 + HL_term2 * std::log( seg_nradl[i] ) ) * sqrt( seg_nradl[i] );
229  const double rms = sqrt( 2.0*( tH0 * tH0 + theta0x * theta0x ) );
230  if (rms==0.0) {
231  std::cout << " Error : RMS cannot be zero ! " << std::endl;
233  }
234  const double arg = dthetaij[i]/rms;
235  result += ( std::log( rms ) + 0.5 * arg * arg + fixedterm);
236  // if (print && fwd==true) cout << "TrajectoryMCSFitter pij=" << pij << " dthetaij[i]=" << dthetaij[i] << " tH0=" << tH0 << " rms=" << rms << " prob=" << ( std::log( rms ) + 0.5 * arg * arg + fixedterm) << " const=" << (momDepConst ? MomentumDependentConstant(pij) : tuned_HL_term1) << " beta=" << beta << " red_length=" << seg_nradl[i] << endl;
237  }
238  return result;
239 }
STL namespace.
Double_t beta
Int_t max
Definition: plot.C:27
double mass(int pid) const
double GetE(const double initial_E, const double length_travelled, const double mass) const
double MomentumDependentConstant(const double p) const
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
double trkf::TrajectoryMCSFitter::MomentumDependentConstant ( const double  p) const
inline

Definition at line 132 of file TrajectoryMCSFitter.h.

132  {
133  //these are from https://arxiv.org/abs/1703.06187
134  constexpr double a = 0.1049;
135  constexpr double c = 11.0038;
136  return (a/(p*p)) + c;
137  }

Member Data Documentation

double trkf::TrajectoryMCSFitter::angResol_
private

Definition at line 160 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::eLossMode_
private

Definition at line 156 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::minHitsPerSegment_
private

Definition at line 154 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::minNSegs_
private

Definition at line 152 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::nElossSteps_
private

Definition at line 155 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

int trkf::TrajectoryMCSFitter::pIdHyp_
private

Definition at line 151 of file TrajectoryMCSFitter.h.

Referenced by fitMcs(), and TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pMax_
private

Definition at line 158 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pMin_
private

Definition at line 157 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::pStep_
private

Definition at line 159 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().

double trkf::TrajectoryMCSFitter::segLen_
private

Definition at line 153 of file TrajectoryMCSFitter.h.

Referenced by TrajectoryMCSFitter().


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