LCFIVertex  0.7.2
vertexfitterlsm.cpp
1 /*
2  Simple chi2 vertex fitter for zvtop.
3 
4  03/Nov/05
5  Currently no use of errors at all, also error handling needs to be improved.
6  Could do with some doxygen comments for improved documentation too.
7  01/06/06
8  Error included from FORTRAN for now- B Jeffery
9 
10  mark.grimes@bristol.ac.uk
11 */
12 
13 #include <map>
14 #include "../include/vertexfitterlsm.h"
15 
16 #include "../include/maxminfinder.h"
17 #include "../include/candidatevertex.h"
18 #include "../../inc/trackstate.h"
19 #include "../../inc/track.h"
20 #include "../include/interactionpoint.h"
21 #include "../../util/inc/matrix.h"
22 #include "../../util/inc/vector3.h"
23 
24 namespace vertex_lcfi { namespace ZVTOP
25 {
26  VertexFitterLSM::VertexFitterLSM()
27  :_UseManualSeed(false),_InitialStep(100.0/1000.0)
28  {
29  }
30  void VertexFitterLSM::fitVertex(const std::vector<TrackState*> & Tracks, InteractionPoint* IP, Vector3 & Result)
31  {
32  _trackStateList=Tracks; // so that valueAt() can access them without having to pass as parameter
33  _ip=IP;
34 
35  //TODO - Is this the optimal seed? - For Vertexes we have added or removed tracks from the old fit maybe a good start....
36  //TODO Throw something if we have <2 objects to fit?
37  //Find seed position, we do a load of 2 prong fits and take the average
38  Vector3 Seed(0,0,0);
39  if (_UseManualSeed)
40  {
41  Seed = _ManualSeed;
42  }
43  else
44  {
45  if (Tracks.size()>1)
46  {
47  int NumSeeds = 0;
48  int N = Tracks.size();
49  for (int OuterIndex=0;OuterIndex < N-1;++OuterIndex)
50  {
51  for (int InnerIndex=OuterIndex+1;InnerIndex < N;++InnerIndex)
52  {
53  Tracks[OuterIndex]->swimToStateNearest(Tracks[InnerIndex]);
54  Tracks[InnerIndex]->swimToStateNearest(Tracks[OuterIndex]->position());
55  Vector3 LeftPoint = Tracks[OuterIndex]->position();
56  Vector3 MidPoint = Tracks[OuterIndex]->position() + ((Tracks[InnerIndex]->position() - Tracks[OuterIndex]->position()) * 0.5);
57  Vector3 RightPoint = Tracks[InnerIndex]->position();
58  //std::cout << "lmr " << LeftPoint << MidPoint << RightPoint << std::endl;
59  double p = (RightPoint.subtract(LeftPoint)).mag();
60  if (p > 0.0001/1000.0)
61  {
62  //std::cout << "p" << RightPoint.subtract(LeftPoint) << std::endl;
63  //std::cout << "CL" << std::endl;
64  double ChiLeft = _chi2Contribution( MidPoint, Tracks[OuterIndex]);
65  //std::cout << "CR" << std::endl;
66  double ChiRight = _chi2Contribution( MidPoint, Tracks[InnerIndex]);
67  //std::cout << "chi " << ChiLeft << " " << ChiRight;
68  double Sigma2Left = ((0.5*p)*(0.5*p))/ChiLeft;
69  double Sigma2Right = ((-0.5*p)*(-0.5*p))/ChiRight;
70  double Frac = Sigma2Left/(Sigma2Left+Sigma2Right);
71  Seed = Seed + (LeftPoint + ((RightPoint-LeftPoint)*Frac));
72  }
73  else
74  Seed = Seed + LeftPoint;
75  ++NumSeeds;
76  }
77  }
78  Seed = Seed/NumSeeds; //Take the average of the 2-prong seeds
79  }
80  if (Tracks.size()==1)
81  {
82  if (IP)
83  {
84  Tracks[0]->swimToStateNearest(IP->position());
85  Vector3 LeftPoint = IP->position();
86  Vector3 MidPoint = IP->position() + ((Tracks[0]->position() - IP->position()) * 0.5);
87  Vector3 RightPoint = Tracks[0]->position();
88  double p = (RightPoint.subtract(LeftPoint)).mag();
89  if (p > 0.0001/1000.0)
90  {
91  double ChiLeft = _chi2Contribution( MidPoint, IP);
92  double ChiRight = _chi2Contribution( MidPoint, Tracks[0]);
93  double Sigma2Left = ((0.5*p)*(0.5*p))/ChiLeft;
94  double Sigma2Right = ((-0.5*p)*(-0.5*p))/ChiRight;
95  double Frac = Sigma2Left/(Sigma2Left+Sigma2Right);
96  Seed = Seed + (LeftPoint + ((RightPoint-LeftPoint)*Frac));
97  }
98  else
99  Seed = Seed + LeftPoint;
100  }
101  }
102  if (Tracks.size()==0)
103  {
104  if (IP)
105  {
106  Seed = IP->position();
107  }
108  }
109  }
110  //TODO Minimisation needed qfor two object fits? Check if seed is close enough
111  //std::cout << "Seed " << Seed <<std::endl;
112  //Minimise the chi squared if we havee more than 2 objects to fit
113  if (Tracks.size()>2 || (Tracks.size()>1 && IP))
114  //if (Tracks.size()>1)
115  {
116  //for (double mu=1;mu<10000;mu=mu+((10000-1)/9))
117  //for (double mu=300;mu<1600;mu=mu+((1600-300)/19))
118  //for (double v=1;v<500;v=v+((500-1)/9))
119  //for (double v=2000;v<6000;v=v+((6000-2000)/19))
120  //Position=minimiser.Minimise(Seed,1400,4000);
121  //Create a minimiser
122  //std::cout << Seed << std::endl;
123  ZVTOP::FunctionMinimiser<ZVTOP::VertexFitterLSM> minimiser( this, _InitialStep, 6 );
124  std::vector<double> result = minimiser.Minimise(Seed.stlVector());//,1,2);
125  Result.x() = result[0];
126  Result.y() = result[1];
127  Result.z() = result[2];
128  //Output distance from seed to position as a function of num tracks and decay length
129  //std::cout << Tracks.size() << " " << Position.mag()*10000 << " " << Position.distanceTo(Seed)*10000 << std::endl;
130  }
131  else
132  Result = Seed;
133  }
134 
135  void VertexFitterLSM::fitVertex(const std::vector<TrackState*> & Tracks, InteractionPoint* IP, Vector3 & Result, double & ChiSquaredOfFit)
136  {
137  this->fitVertex(Tracks,IP,Result);
138  //fill in ChiSquaredOfTrack and of fit
139  ChiSquaredOfFit = 0;
140  for( std::vector<TrackState*>::const_iterator i=Tracks.begin(); i<Tracks.end(); i++ )
141  {
142  ChiSquaredOfFit += _chi2Contribution( Result, (*i) );
143  }
144  //fill in ChiSquaredOfIP if no IP this just gives zero
145  ChiSquaredOfFit += _chi2Contribution( Result, IP );
146 
147  }
148  void VertexFitterLSM::fitVertex(const std::vector<TrackState*> & Tracks, InteractionPoint* IP, Vector3 & Result, double & ChiSquaredOfFit, std::map<TrackState*,double> & ChiSquaredOfTrack,double & ChiSquaredOfIP)
149  {
150 
151  this->fitVertex(Tracks,IP,Result);
152  //fill in ChiSquaredOfTrack and of fit
153  ChiSquaredOfFit = 0;
154  ChiSquaredOfTrack.clear();
155  for( std::vector<TrackState*>::const_iterator i=Tracks.begin(); i<Tracks.end(); i++ )
156  {
157  double chi = _chi2Contribution( Result, (*i) );
158  ChiSquaredOfTrack.insert( std::pair<TrackState*,double>( (*i), chi ) );
159  ChiSquaredOfFit += chi;
160  }
161  //fill in ChiSquaredOfIP if no IP this just gives zero
162  ChiSquaredOfIP=_chi2Contribution( Result, IP );
163  ChiSquaredOfFit += ChiSquaredOfIP;
164  }
165 
166  void VertexFitterLSM::fitVertex(const std::vector<TrackState*> & Tracks, InteractionPoint* IP, Vector3 & Result, Matrix3x3 & ResultError, double & ChiSquaredOfFit, std::map<TrackState*,double> & ChiSquaredOfTrack,double & ChiSquaredOfIP)
167  {
168  this->fitVertex(Tracks,IP,Result,ChiSquaredOfFit,ChiSquaredOfTrack,ChiSquaredOfIP);
169  //set the total of the covariances to zero so that they can be added as we go
170  ResultError.clear();
171  if (Tracks.size()>1)
172  {
173  for( std::vector<TrackState*>::const_iterator i=Tracks.begin(); i<Tracks.end(); i++ )
174  {
175  ResultError += (*i)->vertexErrorContribution(Result);
176  }
177  ResultError = InvertMatrix(ResultError);
178  }
179  else
180  if (IP) ResultError = IP->errorMatrix();
181  }
182 
183  double VertexFitterLSM::valueAt(const Vector3 & point)
184  {
185  //work out the chi2 at the point "point"
186  double chi2=0;
187  //loop over all the tracks
188  for( std::vector<TrackState*>::iterator i=_trackStateList.begin(); i<_trackStateList.end(); i++ )
189  {
190  chi2+=_chi2Contribution( point, (*i) );
191  }
192 
193  //add in the contribution (if any) from the IP if we only have one track!
194  //Line below needed as fortran fitter doesn't use ip?
195  if (_trackStateList.size()<2)
196  chi2+=_chi2Contribution( point, _ip );
197 
198  return chi2;
199  }
200 
201  double VertexFitterLSM::valueAt(const std::vector<double> & point)
202  {
203  return this->valueAt(Vector3(point[0],point[1],point[2]));
204  }
205 
206  void VertexFitterLSM::setSeed(Vector3 Seed)
207  {
208  _UseManualSeed = true;
209  _ManualSeed = Seed;
210  }
211 
212  void VertexFitterLSM::setInitialStep(double Step)
213  {
214  _InitialStep = Step;
215  }
216 
217  double VertexFitterLSM::_chi2Contribution( const Vector3 & point, TrackState* pTrackState )
218  {
219  if( 0==pTrackState ) return 0;
220  else return pTrackState->chi2(point);
221  }
222 
223  double VertexFitterLSM::_chi2Contribution( const Vector3 & point, InteractionPoint* pIP )
224  {
225  if( 0==pIP ) return 0;
226  else return pIP->chi2(point);
227  }
228 }}