BPHNano/plugins/MuonTriggerSelector.cc

0001 /////////////////////////////// MuonTriggerSelector //////////////////////////
0002 /// original authors: G Karathanasis (CERN),  G Melachroinos (NKUA)
0003 // filters muons and produces 3 collections: all muons that pass the selection
0004 // triggering muons and transient tracks of triggering muons
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/stream/EDProducer.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Utilities/interface/StreamID.h"
0014 #include "FWCore/Common/interface/TriggerNames.h"
0015
0016 #include "DataFormats/Common/interface/TriggerResults.h"
0017 #include "DataFormats/PatCandidates/interface/Muon.h"
0018 #include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
0019 #include "DataFormats/PatCandidates/interface/PackedTriggerPrescales.h"
0020 #include "DataFormats/PatCandidates/interface/TriggerPath.h"
0021 #include "DataFormats/PatCandidates/interface/TriggerEvent.h"
0022 #include "DataFormats/PatCandidates/interface/TriggerAlgorithm.h"
0023
0024 #include "MagneticField/Engine/interface/MagneticField.h"
0025 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0026
0027 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0028 #include "helper.h"
0029
0030 using namespace std;
0031
0032 constexpr bool debug = false;
0033
0034 class MuonTriggerSelector : public edm::stream::EDProducer<> {
0035 public:
0036   explicit MuonTriggerSelector(const edm::ParameterSet &iConfig);
0037   ~MuonTriggerSelector() override {};
0038
0039 private:
0040   void produce(edm::Event &, const edm::EventSetup &) override;
0041   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0042   edm::EDGetTokenT<std::vector<pat::Muon>> muonSrc_;
0043   edm::EDGetTokenT<edm::TriggerResults> triggerBits_;
0044   edm::EDGetTokenT<std::vector<pat::TriggerObjectStandAlone>> triggerObjects_;
0045   edm::EDGetTokenT<pat::PackedTriggerPrescales> triggerPrescales_;
0046   const StringCutObjectSelector<pat::Muon> muon_selection_;
0047   // for trigger match
0048   const double maxdR_;
0049   // triggers
0050   std::vector<std::string> HLTPaths_;
0051 };
0052
0053 MuonTriggerSelector::MuonTriggerSelector(const edm::ParameterSet &iConfig)
0054     : bFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0055       muonSrc_(consumes<std::vector<pat::Muon>>(iConfig.getParameter<edm::InputTag>("muonCollection"))),
0056       triggerBits_(consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("bits"))),
0057       triggerObjects_(
0058           consumes<std::vector<pat::TriggerObjectStandAlone>>(iConfig.getParameter<edm::InputTag>("objects"))),
0059       triggerPrescales_(consumes<pat::PackedTriggerPrescales>(iConfig.getParameter<edm::InputTag>("prescales"))),
0060       muon_selection_{iConfig.getParameter<std::string>("muonSelection")},
0061       maxdR_(iConfig.getParameter<double>("maxdR_matching")),
0062       HLTPaths_(iConfig.getParameter<std::vector<std::string>>("HLTPaths"))  // multiple paths with comma
0063 {
0064   // outputs
0065   produces<pat::MuonCollection>("AllMuons");
0066   produces<pat::MuonCollection>("SelectedMuons");
0067   produces<TransientTrackCollection>("SelectedTransientMuons");
0068 }
0069
0070 void MuonTriggerSelector::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0071   const auto &bField = iSetup.getData(bFieldToken_);
0072
0073   edm::Handle<edm::TriggerResults> triggerBits;
0074   iEvent.getByToken(triggerBits_, triggerBits);
0075
0076   std::vector<pat::TriggerObjectStandAlone> triggeringMuons;
0077
0078   //taken from https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookMiniAOD2016#Trigger
0079   edm::Handle<std::vector<pat::TriggerObjectStandAlone>> triggerObjects;
0080   iEvent.getByToken(triggerObjects_, triggerObjects);
0081
0082   std::unique_ptr<pat::MuonCollection> allmuons_out(new pat::MuonCollection);
0083   std::unique_ptr<pat::MuonCollection> muons_out(new pat::MuonCollection);
0084   std::unique_ptr<TransientTrackCollection> trans_muons_out(new TransientTrackCollection);
0085
0086   edm::Handle<std::vector<pat::Muon>> muons;
0087   iEvent.getByToken(muonSrc_, muons);
0088
0089   // check for reco muons matched to triggering muons
0090   std::vector<int> muonIsTrigger(muons->size(), 0);
0091   std::vector<float> muonDR(muons->size(), -1.);
0092   std::vector<float> muonDPT(muons->size(), 10000.);
0093   std::vector<int> loose_id(muons->size(), 0);
0094
0095   std::vector<int> matched_reco_flag(muons->size(), -1);
0096   std::vector<int> matched_trg_index(muons->size(), -1);
0097   std::vector<float> matched_dr(muons->size(), -1.);
0098   std::vector<float> matched_dpt(muons->size(), -10000.);
0099   std::vector<std::vector<int>> fires;
0100   std::vector<std::vector<float>> matcher;
0101   std::vector<std::vector<float>> DR;
0102   std::vector<std::vector<float>> DPT;
0103
0104   for (const pat::Muon &muon : *muons) {
0105     if (debug)
0106       std::cout << "Muon Pt=" << muon.pt() << " Eta=" << muon.eta() << " Phi=" << muon.phi() << endl;
0107
0108     std::vector<int> frs(HLTPaths_.size(), 0);  //path fires for each reco muon
0109     std::vector<float> temp_matched_to(HLTPaths_.size(), 1000.);
0110     std::vector<float> temp_DR(HLTPaths_.size(), 1000.);
0111     std::vector<float> temp_DPT(HLTPaths_.size(), 1000.);
0112     int ipath = -1;
0113
0114     for (const std::string &path : HLTPaths_) {
0115       ipath++;
0116       // the following vectors are used in order to find the minimum DR between a reco muon and all the HLT objects that is matched with it so as a reco muon will be matched with only one HLT object every time so as there is a one-to-one correspondance between the two collection. DPt_rel is not used to create this one-to-one correspondance but only to create a few plots, debugging and be sure thateverything is working fine.
0117       std::vector<float> temp_dr(muon.triggerObjectMatches().size(), 1000.);
0118       std::vector<float> temp_dpt(muon.triggerObjectMatches().size(), 1000.);
0119       std::vector<float> temp_pt(muon.triggerObjectMatches().size(), 1000.);
0120       char cstr[(path + "*").size() + 1];
0121       strcpy(cstr, (path + "*").c_str());
0122       //Here we find all the HLT objects from each HLT path each time that are matched with the reco muon.
0123       if (!muon.triggerObjectMatches().empty()) {
0124         for (size_t i = 0; i < muon.triggerObjectMatches().size(); i++) {
0125           if (muon.triggerObjectMatch(i) != nullptr && muon.triggerObjectMatch(i)->hasPathName(cstr, true, true)) {
0126             frs[ipath] = 1;
0127             float dr = TMath::Sqrt(pow(muon.triggerObjectMatch(i)->eta() - muon.eta(), 2.) +
0128                                    pow(muon.triggerObjectMatch(i)->phi() - muon.phi(), 2.));
0129             float dpt = (muon.triggerObjectMatch(i)->pt() - muon.pt()) / muon.triggerObjectMatch(i)->pt();
0130             temp_dr[i] = dr;
0131             temp_dpt[i] = dpt;
0132             temp_pt[i] = muon.triggerObjectMatch(i)->pt();
0133             if (debug)
0134               std::cout << "Path=" << cstr << endl;
0135             if (debug)
0136               std::cout << "HLT  Pt=" << muon.triggerObjectMatch(i)->pt()
0137                         << " Eta=" << muon.triggerObjectMatch(i)->eta() << " Phi=" << muon.triggerObjectMatch(i)->phi()
0138                         << endl;
0139             if (debug)
0140               std::cout << "Muon Pt=" << muon.pt() << " Eta=" << muon.eta() << " Phi=" << muon.phi() << endl;
0141             if (debug)
0142               std::cout << "DR = " << temp_dr[i] << endl;
0143           }
0144         }
0145         // and now we find the real minimum between the reco muon and all its matched HLT objects.
0146         temp_DR[ipath] = *min_element(temp_dr.begin(), temp_dr.end());
0147         int position = std::min_element(temp_dr.begin(), temp_dr.end()) - temp_dr.begin();
0148         temp_DPT[ipath] = temp_dpt[position];
0149         temp_matched_to[ipath] = temp_pt[position];
0150       }
0151     }
0152     //and now since we have found the minimum DR we save a few variables for plots
0153     fires.push_back(frs);  //This is used in order to see if a reco muon fired a Trigger (1) or not (0).
0154     matcher.push_back(
0155         temp_matched_to);  //This is used in order to see if a reco muon is matched with a HLT object. PT of the reco muon is saved in this vector.
0156     DR.push_back(temp_DR);
0157     DPT.push_back(temp_DPT);
0158   }
0159
0160   //now, check for different reco muons that are matched to the same HLTObject.
0161   for (unsigned int path = 0; path < HLTPaths_.size(); path++) {
0162     for (unsigned int iMuo = 0; iMuo < muons->size(); iMuo++) {
0163       for (unsigned int im = (iMuo + 1); im < muons->size(); im++) {
0164         if (matcher[iMuo][path] != 1000. && matcher[iMuo][path] == matcher[im][path]) {
0165           if (DR[iMuo][path] < DR[im][path]) {  //Keep the one that has the minimum DR with the HLT object
0166             fires[im][path] = 0;
0167             matcher[im][path] = 1000.;
0168             DR[im][path] = 1000.;
0169             DPT[im][path] = 1000.;
0170           } else {
0171             fires[iMuo][path] = 0;
0172             matcher[iMuo][path] = 1000.;
0173             DR[iMuo][path] = 1000.;
0174             DPT[iMuo][path] = 1000.;
0175           }
0176         }
0177       }
0178       if (matcher[iMuo][path] != 1000.) {
0179         muonIsTrigger[iMuo] = 1;
0180         muonDR[iMuo] = DR[iMuo][path];
0181         muonDPT[iMuo] = DPT[iMuo][path];
0182       }
0183     }
0184   }
0185
0186   // Save the reco muon in both collections
0187   for (const pat::Muon &muon : *muons) {
0188     unsigned int iMuo(&muon - &(muons->at(0)));
0189     if (!muon_selection_(muon))
0190       continue;  // selection cuts
0191     const reco::TransientTrack muonTT((*(muon.bestTrack())), &bField);
0192     if (!muonTT.isValid())
0193       continue;
0194
0195     // Save in AllMuons
0196     // Make a copy of the muon and add the trigger flag for filtering Allmuons
0197     pat::Muon muonCopy(muon);
0198     muonCopy.addUserInt("isTriggering", muonIsTrigger[iMuo]);  // 1 if triggered, 0 otherwise
0199     allmuons_out->push_back(muonCopy);
0200
0201     // Save in selectedMuons for triggering muons
0202     if (muonIsTrigger[iMuo] != 1)
0203       continue;
0204
0205     muons_out->emplace_back(muon);
0206     muons_out->back().addUserInt("isTriggering", muonIsTrigger[iMuo]);
0207     muons_out->back().addUserFloat("trgDR", muonDR[iMuo]);
0208     muons_out->back().addUserFloat("trgDPT", muonDPT[iMuo]);
0209     muons_out->back().addUserInt("looseId", loose_id[iMuo]);
0210
0211     for (unsigned int i = 0; i < HLTPaths_.size(); i++)
0212       muons_out->back().addUserInt(HLTPaths_[i], fires[iMuo][i]);
0213
0214     trans_muons_out->emplace_back(muonTT);
0215   }
0216
0217   iEvent.put(std::move(allmuons_out), "AllMuons");
0218   iEvent.put(std::move(muons_out), "SelectedMuons");
0219   iEvent.put(std::move(trans_muons_out), "SelectedTransientMuons");
0220 }
0221
0222 DEFINE_FWK_MODULE(MuonTriggerSelector);