HcalAlCaRecoProducers/plugins/AlCaLowPUHBHEMuonFilter.cc

0001 // system include files
0002 #include <atomic>
0003 #include <memory>
0004 #include <cmath>
0005 #include <iostream>
0006 #include <sstream>
0007 #include <fstream>
0008
0009 // user include files
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/stream/EDFilter.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015 #include "FWCore/Framework/interface/LuminosityBlock.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/Common/interface/TriggerNames.h"
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/Common/interface/Ref.h"
0023 #include "DataFormats/Common/interface/TriggerResults.h"
0024 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0025 #include "DataFormats/MuonReco/interface/Muon.h"
0026 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0027
0028 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0029 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0030
0031 #include "MagneticField/Engine/interface/MagneticField.h"
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0035
0036 //#define EDM_ML_DEBUG
0037 //
0038 // class declaration
0039 //
0040
0041 namespace alCaLowPUHBHEMuonFilter {
0042   struct Counters {
0043     Counters() : nAll_(0), nGood_(0) {}
0044     mutable std::atomic<unsigned int> nAll_, nGood_;
0045   };
0046 }  // namespace alCaLowPUHBHEMuonFilter
0047
0048 class AlCaLowPUHBHEMuonFilter : public edm::stream::EDFilter<edm::GlobalCache<alCaLowPUHBHEMuonFilter::Counters> > {
0049 public:
0050   explicit AlCaLowPUHBHEMuonFilter(edm::ParameterSet const&, const alCaLowPUHBHEMuonFilter::Counters* count);
0051   ~AlCaLowPUHBHEMuonFilter() override;
0052
0053   static std::unique_ptr<alCaLowPUHBHEMuonFilter::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0054     return std::make_unique<alCaLowPUHBHEMuonFilter::Counters>();
0055   }
0056
0057   bool filter(edm::Event&, edm::EventSetup const&) override;
0058   void endStream() override;
0059   static void globalEndJob(const alCaLowPUHBHEMuonFilter::Counters* counters);
0060   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0061
0062 private:
0063   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0064   void endRun(edm::Run const&, edm::EventSetup const&) override;
0065
0066   // ----------member data ---------------------------
0067   HLTConfigProvider hltConfig_;
0068   std::vector<std::string> trigNames_, HLTNames_;
0069   std::string processName_;
0070   bool pfCut_;
0071   double trackIsoCut_, caloIsoCut_, pfIsoCut_, minimumMuonpT_, minimumMuoneta_;
0072   int preScale_;
0073   unsigned int nRun_, nAll_, nGood_;
0074   edm::InputTag triggerResults_, labelMuon_;
0075   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt;
0076   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0077   edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0078   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0079   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0080 };
0081
0082 //
0083 // constants, enums and typedefs
0084 //
0085
0086 //
0087 // static data member definitions
0088 //
0089
0090 //
0091 // constructors and destructor
0092 //
0093 AlCaLowPUHBHEMuonFilter::AlCaLowPUHBHEMuonFilter(edm::ParameterSet const& iConfig,
0094                                                  const alCaLowPUHBHEMuonFilter::Counters* count)
0095     : nRun_(0), nAll_(0), nGood_(0) {
0096   //now do what ever initialization is needed
0097   trigNames_ = iConfig.getParameter<std::vector<std::string> >("triggers");
0098   processName_ = iConfig.getParameter<std::string>("processName");
0099   triggerResults_ = iConfig.getParameter<edm::InputTag>("triggerResultLabel");
0100   labelMuon_ = iConfig.getParameter<edm::InputTag>("muonLabel");
0101   pfIsoCut_ = iConfig.getParameter<double>("pfIsolationCut");
0102   minimumMuonpT_ = iConfig.getParameter<double>("minimumMuonpT");
0103   minimumMuoneta_ = iConfig.getParameter<double>("minimumMuoneta");
0104
0105   // define tokens for access
0106   tok_trigRes_ = consumes<edm::TriggerResults>(triggerResults_);
0107   tok_Muon_ = consumes<reco::MuonCollection>(labelMuon_);
0108
0109   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0110   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0111
0112   edm::LogVerbatim("LowPUHBHEMuon") << "Parameters read from config file \n"
0113                                     << "Process " << processName_ << "  PF Isolation Cuts " << pfIsoCut_
0114                                     << " minimum Muon pT cut " << minimumMuonpT_ << " minimum Muon eta cut "
0115                                     << minimumMuoneta_;
0116   for (unsigned int k = 0; k < trigNames_.size(); ++k)
0117     edm::LogVerbatim("LowPUHBHEMuon") << "Trigger[" << k << "] " << trigNames_[k];
0118 }  // AlCaLowPUHBHEMuonFilter::AlCaLowPUHBHEMuonFilter  constructor
0119
0120 AlCaLowPUHBHEMuonFilter::~AlCaLowPUHBHEMuonFilter() {}
0121
0122 //
0123 // member functions
0124 //
0125
0126 // ------------ method called on each new Event  ------------
0127 bool AlCaLowPUHBHEMuonFilter::filter(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0128   bool accept(false);
0129   ++nAll_;
0130 #ifdef EDM_ML_DEBUG
0131   edm::LogVerbatim("LowPUHBHEMuon") << "AlCaLowPUHBHEMuonFilter::Run " << iEvent.id().run() << " Event "
0132                                     << iEvent.id().event() << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0133                                     << iEvent.bunchCrossing();
0134 #endif
0135   //Step1: Find if the event passes one of the chosen triggers
0136   /////////////////////////////TriggerResults
0137   auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
0138   if (triggerResults.isValid()) {
0139     bool ok(false);
0140     std::vector<std::string> modules;
0141     const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0142     const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0143     for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0144       int hlt = triggerResults->accept(iHLT);
0145 #ifdef EDM_ML_DEBUG
0146       edm::LogVerbatim("LowPUHBHEMuonX") << "trigger names: " << iHLT << " " << triggerNames_[iHLT];
0147 #endif
0148       for (auto const& trigName : trigNames_) {
0149         if (triggerNames_[iHLT].find(trigName) != std::string::npos) {
0150 #ifdef EDM_ML_DEBUG
0151           edm::LogVerbatim("LowPUHBHEMuonX") << "find trigger names: " << trigName;
0152 #endif
0153           if (hlt > 0)
0154             ok = true;
0155 #ifdef EDM_ML_DEBUG
0156           edm::LogVerbatim("LowPUHBHEMuon")
0157               << "AlCaLowPUHBHEMuonFilter::Trigger " << triggerNames_[iHLT] << " Flag " << hlt << ":" << ok;
0158 #endif
0159         }
0160       }
0161     }
0162     if (ok) {
0163       //Step2: Get geometry/B-field information
0164       const MagneticField* bField = &(iSetup.getData(tok_magField_));
0165       const CaloGeometry* geo = &(iSetup.getData(tok_geom_));
0166
0167       // Relevant blocks from iEvent
0168       auto muonHandle = iEvent.getHandle(tok_Muon_);
0169 #ifdef EDM_ML_DEBUG
0170       edm::LogVerbatim("LowPUHBHEMuon") << "AlCaLowPUHBHEMuonFilter::Muon Handle " << muonHandle.isValid();
0171 #endif
0172       if (muonHandle.isValid()) {
0173         for (reco::MuonCollection::const_iterator RecMuon = muonHandle->begin(); RecMuon != muonHandle->end();
0174              ++RecMuon) {
0175 #ifdef EDM_ML_DEBUG
0176           edm::LogVerbatim("LowPUHBHEMuon")
0177               << "AlCaLowPUHBHEMuonFilter::Muon:Track " << RecMuon->track().isNonnull() << " innerTrack "
0178               << RecMuon->innerTrack().isNonnull() << " outerTrack " << RecMuon->outerTrack().isNonnull()
0179               << " globalTrack " << RecMuon->globalTrack().isNonnull();
0180 #endif
0181           if ((RecMuon->pt() < minimumMuonpT_) || std::abs(RecMuon->eta()) < minimumMuoneta_)
0182             continue;
0183           if ((RecMuon->track().isNonnull()) && (RecMuon->innerTrack().isNonnull()) &&
0184               (RecMuon->outerTrack().isNonnull()) && (RecMuon->globalTrack().isNonnull())) {
0185             const reco::Track* pTrack = (RecMuon->innerTrack()).get();
0186             spr::propagatedTrackID trackID = spr::propagateCALO(pTrack, geo, bField, false);
0187 #ifdef EDM_ML_DEBUG
0188             edm::LogVerbatim("LowPUHBHEMuon")
0189                 << "AlCaLowPUHBHEMuonFilter::Propagate: ECAL " << trackID.okECAL << " to HCAL " << trackID.okHCAL;
0190 #endif
0191             double isolR04 =
0192                 ((RecMuon->pfIsolationR04().sumChargedHadronPt +
0193                   std::max(0.,
0194                            RecMuon->pfIsolationR04().sumNeutralHadronEt + RecMuon->pfIsolationR04().sumPhotonEt -
0195                                (0.5 * RecMuon->pfIsolationR04().sumPUPt))) /
0196                  RecMuon->pt());
0197             bool isoCut = (isolR04 < pfIsoCut_);
0198             if ((trackID.okECAL) && (trackID.okHCAL) && isoCut) {
0199               accept = true;
0200               break;
0201             }
0202           }
0203         }
0204       }
0205     }
0206   }
0207   // Step 4:  Return the acceptance flag
0208   if (accept) {
0209     ++nGood_;
0210   }
0211   return accept;
0212
0213 }  // AlCaLowPUHBHEMuonFilter::filter
0214
0215 // ------------ method called once each job just after ending the event loop  ------------
0216 void AlCaLowPUHBHEMuonFilter::endStream() {
0217   globalCache()->nAll_ += nAll_;
0218   globalCache()->nGood_ += nGood_;
0219 }
0220
0221 void AlCaLowPUHBHEMuonFilter::globalEndJob(const alCaLowPUHBHEMuonFilter::Counters* count) {
0222   edm::LogVerbatim("LowPUHBHEMuon") << "Selects " << count->nGood_ << " good events out of " << count->nAll_
0223                                     << " total # of events";
0224 }
0225
0226 // ------------ method called when starting to processes a run  ------------
0227 void AlCaLowPUHBHEMuonFilter::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0228   bool changed(false);
0229   bool flag = hltConfig_.init(iRun, iSetup, processName_, changed);
0230   edm::LogVerbatim("LowPUHBHEMuon") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init " << flag;
0231 }
0232
0233 // ------------ method called when ending the processing of a run  ------------
0234 void AlCaLowPUHBHEMuonFilter::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0235   edm::LogVerbatim("LowPUHBHEMuon") << "endRun[" << nRun_ << "] " << iRun.run();
0236   nRun_++;
0237 }
0238
0239 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0240 void AlCaLowPUHBHEMuonFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0241   edm::ParameterSetDescription desc;
0242   std::vector<std::string> triggers = {"HLT_L1DoubleMu", "HLT_L1SingleMu"};
0243   desc.add<std::string>("processName", "HLT");
0244   desc.add<edm::InputTag>("triggerResultLabel", edm::InputTag("TriggerResults", "", "HLT"));
0245   desc.add<edm::InputTag>("muonLabel", edm::InputTag("muons"));
0246   desc.add<double>("minimumMuonpT", 10.0);
0247   desc.add<double>("minimumMuoneta", 1.305);
0248   desc.add<std::vector<std::string> >("triggers", triggers);
0249   desc.add<double>("pfIsolationCut", 0.15);
0250   descriptions.add("alcaLowPUHBHEMuonFilter", desc);
0251 }
0252
0253 //define this as a plug-in
0254 #include "FWCore/Framework/interface/MakerMacros.h"
0255 DEFINE_FWK_MODULE(AlCaLowPUHBHEMuonFilter);