File indexing completed on 2023-03-17 11:09:42
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;
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
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
0100 edm::Handle<reco::RecoChargedCandidateCollection> L3MuCands;
0101 iEvent.getByToken(L3MuonsToken_, L3MuCands);
0102
0103
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
0111 edm::Handle<reco::RecoChargedCandidateCollection> AllMuCands;
0112 iEvent.getByToken(AllMuonsToken_, AllMuCands);
0113
0114
0115 edm::Handle<reco::TrackCollection> IsoTracks;
0116 iEvent.getByToken(IsoTracksToken_, IsoTracks);
0117
0118 if (AllMuCands->size() >= 3 && L3MuCands->size() >= 2) {
0119
0120
0121 auto AllMuCands_end = AllMuCands->end();
0122 for (auto i = AllMuCands->begin(); i != AllMuCands_end - 2; ++i) {
0123
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
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
0137 if (!(passingPreviousFilter_1 || passingPreviousFilter_2))
0138 continue;
0139 for (auto k = j + 1; k != AllMuCands_end; ++k) {
0140
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
0147 if (!((passingPreviousFilter_1 & passingPreviousFilter_2) ||
0148 (passingPreviousFilter_1 & passingPreviousFilter_3) ||
0149 (passingPreviousFilter_2 & passingPreviousFilter_3)))
0150 continue;
0151
0152
0153 reco::CompositeCandidate tau;
0154
0155
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
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
0174 tau.addDaughter((daughters)[0], "Muon_1");
0175 tau.addDaughter((daughters)[1], "Muon_2");
0176 tau.addDaughter((daughters)[2], "Muon_3");
0177
0178
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());
0187
0188
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
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
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
0219 if ((std::abs(tau.charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0))
0220 continue;
0221
0222
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
0231 Taus->push_back(tau);
0232 }
0233 }
0234 }
0235
0236
0237 std::sort(Taus->begin(), Taus->end(), ptComparer<reco::CompositeCandidate>);
0238
0239
0240 for (const auto& itau : *Taus) {
0241
0242 double sumPt = -itau.pt();
0243
0244
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
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
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