LCFIPlus  0.6.5
VertexFitterSimple.h
Go to the documentation of this file.
1 // VertexFitterLCFI.h
2 
3 #ifndef VertexFitterSimple_h
4 #define VertexFitterSimple_h 1
5 
6 #include "lcfiplus.h"
7 #include "geometry.h"
8 
9 using namespace std;
10 using namespace lcfiplus;
11 
12 namespace lcfiplus {
13 
14 template<class Iterator>
16  public:
17  Vertex* operator() (Iterator tracksBegin, Iterator tracksEnd, Vertex* pointConstraint = 0, bool pointInitialOnly = false) {
18  bool verbose = false;
19 
20  GeometryHandler* gh = GeometryHandler::Instance();
21  if (pointConstraint) {
22  Point* ip = new Point(pointConstraint);
23  vector<PointBase*> tracks;
24  if (!pointInitialOnly)
25  tracks.push_back(ip);
26  int ntracks = 0;
27  for (Iterator it = tracksBegin; it != tracksEnd; it++,ntracks++) {
28  tracks.push_back(new Helix(*it));
29  }
30  if (verbose)
31  cout << "VertexFitterSimple: number of tracks is " << ntracks << endl;
32 
33  TVector3 initial = pointConstraint->getPos();
34  Point* result = new Point;
35  double chi2 = -gh->PointFit(tracks, initial, result);
36 
37  TVector3 vresult = result->GetPos();
38  double cov[6];
39  cov[Vertex::xx] = result->GetErr(0,0);
40  cov[Vertex::xy] = result->GetErr(0,1);
41  cov[Vertex::xz] = result->GetErr(0,2);
42  cov[Vertex::yy] = result->GetErr(1,1);
43  cov[Vertex::yz] = result->GetErr(1,2);
44  cov[Vertex::zz] = result->GetErr(2,2);
45 
46  if (verbose) {
47  cout << "Vertex cov matrix:" << endl;
48  cout << scientific << cov[Vertex::xx] << " ";
49  cout << cov[Vertex::yy] << " ";
50  cout << cov[Vertex::zz] << " ";
51  cout << cov[Vertex::xy] << " ";
52  cout << cov[Vertex::yz] << " ";
53  cout << cov[Vertex::xz] << endl << fixed;
54  }
55 
56  if (verbose) {
57  cout << "VertexFitterSimple: vertex position is " << endl;
58  vresult.Print();
59  }
60 
61  Vertex* vtx = new Vertex(chi2, TMath::Prob(chi2, ntracks*2-3), vresult.x(), vresult.y(), vresult.z(), cov, false);
62  for (Iterator it = tracksBegin; it != tracksEnd; it++,ntracks++) {
63  Helix hel(*it);
64  double ll = hel.LogLikelihood(vresult); // need to incorporate vertex error??
65 
66  if (verbose)
67  cout << "VertexFitterSimple: track loglikelihood is " << ll << endl;
68 
69  vtx->add(*it, -ll);
70  }
71  for (vector<PointBase*>::iterator it = tracks.begin(); it != tracks.end(); it++) {
72  delete *it;
73  }
74 
75  delete result;
76  return vtx;
77  }
78 
79  // without point constraint
80  vector<Helix*> tracks;
81  int ntracks = 0;
82  for (Iterator it = tracksBegin; it != tracksEnd; it++,ntracks++) {
83  tracks.push_back(new Helix(*it));
84  }
85  if (verbose)
86  cout << "VertexFitterSimple: number of tracks is " << ntracks << endl;
87 
88  Point* result = new Point;
89  double chi2 = -gh->HelixPointFit(tracks, result);
90 
91  TVector3 vresult = result->GetPos();
92  double cov[6];
93  cov[Vertex::xx] = result->GetErr(0,0);
94  cov[Vertex::xy] = result->GetErr(0,1);
95  cov[Vertex::xz] = result->GetErr(0,2);
96  cov[Vertex::yy] = result->GetErr(1,1);
97  cov[Vertex::yz] = result->GetErr(1,2);
98  cov[Vertex::zz] = result->GetErr(2,2);
99 
100  if (verbose) {
101  cout << "VertexFitterSimple: vertex position is " << endl;
102  vresult.Print();
103  }
104 
105  Vertex* vtx = new Vertex(chi2, (ntracks > 1 ? TMath::Prob(chi2, ntracks*2-3) : 1), vresult.x(), vresult.y(), vresult.z(), cov, false);
106  for (Iterator it = tracksBegin; it != tracksEnd; it++, ntracks++) {
107  Helix hel(*it);
108  double ll = hel.LogLikelihood(vresult); // need to incorporate vertex error??
109  if (verbose)
110  cout << "VertexFitterSimple: track loglikelihood is " << ll << endl;
111 
112  vtx->add(*it, -ll);
113  }
114  for (vector<Helix*>::iterator it = tracks.begin(); it != tracks.end(); it++) {
115  delete *it;
116  }
117 
118  delete result;
119  return vtx;
120  }
121 
122  double getChi2(const Vertex* vtx, const Track* trk, int mode=1) {
123  // 110510 suehara for IPassoc study
124  if (mode == 0) {
125  // mode 0: no fit at all
126  Helix hel(trk);
127  TVector3 v = vtx->getPos();
128  return -hel.LogLikelihood(v);
129  } else if (mode == 1) {
130  // mode 1: vertex treated as errored point
131 
132  Point ptVtx(vtx);
133  Helix hel(trk);
134 
135  vector<PointBase*> vpt;
136  vpt.push_back(&ptVtx);
137  vpt.push_back(&hel);
138 
139  GeometryHandler* gh = GeometryHandler::Instance();
140  return -gh->PointFit(vpt, vtx->getPos(), 0);
141  } else {
142  // mode 2: vertex treated as track list
143  vector<PointBase*> vpt;
144  for (unsigned int n=0; n<vtx->getTracks().size(); n++) {
145  vpt.push_back(new Helix(vtx->getTracks()[n]));
146  }
147  vpt.push_back(new Helix(trk));
148 
149  GeometryHandler* gh = GeometryHandler::Instance();
150  Point ret;
151  double ll = -gh->PointFit(vpt, vtx->getPos(), &ret);
152 
153  for (unsigned int n=0; n<vpt.size(); n++) {
154  delete vpt[n];
155  }
156  if (mode == 2) {
157  return ll;
158  } else {
159  Helix hel(trk);
160  return -hel.LogLikelihood(ret.GetPos());
161  }
162  }
163  }
164 
165 };
166 
167 
170 }
171 
172 #endif
double GetErr(int i, int j) const
Definition: geometry.h:53
Definition: lcfiplus.h:771
Definition: geometry.h:222
double getChi2(const Vertex *vtx, const Track *trk, int mode=1)
Definition: VertexFitterSimple.h:122
Definition: lcfiplus.h:384
Definition: VertexFitterSimple.h:15
Definition: geometry.h:64
TVector3 getPos() const
Definition: lcfiplus.h:817
VertexFitterSimple< vector< const Track * >::const_iterator > VertexFitterSimple_V
Definition: VertexFitterSimple.h:168
double HelixPointFit(const vector< Helix * > &helices, Point *result=0)
VertexFitterSimple< list< const Track * >::const_iterator > VertexFitterSimple_L
Definition: VertexFitterSimple.h:169
Definition: geometry.h:29
TVector3 GetPos() const
Definition: geometry.h:56
void add(const Track *trk)
Definition: lcfiplus.cc:827
double PointFit(const vector< PointBase * > &points, const TVector3 &initial, Point *result=0)
cov
Definition: lcfiplus.h:26
virtual double LogLikelihood(const TVector3 &p) const
Definition: geometry.h:118
const vector< const Track * > & getTracks() const
Definition: lcfiplus.h:823