SoftLepton/plugins/SoftPFMuonTagInfoProducer.cc

0001 // * Author: Alberto Zucchetta
0002 // * Mail: a.zucchetta@cern.ch
0003 // * January 16, 2015
0004
0005 #include "FWCore/Framework/interface/ModuleFactory.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "DataFormats/Common/interface/ValueMap.h"
0013 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0014 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0015 #include "DataFormats/MuonReco/interface/Muon.h"
0016 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0017 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0018 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0019 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0020 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0021 #include "DataFormats/VertexReco/interface/Vertex.h"
0022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0023
0024 #include "DataFormats/BTauReco/interface/SoftLeptonTagInfo.h"
0025 #include "DataFormats/BTauReco/interface/CandSoftLeptonTagInfo.h"
0026
0027 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0028 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0029 #include "DataFormats/Math/interface/deltaR.h"
0030 #include "DataFormats/Math/interface/Vector3D.h"
0031 #include "DataFormats/Math/interface/Point3D.h"
0032 #include "FWCore/Utilities/interface/EDMException.h"
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0035 #include "DataFormats/JetReco/interface/Jet.h"
0036 #include "DataFormats/PatCandidates/interface/Jet.h"
0037 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0038 #include "DataFormats/PatCandidates/interface/Muon.h"
0039 // Muons
0040 #include "DataFormats/MuonReco/interface/Muon.h"
0041 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0042 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0043 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0044 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0045
0046 #include "DataFormats/BTauReco/interface/SoftLeptonTagInfo.h"
0047
0048 // Transient Track and IP
0049 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0050 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0051 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0052 #include "TrackingTools/IPTools/interface/IPTools.h"
0053 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0054 #include <cmath>
0055 #include <vector>
0056
0057 class SoftPFMuonTagInfoProducer : public edm::global::EDProducer<> {
0058 public:
0059   SoftPFMuonTagInfoProducer(const edm::ParameterSet& conf);
0060
0061   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0062
0063 private:
0064   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const final;
0065   float boostedPPar(const math::XYZVector&, const math::XYZVector&) const;
0066
0067   const edm::EDGetTokenT<edm::View<reco::Jet> > jetToken;
0068   const edm::EDGetTokenT<edm::View<reco::Muon> > muonToken;
0069   const edm::EDGetTokenT<reco::VertexCollection> vertexToken;
0070   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> builderToken;
0071   const edm::EDPutTokenT<reco::CandSoftLeptonTagInfoCollection> putToken;
0072   const float pTcut, SIPsigcut, IPsigcut, ratio1cut, ratio2cut;
0073   const bool useFilter;
0074 };
0075
0076 SoftPFMuonTagInfoProducer::SoftPFMuonTagInfoProducer(const edm::ParameterSet& conf)
0077     : jetToken(consumes(conf.getParameter<edm::InputTag>("jets"))),
0078       muonToken(consumes(conf.getParameter<edm::InputTag>("muons"))),
0079       vertexToken(consumes(conf.getParameter<edm::InputTag>("primaryVertex"))),
0080       builderToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0081       putToken(produces()),
0082       pTcut(conf.getParameter<double>("muonPt")),
0083       SIPsigcut(conf.getParameter<double>("muonSIPsig")),
0084       IPsigcut(conf.getParameter<double>("filterIpsig")),
0085       ratio1cut(conf.getParameter<double>("filterRatio1")),
0086       ratio2cut(conf.getParameter<double>("filterRatio2")),
0087       useFilter(conf.getParameter<bool>("filterPromptMuons")) {}
0088
0089 void SoftPFMuonTagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0090   edm::ParameterSetDescription desc;
0091   desc.add<edm::InputTag>("jets");
0092   desc.add<edm::InputTag>("muons");
0093   desc.add<edm::InputTag>("primaryVertex");
0094   desc.add<double>("muonPt");
0095   desc.add<double>("muonSIPsig");
0096   desc.add<double>("filterIpsig");
0097   desc.add<double>("filterRatio1");
0098   desc.add<double>("filterRatio2");
0099   desc.add<bool>("filterPromptMuons");
0100   descriptions.addDefault(desc);
0101 }
0102
0103 void SoftPFMuonTagInfoProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0104   // Declare produced collection
0105   reco::CandSoftLeptonTagInfoCollection theMuonTagInfo;
0106
0107   // Declare and open Jet collection
0108   edm::Handle<edm::View<reco::Jet> > theJetCollection = iEvent.getHandle(jetToken);
0109
0110   // Declare Muon collection
0111   edm::Handle<edm::View<reco::Muon> > theMuonCollection = iEvent.getHandle(muonToken);
0112
0113   // Declare and open Vertex collection
0114   edm::Handle<reco::VertexCollection> theVertexCollection = iEvent.getHandle(vertexToken);
0115   if (!theVertexCollection.isValid() || theVertexCollection->empty()) {
0116     iEvent.emplace(putToken, std::move(theMuonTagInfo));
0117     return;
0118   }
0119   const reco::Vertex* vertex = &theVertexCollection->front();
0120
0121   // Build TransientTrackBuilder
0122   const TransientTrackBuilder& transientTrackBuilder = iSetup.getData(builderToken);
0123
0124   // Loop on jets
0125   for (unsigned int ij = 0, nj = theJetCollection->size(); ij < nj; ij++) {
0126     edm::RefToBase<reco::Jet> jetRef = theJetCollection->refAt(ij);
0127     // Build TagInfo object
0128     reco::CandSoftLeptonTagInfo tagInfo;
0129     tagInfo.setJetRef(jetRef);
0130     // Loop on jet daughters
0131     for (unsigned int id = 0, nd = jetRef->numberOfDaughters(); id < nd; ++id) {
0132       edm::Ptr<reco::Candidate> lepPtr = jetRef->daughterPtr(id);
0133       if (std::abs(lepPtr->pdgId()) != 13)
0134         continue;
0135
0136       const reco::Muon* muon(nullptr);
0137       // Step 1: try to access the muon from reco::PFCandidate
0138       const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(lepPtr.get());
0139       if (pfcand) {
0140         muon = pfcand->muonRef().get();
0141       }
0142       // If not PFCandidate is available, find a match looping on the muon collection
0143       else {
0144         for (unsigned int im = 0, nm = theMuonCollection->size(); im < nm; ++im) {  // --- Begin loop on muons
0145           const reco::Muon* recomuon = &theMuonCollection->at(im);
0146           const pat::Muon* patmuon = dynamic_cast<const pat::Muon*>(recomuon);
0147           // Step 2: try a match between reco::Candidate
0148           if (patmuon) {
0149             if (patmuon->originalObjectRef() == lepPtr) {
0150               muon = theMuonCollection->refAt(im).get();
0151               break;
0152             }
0153           }
0154           // Step 3: try a match with dR and dpT if pat::Muon casting fails
0155           else {
0156             if (reco::deltaR(*recomuon, *lepPtr) < 0.01 &&
0157                 std::abs(recomuon->pt() - lepPtr->pt()) / lepPtr->pt() < 0.1) {
0158               muon = theMuonCollection->refAt(im).get();
0159               break;
0160             }
0161           }
0162         }  // --- End loop on muons
0163       }
0164       if (!muon || !muon::isLooseMuon(*muon) || muon->pt() < pTcut)
0165         continue;
0166       reco::TrackRef trkRef(muon->innerTrack());
0167       reco::TrackBaseRef trkBaseRef(trkRef);
0168       // Build Transient Track
0169       reco::TransientTrack transientTrack = transientTrackBuilder.build(trkRef);
0170       // Define jet and muon vectors
0171       reco::Candidate::Vector jetvect(jetRef->p4().Vect()), muonvect(muon->p4().Vect());
0172       // Calculate variables
0173       reco::SoftLeptonProperties properties;
0174       Measurement1D ip2d = IPTools::signedTransverseImpactParameter(
0175                                transientTrack, GlobalVector(jetRef->px(), jetRef->py(), jetRef->pz()), *vertex)
0176                                .second;
0177       Measurement1D ip3d = IPTools::signedImpactParameter3D(
0178                                transientTrack, GlobalVector(jetRef->px(), jetRef->py(), jetRef->pz()), *vertex)
0179                                .second;
0180       properties.sip2dsig = ip2d.significance();
0181       properties.sip3dsig = ip3d.significance();
0182       properties.sip2d = ip2d.value();
0183       properties.sip3d = ip3d.value();
0184       properties.deltaR = reco::deltaR(*jetRef, *muon);
0185       properties.ptRel = ((jetvect - muonvect).Cross(muonvect)).R() / jetvect.R();  // | (Pj-Pu) X Pu | / | Pj |
0186       float mag = muonvect.R() * jetvect.R();
0187       float dot = muon->p4().Dot(jetRef->p4());
0188       properties.etaRel = -log((mag - dot) / (mag + dot)) / 2.;
0189       properties.ratio = muon->pt() / jetRef->pt();
0190       properties.ratioRel = muon->p4().Dot(jetRef->p4()) / jetvect.Mag2();
0191       properties.p0Par = boostedPPar(muon->momentum(), jetRef->momentum());
0192       properties.charge = muon->charge();
0193
0194       if (std::abs(properties.sip3dsig) > SIPsigcut)
0195         continue;
0196
0197       // Filter leptons from W, Z decays
0198       if (useFilter &&
0199           ((std::abs(properties.sip3dsig) < IPsigcut && properties.ratio > ratio1cut) || properties.ratio > ratio2cut))
0200         continue;
0201
0202       // Insert lepton properties
0203       tagInfo.insert(lepPtr, properties);
0204
0205     }  // --- End loop on daughters
0206
0207     // Fill the TagInfo collection
0208     theMuonTagInfo.push_back(tagInfo);
0209   }  // --- End loop on jets
0210
0211   // Put the TagInfo collection in the event
0212   iEvent.emplace(putToken, std::move(theMuonTagInfo));
0213 }
0214
0215 // compute the lepton momentum along the jet axis, in the jet rest frame
0216 float SoftPFMuonTagInfoProducer::boostedPPar(const math::XYZVector& vector, const math::XYZVector& axis) const {
0217   static const double lepton_mass = 0.00;  // assume a massless (ultrarelativistic) lepton
0218   static const double jet_mass = 5.279;    // use B±/B0 mass as the jet rest mass [PDG 2007 updates]
0219   ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > lepton(
0220       vector.Dot(axis) / axis.r(), ROOT::Math::VectorUtil::Perp(vector, axis), 0., lepton_mass);
0221   ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > jet(axis.r(), 0., 0., jet_mass);
0222   ROOT::Math::BoostX boost(-jet.Beta());
0223   return boost(lepton).x();
0224 }
0225
0226 DEFINE_FWK_MODULE(SoftPFMuonTagInfoProducer);