L1TNtuples/plugins/L1Muon2RecoTreeProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1Trigger/L1TNtuples
0004 // Class:      L1JetRecoTreeProducer
0005 //
0006 /**\class L1JetRecoTreeProducer L1JetRecoTreeProducer.cc L1Trigger/L1TNtuples/src/L1JetRecoTreeProducer.cc
0007
0008    Description: Produces tree containing reco quantities
0009
0010
0011 */
0012
0013 // system include files
0014 #include <memory>
0015
0016 // framework
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Common/interface/TriggerNames.h"
0026
0027 // cond formats
0028 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0029
0030 // data formats
0031 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0032 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0033 #include "DataFormats/JetReco/interface/JetID.h"
0034 #include "DataFormats/METReco/interface/PFMET.h"
0035 #include "DataFormats/METReco/interface/PFMETCollection.h"
0036
0037 // RECO TRIGGER MATCHING:
0038 #include "DataFormats/Math/interface/deltaR.h"
0039 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0040 #include "DataFormats/Common/interface/TriggerResults.h"
0041 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0042 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0043 #include "TString.h"
0044 #include "TRegexp.h"
0045 #include <utility>
0046
0047 //muons
0048 #include "DataFormats/MuonReco/interface/Muon.h"
0049 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0052 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0053 #include "DataFormats/GeometrySurface/interface/Plane.h"
0054 #include "DataFormats/MuonReco/interface/MuonEnergy.h"
0055 #include "DataFormats/MuonReco/interface/MuonTime.h"
0056 #include "CondFormats/AlignmentRecord/interface/TrackerSurfaceDeformationRcd.h"
0057
0058 //taus
0059 #include "DataFormats/TauReco/interface/PFTau.h"
0060 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0061 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
0062
0063 // ROOT output stuff
0064 #include "FWCore/ServiceRegistry/interface/Service.h"
0065 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0066 #include "TH1.h"
0067 #include "TTree.h"
0068 #include "TF1.h"
0069
0070 //local  data formats
0071 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMuon2.h"
0072 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMet.h"
0073
0074 //vertices
0075 #include "DataFormats/VertexReco/interface/Vertex.h"
0076 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0077 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoVertexDataFormat.h"
0078
0079 using namespace std;
0080
0081 //
0082 // class declaration
0083 //
0084
0085 class L1Muon2RecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources, edm::one::WatchRuns> {
0086 public:
0087   explicit L1Muon2RecoTreeProducer(const edm::ParameterSet &);
0088   ~L1Muon2RecoTreeProducer() override;
0089
0090 private:
0091   void beginJob() override;
0092   void analyze(const edm::Event &, const edm::EventSetup &) override;
0093   void beginRun(const edm::Run &, const edm::EventSetup &) override;
0094   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0095   void endJob() override;
0096
0097 public:
0098   L1Analysis::L1AnalysisRecoMuon2 *muon;
0099
0100   L1Analysis::L1AnalysisRecoMuon2DataFormat *muon_data;
0101
0102   double match_trigger(std::vector<int> &trigIndices,
0103                        const trigger::TriggerObjectCollection &trigObjs,
0104                        edm::Handle<trigger::TriggerEvent> &triggerEvent,
0105                        const reco::Muon &mu);
0106   void empty_hlt();
0107
0108 private:
0109   // output file
0110   edm::Service<TFileService> fs_;
0111
0112   // tree
0113   TTree *tree_;
0114
0115   // EDM input tags
0116   edm::EDGetTokenT<reco::MuonCollection> MuonToken_;
0117   edm::EDGetTokenT<edm::TriggerResults> TriggerResultsToken_;
0118   edm::EDGetTokenT<trigger::TriggerEvent> triggerSummaryLabelToken_;
0119   edm::EDGetTokenT<reco::VertexCollection> VtxToken_;
0120   edm::EDGetTokenT<reco::PFMETCollection> metToken_;
0121
0122   // bool triggerMatching_;
0123   // edm::InputTag triggerSummaryLabel_;
0124   // std::string triggerProcessLabel_;
0125   // std::vector<std::string> isoTriggerNames_;
0126   // std::vector<std::string> triggerNames_;
0127
0128   edm::Handle<edm::TriggerResults> isoTriggerToken_;
0129   edm::Handle<std::vector<std::string>> isoTriggerNamesToken_;
0130   // edm::Handle<edm::TriggerResults>     TriggerToken_;
0131
0132   HLTConfigProvider hltConfig_;
0133
0134   // debug stuff
0135   unsigned int maxMuon_;
0136   double triggerMaxDeltaR_;
0137   bool triggerMatching_;
0138   std::string triggerProcessLabel_;
0139   std::vector<std::string> isoTriggerNames_;
0140   std::vector<std::string> triggerNames_;
0141   std::vector<int> isoTriggerIndices_;
0142   std::vector<int> triggerIndices_;
0143 };
0144
0145 L1Muon2RecoTreeProducer::L1Muon2RecoTreeProducer(const edm::ParameterSet &iConfig) {
0146   usesResource(TFileService::kSharedResource);
0147
0148   maxMuon_ = iConfig.getParameter<unsigned int>("maxMuon");
0149   isoTriggerNames_ = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
0150   // isoTriggerToken_         = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
0151   // TriggerToken_         = iConfig.getParameter<std::vector<std::string>>("triggerNames");
0152   MuonToken_ = consumes<reco::MuonCollection>(iConfig.getUntrackedParameter("MuonToken", edm::InputTag("muons")));
0153   VtxToken_ = consumes<reco::VertexCollection>(
0154       iConfig.getUntrackedParameter("VertexToken", edm::InputTag("offlinePrimaryVertices")));
0155
0156   TriggerResultsToken_ = consumes<edm::TriggerResults>(
0157       iConfig.getUntrackedParameter("TriggerResultsToken", edm::InputTag("TriggerResults", "", "HLT")));
0158   triggerSummaryLabelToken_ = consumes<trigger::TriggerEvent>(
0159       iConfig.getUntrackedParameter("triggerSummaryLabelToken", edm::InputTag("hltTriggerSummaryAOD", "", "HLT")));
0160   // TriggerResultsToken_      = consumes<edm::TriggerResults>(iConfig.getUntrackedParameter("TriggerResultsToken",edm::InputTag("triggerSummaryLabel")));
0161   // triggerSummaryLabelToken_ = consumes<trigger::TriggerEvent>(iConfig.getUntrackedParameter("triggerSummaryLabelToken",edm::InputTag("triggerSummaryLabel")));
0162   //iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
0163   metToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("metToken", edm::InputTag("pfMet")));
0164
0165   muon = new L1Analysis::L1AnalysisRecoMuon2(iConfig, consumesCollector());
0166   muon_data = muon->getData();
0167
0168   tree_ = fs_->make<TTree>("Muon2RecoTree", "Muon2RecoTree");
0169   tree_->Branch("Muon", "L1Analysis::L1AnalysisRecoMuon2DataFormat", &muon_data, 32000, 3);
0170
0171   triggerMaxDeltaR_ = iConfig.getParameter<double>("triggerMaxDeltaR");
0172   triggerMatching_ = iConfig.getUntrackedParameter<bool>("triggerMatching");
0173   triggerProcessLabel_ = iConfig.getUntrackedParameter<std::string>("triggerProcessLabel");
0174   isoTriggerNames_ = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
0175   triggerNames_ = iConfig.getParameter<std::vector<std::string>>("triggerNames");
0176
0177   // triggerMatching_     = iConfig.getUntrackedParameter<bool>("triggerMatching");
0178   // _ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
0179   // triggerProcessLabel_ = iConfig.getUntrackedParameter<std::string>("triggerProcessLabel");
0180   // triggerNames_        = iConfig.getParameter<std::vector<std::string> > ("triggerNames");
0181   // isoTriggerNames_     = iConfig.getParameter<std::vector<std::string> > ("isoTriggerNames");
0182   // triggerMaxDeltaR_    = iConfig.getParameter<double> ("triggerMaxDeltaR");
0183 }
0184
0185 L1Muon2RecoTreeProducer::~L1Muon2RecoTreeProducer() {
0186   // do anything here that needs to be done at desctruction time
0187   // (e.g. close files, deallocate resources etc.)
0188 }
0189
0190 //
0191 // member functions
0192 //
0193 // ------------ method called to for each event  ------------
0194 void L1Muon2RecoTreeProducer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0195   muon->init(iSetup);
0196
0197   muon->Reset();
0198   edm::Handle<reco::MuonCollection> recoMuons;
0199   iEvent.getByToken(MuonToken_, recoMuons);
0200
0201   edm::Handle<reco::VertexCollection> vertices;
0202   iEvent.getByToken(VtxToken_, vertices);
0203
0204   edm::Handle<edm::TriggerResults> triggerResults;
0205   iEvent.getByToken(TriggerResultsToken_, triggerResults);
0206
0207   edm::Handle<trigger::TriggerEvent> triggerSummaryLabel_;
0208   iEvent.getByToken(triggerSummaryLabelToken_, triggerSummaryLabel_);
0209
0210   edm::Handle<reco::PFMETCollection> metLabel_;
0211   iEvent.getByToken(metToken_, metLabel_);
0212
0213   int counter_met = 0;
0214
0215   double METx = 0.;
0216   double METy = 0.;
0217
0218   for (reco::PFMETCollection::const_iterator imet = metLabel_->begin();
0219        imet != metLabel_->end() && (unsigned)counter_met < 1;
0220        imet++) {
0221     METx = imet->px();
0222     METy = imet->py();
0223   }
0224
0225   if (recoMuons.isValid()) {
0226     muon->SetMuon(iEvent, iSetup, recoMuons, vertices, METx, METy, maxMuon_);
0227   } else {
0228   }
0229
0230   int counter_mu = 0;
0231   for (reco::MuonCollection::const_iterator imu = recoMuons->begin();
0232        imu != recoMuons->end() && (unsigned)counter_mu < maxMuon_;
0233        imu++) {
0234     //---------------------------------------------------------------------
0235     // TRIGGER MATCHING:
0236     // if specified the reconstructed muons are matched to a trigger
0237     //---------------------------------------------------------------------
0238     if (triggerMatching_) {
0239       double isoMatchDeltaR = 9999.;
0240       double matchDeltaR = 9999.;
0241       int hasIsoTriggered = 0;
0242       int hasTriggered = 0;
0243
0244       int passesSingleMuonFlag = 0;
0245
0246       // first check if the trigger results are valid:
0247       if (triggerResults.isValid()) {
0248         if (triggerSummaryLabel_.isValid()) {
0249           const edm::TriggerNames &trigNames = iEvent.triggerNames(*triggerResults);
0250           // for(UInt_t iPath = 0 ; iPath < trigNames.size() ; ++iPath)
0251           //   {
0252           //     cout<<iPath<<": "<<trigNames.triggerName(iPath)<<endl;
0253           //   }
0254
0255           for (UInt_t iPath = 0; iPath < isoTriggerNames_.size(); ++iPath) {
0256             if (passesSingleMuonFlag == 1)
0257               continue;
0258             std::string pathName = isoTriggerNames_.at(iPath);
0259
0260             bool passTrig = false;
0261             //cout<<"testing pathName                         = "<<pathName<<endl;
0262             //cout<<"trigNames.triggerIndex(pathName) = "<<trigNames.triggerIndex(pathName)<<endl;
0263             // cout<<"trigNames.triggerIndex(pathName)<trigNames.size()= "<<(trigNames.triggerIndex(pathName)<trigNames.size())<<endl;
0264
0265             if (trigNames.triggerIndex(pathName) < trigNames.size())
0266               passTrig = triggerResults->accept(trigNames.triggerIndex(pathName));
0267             // cout<<"pass = "<<passTrig<<endl;
0268             if (passTrig)
0269               passesSingleMuonFlag = 1;
0270           }
0271
0272           if (triggerSummaryLabel_.isValid()) {
0273             // get trigger objects:
0274             const trigger::TriggerObjectCollection triggerObjects = triggerSummaryLabel_->getObjects();
0275
0276             matchDeltaR = match_trigger(triggerIndices_, triggerObjects, triggerSummaryLabel_, (*imu));
0277             if (matchDeltaR < triggerMaxDeltaR_)
0278               hasTriggered = 1;
0279
0280             // cout<<"isoTriggerIndices_.size() = "<<isoTriggerIndices_.size()<<endl;
0281             // cout<<"triggerObjects.size() = "<<triggerObjects.size()<<endl;
0282             // cout<<"imu->eta() = "<<imu->eta()<<endl;
0283
0284             isoMatchDeltaR = match_trigger(isoTriggerIndices_, triggerObjects, triggerSummaryLabel_, (*imu));
0285             if (isoMatchDeltaR < triggerMaxDeltaR_)
0286               hasIsoTriggered = 1;
0287
0288           }  // end if (triggerEvent.isValid())
0289
0290         }  // end if (triggerResults.isValid())
0291
0292       }  // end if (triggerResults.isValid())
0293
0294       // fill trigger matching variables:
0295       muon_data->hlt_isomu.push_back(hasIsoTriggered);
0296       muon_data->hlt_mu.push_back(hasTriggered);
0297       muon_data->hlt_isoDeltaR.push_back(isoMatchDeltaR);
0298       muon_data->hlt_deltaR.push_back(matchDeltaR);
0299       muon_data->passesSingleMuon.push_back(passesSingleMuonFlag);
0300
0301     } else {
0302       empty_hlt();
0303     }  // end if (triggerMatching_)
0304   }
0305
0306   tree_->Fill();
0307 }
0308
0309 // ------------ method called once each job just before starting event loop  ------------
0310 void L1Muon2RecoTreeProducer::beginJob(void) {}
0311
0312 // ------------ method called once each job just after ending the event loop  ------------
0313 void L1Muon2RecoTreeProducer::endJob() {}
0314
0315 void L1Muon2RecoTreeProducer::empty_hlt() {
0316   muon_data->hlt_isomu.push_back(-9999);
0317   muon_data->hlt_mu.push_back(-9999);
0318   muon_data->hlt_isoDeltaR.push_back(-9999);
0319   muon_data->hlt_deltaR.push_back(-9999);
0320 }
0321
0322 double L1Muon2RecoTreeProducer::match_trigger(std::vector<int> &trigIndices,
0323                                               const trigger::TriggerObjectCollection &trigObjs,
0324                                               edm::Handle<trigger::TriggerEvent> &triggerEvent,
0325                                               const reco::Muon &mu) {
0326   double matchDeltaR = 9999;
0327
0328   for (size_t iTrigIndex = 0; iTrigIndex < trigIndices.size(); ++iTrigIndex) {
0329     int triggerIndex = trigIndices[iTrigIndex];
0330     const std::vector<std::string> moduleLabels(hltConfig_.moduleLabels(triggerIndex));
0331     // find index of the last module:
0332     const unsigned moduleIndex = hltConfig_.size(triggerIndex) - 2;
0333     // find index of HLT trigger name:
0334     const unsigned hltFilterIndex =
0335         triggerEvent->filterIndex(edm::InputTag(moduleLabels[moduleIndex], "", triggerProcessLabel_));
0336
0337     if (hltFilterIndex < triggerEvent->sizeFilters()) {
0338       const trigger::Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
0339       const trigger::Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));
0340
0341       const unsigned nTriggers = triggerVids.size();
0342       for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
0343         // loop over all trigger objects:
0344         const trigger::TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];
0345
0346         double dRtmp = deltaR(mu, trigObject);
0347
0348         if (dRtmp < matchDeltaR) {
0349           matchDeltaR = dRtmp;
0350         }
0351
0352       }  // loop over different trigger objects
0353     }  // if trigger is in event (should apply hltFilter with used trigger...)
0354   }  // loop over muon candidates
0355
0356   return matchDeltaR;
0357 }
0358
0359 void L1Muon2RecoTreeProducer::beginRun(const edm::Run &run, const edm::EventSetup &eventSetup) {
0360   // Prepare for trigger matching for each new run:
0361   // Look up triggetIndices in the HLT config for the different paths
0362   if (triggerMatching_) {
0363     bool changed = true;
0364     if (!hltConfig_.init(run, eventSetup, triggerProcessLabel_, changed)) {
0365       // if you can't initialize hlt configuration, crash!
0366       std::cout << "Error: didn't find process" << triggerProcessLabel_ << std::endl;
0367       assert(false);
0368     }
0369
0370     bool enableWildcard = true;
0371     for (size_t iTrig = 0; iTrig < triggerNames_.size(); ++iTrig) {
0372       // prepare for regular expression (with wildcards) functionality:
0373       TString tNameTmp = TString(triggerNames_[iTrig]);
0374       TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
0375       int tIndex = -1;
0376       // find the trigger index:
0377       for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
0378         // use TString since it provides reg exp functionality:
0379         TString tmpName = TString(hltConfig_.triggerName(ipath));
0380         if (tmpName.Contains(tNamePattern)) {
0381           tIndex = int(ipath);
0382           triggerIndices_.push_back(tIndex);
0383         }
0384       }
0385       if (tIndex < 0) {  // if can't find trigger path at all, give warning:
0386         std::cout << "Warning: Could not find trigger" << triggerNames_[iTrig] << std::endl;
0387         //assert(false);
0388       }
0389     }  // end for triggerNames
0390     for (size_t iTrig = 0; iTrig < isoTriggerNames_.size(); ++iTrig) {
0391       // prepare for regular expression functionality:
0392       TString tNameTmp = TString(isoTriggerNames_[iTrig]);
0393       TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
0394       int tIndex = -1;
0395       // find the trigger index:
0396       for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
0397         // use TString since it provides reg exp functionality:
0398         TString tmpName = TString(hltConfig_.triggerName(ipath));
0399         if (tmpName.Contains(tNamePattern)) {
0400           tIndex = int(ipath);
0401           isoTriggerIndices_.push_back(tIndex);
0402         }
0403       }
0404       if (tIndex < 0) {  // if can't find trigger path at all, give warning:
0405         std::cout << "Warning: Could not find trigger" << isoTriggerNames_[iTrig] << std::endl;
0406         //assert(false);
0407       }
0408     }  // end for isoTriggerNames
0409   }  // end if (triggerMatching_)
0410 }
0411
0412 //define this as a plug-in
0413 DEFINE_FWK_MODULE(L1Muon2RecoTreeProducer);