From e0b3aa81891a4b3283be4063154a973b06f5babf Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Wed, 12 Jul 2017 09:12:21 +0200 Subject: [PATCH 1/7] add tau3mu HLT producer --- HLTrigger/Muon/BuildFile.xml | 1 + .../Muon/interface/HLTTriMuonIsolation.h | 230 ++++++++++++++++++ HLTrigger/Muon/src/SealModule.cc | 3 + 3 files changed, 234 insertions(+) create mode 100644 HLTrigger/Muon/interface/HLTTriMuonIsolation.h 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..ba5619cbf604d --- /dev/null +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -0,0 +1,230 @@ +#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 MinTriMuonMass_ ; + const double MaxTriMuonMass_ ; + 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" )) , + MinTriMuonMass_ (iConfig.getParameter ("MinTriMuonMass" )) , + MaxTriMuonMass_ (iConfig.getParameter ("MaxTriMuonMass" )) , + 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){ + // Create the 3-muon candidates + // loop over L3/Trk muons and create all combinations + for (unsigned int i = 0; i != AllMuCands->size(); ++i){ + const reco::TrackRef &tk_i = (*AllMuCands)[i].track(); + for (unsigned int j = i+1; j != AllMuCands->size(); ++j){ + const reco::TrackRef &tk_j = (*AllMuCands)[j].track(); + for (unsigned int k = j+1; k != AllMuCands->size(); ++k){ + const reco::TrackRef &tk_k = (*AllMuCands)[k].track(); + + int passingPreviousFilter_1 = 0; + int passingPreviousFilter_2 = 0; + int passingPreviousFilter_3 = 0; + + for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ + reco::TrackRef candTrkRef = (*imu)->get(); + if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1++; + if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2++; + if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3++; + } + // at least two muons must have passed the previous di-muon filter + if (((passingPreviousFilter_1 < 1) & (passingPreviousFilter_2 < 1)) || + ((passingPreviousFilter_1 < 1) & (passingPreviousFilter_3 < 1)) || + ((passingPreviousFilter_2 < 1) & (passingPreviousFilter_3 < 1))) 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.push_back((*AllMuCands)[i]); + Daughters.push_back((*AllMuCands)[j]); + Daughters.push_back((*AllMuCands)[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 + + Taus->push_back(Tau); + } + } + } + + // Loop over taus and further select + for (reco::CompositeCandidateCollection::const_iterator itau = Taus->begin(); itau != Taus->end(); ++itau){ + if ( itau->pt() < TriMuonPtCut_ ) continue; + if ( itau->mass() < MinTriMuonMass_) continue; + if ( itau->mass() > MaxTriMuonMass_) continue; + if (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 ((abs(itau->charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue; + if (abs(itau->daughter(0)->vz() - itau->vz()) > MaxDZ_) continue; + if (abs(itau->daughter(1)->vz() - itau->vz()) > MaxDZ_) continue; + if (abs(itau->daughter(2)->vz() - itau->vz()) > MaxDZ_) continue; + + double sumPt = -itau->pt(); // remove the candidate pt from the iso sum + + // compute iso sum pT + for (reco::TrackCollection::const_iterator itrk = IsoTracks->begin(); itrk != IsoTracks->end(); ++itrk){ + if (reco::deltaR2(itrk->momentum(), itau->p4()) > IsoConeSize_*IsoConeSize_) continue; + if (abs(itrk->vz() - itau->vz()) > MaxDZ_) continue; + sumPt += itrk->pt(); + } + + // apply the isolation cut + if ((std::max(0., sumPt) > (EnableAbsIso_ * ChargedAbsIsoCut_)) || + (std::max(0., 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("MinTriMuonMass" , 0.5 ); + desc.add("MaxTriMuonMass" , 2.8 ); + 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); From 3a1a02a1dfa1d167dcdb9b93c146c7a2aa6cf90d Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Fri, 14 Jul 2017 11:42:44 +0200 Subject: [PATCH 2/7] improve loop logic, anticipate dZ cut --- .../Muon/interface/HLTTriMuonIsolation.h | 80 +++++++++++-------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index ba5619cbf604d..bdcdbdad171eb 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -3,6 +3,7 @@ #include #include +#include #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/global/EDProducer.h" @@ -110,30 +111,38 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS edm::Handle IsoTracks; iEvent.getByToken(IsoTracksToken_, IsoTracks); - if (AllMuCands->size() >= 3){ + if (AllMuCands->size() >= 3 && L3MuCands->size() >= 2){ // Create the 3-muon candidates // loop over L3/Trk muons and create all combinations - for (unsigned int i = 0; i != AllMuCands->size(); ++i){ - const reco::TrackRef &tk_i = (*AllMuCands)[i].track(); - for (unsigned int j = i+1; j != AllMuCands->size(); ++j){ - const reco::TrackRef &tk_j = (*AllMuCands)[j].track(); - for (unsigned int k = j+1; k != AllMuCands->size(); ++k){ - const reco::TrackRef &tk_k = (*AllMuCands)[k].track(); - - int passingPreviousFilter_1 = 0; - int passingPreviousFilter_2 = 0; - int passingPreviousFilter_3 = 0; - + auto AllMuCands_end = AllMuCands->end(); + for (auto i = AllMuCands->begin(); i != AllMuCands_end; ++i) { + // check that muon_i passes the previous filter + bool passingPreviousFilter_1 = false; + const reco::TrackRef &tk_i = i->track(); + for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ + reco::TrackRef candTrkRef = (*imu)->get(); + if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1 = true; + } + for (auto j = i+1; j != AllMuCands_end; ++j) { + // check that muon_j passes the previous filter + bool passingPreviousFilter_2 = false; + const reco::TrackRef &tk_j = j->track(); + for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ + reco::TrackRef candTrkRef = (*imu)->get(); + if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2 = true; + } + for (auto k = j+1; k != AllMuCands_end; ++k){ + // check that muon_k passes the previous filter + bool passingPreviousFilter_3 = false; + const reco::TrackRef &tk_k = k->track(); for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1++; - if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2++; - if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3++; + if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3 = true; } // at least two muons must have passed the previous di-muon filter - if (((passingPreviousFilter_1 < 1) & (passingPreviousFilter_2 < 1)) || - ((passingPreviousFilter_1 < 1) & (passingPreviousFilter_3 < 1)) || - ((passingPreviousFilter_2 < 1) & (passingPreviousFilter_3 < 1))) continue; + 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; @@ -141,9 +150,9 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS // sort the muons by pt and add them to the tau reco::RecoChargedCandidateCollection Daughters; - Daughters.push_back((*AllMuCands)[i]); - Daughters.push_back((*AllMuCands)[j]); - Daughters.push_back((*AllMuCands)[k]); + Daughters.push_back(*i); + Daughters.push_back(*j); + Daughters.push_back(*k); std::sort(Daughters.begin(), Daughters.end(), ptComparer); @@ -160,6 +169,11 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS 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 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; Taus->push_back(Tau); } @@ -168,24 +182,22 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS // Loop over taus and further select for (reco::CompositeCandidateCollection::const_iterator itau = Taus->begin(); itau != Taus->end(); ++itau){ - if ( itau->pt() < TriMuonPtCut_ ) continue; - if ( itau->mass() < MinTriMuonMass_) continue; - if ( itau->mass() > MaxTriMuonMass_) continue; - if (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 ((abs(itau->charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue; - if (abs(itau->daughter(0)->vz() - itau->vz()) > MaxDZ_) continue; - if (abs(itau->daughter(1)->vz() - itau->vz()) > MaxDZ_) continue; - if (abs(itau->daughter(2)->vz() - itau->vz()) > MaxDZ_) continue; + 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; - double sumPt = -itau->pt(); // remove the candidate pt from the iso sum + // remove the candidate pt from the iso sum + double sumPt = -itau->pt(); // compute iso sum pT for (reco::TrackCollection::const_iterator itrk = IsoTracks->begin(); itrk != IsoTracks->end(); ++itrk){ if (reco::deltaR2(itrk->momentum(), itau->p4()) > IsoConeSize_*IsoConeSize_) continue; - if (abs(itrk->vz() - itau->vz()) > MaxDZ_) continue; + if (std::abs(itrk->vz() - itau->vz()) > MaxDZ_) continue; sumPt += itrk->pt(); } From 049a2d38444a8f8cc6227f6ee8a0c5d207427023 Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Fri, 14 Jul 2017 12:09:43 +0200 Subject: [PATCH 3/7] spare a loop if you can --- HLTrigger/Muon/interface/HLTTriMuonIsolation.h | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index bdcdbdad171eb..944af835283d4 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -51,6 +51,7 @@ class HLTTriMuonIsolation : public edm::global::EDProducer<> { const double ChargedRelIsoCut_; const double ChargedAbsIsoCut_; const double IsoConeSize_ ; + const double MatchingConeSize_; const double MinTriMuonMass_ ; const double MaxTriMuonMass_ ; const int TriMuonAbsCharge_; @@ -72,6 +73,7 @@ HLTTriMuonIsolation::HLTTriMuonIsolation(const edm::ParameterSet& iConfig): ChargedRelIsoCut_ (iConfig.getParameter ("ChargedRelIsoCut" )) , ChargedAbsIsoCut_ (iConfig.getParameter ("ChargedAbsIsoCut" )) , IsoConeSize_ (iConfig.getParameter ("IsoConeSize" )) , + MatchingConeSize_ (iConfig.getParameter ("MatchingConeSize" )) , MinTriMuonMass_ (iConfig.getParameter ("MinTriMuonMass" )) , MaxTriMuonMass_ (iConfig.getParameter ("MaxTriMuonMass" )) , TriMuonAbsCharge_ (iConfig.getParameter ("TriMuonAbsCharge" )) , @@ -121,7 +123,7 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS const reco::TrackRef &tk_i = i->track(); for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1 = true; + if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_1 = true; } for (auto j = i+1; j != AllMuCands_end; ++j) { // check that muon_j passes the previous filter @@ -129,16 +131,18 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS const reco::TrackRef &tk_j = j->track(); for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2 = true; + if (reco::deltaR2(tk_j->momentum(), candTrkRef->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; const reco::TrackRef &tk_k = k->track(); for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3 = true; - } + if (reco::deltaR2(tk_k->momentum(), candTrkRef->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 ) || @@ -170,7 +174,7 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS Tau.setPdgId(tauPdgId); Tau.setVertex((Daughters)[0].vertex()); // assign the leading muon vertex as tau vertex - // the three muons must be close in Z + // 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; @@ -230,6 +234,7 @@ HLTTriMuonIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptio 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 ("TriMuonAbsCharge", -1 ); From 32439c58371d6fcce887fb03ca9a9754d0356629 Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Fri, 14 Jul 2017 12:12:36 +0200 Subject: [PATCH 4/7] remove unneeded include from boost lib --- HLTrigger/Muon/interface/HLTTriMuonIsolation.h | 1 - 1 file changed, 1 deletion(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index 944af835283d4..f0f65c65850af 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -3,7 +3,6 @@ #include #include -#include #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/global/EDProducer.h" From a3a2fd1177a6cb9d5ee514b2989fa1dab51393a4 Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Fri, 14 Jul 2017 16:05:28 +0200 Subject: [PATCH 5/7] add option to require collimated muons --- HLTrigger/Muon/interface/HLTTriMuonIsolation.h | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index f0f65c65850af..ddd6785f25ea5 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -53,6 +53,7 @@ class HLTTriMuonIsolation : public edm::global::EDProducer<> { const double MatchingConeSize_; const double MinTriMuonMass_ ; const double MaxTriMuonMass_ ; + const double MaxTriMuonRadius_; const int TriMuonAbsCharge_; const double MaxDZ_ ; const bool EnableRelIso_ ; @@ -75,6 +76,7 @@ HLTTriMuonIsolation::HLTTriMuonIsolation(const edm::ParameterSet& iConfig): 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" )) , @@ -177,7 +179,19 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS 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); } } @@ -236,6 +250,7 @@ HLTTriMuonIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptio 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); From 97140bd7050e9463871ca8f2a2f3e2ec772251c5 Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Mon, 17 Jul 2017 17:11:50 +0200 Subject: [PATCH 6/7] pass by reference, improve loops, stick to CMS' variable naming style --- .../Muon/interface/HLTTriMuonIsolation.h | 102 +++++++++--------- 1 file changed, 50 insertions(+), 52 deletions(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index ddd6785f25ea5..d42e51f34d7c6 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -40,7 +40,7 @@ class HLTTriMuonIsolation : public edm::global::EDProducer<> { 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(); } + 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_ ; @@ -121,69 +121,64 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS for (auto i = AllMuCands->begin(); i != AllMuCands_end; ++i) { // check that muon_i passes the previous filter bool passingPreviousFilter_1 = false; - const reco::TrackRef &tk_i = i->track(); - for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ - reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_1 = true; + 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; ++j) { // check that muon_j passes the previous filter bool passingPreviousFilter_2 = false; - const reco::TrackRef &tk_j = j->track(); - for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ - reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_2 = true; + 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; - const reco::TrackRef &tk_k = k->track(); - for (std::vector::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){ - reco::TrackRef candTrkRef = (*imu)->get(); - if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_3 = true; - } + 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; + reco::CompositeCandidate tau; // sort the muons by pt and add them to the tau - reco::RecoChargedCandidateCollection Daughters; + reco::RecoChargedCandidateCollection daughters; + daughters.reserve(3); - Daughters.push_back(*i); - Daughters.push_back(*j); - Daughters.push_back(*k); + daughters.push_back(*i); + daughters.push_back(*j); + daughters.push_back(*k); - std::sort(Daughters.begin(), Daughters.end(), ptComparer); + 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"); + 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 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 + 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; + 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_) { + for (auto const &idau : daughters){ + if (reco::deltaR2(tau.p4(), idau.p4()) > MaxTriMuonRadius_*MaxTriMuonRadius_) { collimated = false; break; } @@ -192,37 +187,40 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS if (!collimated) continue; // a good tau, at last - Taus->push_back(Tau); + Taus->push_back(tau); } } } // Loop over taus and further select - for (reco::CompositeCandidateCollection::const_iterator itau = Taus->begin(); itau != Taus->end(); ++itau){ - 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; - + 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(); + double sumPt = -itau.pt(); // compute iso sum pT - for (reco::TrackCollection::const_iterator itrk = IsoTracks->begin(); itrk != IsoTracks->end(); ++itrk){ - if (reco::deltaR2(itrk->momentum(), itau->p4()) > IsoConeSize_*IsoConeSize_) continue; - if (std::abs(itrk->vz() - itau->vz()) > MaxDZ_) continue; - sumPt += itrk->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 ((std::max(0., sumPt) > (EnableAbsIso_ * ChargedAbsIsoCut_)) || - (std::max(0., sumPt) > (EnableRelIso_ * ChargedRelIsoCut_ * itau->pt()))) continue; + (std::max(0., sumPt) > (EnableRelIso_ * ChargedRelIsoCut_ * itau.pt()))) continue; - SelectedTaus->push_back(*itau); + SelectedTaus->push_back(itau); } } From dc4fe754b6367cd8f57479adfe69a46d8873107f Mon Sep 17 00:00:00 2001 From: Riccardo Manzoni Date: Mon, 17 Jul 2017 17:25:45 +0200 Subject: [PATCH 7/7] shorten loops when muons run out --- HLTrigger/Muon/interface/HLTTriMuonIsolation.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h index d42e51f34d7c6..32277e6bde769 100644 --- a/HLTrigger/Muon/interface/HLTTriMuonIsolation.h +++ b/HLTrigger/Muon/interface/HLTTriMuonIsolation.h @@ -118,13 +118,13 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS // 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; ++i) { + 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; ++j) { + 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){ @@ -217,8 +217,8 @@ HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventS } // apply the isolation cut - if ((std::max(0., sumPt) > (EnableAbsIso_ * ChargedAbsIsoCut_)) || - (std::max(0., sumPt) > (EnableRelIso_ * ChargedRelIsoCut_ * itau.pt()))) continue; + if ((sumPt > (EnableAbsIso_ * ChargedAbsIsoCut_)) || + (sumPt > (EnableRelIso_ * ChargedRelIsoCut_ * itau.pt()))) continue; SelectedTaus->push_back(itau); }