CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFDisplacedVertexCandidateFinder.cc
Go to the documentation of this file.
2 
6 
8 
9 using namespace std;
10 using namespace reco;
11 
12 //for debug only
13 //#define PFLOW_DEBUG
14 
16  vertexCandidates_( new PFDisplacedVertexCandidateCollection ),
17  dcaCut_(1000),
18  primaryVertexCut2_(0.0),
19  dcaPInnerHitCut2_(1000.0),
20  vertexCandidatesSize_(50),
21  debug_(false) {
22 
24  theMinimum_ = theMinimum;
25 }
26 
27 
29 
30 #ifdef PFLOW_DEBUG
31  if(debug_)
32  cout<<"~PFDisplacedVertexCandidateFinder - number of remaining elements: "
33  <<eventTracks_.size()<<endl;
34 #endif
35 
36 }
37 
38 
40  edm::Handle< reco::VertexCollection > mainVertexHandle,
41  edm::Handle< reco::BeamSpot > beamSpotHandle){
42 
43  const math::XYZPoint beamSpot = beamSpotHandle.isValid() ?
44  math::XYZPoint(beamSpotHandle->x0(), beamSpotHandle->y0(), beamSpotHandle->z0()) :
45  math::XYZPoint(0, 0, 0);
46 
47  // The primary vertex is taken from the refitted list,
48  // if does not exist from the average offline beam spot position
49  // if does not exist (0,0,0) is used
50  pvtx_ = mainVertexHandle.isValid() ?
51  math::XYZPoint(mainVertexHandle->begin()->x(),
52  mainVertexHandle->begin()->y(),
53  mainVertexHandle->begin()->z()) :
54  beamSpot;
55 }
56 
57 // Set the imput collection of tracks and calculate their
58 // trajectory parameters the Global Trajectory Parameters
59 void
61  const MagneticField* magField) {
62 
63  magField_ = magField;
64  trackMask_.clear();
65  trackMask_.resize(trackh->size());
67  eventTrackTrajectories_.resize(trackh->size());
68 
69  for(unsigned i=0;i<trackMask_.size(); i++) trackMask_[i] = true;
70 
71  eventTracks_.clear();
72  if(trackh.isValid()) {
73  for(unsigned i=0;i<trackh->size(); i++) {
74 
75  TrackRef tref( trackh, i);
76  TrackBaseRef tbref(tref);
77 
78  if( !isSelected( tbref ) ) {
79  trackMask_[i] = false;
80  continue;
81  }
82 
83  const Track* trk = tref.get();
84  eventTracks_.push_back( tbref );
86  }
87  }
88 }
89 
90 
91 // -------- Main function which find vertices -------- //
92 
93 void
95 
96  if (debug_) cout << "========= Start Finding Displaced Vertex Candidates =========" << endl;
97 
98  // The vertexCandidates have not been passed to the event, and need to be cleared
99  if(vertexCandidates_.get() )vertexCandidates_->clear();
100  else
102 
104  for(IE ie = eventTracks_.begin();
105  ie != eventTracks_.end();) {
106 
107  // Run the recursive procedure to find all tracks link together
108  // In one blob called Candidate
109 
110  PFDisplacedVertexCandidate tempVertexCandidate;
111 
112  ie = associate( eventTracks_.end() , ie, tempVertexCandidate);
113 
114  // Build remaining links in current block
115  if(tempVertexCandidate.isValid()) {
116  packLinks( tempVertexCandidate );
117  vertexCandidates_->push_back( tempVertexCandidate );
118  }
119  }
120 
121  if(debug_) cout << "========= End Finding Displaced Vertex Candidates =========" << endl;
122 
123 }
124 
125 
126 
127 
130  PFDisplacedVertexCandidate& tempVertexCandidate) {
131 
132 
133 #ifdef PFLOW_DEBUG
134  if(debug_ ) cout<<"== Start the association procedure =="<<endl;
135 #endif
136 
137  if( last!= eventTracks_.end() ) {
138  double dist = -1;
139  GlobalPoint P(0,0,0);
141  link( (*last), (*next), dist, P, linktest);
142 
143 
144  if(dist<-0.5) {
145 #ifdef PFLOW_DEBUG
146  if(debug_ ) cout<<"link failed"<<endl;
147 #endif
148  return ++next; // association failed
149  }
150  else {
151  // add next element to the current pflowblock
152  tempVertexCandidate.addElement( (*next));
153  trackMask_[(*next).key()] = false;
154 #ifdef PFLOW_DEBUG
155  if(debug_ )
156  cout<<"link parameters "
157  << " *next = " << (*next).key()
158  << " *last = " << (*last).key()
159  << " dist = " << dist
160  << " P.x = " << P.x()
161  << " P.y = " << P.y()
162  << " P.z = " << P.z()
163  << endl;
164 #endif
165  }
166  }
167  else {
168  // add next element to this eflowblock
169 #ifdef PFLOW_DEBUG
170  if(debug_ ) cout<<"adding to block element "
171  << (*next).key()
172  <<endl;
173 #endif
174  tempVertexCandidate.addElement( (*next));
175  trackMask_[(*next).key()] = false;
176 
177 
178  }
179 
180  // recursive call: associate next and other unused elements
181 #ifdef PFLOW_DEBUG
182  if(debug_ ) {
183  for(unsigned i=0; i<trackMask_.size(); i++)
184  cout << " Mask[" << i << "] = " << trackMask_[i];
185  cout << "" << endl;
186  }
187 #endif
188 
189  for(IE ie = eventTracks_.begin();
190  ie != eventTracks_.end();) {
191 
192  if( ie == last || ie == next ) { ++ie; continue;}
193 
194  // *ie already included to a block
195  if( !trackMask_[(*ie).key()] ) { ++ie; continue;}
196 
197 #ifdef PFLOW_DEBUG
198  if(debug_ ) cout<<"calling associate "
199  << (*next).key()
200  <<" & "
201  << (*ie).key()
202  <<endl;
203 #endif
204  ie = associate(next, ie, tempVertexCandidate);
205 
206  }
207 
208 #ifdef PFLOW_DEBUG
209  if(debug_ ) {
210  cout<<"**** removing element "<<endl;
211  }
212 #endif
213 
214  IE iteratorToNextFreeElement = eventTracks_.erase( next );
215 
216 #ifdef PFLOW_DEBUG
217  if(debug_ ) cout<< "== End the association procedure ==" <<endl;
218 #endif
219 
220  return iteratorToNextFreeElement;
221 }
222 
223 
224 
225 void
227  const TrackBaseRef& el2,
228  double& dist, GlobalPoint& P,
230  ){
231 
232  if ( fabs(el1->eta()-el2->eta()) > 1) {dist = -1; return;}
233  if ( el1->pt()>2 && el2->pt()>2 && fabs(el1->phi()-el2->phi()) > 1) {dist = -1; return;}
234 
237 
238 
239  // Closest approach
240  theMinimum_.calculate(gt1,gt2);
241 
242  // Fill the parameters
243  dist = theMinimum_.distance();
245 
246  vertexLinkTest = PFDisplacedVertexCandidate::LINKTEST_DCA; //rechit by default
247 
248 
249  // Check if the link test fails
250  if (dist > dcaCut_) {dist = -1; return;}
251 
252  // Check if the closses approach point is too close to the primary vertex/beam pipe
253  double rho2 = P.x()*P.x()+P.y()*P.y();
254 
255  if ( rho2 < primaryVertexCut2_) {dist = -1; return;}
256 
257  // Check if the closses approach point is out of range of the tracker
258 
259  double tob_rho_limit2 = 10000;
260  double tec_z_limit = 270;
261 
262  if ( rho2 > tob_rho_limit2) {dist = -1; return;}
263  if ( fabs(P.z()) > tec_z_limit) {dist = -1; return;}
264 
265  return;
266 
267 
268  /*
269  // Check if the inner hit of one of the tracks is not too far away from the vertex
270  double dx = el1->innerPosition().x()-P.x();
271  double dy = el1->innerPosition().y()-P.y();
272  double dz = el1->innerPosition().z()-P.z();
273  double dist2 = dx*dx+dy*dy+dz*dz;
274 
275  if (dist2 > dcaPInnerHitCut2_) {
276  dist = -1;
277  if (debug_) cout << "track " << el1.key() << " dist to vertex " << sqrt(dist2) << endl;
278  return;
279  }
280 
281  dx = el2->innerPosition().x()-P.x();
282  dy = el2->innerPosition().y()-P.y();
283  dz = el2->innerPosition().z()-P.z();
284  dist2 = dx*dx+dy*dy+dz*dz;
285 
286  if (dist2 > dcaPInnerHitCut2_) {
287  dist = -1;
288  if (debug_) cout << "track " << el2.key() << " dist to vertex " << sqrt(dist2) << endl;
289  return;
290  }
291  */
292 
293 }
294 
295 
296 // Build up a matrix of all the links between different tracks in
297 // In the Candidate
298 void
300 
301 
302  const vector < TrackBaseRef >& els = vertexCandidate.elements();
303 
304  //First Loop: update all link data
305  for( unsigned i1=0; i1<els.size(); i1++ ) {
306  for( unsigned i2=i1+1; i2<els.size(); i2++ ) {
307 
308  // no reflexive link
309  if( i1==i2 ) continue;
310 
311  double dist = -1;
312  GlobalPoint P(0,0,0);
314 
315  link( els[i1], els[i2], dist, P, linktest);
316 
317 
318 #ifdef PFLOW_DEBUG
319  if( debug_ )
320  cout << "Setting link between elements " << i1 << " key " << els[i1].key()
321  << " and " << i2 << " key " << els[i2].key()
322  << " of dist =" << dist << " computed from link test "
323  << linktest << endl;
324 #endif
325 
326  if(dist >-0.5) vertexCandidate.setLink( i1, i2, dist, P, linktest );
327  }
328  }
329 
330 }
331 
332 
333 
334 // --------------- TOOLS -------------- //
335 
336 
337 
338 // This tool is a copy from VZeroFinder
341 (const Track* track) const
342 {
343 
344  const GlobalPoint position(track->vx(),
345  track->vy(),
346  track->vz());
347 
348  const GlobalVector momentum(track->momentum().x(),
349  track->momentum().y(),
350  track->momentum().z());
351 
353  track->charge(),magField_);
354 
355  return gtp;
356 }
357 
358 
359 // This tool allow to pre-select the track for the displaced vertex finder
360 bool
362 
363  double nChi2 = trackref->normalizedChi2();
364  double pt = trackref->pt();
365  double dpt = trackref->ptError();
366  double dxy = trackref->dxy(pvtx_);
367 
368  double pt_error = dpt/pt*100;
369 
370  if (debug_) cout << " PFDisplacedVertexFinder: PFrecTrack->Track Pt= "
371  << pt << " dPt/Pt = " << pt_error << "% nChi2 = " << nChi2 << endl;
372  if (nChi2 > nChi2_max_ || pt < pt_min_){
373 
374  if (debug_) cout << " PFBlockAlgo: skip badly measured or low pt track"
375  << " nChi2_cut = " << 5
376  << " pt_cut = " << 0.2 << endl;
377  return false;
378  }
379  // cout << "dxy = " << dxy << endl;
380  if (fabs(dxy) < dxy_ && pt < pt_min_prim_) return false;
381  // if (fabs(dxy) < 0.2 && pt < 0.8) return false;
382 
383 
384 
385  return true;
386 }
387 
388 
389 
390 
391 
392 
393 
394 
395 
396 
397 
398 ostream& operator<<(std::ostream& out, const PFDisplacedVertexCandidateFinder& a) {
399  if(! out) return out;
400 
401  out<<"====== Particle Flow Block Algorithm ======= ";
402  out<<endl;
403  out<<"number of unassociated elements : "<<a.eventTracks_.size()<<endl;
404  out<<endl;
405 
406  out << " Tracks selection based on " << std::endl;
407  out << " pvtx_ = " << a.pvtx_ << std::endl;
408  out << " fabs(dxy) < " << a.dxy_ << " and pt < "<< a.pt_min_prim_ << std::endl;
409  out << " nChi2 < " << a.nChi2_max_ << " and pt < "<< a.pt_min_ << std::endl;
410 
411  out<<endl;
412 
413 
415  ie != a.eventTracks_.end(); ie++) {
416 
417  double pt = (*ie).get()->pt();
418 
419  math::XYZPoint Pi = (*ie).get()->innerPosition();
420  math::XYZPoint Po = (*ie).get()->outerPosition();
421 
422  double innermost_radius = sqrt(Pi.x()*Pi.x() + Pi.y()*Pi.y() + Pi.z()*Pi.z());
423  double outermost_radius = sqrt(Po.x()*Po.x() + Po.y()*Po.y() + Po.z()*Po.z());
424  double innermost_rho = sqrt(Pi.x()*Pi.x() + Pi.y()*Pi.y());
425  double outermost_rho = sqrt(Po.x()*Po.x() + Po.y()*Po.y());
426 
427  out<<"ie = " << (*ie).key()
428  <<" pt = " << pt
429  <<" innermost hit radius = " << innermost_radius << " rho = " << innermost_rho
430  <<" outermost hit radius = " << outermost_radius << " rho = " << outermost_rho
431  <<endl;
432  }
433 
434 
435  const std::auto_ptr< reco::PFDisplacedVertexCandidateCollection >& vertexCandidates
436  = a.vertexCandidates();
437 
438  if(!vertexCandidates.get() ) {
439  out<<"vertexCandidates already transfered"<<endl;
440  }
441  else {
442  out<<"number of vertexCandidates : "<<vertexCandidates->size()<<endl;
443  out<<endl;
444 
445 
446  for(PFDisplacedVertexCandidateFinder::IBC ib=vertexCandidates->begin();
447  ib != vertexCandidates->end(); ib++)
448  ib->Dump();
449 
450 
451 
452  }
453 
454  return out;
455 }
virtual float distance() const
void setInput(const edm::Handle< reco::TrackCollection > &trackh, const MagneticField *magField)
const double Pi
void addElement(const TrackBaseRef)
add a track Reference to the current Candidate
void link(const reco::TrackBaseRef &el1, const reco::TrackBaseRef &el2, double &dist, GlobalPoint &P, reco::PFDisplacedVertexCandidate::VertexLinkTest &linktest)
Check whether 2 elements are linked and fill the link parameters.
int i
Definition: DBlmapReader.cc:9
virtual bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb)
std::vector< PFDisplacedVertexCandidate > PFDisplacedVertexCandidateCollection
collection of PFDisplacedVertexCandidate objects
int ib
Definition: cuy.py:660
A block of tracks linked together.
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:609
GlobalTrajectoryParameters getGlobalTrajectoryParameters(const reco::Track *) const
--—— TOOLS --—— //
virtual GlobalPoint crossingPoint() const
IE associate(IE next, IE last, reco::PFDisplacedVertexCandidate &tempVertexCandidate)
--—— Different steps of the finder algorithm --—— ///
void setPrimaryVertex(edm::Handle< reco::VertexCollection > mainVertexHandle, edm::Handle< reco::BeamSpot > beamSpotHandle)
T y() const
Definition: PV3DBase.h:63
std::vector< GlobalTrajectoryParameters > eventTrackTrajectories_
#define P
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:693
tuple els
Definition: asciidump.py:420
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:723
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
std::auto_ptr< reco::PFDisplacedVertexCandidateCollection > vertexCandidates_
-----— Members -----— ///
std::list< reco::TrackBaseRef >::const_iterator IEC
void setLink(unsigned i1, unsigned i2, const float dist, const GlobalPoint &dcaPoint, const VertexLinkTest test=LINKTEST_DCA)
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:699
size_t key() const
Definition: RefToBase.h:229
T sqrt(T t)
Definition: SSEVec.h:48
double dcaCut_
–— Algo parameters for the vertex finder -— ///
double pt() const
track transverse momentum
Definition: TrackBase.h:669
T z() const
Definition: PV3DBase.h:64
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
Definition: TrackBase.h:811
bool isValid() const
A Vertex Candidate is valid if it has at least two tracks.
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
bool isValid() const
Definition: HandleBase.h:76
tuple out
Definition: dbtoconf.py:99
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:717
reco::PFDisplacedVertexCandidateCollection::const_iterator IBC
const std::auto_ptr< reco::PFDisplacedVertexCandidateCollection > & vertexCandidates() const
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
bool goodPtResolution(const reco::TrackBaseRef &trackref) const
Quality Criterion on the Pt resolution to select a Track.
bool isSelected(const reco::TrackBaseRef &trackref)
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:711
double a
Definition: hdecay.h:121
static int position[264][3]
Definition: ReadPGInfo.cc:509
void packLinks(reco::PFDisplacedVertexCandidate &vertexCandidate)
void findDisplacedVertexCandidates()
-----— Main function which find vertices -----— ///
double primaryVertexCut2_
Do not reconstruct vertices wich are too close to the beam pipe.
std::list< reco::TrackBaseRef > eventTracks_
The track refs.
tuple cout
Definition: gather_cfg.py:121
int charge() const
track electric charge
Definition: TrackBase.h:615
volatile std::atomic< bool > shutdown_flag false
const std::vector< TrackBaseRef > & elements() const
std::list< reco::TrackBaseRef >::iterator IE
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:639
T x() const
Definition: PV3DBase.h:62
bool debug_
if true, debug printouts activated
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:705