1 #include "../include/vertexfinderclassic.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 "../../util/inc/memorymanager.h"
13 namespace vertex_lcfi {
namespace ZVTOP
15 VertexFinderClassic::VertexFinderClassic(
const std::vector<Track*> &Tracks, InteractionPoint* IP,
const Vector3 &JetAxis,
double Kip,
double Kalpha,
double TwoProngCut,
double TrackTrimCut,
double ResolverCutOff)
16 : _TrackList(Tracks),_IP(IP),_Kip(Kip),_Kalpha(Kalpha),_JetAxis(JetAxis),_TwoProngCut(TwoProngCut),_TrackTrimCut(TrackTrimCut),_ResolverCutOff(ResolverCutOff)
21 void VertexFinderClassic::addTrack(Track*
const Track)
23 _TrackList.push_back(Track);
26 void VertexFinderClassic::setIP(InteractionPoint*
const IP)
31 bool VertexFinderClassic::removeTrack(Track*
const TrackToRemove)
33 std::vector<Track*>::iterator iTrack;
35 for (iTrack = _TrackList.begin();((*iTrack) != TrackToRemove) && iTrack!=_TrackList.end();++iTrack)
38 if (iTrack!=_TrackList.end())
41 _TrackList.erase(iTrack);
48 bool VertexFinderClassic::clearIP()
101 return rpStart->
position().distanceTo2(IP->position()) < rpEnd->
position().distanceTo2(IP->position());
110 std::list<CandidateVertex*> VertexFinderClassic::findVertices()
112 using std::cout;
using std::endl;clock_t start,pstart;
int debug=0;
114 if (debug) {cout <<
"Constructing Vertex Function....."; cout.flush();pstart=clock();}start=clock();
117 if (debug) {cout <<
"\tdone!\t\t\t" << ((double)clock()-(double)pstart)*1000.0/CLOCKS_PER_SEC <<
"ms" << endl; cout.flush();}
121 std::list<CandidateVertex*> CVList;
122 if (debug) {cout <<
"2-Prong Creation....."; cout.flush();pstart=clock();}
123 if (debug<0) std::cout << std::endl;
125 int N = _TrackList.size();
127 std::vector<TrackState*> TrackStates;
128 for (std::vector<Track*>::iterator iTrack = _TrackList.begin();iTrack != _TrackList.end();++iTrack)
131 TrackStates.push_back(Track);
133 for (
int OuterIndex=0;OuterIndex < N-1;++OuterIndex)
135 for (
int InnerIndex=OuterIndex+1;InnerIndex < N;++InnerIndex)
137 std::vector<TrackState*> Tracks;
138 Tracks.push_back(TrackStates[OuterIndex]);
139 Tracks.push_back(TrackStates[InnerIndex]);
154 CVList.push_back(CV);
155 if (debug<0) std::cout <<
"2-prong of:" << TrackStates[OuterIndex]->parentTrack()->trackingNum() <<
"," << TrackStates[InnerIndex]->parentTrack()->trackingNum() <<
" @ " << CV->
position() << std::endl;
159 if (debug<0) std::cout <<
"-";
169 if (debug) {cout <<
"Track+IP....."; cout.flush();}
170 if (debug<0) std::cout << std::endl;
173 for (
int Index=0;Index < N;++Index)
175 std::vector<TrackState*> Tracks;
176 Tracks.push_back(TrackStates[Index]);
190 CVList.push_back(CV);
191 if (debug<0) std::cout <<
"2-prong of:" <<
"IP" <<
"," << TrackStates[Index]->parentTrack()->trackingNum() <<
" @ " << CV->
position() << std::endl;
195 if (debug<0) std::cout <<
"-";
209 if (debug) {cout <<
"\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
210 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
214 if (debug) {cout <<
"Find V(r) max....."; cout.flush();pstart=clock();}
215 for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV)
217 (*iCV)->findVertexFuncMax();
219 if (debug) {cout <<
"\t\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" <<endl; cout.flush();}
221 if (debug) {cout <<
"Track removal based on clustering....."; cout.flush();pstart=clock();}
224 std::vector<CandidateVertex*> RemoveFrom;
225 std::vector<Track*> TrackToRemove;
227 for (std::vector<Track*>::iterator iTrack = _TrackList.begin();iTrack != _TrackList.end();++iTrack)
230 std::list<CandidateVertex*> AssocCVs;
231 for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV)
233 if ((*iCV)->hasTrack(*iTrack)) AssocCVs.push_back(*iCV);
235 if (AssocCVs.empty())
244 std::list<CandidateVertex*> TempAssocCVs = AssocCVs;
245 for (std::list<CandidateVertex*>::iterator iCV = TempAssocCVs.begin();iCV != TempAssocCVs.end();++iCV)
249 RemoveFrom.push_back(*iCV);
250 TrackToRemove.push_back(*iTrack);
251 AssocCVs.remove(*iCV);
255 std::list<CandidateVertex*> RetainedCVs;
259 for (std::list<CandidateVertex*>::iterator iCV = AssocCVs.begin();iCV != AssocCVs.end();++iCV)
261 bool resolved =
true;
262 for (std::list<CandidateVertex*>::iterator iRetainedCV = RetainedCVs.begin();iRetainedCV != RetainedCVs.end();++iRetainedCV)
264 resolved = (*iCV)->isResolvedFrom((*iRetainedCV),_ResolverCutOff, CandidateVertex::FittedPosition);
268 RemoveFrom.push_back(*iCV);
269 TrackToRemove.push_back(*iTrack);
275 RetainedCVs.push_back(*iCV);
280 std::vector<Track*>::iterator iTrack = TrackToRemove.begin();
281 for (std::vector<CandidateVertex*>::iterator iCV = RemoveFrom.begin();iCV != RemoveFrom.end();++iCV)
283 (*iCV)->removeTrack(*iTrack);
287 if (debug) {cout <<
"\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
288 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
289 if (debug) {cout <<
"Remove IP from all but highest V(r)....."; cout.flush();pstart=clock();}
293 for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV)
296 if ((*iCV)->interactionPoint() && ((*iCV)->vertexFuncValue()>HighValue))
303 for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV)
304 if ((*iCV)->interactionPoint() && ((*iCV) != HighVertex))
309 std::vector<CandidateVertex*> removed;
310 for (std::list<CandidateVertex*>::iterator iVertex=CVList.begin();iVertex != CVList.end();++iVertex)
312 if ((*iVertex)->trackStateList().empty() && !(*iVertex)->interactionPoint())
314 removed.push_back(*iVertex);
317 for (std::vector<CandidateVertex*>::iterator iCVertex=removed.begin();iCVertex != removed.end();++iCVertex)
319 CVList.remove(*iCVertex);
322 if (debug) {cout <<
"done!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
323 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
324 if (debug) {cout <<
"Sort in V(r)max....."; cout.flush();pstart=clock();}
328 if (debug) {cout <<
"\t\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
330 if (debug) {cout <<
"Clustering by V(r)Max resolution....."; cout.flush();pstart=clock();}
334 std::list<std::list<CandidateVertex*> > ClusterLists;
340 CVList.erase(CVList.begin());
342 std::list<CandidateVertex*> Cluster;
343 Cluster.push_back(Seed);
349 std::list<CandidateVertex*>::iterator iCVertex;
350 for (iCVertex=Cluster.begin();iCVertex != Cluster.end();++iCVertex)
353 std::list<CandidateVertex*>::iterator iVertex;
354 for (iVertex=CVList.begin();iVertex != CVList.end();++iVertex)
357 if (!(*iCVertex)->isResolvedFrom(*iVertex,_ResolverCutOff, CandidateVertex::NearestMaximum))
360 Cluster.push_back(*iVertex);
365 for (std::list<CandidateVertex*>::iterator iCVertex2=Cluster.begin();iCVertex2 != Cluster.end();++iCVertex2)
367 CVList.remove(*iCVertex2);
374 ClusterLists.push_back(Cluster);
376 while (!CVList.empty());
377 if (debug) {cout <<
"\tdone!" <<
"\t\t\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
378 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
379 if (debug) {cout <<
"Merging....."; cout.flush();pstart=clock();}
381 std::list<std::list<CandidateVertex*> >::iterator iList;
382 for (iList=ClusterLists.begin();iList != ClusterLists.end();++iList)
384 std::list<CandidateVertex*>::iterator iVertex;
385 for (iVertex=(++((*iList).begin()));iVertex != (*iList).end(); ++iVertex)
387 (*((*iList).begin()))->mergeCandidateVertex(*iVertex);
390 CVList.push_back(*((*iList).begin()));
394 _ifNoIPAddIP(&CVList);
396 if (debug) {cout <<
"\t\t\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
397 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
398 if (debug) {cout <<
"Sort in V(r)Max....."; cout.flush();pstart=clock();}
402 if (debug) {cout <<
"\t\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
404 if (debug) {cout <<
"Chi Squared Track Cut....."; cout.flush();pstart=clock();}
406 for (std::list<CandidateVertex*>::iterator iVertex=CVList.begin();iVertex != CVList.end();++iVertex)
408 (*iVertex)->trimByChi2(_TrackTrimCut);
410 if (debug) {cout <<
"\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
412 if (debug) {cout <<
"Discard <2 Track Vertices....."; cout.flush();}pstart=clock();
414 _removeOneTrackNoIPVertices(&CVList);
415 if (debug) {cout <<
"\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
416 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
417 if (debug) {cout <<
"Claim tracks by V(r)....."; cout.flush();pstart=clock();}
419 std::list<CandidateVertex*> LosingList;
420 std::list<CandidateVertex*> LeftToDo;
421 for (std::list<CandidateVertex*>::iterator iCVertex=CVList.begin();iCVertex != CVList.end();++iCVertex)
423 LosingList.push_back(*iCVertex);
424 LeftToDo.push_back(*iCVertex);
426 while (!LeftToDo.empty())
429 LosingList.erase(LosingList.begin());
431 (*LeftToDo.begin())->claimTracksFrom(LosingList);
432 LeftToDo.erase(LeftToDo.begin());
434 std::vector<CandidateVertex*> removed = _removeOneTrackNoIPVertices(&CVList);
436 for (std::vector<CandidateVertex*>::iterator iCVertex=removed.begin();iCVertex != removed.end();++iCVertex)
438 LosingList.remove(*iCVertex);
439 LeftToDo.remove(*iCVertex);
443 if (debug) {cout <<
"\t\tdone!" <<
" "<< CVList.size() <<
" Vertices" <<
"\t" << ((double(clock())-double(pstart))/CLOCKS_PER_SEC)*1000 <<
"ms" << endl; cout.flush();}
444 if (debug>1) {
for (std::list<CandidateVertex*>::iterator iCV = CVList.begin();iCV != CVList.end();++iCV) cout << **iCV <<endl;}
445 if (debug) {cout <<
"Vertex find complete" <<
"\t" << (double(clock())-double(start))/CLOCKS_PER_SEC <<
"s" <<
" " << ((
double(clock())-
double(start))*1000)/(CLOCKS_PER_SEC*_TrackList.size()*_TrackList.size()) <<
"ms per track^2" << endl; cout.flush();}
454 std::vector<CandidateVertex*> VertexFinderClassic::_removeOneTrackNoIPVertices(std::list<CandidateVertex*>* CVList)
457 std::vector<CandidateVertex*> removed;
458 for (std::list<CandidateVertex*>::iterator iVertex=CVList->begin();iVertex != CVList->end();++iVertex)
460 if ((*iVertex)->trackStateList().size() < 2)
463 if (!(*iVertex)->interactionPoint())
464 removed.push_back(*iVertex);
467 for (std::vector<CandidateVertex*>::iterator iCVertex=removed.begin();iCVertex != removed.end();++iCVertex)
469 CVList->remove(*iCVertex);
476 void VertexFinderClassic::_ifNoIPAddIP(std::list<CandidateVertex*>* CVList)
479 for (std::list<CandidateVertex*>::iterator iVertex=CVList->begin();iVertex != CVList->end();++iVertex)
481 if ((*iVertex)->interactionPoint())
return;
484 std::vector<TrackState*> Tracks;
485 CandidateVertex* CV =
new CandidateVertex(Tracks,_IP,_VF);
486 MemoryManager<CandidateVertex>::Event()->registerObject(CV);
487 CVList->push_back(CV);
double vertexFuncMaxValue() const
Return the value of the vertex function at this vertexes local maximum.
void registerObject(T *pointer)
Register an object for memory management.
VertexFunction as in ZVTOP paper.
Interaction Point representation.
double maxChiSquaredOfTrackIP() const
Return the chi squared contribution of the Track or IP with the highest chi square contribution...
A collection of TrackState objects with a fit and vertex function maximum.
double vertexFuncValue() const
Return the value of the vertex function at the vertices position.
const Vector3 & position() const
Return the fitted position of this Vertex.
Spatial point on a Track.
bool removeIP()
Remove this vertices InteractionPoint.