PatExamples/test/MuonVIDCITKAnalyzer.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "DataFormats/Common/interface/ValueMap.h"
0008 #include "DataFormats/PatCandidates/interface/Muon.h"
0009 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0010 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0011 #include "DataFormats/VertexReco/interface/Vertex.h"
0012 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include "TH1F.h"
0018 #include "TH2F.h"
0019
0020 #include <sstream>
0021
0022 using namespace std;
0023 using namespace edm;
0024 //
0025 // class decleration
0026 //
0027 class MuonVIDCITKAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0028 public:
0029   MuonVIDCITKAnalyzer(const edm::ParameterSet&);
0030   void analyze(const edm::Event&, const edm::EventSetup&) override;
0031
0032 private:
0033   edm::EDGetTokenT<pat::MuonCollection> muonToken_;
0034   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0035
0036   // VID
0037   typedef edm::ValueMap<bool> VIDMap;
0038   edm::EDGetTokenT<VIDMap> muonLooseVIDToken_;
0039   edm::EDGetTokenT<VIDMap> muonMediumVIDToken_;
0040   edm::EDGetTokenT<VIDMap> muonTightVIDToken_;
0041   edm::EDGetTokenT<VIDMap> muonSoftVIDToken_;
0042   edm::EDGetTokenT<VIDMap> muonHighPtVIDToken_;
0043
0044   // CITK
0045   typedef edm::ValueMap<float> CITKMap;
0046   edm::EDGetTokenT<CITKMap> muonChIsoCITKToken_;
0047   edm::EDGetTokenT<CITKMap> muonNhIsoCITKToken_;
0048   edm::EDGetTokenT<CITKMap> muonPhIsoCITKToken_;
0049   edm::EDGetTokenT<CITKMap> muonPuIsoCITKToken_;
0050
0051   // Histograms
0052   TH2F *h2Ch_, *h2Nh_, *h2Ph_, *h2Pu_;
0053   TH1F *hIsoDiffCh_, *hIsoDiffNh_, *hIsoDiffPh_, *hIsoDiffPu_;
0054 };
0055
0056 MuonVIDCITKAnalyzer::MuonVIDCITKAnalyzer(const edm::ParameterSet& iConfig) {
0057   muonToken_ = consumes<pat::MuonCollection>(iConfig.getParameter<edm::InputTag>("muon"));
0058   vertexToken_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertex"));
0059
0060   string vidPrefix = "muoMuonIDs:cutBasedMuonId-MuonPOG-V0-";
0061   muonLooseVIDToken_ = consumes<VIDMap>(edm::InputTag(vidPrefix + "loose"));
0062   muonMediumVIDToken_ = consumes<VIDMap>(edm::InputTag(vidPrefix + "medium"));
0063   muonTightVIDToken_ = consumes<VIDMap>(edm::InputTag(vidPrefix + "tight"));
0064   muonSoftVIDToken_ = consumes<VIDMap>(edm::InputTag(vidPrefix + "soft"));
0065   muonHighPtVIDToken_ = consumes<VIDMap>(edm::InputTag(vidPrefix + "highpt"));
0066
0067   muonChIsoCITKToken_ =
0068       consumes<CITKMap>(edm::InputTag("muonPFNoPileUpIsolation:h+-DR040-ThresholdVeto000-ConeVeto000"));
0069   muonNhIsoCITKToken_ =
0070       consumes<CITKMap>(edm::InputTag("muonPFNoPileUpIsolation:h0-DR040-ThresholdVeto050-ConeVeto001"));
0071   muonPhIsoCITKToken_ =
0072       consumes<CITKMap>(edm::InputTag("muonPFNoPileUpIsolation:gamma-DR040-ThresholdVeto050-ConeVeto001"));
0073   muonPuIsoCITKToken_ = consumes<CITKMap>(edm::InputTag("muonPFPileUpIsolation:h+-DR040-ThresholdVeto050-ConeVeto001"));
0074
0075   usesResource("TFileService");
0076   edm::Service<TFileService> fs;
0077   h2Ch_ = fs->make<TH2F>("h2Ch", "ChIso;IsoValue ChIso;CITK ChIso", 100, 0, 10, 100, 0, 10);
0078   h2Nh_ = fs->make<TH2F>("h2Nh", "NhIso;IsoValue NhIso;CITK NhIso", 100, 0, 10, 100, 0, 10);
0079   h2Ph_ = fs->make<TH2F>("h2Ph", "PhIso;IsoValue PhIso;CITK PhIso", 100, 0, 10, 100, 0, 10);
0080   h2Pu_ = fs->make<TH2F>("h2Pu", "PuIso;IsoValue PuIso;CITK PuIso", 100, 0, 10, 100, 0, 10);
0081
0082   hIsoDiffCh_ = fs->make<TH1F>("hIsoDiffCh", "Diff. ChIso;CITK-IsoValue", 1000, -100, 100);
0083   hIsoDiffNh_ = fs->make<TH1F>("hIsoDiffNh", "Diff. NhIso;CITK-IsoValue", 1000, -100, 100);
0084   hIsoDiffPh_ = fs->make<TH1F>("hIsoDiffPh", "Diff. PhIso;CITK-IsoValue", 1000, -100, 100);
0085   hIsoDiffPu_ = fs->make<TH1F>("hIsoDiffPu", "Diff. PuIso;CITK-IsoValue", 1000, -100, 100);
0086 }
0087
0088 void MuonVIDCITKAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& eSetup) {
0089   edm::Handle<pat::MuonCollection> muonHandle;
0090   event.getByToken(muonToken_, muonHandle);
0091
0092   edm::Handle<reco::VertexCollection> vertexHandle;
0093   event.getByToken(vertexToken_, vertexHandle);
0094   reco::Vertex vertex = vertexHandle->at(0);
0095
0096   // VID
0097   edm::Handle<VIDMap> muonLooseVIDHandle;
0098   edm::Handle<VIDMap> muonMediumVIDHandle;
0099   edm::Handle<VIDMap> muonTightVIDHandle;
0100   edm::Handle<VIDMap> muonSoftVIDHandle;
0101   edm::Handle<VIDMap> muonHighPtVIDHandle;
0102   event.getByToken(muonLooseVIDToken_, muonLooseVIDHandle);
0103   event.getByToken(muonMediumVIDToken_, muonMediumVIDHandle);
0104   event.getByToken(muonTightVIDToken_, muonTightVIDHandle);
0105   event.getByToken(muonSoftVIDToken_, muonSoftVIDHandle);
0106   event.getByToken(muonHighPtVIDToken_, muonHighPtVIDHandle);
0107
0108   // CITK
0109   edm::Handle<CITKMap> muonChIsoCITKHandle;
0110   edm::Handle<CITKMap> muonNhIsoCITKHandle;
0111   edm::Handle<CITKMap> muonPhIsoCITKHandle;
0112   edm::Handle<CITKMap> muonPuIsoCITKHandle;
0113   event.getByToken(muonChIsoCITKToken_, muonChIsoCITKHandle);
0114   event.getByToken(muonNhIsoCITKToken_, muonNhIsoCITKHandle);
0115   event.getByToken(muonPhIsoCITKToken_, muonPhIsoCITKHandle);
0116   event.getByToken(muonPuIsoCITKToken_, muonPuIsoCITKHandle);
0117
0118   for (size_t i = 0, n = muonHandle->size(); i < n; ++i) {
0119     const auto& mu = dynamic_cast<const pat::Muon&>(muonHandle->at(i));
0120     const auto& muRef = pat::MuonRef(muonHandle, i);
0121
0122     stringstream sout;
0123
0124     // Check standard ID vs VID
0125     const bool vidLoose = (*muonLooseVIDHandle)[muRef];
0126     const bool vidMedium = (*muonMediumVIDHandle)[muRef];
0127     const bool vidTight = (*muonTightVIDHandle)[muRef];
0128     const bool vidSoft = (*muonSoftVIDHandle)[muRef];
0129     const bool vidHighPt = (*muonHighPtVIDHandle)[muRef];
0130
0131     const bool isLoose = muon::isLooseMuon(mu);
0132     const bool isMedium = muon::isMediumMuon(mu);
0133     const bool isTight = muon::isTightMuon(mu, vertex);
0134     const bool isSoft = muon::isSoftMuon(mu, vertex);
0135     const bool isHighPt = muon::isHighPtMuon(mu, vertex);
0136
0137     if (vidLoose != isLoose) {
0138       sout << " isLoose " << vidLoose << ' ' << isLoose << endl;
0139     }
0140     if (vidMedium != isMedium) {
0141       sout << " isMedium" << vidMedium << ' ' << isMedium << endl;
0142     }
0143     if (vidTight != isTight) {
0144       sout << " isTight " << vidTight << ' ' << isTight << endl;
0145     }
0146     if (vidSoft != isSoft) {
0147       sout << " isSoft  " << vidSoft << ' ' << isSoft << endl;
0148     }
0149     if (vidHighPt != isHighPt) {
0150       sout << " isHighPt" << vidHighPt << ' ' << isHighPt << endl;
0151     }
0152
0153     // Check standard IsoDeposit vs CITK
0154     const double citkChIso = (*muonChIsoCITKHandle)[muRef];
0155     const double citkNhIso = (*muonNhIsoCITKHandle)[muRef];
0156     const double citkPhIso = (*muonPhIsoCITKHandle)[muRef];
0157     const double citkPuIso = (*muonPuIsoCITKHandle)[muRef];
0158
0159     const double patChIso = mu.chargedHadronIso();
0160     const double patNhIso = mu.neutralHadronIso();
0161     const double patPhIso = mu.photonIso();
0162     const double patPuIso = mu.puChargedHadronIso();
0163
0164     h2Ch_->Fill(patChIso, citkChIso);
0165     h2Nh_->Fill(patNhIso, citkNhIso);
0166     h2Ph_->Fill(patPhIso, citkPhIso);
0167     h2Pu_->Fill(patPuIso, citkPuIso);
0168
0169     hIsoDiffCh_->Fill(citkChIso - patChIso);
0170     hIsoDiffNh_->Fill(citkNhIso - patNhIso);
0171     hIsoDiffPh_->Fill(citkPhIso - patPhIso);
0172     hIsoDiffPu_->Fill(citkPuIso - patPuIso);
0173
0174     if (std::abs(citkChIso - patChIso) >= 1e-4) {
0175       sout << " ChIso citk=" << citkChIso << " isodep=" << patChIso << endl;
0176     }
0177     if (std::abs(citkNhIso - patNhIso) >= 1e-4) {
0178       sout << " NhIso citk=" << citkNhIso << " isodep=" << patNhIso << endl;
0179     }
0180     if (std::abs(citkPhIso - patPhIso) >= 1e-4) {
0181       sout << " PhIso citk=" << citkPhIso << " isodep=" << patPhIso << endl;
0182     }
0183     if (std::abs(citkPuIso - patPuIso) >= 1e-4) {
0184       sout << " PuIso citk=" << citkPuIso << " isodep=" << patPuIso << endl;
0185     }
0186
0187     if (!sout.str().empty()) {
0188       cout << event.id().event() << " mu" << i << " pt=" << muRef->pt() << " eta=" << muRef->eta()
0189            << " phi=" << muRef->phi() << endl;
0190       cout << sout.rdbuf();
0191     }
0192   }
0193 }
0194
0195 //define this as a plug-in
0196 #include "FWCore/Framework/interface/MakerMacros.h"
0197 DEFINE_FWK_MODULE(MuonVIDCITKAnalyzer);