diff --git a/HLTrigger/Muon/BuildFile.xml b/HLTrigger/Muon/BuildFile.xml index 4efaf0ad087d0..d95974d4cf285 100644 --- a/HLTrigger/Muon/BuildFile.xml +++ b/HLTrigger/Muon/BuildFile.xml @@ -1,3 +1,4 @@ + diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h new file mode 100644 index 0000000000000..32277e6bde769 --- /dev/null +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -0,0 +1,259 @@ +#ifndef HLTrigger_Muon_HLTTriMuonIsolation_h +#define HLTrigger_Muon_HLTTriMuonIsolation_h + +#include +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/Candidate/interface/CompositeCandidate.h" +#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/LorentzVector.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.h" +#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h" +#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h" + + +class HLTTriMuonIsolation : public edm::global::EDProducer<> { + public: + explicit HLTTriMuonIsolation(const edm::ParameterSet& iConfig); + ~HLTTriMuonIsolation(); + virtual void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + const edm::EDGetTokenT L3MuonsToken_ ; + const edm::EDGetTokenT AllMuonsToken_ ; + const edm::EDGetTokenT L3DiMuonsFilterToken_; + const edm::EDGetTokenT IsoTracksToken_ ; + + edm::Handle L3MuCands ; + edm::Handle L3DiMuonsFilterCands; + edm::Handle PassedL3Muons ; + edm::Handle AllMuCands ; + edm::Handle IsoTracks ; + + static bool ptComparer(const reco::RecoChargedCandidate & mu_1, const reco::RecoChargedCandidate & mu_2) { return mu_1.pt() > mu_2.pt(); } + + const double Muon1PtCut_ ; + const double Muon2PtCut_ ; + const double Muon3PtCut_ ; + const double TriMuonPtCut_ ; + const double TriMuonEtaCut_ ; + const double ChargedRelIsoCut_; + const double ChargedAbsIsoCut_; + const double IsoConeSize_ ; + const double MatchingConeSize_; + const double MinTriMuonMass_ ; + const double MaxTriMuonMass_ ; + const double MaxTriMuonRadius_; + const int TriMuonAbsCharge_; + const double MaxDZ_ ; + const bool EnableRelIso_ ; + const bool EnableAbsIso_ ; +}; + +HLTTriMuonIsolation::HLTTriMuonIsolation(const edm::ParameterSet& iConfig): + L3MuonsToken_ (consumes (iConfig.getParameter("L3MuonsSrc" ))), + AllMuonsToken_ (consumes (iConfig.getParameter("AllMuonsSrc" ))), + L3DiMuonsFilterToken_(consumes (iConfig.getParameter("L3DiMuonsFilterSrc"))), + IsoTracksToken_ (consumes (iConfig.getParameter("IsoTracksSrc" ))), + Muon1PtCut_ (iConfig.getParameter ("Muon1PtCut" )) , + Muon2PtCut_ (iConfig.getParameter ("Muon2PtCut" )) , + Muon3PtCut_ (iConfig.getParameter ("Muon3PtCut" )) , + TriMuonPtCut_ (iConfig.getParameter ("TriMuonPtCut" )) , + TriMuonEtaCut_ (iConfig.getParameter ("TriMuonEtaCut" )) , + ChargedRelIsoCut_ (iConfig.getParameter ("ChargedRelIsoCut" )) , + ChargedAbsIsoCut_ (iConfig.getParameter ("ChargedAbsIsoCut" )) , + IsoConeSize_ (iConfig.getParameter ("IsoConeSize" )) , + MatchingConeSize_ (iConfig.getParameter ("MatchingConeSize" )) , + MinTriMuonMass_ (iConfig.getParameter ("MinTriMuonMass" )) , + MaxTriMuonMass_ (iConfig.getParameter ("MaxTriMuonMass" )) , + MaxTriMuonRadius_ (iConfig.getParameter ("MaxTriMuonRadius" )) , + TriMuonAbsCharge_ (iConfig.getParameter ("TriMuonAbsCharge" )) , + MaxDZ_ (iConfig.getParameter ("MaxDZ" )) , + EnableRelIso_ (iConfig.getParameter ("EnableRelIso" )) , + EnableAbsIso_ (iConfig.getParameter ("EnableAbsIso" )) +{ + //register products + produces("Taus"); + produces("SelectedTaus"); +} + +HLTTriMuonIsolation::~HLTTriMuonIsolation(){ } + +void +HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventSetup const & iSetup) const +{ + std::unique_ptr Taus (new reco::CompositeCandidateCollection); + std::unique_ptr SelectedTaus(new reco::CompositeCandidateCollection); + + // Get the L3 muon candidates + edm::Handle L3MuCands; + iEvent.getByToken(L3MuonsToken_, L3MuCands); + + // Get the L3 muon candidates that passed the filter + edm::Handle L3DiMuonsFilterCands; + iEvent.getByToken(L3DiMuonsFilterToken_, L3DiMuonsFilterCands); + + std::vector PassedL3Muons; + L3DiMuonsFilterCands->getObjects(trigger::TriggerMuon, PassedL3Muons); + + // Get the Trk + L3 muon candidates (after merging) + edm::Handle AllMuCands; + iEvent.getByToken(AllMuonsToken_, AllMuCands); + + // Get iso tracks + edm::Handle IsoTracks; + iEvent.getByToken(IsoTracksToken_, IsoTracks); + + if (AllMuCands->size() >= 3 && L3MuCands->size() >= 2){ + // Create the 3-muon candidates + // loop over L3/Trk muons and create all combinations + auto AllMuCands_end = AllMuCands->end(); + for (auto i = AllMuCands->begin(); i != AllMuCands_end-2; ++i) { + // check that muon_i passes the previous filter + bool passingPreviousFilter_1 = false; + for (const auto & imu : PassedL3Muons){ + if (reco::deltaR2(i->momentum(), imu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_1 = true; + } + for (auto j = i+1; j != AllMuCands_end-1; ++j) { + // check that muon_j passes the previous filter + bool passingPreviousFilter_2 = false; + for (const auto & jmu : PassedL3Muons){ + if (reco::deltaR2(j->momentum(), jmu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_2 = true; + } + // if, at this point, no muons passed the previous filter just skip to the next iteration + if (!(passingPreviousFilter_1 || passingPreviousFilter_2)) continue; + for (auto k = j+1; k != AllMuCands_end; ++k){ + // check that muon_k passes the previous filter + bool passingPreviousFilter_3 = false; + for (const auto & kmu : PassedL3Muons){ + if (reco::deltaR2(k->momentum(), kmu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_3 = true; + } + // at least two muons must have passed the previous di-muon filter + if (!( (passingPreviousFilter_1 & passingPreviousFilter_2 ) || + (passingPreviousFilter_1 & passingPreviousFilter_3 ) || + (passingPreviousFilter_2 & passingPreviousFilter_3 ) )) continue; + + // Create a composite candidate to be a tau + reco::CompositeCandidate tau; + + // sort the muons by pt and add them to the tau + reco::RecoChargedCandidateCollection daughters; + daughters.reserve(3); + + daughters.push_back(*i); + daughters.push_back(*j); + daughters.push_back(*k); + + std::sort(daughters.begin(), daughters.end(), ptComparer); + + tau.addDaughter((daughters)[0], "Muon_1"); + tau.addDaughter((daughters)[1], "Muon_2"); + tau.addDaughter((daughters)[2], "Muon_3"); + + // start building the tau + int charge = daughters[0].charge() + daughters[1].charge() + daughters[2].charge(); + math::XYZTLorentzVectorD taup4 = daughters[0].p4() + daughters[1].p4() + daughters[2].p4() ; + int tauPdgId = charge > 0? 15 : -15; + + tau.setP4(taup4); + tau.setCharge(charge); + tau.setPdgId(tauPdgId); + tau.setVertex((daughters)[0].vertex()); // assign the leading muon vertex as tau vertex + + // the three muons must be close to each other in Z + if (std::abs(tau.daughter(0)->vz() - tau.vz()) > MaxDZ_) continue; + if (std::abs(tau.daughter(1)->vz() - tau.vz()) > MaxDZ_) continue; + if (std::abs(tau.daughter(2)->vz() - tau.vz()) > MaxDZ_) continue; + + // require muons to be collimated + bool collimated = true; + for (auto const &idau : daughters){ + if (reco::deltaR2(tau.p4(), idau.p4()) > MaxTriMuonRadius_*MaxTriMuonRadius_) { + collimated = false; + break; + } + } + + if (!collimated) continue; + + // a good tau, at last + Taus->push_back(tau); + } + } + } + + // Loop over taus and further select + for (const auto & itau : *Taus){ + if ( itau.pt() < TriMuonPtCut_ ) continue; + if ( itau.mass() < MinTriMuonMass_) continue; + if ( itau.mass() > MaxTriMuonMass_) continue; + if (std::abs(itau.eta()) > TriMuonEtaCut_ ) continue; + if (itau.daughter(0)->pt() < Muon1PtCut_ ) continue; + if (itau.daughter(1)->pt() < Muon2PtCut_ ) continue; + if (itau.daughter(2)->pt() < Muon3PtCut_ ) continue; + if ((std::abs(itau.charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue; + if (std::abs(itau.daughter(0)->vz() - itau.vz()) > MaxDZ_) continue; + if (std::abs(itau.daughter(1)->vz() - itau.vz()) > MaxDZ_) continue; + if (std::abs(itau.daughter(2)->vz() - itau.vz()) > MaxDZ_) continue; + + // remove the candidate pt from the iso sum + double sumPt = -itau.pt(); + + // compute iso sum pT + for (const auto & itrk : *IsoTracks){ + if (reco::deltaR2(itrk.momentum(), itau.p4()) > IsoConeSize_*IsoConeSize_) continue; + if (std::abs(itrk.vz() - itau.vz()) > MaxDZ_) continue; + sumPt += itrk.pt(); + } + + // apply the isolation cut + if ((sumPt > (EnableAbsIso_ * ChargedAbsIsoCut_)) || + (sumPt > (EnableRelIso_ * ChargedRelIsoCut_ * itau.pt()))) continue; + + SelectedTaus->push_back(itau); + } + } + + // finally put the vector of 3-muon candidates in the event + iEvent.put(std::move(Taus) , "Taus" ); + iEvent.put(std::move(SelectedTaus), "SelectedTaus"); +} + +void +HLTTriMuonIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) +{ + edm::ParameterSetDescription desc; + desc.add("L3MuonsSrc" , edm::InputTag("hltIterL3FromL2MuonCandidates" )); + desc.add("AllMuonsSrc" , edm::InputTag("hltGlbTrkMuonCands" )); + desc.add("L3DiMuonsFilterSrc", edm::InputTag("hltDiMuonForTau3MuDzFiltered0p3")); + desc.add("IsoTracksSrc" , edm::InputTag("hltIter2L3FromL2MuonMerged" )); + desc.add("Muon1PtCut" , 5. ); + desc.add("Muon2PtCut" , 3. ); + desc.add("Muon3PtCut" , 0. ); + desc.add("TriMuonPtCut" , 8. ); + desc.add("TriMuonEtaCut" , 2.5 ); + desc.add("ChargedAbsIsoCut", 3.0 ); + desc.add("ChargedRelIsoCut", 0.1 ); + desc.add("IsoConeSize" , 0.5 ); + desc.add("MatchingConeSize", 0.03 ); + desc.add("MinTriMuonMass" , 0.5 ); + desc.add("MaxTriMuonMass" , 2.8 ); + desc.add("MaxTriMuonRadius", 0.6 ); + desc.add ("TriMuonAbsCharge", -1 ); + desc.add("MaxDZ" , 0.3 ); + desc.add ("EnableRelIso" , false); + desc.add ("EnableAbsIso" , true ); + descriptions.add("hltTriMuonIsolationProducer",desc); +} + +#endif diff --git a/HLTrigger/Muon/src/SealModule.cc b/HLTrigger/Muon/src/SealModule.cc index 5a25ae228768e..7d3b264c40c62 100644 --- a/HLTrigger/Muon/src/SealModule.cc +++ b/HLTrigger/Muon/src/SealModule.cc @@ -45,3 +45,6 @@ DEFINE_FWK_MODULE(HLTL1MuonSelector); DEFINE_FWK_MODULE(HLTL1TMuonSelector); DEFINE_FWK_MODULE(HLTScoutingMuonProducer); DEFINE_FWK_MODULE(HLTL1MuonNoL2Selector); + +#include "HLTrigger/Muon/interface/HLTTriMuonIsolation.h" +DEFINE_FWK_MODULE(HLTTriMuonIsolation);