Muon/plugins/HLTTriMuonIsolation.h

0001 #ifndef HLTrigger_Muon_HLTTriMuonIsolation_h
0002 #define HLTrigger_Muon_HLTTriMuonIsolation_h
0003
0004 #include <iostream>
0005 #include <string>
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013
0014 #include "DataFormats/Candidate/interface/CompositeCandidate.h"
0015 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0016 #include "DataFormats/Math/interface/deltaR.h"
0017 #include "DataFormats/Math/interface/LorentzVector.h"
0018 #include "DataFormats/MuonReco/interface/Muon.h"
0019 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0020 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0021 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0022
0023 class HLTTriMuonIsolation : public edm::global::EDProducer<> {
0024 public:
0025   explicit HLTTriMuonIsolation(const edm::ParameterSet& iConfig);
0026   ~HLTTriMuonIsolation() override = default;
0027   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0028   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0029
0030 private:
0031   const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> L3MuonsToken_;
0032   const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> AllMuonsToken_;
0033   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterToken_;
0034   const edm::EDGetTokenT<reco::TrackCollection> IsoTracksToken_;
0035
0036   edm::Handle<reco::RecoChargedCandidateCollection> L3MuCands;
0037   edm::Handle<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterCands;
0038   edm::Handle<reco::RecoChargedCandidateRef> PassedL3Muons;
0039   edm::Handle<reco::RecoChargedCandidateCollection> AllMuCands;
0040   edm::Handle<reco::TrackCollection> IsoTracks;
0041
0042   template <typename T>
0043   static bool ptComparer(const T& cand_1, const T& cand_2) {
0044     return cand_1.pt() > cand_2.pt();
0045   }
0046
0047   const double TwiceMuonMass_ = 2. * 0.1056583715;  // in GeV
0048
0049   const double Muon1PtCut_;
0050   const double Muon2PtCut_;
0051   const double Muon3PtCut_;
0052   const double TriMuonPtCut_;
0053   const double TriMuonEtaCut_;
0054   const double ChargedRelIsoCut_;
0055   const double ChargedAbsIsoCut_;
0056   const double IsoConeSize_;
0057   const double MatchingConeSize_;
0058   const double MinTriMuonMass_;
0059   const double MaxTriMuonMass_;
0060   const double MaxTriMuonRadius_;
0061   const int TriMuonAbsCharge_;
0062   const double MaxDZ_;
0063   const bool EnableRelIso_;
0064   const bool EnableAbsIso_;
0065 };
0066
0067 inline HLTTriMuonIsolation::HLTTriMuonIsolation(const edm::ParameterSet& iConfig)
0068     : L3MuonsToken_(consumes<reco::RecoChargedCandidateCollection>(iConfig.getParameter<edm::InputTag>("L3MuonsSrc"))),
0069       AllMuonsToken_(
0070           consumes<reco::RecoChargedCandidateCollection>(iConfig.getParameter<edm::InputTag>("AllMuonsSrc"))),
0071       L3DiMuonsFilterToken_(
0072           consumes<trigger::TriggerFilterObjectWithRefs>(iConfig.getParameter<edm::InputTag>("L3DiMuonsFilterSrc"))),
0073       IsoTracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("IsoTracksSrc"))),
0074       Muon1PtCut_(iConfig.getParameter<double>("Muon1PtCut")),
0075       Muon2PtCut_(iConfig.getParameter<double>("Muon2PtCut")),
0076       Muon3PtCut_(iConfig.getParameter<double>("Muon3PtCut")),
0077       TriMuonPtCut_(iConfig.getParameter<double>("TriMuonPtCut")),
0078       TriMuonEtaCut_(iConfig.getParameter<double>("TriMuonEtaCut")),
0079       ChargedRelIsoCut_(iConfig.getParameter<double>("ChargedRelIsoCut")),
0080       ChargedAbsIsoCut_(iConfig.getParameter<double>("ChargedAbsIsoCut")),
0081       IsoConeSize_(iConfig.getParameter<double>("IsoConeSize")),
0082       MatchingConeSize_(iConfig.getParameter<double>("MatchingConeSize")),
0083       MinTriMuonMass_(iConfig.getParameter<double>("MinTriMuonMass")),
0084       MaxTriMuonMass_(iConfig.getParameter<double>("MaxTriMuonMass")),
0085       MaxTriMuonRadius_(iConfig.getParameter<double>("MaxTriMuonRadius")),
0086       TriMuonAbsCharge_(iConfig.getParameter<int>("TriMuonAbsCharge")),
0087       MaxDZ_(iConfig.getParameter<double>("MaxDZ")),
0088       EnableRelIso_(iConfig.getParameter<bool>("EnableRelIso")),
0089       EnableAbsIso_(iConfig.getParameter<bool>("EnableAbsIso")) {
0090   //register products
0091   produces<reco::CompositeCandidateCollection>("Taus");
0092   produces<reco::CompositeCandidateCollection>("SelectedTaus");
0093 }
0094
0095 inline void HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
0096   std::unique_ptr<reco::CompositeCandidateCollection> Taus(new reco::CompositeCandidateCollection);
0097   std::unique_ptr<reco::CompositeCandidateCollection> SelectedTaus(new reco::CompositeCandidateCollection);
0098
0099   // Get the L3 muon candidates
0100   edm::Handle<reco::RecoChargedCandidateCollection> L3MuCands;
0101   iEvent.getByToken(L3MuonsToken_, L3MuCands);
0102
0103   // Get the L3 muon candidates that passed the filter
0104   edm::Handle<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterCands;
0105   iEvent.getByToken(L3DiMuonsFilterToken_, L3DiMuonsFilterCands);
0106
0107   std::vector<reco::RecoChargedCandidateRef> PassedL3Muons;
0108   L3DiMuonsFilterCands->getObjects(trigger::TriggerMuon, PassedL3Muons);
0109
0110   // Get the Trk + L3 muon candidates (after merging)
0111   edm::Handle<reco::RecoChargedCandidateCollection> AllMuCands;
0112   iEvent.getByToken(AllMuonsToken_, AllMuCands);
0113
0114   // Get iso tracks
0115   edm::Handle<reco::TrackCollection> IsoTracks;
0116   iEvent.getByToken(IsoTracksToken_, IsoTracks);
0117
0118   if (AllMuCands->size() >= 3 && L3MuCands->size() >= 2) {
0119     // Create the 3-muon candidates
0120     // loop over L3/Trk muons and create all combinations
0121     auto AllMuCands_end = AllMuCands->end();
0122     for (auto i = AllMuCands->begin(); i != AllMuCands_end - 2; ++i) {
0123       // check that muon_i passes the previous filter
0124       bool passingPreviousFilter_1 = false;
0125       for (const auto& imu : PassedL3Muons) {
0126         if (reco::deltaR2(i->momentum(), imu->momentum()) < (MatchingConeSize_ * MatchingConeSize_))
0127           passingPreviousFilter_1 = true;
0128       }
0129       for (auto j = i + 1; j != AllMuCands_end - 1; ++j) {
0130         // check that muon_j passes the previous filter
0131         bool passingPreviousFilter_2 = false;
0132         for (const auto& jmu : PassedL3Muons) {
0133           if (reco::deltaR2(j->momentum(), jmu->momentum()) < (MatchingConeSize_ * MatchingConeSize_))
0134             passingPreviousFilter_2 = true;
0135         }
0136         // if, at this point, no muons passed the previous filter just skip to the next iteration
0137         if (!(passingPreviousFilter_1 || passingPreviousFilter_2))
0138           continue;
0139         for (auto k = j + 1; k != AllMuCands_end; ++k) {
0140           // check that muon_k passes the previous filter
0141           bool passingPreviousFilter_3 = false;
0142           for (const auto& kmu : PassedL3Muons) {
0143             if (reco::deltaR2(k->momentum(), kmu->momentum()) < (MatchingConeSize_ * MatchingConeSize_))
0144               passingPreviousFilter_3 = true;
0145           }
0146           // at least two muons must have passed the previous di-muon filter
0147           if (!((passingPreviousFilter_1 & passingPreviousFilter_2) ||
0148                 (passingPreviousFilter_1 & passingPreviousFilter_3) ||
0149                 (passingPreviousFilter_2 & passingPreviousFilter_3)))
0150             continue;
0151
0152           // Create a composite candidate to be a tau
0153           reco::CompositeCandidate tau;
0154
0155           // sort the muons by pt and add them to the tau
0156           reco::RecoChargedCandidateCollection daughters;
0157           daughters.reserve(3);
0158
0159           daughters.push_back(*i);
0160           daughters.push_back(*j);
0161           daughters.push_back(*k);
0162
0163           std::sort(daughters.begin(), daughters.end(), ptComparer<reco::RecoChargedCandidate>);
0164
0165           // Muon kinematic selections
0166           if (daughters[0].pt() < Muon1PtCut_)
0167             continue;
0168           if (daughters[1].pt() < Muon2PtCut_)
0169             continue;
0170           if (daughters[2].pt() < Muon3PtCut_)
0171             continue;
0172
0173           // assign the tau its daughters
0174           tau.addDaughter((daughters)[0], "Muon_1");
0175           tau.addDaughter((daughters)[1], "Muon_2");
0176           tau.addDaughter((daughters)[2], "Muon_3");
0177
0178           // start building the tau
0179           int charge = daughters[0].charge() + daughters[1].charge() + daughters[2].charge();
0180           math::XYZTLorentzVectorD taup4 = daughters[0].p4() + daughters[1].p4() + daughters[2].p4();
0181           int tauPdgId = charge > 0 ? 15 : -15;
0182
0183           tau.setP4(taup4);
0184           tau.setCharge(charge);
0185           tau.setPdgId(tauPdgId);
0186           tau.setVertex((daughters)[0].vertex());  // assign the leading muon vertex as tau vertex
0187
0188           // the three muons must be close to each other in Z
0189           if (std::abs(tau.daughter(0)->vz() - tau.vz()) > MaxDZ_)
0190             continue;
0191           if (std::abs(tau.daughter(1)->vz() - tau.vz()) > MaxDZ_)
0192             continue;
0193           if (std::abs(tau.daughter(2)->vz() - tau.vz()) > MaxDZ_)
0194             continue;
0195
0196           // require muons to be collimated
0197           bool collimated = true;
0198           for (auto const& idau : daughters) {
0199             if (reco::deltaR2(tau.p4(), idau.p4()) > MaxTriMuonRadius_ * MaxTriMuonRadius_) {
0200               collimated = false;
0201               break;
0202             }
0203           }
0204
0205           if (!collimated)
0206             continue;
0207
0208           // Tau kinematic selections
0209           if (tau.pt() < TriMuonPtCut_)
0210             continue;
0211           if (tau.mass() < MinTriMuonMass_)
0212             continue;
0213           if (tau.mass() > MaxTriMuonMass_)
0214             continue;
0215           if (std::abs(tau.eta()) > TriMuonEtaCut_)
0216             continue;
0217
0218           // Tau charge selection
0219           if ((std::abs(tau.charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0))
0220             continue;
0221
0222           // Sanity check against duplicates, di-muon masses must be > 2 * mass_mu
0223           if ((tau.daughter(0)->p4() + tau.daughter(1)->p4()).mass() < TwiceMuonMass_)
0224             continue;
0225           if ((tau.daughter(0)->p4() + tau.daughter(2)->p4()).mass() < TwiceMuonMass_)
0226             continue;
0227           if ((tau.daughter(1)->p4() + tau.daughter(2)->p4()).mass() < TwiceMuonMass_)
0228             continue;
0229
0230           // a good tau, at last
0231           Taus->push_back(tau);
0232         }
0233       }
0234     }
0235
0236     // Sort taus by pt
0237     std::sort(Taus->begin(), Taus->end(), ptComparer<reco::CompositeCandidate>);
0238
0239     // Loop over taus and further select by isolation
0240     for (const auto& itau : *Taus) {
0241       // remove the candidate pt from the iso sum
0242       double sumPt = -itau.pt();
0243
0244       // compute iso sum pT
0245       for (const auto& itrk : *IsoTracks) {
0246         if (reco::deltaR2(itrk.momentum(), itau.p4()) > IsoConeSize_ * IsoConeSize_)
0247           continue;
0248         if (std::abs(itrk.vz() - itau.vz()) > MaxDZ_)
0249           continue;
0250         sumPt += itrk.pt();
0251       }
0252
0253       // apply the isolation cut
0254       double chAbsIsoCut = EnableAbsIso_ ? ChargedAbsIsoCut_ : std::numeric_limits<double>::infinity();
0255       double chRelIsoCut = EnableRelIso_ ? ChargedRelIsoCut_ * itau.pt() : std::numeric_limits<double>::infinity();
0256
0257       if (!((sumPt < chAbsIsoCut) || (sumPt < chRelIsoCut)))
0258         continue;
0259
0260       SelectedTaus->push_back(itau);
0261     }
0262   }
0263
0264   // finally put the vector of 3-muon candidates in the event
0265   iEvent.put(std::move(Taus), "Taus");
0266   iEvent.put(std::move(SelectedTaus), "SelectedTaus");
0267 }
0268
0269 inline void HLTTriMuonIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0270   edm::ParameterSetDescription desc;
0271   desc.add<edm::InputTag>("L3MuonsSrc", edm::InputTag("hltIterL3FromL2MuonCandidates"));
0272   desc.add<edm::InputTag>("AllMuonsSrc", edm::InputTag("hltGlbTrkMuonCands"));
0273   desc.add<edm::InputTag>("L3DiMuonsFilterSrc", edm::InputTag("hltDiMuonForTau3MuDzFiltered0p3"));
0274   desc.add<edm::InputTag>("IsoTracksSrc", edm::InputTag("hltIter2L3FromL2MuonMerged"));
0275   desc.add<double>("Muon1PtCut", 5.);
0276   desc.add<double>("Muon2PtCut", 3.);
0277   desc.add<double>("Muon3PtCut", 0.);
0278   desc.add<double>("TriMuonPtCut", 8.);
0279   desc.add<double>("TriMuonEtaCut", 2.5);
0280   desc.add<double>("ChargedAbsIsoCut", 3.0);
0281   desc.add<double>("ChargedRelIsoCut", 0.1);
0282   desc.add<double>("IsoConeSize", 0.5);
0283   desc.add<double>("MatchingConeSize", 0.03);
0284   desc.add<double>("MinTriMuonMass", 0.5);
0285   desc.add<double>("MaxTriMuonMass", 2.8);
0286   desc.add<double>("MaxTriMuonRadius", 0.6);
0287   desc.add<int>("TriMuonAbsCharge", -1);
0288   desc.add<double>("MaxDZ", 0.3);
0289   desc.add<bool>("EnableRelIso", false);
0290   desc.add<bool>("EnableAbsIso", true);
0291   descriptions.add("hltTriMuonIsolationProducer", desc);
0292 }
0293
0294 #endif