1 #include "../include/vertexfinderghost.h"
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"
17 namespace vertex_lcfi {
namespace ZVTOP
21 struct IPDistAscending
25 IPDistAscending(InteractionPoint* _IP) : IP(_IP) {}
27 bool operator()(CandidateVertex*& rpStart, CandidateVertex*& rpEnd)
31 return rpStart->position().distanceTo2(IP->
position()) < rpEnd->position().distanceTo2(IP->
position());
35 return rpStart->position().distanceTo2(
Vector3(0,0,0)) < rpEnd->position().distanceTo2(
Vector3(0,0,0));
40 VertexFinderGhost::VertexFinderGhost(
const std::vector<Track*> &Tracks, InteractionPoint* IP)
41 :_TrackList(Tracks),_IP(IP)
46 void VertexFinderGhost::addTrack(Track*
const Track)
48 _TrackList.push_back(Track);
51 void VertexFinderGhost::setIP(InteractionPoint*
const IP)
56 bool VertexFinderGhost::removeTrack(Track*
const TrackToRemove)
58 std::vector<Track*>::iterator iTrack;
60 for (iTrack = _TrackList.begin();((*iTrack) != TrackToRemove) && iTrack!=_TrackList.end();++iTrack)
63 if (iTrack!=_TrackList.end())
66 _TrackList.erase(iTrack);
73 bool VertexFinderGhost::clearIP()
80 std::list<CandidateVertex*> VertexFinderGhost::findVertices()
87 using std::cout;
using std::endl;clock_t start,pstart;
int debug=0;
89 if (_TrackList.empty())
91 if (debug) std::cout <<
"Ghost: No Tracks" << std::endl;
93 return std::list<CandidateVertex*>();
96 std::list<CandidateVertex*> ret;
97 CandidateVertex* CV =
new CandidateVertex(std::vector<TrackState*>(),_IP,0);
103 if (debug) {cout <<
"Initial Ghost direction finding and resize...."; cout.flush();pstart=clock();start=clock();}
105 Track* GhostTrack = GhostFinderStage1().findGhost(_InitialGhostWidth,_MaxChi2Allowed,_SeedDirection,_TrackList,_IP);
106 if (debug) {cout <<
"\tdone!\t\t\t" << ((double)clock()-(double)pstart)*1000.0/CLOCKS_PER_SEC <<
"ms" << endl; cout.flush();}
109 std::vector<TrackState*> TrackStates;
110 for (std::vector<Track*>::iterator iTrack = _TrackList.begin();iTrack != _TrackList.end();++iTrack)
112 TrackState* Track = (*iTrack)->makeState();
113 TrackStates.push_back(Track);
116 TrackState* GhostTrackState = GhostTrack->makeState();
118 if (debug) {cout <<
"Generating 1-prong Candidate Vertices + IP...."; cout.flush();start=clock();}
120 std::list<CandidateVertex*> Candidates;
122 for( std::vector<TrackState*>::const_iterator iTrack=TrackStates.begin(); iTrack != TrackStates.end(); ++iTrack)
125 std::vector<TrackState*> Tracks;
126 Tracks.push_back(*iTrack);
127 Tracks.push_back(GhostTrackState);
128 CandidateVertex* CV =
new CandidateVertex(Tracks,(InteractionPoint*)0,0);
130 Candidates.push_back(CV);
134 std::vector<TrackState*> Tracks;
135 CandidateVertex* CV =
new CandidateVertex(Tracks,_IP,0);
137 Candidates.push_back(CV);
139 if (debug) {cout <<
"\tdone!" <<
" "<< Candidates.size() <<
" Candidates" <<
"\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 <<
"ms" <<endl; cout.flush();}
140 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = Candidates.begin();iCV != Candidates.end();++iCV) cout << *iCV <<endl;}
142 if (debug) {cout <<
"Generating Trial Merged Vertices...."; cout.flush();start=clock();}
144 std::list<CandidateVertex*> TrialMergedCandidates;
146 std::map<CandidateVertex*,std::vector<CandidateVertex*> > VerticesContainedIn;
148 for (std::list<CandidateVertex*>::const_iterator iOuterCV = Candidates.begin();iOuterCV != --(Candidates.end());++iOuterCV)
150 std::list<CandidateVertex*>::const_iterator CopyiOuterCV = iOuterCV;
151 for (std::list<CandidateVertex*>::const_iterator iInnerCV = (++CopyiOuterCV);iInnerCV != Candidates.end();++iInnerCV)
153 std::vector<CandidateVertex*> ToMerge;
154 ToMerge.push_back(*iOuterCV);
155 ToMerge.push_back(*iInnerCV);
157 CandidateVertex* Merged =
new CandidateVertex(ToMerge);
160 if (Merged->hasTrack(GhostTrack) && Merged->interactionPoint())
162 Merged->removeTrack(GhostTrack);
165 TrialMergedCandidates.push_back(Merged);
167 VerticesContainedIn[Merged] = ToMerge;
170 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 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = TrialMergedCandidates.begin();iCV != TrialMergedCandidates.end();++iCV) cout << *iCV <<endl;}
174 if (debug) {cout <<
"Find most probable trial vertex..."; cout.flush();start=clock();}
176 double HighestProb = -1;
177 CandidateVertex* MostProbableVertex = 0;
178 for (std::list<CandidateVertex*>::const_iterator iCV = TrialMergedCandidates.begin();iCV != TrialMergedCandidates.end();++iCV)
180 int DegreesOfFreedom;
181 if ((*iCV)->interactionPoint())
184 DegreesOfFreedom = 2 * (*iCV)->trackStateList().size();
189 DegreesOfFreedom = 2 * ((*iCV)->trackStateList().size() - 1) - 2;
192 std::cout << (*iCV)->chiSquaredOfFit() <<
":"<<
util::prob((*iCV)->chiSquaredOfFit(),DegreesOfFreedom) <<
" ";
193 if(
util::prob((*iCV)->chiSquaredOfFit(),DegreesOfFreedom) > HighestProb)
195 HighestProb =
util::prob((*iCV)->chiSquaredOfFit() , DegreesOfFreedom);
196 MostProbableVertex = (*iCV);
200 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();}
202 if (0) {cout <<
"Promote trial to candidate and add new trials..."; cout.flush();start=clock();}
203 if (HighestProb > _MinimumProbability)
206 TrialMergedCandidates.remove(MostProbableVertex);
207 Candidates.push_back(MostProbableVertex);
211 for (std::vector<CandidateVertex*>::const_iterator iOriginalCV = VerticesContainedIn[MostProbableVertex].begin();iOriginalCV != VerticesContainedIn[MostProbableVertex].end();++iOriginalCV)
213 Candidates.remove(*iOriginalCV);
214 for (std::list<CandidateVertex*>::iterator iTrialCV = TrialMergedCandidates.begin();iTrialCV != TrialMergedCandidates.end();)
217 if(find(VerticesContainedIn[*iTrialCV].begin(), VerticesContainedIn[*iTrialCV].end(), *iOriginalCV) != VerticesContainedIn[*iTrialCV].end())
220 std::list<CandidateVertex*>::iterator CopyiTrialCV = iTrialCV;
222 TrialMergedCandidates.erase(CopyiTrialCV);
230 for (std::list<CandidateVertex*>::const_iterator iCV = Candidates.begin();iCV != Candidates.end();++iCV)
232 if (*iCV != MostProbableVertex)
234 std::vector<CandidateVertex*> ToMerge;
235 ToMerge.push_back(MostProbableVertex);
236 ToMerge.push_back(*iCV);
238 CandidateVertex* Merged =
new CandidateVertex(ToMerge);
240 TrialMergedCandidates.push_back(Merged);
242 VerticesContainedIn[Merged] = ToMerge;
250 if (0) {cout <<
"done!" <<
" "<< TrialMergedCandidates.size() <<
"T " << Candidates.size() <<
"C Verts"<<
"\t" << ((double(clock())-double(start))/CLOCKS_PER_SEC)*1000 <<
"ms" <<endl; cout.flush();}
253 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 if (debug) {cout <<
"Vertex find complete" <<
"\t" << (double(clock())-double(pstart))/CLOCKS_PER_SEC <<
"s" << endl; cout.flush();}
256 Candidates.sort(IPDistAscending(_IP));
257 _LastGhost = GhostTrack;
262 Track* VertexFinderGhost::lastGhost()
const
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.