PFTracking/src/PFDisplacedVertexCandidateFinder.cc

0001 #include <memory>
0002
0003 #include "RecoParticleFlow/PFTracking/interface/PFDisplacedVertexCandidateFinder.h"
0004
0005 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011
0012 using namespace std;
0013 using namespace reco;
0014
0015 //for debug only
0016 //#define PFLOW_DEBUG
0017
0018 PFDisplacedVertexCandidateFinder::PFDisplacedVertexCandidateFinder()
0019     : vertexCandidates_(new PFDisplacedVertexCandidateCollection),
0020       dcaCut_(1000),
0021       primaryVertexCut2_(0.0),
0022       dcaPInnerHitCut2_(1000.0),
0023       vertexCandidatesSize_(50),
0024       theMinimum_(TwoTrackMinimumDistance::SlowMode),
0025       debug_(false),
0026       magField_(nullptr) {}
0027
0028 PFDisplacedVertexCandidateFinder::~PFDisplacedVertexCandidateFinder() {
0029   LogDebug("PFDisplacedVertexCandidateFinder")
0030       << "~PFDisplacedVertexCandidateFinder - number of remaining elements: " << eventTracks_.size();
0031 }
0032
0033 void PFDisplacedVertexCandidateFinder::setPrimaryVertex(edm::Handle<reco::VertexCollection> mainVertexHandle,
0034                                                         edm::Handle<reco::BeamSpot> beamSpotHandle) {
0035   const math::XYZPoint beamSpot = beamSpotHandle.isValid()
0036                                       ? math::XYZPoint(beamSpotHandle->x0(), beamSpotHandle->y0(), beamSpotHandle->z0())
0037                                       : math::XYZPoint(0, 0, 0);
0038
0039   // The primary vertex is taken from the refitted list,
0040   // if does not exist from the average offline beam spot position
0041   // if does not exist (0,0,0) is used
0042   pvtx_ = mainVertexHandle.isValid()
0043               ? math::XYZPoint(
0044                     mainVertexHandle->begin()->x(), mainVertexHandle->begin()->y(), mainVertexHandle->begin()->z())
0045               : beamSpot;
0046 }
0047
0048 // Set the imput collection of tracks and calculate their
0049 // trajectory parameters the Global Trajectory Parameters
0050 void PFDisplacedVertexCandidateFinder::setInput(const edm::Handle<TrackCollection>& trackh,
0051                                                 const MagneticField* magField) {
0052   magField_ = magField;
0053   trackMask_.clear();
0054   trackMask_.resize(trackh->size());
0055   eventTrackTrajectories_.clear();
0056   eventTrackTrajectories_.resize(trackh->size());
0057
0058   for (unsigned i = 0; i < trackMask_.size(); i++)
0059     trackMask_[i] = true;
0060
0061   eventTracks_.clear();
0062   if (trackh.isValid()) {
0063     for (unsigned i = 0; i < trackh->size(); i++) {
0064       TrackRef tref(trackh, i);
0065       TrackBaseRef tbref(tref);
0066
0067       if (!isSelected(tbref)) {
0068         trackMask_[i] = false;
0069         continue;
0070       }
0071
0072       const Track* trk = tref.get();
0073       eventTracks_.push_back(tbref);
0074       eventTrackTrajectories_[i] = getGlobalTrajectoryParameters(trk);
0075     }
0076   }
0077   track_table_ = edm::soa::PtEtaPhiTable{*trackh};
0078 }
0079
0080 // -------- Main function which find vertices -------- //
0081
0082 void PFDisplacedVertexCandidateFinder::findDisplacedVertexCandidates() {
0083   LogDebug("PFDisplacedVertexCandidateFinder") << "========= Start Finding Displaced Vertex Candidates =========";
0084
0085   // The vertexCandidates have not been passed to the event, and need to be cleared
0086   if (vertexCandidates_.get())
0087     vertexCandidates_->clear();
0088   else
0089     vertexCandidates_ = std::make_unique<PFDisplacedVertexCandidateCollection>();
0090
0091   vertexCandidates_->reserve(vertexCandidatesSize_);
0092   for (IE ie = eventTracks_.begin(); ie != eventTracks_.end();) {
0093     // Run the recursive procedure to find all tracks link together
0094     // In one blob called Candidate
0095
0096     PFDisplacedVertexCandidate tempVertexCandidate;
0097
0098     ie = associate(eventTracks_.end(), ie, tempVertexCandidate);
0099     LogDebug("PFDisplacedVertexCandidateFinder") << "elements=" << tempVertexCandidate.elements().size();
0100
0101     // Build remaining links in current block
0102     if (tempVertexCandidate.isValid()) {
0103       packLinks(tempVertexCandidate);
0104       vertexCandidates_->push_back(tempVertexCandidate);
0105     }
0106   }
0107
0108   LogDebug("PFDisplacedVertexCandidateFinder") << "========= End Finding Displaced Vertex Candidates =========";
0109 }
0110
0111 PFDisplacedVertexCandidateFinder::IE PFDisplacedVertexCandidateFinder::associate(
0112     IE last, IE next, PFDisplacedVertexCandidate& tempVertexCandidate) {
0113   LogDebug("PFDisplacedVertexCandidateFinder") << "== Start the association procedure ==";
0114
0115   if (last != eventTracks_.end()) {
0116     double dist = -1;
0117     GlobalPoint crossingPoint(0, 0, 0);
0118     PFDisplacedVertexCandidate::VertexLinkTest linktest;
0119     link(*last, *next, dist, crossingPoint, linktest);
0120
0121     if (dist < -0.5) {
0122       LogDebug("PFDisplacedVertexCandidateFinder") << "link failed";
0123
0124       return ++next;  // association failed
0125     } else {
0126       // add next element to the current pflowblock
0127       LogDebug("PFDisplacedVertexCandidateFinder")
0128           << "adding element " << (*next).key() << " to PFDisplacedVertexCandidate with "
0129           << tempVertexCandidate.elements().size() << " elements";
0130       tempVertexCandidate.addElement((*next));
0131       trackMask_[(*next).key()] = false;
0132
0133       LogDebug("PFDisplacedVertexCandidateFinder")
0134           << "link parameters "
0135           << " *next = " << (*next).key() << " *last = " << (*last).key() << "  dist = " << crossingPoint
0136           << " P.x = " << crossingPoint.x() << " P.y = " << crossingPoint.y() << " P.z = " << crossingPoint.z();
0137     }
0138   } else {
0139     // add next element to this eflowblock
0140     LogDebug("PFDisplacedVertexCandidateFinder") << "adding to block element " << (*next).key();
0141     tempVertexCandidate.addElement((*next));
0142     trackMask_[(*next).key()] = false;
0143   }
0144
0145   // recursive call: associate next and other unused elements
0146 #ifdef EDM_ML_DEBUG
0147   for (unsigned i = 0; i < trackMask_.size(); i++)
0148     LogTrace("PFDisplacedVertexCandidateFinder") << " Mask[" << i << "] = " << trackMask_[i];
0149 #endif
0150
0151   for (IE ie = eventTracks_.begin(); ie != eventTracks_.end();) {
0152     if (ie == last || ie == next) {
0153       ++ie;
0154       continue;
0155     }
0156
0157     // *ie already included to a block
0158     if (!trackMask_[(*ie).key()]) {
0159       ++ie;
0160       continue;
0161     }
0162
0163     LogDebug("PFDisplacedVertexCandidateFinder") << "calling associate " << (*next).key() << " & " << (*ie).key();
0164     ie = associate(next, ie, tempVertexCandidate);
0165   }
0166
0167   LogDebug("PFDisplacedVertexCandidateFinder") << "**** removing element ";
0168
0169   IE iteratorToNextFreeElement = eventTracks_.erase(next);
0170
0171   LogDebug("PFDisplacedVertexCandidateFinder") << "== End the association procedure ==";
0172
0173   return iteratorToNextFreeElement;
0174 }
0175
0176 void PFDisplacedVertexCandidateFinder::link(const TrackBaseRef& el1,
0177                                             const TrackBaseRef& el2,
0178                                             double& dist,
0179                                             GlobalPoint& crossingPoint,
0180                                             PFDisplacedVertexCandidate::VertexLinkTest& vertexLinkTest) {
0181   using namespace edm::soa::col;
0182   const auto iel1 = el1.key();
0183   const auto iel2 = el2.key();
0184   const auto pt1 = track_table_.get<Pt>(iel1);
0185   const auto pt2 = track_table_.get<Pt>(iel2);
0186   const auto eta1 = track_table_.get<Eta>(iel1);
0187   const auto eta2 = track_table_.get<Eta>(iel2);
0188   const auto phi1 = track_table_.get<Phi>(iel1);
0189   const auto phi2 = track_table_.get<Phi>(iel2);
0190
0191   if (std::abs(eta1 - eta2) > 1) {
0192     dist = -1;
0193     return;
0194   }
0195   if (pt1 > 2 && pt2 > 2 && std::abs(::deltaPhi(phi1, phi2)) > 1) {
0196     dist = -1;
0197     return;
0198   }
0199
0200   const GlobalTrajectoryParameters& gt1 = eventTrackTrajectories_[iel1];
0201   const GlobalTrajectoryParameters& gt2 = eventTrackTrajectories_[iel2];
0202
0203   // Closest approach
0204   theMinimum_.calculate(gt1, gt2);
0205
0206   // Fill the parameters
0207   dist = theMinimum_.distance();
0208   LogDebug("PFDisplacedVertexCandidateFinder") << "link iel1=" << iel1 << " iel2=" << iel2 << " dist=" << dist;
0209   crossingPoint = theMinimum_.crossingPoint();
0210
0211   vertexLinkTest = PFDisplacedVertexCandidate::LINKTEST_DCA;  //rechit by default
0212
0213   // Check if the link test fails
0214   if (dist > dcaCut_) {
0215     dist = -1;
0216     return;
0217   }
0218
0219   // Check if the closses approach point is too close to the primary vertex/beam pipe
0220   double rho2 = crossingPoint.x() * crossingPoint.x() + crossingPoint.y() * crossingPoint.y();
0221
0222   if (rho2 < primaryVertexCut2_) {
0223     dist = -1;
0224     return;
0225   }
0226
0227   // Check if the closses approach point is out of range of the tracker
0228
0229   double tob_rho_limit2 = 10000;
0230   double tec_z_limit = 270;
0231
0232   if (rho2 > tob_rho_limit2) {
0233     dist = -1;
0234     return;
0235   }
0236   if (std::abs(crossingPoint.z()) > tec_z_limit) {
0237     dist = -1;
0238     return;
0239   }
0240
0241   return;
0242
0243   /*
0244   // Check if the inner hit of one of the tracks is not too far away from the vertex
0245   double dx = el1->innerPosition().x()-P.x();
0246   double dy = el1->innerPosition().y()-P.y();
0247   double dz = el1->innerPosition().z()-P.z();
0248   double dist2 = dx*dx+dy*dy+dz*dz;
0249
0250   if (dist2 > dcaPInnerHitCut2_) {
0251     dist = -1;
0252     if (debug_) cout << "track " << el1.key() << " dist to vertex " << sqrt(dist2) << endl;
0253     return;
0254   }
0255
0256   dx = el2->innerPosition().x()-P.x();
0257   dy = el2->innerPosition().y()-P.y();
0258   dz = el2->innerPosition().z()-P.z();
0259   dist2 = dx*dx+dy*dy+dz*dz;
0260
0261   if (dist2 > dcaPInnerHitCut2_) {
0262     dist = -1;
0263     if (debug_) cout << "track " << el2.key() << " dist to vertex " << sqrt(dist2) << endl;
0264     return;
0265   }
0266   */
0267 }
0268
0269 // Build up a matrix of all the links between different tracks in
0270 // In the Candidate
0271 void PFDisplacedVertexCandidateFinder::packLinks(PFDisplacedVertexCandidate& vertexCandidate) {
0272   const vector<TrackBaseRef>& els = vertexCandidate.elements();
0273
0274   //First Loop: update all link data
0275   for (unsigned i1 = 0; i1 < els.size(); i1++) {
0276     for (unsigned i2 = i1 + 1; i2 < els.size(); i2++) {
0277       // no reflexive link
0278       if (i1 == i2)
0279         continue;
0280
0281       double dist = -1;
0282       GlobalPoint crossingPoint(0, 0, 0);
0283       PFDisplacedVertexCandidate::VertexLinkTest linktest;
0284
0285       link(els[i1], els[i2], dist, crossingPoint, linktest);
0286       LogDebug("PFDisplacedVertexCandidateFinder")
0287           << "Setting link between elements " << i1 << " key " << els[i1].key() << " and " << i2 << " key "
0288           << els[i2].key() << " of dist =" << dist << " computed from link test " << linktest;
0289
0290       if (dist > -0.5)
0291         vertexCandidate.setLink(i1, i2, dist, crossingPoint, linktest);
0292     }
0293   }
0294 }
0295
0296 // --------------- TOOLS -------------- //
0297
0298 // This tool is a copy from VZeroFinder
0299 GlobalTrajectoryParameters PFDisplacedVertexCandidateFinder::getGlobalTrajectoryParameters(const Track* track) const {
0300   const GlobalPoint position(track->vx(), track->vy(), track->vz());
0301
0302   const GlobalVector momentum(track->momentum().x(), track->momentum().y(), track->momentum().z());
0303
0304   GlobalTrajectoryParameters gtp(position, momentum, track->charge(), magField_);
0305
0306   return gtp;
0307 }
0308
0309 // This tool allow to pre-select the track for the displaced vertex finder
0310 bool PFDisplacedVertexCandidateFinder::goodPtResolution(const TrackBaseRef& trackref) const {
0311   double nChi2 = trackref->normalizedChi2();
0312   double pt = trackref->pt();
0313   double dpt = trackref->ptError();
0314   double dxy = trackref->dxy(pvtx_);
0315
0316   double pt_error = dpt / pt * 100;
0317
0318   LogDebug("PFDisplacedVertexCandidateFinder")
0319       << " PFDisplacedVertexFinder: PFrecTrack->Track Pt= " << pt << " dPt/Pt = " << pt_error << "% nChi2 = " << nChi2;
0320   if (nChi2 > nChi2_max_ || pt < pt_min_) {
0321     LogDebug("PFDisplacedVertexCandidateFinder") << " PFBlockAlgo: skip badly measured or low pt track"
0322                                                  << " nChi2_cut = " << 5 << " pt_cut = " << 0.2;
0323     return false;
0324   }
0325   //  cout << "dxy = " << dxy << endl;
0326   if (std::abs(dxy) < dxy_ && pt < pt_min_prim_)
0327     return false;
0328   //  if (fabs(dxy) < 0.2 && pt < 0.8) return false;
0329
0330   return true;
0331 }
0332
0333 ostream& operator<<(std::ostream& out, const PFDisplacedVertexCandidateFinder& a) {
0334   if (!out)
0335     return out;
0336
0337   out << "====== Particle Flow Block Algorithm ======= ";
0338   out << endl;
0339   out << "number of unassociated elements : " << a.eventTracks_.size() << endl;
0340   out << endl;
0341
0342   out << " Tracks selection based on " << std::endl;
0343   out << " pvtx_ = " << a.pvtx_ << std::endl;
0344   out << " std::abs(dxy) < " << a.dxy_ << " and pt < " << a.pt_min_prim_ << std::endl;
0345   out << " nChi2 < " << a.nChi2_max_ << " and pt < " << a.pt_min_ << std::endl;
0346
0347   out << endl;
0348
0349   for (PFDisplacedVertexCandidateFinder::IEC ie = a.eventTracks_.begin(); ie != a.eventTracks_.end(); ie++) {
0350     double pt = (*ie).get()->pt();
0351
0352     math::XYZPoint Pi = (*ie).get()->innerPosition();
0353     math::XYZPoint Po = (*ie).get()->outerPosition();
0354
0355     double innermost_radius = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y() + Pi.z() * Pi.z());
0356     double outermost_radius = sqrt(Po.x() * Po.x() + Po.y() * Po.y() + Po.z() * Po.z());
0357     double innermost_rho = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y());
0358     double outermost_rho = sqrt(Po.x() * Po.x() + Po.y() * Po.y());
0359
0360     out << "ie = " << (*ie).key() << " pt = " << pt << " innermost hit radius = " << innermost_radius
0361         << " rho = " << innermost_rho << " outermost hit radius = " << outermost_radius << " rho = " << outermost_rho
0362         << endl;
0363   }
0364
0365   const std::unique_ptr<reco::PFDisplacedVertexCandidateCollection>& vertexCandidates = std::move(a.vertexCandidates());
0366
0367   if (!vertexCandidates.get()) {
0368     out << "vertexCandidates already transfered" << endl;
0369   } else {
0370     out << "number of vertexCandidates : " << vertexCandidates->size() << endl;
0371     out << endl;
0372
0373     for (PFDisplacedVertexCandidateFinder::IBC ib = vertexCandidates->begin(); ib != vertexCandidates->end(); ib++)
0374       ib->Dump();
0375   }
0376
0377   return out;
0378 }