All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Types | Public Member Functions | Private Attributes | List of all members
trkf::sbn::TrajectoryMCSFitter Class Reference

#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

Definition at line 33 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

Member Typedef Documentation

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

Definition at line 91 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

Constructor & Destructor Documentation

trkf::sbn::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 sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

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::sbn::TrajectoryMCSFitter::TrajectoryMCSFitter ( const Parameters &  p)
inlineexplicit

Definition at line 105 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

106  : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {}
p
pdgs p
Definition: selectors.fcl:22
trkf::sbn::TrajectoryMCSFitter::TrajectoryMCSFitter
TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol)
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h:93

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 sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
recob::TrackTrajectory::LastValidPoint
size_t LastValidPoint() const
Returns the index of the last valid point in the trajectory.
Definition: TrackTrajectory.h:235
recob::TrackTrajectory::Length
double Length(size_t startAt=0) const
Returns the approximate length of the trajectory.
Definition: TrackTrajectory.cxx:70
recob::TrackTrajectory::LocationAtPoint
T LocationAtPoint(unsigned int p) const
Position at point p. Use e.g. as:
Definition: TrackTrajectory.h:547
recob::TrackTrajectory::FirstValidPoint
size_t FirstValidPoint() const
Returns the index of the first valid point in the trajectory.
Definition: TrackTrajectory.h:197
recob::TrackTrajectory::InvalidIndex
static constexpr size_t InvalidIndex
Value returned on failed index queries.
Definition: TrackTrajectory.h:661
recob::TrackTrajectory::NextValidPoint
size_t NextValidPoint(size_t index) const
Returns the index of the next valid point in the trajectory.
Definition: TrackTrajectory.h:211
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 sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
trkf::sbn::TrajectoryMCSFitter::mcsLikelihood
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
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:188
double TrajectoryMCSFitter::energyLossBetheBloch ( const double  mass,
const double  e2 
) const

Definition at line 259 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
E
process_name E
Definition: rawtorootconvert.fcl:26
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
double TrajectoryMCSFitter::energyLossLandau ( const double  mass2,
const double  E2,
const double  x 
) const

Definition at line 241 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
recob::MCSFitResult trkf::sbn::TrajectoryMCSFitter::fitMcs ( const recob::TrackTrajectory &  traj,
bool  momDepConst = true 
) const
inline

Definition at line 108 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

108 { return fitMcs(traj,pIdHyp_,momDepConst); }
trkf::sbn::TrajectoryMCSFitter::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h:108
recob::MCSFitResult trkf::sbn::TrajectoryMCSFitter::fitMcs ( const recob::Track &  track,
bool  momDepConst = true 
) const
inline

Definition at line 109 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

109 { return fitMcs(track,pIdHyp_,momDepConst); }
trkf::sbn::TrajectoryMCSFitter::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h:108
recob::MCSFitResult trkf::sbn::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory &  traj,
bool  momDepConst = true 
) const
inline

Definition at line 110 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

110 { return fitMcs(traj,pIdHyp_,momDepConst); }
trkf::sbn::TrajectoryMCSFitter::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h:108
recob::MCSFitResult TrajectoryMCSFitter::fitMcs ( const recob::TrackTrajectory &  traj,
int  pid,
bool  momDepConst = true 
) const

Definition at line 11 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
trkf::sbn::TrajectoryMCSFitter::breakTrajInSegments
void breakTrajInSegments(const recob::TrackTrajectory &traj, std::vector< size_t > &breakpoints, std::vector< float > &segradlengths, std::vector< float > &cumseglens) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:59
p
pdgs p
Definition: selectors.fcl:22
trkf::sbn::TrajectoryMCSFitter::linearRegression
void linearRegression(const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:136
trkf::sbn::TrajectoryMCSFitter::doLikelihoodScan
const ScanResult doLikelihoodScan(std::vector< float > &dtheta, std::vector< float > &seg_nradlengths, std::vector< float > &cumLen, bool fwdFit, bool momDepConst, int pid) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:98
recob::tracking::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space. See recob::tracking::Coord_t for more details on the ...
Definition: TrackingTypes.h:29
recob::MCSFitResult
Class storing the result of the Maximum Likelihood fit of Multiple Coulomb Scattering angles between ...
Definition: MCSFitResult.h:19
recob::MCSFitResult trkf::sbn::TrajectoryMCSFitter::fitMcs ( const recob::Track &  track,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 113 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

113 { return fitMcs(track.Trajectory(),pid,momDepConst); }
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: bj/lardataobj/RecoBase/Track.h:98
trkf::sbn::TrajectoryMCSFitter::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h:108
recob::MCSFitResult trkf::sbn::TrajectoryMCSFitter::fitMcs ( const recob::Trajectory &  traj,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 114 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

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

Definition at line 283 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
trkf::sbn::TrajectoryMCSFitter::energyLossLandau
double energyLossLandau(const double mass2, const double E2, const double x) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:241
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
trkf::sbn::TrajectoryMCSFitter::energyLossBetheBloch
double energyLossBetheBloch(const double mass, const double e2) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:259
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 sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 }
recob::TrackTrajectory::DirectionAtPoint
T DirectionAtPoint(unsigned int p) const
Direction at point p. Use e.g. as:
Definition: TrackTrajectory.h:559
trkf::Vector_t
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
geo::vect::MiddlePointAccumulatorDim
Helper class to compute the middle point in a point set.
Definition: geo_vectors_utils.h:1246
p
pdgs p
Definition: selectors.fcl:22
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
geo::vect::MiddlePointAccumulatorDim::add
void add(Point const &p)
Accumulates a point.
Definition: geo_vectors_utils.h:1316
recob::TrackTrajectory::LocationAtPoint
T LocationAtPoint(unsigned int p) const
Position at point p. Use e.g. as:
Definition: TrackTrajectory.h:547
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:1202
recob::TrackTrajectory::NextValidPoint
size_t NextValidPoint(size_t index) const
Returns the index of the next valid point in the trajectory.
Definition: TrackTrajectory.h:211
double trkf::sbn::TrajectoryMCSFitter::mass ( int  pid) const
inline

Definition at line 138 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

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  }
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
util::kBogusD
constexpr double kBogusD
obviously bogus double value
Definition: PhysicalConstants.h:71
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 sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx.

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 ) {
221  result = std::numeric_limits<double>::max();
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;
232  return std::numeric_limits<double>::max();
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 }
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
p
pdgs p
Definition: selectors.fcl:22
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
trkf::sbn::TrajectoryMCSFitter::GetE
double GetE(const double initial_E, const double length_travelled, const double mass) const
Definition: sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx:283
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
double trkf::sbn::TrajectoryMCSFitter::MomentumDependentConstant ( const double  p) const
inline

Definition at line 132 of file sbncode/sbncode/LArRecoProducer/LArReco/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  }
p
pdgs p
Definition: selectors.fcl:22
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37

Member Data Documentation

double trkf::sbn::TrajectoryMCSFitter::angResol_
private

Definition at line 160 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

int trkf::sbn::TrajectoryMCSFitter::eLossMode_
private

Definition at line 156 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

int trkf::sbn::TrajectoryMCSFitter::minHitsPerSegment_
private

Definition at line 154 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

int trkf::sbn::TrajectoryMCSFitter::minNSegs_
private

Definition at line 152 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

int trkf::sbn::TrajectoryMCSFitter::nElossSteps_
private

Definition at line 155 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

int trkf::sbn::TrajectoryMCSFitter::pIdHyp_
private

Definition at line 151 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

double trkf::sbn::TrajectoryMCSFitter::pMax_
private

Definition at line 158 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

double trkf::sbn::TrajectoryMCSFitter::pMin_
private

Definition at line 157 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

double trkf::sbn::TrajectoryMCSFitter::pStep_
private

Definition at line 159 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

double trkf::sbn::TrajectoryMCSFitter::segLen_
private

Definition at line 153 of file sbncode/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5