MuonDPG/src/BaseTnPEfficiencyTask.cc

0001 /*
0002  * \file BaseTnPEfficiencyTask.cc
0003  *
0004  * \author L. Lunerti - INFN Bologna
0005  *
0006  */
0007
0008 #include "DQMOffline/MuonDPG/interface/BaseTnPEfficiencyTask.h"
0009
0010 // Framework
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012
0013 // Geometry
0014 #include "DataFormats/GeometryVector/interface/Pi.h"
0015
0016 //Math
0017 #include "DataFormats/Math/interface/deltaR.h"
0018
0019 //Root
0020 #include "TRegexp.h"
0021
0022 #include <tuple>
0023 #include <algorithm>
0024
0025 BaseTnPEfficiencyTask::BaseTnPEfficiencyTask(const edm::ParameterSet& config)
0026     : m_nEvents(0),
0027       m_muToken(consumes<reco::MuonCollection>(config.getUntrackedParameter<edm::InputTag>("inputTagMuons"))),
0028       m_borderCut(config.getUntrackedParameter<double>("borderCut")),
0029       m_dxCut(config.getUntrackedParameter<double>("dx_cut")),
0030       m_detailedAnalysis(config.getUntrackedParameter<bool>("detailedAnalysis")),
0031       m_primaryVerticesToken(
0032           consumes<std::vector<reco::Vertex>>(config.getUntrackedParameter<edm::InputTag>("inputTagPrimaryVertices"))),
0033       m_triggerResultsToken(
0034           consumes<edm::TriggerResults>(config.getUntrackedParameter<edm::InputTag>("trigResultsTag"))),
0035       m_triggerEventToken(consumes<trigger::TriggerEvent>(config.getUntrackedParameter<edm::InputTag>("trigEventTag"))),
0036       m_trigName(config.getUntrackedParameter<std::string>("trigName")),
0037       m_probeSelector(config.getUntrackedParameter<std::string>("probeCut")),
0038       m_dxyCut(config.getUntrackedParameter<double>("probeDxyCut")),
0039       m_dzCut(config.getUntrackedParameter<double>("probeDzCut")),
0040       m_tagSelector(config.getUntrackedParameter<std::string>("tagCut")),
0041       m_lowPairMassCut(config.getUntrackedParameter<double>("lowPairMassCut")),
0042       m_highPairMassCut(config.getUntrackedParameter<double>("highPairMassCut")) {
0043   LogTrace("DQMOffline|MuonDPG|BaseTnPEfficiencyTask") << "[BaseTnPEfficiencyTask]: Constructor" << std::endl;
0044 }
0045
0046 BaseTnPEfficiencyTask::~BaseTnPEfficiencyTask() {
0047   LogTrace("DQMOffline|MuonDPG|BaseTnPEfficiencyTask")
0048       << "[BaseTnPEfficiencyTask]: analyzed " << m_nEvents << " events" << std::endl;
0049 }
0050
0051 void BaseTnPEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& context) {
0052   bool changed = true;
0053   m_hltConfig.init(run, context, "HLT", changed);
0054
0055   bool enableWildCard = true;
0056
0057   TString tName = TString(m_trigName);
0058   TRegexp tNamePattern = TRegexp(tName, enableWildCard);
0059
0060   for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
0061     TString pathName = TString(m_hltConfig.triggerName(iPath));
0062     if (pathName.Contains(tNamePattern)) {
0063       m_trigIndices.push_back(static_cast<int>(iPath));
0064     }
0065   }
0066 }
0067
0068 void BaseTnPEfficiencyTask::bookHistograms(DQMStore::IBooker& iBooker,
0069                                            edm::Run const& run,
0070                                            edm::EventSetup const& context) {
0071   LogTrace("DQMOffline|MuonDPG|BaseTnPEfficiencyTask") << "[BaseTnPEfficiencyTask]: bookHistograms" << std::endl;
0072
0073   if (m_detailedAnalysis) {
0074     std::string baseDir = topFolder() + "/detailed/";
0075     iBooker.setCurrentFolder(baseDir);
0076
0077     LogTrace("DQMOffline|MuonDPG|BaseTnPEfficiencyTask")
0078         << "[BaseTnPEfficiencyTask]: booking histos in " << baseDir << std::endl;
0079
0080     m_histos["probePt"] = iBooker.book1D("probePt", "probePt;probe p_{T} [GeV];Events", 125, 0., 250.);
0081     m_histos["probeEta"] = iBooker.book1D("probeEta", "probeEta;probe #eta;Events", 24, -2.4, 2.4);
0082     m_histos["probePhi"] = iBooker.book1D("probePhi", "probePhi;probe #phi; Events", 36, -TMath::Pi(), TMath::Pi());
0083     m_histos["probeNumberOfMatchedStations"] = iBooker.book1D(
0084         "probeNumberOfMatchedStations", "probeNumberOfMatchedStations;Number of matched stations;Events", 5, 0., 5);
0085     m_histos["pairMass"] = iBooker.book1D("pairMass", "pairMass", 25, 70., 120.);
0086   }
0087 }
0088
0089 void BaseTnPEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& context) {
0090   ++m_nEvents;
0091
0092   edm::Handle<reco::MuonCollection> muons;
0093   if (!event.getByToken(m_muToken, muons))
0094     return;
0095
0096   edm::Handle<std::vector<reco::Vertex>> vtxs;
0097   if (!event.getByToken(m_primaryVerticesToken, vtxs))
0098     return;
0099   const reco::Vertex& vertex = vtxs->at(0);
0100
0101   edm::Handle<edm::TriggerResults> triggerResults;
0102   if (!event.getByToken(m_triggerResultsToken, triggerResults))
0103     return;
0104
0105   edm::Handle<trigger::TriggerEvent> triggerEvent;
0106   if (!event.getByToken(m_triggerEventToken, triggerEvent))
0107     return;
0108
0109   //common tnp variables
0110   std::vector<unsigned> preSel_tag_indices;
0111   std::vector<unsigned> tag_indices;
0112   std::vector<unsigned> preSel_probe_indices;
0113   std::vector<unsigned> probe_indices;
0114
0115   if (muons.isValid() && vtxs.isValid()) {
0116     //Is there a better way to initialize two different type variables?
0117     for (auto [muon, muonColl_index] = std::tuple{std::vector<reco::Muon>::const_iterator{(*muons).begin()}, 0};
0118          muon != (*muons).end();
0119          ++muon, ++muonColl_index) {
0120       bool trigMatch = false;
0121
0122       //Getting trigger infos for tag selection
0123       if (triggerResults.isValid() && triggerEvent.isValid()) {
0124         const trigger::TriggerObjectCollection trigObjColl = triggerEvent->getObjects();
0125         trigMatch = hasTrigger(m_trigIndices, trigObjColl, triggerEvent, *muon);
0126       }
0127
0128       //Probe selection
0129       if (m_probeSelector(*muon) && (std::abs(muon->muonBestTrack()->dxy(vertex.position())) < m_dxyCut) &&
0130           (std::abs(muon->muonBestTrack()->dz(vertex.position())) < m_dzCut)) {
0131         preSel_probe_indices.push_back(muonColl_index);
0132       }
0133       //Tag selection
0134       if (m_tagSelector(*muon) && trigMatch) {
0135         preSel_tag_indices.push_back(muonColl_index);
0136       }
0137
0138     }  //loop over muons
0139   }
0140
0141   //Probe selection
0142   for (const auto i_tag : preSel_tag_indices) {
0143     reco::Muon tag = (*muons).at(i_tag);
0144     float pt_max = 0.;
0145     int max_pt_idx;
0146     bool pair_found = false;
0147
0148     for (const auto i_probe : preSel_probe_indices) {
0149       //Prevent tag and probe to be the same object
0150       if (i_probe == i_tag)
0151         continue;
0152
0153       reco::Muon preSel_probe = (*muons).at(i_probe);
0154
0155       int pair_charge_product = tag.charge() * preSel_probe.charge();
0156
0157       //check if tag+probe pair is compatible with Z decay
0158       if (pair_charge_product > 0)
0159         continue;
0160
0161       float pair_mass = (tag.polarP4() + preSel_probe.polarP4()).M();
0162       m_histos.find("pairMass")->second->Fill(pair_mass);
0163
0164       if (pair_mass < m_lowPairMassCut || pair_mass > m_highPairMassCut)
0165         continue;
0166
0167       float pair_pt = (tag.polarP4() + preSel_probe.polarP4()).Pt();
0168       if (pair_pt > pt_max) {
0169         pair_found = true;
0170         pt_max = pair_pt;
0171         max_pt_idx = i_probe;
0172       }
0173     }
0174     if (pair_found) {
0175       probe_indices.push_back(max_pt_idx);
0176       tag_indices.push_back(i_tag);
0177     }
0178   }
0179
0180   m_probeIndices.push_back(probe_indices);
0181   m_tagIndices.push_back(tag_indices);
0182 }
0183
0184 bool BaseTnPEfficiencyTask::hasTrigger(std::vector<int>& trigIndices,
0185                                        const trigger::TriggerObjectCollection& trigObjs,
0186                                        edm::Handle<trigger::TriggerEvent>& trigEvent,
0187                                        const reco::Muon& muon) {
0188   float dR2match = 999.;
0189   for (int trigIdx : trigIndices) {
0190     const std::vector<std::string> trigModuleLabels = m_hltConfig.moduleLabels(trigIdx);
0191
0192     const unsigned trigModuleIndex =
0193         std::find(trigModuleLabels.begin(), trigModuleLabels.end(), "hltBoolEnd") - trigModuleLabels.begin() - 1;
0194     const unsigned hltFilterIndex = trigEvent->filterIndex(edm::InputTag(trigModuleLabels[trigModuleIndex], "", "HLT"));
0195     if (hltFilterIndex < trigEvent->sizeFilters()) {
0196       const trigger::Keys keys = trigEvent->filterKeys(hltFilterIndex);
0197       const trigger::Vids vids = trigEvent->filterIds(hltFilterIndex);
0198       const unsigned nTriggers = vids.size();
0199
0200       for (unsigned iTrig = 0; iTrig < nTriggers; ++iTrig) {
0201         trigger::TriggerObject trigObj = trigObjs[keys[iTrig]];
0202         float dR2 = deltaR2(muon, trigObj);
0203         if (dR2 < dR2match)
0204           dR2match = dR2;
0205       }
0206     }
0207   }
0208
0209   return dR2match < 0.01;
0210 }