CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackerOnlyConversionProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackerOnlyConversionProducer
4 // Class: TrackerOnlyConversionProducer
5 //
13 //
14 // Original Author: Hongliang Liu
15 // Created: Thu Mar 13 17:40:48 CDT 2008
16 // $Id: TrackerOnlyConversionProducer.cc,v 1.40 2010/11/22 02:02:08 bendavid Exp $
17 //
18 //
19 
20 
21 // system include files
22 #include <memory>
23 #include <map>
24 
25 
31 
36 
39 
40 
45 
49 
52 
59 
60 
61 //Kinematic constraint vertex fitter
71 
72 
73 
75  theVertexFinder_(0)
76 
77 {
78  algoName_ = iConfig.getParameter<std::string>( "AlgorithmName" );
79 
80  src_ = iConfig.getParameter<edm::InputTag>("src");
81  maxNumOfTrackInPU_ = iConfig.getParameter<int>("maxNumOfTrackInPU");
82  allowTrackBC_ = iConfig.getParameter<bool>("AllowTrackBC");
83  allowD0_ = iConfig.getParameter<bool>("AllowD0");
84  allowDeltaPhi_ = iConfig.getParameter<bool>("AllowDeltaPhi");
85  allowDeltaCot_ = iConfig.getParameter<bool>("AllowDeltaCot");
86  allowMinApproach_ = iConfig.getParameter<bool>("AllowMinApproach");
87  allowOppCharge_ = iConfig.getParameter<bool>("AllowOppCharge");
88 
89  allowVertex_ = iConfig.getParameter<bool>("AllowVertex");
90 
91  halfWayEta_ = iConfig.getParameter<double>("HalfwayEta");//open angle to search track matches with BC
92 
93  if (allowD0_)
94  d0Cut_ = iConfig.getParameter<double>("d0");
95 
96  usePvtx_ = iConfig.getParameter<bool>("UsePvtx");//if use primary vertices
97 
98  if (usePvtx_){
99  vertexProducer_ = iConfig.getParameter<std::string>("primaryVertexProducer");
100  }
101 
102  if (allowTrackBC_) {
103  //Track-cluster matching eta and phi cuts
104  dEtaTkBC_ = iConfig.getParameter<double>("dEtaTrackBC");//TODO research on cut endcap/barrel
105  dPhiTkBC_ = iConfig.getParameter<double>("dPhiTrackBC");
106 
107  bcBarrelCollection_ = iConfig.getParameter<edm::InputTag>("bcBarrelCollection");
108  bcEndcapCollection_ = iConfig.getParameter<edm::InputTag>("bcEndcapCollection");
109 
110  energyBC_ = iConfig.getParameter<double>("EnergyBC");//BC energy cut
111  energyTotalBC_ = iConfig.getParameter<double>("EnergyTotalBC");//BC energy cut
112 
113  }
114 
115  //Track cuts on left right track: at least one leg reaches ECAL
116  //Left track: must exist, must reach Ecal and match BC, so loose cut on Chi2 and tight on hits
117  //Right track: not necessary to exist (if allowSingleLeg_), not necessary to reach ECAL or match BC, so tight cut on Chi2 and loose on hits
118  maxChi2Left_ = iConfig.getParameter<double>("MaxChi2Left");
119  maxChi2Right_ = iConfig.getParameter<double>("MaxChi2Right");
120  minHitsLeft_ = iConfig.getParameter<int>("MinHitsLeft");
121  minHitsRight_ = iConfig.getParameter<int>("MinHitsRight");
122 
123  //Track Open angle cut on delta cot(theta) and delta phi
124  if (allowDeltaCot_)
125  deltaCotTheta_ = iConfig.getParameter<double>("DeltaCotTheta");
126  if (allowDeltaPhi_)
127  deltaPhi_ = iConfig.getParameter<double>("DeltaPhi");
128  if (allowMinApproach_){
129  minApproachLow_ = iConfig.getParameter<double>("MinApproachLow");
130  minApproachHigh_ = iConfig.getParameter<double>("MinApproachHigh");
131  }
132 
133  // if allow single track collection, by default False
134  allowSingleLeg_ = iConfig.getParameter<bool>("AllowSingleLeg");
135  rightBC_ = iConfig.getParameter<bool>("AllowRightBC");
136 
137  //track inner position dz cut, need RECO
138  dzCut_ = iConfig.getParameter<double>("dz");
139  //track analytical cross cut
140  r_cut = iConfig.getParameter<double>("rCut");
141  vtxChi2_ = iConfig.getParameter<double>("vtxChi2");
142 
143 
144  theVertexFinder_ = new ConversionVertexFinder ( iConfig );
145 
146  thettbuilder_ = 0;
147 
148  //output
149  ConvertedPhotonCollection_ = iConfig.getParameter<std::string>("convertedPhotonCollection");
150 
151  produces< reco::ConversionCollection >(ConvertedPhotonCollection_);
152 
153 }
154 
155 
157 {
158 
159  // do anything here that needs to be done at desctruction time
160  // (e.g. close files, deallocate resources etc.)
161  delete theVertexFinder_;
162 }
163 
164 
165 // ------------ method called to produce the data ------------
166  void
168 {
169  using namespace edm;
170 
171  reco::ConversionCollection outputConvPhotonCollection;
172  std::auto_ptr<reco::ConversionCollection> outputConvPhotonCollection_p(new reco::ConversionCollection);
173 
174  //std::cout << " TrackerOnlyConversionProducer::produce " << std::endl;
175  //Read multiple track input collections
176 
177  edm::Handle<reco::ConversionTrackCollection> trackCollectionHandle;
178  iEvent.getByLabel(src_,trackCollectionHandle);
179 
180  edm::Handle<edm::View<reco::CaloCluster> > bcBarrelHandle;//TODO error handling if no collection found
181  edm::Handle<edm::View<reco::CaloCluster> > bcEndcapHandle;//TODO check cluster type if BasicCluster or PFCluster
182  if (allowTrackBC_){
183  iEvent.getByLabel( bcBarrelCollection_, bcBarrelHandle);
184  iEvent.getByLabel( bcEndcapCollection_, bcEndcapHandle);
185  }
186 
189  if (usePvtx_){
190  iEvent.getByLabel(vertexProducer_, vertexHandle);
191  if (!vertexHandle.isValid()) {
192  edm::LogError("TrackerOnlyConversionProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
193  usePvtx_ = false;
194  }
195  if (usePvtx_)
196  vertexCollection = *(vertexHandle.product());
197  }
198 
199  edm::ESHandle<TransientTrackBuilder> hTransientTrackBuilder;
200  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",hTransientTrackBuilder);
201  thettbuilder_ = hTransientTrackBuilder.product();
202 
203  reco::Vertex the_pvtx;
204  //because the priamry vertex is sorted by quality, the first one is the best
205  if (!vertexCollection.empty())
206  the_pvtx = *(vertexCollection.begin());
207 
208  if (trackCollectionHandle->size()> maxNumOfTrackInPU_ ){
209  iEvent.put( outputConvPhotonCollection_p, ConvertedPhotonCollection_);
210  return;
211  }
212 
213  //2. select track by propagating
214  //2.0 build Basic cluster geometry map to search in eta bounds for clusters
215  std::multimap<double, reco::CaloClusterPtr> basicClusterPtrs;
216  edm::Handle<edm::View<reco::CaloCluster> > bcHandle = bcBarrelHandle;
217 
218  if (allowTrackBC_){
219  for (unsigned jj = 0; jj < 2; ++jj ){
220  for (unsigned ii = 0; ii < bcHandle->size(); ++ii ) {
221  if (bcHandle->ptrAt(ii)->energy()>energyBC_)
222  basicClusterPtrs.insert(std::make_pair(bcHandle->ptrAt(ii)->position().eta(), bcHandle->ptrAt(ii)));
223  }
224  bcHandle = bcEndcapHandle;
225  }
226  }
227 
228  buildCollection( iEvent, iSetup, *trackCollectionHandle.product(), basicClusterPtrs, the_pvtx, outputConvPhotonCollection);//allow empty basicClusterPtrs
229 
230  outputConvPhotonCollection_p->assign(outputConvPhotonCollection.begin(), outputConvPhotonCollection.end());
231  iEvent.put( outputConvPhotonCollection_p, ConvertedPhotonCollection_);
232 
233 }
234 
235 
236 
237 
240  const std::multimap<double, reco::CaloClusterPtr>& basicClusterPtrs,
241  const reco::Vertex& the_pvtx,
242  reco::ConversionCollection & outputConvPhotonCollection){
243 
244 
245  edm::ESHandle<TrackerGeometry> trackerGeomHandle;
246  edm::ESHandle<MagneticField> magFieldHandle;
247  iSetup.get<TrackerDigiGeometryRecord>().get( trackerGeomHandle );
248  iSetup.get<IdealMagneticFieldRecord>().get( magFieldHandle );
249 
250  const TrackerGeometry* trackerGeom = trackerGeomHandle.product();
251  const MagneticField* magField = magFieldHandle.product();;
252 
253  std::vector<math::XYZPoint> trackImpactPosition;
254  trackImpactPosition.reserve(allTracks.size());//track impact position at ECAL
255  std::vector<bool> trackValidECAL;//Does this track reach ECAL basic cluster (reach ECAL && match with BC)
256  trackValidECAL.assign(allTracks.size(), false);
257 
258  std::vector<reco::CaloClusterPtr> trackMatchedBC;
259  reco::CaloClusterPtr empty_bc;
260  trackMatchedBC.assign(allTracks.size(), empty_bc);//TODO find a better way to avoid copy constructor
261 
262  std::vector<int> bcHandleId;//the associated BC handle id, -1 invalid, 0 barrel 1 endcap
263  bcHandleId.assign(allTracks.size(), -1);
264 
265  std::multimap<double, int> trackInnerEta;//Track innermost state Eta map to TrackRef index, to be used in track pair sorting
266 
267 
268  ConversionHitChecker hitChecker;
269 
270 
271  //2 propagate all tracks into ECAL, record its eta and phi
272  for (reco::ConversionTrackCollection::const_iterator tk_ref = allTracks.begin(); tk_ref != allTracks.end(); ++tk_ref ){
273  const reco::Track* tk = tk_ref->trackRef().get() ;
274 
275  //map TrackRef to Eta
276  trackInnerEta.insert(std::make_pair(tk->innerMomentum().eta(), tk_ref-allTracks.begin()));
277 
278  if (allowTrackBC_){
279  //check impact position then match with BC
280  math::XYZPoint ew;
281  if ( getTrackImpactPosition(tk, trackerGeom, magField, ew) ){
282  trackImpactPosition.push_back(ew);
283 
284  reco::CaloClusterPtr closest_bc;//the closest matching BC to track
285 
286  if ( getMatchedBC(basicClusterPtrs, ew, closest_bc) ){
287  trackMatchedBC[tk_ref-allTracks.begin()] = closest_bc;
288  trackValidECAL[tk_ref-allTracks.begin()] = true;
289  bcHandleId[tk_ref-allTracks.begin()] = (fabs(closest_bc->position().eta())>1.479)?1:0;
290  }
291  } else {
292  trackImpactPosition.push_back(ew);
293  continue;
294  }
295 
296  }
297  }
298 
299 
300 
301 
302  //3. pair up : to select one track from matched EBC, and select the other track to fit cot theta and phi open angle cut
303  //TODO it is k-Closest pair of point problem
304  //std::cout << " allTracks.size() " << allTracks.size() << std::endl;
305  for(reco::ConversionTrackCollection::const_iterator ll = allTracks.begin(); ll != allTracks.end(); ++ ll ) {
306  bool track1HighPurity=true;
307  //std::cout << " Loop on allTracks " << std::endl;
308  const edm::RefToBase<reco::Track> & left = ll->trackRef();
309 
310  //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
311  //(Note that the TrackRef and GsfTrackRef versions of the constructor are needed
312  // to properly get refit tracks in the output vertex)
313  reco::TransientTrack ttk_l;
314  if (dynamic_cast<const reco::GsfTrack*>(left.get())) {
315  ttk_l = thettbuilder_->build(left.castTo<reco::GsfTrackRef>());
316  }
317  else {
318  ttk_l = thettbuilder_->build(left.castTo<reco::TrackRef>());
319  }
320 
321 
322  if ((allowTrackBC_ && !trackValidECAL[ll-allTracks.begin()]) )//this Left leg should have valid BC
323  continue;
324 
325  if (the_pvtx.isValid()){
326  if (!(trackD0Cut(left, the_pvtx))) track1HighPurity=false;
327  } else {
328  if (!(trackD0Cut(left))) track1HighPurity=false;
329  }
330 
331  std::vector<int> right_candidates;//store all right legs passed the cut (theta/approach and ref pair)
332  std::vector<double> right_candidate_theta, right_candidate_approach;
333  std::vector<std::pair<bool, reco::Vertex> > vertex_candidates;
334 
335 
336  for (reco::ConversionTrackCollection::const_iterator rr = ll+1; rr != allTracks.end(); ++ rr ) {
337  bool track2HighPurity = true;
338  bool highPurityPair = true;
339 
340  const edm::RefToBase<reco::Track> & right = rr->trackRef();
341 
342  //TODO: This is a workaround, should be fixed with a proper function in the TTBuilder
343  reco::TransientTrack ttk_r;
344  if (dynamic_cast<const reco::GsfTrack*>(right.get())) {
345  ttk_r = thettbuilder_->build(right.castTo<reco::GsfTrackRef>());
346  }
347  else {
348  ttk_r = thettbuilder_->build(right.castTo<reco::TrackRef>());
349  }
350  //std::cout << " This track is " << right->algoName() << std::endl;
351 
352 
353  //all vertexing preselection should go here
354 
355  //check for opposite charge
356  if ( allowOppCharge_ && (left->charge()*right->charge() > 0) )
357  continue; //same sign, reject pair
358 
359  if ( (allowTrackBC_ && !trackValidECAL[rr-allTracks.begin()] && rightBC_) )// if right track matches ECAL
360  continue;
361 
362 
363  double approachDist = -999.;
364  //apply preselection to track pair, unless one or both tracks are gsf
365  if (!preselectTrackPair(ttk_l,ttk_r, approachDist) && left->algo()!=29 && right->algo()!=29) {
366  continue;
367  }
368 
369  //do the actual vertex fit
370  reco::Vertex theConversionVertex;//by default it is invalid
371  bool goodVertex = checkVertex(ttk_l, ttk_r, magField, theConversionVertex);
372 
373  //bail as early as possible in case the fit didn't return a good vertex
374  if (!goodVertex) {
375  continue;
376  }
377 
378 
379 
380  //track pair pass the quality cut
381  if ( !( (trackQualityFilter(left, true) && trackQualityFilter(right, false))
382  || (trackQualityFilter(left, false) && trackQualityFilter(right, true)) ) ) {
383  highPurityPair=false;
384  }
385 
386  if (the_pvtx.isValid()){
387  if (!(trackD0Cut(right, the_pvtx))) track2HighPurity=false;
388  } else {
389  if (!(trackD0Cut(right))) track2HighPurity=false;
390  }
391 
392 
393  const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr> the_left = std::make_pair(left, trackMatchedBC[ll-allTracks.begin()]);
394  const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr> the_right = std::make_pair(right, trackMatchedBC[rr-allTracks.begin()]);
395 
396 
397 
398  //signature cuts, then check if vertex, then post-selection cuts
399  highPurityPair= highPurityPair && track1HighPurity && track2HighPurity && checkTrackPair(the_left, the_right) ;
400  highPurityPair = highPurityPair && goodVertex && checkPhi(left, right, trackerGeom, magField, theConversionVertex) ;
401 
402  //if all cuts passed, go ahead to make conversion candidates
403  std::vector<edm::RefToBase<reco::Track> > trackPairRef;
404  trackPairRef.push_back(left);//left track
405  trackPairRef.push_back(right);//right track
406 
407  std::vector<math::XYZVector> trackPin;
408  std::vector<math::XYZVector> trackPout;
409  std::vector<math::XYZPoint> trackInnPos;
410  std::vector<uint8_t> nHitsBeforeVtx;
411  std::vector<Measurement1DFloat> dlClosestHitToVtx;
412 
413  if (left->extra().isNonnull() && right->extra().isNonnull()){//only available on TrackExtra
414  trackInnPos.push_back( left->innerPosition());
415  trackInnPos.push_back( right->innerPosition());
416  trackPin.push_back(left->innerMomentum());
417  trackPin.push_back(right->innerMomentum());
418  trackPout.push_back(left->outerMomentum());
419  trackPout.push_back(right->outerMomentum());
420  }
421 
422  if (ll->trajRef().isNonnull() && rr->trajRef().isNonnull()) {
423  std::pair<uint8_t,Measurement1DFloat> leftWrongHits = hitChecker.nHitsBeforeVtx(*ll->trajRef().get(),theConversionVertex);
424  std::pair<uint8_t,Measurement1DFloat> rightWrongHits = hitChecker.nHitsBeforeVtx(*rr->trajRef().get(),theConversionVertex);
425  nHitsBeforeVtx.push_back(leftWrongHits.first);
426  nHitsBeforeVtx.push_back(rightWrongHits.first);
427  dlClosestHitToVtx.push_back(leftWrongHits.second);
428  dlClosestHitToVtx.push_back(rightWrongHits.second);
429  }
430 
431  uint8_t nSharedHits = hitChecker.nSharedHits(*left.get(),*right.get());
432 
434  //if using kinematic fit, check with chi2 post cut
435  if (theConversionVertex.isValid()){
436  const float chi2Prob = ChiSquaredProbability(theConversionVertex.chi2(), theConversionVertex.ndof());
437  if (chi2Prob<vtxChi2_) highPurityPair=false;
438  }
439 
440  //std::cout << " highPurityPair after vertex cut " << highPurityPair << std::endl;
441  std::vector<math::XYZPoint> trkPositionAtEcal;
442  std::vector<reco::CaloClusterPtr> matchingBC;
443 
444  if (allowTrackBC_){//TODO find out the BC ptrs if not doing matching, otherwise, leave it empty
445  const int lbc_handle = bcHandleId[ll-allTracks.begin()],
446  rbc_handle = bcHandleId[rr-allTracks.begin()];
447 
448  trkPositionAtEcal.push_back(trackImpactPosition[ll-allTracks.begin()]);//left track
449  if (trackValidECAL[rr-allTracks.begin()])//second track ECAL position may be invalid
450  trkPositionAtEcal.push_back(trackImpactPosition[rr-allTracks.begin()]);
451 
452  matchingBC.push_back(trackMatchedBC[ll-allTracks.begin()]);//left track
453  scPtrVec.push_back(trackMatchedBC[ll-allTracks.begin()]);//left track
454  if (trackValidECAL[rr-allTracks.begin()]){//second track ECAL position may be invalid
455  matchingBC.push_back(trackMatchedBC[rr-allTracks.begin()]);
456  if (!(trackMatchedBC[rr-allTracks.begin()] == trackMatchedBC[ll-allTracks.begin()])//avoid 1 bc 2 tk
457  && lbc_handle == rbc_handle )//avoid ptr from different collection
458  scPtrVec.push_back(trackMatchedBC[rr-allTracks.begin()]);
459  }
460  }
461  const float minAppDist = approachDist;
462 
464  float dummy=0;
465  reco::Conversion newCandidate(scPtrVec, trackPairRef, trkPositionAtEcal, theConversionVertex, matchingBC, minAppDist, trackInnPos, trackPin, trackPout, nHitsBeforeVtx, dlClosestHitToVtx, nSharedHits, dummy, algo );
466  newCandidate.setQuality(reco::Conversion::highPurity, highPurityPair);
467 
468  bool generalTracksOnly = ll->isTrackerOnly() && rr->isTrackerOnly() && !dynamic_cast<const reco::GsfTrack*>(ll->trackRef().get()) && !dynamic_cast<const reco::GsfTrack*>(rr->trackRef().get());
469  bool arbitratedEcalSeeded = ll->isArbitratedEcalSeeded() && rr->isArbitratedEcalSeeded();
470  bool arbitratedMerged = ll->isArbitratedMerged() && rr->isArbitratedMerged();
471  bool arbitratedMergedEcalGeneral = ll->isArbitratedMergedEcalGeneral() && rr->isArbitratedMergedEcalGeneral();
472 
473  newCandidate.setQuality(reco::Conversion::generalTracksOnly, generalTracksOnly);
474  newCandidate.setQuality(reco::Conversion::arbitratedEcalSeeded, arbitratedEcalSeeded);
475  newCandidate.setQuality(reco::Conversion::arbitratedMerged, arbitratedMerged);
476  newCandidate.setQuality(reco::Conversion::arbitratedMergedEcalGeneral, arbitratedMergedEcalGeneral);
477 
478  outputConvPhotonCollection.push_back(newCandidate);
479 
480  }
481 
482  }
483 
484 
485 
486 
487 
488 
489 }
490 
491 
492 
493 
494 
495 //
496 // member functions
497 //
498 
500  bool pass = true;
501  if (isLeft){
502  pass = (ref->normalizedChi2() < maxChi2Left_ && ref->found() >= minHitsLeft_);
503  } else {
504  pass = (ref->normalizedChi2() < maxChi2Right_ && ref->found() >= minHitsRight_);
505  }
506 
507  return pass;
508 }
509 
511  //NOTE if not allow d0 cut, always true
512  return ((!allowD0_) || !(ref->d0()*ref->charge()/ref->d0Error()<d0Cut_));
513 }
514 
516  //
517  return ((!allowD0_) || !(-ref->dxy(the_pvtx.position())*ref->charge()/ref->dxyError()<d0Cut_));
518 }
519 
520 
522  const TrackerGeometry* trackerGeom, const MagneticField* magField,
523  math::XYZPoint& ew){
524 
525  PropagatorWithMaterial propag( alongMomentum, 0.000511, magField );
526  TrajectoryStateTransform transformer;
528  new BoundCylinder( GlobalPoint(0.,0.,0.), TkRotation<float>(),
529  SimpleCylinderBounds( 129, 129, -320.5, 320.5 ) ) );
530  const float epsilon = 0.001;
531  Surface::RotationType rot; // unit rotation matrix
532  const float barrelRadius = 129.f;
533  const float barrelHalfLength = 270.9f;
534  const float endcapRadius = 171.1f;
535  const float endcapZ = 320.5f;
537  SimpleCylinderBounds( barrelRadius-epsilon, barrelRadius+epsilon, -barrelHalfLength, barrelHalfLength)));
538  ReferenceCountingPointer<BoundDisk> theNegativeEtaEndcap_(
539  new BoundDisk( Surface::PositionType( 0, 0, -endcapZ), rot,
540  SimpleDiskBounds( 0, endcapRadius, -epsilon, epsilon)));
541  ReferenceCountingPointer<BoundDisk> thePositiveEtaEndcap_(
542  new BoundDisk( Surface::PositionType( 0, 0, endcapZ), rot,
543  SimpleDiskBounds( 0, endcapRadius, -epsilon, epsilon)));
544 
545  //const TrajectoryStateOnSurface myTSOS = transformer.innerStateOnSurface(*(*ref), *trackerGeom, magField);
546  const TrajectoryStateOnSurface myTSOS = transformer.outerStateOnSurface(*tk_ref, *trackerGeom, magField);
547  TrajectoryStateOnSurface stateAtECAL;
548  stateAtECAL = propag.propagate(myTSOS, *theBarrel_);
549  if (!stateAtECAL.isValid() || ( stateAtECAL.isValid() && fabs(stateAtECAL.globalPosition().eta() ) >1.479 ) ) {
550  //endcap propagator
551  if (myTSOS.globalPosition().eta() > 0.) {
552  stateAtECAL = propag.propagate(myTSOS, *thePositiveEtaEndcap_);
553  } else {
554  stateAtECAL = propag.propagate(myTSOS, *theNegativeEtaEndcap_);
555  }
556  }
557  if (stateAtECAL.isValid()){
558  ew = stateAtECAL.globalPosition();
559  return true;
560  }
561  else
562  return false;
563 }
564 
565 bool TrackerOnlyConversionProducer::getMatchedBC(const std::multimap<double, reco::CaloClusterPtr>& bcMap,
566  const math::XYZPoint& trackImpactPosition,
567  reco::CaloClusterPtr& closestBC){
568  const double track_eta = trackImpactPosition.eta();
569  const double track_phi = trackImpactPosition.phi();
570 
571  double min_eta = 999., min_phi = 999.;
572  reco::CaloClusterPtr closest_bc;
573  for (std::multimap<double, reco::CaloClusterPtr>::const_iterator bc = bcMap.lower_bound(track_eta - halfWayEta_);
574  bc != bcMap.upper_bound(track_eta + halfWayEta_); ++bc){//use eta map to select possible BC collection then loop in
575  const reco::CaloClusterPtr& ebc = bc->second;
576  const double delta_eta = track_eta-(ebc->position().eta());
577  const double delta_phi = reco::deltaPhi(track_phi, (ebc->position().phi()));
578  if (fabs(delta_eta)<dEtaTkBC_ && fabs(delta_phi)<dPhiTkBC_){
579  if (fabs(min_eta)>fabs(delta_eta) && fabs(min_phi)>fabs(delta_phi)){//take the closest to track BC
580  min_eta = delta_eta;
581  min_phi = delta_phi;
582  closest_bc = bc->second;
583  //TODO check if min_eta>delta_eta but min_phi<delta_phi
584  }
585  }
586  }
587 
588  if (min_eta < 999.){
589  closestBC = closest_bc;
590  return true;
591  } else
592  return false;
593 }
594 
596  const math::XYZPoint& trackImpactPosition,
597  reco::CaloClusterPtr& closestBC){
598  const double track_eta = trackImpactPosition.eta();
599  const double track_phi = trackImpactPosition.phi();
600 
601  double min_eta = 999., min_phi = 999.;
602  reco::CaloClusterPtr closest_bc;
604  bc != bcMap.end(); ++bc){//use eta map to select possible BC collection then loop in
605  const reco::CaloClusterPtr& ebc = (*bc);
606  const double delta_eta = track_eta-(ebc->position().eta());
607  const double delta_phi = reco::deltaPhi(track_phi, (ebc->position().phi()));
608  if (fabs(delta_eta)<dEtaTkBC_ && fabs(delta_phi)<dPhiTkBC_){
609  if (fabs(min_eta)>fabs(delta_eta) && fabs(min_phi)>fabs(delta_phi)){//take the closest to track BC
610  min_eta = delta_eta;
611  min_phi = delta_phi;
612  closest_bc = ebc;
613  //TODO check if min_eta>delta_eta but min_phi<delta_phi
614  }
615  }
616  }
617 
618  if (min_eta < 999.){
619  closestBC = closest_bc;
620  return true;
621  } else
622  return false;
623 }
624 
625 //check track open angle of phi at vertex
627  const TrackerGeometry* trackerGeom, const MagneticField* magField,
628  const reco::Vertex& vtx){
629  if (!allowDeltaPhi_)
630  return true;
631  //if track has innermost momentum, check with innermost phi
632  //if track also has valid vertex, propagate to vertex then calculate phi there
633  //if track has no innermost momentum, just return true, because track->phi() makes no sense
634  if (tk_l->extra().isNonnull() && tk_r->extra().isNonnull()){
635  double iphi1 = tk_l->innerMomentum().phi(), iphi2 = tk_r->innerMomentum().phi();
636  if (vtx.isValid()){
637  PropagatorWithMaterial propag( anyDirection, 0.000511, magField );
638  TrajectoryStateTransform transformer;
639  double recoPhoR = vtx.position().Rho();
642  SimpleCylinderBounds( recoPhoR-0.001, recoPhoR+0.001, -fabs(vtx.position().z()), fabs(vtx.position().z()))));
644  new BoundDisk( Surface::PositionType( 0, 0, vtx.position().z()), rot,
645  SimpleDiskBounds( 0, recoPhoR, -0.001, 0.001)));
646 
647  const TrajectoryStateOnSurface myTSOS1 = transformer.innerStateOnSurface(*tk_l, *trackerGeom, magField);
648  const TrajectoryStateOnSurface myTSOS2 = transformer.innerStateOnSurface(*tk_r, *trackerGeom, magField);
649  TrajectoryStateOnSurface stateAtVtx1, stateAtVtx2;
650  stateAtVtx1 = propag.propagate(myTSOS1, *theBarrel_);
651  if (!stateAtVtx1.isValid() ) {
652  stateAtVtx1 = propag.propagate(myTSOS1, *theDisk_);
653  }
654  if (stateAtVtx1.isValid()){
655  iphi1 = stateAtVtx1.globalDirection().phi();
656  }
657  stateAtVtx2 = propag.propagate(myTSOS2, *theBarrel_);
658  if (!stateAtVtx2.isValid() ) {
659  stateAtVtx2 = propag.propagate(myTSOS2, *theDisk_);
660  }
661  if (stateAtVtx2.isValid()){
662  iphi2 = stateAtVtx2.globalDirection().phi();
663  }
664  }
665  const double dPhi = reco::deltaPhi(iphi1, iphi2);
666  return (fabs(dPhi) < deltaPhi_);
667  } else {
668  return true;
669  }
670 }
671 
673  double& appDist) {
674 
675 
676  double dCotTheta = 1./tan(ttk_l.track().innerMomentum().theta()) - 1./tan(ttk_r.track().innerMomentum().theta());
677  if (allowDeltaCot_ && (std::abs(dCotTheta) > deltaCotTheta_)) {
678  return false;
679  }
680 
681  //non-conversion hypothesis, reject prompt track pairs
682  ClosestApproachInRPhi closest;
684  if (!closest.status()) {
685  return false;
686  }
687 
688  if (closest.crossingPoint().perp() < r_cut) {
689  return false;
690  }
691 
692 
693  //compute tangent point btw tracks (conversion hypothesis)
694  TangentApproachInRPhi tangent;
696  if (!tangent.status()) {
697  return false;
698  }
699 
700  GlobalPoint tangentPoint = tangent.crossingPoint();
701  double rho = tangentPoint.perp();
702 
703  //reject candidates well outside of tracker bounds
704  if (rho > 120.0) {
705  return false;
706  }
707 
708  if (std::abs(tangentPoint.z()) > 300.0) {
709  return false;
710  }
711 
712  std::pair<GlobalTrajectoryParameters,GlobalTrajectoryParameters> trajs = tangent.trajectoryParameters();
713 
714  //very large separation in z, no hope
715  if (std::abs(trajs.first.position().z() - trajs.second.position().z()) > dzCut_) {
716  return false;
717  }
718 
719 
720  float minApproach = tangent.perpdist();
721  appDist = minApproach;
722 
723  if (allowMinApproach_ && (minApproach < minApproachLow_ || minApproach > minApproachHigh_) ) {
724  return false;
725  }
726 
727  return true;
728 
729 
730 }
731 
733  const std::pair<edm::RefToBase<reco::Track>, reco::CaloClusterPtr>& rr){
734 
735  const reco::CaloClusterPtr& bc_l = ll.second;//can be null, so check isNonnull()
736  const reco::CaloClusterPtr& bc_r = rr.second;
737 
738  //The cuts should be ordered by considering if takes time and if cuts off many fakes
739  if (allowTrackBC_){
740  //check energy of BC
741  double total_e_bc = 0;
742  if (bc_l.isNonnull()) total_e_bc += bc_l->energy();
743  if (rightBC_)
744  if (bc_r.isNonnull())
745  total_e_bc += bc_r->energy();
746 
747  if (total_e_bc < energyTotalBC_) return false;
748  }
749 
750  return true;
751 }
752 
753 
754 
755 //because reco::vertex uses track ref, so have to keep them
757  const MagneticField* magField,
758  reco::Vertex& the_vertex){
759  bool found = false;
760 
761  std::vector<reco::TransientTrack> pair;
762  pair.push_back(ttk_l);
763  pair.push_back(ttk_r);
764 
765  found = theVertexFinder_->run(pair, the_vertex);
766 
767 
768 
769  return found;
770 }
771 
772 
T getParameter(std::string const &) const
double d0Error() const
error on d0
Definition: TrackBase.h:204
double delta_eta(double eta1, double eta2)
Definition: AnglesUtil.h:98
T perp() const
Definition: PV3DBase.h:66
double d0() const
dxy parameter in perigee convention (d0 = - dxy)
Definition: TrackBase.h:122
const TrackExtraRef & extra() const
reference to &quot;extra&quot; object
Definition: Track.h:97
std::pair< GlobalTrajectoryParameters, GlobalTrajectoryParameters > trajectoryParameters() const
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:110
double dxyError() const
error on dxy
Definition: TrackBase.h:202
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:135
std::pair< uint8_t, Measurement1DFloat > nHitsBeforeVtx(const Trajectory &traj, const reco::Vertex &vtx, double sigmaTolerance=3.0) const
void setQuality(ConversionQuality q, bool b)
Definition: Conversion.h:233
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:61
Definition: DDAxes.h:10
bool checkPhi(const edm::RefToBase< reco::Track > &tk_l, const edm::RefToBase< reco::Track > &tk_r, const TrackerGeometry *trackerGeom, const MagneticField *magField, const reco::Vertex &the_vertex)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
reco::TransientTrack build(const reco::Track *p) const
#define abs(x)
Definition: mlp_lapack.h:159
static ConversionAlgorithm algoByName(const std::string &name)
Definition: Conversion.cc:167
GlobalPoint globalPosition() const
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
const Point & position() const
position
Definition: Vertex.h:93
std::vector< Conversion > ConversionCollection
collectin of Conversion objects
Definition: ConversionFwd.h:9
bool trackQualityFilter(const edm::RefToBase< reco::Track > &ref, bool isLeft)
TrajectoryStateOnSurface innermostMeasurementState() const
virtual GlobalPoint crossingPoint() const
virtual GlobalPoint crossingPoint() const
const_iterator begin() const
Definition: PtrVector.h:124
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
TrackAlgorithm algo() const
Definition: TrackBase.h:303
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:248
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:156
std::vector< ConversionTrack > ConversionTrackCollection
collection of ConversionTracks
int iEvent
Definition: GenABIO.cc:243
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
TransientVertex run(std::vector< reco::TransientTrack > pair)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
T z() const
Definition: PV3DBase.h:58
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
const_iterator end() const
Definition: PtrVector.h:129
double chi2() const
chi-squares
Definition: Vertex.h:82
TrajectoryStateOnSurface outerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field) const
float ChiSquaredProbability(double chiSquared, double nrDOF)
virtual void produce(edm::Event &, const edm::EventSetup &)
TrackerOnlyConversionProducer(const edm::ParameterSet &)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
uint8_t nSharedHits(const reco::Track &trk1, const reco::Track &trk2) const
double ndof() const
Definition: Vertex.h:89
double delta_phi(double ph11, double phi2)
Definition: AnglesUtil.h:91
virtual TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface &tsos, const Plane &plane) const
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
bool getTrackImpactPosition(const reco::Track *tk_ref, const TrackerGeometry *trackerGeom, const MagneticField *magField, math::XYZPoint &ew)
bool trackD0Cut(const edm::RefToBase< reco::Track > &ref)
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:49
virtual bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb)
const Track & track() const
const TransientTrackBuilder * thettbuilder_
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
const T & get() const
Definition: EventSetup.h:55
REF castTo() const
cast to a concrete type
Definition: RefToBase.h:236
T const * product() const
Definition: ESHandle.h:62
bool getMatchedBC(const std::multimap< double, reco::CaloClusterPtr > &bcMap, const math::XYZPoint &trackImpactPosition, reco::CaloClusterPtr &closestBC)
bool preselectTrackPair(const reco::TransientTrack &ttk_l, const reco::TransientTrack &ttk_r, double &appDist)
TrajectoryStateOnSurface innerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field) const
T eta() const
Definition: PV3DBase.h:70
virtual bool status() const
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:45
bool checkTrackPair(const std::pair< edm::RefToBase< reco::Track >, reco::CaloClusterPtr > &ll, const std::pair< edm::RefToBase< reco::Track >, reco::CaloClusterPtr > &rr)
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:100
int charge() const
track electric charge
Definition: TrackBase.h:112
bool checkVertex(const reco::TransientTrack &ttk_l, const reco::TransientTrack &ttk_r, const MagneticField *magField, reco::Vertex &the_vertex)
const double epsilon
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:120
void buildCollection(edm::Event &iEvent, const edm::EventSetup &iSetup, const reco::ConversionTrackCollection &allTracks, const std::multimap< double, reco::CaloClusterPtr > &basicClusterPtrs, const reco::Vertex &the_pvtx, reco::ConversionCollection &outputConvPhotonCollection)
virtual bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb)
value_type const * get() const
Definition: RefToBase.h:207
virtual bool status() const
GlobalVector globalDirection() const