Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:06

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_ = new RecoTauQualityCuts(pset.getParameterSet("vxAssocQualityCuts"));
0176       } else {
0177         //std::cout << " reading 'signalQualityCuts'" << std::endl;
0178         qcuts_ = new 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_ = new 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     RecoTauVertexAssociator::~RecoTauVertexAssociator() {
0237       delete vertexSelector_;
0238       delete qcuts_;
0239       delete jetToVertexAssociation_;
0240     }
0241 
0242     void RecoTauVertexAssociator::setEvent(const edm::Event& evt) {
0243       edm::Handle<reco::VertexCollection> vertices;
0244       evt.getByToken(vxToken_, vertices);
0245       selectedVertices_.clear();
0246       selectedVertices_.reserve(vertices->size());
0247       for (size_t idxVertex = 0; idxVertex < vertices->size(); ++idxVertex) {
0248         reco::VertexRef vertex(vertices, idxVertex);
0249         if (vertexSelector_ && !(*vertexSelector_)(*vertex))
0250           continue;
0251         selectedVertices_.push_back(vertex);
0252       }
0253       if (!selectedVertices_.empty()) {
0254         qcuts_->setPV(selectedVertices_[0]);
0255       }
0256       edm::EventNumber_t currentEvent = evt.id().event();
0257       if (currentEvent != lastEvent_ || !jetToVertexAssociation_) {
0258         if (!jetToVertexAssociation_)
0259           jetToVertexAssociation_ = new std::map<const reco::Jet*, reco::VertexRef>;
0260         else
0261           jetToVertexAssociation_->clear();
0262         lastEvent_ = currentEvent;
0263       }
0264     }
0265 
0266     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const PFTau& tau, bool useJet) const {
0267       if (!useJet) {
0268         if (tau.leadChargedHadrCand().isNonnull()) {
0269           const reco::Track* track = getTrack(*tau.leadChargedHadrCand());
0270           if (track != nullptr)
0271             return associatedVertex(track);
0272         }
0273       }
0274       // MB: use vertex associated to a given jet if explicitely requested or in case of missing leading track
0275       reco::JetBaseRef jetRef = tau.jetRef();
0276       // FIXME workaround for HLT which does not use updated data format
0277       if (jetRef.isNull())
0278         jetRef = tau.pfTauTagInfoRef()->pfjetRef();
0279       return associatedVertex(*jetRef);
0280     }
0281 
0282     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const Track* track) const {
0283       reco::VertexRef trkVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0284 
0285       // algos kHighestWeigtForLeadTrack and kCombined not supported
0286       if (algo_ == kHighestPtInEvent) {
0287         if (!selectedVertices_.empty())
0288           trkVertex = selectedVertices_[0];
0289       } else if (algo_ == kClosestDeltaZ) {
0290         if (track) {
0291           double closestDistance = 1.e+6;
0292           // Find the vertex that has the lowest dZ to the track
0293           int idxVertex = 0;
0294           for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0295                selectedVertex != selectedVertices_.end();
0296                ++selectedVertex) {
0297             double dZ = dzToTrack(*selectedVertex, track);
0298             if (verbosity_) {
0299               std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0300                         << ", y = " << (*selectedVertex)->position().y()
0301                         << ", z = " << (*selectedVertex)->position().z() << " --> dZ = " << dZ << std::endl;
0302             }
0303             if (dZ < closestDistance) {
0304               trkVertex = (*selectedVertex);
0305               closestDistance = dZ;
0306             }
0307             ++idxVertex;
0308           }
0309         }
0310       }
0311 
0312       if (verbosity_ >= 1) {
0313         std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y()
0314                   << ", z = " << trkVertex->position().z() << std::endl;
0315       }
0316 
0317       return trkVertex;
0318     }
0319 
0320     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const TrackBaseRef& track) const {
0321       reco::VertexRef trkVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0322 
0323       if (algo_ == kHighestWeigtForLeadTrack || algo_ == kCombined) {
0324         if (track.isNonnull()) {
0325           double largestWeight = -1.;
0326           // Find the vertex that has the highest association probability to the track
0327           int idxVertex = 0;
0328           for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0329                selectedVertex != selectedVertices_.end();
0330                ++selectedVertex) {
0331             double weight = trackWeightInVertex(*selectedVertex, track);
0332             if (verbosity_) {
0333               std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0334                         << ", y = " << (*selectedVertex)->position().y()
0335                         << ", z = " << (*selectedVertex)->position().z() << " --> weight = " << weight << std::endl;
0336             }
0337             if (weight > largestWeight) {
0338               trkVertex = (*selectedVertex);
0339               largestWeight = weight;
0340             }
0341             ++idxVertex;
0342           }
0343           // the weight was never larger than zero
0344           if (algo_ == kCombined && largestWeight < 1.e-7) {
0345             if (verbosity_) {
0346               std::cout << "No vertex had positive weight! Trying dZ instead... " << std::endl;
0347             }
0348             double closestDistance = 1.e+6;
0349             // Find the vertex that has the lowest dZ to the leading track
0350             int idxVertex = 0;
0351             for (std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
0352                  selectedVertex != selectedVertices_.end();
0353                  ++selectedVertex) {
0354               double dZ = dzToTrack(*selectedVertex, track.get());
0355               if (verbosity_) {
0356                 std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x()
0357                           << ", y = " << (*selectedVertex)->position().y()
0358                           << ", z = " << (*selectedVertex)->position().z() << " --> dZ = " << dZ << std::endl;
0359               }
0360               if (dZ < closestDistance) {
0361                 trkVertex = (*selectedVertex);
0362                 closestDistance = dZ;
0363               }
0364               ++idxVertex;
0365             }
0366           }
0367         }
0368       }
0369 
0370       if (verbosity_ >= 1) {
0371         std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y()
0372                   << ", z = " << trkVertex->position().z() << std::endl;
0373       }
0374 
0375       return trkVertex;
0376     }
0377 
0378     reco::VertexRef RecoTauVertexAssociator::associatedVertex(const Jet& jet) const {
0379       if (verbosity_ >= 1) {
0380         std::cout << "<RecoTauVertexAssociator::associatedVertex>:" << std::endl;
0381         std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
0382         std::cout << " num. Vertices = " << selectedVertices_.size() << std::endl;
0383         std::cout << " size(jetToVertexAssociation) = " << jetToVertexAssociation_->size() << std::endl;
0384         std::cout << " vertexTag = " << vertexTag_ << std::endl;
0385         std::cout << " algorithm = " << algorithm_ << std::endl;
0386         std::cout << " recoverLeadingTrk = " << recoverLeadingTrk_ << std::endl;
0387       }
0388 
0389       reco::VertexRef jetVertex = (!selectedVertices_.empty()) ? selectedVertices_[0] : reco::VertexRef();
0390       const Jet* jetPtr = &jet;
0391 
0392       // check if jet-vertex association has been determined for this jet before
0393       std::map<const reco::Jet*, reco::VertexRef>::iterator vertexPtr = jetToVertexAssociation_->find(jetPtr);
0394       if (vertexPtr != jetToVertexAssociation_->end()) {
0395         jetVertex = vertexPtr->second;
0396       } else {
0397         // no jet-vertex association exists for this jet yet, compute it!
0398         if (algo_ == kHighestPtInEvent) {
0399           if (!selectedVertices_.empty())
0400             jetVertex = selectedVertices_[0];
0401         } else if (algo_ == kClosestDeltaZ) {
0402           // find "leading" (highest Pt) track in jet
0403           const reco::Track* leadTrack = getLeadTrack(jet);
0404           if (verbosity_) {
0405             if (leadTrack != nullptr)
0406               std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta()
0407                         << ", phi = " << leadTrack->phi() << std::endl;
0408             else
0409               std::cout << "leadTrack: N/A" << std::endl;
0410           }
0411           if (leadTrack != nullptr) {
0412             jetVertex = associatedVertex(leadTrack);
0413           }
0414         } else if (algo_ == kHighestWeigtForLeadTrack || algo_ == kCombined) {
0415           // find "leading" (highest Pt) track in jet
0416           const reco::TrackBaseRef leadTrack = getLeadTrackRef(jet);
0417           if (verbosity_) {
0418             if (leadTrack.isNonnull())
0419               std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta()
0420                         << ", phi = " << leadTrack->phi() << std::endl;
0421             else
0422               std::cout << "leadTrack: N/A" << std::endl;
0423           }
0424           if (leadTrack.isNonnull()) {
0425             jetVertex = associatedVertex(leadTrack);
0426           }
0427         }
0428 
0429         jetToVertexAssociation_->insert(std::pair<const Jet*, reco::VertexRef>(jetPtr, jetVertex));
0430       }
0431 
0432       if (verbosity_ >= 1) {
0433         std::cout << "--> returning vertex: x = " << jetVertex->position().x() << ", y = " << jetVertex->position().y()
0434                   << ", z = " << jetVertex->position().z() << std::endl;
0435       }
0436 
0437       return jetVertex;
0438     }
0439 
0440   }  // namespace tau
0441 }  // namespace reco