TimingFitAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file TimingFitAlg.cxx
3 // \brief Algorithm to determine the direction of tracks using timing
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
10 #include "RecoBase/RecoHit.h"
11 #include "RecoBase/Track.h"
12 #include "Utilities/func/MathUtil.h"
13 
15 
16 #include <cmath>
17 
18 #include "TError.h"
19 #include "TGraphErrors.h"
20 #include "TH2.h"
21 #include "TPad.h"
22 
23 #include "Minuit2/FCNBase.h"
24 #include "Minuit2/FunctionMinimum.h"
25 #include "Minuit2/MnMigrad.h"
26 #include "Minuit2/MnUserParameters.h"
27 
28 namespace tf
29 {
30  //......................................................................
32  {
33  }
34 
35  class HoughScore: public ROOT::Minuit2::FCNBase
36  {
37  public:
38  HoughScore(const std::vector<Pt>& pts) : fPts(pts) {}
39 
40  double Eval(double t0, double cinv) const
41  {
42  double ret = 0;
43  const unsigned int I = fPts.size();
44  for(unsigned int i = 0; i < I; ++i){
45  const Pt& pt = fPts[i];
46 
47  const double dt = pt.t-(t0+pt.L*cinv);
48  const double w = exp(-util::sqr(dt/pt.err)/2)/pt.err;
49  ret += w;
50  }
51  return ret;
52  }
53 
54  virtual double operator()(const std::vector<double>& xs) const
55  {
56  // Minuit is trying to minimize, so use a minus sign to maximize score
57  return -Eval(xs[0], xs[1]);
58  }
59 
60  virtual double Up() const {return 1;}
61  protected:
62  const std::vector<Pt>& fPts;
63  };
64 
65  ErrorHandlerFunc_t gOldHandler = 0;
66 
67  //......................................................................
68  // For some reason gErrorIgnoreLevel isn't quite working, so we have to bust
69  // out the nuclear option.
70  void TimingFitErrorHandler(int level, bool b_abort,
71  const char* loc, const char* msg)
72  {
73  // If it's a serious message, do what we would normally do with it
74  if(level > kInfo) (*gOldHandler)(level, b_abort, loc, msg);
75  // Otherwise surppress it
76  }
77 
78  //......................................................................
79  TimingFitResult TimingFitAlg::HoughFitPts(const std::vector<Pt>& pts) const
80  {
82  if(pts.size() < 2){
83  ret.valid = false;
84  return ret;
85  }
86 
87  gOldHandler = SetErrorHandler(TimingFitErrorHandler);
88 
89  HoughScore hs(pts);
90 
91  ROOT::Minuit2::MnUserParameters mnPars;
92  mnPars.Add("t0", 0, 100);
93  mnPars.Add("cinv", 0, 1./30);
94 
95  // Best fit
96  ROOT::Minuit2::MnMigrad mnMigrad(hs, mnPars);
97  ROOT::Minuit2::FunctionMinimum funcMin = mnMigrad();
98  ret.bestInvSpeed = mnMigrad.Parameter(1).Value();
99  ret.bestScore = -funcMin.Fval()/pts.size();
100  ret.valid = funcMin.IsValid();
101 
102  // Degugging code to look at fit failures
103  /*
104  if(!ret.valid){
105  TH2F h("", ";t_{0} (ns);1/c (ns/cm)", 100, -300, 300, 100, -.2, +.2);
106 
107  for(int x = 0; x < h.GetNbinsX()+2; ++x){
108  const double t0 = h.GetXaxis()->GetBinCenter(x);
109  for(int y = 0; y < h.GetNbinsY()+2; ++y){
110  const double cinv = h.GetYaxis()->GetBinCenter(y);
111 
112  h.SetBinContent(x, y, hs.Eval(t0, cinv));
113  } // end for y (cinv)
114  } // end for x (t0)
115 
116  h.Draw("colz");
117  gPad->Print("map.eps");
118 
119  TGraphErrors* g = new TGraphErrors;
120  for(unsigned int i = 0; i < pts.size(); ++i){
121  g->SetPoint(g->GetN(), pts[i].L, pts[i].t);
122  g->SetPointError(g->GetN()-1, 0, pts[i].err);
123  }
124 
125  g->SetMarkerStyle(kFullCircle);
126  g->Draw("ap");
127  gPad->Print("t_vs_z.eps");
128 
129  exit(0);
130  }
131  */
132 
133  // Best fit constrained to +c
134  mnMigrad.SetValue(0u, 0.);
135  mnMigrad.SetValue(1, +1./30);
136  mnMigrad.Fix(1);
137  funcMin = mnMigrad();
138  ret.fwdScore = -funcMin.Fval()/pts.size();
139 
140  // Best fit constrained to -c
141  mnMigrad.SetValue(0u, 0.);
142  mnMigrad.SetValue(1, -1./30);
143  mnMigrad.Fix(1);
144  funcMin = mnMigrad();
145  ret.revScore = -funcMin.Fval()/pts.size();
146 
147  SetErrorHandler(gOldHandler);
148 
149  return ret;
150  }
151 
152  //......................................................................
154  {
155  std::vector<Pt> pts;
156 
158 
159  const unsigned int cellMax = track->NCell();
160  const TVector3 dir = track->Dir();
161  // Center Pts around (0,0). Makes natural choice of MINUIT seeding
162  const TVector3 origin = track->MeanXYZ();
163  const double t0 = track->MeanTNS();
164  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
165  const art::Ptr<rb::CellHit> chit = track->Cell(cellIdx);
166 
167  const rb::RecoHit rhit = track->RecoHit(chit);
168  const double t = rhit.T();
169 
170  const TVector3 pos(rhit.X(), rhit.Y(), rhit.Z());
171 
172  Pt pt((pos-origin).Dot(dir), t-t0, cal->GetTimeRes(*chit));
173  pts.push_back(pt);
174  } // end for cellIdx
175 
176  return HoughFitPts(pts);
177  }
178 
179 } // namespace
float T() const
Definition: RecoHit.h:39
void TimingFitErrorHandler(int level, bool b_abort, const char *loc, const char *msg)
double GetTimeRes(rb::CellHit const &cellhit)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
Represent one distance/time pair on a track.
Definition: TimingFitAlg.h:19
ErrorHandlerFunc_t gOldHandler
double Eval(double t0, double cinv) const
float Dot(const Proxy &v) const
Store results of TimingFitAlg.
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
Timing fit.
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
TimingFitResult HoughFit(const rb::Track *track) const
HoughScore(const std::vector< Pt > &pts)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual double Up() const
Definition: Cand.cxx:23
TimingFitResult HoughFitPts(const std::vector< Pt > &pts) const
Helper for HoughFit. May be useful for more advanced users.
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
double L
Definition: TimingFitAlg.h:22
virtual double operator()(const std::vector< double > &xs) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
TDirectory * dir
Definition: macro.C:5
bool valid
As reported by MINUIT.
const int cellMax
float X() const
Definition: RecoHit.h:36
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
double t
Definition: TimingFitAlg.h:22
float Y() const
Definition: RecoHit.h:37
double err
Definition: TimingFitAlg.h:22
double bestInvSpeed
ns/cm
Float_t w
Definition: plot.C:20
const std::vector< Pt > & fPts