RecoTau/src/PFTauPrimaryVertexProducerBase.cc

0001 #include "RecoTauTag/RecoTau/interface/PFTauPrimaryVertexProducerBase.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005
0006 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0007 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0008
0009 #include "DataFormats/Common/interface/AssociationVector.h"
0010 #include "DataFormats/Common/interface/RefProd.h"
0011 #include "DataFormats/Common/interface/RefToPtr.h"
0012
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014
0015 #include <memory>
0016
0017 PFTauPrimaryVertexProducerBase::PFTauPrimaryVertexProducerBase(const edm::ParameterSet& iConfig)
0018     : pftauToken_(consumes<std::vector<reco::PFTau>>(iConfig.getParameter<edm::InputTag>("PFTauTag"))),
0019       electronToken_(consumes<edm::View<reco::Electron>>(iConfig.getParameter<edm::InputTag>("ElectronTag"))),
0020       muonToken_(consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("MuonTag"))),
0021       pvToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PVTag"))),
0022       beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))),
0023       transTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0024       algorithm_(iConfig.getParameter<int>("Algorithm")),
0025       qualityCutsPSet_(iConfig.getParameter<edm::ParameterSet>("qualityCuts")),
0026       useBeamSpot_(iConfig.getParameter<bool>("useBeamSpot")),
0027       useSelectedTaus_(iConfig.getParameter<bool>("useSelectedTaus")),
0028       removeMuonTracks_(iConfig.getParameter<bool>("RemoveMuonTracks")),
0029       removeElectronTracks_(iConfig.getParameter<bool>("RemoveElectronTracks")) {
0030   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0031   std::vector<edm::ParameterSet> discriminators =
0032       iConfig.getParameter<std::vector<edm::ParameterSet>>("discriminators");
0033   // Build each of our cuts
0034   for (auto const& pset : discriminators) {
0035     DiscCutPair* newCut = new DiscCutPair();
0036     newCut->inputToken_ = consumes<reco::PFTauDiscriminator>(pset.getParameter<edm::InputTag>("discriminator"));
0037
0038     if (pset.existsAs<std::string>("selectionCut"))
0039       newCut->cutFormula_ = new TFormula("selectionCut", pset.getParameter<std::string>("selectionCut").data());
0040     else
0041       newCut->cut_ = pset.getParameter<double>("selectionCut");
0042     discriminators_.push_back(newCut);
0043   }
0044   // Build a string cut if desired
0045   if (iConfig.exists("cut"))
0046     cut_ = std::make_unique<StringCutObjectSelector<reco::PFTau>>(iConfig.getParameter<std::string>("cut"));
0047   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0048   produces<edm::AssociationVector<reco::PFTauRefProd, std::vector<reco::VertexRef>>>();
0049   produces<reco::VertexCollection>("PFTauPrimaryVertices");
0050
0051   vertexAssociator_ = std::make_unique<reco::tau::RecoTauVertexAssociator>(qualityCutsPSet_, consumesCollector());
0052 }
0053
0054 PFTauPrimaryVertexProducerBase::~PFTauPrimaryVertexProducerBase() {}
0055
0056 namespace {
0057   edm::Ptr<reco::TrackBase> getTrack(const reco::Candidate& cand) {
0058     const reco::PFCandidate* pfCandPtr = dynamic_cast<const reco::PFCandidate*>(&cand);
0059     if (pfCandPtr) {
0060       if (pfCandPtr->trackRef().isNonnull())
0061         return edm::refToPtr(pfCandPtr->trackRef());
0062       else if (pfCandPtr->gsfTrackRef().isNonnull())
0063         return edm::refToPtr(pfCandPtr->gsfTrackRef());
0064       else
0065         return edm::Ptr<reco::TrackBase>();
0066     }
0067     const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0068     if (pCand && pCand->hasTrackDetails()) {
0069       const reco::TrackBase* trkPtr = pCand->bestTrack();
0070       return edm::Ptr<reco::TrackBase>(trkPtr, 0);
0071     }
0072     return edm::Ptr<reco::TrackBase>();
0073   }
0074 }  // namespace
0075
0076 void PFTauPrimaryVertexProducerBase::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0077   beginEvent(iEvent, iSetup);
0078
0079   // Obtain Collections
0080   auto const& transTrackBuilder = iSetup.getData(transTrackBuilderToken_);
0081
0082   edm::Handle<std::vector<reco::PFTau>> pfTaus;
0083   iEvent.getByToken(pftauToken_, pfTaus);
0084
0085   edm::Handle<edm::View<reco::Electron>> electrons;
0086   if (removeElectronTracks_)
0087     iEvent.getByToken(electronToken_, electrons);
0088
0089   edm::Handle<edm::View<reco::Muon>> muons;
0090   if (removeMuonTracks_)
0091     iEvent.getByToken(muonToken_, muons);
0092
0093   edm::Handle<reco::VertexCollection> vertices;
0094   iEvent.getByToken(pvToken_, vertices);
0095
0096   edm::Handle<reco::BeamSpot> beamSpot;
0097   if (useBeamSpot_)
0098     iEvent.getByToken(beamSpotToken_, beamSpot);
0099
0100   // Set Association Map
0101   auto avPFTauPV = std::make_unique<edm::AssociationVector<reco::PFTauRefProd, std::vector<reco::VertexRef>>>(
0102       reco::PFTauRefProd(pfTaus));
0103   auto vertexCollection_out = std::make_unique<reco::VertexCollection>();
0104   reco::VertexRefProd vertexRefProd_out = iEvent.getRefBeforePut<reco::VertexCollection>("PFTauPrimaryVertices");
0105
0106   // Load each discriminator
0107   for (auto& disc : discriminators_) {
0108     edm::Handle<reco::PFTauDiscriminator> discr;
0109     iEvent.getByToken(disc->inputToken_, discr);
0110     disc->discr_ = &(*discr);
0111   }
0112
0113   // Set event for VerexAssociator if needed
0114   if (useInputPV == algorithm_)
0115     vertexAssociator_->setEvent(iEvent);
0116
0117   // For each Tau Run Algorithim
0118   if (pfTaus.isValid()) {
0119     for (reco::PFTauCollection::size_type iPFTau = 0; iPFTau < pfTaus->size(); iPFTau++) {
0120       reco::PFTauRef tau(pfTaus, iPFTau);
0121       reco::VertexRef thePVRef;
0122       if (useInputPV == algorithm_) {
0123         thePVRef = vertexAssociator_->associatedVertex(*tau);
0124       } else if (useFrontPV == algorithm_) {
0125         thePVRef = reco::VertexRef(vertices, 0);
0126       }
0127       reco::Vertex thePV = *thePVRef;
0128       ///////////////////////
0129       // Check if it passed all the discrimiantors
0130       bool passed(true);
0131       for (auto const& disc : discriminators_) {
0132         // Check this discriminator passes
0133         bool passedDisc = true;
0134         if (disc->cutFormula_)
0135           passedDisc = (disc->cutFormula_->Eval((*disc->discr_)[tau]) > 0.5);
0136         else
0137           passedDisc = ((*disc->discr_)[tau] > disc->cut_);
0138         if (!passedDisc) {
0139           passed = false;
0140           break;
0141         }
0142       }
0143       if (passed && cut_.get()) {
0144         passed = (*cut_)(*tau);
0145       }
0146       if (passed) {
0147         std::vector<edm::Ptr<reco::TrackBase>> signalTracks;
0148         for (reco::PFTauCollection::size_type jPFTau = 0; jPFTau < pfTaus->size(); jPFTau++) {
0149           if (useSelectedTaus_ || iPFTau == jPFTau) {
0150             reco::PFTauRef pfTauRef(pfTaus, jPFTau);
0151             ///////////////////////////////////////////////////////////////////////////////////////////////
0152             // Get tracks from PFTau daughters
0153             for (const auto& pfcand : pfTauRef->signalChargedHadrCands()) {
0154               if (pfcand.isNull())
0155                 continue;
0156               const edm::Ptr<reco::TrackBase>& trackPtr = getTrack(*pfcand);
0157               if (trackPtr.isNonnull())
0158                 signalTracks.push_back(trackPtr);
0159             }
0160           }
0161         }
0162         // Get Muon tracks
0163         if (removeMuonTracks_) {
0164           if (muons.isValid()) {
0165             for (const auto& muon : *muons) {
0166               if (muon.track().isNonnull())
0167                 signalTracks.push_back(edm::refToPtr(muon.track()));
0168             }
0169           }
0170         }
0171         // Get Electron Tracks
0172         if (removeElectronTracks_) {
0173           if (electrons.isValid()) {
0174             for (const auto& electron : *electrons) {
0175               if (electron.track().isNonnull())
0176                 signalTracks.push_back(edm::refToPtr(electron.track()));
0177               if (electron.gsfTrack().isNonnull())
0178                 signalTracks.push_back(edm::refToPtr(electron.gsfTrack()));
0179             }
0180           }
0181         }
0182         ///////////////////////////////////////////////////////////////////////////////////////////////
0183         // Get Non-Tau tracks
0184         std::vector<const reco::Track*> nonTauTracks;
0185         nonTauTracksInPV(thePVRef, signalTracks, nonTauTracks);
0186
0187         ///////////////////////////////////////////////////////////////////////////////////////////////
0188         // Refit the vertex
0189         TransientVertex transVtx;
0190         std::vector<reco::TransientTrack> transTracks;
0191         transTracks.reserve(nonTauTracks.size());
0192         for (const auto track : nonTauTracks) {
0193           transTracks.push_back(transTrackBuilder.build(*track));
0194         }
0195         bool fitOK(true);
0196         if (transTracks.size() >= 2) {
0197           AdaptiveVertexFitter avf;
0198           avf.setWeightThreshold(0.1);  //weight per track. allow almost every fit, else --> exception
0199           if (!useBeamSpot_) {
0200             transVtx = avf.vertex(transTracks);
0201           } else {
0202             transVtx = avf.vertex(transTracks, *beamSpot);
0203           }
0204           if (!transVtx.isValid()) {
0205             fitOK = false;
0206           } else {
0207             //MB: protect against rare cases when transVtx is valid but its position is ill-defined
0208             if (!std::isfinite(transVtx.position().z()))  //MB: it is enough to check one coordinate (?)
0209               fitOK = false;
0210           }
0211         } else
0212           fitOK = false;
0213         if (fitOK)
0214           thePV = transVtx;
0215       }
0216       reco::VertexRef vtxRef = reco::VertexRef(vertexRefProd_out, vertexCollection_out->size());
0217       vertexCollection_out->push_back(thePV);
0218       avPFTauPV->setValue(iPFTau, vtxRef);
0219     }
0220   }
0221   iEvent.put(std::move(vertexCollection_out), "PFTauPrimaryVertices");
0222   iEvent.put(std::move(avPFTauPV));
0223 }
0224
0225 edm::ParameterSetDescription PFTauPrimaryVertexProducerBase::getDescriptionsBase() {
0226   // PFTauPrimaryVertexProducerBase
0227   edm::ParameterSetDescription desc;
0228
0229   {
0230     edm::ParameterSetDescription vpsd1;
0231     vpsd1.add<edm::InputTag>("discriminator");
0232     vpsd1.add<double>("selectionCut");
0233     desc.addVPSet("discriminators", vpsd1);
0234   }
0235
0236   edm::ParameterSetDescription desc_qualityCuts;
0237   reco::tau::RecoTauQualityCuts::fillDescriptions(desc_qualityCuts);
0238   desc.add<edm::ParameterSetDescription>("qualityCuts", desc_qualityCuts);
0239
0240   desc.add<std::string>("cut", "pt > 18.0 & abs(eta)<2.3");
0241   desc.add<int>("Algorithm", 0);
0242   desc.add<bool>("RemoveElectronTracks", false);
0243   desc.add<bool>("RemoveMuonTracks", false);
0244   desc.add<bool>("useBeamSpot", true);
0245   desc.add<bool>("useSelectedTaus", false);
0246   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0247   desc.add<edm::InputTag>("ElectronTag", edm::InputTag("MyElectrons"));
0248   desc.add<edm::InputTag>("PFTauTag", edm::InputTag("hpsPFTauProducer"));
0249   desc.add<edm::InputTag>("MuonTag", edm::InputTag("MyMuons"));
0250   desc.add<edm::InputTag>("PVTag", edm::InputTag("offlinePrimaryVertices"));
0251
0252   return desc;
0253 }