Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:09:37

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Common/interface/TriggerNames.h"
0003 #include "FWCore/Utilities/interface/RegexMatch.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/Common/interface/TriggerResults.h"
0006 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/MuonReco/interface/Muon.h"
0009 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0010 #include "DataFormats/Math/interface/deltaR.h"
0011 
0012 #include <TLorentzVector.h>
0013 
0014 #include <memory>
0015 
0016 #include "DQMOffline/Lumi/plugins/ZCounting.h"
0017 
0018 //
0019 // -------------------------------------- Constructor --------------------------------------------
0020 //
0021 ZCounting::ZCounting(const edm::ParameterSet& iConfig)
0022     : triggerResultsInputTag_(iConfig.getParameter<edm::InputTag>("TriggerResults")),
0023       fPVName_token(consumes<reco::VertexCollection>(
0024           iConfig.getUntrackedParameter<std::string>("edmPVName", "offlinePrimaryVertices"))),
0025       fMuonName_token(consumes<reco::MuonCollection>(iConfig.getUntrackedParameter<std::string>("edmName", "muons"))),
0026       fStandaloneRegName_token(consumes<reco::TrackCollection>(
0027           iConfig.getUntrackedParameter<std::string>("StandaloneReg", "standAloneMuons"))),
0028       fStandaloneUpdName_token(consumes<reco::TrackCollection>(
0029           iConfig.getUntrackedParameter<std::string>("StandaloneUpd", "standAloneMuons:UpdatedAtVtx"))),
0030       fTrackName_token(
0031           consumes<reco::TrackCollection>(iConfig.getUntrackedParameter<std::string>("edmTrackName", "generalTracks"))),
0032 
0033       PtCutL1_(iConfig.getUntrackedParameter<double>("PtCutL1")),
0034       PtCutL2_(iConfig.getUntrackedParameter<double>("PtCutL2")),
0035       EtaCutL1_(iConfig.getUntrackedParameter<double>("EtaCutL1")),
0036       EtaCutL2_(iConfig.getUntrackedParameter<double>("EtaCutL2")),
0037 
0038       MassBin_(iConfig.getUntrackedParameter<int>("MassBin")),
0039       MassMin_(iConfig.getUntrackedParameter<double>("MassMin")),
0040       MassMax_(iConfig.getUntrackedParameter<double>("MassMax")),
0041 
0042       LumiBin_(iConfig.getUntrackedParameter<int>("LumiBin")),
0043       LumiMin_(iConfig.getUntrackedParameter<double>("LumiMin")),
0044       LumiMax_(iConfig.getUntrackedParameter<double>("LumiMax")),
0045 
0046       PVBin_(iConfig.getUntrackedParameter<int>("PVBin")),
0047       PVMin_(iConfig.getUntrackedParameter<double>("PVMin")),
0048       PVMax_(iConfig.getUntrackedParameter<double>("PVMax")),
0049 
0050       VtxNTracksFitCut_(iConfig.getUntrackedParameter<double>("VtxNTracksFitMin")),
0051       VtxNdofCut_(iConfig.getUntrackedParameter<double>("VtxNdofMin")),
0052       VtxAbsZCut_(iConfig.getUntrackedParameter<double>("VtxAbsZMax")),
0053       VtxRhoCut_(iConfig.getUntrackedParameter<double>("VtxRhoMax")),
0054 
0055       IDTypestr_(iConfig.getUntrackedParameter<std::string>("IDType")),
0056       IsoTypestr_(iConfig.getUntrackedParameter<std::string>("IsoType")),
0057       IsoCut_(iConfig.getUntrackedParameter<double>("IsoCut")) {
0058   edm::LogInfo("ZCounting") << "Constructor  ZCounting::ZCounting " << std::endl;
0059 
0060   // Trigger settings
0061   triggers = new TriggerTools();
0062   triggers->setTriggerResultsToken(consumes<edm::TriggerResults>(triggerResultsInputTag_));
0063   triggers->setTriggerEventToken(consumes<trigger::TriggerEvent>(iConfig.getParameter<edm::InputTag>("TriggerEvent")));
0064   triggers->setDRMAX(DRMAX_HLT);
0065 
0066   edm::LogVerbatim("ZCounting") << "ZCounting::ZCounting set trigger names";
0067   const std::vector<std::string> patterns_ = iConfig.getParameter<std::vector<std::string>>("MuonTriggerNames");
0068   for (const std::string& pattern_ : patterns_) {
0069     triggers->addTriggerRecord(pattern_);
0070   }
0071 
0072   if (IDTypestr_ == "Loose")
0073     IDType_ = LooseID;
0074   else if (IDTypestr_ == "Medium")
0075     IDType_ = MediumID;
0076   else if (IDTypestr_ == "Tight")
0077     IDType_ = TightID;
0078   else if (IDTypestr_ == "CustomTight")
0079     IDType_ = CustomTightID;
0080   else
0081     IDType_ = NoneID;
0082 
0083   if (IsoTypestr_ == "Tracker-based")
0084     IsoType_ = TrackerIso;
0085   else if (IsoTypestr_ == "PF-based")
0086     IsoType_ = PFIso;
0087   else
0088     IsoType_ = NoneIso;
0089 }
0090 
0091 //
0092 //  -------------------------------------- Destructor --------------------------------------------
0093 //
0094 ZCounting::~ZCounting() { edm::LogInfo("ZCounting") << "Destructor ZCounting::~ZCounting " << std::endl; }
0095 
0096 //
0097 // -------------------------------------- beginRun --------------------------------------------
0098 //
0099 void ZCounting::dqmBeginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0100   edm::LogInfo("ZCounting") << "ZCounting::beginRun" << std::endl;
0101 
0102   // initialize triggers
0103 
0104   edm::LogVerbatim("ZCounting") << "ZCounting::dqmBeginRun now at " << iRun.id();
0105   bool hltChanged_ = true;
0106   if (hltConfigProvider_.init(iRun, iSetup, triggerResultsInputTag_.process(), hltChanged_)) {
0107     edm::LogVerbatim("ZCounting") << "ZCounting::dqmBeginRun [TriggerObjMatchValueMapsProducer::beginRun] "
0108                                      "HLTConfigProvider initialized [processName() = \""
0109                                   << hltConfigProvider_.processName() << "\", tableName() = \""
0110                                   << hltConfigProvider_.tableName() << "\", size() = " << hltConfigProvider_.size()
0111                                   << "]";
0112   } else {
0113     edm::LogError("ZCounting") << "ZCounting::dqmBeginRun Initialization of HLTConfigProvider failed for Run="
0114                                << iRun.id() << " (process=\"" << triggerResultsInputTag_.process()
0115                                << "\") -> plugin will not produce outputs for this Run";
0116     return;
0117   }
0118 
0119   triggers->initHLTObjects(hltConfigProvider_);
0120 }
0121 
0122 //
0123 // -------------------------------------- bookHistos --------------------------------------------
0124 //
0125 void ZCounting::bookHistograms(DQMStore::IBooker& ibooker_, edm::Run const&, edm::EventSetup const&) {
0126   edm::LogInfo("ZCounting") << "ZCounting::bookHistograms" << std::endl;
0127   ibooker_.cd();
0128   ibooker_.setCurrentFolder("ZCounting/Histograms");
0129 
0130   // Muon histograms
0131   h_mass_2HLT_BB = ibooker_.book2D("h_mass_2HLT_BB",
0132                                    "Both muon pass HLT in barrel-barrel",
0133                                    LumiBin_,
0134                                    LumiMin_,
0135                                    LumiMax_,
0136                                    MassBin_,
0137                                    MassMin_,
0138                                    MassMax_);
0139   h_mass_2HLT_BE = ibooker_.book2D("h_mass_2HLT_BE",
0140                                    "Both muon pass HLT passing in barrel-endcap",
0141                                    LumiBin_,
0142                                    LumiMin_,
0143                                    LumiMax_,
0144                                    MassBin_,
0145                                    MassMin_,
0146                                    MassMax_);
0147   h_mass_2HLT_EE = ibooker_.book2D("h_mass_2HLT_EE",
0148                                    "Both muon pass HLT passing in endcap-endcap",
0149                                    LumiBin_,
0150                                    LumiMin_,
0151                                    LumiMax_,
0152                                    MassBin_,
0153                                    MassMin_,
0154                                    MassMax_);
0155   h_mass_1HLT_BB = ibooker_.book2D("h_mass_1HLT_BB",
0156                                    "One muon pass HLT in barrel-barrel",
0157                                    LumiBin_,
0158                                    LumiMin_,
0159                                    LumiMax_,
0160                                    MassBin_,
0161                                    MassMin_,
0162                                    MassMax_);
0163   h_mass_1HLT_BE = ibooker_.book2D("h_mass_1HLT_BE",
0164                                    "One muon pass HLT passing in barrel-endcap",
0165                                    LumiBin_,
0166                                    LumiMin_,
0167                                    LumiMax_,
0168                                    MassBin_,
0169                                    MassMin_,
0170                                    MassMax_);
0171   h_mass_1HLT_EE = ibooker_.book2D("h_mass_1HLT_EE",
0172                                    "One muon pass HLT passing in endcap-endcap",
0173                                    LumiBin_,
0174                                    LumiMin_,
0175                                    LumiMax_,
0176                                    MassBin_,
0177                                    MassMin_,
0178                                    MassMax_);
0179 
0180   h_mass_ID_fail_BB = ibooker_.book2D(
0181       "h_mass_ID_fail_BB", "Muon ID failing barrel-barrel", LumiBin_, LumiMin_, LumiMax_, MassBin_, MassMin_, MassMax_);
0182   h_mass_ID_fail_BE = ibooker_.book2D(
0183       "h_mass_ID_fail_BE", "Muon ID failing barrel-endcap", LumiBin_, LumiMin_, LumiMax_, MassBin_, MassMin_, MassMax_);
0184 
0185   h_mass_ID_fail_EE = ibooker_.book2D(
0186       "h_mass_ID_fail_EE", "Muon ID failing endcap-endcap", LumiBin_, LumiMin_, LumiMax_, MassBin_, MassMin_, MassMax_);
0187 
0188   h_mass_Glo_pass_BB = ibooker_.book2D("h_mass_Glo_pass_BB",
0189                                        "Muon Glo passing barrel-barrel",
0190                                        LumiBin_,
0191                                        LumiMin_,
0192                                        LumiMax_,
0193                                        MassBin_,
0194                                        MassMin_,
0195                                        MassMax_);
0196   h_mass_Glo_pass_BE = ibooker_.book2D("h_mass_Glo_pass_BE",
0197                                        "Muon Glo passing barrel-endcap",
0198                                        LumiBin_,
0199                                        LumiMin_,
0200                                        LumiMax_,
0201                                        MassBin_,
0202                                        MassMin_,
0203                                        MassMax_);
0204 
0205   h_mass_Glo_pass_EE = ibooker_.book2D("h_mass_Glo_pass_EE",
0206                                        "Muon Glo passing endcap-endcap",
0207                                        LumiBin_,
0208                                        LumiMin_,
0209                                        LumiMax_,
0210                                        MassBin_,
0211                                        MassMin_,
0212                                        MassMax_);
0213 
0214   h_mass_Glo_fail_BB = ibooker_.book2D("h_mass_Glo_fail_BB",
0215                                        "Muon Glo failing barrel-barrel",
0216                                        LumiBin_,
0217                                        LumiMin_,
0218                                        LumiMax_,
0219                                        MassBin_,
0220                                        MassMin_,
0221                                        MassMax_);
0222   h_mass_Glo_fail_BE = ibooker_.book2D("h_mass_Glo_fail_BE",
0223                                        "Muon Glo failing barrel-endcap",
0224                                        LumiBin_,
0225                                        LumiMin_,
0226                                        LumiMax_,
0227                                        MassBin_,
0228                                        MassMin_,
0229                                        MassMax_);
0230 
0231   h_mass_Glo_fail_EE = ibooker_.book2D("h_mass_Glo_fail_EE",
0232                                        "Muon Glo failing endcap-endcap",
0233                                        LumiBin_,
0234                                        LumiMin_,
0235                                        LumiMax_,
0236                                        MassBin_,
0237                                        MassMin_,
0238                                        MassMax_);
0239 
0240   h_mass_Sta_pass_BB = ibooker_.book2D("h_mass_Sta_pass_BB",
0241                                        "Muon Sta passing barrel-barrel",
0242                                        LumiBin_,
0243                                        LumiMin_,
0244                                        LumiMax_,
0245                                        MassBin_,
0246                                        MassMin_,
0247                                        MassMax_);
0248   h_mass_Sta_pass_BE = ibooker_.book2D("h_mass_Sta_pass_BE",
0249                                        "Muon Sta passing barrel-endcap",
0250                                        LumiBin_,
0251                                        LumiMin_,
0252                                        LumiMax_,
0253                                        MassBin_,
0254                                        MassMin_,
0255                                        MassMax_);
0256 
0257   h_mass_Sta_pass_EE = ibooker_.book2D("h_mass_Sta_pass_EE",
0258                                        "Muon Sta passing endcap-endcap",
0259                                        LumiBin_,
0260                                        LumiMin_,
0261                                        LumiMax_,
0262                                        MassBin_,
0263                                        MassMin_,
0264                                        MassMax_);
0265 
0266   h_mass_Sta_fail_BB = ibooker_.book2D("h_mass_Sta_fail_BB",
0267                                        "Muon Sta failing barrel-barrel",
0268                                        LumiBin_,
0269                                        LumiMin_,
0270                                        LumiMax_,
0271                                        MassBin_,
0272                                        MassMin_,
0273                                        MassMax_);
0274   h_mass_Sta_fail_BE = ibooker_.book2D("h_mass_Sta_fail_BE",
0275                                        "Muon Sta failing barrel-endcap",
0276                                        LumiBin_,
0277                                        LumiMin_,
0278                                        LumiMax_,
0279                                        MassBin_,
0280                                        MassMin_,
0281                                        MassMax_);
0282 
0283   h_mass_Sta_fail_EE = ibooker_.book2D("h_mass_Sta_fail_EE",
0284                                        "Muon Sta failing endcap-endcap",
0285                                        LumiBin_,
0286                                        LumiMin_,
0287                                        LumiMax_,
0288                                        MassBin_,
0289                                        MassMin_,
0290                                        MassMax_);
0291 
0292   h_npv = ibooker_.book2D(
0293       "h_npv", "Events with valid primary vertex", LumiBin_, LumiMin_, LumiMax_, PVBin_, PVMin_, PVMax_);
0294 
0295   // Axis titles
0296   h_mass_2HLT_BB->setAxisTitle("luminosity section", 1);
0297   h_mass_2HLT_BE->setAxisTitle("luminosity section", 1);
0298   h_mass_2HLT_EE->setAxisTitle("luminosity section", 1);
0299   h_mass_1HLT_BB->setAxisTitle("luminosity section", 1);
0300   h_mass_1HLT_BE->setAxisTitle("luminosity section", 1);
0301   h_mass_1HLT_EE->setAxisTitle("luminosity section", 1);
0302   h_mass_ID_fail_BB->setAxisTitle("luminosity section", 1);
0303   h_mass_ID_fail_BE->setAxisTitle("luminosity section", 1);
0304   h_mass_ID_fail_EE->setAxisTitle("luminosity section", 1);
0305   h_mass_Glo_pass_BB->setAxisTitle("luminosity section", 1);
0306   h_mass_Glo_pass_BE->setAxisTitle("luminosity section", 1);
0307   h_mass_Glo_pass_EE->setAxisTitle("luminosity section", 1);
0308   h_mass_Glo_fail_BB->setAxisTitle("luminosity section", 1);
0309   h_mass_Glo_fail_BE->setAxisTitle("luminosity section", 1);
0310   h_mass_Glo_fail_EE->setAxisTitle("luminosity section", 1);
0311   h_mass_Sta_pass_BB->setAxisTitle("luminosity section", 1);
0312   h_mass_Sta_pass_BE->setAxisTitle("luminosity section", 1);
0313   h_mass_Sta_pass_EE->setAxisTitle("luminosity section", 1);
0314   h_mass_Sta_fail_BB->setAxisTitle("luminosity section", 1);
0315   h_mass_Sta_fail_BE->setAxisTitle("luminosity section", 1);
0316   h_mass_Sta_fail_EE->setAxisTitle("luminosity section", 1);
0317   h_mass_2HLT_BB->setAxisTitle("tag and probe mass", 2);
0318   h_mass_2HLT_BE->setAxisTitle("tag and probe mass", 2);
0319   h_mass_2HLT_EE->setAxisTitle("tag and probe mass", 2);
0320   h_mass_1HLT_BB->setAxisTitle("tag and probe mass", 2);
0321   h_mass_1HLT_BE->setAxisTitle("tag and probe mass", 2);
0322   h_mass_1HLT_EE->setAxisTitle("tag and probe mass", 2);
0323   h_mass_ID_fail_BB->setAxisTitle("tag and probe mass", 2);
0324   h_mass_ID_fail_BE->setAxisTitle("tag and probe mass", 2);
0325   h_mass_ID_fail_EE->setAxisTitle("tag and probe mass", 2);
0326   h_mass_Glo_pass_BB->setAxisTitle("tag and probe mass", 2);
0327   h_mass_Glo_pass_BE->setAxisTitle("tag and probe mass", 2);
0328   h_mass_Glo_pass_EE->setAxisTitle("tag and probe mass", 2);
0329   h_mass_Glo_fail_BB->setAxisTitle("tag and probe mass", 2);
0330   h_mass_Glo_fail_BE->setAxisTitle("tag and probe mass", 2);
0331   h_mass_Glo_fail_EE->setAxisTitle("tag and probe mass", 2);
0332   h_mass_Sta_pass_BB->setAxisTitle("tag and probe mass", 2);
0333   h_mass_Sta_pass_BE->setAxisTitle("tag and probe mass", 2);
0334   h_mass_Sta_pass_EE->setAxisTitle("tag and probe mass", 2);
0335   h_mass_Sta_fail_BB->setAxisTitle("tag and probe mass", 2);
0336   h_mass_Sta_fail_BE->setAxisTitle("tag and probe mass", 2);
0337   h_mass_Sta_fail_EE->setAxisTitle("tag and probe mass", 2);
0338   h_npv->setAxisTitle("luminosity section", 1);
0339   h_npv->setAxisTitle("number of primary vertices", 2);
0340 }
0341 
0342 //
0343 // -------------------------------------- Analyze --------------------------------------------
0344 //
0345 //--------------------------------------------------------------------------------------------------
0346 void ZCounting::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {  // Fill event tree on the fly
0347   edm::LogInfo("ZCounting") << "ZCounting::analyze" << std::endl;
0348 
0349   //-------------------------------
0350   //--- Vertex
0351   //-------------------------------
0352   edm::Handle<reco::VertexCollection> hVertexProduct;
0353   iEvent.getByToken(fPVName_token, hVertexProduct);
0354   if (!hVertexProduct.isValid()) {
0355     edm::LogWarning("ZCounting") << "ZCounting::analyze - no valid primary vertex product found" << std::endl;
0356     return;
0357   }
0358 
0359   const reco::Vertex* pv = nullptr;
0360   int nvtx = 0;
0361 
0362   for (auto const& itVtx : *hVertexProduct) {
0363     if (itVtx.isFake())
0364       continue;
0365     if (itVtx.tracksSize() < VtxNTracksFitCut_)
0366       continue;
0367     if (itVtx.ndof() < VtxNdofCut_)
0368       continue;
0369     if (fabs(itVtx.z()) > VtxAbsZCut_)
0370       continue;
0371     if (itVtx.position().Rho() > VtxRhoCut_)
0372       continue;
0373 
0374     if (nvtx == 0) {
0375       pv = &itVtx;
0376     }
0377     nvtx++;
0378   }
0379 
0380   h_npv->Fill(iEvent.luminosityBlock(), nvtx);
0381 
0382   //-------------------------------
0383   //--- Trigger
0384   //-------------------------------
0385   triggers->readEvent(iEvent);
0386 
0387   // Trigger requirement
0388   if (!triggers->pass())
0389     return;
0390 
0391   //-------------------------------
0392   //--- Muon and Track collections
0393   //-------------------------------
0394   edm::Handle<reco::MuonCollection> hMuonProduct;
0395   iEvent.getByToken(fMuonName_token, hMuonProduct);
0396   if (!hMuonProduct.isValid()) {
0397     edm::LogWarning("ZCounting") << "ZCounting::analyze - no valid hMuonProduct found" << std::endl;
0398     return;
0399   }
0400 
0401   edm::Handle<reco::TrackCollection> hTrackProduct;
0402   iEvent.getByToken(fTrackName_token, hTrackProduct);
0403   if (!hTrackProduct.isValid()) {
0404     edm::LogWarning("ZCounting") << "ZCounting::analyze - no valid hTrackProduct found" << std::endl;
0405     return;
0406   }
0407 
0408   //-------------------------------
0409   //--- Merged standalone muon collections
0410   //--- The muon collection contains duplicates (from standAloneMuons and standAloneMuons:UpdatedAtVtx collections) and missing standAloneMuons
0411   //--- We need to produce a merged standalone muon collection to reproduce the decision in the global muon producer
0412   //-------------------------------
0413   edm::Handle<reco::TrackCollection> tracksStandAlone;
0414   iEvent.getByToken(fStandaloneRegName_token, tracksStandAlone);
0415   if (!tracksStandAlone.isValid()) {
0416     edm::LogWarning("ZCounting") << "ZCounting::analyze - no valid tracksStandAlone found" << std::endl;
0417     return;
0418   }
0419 
0420   edm::Handle<reco::TrackCollection> tracksStandAloneUpdatedAtVtx;
0421   iEvent.getByToken(fStandaloneUpdName_token, tracksStandAloneUpdatedAtVtx);
0422   if (!tracksStandAloneUpdatedAtVtx.isValid()) {
0423     edm::LogWarning("ZCounting") << "ZCounting::analyze - no valid tracksStandAloneUpdatedAtVtx found" << std::endl;
0424     return;
0425   }
0426 
0427   std::vector<const reco::Track*> hStandaloneProduct;
0428   std::vector<bool> passGlobalMuonMap;
0429 
0430   for (auto const& standAlone : *tracksStandAlone) {
0431     auto const extraIdx = standAlone.extra().key();
0432 
0433     const reco::Track* track = &standAlone;
0434 
0435     // replicate logic in GlobalMuonProducer, take the updatedAtVtx track if it exists and has
0436     // the same eta sign as the original, otherwise take the original
0437     for (auto const& standAloneUpdatedAtVtx : *tracksStandAloneUpdatedAtVtx) {
0438       if (standAloneUpdatedAtVtx.extra().key() == extraIdx) {
0439         const bool etaFlip1 = (standAloneUpdatedAtVtx.eta() * standAlone.eta()) >= 0;
0440         if (etaFlip1) {
0441           track = &standAloneUpdatedAtVtx;
0442         }
0443         break;
0444       }
0445     }
0446 
0447     // kinematic cuts
0448     if (track->pt() < MIN_PT_STA)
0449       continue;
0450     if (fabs(track->eta()) > MAX_ETA_STA)
0451       continue;
0452     // require minimum number of valid hits (mainly to reduce background)
0453     if (track->numberOfValidHits() < N_STA_HITS)
0454       continue;
0455 
0456     // look for corresponding muon object to check if the standalone muon is global
0457     bool isGlobalMuon = false;
0458     for (auto const& itMu2 : *hMuonProduct) {
0459       if (itMu2.standAloneMuon().isNull())
0460         continue;
0461 
0462       auto const& muonStandAlone = *itMu2.standAloneMuon();
0463 
0464       if (track->extra().key() == muonStandAlone.extra().key()) {
0465         // we found a corresponding muon object
0466         if (muonStandAlone.pt() == track->pt() && muonStandAlone.eta() == track->eta() &&
0467             muonStandAlone.phi() == track->phi()) {
0468           // the corresponding muon object uses the same standalone muon track
0469           // check if is a global muon
0470           isGlobalMuon = passGlobalMuon(itMu2);
0471         }
0472         break;
0473       }
0474     }
0475 
0476     passGlobalMuonMap.push_back(isGlobalMuon);
0477     hStandaloneProduct.push_back(track);
0478   }
0479 
0480   TLorentzVector vTag(0., 0., 0., 0.);
0481   TLorentzVector vProbe(0., 0., 0., 0.);
0482   TLorentzVector vTrack(0., 0., 0., 0.);
0483 
0484   // Tag loop
0485   for (auto const& itMu1 : *hMuonProduct) {
0486     const float pt1 = itMu1.muonBestTrack()->pt();
0487     const float eta1 = itMu1.muonBestTrack()->eta();
0488     const float phi1 = itMu1.muonBestTrack()->phi();
0489     const float q1 = itMu1.muonBestTrack()->charge();
0490 
0491     // Tag selection: kinematic cuts, lepton selection and trigger matching
0492     if (pt1 < PtCutL1_)
0493       continue;
0494     if (fabs(eta1) > EtaCutL1_)
0495       continue;
0496     if (!(passGlobalMuon(itMu1) && passMuonID(itMu1, pv) && passMuonIso(itMu1)))
0497       continue;
0498     if (!triggers->passObj(eta1, phi1))
0499       continue;
0500 
0501     vTag.SetPtEtaPhiM(pt1, eta1, phi1, MUON_MASS);
0502 
0503     bool isTagCentral = false;
0504     if (fabs(eta1) < MUON_BOUND)
0505       isTagCentral = true;
0506 
0507     // Probe loop over muons
0508     for (auto const& itMu2 : *hMuonProduct) {
0509       if (&itMu2 == &itMu1)
0510         continue;
0511 
0512       const float pt2 = itMu2.muonBestTrack()->pt();
0513       const float eta2 = itMu2.muonBestTrack()->eta();
0514       const float phi2 = itMu2.muonBestTrack()->phi();
0515       const float q2 = itMu2.muonBestTrack()->charge();
0516 
0517       // Probe selection: kinematic cuts and opposite charge requirement
0518       if (pt2 < PtCutL2_)
0519         continue;
0520       if (fabs(eta2) > EtaCutL2_)
0521         continue;
0522       if (q1 == q2)
0523         continue;
0524 
0525       vProbe.SetPtEtaPhiM(pt2, eta2, phi2, MUON_MASS);
0526 
0527       // Mass window
0528       TLorentzVector vDilep = vTag + vProbe;
0529       float dilepMass = vDilep.M();
0530       if ((dilepMass < MassMin_) || (dilepMass > MassMax_))
0531         continue;
0532 
0533       bool isProbeCentral = fabs(eta2) < MUON_BOUND;
0534 
0535       // Determine event category for efficiency calculation
0536       if (passGlobalMuon(itMu2) && passMuonID(itMu2, pv) && passMuonIso(itMu2)) {
0537         if (triggers->passObj(eta2, phi2)) {
0538           // category 2HLT: both muons passing trigger requirements
0539           if (&itMu1 > &itMu2)
0540             continue;  // make sure we don't double count MuMu2HLT category
0541 
0542           if (isTagCentral && isProbeCentral) {
0543             h_mass_2HLT_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0544           } else if (!isTagCentral && !isProbeCentral) {
0545             h_mass_2HLT_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0546           } else {
0547             h_mass_2HLT_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0548           }
0549         } else {
0550           // category 1HLT: only one muon passes trigger
0551           if (isTagCentral && isProbeCentral) {
0552             h_mass_1HLT_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0553           } else if (!isTagCentral && !isProbeCentral) {
0554             h_mass_1HLT_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0555           } else {
0556             h_mass_1HLT_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0557           }
0558         }
0559       } else if (passGlobalMuon(itMu2)) {
0560         // category Glo: probe is a Global muon but failing selection
0561         if (isTagCentral && isProbeCentral) {
0562           h_mass_ID_fail_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0563         } else if (!isTagCentral && !isProbeCentral) {
0564           h_mass_ID_fail_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0565         } else {
0566           h_mass_ID_fail_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0567         }
0568       }
0569     }  // End of probe loop over muons
0570 
0571     // Probe loop over standalone muons, for global muon efficiency calculation
0572     for (std::vector<reco::Track>::size_type idx = 0; idx < hStandaloneProduct.size(); idx++) {
0573       const reco::Track* itSta = hStandaloneProduct[idx];
0574 
0575       // standalone muon kinematics
0576       const float pt2 = itSta->pt();
0577       const float eta2 = itSta->eta();
0578       const float phi2 = itSta->phi();
0579 
0580       // kinematic cuts
0581       if (pt2 < PtCutL2_)
0582         continue;
0583       if (fabs(eta2) > EtaCutL2_)
0584         continue;
0585 
0586       vProbe.SetPtEtaPhiM(pt2, eta2, phi2, MUON_MASS);
0587 
0588       // Mass window
0589       TLorentzVector vDilep = vTag + vProbe;
0590       float dilepMass = vDilep.M();
0591       if ((dilepMass < MassMin_) || (dilepMass > MassMax_))
0592         continue;
0593 
0594       const bool isProbeCentral = fabs(eta2) < MUON_BOUND;
0595 
0596       if (passGlobalMuonMap[idx]) {
0597         if (isTagCentral && isProbeCentral) {
0598           h_mass_Glo_pass_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0599         } else if (!isTagCentral && !isProbeCentral) {
0600           h_mass_Glo_pass_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0601         } else {
0602           h_mass_Glo_pass_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0603         }
0604       } else {
0605         if (isTagCentral && isProbeCentral) {
0606           h_mass_Glo_fail_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0607         } else if (!isTagCentral && !isProbeCentral) {
0608           h_mass_Glo_fail_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0609         } else {
0610           h_mass_Glo_fail_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0611         }
0612       }
0613     }
0614 
0615     // Probe loop over tracks, only for standalone efficiency calculation
0616     for (auto const& itTrk : *hTrackProduct) {
0617       const float pt2 = itTrk.pt();
0618       const float eta2 = itTrk.eta();
0619       const float phi2 = itTrk.phi();
0620       const float q2 = itTrk.charge();
0621 
0622       // Probe selection:  kinematic cuts and opposite charge requirement
0623       if (pt2 < PtCutL2_)
0624         continue;
0625       if (fabs(eta2) > EtaCutL2_)
0626         continue;
0627       if (q1 == q2)
0628         continue;
0629       if (!passTrack(itTrk))
0630         continue;
0631 
0632       vTrack.SetPtEtaPhiM(pt2, eta2, phi2, MUON_MASS);
0633 
0634       TLorentzVector vDilep = vTag + vTrack;
0635       float dilepMass = vDilep.M();
0636       if ((dilepMass < MassMin_) || (dilepMass > MassMax_))
0637         continue;
0638 
0639       // check if track is matched to standalone muon
0640       bool isStandalone = false;
0641       for (const reco::Track* itSta : hStandaloneProduct) {
0642         if (reco::deltaR2(itSta->eta(), itSta->phi(), eta2, phi2) < DRMAX_IO) {
0643           isStandalone = true;
0644           break;
0645         }
0646       }
0647 
0648       const bool isTrackCentral = fabs(eta2) < MUON_BOUND;
0649 
0650       if (isStandalone) {
0651         if (isTagCentral && isTrackCentral) {
0652           h_mass_Sta_pass_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0653         } else if (!isTagCentral && !isTrackCentral) {
0654           h_mass_Sta_pass_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0655         } else {
0656           h_mass_Sta_pass_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0657         }
0658       } else {
0659         if (isTagCentral && isTrackCentral) {
0660           h_mass_Sta_fail_BB->Fill(iEvent.luminosityBlock(), dilepMass);
0661         } else if (!isTagCentral && !isTrackCentral) {
0662           h_mass_Sta_fail_EE->Fill(iEvent.luminosityBlock(), dilepMass);
0663         } else {
0664           h_mass_Sta_fail_BE->Fill(iEvent.luminosityBlock(), dilepMass);
0665         }
0666       }
0667     }  //End of probe loop over tracks
0668   }    //End of tag loop
0669 }
0670 
0671 //
0672 // -------------------------------------- functions --------------------------------------------
0673 //
0674 
0675 //--------------------------------------------------------------------------------------------------
0676 // Definition of the CustomTightID function
0677 bool ZCounting::isCustomTightMuon(const reco::Muon& muon) {
0678   // tight POG cut based ID w/o impact parameter cuts
0679   return muon.isGlobalMuon() && muon.isPFMuon() && muon.globalTrack()->normalizedChi2() < 10. &&
0680          muon.globalTrack()->hitPattern().numberOfValidMuonHits() > 0 && muon.numberOfMatchedStations() > 1 &&
0681          muon.innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5 &&
0682          muon.innerTrack()->hitPattern().numberOfValidPixelHits() > 0;
0683 }
0684 
0685 //--------------------------------------------------------------------------------------------------
0686 bool ZCounting::passMuonID(const reco::Muon& muon, const reco::Vertex* vtx) {
0687   // Muon ID selection, using internal function "DataFormats/MuonReco/src/MuonSelectors.cc
0688   switch (IDType_) {
0689     case LooseID:
0690       return muon::isLooseMuon(muon);
0691     case MediumID:
0692       return muon::isMediumMuon(muon);
0693     case CustomTightID:
0694       return isCustomTightMuon(muon);
0695     case TightID:
0696       return vtx != nullptr && muon::isTightMuon(muon, *vtx);
0697     case NoneID:
0698       return true;
0699   }
0700   return false;
0701 }
0702 
0703 //--------------------------------------------------------------------------------------------------
0704 bool ZCounting::passGlobalMuon(const reco::Muon& muon) {
0705   // Global muon selection:
0706   // - standard global muon criterium,
0707   // - requirements on inner and outer track pT>15 and |eta|
0708   // - requirements on deltaR(inner track, outer track)
0709 
0710   return muon.isGlobalMuon() && muon.outerTrack()->numberOfValidHits() >= N_STA_HITS &&
0711          muon.innerTrack()->pt() > MIN_PT_TRK && std::abs(muon.innerTrack()->eta()) < MAX_ETA_TRK &&
0712          muon.outerTrack()->pt() > MIN_PT_STA && std::abs(muon.outerTrack()->eta()) < MAX_ETA_STA &&
0713          reco::deltaR2(
0714              muon.outerTrack()->eta(), muon.outerTrack()->phi(), muon.innerTrack()->eta(), muon.innerTrack()->phi()) <
0715              DRMAX_IO;
0716 }
0717 
0718 //--------------------------------------------------------------------------------------------------
0719 bool ZCounting::passTrack(const reco::Track& track) {
0720   return track.hitPattern().trackerLayersWithMeasurement() >= 6 && track.hitPattern().numberOfValidPixelHits() >= 1 &&
0721          track.originalAlgo() != 13      // reject muon seeded tracks - InOut
0722          && track.originalAlgo() != 14;  // reject muon seeded tracks - OutIn
0723 }
0724 
0725 //--------------------------------------------------------------------------------------------------
0726 bool ZCounting::passMuonIso(const reco::Muon& muon) {
0727   //Muon isolation selection, up-to-date with MUO POG recommendation
0728   switch (IsoType_) {
0729     case TrackerIso:
0730       return muon.isolationR03().sumPt < IsoCut_;
0731     case PFIso:
0732       return muon.pfIsolationR04().sumChargedHadronPt +
0733                  std::max(0.,
0734                           muon.pfIsolationR04().sumNeutralHadronEt + muon.pfIsolationR04().sumPhotonEt -
0735                               0.5 * muon.pfIsolationR04().sumPUPt) <
0736              IsoCut_;
0737     case NoneIso:
0738       return true;
0739   }
0740 
0741   return false;
0742 }
0743 
0744 DEFINE_FWK_MODULE(ZCounting);