RecoTau/src/RecoTauVertexAssociator.cc

0001 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0002
0003 #include <functional>
0004
0005 #include "DataFormats/TauReco/interface/PFTau.h"
0006 #include "DataFormats/VertexReco/interface/Vertex.h"
0007 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0012 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014
0015 namespace reco {
0016   namespace tau {
0017
0018     namespace {
0019       inline const reco::Track* getTrack(const Candidate& cand) {
0020         const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0021         if (pfCandPtr != nullptr) {
0022           if (pfCandPtr->trackRef().isNonnull())
0023             return pfCandPtr->trackRef().get();
0024           else if (pfCandPtr->gsfTrackRef().isNonnull())
0025             return pfCandPtr->gsfTrackRef().get();
0026           else
0027             return nullptr;
0028         }
0029         const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0030         if (packedCand != nullptr && packedCand->hasTrackDetails())
0031           return &packedCand->pseudoTrack();
0032
0033         return nullptr;
0034       }
0035
0036       inline const reco::TrackBaseRef getTrackRef(const Candidate& cand) {
0037         // TauReco@MiniAOD: This version does not work on top of MiniAOD, however,
0038         // it is only used for non-default track-vertex associations
0039         const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0040         if (pfCandPtr != nullptr) {
0041           if (pfCandPtr->trackRef().isNonnull())
0042             return reco::TrackBaseRef(pfCandPtr->trackRef());
0043           else if (pfCandPtr->gsfTrackRef().isNonnull())
0044             return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
0045           else
0046             return reco::TrackBaseRef();
0047         }
0048
0049         return reco::TrackBaseRef();
0050       }
0051     }  // namespace
0052
0053     // Get the highest pt track in a jet.
0054     // Get the KF track if it exists.  Otherwise, see if it has a GSF track.
0055     const reco::CandidatePtr RecoTauVertexAssociator::getLeadCand(const Jet& jet) const {
0056       std::vector<CandidatePtr> chargedPFCands = pfChargedCands(jet, true);
0057       if (verbosity_ >= 1) {
0058         std::cout << "<RecoTauVertexAssociator::getLeadTrack>:" << std::endl;
0059         std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
0060         std::cout << " num. chargedPFCands = " << chargedPFCands.size() << std::endl;
0061         std::cout << " vxTrkFiltering = " << vxTrkFiltering_ << std::endl;
0062       }
0063
0064       if (chargedPFCands.empty()) {
0065         return reco::CandidatePtr(nullptr, 0);
0066       }
0067
0068       std::vector<CandidatePtr> selectedPFCands;
0069       if (vxTrkFiltering_) {
0070         selectedPFCands = qcuts_->filterCandRefs(chargedPFCands);
0071       } else {
0072         selectedPFCands = chargedPFCands;
0073       }
0074       if (verbosity_ >= 1) {
0075         std::cout << " num. selectedPFCands = " << selectedPFCands.size() << std::endl;
0076       }
0077
0078       CandidatePtr leadCand;
0079       if (!selectedPFCands.empty()) {
0080         double leadTrackPt = 0.;
0081         if (leadingTrkOrPFCandOption_ == kFirstTrack) {
0082           leadCand = selectedPFCands[0];
0083         } else {
0084           for (std::vector<CandidatePtr>::const_iterator pfCand = selectedPFCands.begin();
0085                pfCand != selectedPFCands.end();
0086                ++pfCand) {
0087             const reco::Track* track = getTrack(**pfCand);
0088             double actualTrackPt = 0., actualTrackPtError = 0.;
0089             if (track != nullptr) {
0090               actualTrackPt = track->pt();
0091               actualTrackPtError = track->ptError();
0092             }
0093             double trackPt = 0.;
0094             if (leadingTrkOrPFCandOption_ == kLeadTrack) {
0095               trackPt = actualTrackPt - 2. * actualTrackPtError;
0096             } else if (leadingTrkOrPFCandOption_ == kLeadPFCand) {
0097               trackPt = (*pfCand)->pt();
0098             } else if (leadingTrkOrPFCandOption_ == kMinLeadTrackOrPFCand) {
0099               trackPt = std::min(actualTrackPt, (double)(*pfCand)->pt());
0100             } else
0101               assert(0);
0102             if (trackPt > leadTrackPt) {
0103               leadCand = (*pfCand);
0104               leadTrackPt = trackPt;
0105             }
0106           }
0107         }
0108       }
0109       if (leadCand.isNull()) {
0110         if (recoverLeadingTrk_) {
0111           leadCand = chargedPFCands[0];
0112         } else {
0113           return reco::CandidatePtr(nullptr, 0);
0114         }
0115       }
0116       if (verbosity_ >= 1) {
0117         std::cout << "leadCand: Pt = " << leadCand->pt() << ", eta = " << leadCand->eta()
0118                   << ", phi = " << leadCand->phi() << std::endl;
0119       }
0120       return leadCand;
0121     }
0122
0123     const reco::Track* RecoTauVertexAssociator::getLeadTrack(const Jet& jet) const {
0124       auto leadCand = getLeadCand(jet);
0125       if (leadCand.isNull())
0126         return nullptr;
0127       const reco::Track* track = getTrack(*leadCand);
0128       return track;
0129     }
0130
0131     const reco::TrackBaseRef RecoTauVertexAssociator::getLeadTrackRef(const Jet& jet) const {
0132       auto leadCand = getLeadCand(jet);
0133       if (leadCand.isNull())
0134         return reco::TrackBaseRef();
0135       return getTrackRef(*leadCand);
0136     }
0137
0138     namespace {
0139       // Define functors which extract the relevant information from a collection of
0140       // vertices.
0141       double dzToTrack(const reco::VertexRef& vtx, const reco::Track* trk) {
0142         if (!trk || !vtx) {
0143           return std::numeric_limits<double>::infinity();
0144         }
0145         return std::abs(trk->dz(vtx->position()));
0146       }
0147
0148       double trackWeightInVertex(const reco::VertexRef& vtx, const reco::TrackBaseRef& trk) {
0149         if (!trk || !vtx) {
0150           return 0.0;
0151         }
0152         return vtx->trackWeight(trk);
0153       }
0154     }  // namespace
0155
0156     RecoTauVertexAssociator::RecoTauVertexAssociator(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0157         : vertexSelector_(nullptr), qcuts_(nullptr), jetToVertexAssociation_(nullptr), lastEvent_(0) {
0158       //std::cout << "<RecoTauVertexAssociator::RecoTauVertexAssociator>:" << std::endl;
0159
0160       vertexTag_ = edm::InputTag("offlinePrimaryVertices", "");
0161       algorithm_ = "highestPtInEvent";
0162       // Sanity check, will remove once HLT module configs are updated.
0163       if (!pset.exists("primaryVertexSrc") || !pset.exists("pvFindingAlgo")) {
0164         edm::LogWarning("NoVertexFindingMethodSpecified")
0165             << "The PSet passed to the RecoTauVertexAssociator was incorrectly configured."
0166             << " The vertex will be taken as the highest Pt vertex from the 'offlinePrimaryVertices' collection."
0167             << std::endl;
0168       } else {
0169         vertexTag_ = pset.getParameter<edm::InputTag>("primaryVertexSrc");
0170         algorithm_ = pset.getParameter<std::string>("pvFindingAlgo");
0171       }
0172
0173       if (pset.exists("vxAssocQualityCuts")) {
0174         //std::cout << " reading 'vxAssocQualityCuts'" << std::endl;
0175         qcuts_ = std::make_unique<RecoTauQualityCuts>(pset.getParameterSet("vxAssocQualityCuts"));
0176       } else {
0177         //std::cout << " reading 'signalQualityCuts'" << std::endl;
0178         qcuts_ = std::make_unique<RecoTauQualityCuts>(pset.getParameterSet("signalQualityCuts"));
0179       }
0180       assert(qcuts_);
0181
0182       vxTrkFiltering_ = false;
0183       if (!pset.exists("vertexTrackFiltering") && pset.exists("vxAssocQualityCuts")) {
0184         edm::LogWarning("NoVertexTrackFilteringSpecified")
0185             << "The PSet passed to the RecoTauVertexAssociator was incorrectly configured."
0186             << " Please define vertexTrackFiltering in config file."
0187             << " No filtering of tracks to vertices will be applied." << std::endl;
0188       } else {
0189         vxTrkFiltering_ = pset.exists("vertexTrackFiltering") ? pset.getParameter<bool>("vertexTrackFiltering") : false;
0190       }
0191       if (pset.exists("vertexSelection")) {
0192         std::string vertexSelection = pset.getParameter<std::string>("vertexSelection");
0193         if (!vertexSelection.empty()) {
0194           vertexSelector_ = std::make_unique<StringCutObjectSelector<reco::Vertex>>(vertexSelection);
0195         }
0196       }
0197
0198       if (algorithm_ == "highestPtInEvent") {
0199         algo_ = kHighestPtInEvent;
0200       } else if (algorithm_ == "closestInDeltaZ") {
0201         algo_ = kClosestDeltaZ;
0202       } else if (algorithm_ == "highestWeightForLeadTrack") {
0203         algo_ = kHighestWeigtForLeadTrack;
0204       } else if (algorithm_ == "combined") {
0205         algo_ = kCombined;
0206       } else {
0207         throw cms::Exception("BadVertexAssociatorConfig")
0208             << "Invalid Configuration parameter 'algorithm' " << algorithm_ << "."
0209             << " Valid options are: 'highestPtInEvent', 'closestInDeltaZ', 'highestWeightForLeadTrack' and "
0210                "'combined'.\n";
0211       }
0212
0213       vxToken_ = iC.consumes<reco::VertexCollection>(vertexTag_);
0214
0215       recoverLeadingTrk_ = pset.exists("recoverLeadingTrk") ? pset.getParameter<bool>("recoverLeadingTrk") : false;
0216
0217       std::string leadingTrkOrPFCandOption_string = pset.exists("leadingTrkOrPFCandOption")
0218                                                         ? pset.getParameter<std::string>("leadingTrkOrPFCandOption")
0219                                                         : "firstTrack";
0220       if (leadingTrkOrPFCandOption_string == "leadTrack")
0221         leadingTrkOrPFCandOption_ = kLeadTrack;
0222       else if (leadingTrkOrPFCandOption_string == "leadPFCand")
0223         leadingTrkOrPFCandOption_ = kLeadPFCand;
0224       else if (leadingTrkOrPFCandOption_string == "minLeadTrackOrPFCand")
0225         leadingTrkOrPFCandOption_ = kMinLeadTrackOrPFCand;
0226       else if (leadingTrkOrPFCandOption_string == "firstTrack")
0227         leadingTrkOrPFCandOption_ = kFirstTrack;
0228       else
0229         throw cms::Exception("BadVertexAssociatorConfig")
0230             << "Invalid Configuration parameter 'leadingTrkOrPFCandOption' " << leadingTrkOrPFCandOption_string << "."
0231             << " Valid options are: 'leadTrack', 'leadPFCand', 'firstTrack'.\n";
0232
0233       verbosity_ = (pset.exists("verbosity")) ? pset.getParameter<int>("verbosity") : 0;
0234     }
0235
0236     void RecoTauVertexAssociator::setEvent(const edm::Event& evt) {
0237       edm::Handle<reco::VertexCollection> vertices;
0238       evt.getByToken(vxToken_, vertices);
0239       selectedVertices_.clear();
0240       selectedVertices_.reserve(vertices->size());
0241       for (size_t idxVertex = 0; idxVertex < vertices->size(); ++idxVertex) {
0242         reco::VertexRef vertex(vertices, idxVertex);
0243         if (vertexSelector_ && !(*vertexSelector_)(*vertex))
0244           continue;
0245         selectedVertices_.push_back(vertex);
0246       }
0247       if (!selectedVertices_.empty()) {
0248         qcuts_->setPV(selectedVertices_[0]);
0249       }
0250       edm::EventNumber_t currentEvent = evt.id().event();
0251       if (currentEvent != lastEvent_ || !jetToVertexAssociation_) {
0252         if (!jetToVertexAssociation_)
0253           jetToVertexAssociation_ = std::make_unique<JetToVtxAssoc>();
0254         else
0255           jetToVertexAssociation_->clear();
0256         lastEvent_ = currentEvent;
0257       }
0258     }
0259
0260     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const PFTau& tau, bool useJet) const {
0261       if (!useJet) {
0262         if (tau.leadChargedHadrCand().isNonnull()) {
0263           const reco::Track* track = getTrack(*tau.leadChargedHadrCand());
0264           if (track != nullptr)
0265             return associatedVertex(track);
0266         }
0267       }
0268       // MB: use vertex associated to a given jet if explicitely requested or in case of missing leading track
0269       reco::JetBaseRef jetRef = tau.jetRef();
0270       // FIXME workaround for HLT which does not use updated data format
0271       if (jetRef.isNull())
0272         jetRef = tau.pfTauTagInfoRef()->pfjetRef();
0273       return associatedVertex(*jetRef);
0274     }
0275
0276     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const Track* track) const {
0277       reco::VertexRef trkVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0278
0279       // algos kHighestWeigtForLeadTrack and kCombined not supported
0280       if (algo_ == kHighestPtInEvent) {
0281         if (!selectedVertices_.empty())
0282           trkVertex = selectedVertices_[0];
0283       } else if (algo_ == kClosestDeltaZ) {
0284         if (track) {
0285           double closestDistance = 1.e+6;
0286           // Find the vertex that has the lowest dZ to the track
0287           int idxVertex = 0;
0288           for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0289                selectedVertex != selectedVertices_.end();
0290                ++selectedVertex) {
0291             double dZ = dzToTrack(*selectedVertex, track);
0292             if (verbosity_) {
0293               std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0294                         << ", y = " << (*selectedVertex)->position().y()
0295                         << ", z = " << (*selectedVertex)->position().z() << " --> dZ = " << dZ << std::endl;
0296             }
0297             if (dZ < closestDistance) {
0298               trkVertex = (*selectedVertex);
0299               closestDistance = dZ;
0300             }
0301             ++idxVertex;
0302           }
0303         }
0304       }
0305
0306       if (verbosity_ >= 1) {
0307         std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y()
0308                   << ", z = " << trkVertex->position().z() << std::endl;
0309       }
0310
0311       return trkVertex;
0312     }
0313
0314     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const TrackBaseRef& track) const {
0315       reco::VertexRef trkVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0316
0317       if (algo_ == kHighestWeigtForLeadTrack || algo_ == kCombined) {
0318         if (track.isNonnull()) {
0319           double largestWeight = -1.;
0320           // Find the vertex that has the highest association probability to the track
0321           int idxVertex = 0;
0322           for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0323                selectedVertex != selectedVertices_.end();
0324                ++selectedVertex) {
0325             double weight = trackWeightInVertex(*selectedVertex, track);
0326             if (verbosity_) {
0327               std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0328                         << ", y = " << (*selectedVertex)->position().y()
0329                         << ", z = " << (*selectedVertex)->position().z() << " --> weight = " << weight << std::endl;
0330             }
0331             if (weight > largestWeight) {
0332               trkVertex = (*selectedVertex);
0333               largestWeight = weight;
0334             }
0335             ++idxVertex;
0336           }
0337           // the weight was never larger than zero
0338           if (algo_ == kCombined && largestWeight < 1.e-7) {
0339             if (verbosity_) {
0340               std::cout << "No vertex had positive weight! Trying dZ instead... " << std::endl;
0341             }
0342             double closestDistance = 1.e+6;
0343             // Find the vertex that has the lowest dZ to the leading track
0344             int idxVertex = 0;
0345             for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0346                  selectedVertex != selectedVertices_.end();
0347                  ++selectedVertex) {
0348               double dZ = dzToTrack(*selectedVertex, track.get());
0349               if (verbosity_) {
0350                 std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0351                           << ", y = " << (*selectedVertex)->position().y()
0352                           << ", z = " << (*selectedVertex)->position().z() << " --> dZ = " << dZ << std::endl;
0353               }
0354               if (dZ < closestDistance) {
0355                 trkVertex = (*selectedVertex);
0356                 closestDistance = dZ;
0357               }
0358               ++idxVertex;
0359             }
0360           }
0361         }
0362       }
0363
0364       if (verbosity_ >= 1) {
0365         std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y()
0366                   << ", z = " << trkVertex->position().z() << std::endl;
0367       }
0368
0369       return trkVertex;
0370     }
0371
0372     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const Jet& jet) const {
0373       if (verbosity_ >= 1) {
0374         std::cout << "<RecoTauVertexAssociator::associatedVertex>:" << std::endl;
0375         std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
0376         std::cout << " num. Vertices = " << selectedVertices_.size() << std::endl;
0377         std::cout << " size(jetToVertexAssociation) = " << jetToVertexAssociation_->size() << std::endl;
0378         std::cout << " vertexTag = " << vertexTag_ << std::endl;
0379         std::cout << " algorithm = " << algorithm_ << std::endl;
0380         std::cout << " recoverLeadingTrk = " << recoverLeadingTrk_ << std::endl;
0381       }
0382
0383       reco::VertexRef jetVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0384       const Jet* jetPtr = &jet;
0385
0386       // check if jet-vertex association has been determined for this jet before
0387       auto vertexPtr = jetToVertexAssociation_->find(jetPtr);
0388       if (vertexPtr != jetToVertexAssociation_->end()) {
0389         jetVertex = vertexPtr->second;
0390       } else {
0391         // no jet-vertex association exists for this jet yet, compute it!
0392         if (algo_ == kHighestPtInEvent) {
0393           if (!selectedVertices_.empty())
0394             jetVertex = selectedVertices_[0];
0395         } else if (algo_ == kClosestDeltaZ) {
0396           // find "leading" (highest Pt) track in jet
0397           const reco::Track* leadTrack = getLeadTrack(jet);
0398           if (verbosity_) {
0399             if (leadTrack != nullptr)
0400               std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta()
0401                         << ", phi = " << leadTrack->phi() << std::endl;
0402             else
0403               std::cout << "leadTrack: N/A" << std::endl;
0404           }
0405           if (leadTrack != nullptr) {
0406             jetVertex = associatedVertex(leadTrack);
0407           }
0408         } else if (algo_ == kHighestWeigtForLeadTrack || algo_ == kCombined) {
0409           // find "leading" (highest Pt) track in jet
0410           const reco::TrackBaseRef leadTrack = getLeadTrackRef(jet);
0411           if (verbosity_) {
0412             if (leadTrack.isNonnull())
0413               std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta()
0414                         << ", phi = " << leadTrack->phi() << std::endl;
0415             else
0416               std::cout << "leadTrack: N/A" << std::endl;
0417           }
0418           if (leadTrack.isNonnull()) {
0419             jetVertex = associatedVertex(leadTrack);
0420           }
0421         }
0422
0423         jetToVertexAssociation_->insert({jetPtr, jetVertex});
0424       }
0425
0426       if (verbosity_ >= 1) {
0427         std::cout << "--> returning vertex: x = " << jetVertex->position().x() << ", y = " << jetVertex->position().y()
0428                   << ", z = " << jetVertex->position().z() << std::endl;
0429       }
0430
0431       return jetVertex;
0432     }
0433
0434   }  // namespace tau
0435 }  // namespace reco