LCFIVertex  0.7.2
vertexfinderghost.cpp
1 #include "../include/vertexfinderghost.h"
2 
3 #include "../../inc/track.h"
4 #include "../include/candidatevertex.h"
5 #include "../include/interactionpoint.h"
6 #include "../include/vertexfunction.h"
7 #include "../include/vertexfunctionclassic.h"
8 #include "../../inc/trackstate.h"
9 #include "../include/ghostfinderstage1.h"
10 #include "../../util/inc/matrix.h"
11 #include "../../util/inc/memorymanager.h"
12 #include "../../util/inc/util.h"
13 #include <vector>
14 #include <list>
15 #include <ctime>
16 
17 namespace vertex_lcfi { namespace ZVTOP
18 {
19 
20 // Ascending distance from IP func
21 struct IPDistAscending
22 {
23  InteractionPoint* IP;
24 
25  IPDistAscending(InteractionPoint* _IP) : IP(_IP) {}
26 
27  bool operator()(CandidateVertex*& rpStart, CandidateVertex*& rpEnd)
28  {
29  if (IP)
30  {
31  return rpStart->position().distanceTo2(IP->position()) < rpEnd->position().distanceTo2(IP->position());
32  }
33  else
34  {
35  return rpStart->position().distanceTo2(Vector3(0,0,0)) < rpEnd->position().distanceTo2(Vector3(0,0,0));
36  }
37  }
38 };
39 
40 VertexFinderGhost::VertexFinderGhost(const std::vector<Track*> &Tracks, InteractionPoint* IP)
41 :_TrackList(Tracks),_IP(IP)
42 {
43 }
44 
45 
46 void VertexFinderGhost::addTrack(Track* const Track)
47 {
48  _TrackList.push_back(Track);
49 }
50 
51 void VertexFinderGhost::setIP(InteractionPoint* const IP)
52 {
53  _IP=IP;
54 }
55 // returns true if track was in jet and removed
56 bool VertexFinderGhost::removeTrack(Track* const TrackToRemove)
57 {
58  std::vector<Track*>::iterator iTrack;
59 
60  for (iTrack = _TrackList.begin();((*iTrack) != TrackToRemove) && iTrack!=_TrackList.end();++iTrack)
61  ;
62  //iTrack now points to end or the track we want to remove
63  if (iTrack!=_TrackList.end())
64  {
65 
66  _TrackList.erase(iTrack);
67  return 1;
68  }
69  else
70  return 0;
71 }
72 
73 bool VertexFinderGhost::clearIP()
74 {
75  _IP=0;
76  return true;
77 }
78 
79 
80 std::list<CandidateVertex*> VertexFinderGhost::findVertices()
81 {
82  //TODO Momentum factor in fortran?
83  //TODO Cope with no tracks etc.
84  //TODO Check for existance of IP
85  //TODO Create candiates without vf (new constructor)
86  //TODO Change this next line to a parameter
87  using std::cout;using std::endl;clock_t start,pstart;int debug=0;
88  //If we have no tracks return the IP
89  if (_TrackList.empty())
90  {
91  if (debug) std::cout << "Ghost: No Tracks" << std::endl;
92  if(!_IP)
93  return std::list<CandidateVertex*>();
94  else
95  {
96  std::list<CandidateVertex*> ret;
97  CandidateVertex* CV = new CandidateVertex(std::vector<TrackState*>(),_IP,0);
98  MemoryManager<CandidateVertex>::Event()->registerObject(CV);
99  ret.push_back(CV);
100  return ret;
101  }
102  }
103  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "Initial Ghost direction finding and resize...."; cout.flush();pstart=clock();start=clock();}
104  //First we find the ghost
105  Track* GhostTrack = GhostFinderStage1().findGhost(_InitialGhostWidth,_MaxChi2Allowed,_SeedDirection,_TrackList,_IP);
106  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "\tdone!\t\t\t" << ((double)clock()-(double)pstart)*1000.0/CLOCKS_PER_SEC << "ms" << endl; cout.flush();}
107 
108  //Make trackstates of the tracks
109  std::vector<TrackState*> TrackStates;
110  for (std::vector<Track*>::iterator iTrack = _TrackList.begin();iTrack != _TrackList.end();++iTrack)
111  {
112  TrackState* Track = (*iTrack)->makeState();
113  TrackStates.push_back(Track);
114  }
115  //And of the ghost
116  TrackState* GhostTrackState = GhostTrack->makeState();
117 
118  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "Generating 1-prong Candidate Vertices + IP...."; cout.flush();start=clock();}
119  //Then we make a list of all the 1 prong (Ghost + 1 track) vertices
120  std::list<CandidateVertex*> Candidates;
121  //Loop over Jet Tracks
122  for( std::vector<TrackState*>::const_iterator iTrack=TrackStates.begin(); iTrack != TrackStates.end(); ++iTrack)
123  {
124  //Make a CV of this track and the ghost and add it to the list
125  std::vector<TrackState*> Tracks;
126  Tracks.push_back(*iTrack);
127  Tracks.push_back(GhostTrackState);
128  CandidateVertex* CV = new CandidateVertex(Tracks,(InteractionPoint*)0,0);
129  MemoryManager<CandidateVertex>::Event()->registerObject(CV);
130  Candidates.push_back(CV);
131  }
132  //And add a CV with just the IP
133  {
134  std::vector<TrackState*> Tracks;
135  CandidateVertex* CV = new CandidateVertex(Tracks,_IP,0);
136  MemoryManager<CandidateVertex>::Event()->registerObject(CV);
137  Candidates.push_back(CV);
138  }
139  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "\tdone!" << " "<< Candidates.size() << " Candidates" << "\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 << "ms" <<endl; cout.flush();}
140  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug>1) {for (std::list<CandidateVertex*>::iterator iCV = Candidates.begin();iCV != Candidates.end();++iCV) cout << *iCV <<endl;}
141 
142  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "Generating Trial Merged Vertices...."; cout.flush();start=clock();}
143  //We make a list of vertices that are all merges of all possible combinations of the current candidates
144  std::list<CandidateVertex*> TrialMergedCandidates;
145  //And we remember which vertices were merged to make each trial
146  std::map<CandidateVertex*,std::vector<CandidateVertex*> > VerticesContainedIn;
147 
148  for (std::list<CandidateVertex*>::const_iterator iOuterCV = Candidates.begin();iOuterCV != --(Candidates.end());++iOuterCV)
149  {
150  std::list<CandidateVertex*>::const_iterator CopyiOuterCV = iOuterCV;
151  for (std::list<CandidateVertex*>::const_iterator iInnerCV = (++CopyiOuterCV);iInnerCV != Candidates.end();++iInnerCV)
152  {
153  std::vector<CandidateVertex*> ToMerge;
154  ToMerge.push_back(*iOuterCV);
155  ToMerge.push_back(*iInnerCV);
156 
157  CandidateVertex* Merged = new CandidateVertex(ToMerge);
158  MemoryManager<CandidateVertex>::Event()->registerObject(Merged);
159  //If we merged the ghost and ip, just keep the IP
160  if (Merged->hasTrack(GhostTrack) && Merged->interactionPoint())
161  {
162  Merged->removeTrack(GhostTrack);
163  }
164 
165  TrialMergedCandidates.push_back(Merged);
166 
167  VerticesContainedIn[Merged] = ToMerge;
168  }
169  }
170  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "\t\tdone!" << " "<< TrialMergedCandidates.size() << "T " << Candidates.size() << "C Verts"<< "\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 << "ms" <<endl; cout.flush();}
171  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug>1) {for (std::list<CandidateVertex*>::iterator iCV = TrialMergedCandidates.begin();iCV != TrialMergedCandidates.end();++iCV) cout << *iCV <<endl;}
172  do
173  {
174  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "Find most probable trial vertex..."; cout.flush();start=clock();}
175  //Loop over the merged vertices to find the highest probability
176  double HighestProb = -1;
177  CandidateVertex* MostProbableVertex = 0;
178  for (std::list<CandidateVertex*>::const_iterator iCV = TrialMergedCandidates.begin();iCV != TrialMergedCandidates.end();++iCV)
179  {
180  int DegreesOfFreedom;
181  if ((*iCV)->interactionPoint())
182  {
183  //DOF = 2N when fitting N tracks with IP
184  DegreesOfFreedom = 2 * (*iCV)->trackStateList().size();
185  }
186  else
187  {
188  //DOF = 2N-2 when fitting N tracks with IP = 2(N-1)-1 as ghost is in our N
189  DegreesOfFreedom = 2 * ((*iCV)->trackStateList().size() - 1) - 2;
190  }
191  //TODO check calc above
192  std::cout << (*iCV)->chiSquaredOfFit() <<":"<<util::prob((*iCV)->chiSquaredOfFit(),DegreesOfFreedom) << " ";
193  if(util::prob((*iCV)->chiSquaredOfFit(),DegreesOfFreedom) > HighestProb)
194  {
195  HighestProb = util::prob((*iCV)->chiSquaredOfFit() , DegreesOfFreedom);
196  MostProbableVertex = (*iCV);
197  }
198  }
199  //std::cout << std::endl;
200  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "\t\tdone!" << " "<< TrialMergedCandidates.size() << "T " << Candidates.size() << "C Verts"<< "\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 << "ms" <<endl; cout.flush();}
201  //std::cout << "Highest Prob: " << MostProbableVertex << std::endl;
202  /*////////////////////////////////////////////////////////DEBUGLINE*/if (0) {cout << "Promote trial to candidate and add new trials..."; cout.flush();start=clock();}
203  if (HighestProb > _MinimumProbability)
204  {
205  //So now we need to remove the most probable vertex from the trial list and put it on the candidates list
206  TrialMergedCandidates.remove(MostProbableVertex);
207  Candidates.push_back(MostProbableVertex);
208 
209  //We then need to remove the verticies the most probable one was composed of from the candidates list as they are now represented in the candidates in merged form.
210  //We also need to remove the trial vertices that contain any of the vertices that the most probable one was composed of
211  for (std::vector<CandidateVertex*>::const_iterator iOriginalCV = VerticesContainedIn[MostProbableVertex].begin();iOriginalCV != VerticesContainedIn[MostProbableVertex].end();++iOriginalCV)
212  {
213  Candidates.remove(*iOriginalCV);
214  for (std::list<CandidateVertex*>::iterator iTrialCV = TrialMergedCandidates.begin();iTrialCV != TrialMergedCandidates.end();/*++iTrialCV*/)
215  {
216  //If we find that this Trial Vertex contained the Origial Candidate
217  if(find(VerticesContainedIn[*iTrialCV].begin(), VerticesContainedIn[*iTrialCV].end(), *iOriginalCV) != VerticesContainedIn[*iTrialCV].end())
218  {
219  //Remove it from the trial list (shuffle iterator so it doesn't become invalid)
220  std::list<CandidateVertex*>::iterator CopyiTrialCV = iTrialCV;
221  ++iTrialCV;
222  TrialMergedCandidates.erase(CopyiTrialCV);
223  }
224  else
225  ++iTrialCV;
226  }
227  }
228 
229  //We now need to add new Trials that are composed of our new candidate in combination with all the others
230  for (std::list<CandidateVertex*>::const_iterator iCV = Candidates.begin();iCV != Candidates.end();++iCV)
231  {
232  if (*iCV != MostProbableVertex)
233  {
234  std::vector<CandidateVertex*> ToMerge;
235  ToMerge.push_back(MostProbableVertex);
236  ToMerge.push_back(*iCV);
237 
238  CandidateVertex* Merged = new CandidateVertex(ToMerge);
239  MemoryManager<CandidateVertex>::Event()->registerObject(Merged);
240  TrialMergedCandidates.push_back(Merged);
241 
242  VerticesContainedIn[Merged] = ToMerge;
243  }
244  }
245  }
246  else
247  {
248  break;
249  }
250  /*////////////////////////////////////////////////////////DEBUGLINE*/if (0) {cout << "done!" << " "<< TrialMergedCandidates.size() << "T " << Candidates.size() << "C Verts"<< "\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 << "ms" <<endl; cout.flush();}
251 
252  }while (1); // Loop is broken in at line 216 only
253  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "done!" << " "<< TrialMergedCandidates.size() << "T " << Candidates.size() << "C Verts"<< "\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 << "ms" <<endl; cout.flush();}
254  /*////////////////////////////////////////////////////////DEBUGLINE*/if (debug) {cout << "Vertex find complete" << "\t" << (double(clock())-double(pstart))/CLOCKS_PER_SEC << "s" << endl; cout.flush();}
255 
256  Candidates.sort(IPDistAscending(_IP));
257  _LastGhost = GhostTrack;
258  return Candidates;
259  //Done
260 }
261 
262 Track* VertexFinderGhost::lastGhost() const
263 {
264  return _LastGhost;
265 }
266 
267 
268 }}
269 
270 
271 
272 
double prob(double ChiSquared, double DegreesOfFreedom)
Calculate probability from Chi Squared and Degrees of freedom.
const Vector3 & position() const
Return position of IP.
static MemoryManager< T > * Event()
Returns the Event duration singleton instance of the MemoryManager for type T.