PFTracking/src/PFDisplacedVertexFinder.cc

0001 #include <memory>
0002
0003 #include "RecoParticleFlow/PFTracking/interface/PFDisplacedVertexFinder.h"
0004 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0005
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include "DataFormats/TrackReco/interface/Track.h"
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010
0011 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0012 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0013 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0014
0015 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0016
0017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0018 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0019
0020 #include "MagneticField/Engine/interface/MagneticField.h"
0021
0022 #include "TMath.h"
0023
0024 using namespace std;
0025 using namespace reco;
0026
0027 //for debug only
0028 //#define PFLOW_DEBUG
0029
0030 PFDisplacedVertexFinder::PFDisplacedVertexFinder()
0031     : displacedVertexCandidates_(nullptr),
0032       displacedVertices_(new PFDisplacedVertexCollection),
0033       transvSize_(0.0),
0034       longSize_(0.0),
0035       primaryVertexCut_(0.0),
0036       tobCut_(0.0),
0037       tecCut_(0.0),
0038       minAdaptWeight_(2.0),
0039       debug_(false) {}
0040
0041 PFDisplacedVertexFinder::~PFDisplacedVertexFinder() {}
0042
0043 void PFDisplacedVertexFinder::setInput(
0044     const edm::Handle<reco::PFDisplacedVertexCandidateCollection>& displacedVertexCandidates) {
0045   if (displacedVertexCandidates.isValid()) {
0046     displacedVertexCandidates_ = displacedVertexCandidates.product();
0047   } else {
0048     displacedVertexCandidates_ = nullptr;
0049   }
0050 }
0051
0052 // -------- Main function which find vertices -------- //
0053
0054 void PFDisplacedVertexFinder::findDisplacedVertices() {
0055   if (debug_)
0056     cout << "========= Start Find Displaced Vertices =========" << endl;
0057
0058   // The vertexCandidates have not been passed to the event
0059   // So they need to be cleared is they are not empty
0060   if (displacedVertices_.get())
0061     displacedVertices_->clear();
0062   else
0063     displacedVertices_ = std::make_unique<PFDisplacedVertexCollection>();
0064
0065   if (displacedVertexCandidates_ == nullptr) {
0066     edm::LogInfo("EmptyVertexInput")
0067         << "displacedVertexCandidates are not set or the setInput was called with invalid vertex";
0068     return;
0069   }
0070
0071   // Prepare the collections
0072   PFDisplacedVertexSeedCollection tempDisplacedVertexSeeds;
0073   tempDisplacedVertexSeeds.reserve(4 * displacedVertexCandidates_->size());
0074   PFDisplacedVertexCollection tempDisplacedVertices;
0075   tempDisplacedVertices.reserve(4 * displacedVertexCandidates_->size());
0076
0077   if (debug_)
0078     cout << "1) Parsing displacedVertexCandidates into displacedVertexSeeds" << endl;
0079
0080   // 1) Parsing displacedVertexCandidates into displacedVertexSeeds which would
0081   // be later used for vertex fitting
0082
0083   int i = -1;
0084
0085   for (auto const& idvc : *displacedVertexCandidates_) {
0086     i++;
0087     if (debug_) {
0088       cout << "Analyse Vertex Candidate " << i << endl;
0089     }
0090
0091     findSeedsFromCandidate(idvc, tempDisplacedVertexSeeds);
0092   }
0093
0094   if (debug_)
0095     cout << "2) Merging Vertex Seeds" << endl;
0096
0097   // 2) Some displacedVertexSeeds coming from different displacedVertexCandidates
0098   // may be closed enough to be merged together. bLocked is an array which keeps the
0099   // information about the seeds which are desabled.
0100   vector<bool> bLockedSeeds;
0101   bLockedSeeds.resize(tempDisplacedVertexSeeds.size());
0102   mergeSeeds(tempDisplacedVertexSeeds, bLockedSeeds);
0103
0104   if (debug_)
0105     cout << "3) Fitting Vertices From Seeds" << endl;
0106
0107   // 3) Fit displacedVertices from displacedVertexSeeds
0108   for (unsigned idv = 0; idv < tempDisplacedVertexSeeds.size(); idv++) {
0109     if (!tempDisplacedVertexSeeds[idv].isEmpty() && !bLockedSeeds[idv]) {
0110       PFDisplacedVertex displacedVertex;
0111       bLockedSeeds[idv] = fitVertexFromSeed(tempDisplacedVertexSeeds[idv], displacedVertex);
0112       if (!bLockedSeeds[idv])
0113         tempDisplacedVertices.emplace_back(displacedVertex);
0114     }
0115   }
0116
0117   if (debug_)
0118     cout << "4) Rejecting Bad Vertices and label them" << endl;
0119
0120   // 4) Reject displaced vertices which may be considered as fakes
0121   vector<bool> bLocked;
0122   bLocked.resize(tempDisplacedVertices.size());
0123   selectAndLabelVertices(tempDisplacedVertices, bLocked);
0124
0125   if (debug_)
0126     cout << "5) Fill the Displaced Vertices" << endl;
0127
0128   // 5) Fill the displacedVertex_ which would be transfered to the producer
0129   displacedVertices_->reserve(tempDisplacedVertices.size());
0130
0131   for (unsigned idv = 0; idv < tempDisplacedVertices.size(); idv++)
0132     if (!bLocked[idv])
0133       displacedVertices_->push_back(tempDisplacedVertices[idv]);
0134
0135   if (debug_)
0136     cout << "========= End Find Displaced Vertices =========" << endl;
0137 }
0138
0139 // -------- Different steps of the finder algorithm -------- //
0140
0141 void PFDisplacedVertexFinder::findSeedsFromCandidate(const PFDisplacedVertexCandidate& vertexCandidate,
0142                                                      PFDisplacedVertexSeedCollection& tempDisplacedVertexSeeds) {
0143   const PFDisplacedVertexCandidate::DistMap r2Map = vertexCandidate.r2Map();
0144   bool bNeedNewCandidate = false;
0145
0146   tempDisplacedVertexSeeds.push_back(PFDisplacedVertexSeed());
0147
0148   IDVS idvc_current;
0149
0150   for (PFDisplacedVertexCandidate::DistMap::const_iterator imap = r2Map.begin(); imap != r2Map.end(); imap++) {
0151     unsigned ie1 = (*imap).second.first;
0152     unsigned ie2 = (*imap).second.second;
0153
0154     if (debug_)
0155       cout << "ie1 = " << ie1 << " ie2 = " << ie2 << " radius = " << sqrt((*imap).first) << endl;
0156
0157     GlobalPoint dcaPoint = vertexCandidate.dcaPoint(ie1, ie2);
0158     if (fabs(dcaPoint.x()) > 1e9)
0159       continue;
0160
0161     bNeedNewCandidate = true;
0162     for (idvc_current = tempDisplacedVertexSeeds.begin(); idvc_current != tempDisplacedVertexSeeds.end();
0163          idvc_current++) {
0164       if ((*idvc_current).isEmpty()) {
0165         bNeedNewCandidate = false;
0166         break;
0167       }
0168       const GlobalPoint vertexPoint = (*idvc_current).seedPoint();
0169       std::pair<float, float> diffs = getTransvLongDiff(vertexPoint, dcaPoint);
0170       if (diffs.second > longSize_)
0171         continue;
0172       if (diffs.first > transvSize_)
0173         continue;
0174       bNeedNewCandidate = false;
0175       break;
0176     }
0177     if (bNeedNewCandidate) {
0178       if (debug_)
0179         cout << "create new displaced vertex" << endl;
0180       tempDisplacedVertexSeeds.push_back(PFDisplacedVertexSeed());
0181       idvc_current = tempDisplacedVertexSeeds.end();
0182       idvc_current--;
0183     }
0184
0185     (*idvc_current).updateSeedPoint(dcaPoint, vertexCandidate.tref(ie1), vertexCandidate.tref(ie2));
0186   }
0187 }
0188
0189 void PFDisplacedVertexFinder::mergeSeeds(PFDisplacedVertexSeedCollection& tempDisplacedVertexSeeds,
0190                                          vector<bool>& bLocked) {
0191   // loop over displaced vertex candidates
0192   // and merge them if they are close to each other
0193   for (unsigned idv_mother = 0; idv_mother < tempDisplacedVertexSeeds.size(); idv_mother++) {
0194     if (!bLocked[idv_mother]) {
0195       for (unsigned idv_daughter = idv_mother + 1; idv_daughter < tempDisplacedVertexSeeds.size(); idv_daughter++) {
0196         if (!bLocked[idv_daughter]) {
0197           if (isCloseTo(tempDisplacedVertexSeeds[idv_mother], tempDisplacedVertexSeeds[idv_daughter])) {
0198             tempDisplacedVertexSeeds[idv_mother].mergeWith(tempDisplacedVertexSeeds[idv_daughter]);
0199             bLocked[idv_daughter] = true;
0200             if (debug_)
0201               cout << "Seeds " << idv_mother << " and " << idv_daughter << " merged" << endl;
0202           }
0203         }
0204       }
0205     }
0206   }
0207 }
0208
0209 bool PFDisplacedVertexFinder::fitVertexFromSeed(const PFDisplacedVertexSeed& displacedVertexSeed,
0210                                                 PFDisplacedVertex& displacedVertex) {
0211   if (debug_)
0212     cout << "== Start vertexing procedure ==" << endl;
0213
0214   // ---- Prepare transient track list ----
0215
0216   auto const& tracksToFit = displacedVertexSeed.elements();
0217   const GlobalPoint& seedPoint = displacedVertexSeed.seedPoint();
0218
0219   vector<TransientTrack> transTracks;
0220   vector<TransientTrack> transTracksRaw;
0221   vector<TrackBaseRef> transTracksRef;
0222
0223   transTracks.reserve(tracksToFit.size());
0224   transTracksRaw.reserve(tracksToFit.size());
0225   transTracksRef.reserve(tracksToFit.size());
0226
0227   TransientVertex theVertexAdaptiveRaw;
0228   TransientVertex theRecoVertex;
0229
0230   // ---- 1) Clean for potentially poor seeds ------- //
0231   // --------------------------------------------- //
0232
0233   if (tracksToFit.size() < 2) {
0234     if (debug_)
0235       cout << "Only one to Fit Track" << endl;
0236     return true;
0237   }
0238
0239   double rho = sqrt(seedPoint.x() * seedPoint.x() + seedPoint.y() * seedPoint.y());
0240   double z = seedPoint.z();
0241
0242   if (rho > tobCut_ || fabs(z) > tecCut_) {
0243     if (debug_)
0244       cout << "Seed Point out of the tracker rho = " << rho << " z = " << z << " nTracks = " << tracksToFit.size()
0245            << endl;
0246     return true;
0247   }
0248
0249   if (debug_)
0250     displacedVertexSeed.Dump();
0251
0252   int nStep45 = 0;
0253   int nNotIterative = 0;
0254
0255   // Fill vectors of TransientTracks and TrackRefs after applying preselection cuts.
0256   for (auto const& ie : tracksToFit) {
0257     TransientTrack tmpTk(*(ie.get()), magField_, globTkGeomHandle_);
0258     transTracksRaw.emplace_back(tmpTk);
0259     bool nonIt = PFTrackAlgoTools::nonIterative((ie)->algo());
0260     bool step45 = PFTrackAlgoTools::step45((ie)->algo());
0261     bool highQ = PFTrackAlgoTools::highQuality((ie)->algo());
0262     if (step45)
0263       nStep45++;
0264     else if (nonIt)
0265       nNotIterative++;
0266     else if (!highQ) {
0267       nNotIterative++;
0268       nStep45++;
0269     }
0270   }
0271
0272   if (rho > 25 && nStep45 + nNotIterative < 1) {
0273     if (debug_)
0274       cout << "Seed point at rho > 25 cm but no step 4-5 tracks" << endl;
0275     return true;
0276   }
0277
0278   // ----------------------------------------------- //
0279   // ---- PRESELECT GOOD TRACKS FOR FINAL VERTEX --- //
0280   // ----------------------------------------------- //
0281
0282   // 1)If only two track are found do not prefit
0283
0284   if (transTracksRaw.size() == 2) {
0285     if (debug_)
0286       cout << "No raw fit done" << endl;
0287     if (switchOff2TrackVertex_) {
0288       if (debug_)
0289         cout << "Due to probably high pile-up conditions 2 track vertices switched off" << endl;
0290       return true;
0291     }
0292     GlobalError globalError;
0293
0294     theVertexAdaptiveRaw = TransientVertex(seedPoint, globalError, transTracksRaw, 1.);
0295
0296   } else {  //branch with transTracksRaw.size of at least 3
0297
0298     if (debug_)
0299       cout << "Raw fit done." << endl;
0300
0301     AdaptiveVertexFitter theAdaptiveFitterRaw(GeometricAnnealing(sigmacut_, t_ini_, ratio_),
0302                                               DefaultLinearizationPointFinder(),
0303                                               KalmanVertexUpdator<5>(),
0304                                               KalmanVertexTrackCompatibilityEstimator<5>(),
0305                                               KalmanVertexSmoother());
0306
0307     if (transTracksRaw.size() == 3) {
0308       theVertexAdaptiveRaw = theAdaptiveFitterRaw.vertex(transTracksRaw, seedPoint);
0309
0310     } else if (transTracksRaw.size() < 1000) {
0311       /// This prefit procedure allow to reduce the Warning rate from Adaptive Vertex fitter
0312       /// It reject also many fake tracks
0313
0314       if (debug_)
0315         cout << "First test with KFT" << endl;
0316
0317       KalmanVertexFitter theKalmanFitter(true);
0318       theVertexAdaptiveRaw = theKalmanFitter.vertex(transTracksRaw, seedPoint);
0319
0320       if (!theVertexAdaptiveRaw.isValid() || theVertexAdaptiveRaw.totalChiSquared() < 0.) {
0321         if (debug_)
0322           cout << "Prefit failed : valid? " << theVertexAdaptiveRaw.isValid()
0323                << " totalChi2 = " << theVertexAdaptiveRaw.totalChiSquared() << endl;
0324         return true;
0325       }
0326
0327       if (debug_)
0328         cout << "We use KFT instead of seed point to set up a point for AVF "
0329              << " x = " << theVertexAdaptiveRaw.position().x() << " y = " << theVertexAdaptiveRaw.position().y()
0330              << " z = " << theVertexAdaptiveRaw.position().z() << endl;
0331
0332       // To save time: reject the Displaced vertex if it is too close to the beam pipe.
0333       // Frequently it is very big vertices, with some dosens of tracks
0334
0335       Vertex vtx = theVertexAdaptiveRaw;
0336       rho = vtx.position().rho();
0337
0338       //   cout << "primary vertex cut  = " << primaryVertexCut_ << endl;
0339
0340       if (rho < primaryVertexCut_ || rho > 100) {
0341         if (debug_)
0342           cout << "KFT Vertex geometrically rejected with  tracks #rho = " << rho << endl;
0343         return true;
0344       }
0345
0346       //     cout << "primary vertex cut  = " << primaryVertexCut_ << " rho = " << rho << endl;
0347
0348       theVertexAdaptiveRaw = theAdaptiveFitterRaw.vertex(transTracksRaw, theVertexAdaptiveRaw.position());
0349
0350     } else {
0351       edm::LogWarning("TooManyPFDVCandidates") << "gave up vertex reco for " << transTracksRaw.size() << " tracks";
0352     }
0353
0354     if (!theVertexAdaptiveRaw.isValid() || theVertexAdaptiveRaw.totalChiSquared() < 0.) {
0355       if (debug_)
0356         cout << "Fit failed : valid? " << theVertexAdaptiveRaw.isValid()
0357              << " totalChi2 = " << theVertexAdaptiveRaw.totalChiSquared() << endl;
0358       return true;
0359     }
0360
0361     // To save time: reject the Displaced vertex if it is too close to the beam pipe.
0362     // Frequently it is very big vertices, with some dosens of tracks
0363
0364     Vertex vtx = theVertexAdaptiveRaw;
0365     rho = vtx.position().rho();
0366
0367     if (rho < primaryVertexCut_) {
0368       if (debug_)
0369         cout << "Vertex "
0370              << " geometrically rejected with " << transTracksRaw.size() << " tracks #rho = " << rho << endl;
0371       return true;
0372     }
0373   }
0374
0375   // ---- Remove tracks with small weight or
0376   //      big first (last) hit_to_vertex distance
0377   //      and then refit                          ---- //
0378
0379   for (unsigned i = 0; i < transTracksRaw.size(); i++) {
0380     if (debug_)
0381       cout << "Raw Weight track " << i << " = " << theVertexAdaptiveRaw.trackWeight(transTracksRaw[i]) << endl;
0382
0383     if (theVertexAdaptiveRaw.trackWeight(transTracksRaw[i]) > minAdaptWeight_) {
0384       PFTrackHitFullInfo pattern = hitPattern_.analyze(tkerTopo_, tkerGeom_, tracksToFit[i], theVertexAdaptiveRaw);
0385
0386       PFDisplacedVertex::VertexTrackType vertexTrackType = getVertexTrackType(pattern);
0387
0388       if (vertexTrackType != PFDisplacedVertex::T_NOT_FROM_VERTEX) {
0389         bool bGoodTrack = helper_.isTrackSelected(transTracksRaw[i].track(), vertexTrackType);
0390
0391         if (bGoodTrack) {
0392           transTracks.push_back(transTracksRaw[i]);
0393           transTracksRef.push_back(tracksToFit[i]);
0394         } else {
0395           if (debug_)
0396             cout << "Track rejected nChi2 = " << transTracksRaw[i].track().normalizedChi2()
0397                  << " pt = " << transTracksRaw[i].track().pt()
0398                  << " dxy (wrt (0,0,0)) = " << transTracksRaw[i].track().dxy()
0399                  << " nHits = " << transTracksRaw[i].track().numberOfValidHits() << " nOuterHits = "
0400                  << transTracksRaw[i].track().hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS) << endl;
0401         }
0402       } else {
0403         if (debug_) {
0404           cout << "Remove track because too far away from the vertex:" << endl;
0405         }
0406       }
0407     }
0408   }
0409
0410   if (debug_)
0411     cout << "All Tracks " << transTracksRaw.size() << " with good weight " << transTracks.size() << endl;
0412
0413   // ---- Refit ---- //
0414   FitterType vtxFitter = F_NOTDEFINED;
0415
0416   if (transTracks.size() < 2)
0417     return true;
0418   else if (transTracks.size() == 2) {
0419     if (switchOff2TrackVertex_) {
0420       if (debug_)
0421         cout << "Due to probably high pile-up conditions 2 track vertices switched off" << endl;
0422       return true;
0423     }
0424     vtxFitter = F_KALMAN;
0425   } else if (transTracks.size() > 2 && transTracksRaw.size() > transTracks.size())
0426     vtxFitter = F_ADAPTIVE;
0427   else if (transTracks.size() > 2 && transTracksRaw.size() == transTracks.size())
0428     vtxFitter = F_DONOTREFIT;
0429   else
0430     return true;
0431
0432   if (debug_)
0433     cout << "Vertex Fitter " << vtxFitter << endl;
0434
0435   if (vtxFitter == F_KALMAN) {
0436     KalmanVertexFitter theKalmanFitter(true);
0437     theRecoVertex = theKalmanFitter.vertex(transTracks, seedPoint);
0438
0439   } else if (vtxFitter == F_ADAPTIVE) {
0440     AdaptiveVertexFitter theAdaptiveFitter(GeometricAnnealing(sigmacut_, t_ini_, ratio_),
0441                                            DefaultLinearizationPointFinder(),
0442                                            KalmanVertexUpdator<5>(),
0443                                            KalmanVertexTrackCompatibilityEstimator<5>(),
0444                                            KalmanVertexSmoother());
0445
0446     theRecoVertex = theAdaptiveFitter.vertex(transTracks, seedPoint);
0447
0448   } else if (vtxFitter == F_DONOTREFIT) {
0449     theRecoVertex = theVertexAdaptiveRaw;
0450   } else {
0451     return true;
0452   }
0453
0454   // ---- Check if the fitted vertex is valid ---- //
0455
0456   if (!theRecoVertex.isValid() || theRecoVertex.totalChiSquared() < 0.) {
0457     if (debug_)
0458       cout << "Refit failed : valid? " << theRecoVertex.isValid() << " totalChi2 = " << theRecoVertex.totalChiSquared()
0459            << endl;
0460     return true;
0461   }
0462
0463   // ---- Create vertex ---- //
0464
0465   Vertex theRecoVtx = theRecoVertex;
0466
0467   double chi2 = theRecoVtx.chi2();
0468   double ndf = theRecoVtx.ndof();
0469
0470   if (chi2 > TMath::ChisquareQuantile(0.95, ndf)) {
0471     if (debug_)
0472       cout << "Rejected because of chi2 = " << chi2 << " ndf = " << ndf
0473            << " confid. level: " << TMath::ChisquareQuantile(0.95, ndf) << endl;
0474     return true;
0475   }
0476
0477   // ---- Create and fill vector of refitted TransientTracks ---- //
0478
0479   // -----------------------------------------------//
0480   // ---- Prepare and Fill the Displaced Vertex ----//
0481   // -----------------------------------------------//
0482
0483   displacedVertex = theRecoVtx;
0484   displacedVertex.removeTracks();
0485
0486   for (unsigned i = 0; i < transTracks.size(); i++) {
0487     PFTrackHitFullInfo pattern = hitPattern_.analyze(tkerTopo_, tkerGeom_, transTracksRef[i], theRecoVertex);
0488
0489     PFDisplacedVertex::VertexTrackType vertexTrackType = getVertexTrackType(pattern);
0490
0491     Track refittedTrack;
0492     float weight = theRecoVertex.trackWeight(transTracks[i]);
0493
0494     if (weight < minAdaptWeight_)
0495       continue;
0496
0497     try {
0498       refittedTrack = theRecoVertex.refittedTrack(transTracks[i]).track();
0499     } catch (cms::Exception& exception) {
0500       continue;
0501     }
0502
0503     if (debug_) {
0504       cout << "Vertex Track Type = " << vertexTrackType << endl;
0505
0506       cout << "nHitBeforeVertex = " << pattern.first.first << " nHitAfterVertex = " << pattern.second.first
0507            << " nMissHitBeforeVertex = " << pattern.first.second << " nMissHitAfterVertex = " << pattern.second.second
0508            << " Weight = " << weight << endl;
0509     }
0510
0511     displacedVertex.addElement(transTracksRef[i], refittedTrack, pattern, vertexTrackType, weight);
0512   }
0513
0514   displacedVertex.setPrimaryDirection(helper_.primaryVertex());
0515   displacedVertex.calcKinematics();
0516
0517   if (debug_)
0518     cout << "== End vertexing procedure ==" << endl;
0519
0520   return false;
0521 }
0522
0523 void PFDisplacedVertexFinder::selectAndLabelVertices(PFDisplacedVertexCollection& tempDisplacedVertices,
0524                                                      vector<bool>& bLocked) {
0525   if (debug_)
0526     cout << " 4.1) Reject vertices " << endl;
0527
0528   for (unsigned idv = 0; idv < tempDisplacedVertices.size(); idv++) {
0529     // ---- We accept a vertex only if it is not in TOB in the barrel
0530     //      and not in TEC in the end caps
0531     //      and not too much before the first pixel layer            ---- //
0532
0533     const float rho = tempDisplacedVertices[idv].position().rho();
0534     const float z = tempDisplacedVertices[idv].position().z();
0535
0536     if (rho > tobCut_ || rho < primaryVertexCut_ || fabs(z) > tecCut_) {
0537       if (debug_)
0538         cout << "Vertex " << idv << " geometrically rejected #rho = " << rho << " z = " << z << endl;
0539
0540       bLocked[idv] = true;
0541
0542       continue;
0543     }
0544
0545     unsigned nPrimary = tempDisplacedVertices[idv].nPrimaryTracks();
0546     unsigned nMerged = tempDisplacedVertices[idv].nMergedTracks();
0547     unsigned nSecondary = tempDisplacedVertices[idv].nSecondaryTracks();
0548
0549     if (nPrimary + nMerged > 1) {
0550       bLocked[idv] = true;
0551       if (debug_)
0552         cout << "Vertex " << idv << " rejected because two primary or merged tracks" << endl;
0553     }
0554
0555     if (nPrimary + nMerged + nSecondary < 2) {
0556       bLocked[idv] = true;
0557       if (debug_)
0558         cout << "Vertex " << idv << " rejected because only one track related to the vertex" << endl;
0559     }
0560   }
0561
0562   if (debug_)
0563     cout << " 4.2) Check for common vertices" << endl;
0564
0565   // ---- Among the remaining vertex we shall remove one
0566   //      of those which have two common tracks          ---- //
0567
0568   for (unsigned idv_mother = 0; idv_mother < tempDisplacedVertices.size(); idv_mother++) {
0569     for (unsigned idv_daughter = idv_mother + 1; idv_daughter < tempDisplacedVertices.size(); idv_daughter++) {
0570       if (!bLocked[idv_daughter] && !bLocked[idv_mother]) {
0571         const unsigned commonTrks =
0572             commonTracks(tempDisplacedVertices[idv_daughter], tempDisplacedVertices[idv_mother]);
0573
0574         if (commonTrks > 1) {
0575           if (debug_)
0576             cout << "Vertices " << idv_daughter << " and " << idv_mother << " has many common tracks" << endl;
0577
0578           // We keep the vertex vertex which contains the most of the tracks
0579
0580           const int mother_size = tempDisplacedVertices[idv_mother].nTracks();
0581           const int daughter_size = tempDisplacedVertices[idv_daughter].nTracks();
0582
0583           if (mother_size > daughter_size)
0584             bLocked[idv_daughter] = true;
0585           else if (mother_size < daughter_size)
0586             bLocked[idv_mother] = true;
0587           else {
0588             // If they have the same size we keep the vertex with the best normalised chi2
0589
0590             const float mother_normChi2 = tempDisplacedVertices[idv_mother].normalizedChi2();
0591             const float daughter_normChi2 = tempDisplacedVertices[idv_daughter].normalizedChi2();
0592             if (mother_normChi2 < daughter_normChi2)
0593               bLocked[idv_daughter] = true;
0594             else
0595               bLocked[idv_mother] = true;
0596           }
0597         }
0598       }
0599     }
0600   }
0601
0602   for (unsigned idv = 0; idv < tempDisplacedVertices.size(); idv++)
0603     if (!bLocked[idv])
0604       bLocked[idv] = rejectAndLabelVertex(tempDisplacedVertices[idv]);
0605 }
0606
0607 bool PFDisplacedVertexFinder::rejectAndLabelVertex(PFDisplacedVertex& dv) {
0608   PFDisplacedVertex::VertexType vertexType = helper_.identifyVertex(dv);
0609   dv.setVertexType(vertexType);
0610
0611   return dv.isFake();
0612 }
0613
0614 /// -------- Tools -------- ///
0615
0616 bool PFDisplacedVertexFinder::isCloseTo(const PFDisplacedVertexSeed& dv1, const PFDisplacedVertexSeed& dv2) const {
0617   const GlobalPoint& vP1 = dv1.seedPoint();
0618   const GlobalPoint& vP2 = dv2.seedPoint();
0619
0620   std::pair<float, float> diffs = getTransvLongDiff(vP1, vP2);
0621   if (diffs.second > longSize_)
0622     return false;
0623   if (diffs.first > transvSize_)
0624     return false;
0625   //  if (Delta_Long < longSize_ && Delta_Transv < transvSize_) isCloseTo = true;
0626
0627   return true;
0628 }
0629
0630 std::pair<float, float> PFDisplacedVertexFinder::getTransvLongDiff(const GlobalPoint& Ref,
0631                                                                    const GlobalPoint& ToProject) const {
0632   const auto& vRef = Ref.basicVector();
0633   const auto& vToProject = ToProject.basicVector();
0634   float vRefMag2 = vRef.mag2();
0635   float oneOverMag = 1.0f / sqrt(vRefMag2);
0636
0637   return std::make_pair(fabs(vRef.cross(vToProject).mag() * oneOverMag),
0638                         fabs((vRef.dot(vToProject) - vRefMag2) * oneOverMag));
0639 }
0640
0641 reco::PFDisplacedVertex::VertexTrackType PFDisplacedVertexFinder::getVertexTrackType(
0642     PFTrackHitFullInfo& pairTrackHitInfo) const {
0643   unsigned int nHitBeforeVertex = pairTrackHitInfo.first.first;
0644   unsigned int nHitAfterVertex = pairTrackHitInfo.second.first;
0645
0646   unsigned int nMissHitBeforeVertex = pairTrackHitInfo.first.second;
0647   unsigned int nMissHitAfterVertex = pairTrackHitInfo.second.second;
0648
0649   // For the moment those definitions are given apriori a more general study would be useful to refine those criteria
0650
0651   if (nHitBeforeVertex <= 1 && nHitAfterVertex >= 3 && nMissHitAfterVertex <= 1)
0652     return PFDisplacedVertex::T_FROM_VERTEX;
0653   else if (nHitBeforeVertex >= 3 && nHitAfterVertex <= 1 && nMissHitBeforeVertex <= 1)
0654     return PFDisplacedVertex::T_TO_VERTEX;
0655   else if ((nHitBeforeVertex >= 2 && nHitAfterVertex >= 3) || (nHitBeforeVertex >= 3 && nHitAfterVertex >= 2))
0656     return PFDisplacedVertex::T_MERGED;
0657   else
0658     return PFDisplacedVertex::T_NOT_FROM_VERTEX;
0659 }
0660
0661 unsigned PFDisplacedVertexFinder::commonTracks(const PFDisplacedVertex& v1, const PFDisplacedVertex& v2) const {
0662   vector<Track> vt1 = v1.refittedTracks();
0663   vector<Track> vt2 = v2.refittedTracks();
0664
0665   unsigned commonTracks = 0;
0666
0667   for (unsigned il1 = 0; il1 < vt1.size(); il1++) {
0668     unsigned il1_idx = v1.originalTrack(vt1[il1]).key();
0669
0670     for (unsigned il2 = 0; il2 < vt2.size(); il2++)
0671       if (il1_idx == v2.originalTrack(vt2[il2]).key()) {
0672         commonTracks++;
0673         break;
0674       }
0675   }
0676
0677   return commonTracks;
0678 }
0679
0680 std::ostream& operator<<(std::ostream& out, const PFDisplacedVertexFinder& a) {
0681   if (!out)
0682     return out;
0683   out << setprecision(3) << setw(5) << endl;
0684   out << "" << endl;
0685   out << " ====================================== " << endl;
0686   out << " ====== Displaced Vertex Finder ======= " << endl;
0687   out << " ====================================== " << endl;
0688   out << " " << endl;
0689
0690   a.helper_.Dump();
0691   out << "" << endl
0692       << " Adaptive Vertex Fitter parameters are :" << endl
0693       << " sigmacut = " << a.sigmacut_ << " T_ini = " << a.t_ini_ << " ratio = " << a.ratio_ << endl
0694       << endl;
0695
0696   const std::unique_ptr<reco::PFDisplacedVertexCollection>& displacedVertices_ = std::move(a.displacedVertices());
0697
0698   if (!displacedVertices_.get()) {
0699     out << "displacedVertex already transfered" << endl;
0700   } else {
0701     out << "Number of displacedVertices found : " << displacedVertices_->size() << endl << endl;
0702
0703     int i = -1;
0704
0705     for (PFDisplacedVertexFinder::IDV idv = displacedVertices_->begin(); idv != displacedVertices_->end(); idv++) {
0706       i++;
0707       out << i << " ";
0708       idv->Dump();
0709       out << "" << endl;
0710     }
0711   }
0712
0713   return out;
0714 }