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::TrajectoryMCSFitterICARUS Class Reference

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

#include <TrajectoryMCSFitterICARUS.h>

Classes

struct  Config
 
struct  ScanResult
 

Public Types

using Parameters = fhicl::Table< Config >
 

Public Member Functions

 TrajectoryMCSFitterICARUS (int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol)
 
 TrajectoryMCSFitterICARUS (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 findSegmentBarycenter (const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) 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
 
double GetOptimalSegLen (const double guess_p, const int n_points, const int plane, const double length_travelled) const
 
double computeResidual (int i, double &alfa) const
 
void ComputeD3P ()
 
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
 
void set2DHits (std::vector< recob::Hit > h)
 

Private Attributes

int pIdHyp_
 
int minNSegs_
 
double segLen_
 
int minHitsPerSegment_
 
int nElossSteps_
 
int eLossMode_
 
double pMin_
 
double pMax_
 
double pStep_
 
double angResol_
 
std::vector< recob::Hit > hits2d
 
float d3p
 

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 TrajectoryMCSFitterICARUS::Configs

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

Definition at line 34 of file TrajectoryMCSFitterICARUS.h.

Member Typedef Documentation

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

Definition at line 93 of file TrajectoryMCSFitterICARUS.h.

Constructor & Destructor Documentation

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

Definition at line 95 of file TrajectoryMCSFitterICARUS.h.

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

Definition at line 107 of file TrajectoryMCSFitterICARUS.h.

108  : TrajectoryMCSFitterICARUS(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::TrajectoryMCSFitterICARUS::TrajectoryMCSFitterICARUS
TrajectoryMCSFitterICARUS(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol)
Definition: TrajectoryMCSFitterICARUS.h:95

Member Function Documentation

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

Definition at line 112 of file TrajectoryMCSFitterICARUS.cxx.

112  {
113  //
114  const double trajlen = traj.Length();
115  const int nseg = std::max(minNSegs_,int(trajlen/segLen_));
116  const double thisSegLen = trajlen/double(nseg);
117  // std::cout << "track with length=" << trajlen << " broken in nseg=" << nseg << " of length=" << thisSegLen << " where segLen_=" << segLen_ << std::endl;
118  //
119  constexpr double lar_radl_inv = 1./14.0;
120  cumseglens.push_back(0.);//first segment has zero cumulative length from previous segments
121  double thislen = 0.;
122  auto nextValid=traj.FirstValidPoint();
123  breakpoints.push_back(nextValid);
124  auto pos0 = traj.LocationAtPoint(nextValid);
125  nextValid = traj.NextValidPoint(nextValid+1);
126  int npoints = 0;
127  while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
128  auto pos1 = traj.LocationAtPoint(nextValid);
129  thislen += ( (pos1-pos0).R() );
130  pos0=pos1;
131  npoints++;
132  if (thislen>=thisSegLen) {
133  breakpoints.push_back(nextValid);
134  if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
135  else segradlengths.push_back(-999.);
136  cumseglens.push_back(cumseglens.back()+thislen);
137  thislen = 0.;
138  npoints = 0;
139  }
140  nextValid = traj.NextValidPoint(nextValid+1);
141  }
142  //then add last segment
143  if (thislen>0.) {
144  breakpoints.push_back(traj.LastValidPoint()+1);
145  segradlengths.push_back(thislen*lar_radl_inv);
146  cumseglens.push_back(cumseglens.back()+thislen);
147  }
148  return;
149 }
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
void TrajectoryMCSFitterICARUS::ComputeD3P ( )

for each triplet of consecutive hits, save absolute value of residuale

Definition at line 419 of file TrajectoryMCSFitterICARUS.cxx.

419  {
420 {
421  double res;
422  vector<double> h0;
423  vector<double> alf;
424 
425  TH1D* hd3pv=new TH1D("hd3pv","hd3pv",100,-5.,5.);
426 
427  for (unsigned int j=0;j<hits2d.size()-2;j+=3) {
428 
429  double a;
430  res=computeResidual(j,a);
431  //cout << " point " << j << " residual " <<res << endl;
432  //! for each triplet of consecutive hits, save absolute value of residuale
433  if(abs(res)<5) { //mm
434  h0.push_back(res);
435  alf.push_back(a);
436  hd3pv->Fill(res);
437 
438  }}
439 
440  bool writeHisto=true;
441  if(writeHisto) {
442  TFile *f = new TFile("d3phisto.root","UPDATE");
443  hd3pv->Write();
444  f->Close();
445  f->Delete();
446  }
447 
448  if(!h0.size())
449  d3p=0.4;
450  else
451  d3p=hd3pv->GetRMS();
452 
453  //d3pvector=h0;
454  //return d3p;
455 }
456 }
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
trkf::TrajectoryMCSFitterICARUS::computeResidual
double computeResidual(int i, double &alfa) const
Definition: TrajectoryMCSFitterICARUS.cxx:457
double TrajectoryMCSFitterICARUS::computeResidual ( int  i,
double &  alfa 
) const

Definition at line 457 of file TrajectoryMCSFitterICARUS.cxx.

458 {
459  /*const auto p0 = traj.LocationAtPoint(i);
460  const auto p1 = traj.LocationAtPoint(i+1);
461  const auto p2 = traj.LocationAtPoint(i+2);
462 
463  auto x0=p0.X(); auto x1=p1.X(); auto x2=p2.X();
464  auto y0=p0.Y(); auto y1=p1.Y(); auto y2=p2.Y();
465  // auto z0=p0.Z(); auto z1=p1.Z(); auto z2=p2.Z();*/
466 
467 recob::Hit h0=hits2d.at(i);
468 recob::Hit h1=hits2d.at(i+1);
469 recob::Hit h2=hits2d.at(i+2);
470 
471 //std::cout << " PeakTime " << h0.PeakTime() << std::endl;
472 
473 float x0=h0.WireID().Wire*3; auto y0=h0.PeakTime()*0.622;
474 float x1=h1.WireID().Wire*3; auto y1=h1.PeakTime()*0.622;
475 float x2=h2.WireID().Wire*3; auto y2=h2.PeakTime()*0.622;
476 
477 //std::cout << " x0 " << x0 << " x1 " << x1 << " x2 " << x2 << std::endl;
478 //std::cout << " y0 " << y0 << " y1 " << y1 << " y2 " << y2 << std::endl;
479  double ym=y0+(x1-x0)/(x2-x0)*(y2-y0);
480  if(abs(x2-x0)<0.001)
481  return -999;
482  if(abs(x2-x1)<0.001)
483  return -999;
484  if(abs(x1-x0)<0.001)
485  return -999;
486  double K=(x1-x0)/(x2-x0);
487 
488  alfa=1/sqrt(1+K*K+(1-K)*(1-K));
489  double dy=y1-ym;
490  double res=dy*(alfa);
491 
492  //cout << " y1 " << y1 << " ym " << ym << " alfa " << alfa << endl;
493  // cout << " deltafit residual " << res << endl;
494  return res;
495 }
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
const TrajectoryMCSFitterICARUS::ScanResult TrajectoryMCSFitterICARUS::doLikelihoodScan ( std::vector< float > &  dtheta,
std::vector< float > &  seg_nradlengths,
std::vector< float > &  cumLen,
bool  fwdFit,
bool  momDepConst,
int  pid 
) const

Definition at line 151 of file TrajectoryMCSFitterICARUS.cxx.

151  {
152  int best_idx = -1;
153  double best_logL = std::numeric_limits<double>::max();
154  double best_p = -1.0;
155  std::vector<float> vlogL;
156  for (double p_test = pMin_; p_test <= pMax_; p_test+=pStep_) {
157  double logL = mcsLikelihood(p_test, angResol_, dtheta, seg_nradlengths, cumLen, fwdFit, momDepConst, pid);
158  if (logL < best_logL) {
159  best_p = p_test;
160  best_logL = logL;
161  best_idx = vlogL.size();
162  }
163 // std::cout << " ptest " << p_test << " likeli " << logL << " bestp "<< best_p << std::endl;
164 //compute likelihood for MC momentum
165  /* double logL = mcsLikelihood(2., angResol_, dtheta, seg_nradlengths, cumLen, fwdFit, momDepConst, pid);
166  if (logL < best_logL) {
167  best_p = 2.;
168  best_logL = logL;
169  best_idx = vlogL.size();
170  }*/
171  vlogL.push_back(logL);
172  }
173  //
174  //uncertainty from left side scan
175  double lunc = -1.0;
176  if (best_idx>0) {
177  for (int j=best_idx-1;j>=0;j--) {
178  double dLL = vlogL[j]-vlogL[best_idx];
179  if ( dLL<0.5 ) {
180  lunc = (best_idx-j)*pStep_;
181  } else break;
182  }
183  }
184  //uncertainty from right side scan
185  double runc = -1.0;
186  if (best_idx<int(vlogL.size()-1)) {
187  for (unsigned int j=best_idx+1;j<vlogL.size();j++) {
188  double dLL = vlogL[j]-vlogL[best_idx];
189  if ( dLL<0.5 ) {
190  runc = (j-best_idx)*pStep_;
191  } else break;
192  }
193  }
194  return ScanResult(best_p, std::max(lunc,runc), best_logL);
195 }
trkf::TrajectoryMCSFitterICARUS::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: TrajectoryMCSFitterICARUS.cxx:262
double TrajectoryMCSFitterICARUS::energyLossBetheBloch ( const double  mass,
const double  e2 
) const

Definition at line 337 of file TrajectoryMCSFitterICARUS.cxx.

337  {
338  // stolen, mostly, from GFMaterialEffects.
339  constexpr double Iinv = 1./188.E-6;
340  constexpr double matConst = 1.4*18./40.;//density*Z/A
341  constexpr double me = 0.511;
342  constexpr double kappa = 0.307075;
343  //
344  const double beta2 = (e2-mass*mass)/e2;
345  const double gamma2 = 1./(1.0 - beta2);
346  const double massRatio = me/mass;
347  const double argument = (2.*me*gamma2*beta2*Iinv) * std::sqrt(1+2*std::sqrt(gamma2)*massRatio + massRatio*massRatio);
348  //
349  double dedx = kappa*matConst/beta2;
350  //
351  if (mass==0.0) return(0.0);
352  if (argument <= exp(beta2)) {
353  dedx = 0.;
354  } else{
355  dedx *= (log(argument)-beta2)*1.E-3; // Bethe-Bloch, converted to GeV/cm
356  if (dedx<0.) dedx = 0.;
357  }
358  return dedx;
359 }
E
process_name E
Definition: rawtorootconvert.fcl:26
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
double TrajectoryMCSFitterICARUS::energyLossLandau ( const double  mass2,
const double  E2,
const double  x 
) const

Definition at line 319 of file TrajectoryMCSFitterICARUS.cxx.

319  {
320  //
321  // eq. (33.11) in http://pdg.lbl.gov/2016/reviews/rpp2016-rev-passage-particles-matter.pdf (except density correction is ignored)
322  //
323  if (x<=0.) return 0.;
324  constexpr double Iinv2 = 1./(188.E-6*188.E-6);
325  constexpr double matConst = 1.4*18./40.;//density*Z/A
326  constexpr double me = 0.511;
327  constexpr double kappa = 0.307075;
328  constexpr double j = 0.200;
329  //
330  const double beta2 = (e2-mass2)/e2;
331  const double gamma2 = 1./(1.0 - beta2);
332  const double epsilon = 0.5*kappa*x*matConst/beta2;
333  //
334  return 0.001*epsilon*( log(2.*me*beta2*gamma2*epsilon*Iinv2) + j - beta2 );
335 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
void TrajectoryMCSFitterICARUS::findSegmentBarycenter ( const recob::TrackTrajectory &  traj,
const size_t  firstPoint,
const size_t  lastPoint,
recob::tracking::Vector_t &  pcdir 
) const

Definition at line 196 of file TrajectoryMCSFitterICARUS.cxx.

196  {
197  int npoints = 0;
198  geo::vect::MiddlePointAccumulator middlePointCalc;
199  size_t nextValid = firstPoint;
200  while (nextValid<lastPoint) {
201  middlePointCalc.add(traj.LocationAtPoint(nextValid));
202  nextValid = traj.NextValidPoint(nextValid+1);
203  npoints++;
204  }
205  const auto avgpos = middlePointCalc.middlePoint();
206  bary=avgpos;
207  //std::cout << " avgpos " << avgpos << std::endl;
208 }
geo::vect::MiddlePointAccumulatorDim
Helper class to compute the middle point in a point set.
Definition: geo_vectors_utils.h:1246
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
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
recob::MCSFitResult trkf::TrajectoryMCSFitterICARUS::fitMcs ( const recob::TrackTrajectory &  traj,
bool  momDepConst = true 
) const
inline

Definition at line 110 of file TrajectoryMCSFitterICARUS.h.

110 { return fitMcs(traj,pIdHyp_,momDepConst); }
trkf::TrajectoryMCSFitterICARUS::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: TrajectoryMCSFitterICARUS.h:110
recob::MCSFitResult trkf::TrajectoryMCSFitterICARUS::fitMcs ( const recob::Track &  track,
bool  momDepConst = true 
) const
inline

Definition at line 111 of file TrajectoryMCSFitterICARUS.h.

111 { return fitMcs(track,pIdHyp_,momDepConst); }
trkf::TrajectoryMCSFitterICARUS::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: TrajectoryMCSFitterICARUS.h:110
recob::MCSFitResult trkf::TrajectoryMCSFitterICARUS::fitMcs ( const recob::Trajectory &  traj,
bool  momDepConst = true 
) const
inline

Definition at line 112 of file TrajectoryMCSFitterICARUS.h.

112 { return fitMcs(traj,pIdHyp_,momDepConst); }
trkf::TrajectoryMCSFitterICARUS::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: TrajectoryMCSFitterICARUS.h:110
recob::MCSFitResult TrajectoryMCSFitterICARUS::fitMcs ( const recob::TrackTrajectory &  traj,
int  pid,
bool  momDepConst = true 
) const

Definition at line 15 of file TrajectoryMCSFitterICARUS.cxx.

15  {
16 
17  //std::cout << " traj nhits " << traj.NHits() << std::endl;
18  //std::cout << " traj lenght " << traj.Length() << std::endl;
19  GetOptimalSegLen(1000,traj.NPoints(),2,traj.Length());
20 
21  //std::cout << " D3p " << d3p << std::endl;
22  //
23  // Break the trajectory in segments of length approximately equal to segLen_
24  //
25  vector<size_t> breakpoints;
26  vector<float> segradlengths;
27  vector<float> cumseglens;
28  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
29 
30  //std::cout << " n segments " << segradlengths.size() << std::endl;
31  //
32  // Fit segment directions, and get 3D angles between them
33  //resi
34  if (segradlengths.size()<2) return recob::MCSFitResult();
35  vector<float> dtheta;
36  Vector_t pcdir0;
37  Vector_t pcdir1;
38  for (unsigned int p = 0; p<segradlengths.size(); p++) {
39  linearRegression(traj, breakpoints[p], breakpoints[p+1], pcdir1);
40  if (p>0) {
41  if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
42  dtheta.push_back(-999.);
43  } else {
44  const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
45  //assert(std::abs(cosval)<=1);
46  //units are mrad
47  double dt = 1000.*acos(cosval);//should we try to use expansion for small angles?
48  dtheta.push_back(dt);
49  //std::cout << " linearfit angle " << dt << std::endl;
50  }
51  }
52  pcdir0 = pcdir1;
53  }
54 
55  vector<Vector_t> barycenters;
56  Vector_t bary;
57 vector<float> dthetaPoly;
58 for (unsigned int p = 0; p<segradlengths.size(); p++) {
59  findSegmentBarycenter(traj, breakpoints[p], breakpoints[p+1], bary);
60  barycenters.push_back(bary);
61 }
62 for (unsigned int p = 2; p<segradlengths.size(); p++) {
63  if (segradlengths[p]<-100. || segradlengths[p-1]<-100. || segradlengths[p-2]<-100.) {
64  dtheta.push_back(-999.);
65  } else {
66  Vector_t dbcp=barycenters[p]-barycenters[p-1];
67  float norm=sqrt(dbcp.X()*dbcp.X()+dbcp.Y()*dbcp.Y()+dbcp.Z()*dbcp.Z());
68  dbcp/=norm;
69  //std::cout << " dbcp " << dbcp << std::endl;
70  Vector_t dbcm=barycenters[p-1]-barycenters[p-2];
71  norm=sqrt(dbcm.X()*dbcm.X()+dbcm.Y()*dbcm.Y()+dbcm.Z()*dbcm.Z());
72 dbcm/=norm;
73  //std::cout << " dbcm " << dbcm << std::endl;
74  const double cosval = dbcp.X()*dbcm.X()+dbcp.Y()*dbcm.Y()+dbcp.Z()*dbcm.Z();
75  //assert(std::abs(cosval)<=1);
76  //units are mrad
77  //std::cout << " cosval " << cosval << std::endl;
78  double dt = 1000.*acos(cosval);//should we try to use expansion for small angles?
79  dthetaPoly.push_back(dt);
80  //std::cout << " polygonal angle " << dt << std::endl;
81 
82  }
83  }
84  //
85  // Perform likelihood scan in forward and backward directions
86  //
87 
88  // std::ofstream outDtheta("dtheta2gev.out",ios::app);
89  //outDtheta << dtheta[0] << std::endl;
90  vector<float> cumLenFwd;
91  vector<float> cumLenBwd;
92  for (unsigned int i = 0; i<cumseglens.size()-2; i++) {
93  //std::cout << " seglen " << cumseglens[i] << std::endl;
94  cumLenFwd.push_back(cumseglens[i]);
95  cumLenBwd.push_back(cumseglens.back()-cumseglens[i+2]);
96  }
97  const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true, momDepConst, pid);
98  const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid);
99  //std::cout << " fwdResult " << fwdResult.p << " bwdResult " << bwdResult.p << std::endl;
100  //const ScanResult fwdResultPoly = doLikelihoodScan(dthetaPoly, segradlengths, cumLenFwd, true, momDepConst, pid);
101  //const ScanResult bwdResultPoly = doLikelihoodScan(dthetaPoly, segradlengths, cumLenBwd, false, momDepConst, pid);
102  //std::cout << " fwdResultPoly " << fwdResultPoly.p << " bwdResultPoly " << bwdResultPoly.p << std::endl;
103 //
104 //for(unsigned int j=0;j<segradlengths.size();j++)
105  //std::cout << " dtheta " << dtheta.at(j) << std::endl;
106  return recob::MCSFitResult(pid,
107  fwdResult.p,fwdResult.pUnc,fwdResult.logL,
108  bwdResult.p,bwdResult.pUnc,bwdResult.logL,
109  segradlengths,dtheta);
110 }
trkf::TrajectoryMCSFitterICARUS::linearRegression
void linearRegression(const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) const
Definition: TrajectoryMCSFitterICARUS.cxx:210
trkf::Vector_t
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
trkf::TrajectoryMCSFitterICARUS::breakTrajInSegments
void breakTrajInSegments(const recob::TrackTrajectory &traj, std::vector< size_t > &breakpoints, std::vector< float > &segradlengths, std::vector< float > &cumseglens) const
Definition: TrajectoryMCSFitterICARUS.cxx:112
p
pdgs p
Definition: selectors.fcl:22
trkf::TrajectoryMCSFitterICARUS::findSegmentBarycenter
void findSegmentBarycenter(const recob::TrackTrajectory &traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t &pcdir) const
Definition: TrajectoryMCSFitterICARUS.cxx:196
recob::TrackTrajectory::Length
double Length(size_t startAt=0) const
Returns the approximate length of the trajectory.
Definition: TrackTrajectory.cxx:70
recob::Trajectory::NPoints
size_t NPoints() const
Returns the number of stored trajectory points.
Definition: Trajectory.h:167
trkf::TrajectoryMCSFitterICARUS::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: TrajectoryMCSFitterICARUS.cxx:151
trkf::TrajectoryMCSFitterICARUS::GetOptimalSegLen
double GetOptimalSegLen(const double guess_p, const int n_points, const int plane, const double length_travelled) const
Definition: TrajectoryMCSFitterICARUS.cxx:383
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:1202
recob::MCSFitResult
Class storing the result of the Maximum Likelihood fit of Multiple Coulomb Scattering angles between ...
Definition: MCSFitResult.h:19
recob::MCSFitResult trkf::TrajectoryMCSFitterICARUS::fitMcs ( const recob::Track &  track,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 115 of file TrajectoryMCSFitterICARUS.h.

115 { 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::TrajectoryMCSFitterICARUS::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: TrajectoryMCSFitterICARUS.h:110
recob::MCSFitResult trkf::TrajectoryMCSFitterICARUS::fitMcs ( const recob::Trajectory &  traj,
int  pid,
bool  momDepConst = true 
) const
inline

Definition at line 116 of file TrajectoryMCSFitterICARUS.h.

116  {
117  recob::TrackTrajectory::Flags_t flags(traj.NPoints());
118  const recob::TrackTrajectory tt(traj,std::move(flags));
119  return fitMcs(tt,pid,momDepConst);
120  }
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::TrajectoryMCSFitterICARUS::fitMcs
recob::MCSFitResult fitMcs(const recob::TrackTrajectory &traj, bool momDepConst=true) const
Definition: TrajectoryMCSFitterICARUS.h:110
double TrajectoryMCSFitterICARUS::GetE ( const double  initial_E,
const double  length_travelled,
const double  mass 
) const

Definition at line 361 of file TrajectoryMCSFitterICARUS.cxx.

361  {
362  //
363  const double step_size = length_travelled / nElossSteps_;
364  //
365  double current_E = initial_E;
366  const double m2 = m*m;
367  //
368  for (auto i = 0; i < nElossSteps_; ++i) {
369  if (eLossMode_==2) {
370  double dedx = energyLossBetheBloch(m,current_E);
371  current_E -= (dedx * step_size);
372  } else {
373  // MPV of Landau energy loss distribution
374  current_E -= energyLossLandau(m2,current_E*current_E,step_size);
375  }
376  if ( current_E <= m ) {
377  // std::cout<<"WARNING: current_E less than mu mass. it is "<<current_E<<std::endl;
378  return 0.;
379  }
380  }
381  return current_E;
382 }
trkf::TrajectoryMCSFitterICARUS::energyLossBetheBloch
double energyLossBetheBloch(const double mass, const double e2) const
Definition: TrajectoryMCSFitterICARUS.cxx:337
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::TrajectoryMCSFitterICARUS::energyLossLandau
double energyLossLandau(const double mass2, const double E2, const double x) const
Definition: TrajectoryMCSFitterICARUS.cxx:319
double TrajectoryMCSFitterICARUS::GetOptimalSegLen ( const double  guess_p,
const int  n_points,
const int  plane,
const double  length_travelled 
) const

Definition at line 383 of file TrajectoryMCSFitterICARUS.cxx.

383  {
384  //
385 // check units of measurment! (energy, length...)
386 double initial_p=guess_p*-0.211*length_travelled;
387 double MIN_P=1000;
388 if(initial_p<MIN_P) initial_p=MIN_P;
389 
390 double sigma=d3p; //to be replaced by computed D3P when ready!
391 
392 //double cosa=1;
393 double sinb=1; //both to be replaced by angles when ready!
394 double DsTot=length_travelled; //to be checked!
395 double p0=13.6; //MeV! constant in MCS formula
396 double LogTerm=0.038; //log term in MCS formula
397 double X0=140.; //interaction length in mm
398 
399 double xi=(sigma/DsTot)*(initial_p/2.)/p0*sqrt(X0/length_travelled)/sqrt(double(n_points))*sinb;
400 //double K=6*xi*xi;
401 double m=1/sqrt(6*xi);
402 
403 //double sigmaMS=p0/(initial_p/2.)*sqrt(length_travelled/X0/m/cosa)/sqrt(2.)/cosa;
404 //double sigmaMeas=pow(m,1.5)*sqrt(double(6.))*sigma*sinb/cosa/DsTot/sqrt(double(n_points));
405 
406 double alfa=1+LogTerm*log(DsTot/sinb/double(m)/X0);
407 
408 double m_corr=(m==1)?m:sqrt(m/(m-1));
409 
410 double xiCorr=sigma/DsTot*((initial_p/2.)/p0)*sqrt(X0/length_travelled)*sqrt((double)n_points)*sinb/alfa;
411 
412 double m2=m_corr/sqrt(6.*xiCorr);
413 if(m==1) m2=10;
414 m2+=0.5;
415 
416 return m2;
417 
418 }
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
X0
standard_singlep gaussian distribution X0
Definition: multigen.fcl:8
void TrajectoryMCSFitterICARUS::linearRegression ( const recob::TrackTrajectory &  traj,
const size_t  firstPoint,
const size_t  lastPoint,
recob::tracking::Vector_t &  pcdir 
) const

Definition at line 210 of file TrajectoryMCSFitterICARUS.cxx.

210  {
211  //
212  int npoints = 0;
213  geo::vect::MiddlePointAccumulator middlePointCalc;
214  size_t nextValid = firstPoint;
215  while (nextValid<lastPoint) {
216  middlePointCalc.add(traj.LocationAtPoint(nextValid));
217  nextValid = traj.NextValidPoint(nextValid+1);
218  npoints++;
219  }
220  const auto avgpos = middlePointCalc.middlePoint();
221  const double norm = 1./double(npoints);
222  //
223  //assert(npoints>0);
224  //
225  TMatrixDSym m(3);
226  nextValid = firstPoint;
227  while (nextValid<lastPoint) {
228  const auto p = traj.LocationAtPoint(nextValid);
229  const double xxw0 = p.X()-avgpos.X();
230  const double yyw0 = p.Y()-avgpos.Y();
231  const double zzw0 = p.Z()-avgpos.Z();
232  m(0, 0) += xxw0*xxw0*norm;
233  m(0, 1) += xxw0*yyw0*norm;
234  m(0, 2) += xxw0*zzw0*norm;
235  m(1, 0) += yyw0*xxw0*norm;
236  m(1, 1) += yyw0*yyw0*norm;
237  m(1, 2) += yyw0*zzw0*norm;
238  m(2, 0) += zzw0*xxw0*norm;
239  m(2, 1) += zzw0*yyw0*norm;
240  m(2, 2) += zzw0*zzw0*norm;
241  nextValid = traj.NextValidPoint(nextValid+1);
242  }
243  //
244  const TMatrixDSymEigen me(m);
245  const auto& eigenval = me.GetEigenValues();
246  const auto& eigenvec = me.GetEigenVectors();
247  //
248  int maxevalidx = 0;
249  double maxeval = eigenval(0);
250  for (int i=1; i<3; ++i) {
251  if (eigenval(i)>maxeval) {
252  maxevalidx = i;
253  maxeval = eigenval(i);
254  }
255  }
256  //
257  pcdir = Vector_t(eigenvec(0, maxevalidx),eigenvec(1, maxevalidx),eigenvec(2, maxevalidx));
258  if (traj.DirectionAtPoint(firstPoint).Dot(pcdir)<0.) pcdir*=-1.;
259  //
260 }
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::TrajectoryMCSFitterICARUS::mass ( int  pid) const
inline

Definition at line 146 of file TrajectoryMCSFitterICARUS.h.

146  {
147  if (abs(pid)==13) { return mumass; }
148  if (abs(pid)==211) { return pimass; }
149  if (abs(pid)==321) { return kmass; }
150  if (abs(pid)==2212) { return pmass; }
151  return util::kBogusD;
152  }
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 TrajectoryMCSFitterICARUS::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 262 of file TrajectoryMCSFitterICARUS.cxx.

262  {
263  //
264  const int beg = (fwd ? 0 : (dthetaij.size()-1));
265  const int end = (fwd ? dthetaij.size() : -1);
266  const int incr = (fwd ? +1 : -1);
267  //
268 // bool print;
269 // if(p>1.29&&p<1.32) print = true;//(p>1.999 && p<2.001);
270 // else print=false;
271  //
272  const double m = mass(pid);
273  const double m2 = m*m;
274  const double Etot = sqrt(p*p + m2);//Initial energy
275  double Eij2 = 0.;
276  //
277  double const fixedterm = 0.5 * std::log( 2.0 * M_PI );
278  double result = 0;
279  for (int i = beg; i != end; i+=incr ) {
280  if (dthetaij[i]<0) {
281  //cout << "skip segment with too few points" << endl;
282  continue;
283  }
284  //
285  if (eLossMode_==1) {
286  // ELoss mode: MIP (constant)
287  constexpr double kcal = 0.002105;
288  const double Eij = Etot - kcal*cumLen[i];//energy at this segment
289  Eij2 = Eij*Eij;
290  } else {
291  // Non constant energy loss distribution
292  const double Eij = GetE(Etot,cumLen[i],m);
293  Eij2 = Eij*Eij;
294  }
295  //
296  if ( Eij2 <= m2 ) {
297  result = std::numeric_limits<double>::max();
298 
299  break;
300  }
301  const double pij = sqrt(Eij2 - m2);//momentum at this segment
302  const double beta = sqrt( 1. - ((m2)/(pij*pij + m2)) );
303  constexpr double tuned_HL_term1 = 11.0038; // https://arxiv.org/abs/1703.06187
304  constexpr double HL_term2 = 0.038;
305  const double tH0 = ( (momDepConst ? MomentumDependentConstant(pij) : tuned_HL_term1) / (pij*beta) ) * ( 1.0 + HL_term2 * std::log( seg_nradl[i] ) ) * sqrt( seg_nradl[i] );
306  const double rms = sqrt( 2.0*( tH0 * tH0 + theta0x * theta0x ) );
307  if (rms==0.0) {
308  //std::cout << " Error : RMS cannot be zero ! " << std::endl;
309  return std::numeric_limits<double>::max();
310  }
311  const double arg = dthetaij[i]/rms;
312  result += ( std::log( rms ) + 0.5 * arg * arg + fixedterm);
313 // if (print && fwd==true) cout << "TrajectoryMCSFitterICARUS 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] << " result " << result << endl;
314  }
315  //std::cout << " momentum " << p <<" likelihood " << result << std::endl;
316  return result;
317 }
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
trkf::TrajectoryMCSFitterICARUS::GetE
double GetE(const double initial_E, const double length_travelled, const double mass) const
Definition: TrajectoryMCSFitterICARUS.cxx:361
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::TrajectoryMCSFitterICARUS::MomentumDependentConstant
double MomentumDependentConstant(const double p) const
Definition: TrajectoryMCSFitterICARUS.h:140
double trkf::TrajectoryMCSFitterICARUS::MomentumDependentConstant ( const double  p) const
inline

Definition at line 140 of file TrajectoryMCSFitterICARUS.h.

140  {
141  //these are from https://arxiv.org/abs/1703.06187
142  constexpr double a = 0.1049;
143  constexpr double c = 11.0038;
144  return (a/(p*p)) + c;
145  }
p
pdgs p
Definition: selectors.fcl:22
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
void trkf::TrajectoryMCSFitterICARUS::set2DHits ( std::vector< recob::Hit >  h)
inline

Definition at line 157 of file TrajectoryMCSFitterICARUS.h.

157 {hits2d=h;}
h
while getopts h
Definition: getStatistics.sh:13

Member Data Documentation

double trkf::TrajectoryMCSFitterICARUS::angResol_
private

Definition at line 170 of file TrajectoryMCSFitterICARUS.h.

float trkf::TrajectoryMCSFitterICARUS::d3p
private

Definition at line 173 of file TrajectoryMCSFitterICARUS.h.

int trkf::TrajectoryMCSFitterICARUS::eLossMode_
private

Definition at line 166 of file TrajectoryMCSFitterICARUS.h.

std::vector<recob::Hit> trkf::TrajectoryMCSFitterICARUS::hits2d
private

Definition at line 172 of file TrajectoryMCSFitterICARUS.h.

int trkf::TrajectoryMCSFitterICARUS::minHitsPerSegment_
private

Definition at line 164 of file TrajectoryMCSFitterICARUS.h.

int trkf::TrajectoryMCSFitterICARUS::minNSegs_
private

Definition at line 162 of file TrajectoryMCSFitterICARUS.h.

int trkf::TrajectoryMCSFitterICARUS::nElossSteps_
private

Definition at line 165 of file TrajectoryMCSFitterICARUS.h.

int trkf::TrajectoryMCSFitterICARUS::pIdHyp_
private

Definition at line 161 of file TrajectoryMCSFitterICARUS.h.

double trkf::TrajectoryMCSFitterICARUS::pMax_
private

Definition at line 168 of file TrajectoryMCSFitterICARUS.h.

double trkf::TrajectoryMCSFitterICARUS::pMin_
private

Definition at line 167 of file TrajectoryMCSFitterICARUS.h.

double trkf::TrajectoryMCSFitterICARUS::pStep_
private

Definition at line 169 of file TrajectoryMCSFitterICARUS.h.

double trkf::TrajectoryMCSFitterICARUS::segLen_
private

Definition at line 163 of file TrajectoryMCSFitterICARUS.h.

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