Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /** \class MuRecoAnalyzer
0002  *
0003  *  DQM monitoring source for PF muons
0004  *
0005  *  \author C. Battilana - CIEMAT
0006  */
0007 //Base class
0008 #include "DQMOffline/Muon/interface/MuonPFAnalyzer.h"
0009 
0010 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0011 
0012 #include "DataFormats/GeometryVector/interface/Pi.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 
0015 //System included files
0016 #include <memory>
0017 #include <string>
0018 #include <typeinfo>
0019 #include <utility>
0020 
0021 //Root included files
0022 #include "TH1.h"
0023 #include "TH2.h"
0024 #include "TProfile.h"
0025 
0026 //Event framework included files
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 #include "DQMServices/Core/interface/DQMStore.h"
0033 #include "FWCore/ServiceRegistry/interface/Service.h"
0034 
0035 using namespace edm;
0036 using namespace std;
0037 using namespace reco;
0038 
0039 MuonPFAnalyzer::MuonPFAnalyzer(const ParameterSet &pSet) {
0040   LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] Initializing configuration from parameterset.\n";
0041 
0042   theGenLabel_ = consumes<GenParticleCollection>(pSet.getParameter<InputTag>("inputTagGenParticles"));
0043   theRecoLabel_ = consumes<MuonCollection>(pSet.getParameter<InputTag>("inputTagMuonReco"));
0044   theVertexLabel_ = consumes<VertexCollection>(pSet.getParameter<InputTag>("inputTagVertex"));
0045   theBeamSpotLabel_ = consumes<BeamSpot>(pSet.getParameter<InputTag>("inputTagBeamSpot"));
0046 
0047   theHighPtTh = pSet.getParameter<double>("highPtThreshold");
0048   theRecoGenR = pSet.getParameter<double>("recoGenDeltaR");
0049   theIsoCut = pSet.getParameter<double>("relCombIsoCut");
0050   theRunOnMC = pSet.getParameter<bool>("runOnMC");
0051 
0052   theFolder = pSet.getParameter<string>("folder");
0053 
0054   theMuonKinds.push_back("");          // all TUNEP/PF muons
0055   theMuonKinds.push_back("Tight");     // tight TUNEP/PF muons
0056   theMuonKinds.push_back("TightIso");  // tight/iso TUNEP/PF muons
0057 }
0058 
0059 MuonPFAnalyzer::~MuonPFAnalyzer() { LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] Destructor called.\n"; }
0060 
0061 // ------------ method called when starting to processes a run  ------------
0062 void MuonPFAnalyzer::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &, edm::EventSetup const &) {
0063   if (theRunOnMC) {
0064     bookHistos(ibooker, "PF");
0065     bookHistos(ibooker, "PFTight");
0066     bookHistos(ibooker, "PFTightIso");
0067     bookHistos(ibooker, "TUNEP");
0068     bookHistos(ibooker, "TUNEPTight");
0069     bookHistos(ibooker, "TUNEPTightIso");
0070   }
0071 
0072   bookHistos(ibooker, "PFvsTUNEP");
0073   bookHistos(ibooker, "PFvsTUNEPTight");
0074   bookHistos(ibooker, "PFvsTUNEPTightIso");
0075 }
0076 
0077 void MuonPFAnalyzer::analyze(const Event &event, const EventSetup &context) {
0078   Handle<reco::MuonCollection> muons;
0079   event.getByToken(theRecoLabel_, muons);
0080 
0081   Handle<GenParticleCollection> genMuons;
0082   event.getByToken(theGenLabel_, genMuons);
0083 
0084   Handle<BeamSpot> beamSpot;
0085   event.getByToken(theBeamSpotLabel_, beamSpot);
0086 
0087   Handle<VertexCollection> vertex;
0088   event.getByToken(theVertexLabel_, vertex);
0089 
0090   const Vertex primaryVertex = getPrimaryVertex(vertex, beamSpot);
0091 
0092   recoToGenMatch(muons, genMuons);
0093 
0094   RecoGenCollection::const_iterator recoGenIt = theRecoGen.begin();
0095   RecoGenCollection::const_iterator recoGenEnd = theRecoGen.end();
0096 
0097   for (; recoGenIt != recoGenEnd; ++recoGenIt) {
0098     const Muon *muon = recoGenIt->first;
0099     TrackRef tunePTrack = muon->tunePMuonBestTrack();
0100 
0101     const GenParticle *genMuon = recoGenIt->second;
0102 
0103     vector<string>::const_iterator kindIt = theMuonKinds.begin();
0104     vector<string>::const_iterator kindEnd = theMuonKinds.end();
0105 
0106     for (; kindIt != kindEnd; ++kindIt) {
0107       const string &kind = (*kindIt);
0108 
0109       if (kind.find("Tight") != string::npos && !muon::isTightMuon((*muon), primaryVertex))
0110         continue;
0111 
0112       if (kind.find("Iso") != string::npos && combRelIso(muon) > theIsoCut)
0113         continue;
0114 
0115       if (theRunOnMC && genMuon && !muon->innerTrack().isNull())  // has matched gen muon
0116       {
0117         if (!tunePTrack.isNull()) {
0118           string group = "TUNEP" + kind;
0119 
0120           float pt = tunePTrack->pt();
0121           float phi = tunePTrack->phi();
0122           float eta = tunePTrack->eta();
0123 
0124           float genPt = genMuon->pt();
0125           float genPhi = genMuon->p4().phi();
0126           float genEta = genMuon->p4().eta();
0127 
0128           float dPtOverPt = (pt / genPt) - 1;
0129 
0130           if (pt < theHighPtTh) {
0131             fillInRange(getPlot(group, "code"), 1, muonTrackType(muon, false));
0132             fillInRange(getPlot(group, "deltaPtOverPt"), 1, dPtOverPt);
0133           } else {
0134             fillInRange(getPlot(group, "codeHighPt"), 1, muonTrackType(muon, false));
0135             fillInRange(getPlot(group, "deltaPtOverPtHighPt"), 1, dPtOverPt);
0136           }
0137 
0138           fillInRange(getPlot(group, "deltaPt"), 1, (pt - genPt));
0139           fillInRange(getPlot(group, "deltaPhi"), 1, fDeltaPhi(genPhi, phi));
0140           fillInRange(getPlot(group, "deltaEta"), 1, genEta - eta);
0141         }
0142 
0143         if (muon->isPFMuon()) {
0144           string group = "PF" + kind;
0145 
0146           // Assumes that default in muon is PF
0147           float pt = muon->pt();
0148           float phi = muon->p4().phi();
0149           float eta = muon->p4().eta();
0150 
0151           float genPt = genMuon->pt();
0152           float genPhi = genMuon->p4().phi();
0153           float genEta = genMuon->p4().eta();
0154 
0155           float dPtOverPt = (pt / genPt) - 1;
0156 
0157           if (pt < theHighPtTh) {
0158             fillInRange(getPlot(group, "code"), 1, muonTrackType(muon, true));
0159             fillInRange(getPlot(group, "deltaPtOverPt"), 1, dPtOverPt);
0160           } else {
0161             fillInRange(getPlot(group, "codeHighPt"), 1, muonTrackType(muon, true));
0162             fillInRange(getPlot(group, "deltaPtOverPtHighPt"), 1, dPtOverPt);
0163           }
0164 
0165           fillInRange(getPlot(group, "deltaPt"), 1, pt - genPt);
0166           fillInRange(getPlot(group, "deltaPhi"), 1, fDeltaPhi(genPhi, phi));
0167           fillInRange(getPlot(group, "deltaEta"), 1, genEta - eta);
0168         }
0169       }
0170 
0171       if (muon->isPFMuon() && !tunePTrack.isNull() && !muon->innerTrack().isNull())  // Compare PF with TuneP + Tracker
0172       {                                                                              // No gen matching needed
0173 
0174         string group = "PFvsTUNEP" + kind;
0175 
0176         float pt = tunePTrack->pt();
0177         float phi = tunePTrack->phi();
0178         float eta = tunePTrack->eta();
0179 
0180         // Assumes that default in muon is PF
0181         float pfPt = muon->pt();
0182         float pfPhi = muon->p4().phi();
0183         float pfEta = muon->p4().eta();
0184         float dPtOverPt = (pfPt / pt) - 1;  // TUNEP vs PF pt used as denum.
0185 
0186         if (pt < theHighPtTh) {
0187           fillInRange(getPlot(group, "code"), 2, muonTrackType(muon, false), muonTrackType(muon, true));
0188           fillInRange(getPlot(group, "deltaPtOverPt"), 1, dPtOverPt);
0189         } else {
0190           fillInRange(getPlot(group, "codeHighPt"), 2, muonTrackType(muon, false), muonTrackType(muon, true));
0191           fillInRange(getPlot(group, "deltaPtOverPtHighPt"), 1, dPtOverPt);
0192         }
0193 
0194         fillInRange(getPlot(group, "deltaPt"), 1, pfPt - pt);
0195         fillInRange(getPlot(group, "deltaPhi"), 1, fDeltaPhi(pfPhi, phi));
0196         fillInRange(getPlot(group, "deltaEta"), 1, pfEta - eta);
0197 
0198         if (theRunOnMC && genMuon)  // has a matched gen muon
0199 
0200         {
0201           float genPt = genMuon->pt();
0202           float dPtOverPtGen = (pt / genPt) - 1;
0203           float dPtOverPtGenPF = (pfPt / genPt) - 1;
0204 
0205           if (pt < theHighPtTh) {
0206             fillInRange(getPlot(group, "deltaPtOverPtPFvsTUNEP"), 2, dPtOverPtGen, dPtOverPtGenPF);
0207           } else {
0208             fillInRange(getPlot(group, "deltaPtOverPtHighPtPFvsTUNEP"), 2, dPtOverPtGen, dPtOverPtGenPF);
0209           }
0210         }
0211       }
0212     }
0213   }
0214 }
0215 
0216 void MuonPFAnalyzer::bookHistos(DQMStore::IBooker &ibooker, const string &group) {
0217   LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] Booking histos for group :" << group << "\n";
0218 
0219   ibooker.setCurrentFolder(string(theFolder) + group);
0220 
0221   bool isPFvsTUNEP = group.find("PFvsTUNEP") != string::npos;
0222 
0223   string hName;
0224 
0225   hName = "deltaPtOverPt" + group;
0226   thePlots[group]["deltaPtOverPt"] = ibooker.book1D(hName.c_str(), hName.c_str(), 101, -1.01, 1.01);
0227 
0228   hName = "deltaPtOverPtHighPt" + group;
0229   thePlots[group]["deltaPtOverPtHighPt"] = ibooker.book1D(hName.c_str(), hName.c_str(), 101, -1.01, 1.01);
0230 
0231   hName = "deltaPt" + group;
0232   thePlots[group]["deltaPt"] = ibooker.book1D(hName.c_str(), hName.c_str(), 201., -10.25, 10.25);
0233 
0234   hName = "deltaPhi" + group;
0235   thePlots[group]["deltaPhi"] = ibooker.book1D(hName.c_str(), hName.c_str(), 51., 0, .0102);
0236 
0237   hName = "deltaEta" + group;
0238   thePlots[group]["deltaEta"] = ibooker.book1D(hName.c_str(), hName.c_str(), 101., -.00505, .00505);
0239 
0240   if (isPFvsTUNEP) {
0241     hName = "code" + group;
0242     MonitorElement *plot = ibooker.book2D(hName.c_str(), hName.c_str(), 7, -.5, 6.5, 7, -.5, 6.5);
0243     thePlots[group]["code"] = plot;
0244     setCodeLabels(plot, 1);
0245     setCodeLabels(plot, 2);
0246 
0247     hName = "codeHighPt" + group;
0248     plot = ibooker.book2D(hName.c_str(), hName.c_str(), 7, -.5, 6.5, 7, -.5, 6.5);
0249     thePlots[group]["codeHighPt"] = plot;
0250     setCodeLabels(plot, 1);
0251     setCodeLabels(plot, 2);
0252 
0253     if (theRunOnMC) {
0254       hName = "deltaPtOverPtPFvsTUNEP" + group;
0255       thePlots[group]["deltaPtOverPtPFvsTUNEP"] =
0256           ibooker.book2D(hName.c_str(), hName.c_str(), 101, -1.01, 1.01, 101, -1.01, 1.01);
0257 
0258       hName = "deltaPtOverPtHighPtPFvsTUNEP" + group;
0259       thePlots[group]["deltaPtOverPtHighPtPFvsTUNEP"] =
0260           ibooker.book2D(hName.c_str(), hName.c_str(), 101, -1.01, 1.01, 101, -1.01, 1.01);
0261     }
0262   } else {
0263     hName = "code" + group;
0264     MonitorElement *plot = ibooker.book1D(hName.c_str(), hName.c_str(), 7, -.5, 6.5);
0265     thePlots[group]["code"] = plot;
0266     setCodeLabels(plot, 1);
0267 
0268     hName = "codeHighPt" + group;
0269     plot = ibooker.book1D(hName.c_str(), hName.c_str(), 7, -.5, 6.5);
0270     thePlots[group]["codeHighPt"] = plot;
0271     setCodeLabels(plot, 1);
0272   }
0273 }
0274 
0275 MuonPFAnalyzer::MonitorElement *MuonPFAnalyzer::getPlot(const string &group, const string &type) {
0276   map<string, map<string, MonitorElement *> >::iterator groupIt = thePlots.find(group);
0277   if (groupIt == thePlots.end()) {
0278     LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] GROUP : " << group << " is not a valid plot group. Returning 0.\n";
0279     return nullptr;
0280   }
0281 
0282   map<string, MonitorElement *>::iterator typeIt = groupIt->second.find(type);
0283   if (typeIt == groupIt->second.end()) {
0284     LogTrace("MuonPFAnalyzer") << "[MuonPFAnalyzer] TYPE : " << type << " is not a valid type for GROUP : " << group
0285                                << ". Returning 0.\n";
0286     return nullptr;
0287   }
0288 
0289   return typeIt->second;
0290 }
0291 
0292 inline float MuonPFAnalyzer::combRelIso(const reco::Muon *muon) {
0293   MuonIsolation iso = muon->isolationR03();
0294   float combRelIso = (iso.emEt + iso.hadEt + iso.sumPt) / muon->pt();
0295 
0296   return combRelIso;
0297 }
0298 
0299 inline float MuonPFAnalyzer::fDeltaPhi(float phi1, float phi2) {
0300   float fPhiDiff = fabs(acos(cos(phi1 - phi2)));
0301   return fPhiDiff;
0302 }
0303 
0304 void MuonPFAnalyzer::setCodeLabels(MonitorElement *plot, int nAxis) {
0305   TAxis *axis = nullptr;
0306 
0307   TH1 *histo = plot->getTH1();
0308   if (!histo)
0309     return;
0310 
0311   if (nAxis == 1)
0312     axis = histo->GetXaxis();
0313   else if (nAxis == 2)
0314     axis = histo->GetYaxis();
0315 
0316   if (!axis)
0317     return;
0318 
0319   axis->SetBinLabel(1, "Inner Track");
0320   axis->SetBinLabel(2, "Outer Track");
0321   axis->SetBinLabel(3, "Combined");
0322   axis->SetBinLabel(4, "TPFMS");
0323   axis->SetBinLabel(5, "Picky");
0324   axis->SetBinLabel(6, "DYT");
0325   axis->SetBinLabel(7, "None");
0326 }
0327 
0328 void MuonPFAnalyzer::fillInRange(MonitorElement *plot, int nAxis, double x, double y) {
0329   TH1 *histo = plot->getTH1();
0330 
0331   TAxis *axis[2] = {nullptr, nullptr};
0332   axis[0] = histo->GetXaxis();
0333   if (nAxis == 2)
0334     axis[1] = histo->GetYaxis();
0335 
0336   double value[2] = {0, 0};
0337   value[0] = x;
0338   value[1] = y;
0339 
0340   for (int i = 0; i < nAxis; ++i) {
0341     double min = axis[i]->GetXmin();
0342     double max = axis[i]->GetXmax();
0343 
0344     if (value[i] <= min)
0345       value[i] = axis[i]->GetBinCenter(1);
0346 
0347     if (value[i] >= max)
0348       value[i] = axis[i]->GetBinCenter(axis[i]->GetNbins());
0349   }
0350 
0351   if (nAxis == 2)
0352     plot->Fill(value[0], value[1]);
0353   else
0354     plot->Fill(value[0]);
0355 }
0356 
0357 int MuonPFAnalyzer::muonTrackType(const Muon *muon, bool usePF) {
0358   switch (usePF ? muon->muonBestTrackType() : muon->tunePMuonBestTrackType()) {
0359     case Muon::InnerTrack:
0360       return 0;
0361     case Muon::OuterTrack:
0362       return 1;
0363     case Muon::CombinedTrack:
0364       return 2;
0365     case Muon::TPFMS:
0366       return 3;
0367     case Muon::Picky:
0368       return 4;
0369     case Muon::DYT:
0370       return 5;
0371     case Muon::None:
0372       return 6;
0373   }
0374 
0375   return 6;
0376 }
0377 
0378 void MuonPFAnalyzer::recoToGenMatch(Handle<MuonCollection> &muons, Handle<GenParticleCollection> &gens) {
0379   theRecoGen.clear();
0380 
0381   if (muons.isValid()) {
0382     MuonCollection::const_iterator muonIt = muons->begin();
0383     MuonCollection::const_iterator muonEnd = muons->end();
0384 
0385     for (; muonIt != muonEnd; ++muonIt) {
0386       float bestDR = 999.;
0387       const GenParticle *bestGen = nullptr;
0388 
0389       if (theRunOnMC && gens.isValid()) {
0390         GenParticleCollection::const_iterator genIt = gens->begin();
0391         GenParticleCollection::const_iterator genEnd = gens->end();
0392 
0393         for (; genIt != genEnd; ++genIt) {
0394           if (abs(genIt->pdgId()) == 13) {
0395             float muonPhi = muonIt->phi();
0396             float muonEta = muonIt->eta();
0397 
0398             float genPhi = genIt->phi();
0399             float genEta = genIt->eta();
0400 
0401             float dR = deltaR(muonEta, muonPhi, genEta, genPhi);
0402 
0403             if (dR < theRecoGenR && dR < bestDR) {
0404               bestDR = dR;
0405               bestGen = &(*genIt);
0406             }
0407           }
0408         }
0409       }
0410 
0411       theRecoGen.push_back(RecoGenPair(&(*muonIt), bestGen));
0412     }
0413   }
0414 }
0415 
0416 const reco::Vertex MuonPFAnalyzer::getPrimaryVertex(Handle<VertexCollection> &vertex, Handle<BeamSpot> &beamSpot) {
0417   Vertex::Point posVtx;
0418   Vertex::Error errVtx;
0419 
0420   bool hasPrimaryVertex = false;
0421 
0422   if (vertex.isValid()) {
0423     vector<Vertex>::const_iterator vertexIt = vertex->begin();
0424     vector<Vertex>::const_iterator vertexEnd = vertex->end();
0425 
0426     for (; vertexIt != vertexEnd; ++vertexIt) {
0427       if (vertexIt->isValid() && !vertexIt->isFake()) {
0428         posVtx = vertexIt->position();
0429         errVtx = vertexIt->error();
0430         hasPrimaryVertex = true;
0431         break;
0432       }
0433     }
0434   }
0435 
0436   if (!hasPrimaryVertex) {
0437     LogInfo("MuonPFAnalyzer") << "[MuonPFAnalyzer] PrimaryVertex not found, use BeamSpot position instead.\n";
0438 
0439     posVtx = beamSpot->position();
0440     errVtx(0, 0) = beamSpot->BeamWidthX();
0441     errVtx(1, 1) = beamSpot->BeamWidthY();
0442     errVtx(2, 2) = beamSpot->sigmaZ();
0443   }
0444 
0445   const Vertex primaryVertex(posVtx, errVtx);
0446 
0447   return primaryVertex;
0448 }
0449 
0450 //define this as a plug-in
0451 DEFINE_FWK_MODULE(MuonPFAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!