RecoParticleFlow/plugins/OffsetAnalyzerDQM.cc

0001 #include "DQMOffline/PFTau/interface/Matchers.h"
0002 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0003 #include "DQMServices/Core/interface/DQMStore.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0011 #include <SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h>
0012
0013 #include <vector>
0014 #include <map>
0015
0016 class OffsetAnalyzerDQM : public DQMEDAnalyzer {
0017 public:
0018   OffsetAnalyzerDQM(const edm::ParameterSet&);
0019   void analyze(const edm::Event&, const edm::EventSetup&) override;
0020
0021 protected:
0022   //Book histograms
0023   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0024   int getEtaIndex(float eta);
0025
0026 private:
0027   struct Plot1D {
0028     std::string name, title, dir;
0029     int nxbins;
0030     double xlow, xhigh;
0031     MonitorElement* plot;
0032
0033     Plot1D() : name(""), title(""), dir(""), nxbins(0), xlow(0), xhigh(0), plot(nullptr) {}
0034     Plot1D(const std::string& n, const std::string& t, const std::string& d, int nx, double x0, double x1)
0035         : name(n), title(t), dir(d), nxbins(nx), xlow(x0), xhigh(x1), plot(nullptr) {}
0036
0037     virtual void book(DQMStore::IBooker& booker) {
0038       booker.setCurrentFolder(dir);
0039       plot = booker.book1D(name, title, nxbins, xlow, xhigh);
0040       plot->setStatOverflows(kTRUE);
0041     }
0042
0043     virtual void fill(float value) {
0044       assert(plot != nullptr);
0045       plot->Fill(value);
0046     }
0047
0048     virtual ~Plot1D() {}
0049   };
0050
0051   struct PlotProfile : public Plot1D {
0052     std::vector<double> xbins;
0053     int nybins;
0054     double ylow, yhigh;
0055     PlotProfile() : Plot1D(), xbins(0), nybins(0), ylow(0), yhigh(0) {}
0056     PlotProfile(const std::string& n,
0057                 const std::string& t,
0058                 const std::string& d,
0059                 int nx,
0060                 double x0,
0061                 double x1,
0062                 const std::vector<double>& vx,
0063                 int ny,
0064                 double y0,
0065                 double y1)
0066         : Plot1D(n, t, d, nx, x0, x1), xbins(vx), nybins(ny), ylow(y0), yhigh(y1) {}
0067
0068     void book(DQMStore::IBooker& booker) override {
0069       booker.setCurrentFolder(dir);
0070       plot = booker.bookProfile(name, title, xbins.size() - 1, &xbins[0], nybins, ylow, yhigh, " ");
0071     }
0072     //make other booker methods for uniform binning
0073
0074     void fill2D(double value1, double value2) {
0075       assert(plot != nullptr);
0076       plot->Fill(value1, value2);
0077     }
0078   };
0079
0080   std::map<std::string, Plot1D> th1dPlots;
0081   std::map<std::string, PlotProfile> offsetPlots;
0082   std::map<int, std::string> pdgMap;
0083
0084   std::string offsetPlotBaseName;
0085   std::vector<std::string> pftypes;
0086   std::vector<double> etabins;
0087
0088   int muHigh;
0089   int npvHigh;
0090
0091   edm::EDGetTokenT<edm::View<reco::Vertex>> pvToken;
0092   edm::EDGetTokenT<edm::View<PileupSummaryInfo>> muToken;
0093   edm::EDGetTokenT<edm::View<pat::PackedCandidate>> pfToken;
0094 };
0095
0096 OffsetAnalyzerDQM::OffsetAnalyzerDQM(const edm::ParameterSet& iConfig) {
0097   offsetPlotBaseName = iConfig.getParameter<std::string>("offsetPlotBaseName");
0098
0099   pvToken = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvTag"));
0100   muToken = consumes<edm::View<PileupSummaryInfo>>(iConfig.getParameter<edm::InputTag>("muTag"));
0101   pfToken = consumes<edm::View<pat::PackedCandidate>>(iConfig.getParameter<edm::InputTag>("pfTag"));
0102
0103   etabins = iConfig.getParameter<std::vector<double>>("etabins");
0104   pftypes = iConfig.getParameter<std::vector<std::string>>("pftypes");
0105
0106   muHigh = iConfig.getUntrackedParameter<int>("muHigh");
0107   npvHigh = iConfig.getUntrackedParameter<int>("npvHigh");
0108
0109   //initialize offset plots
0110   const auto& offset_psets = iConfig.getParameter<std::vector<edm::ParameterSet>>("offsetPlots");
0111   for (auto& pset : offset_psets) {
0112     std::string name = pset.getParameter<std::string>("name");
0113     std::string title = pset.getParameter<std::string>("title");
0114     std::string dir = pset.getParameter<std::string>("dir");
0115     std::vector<double> vx = pset.getParameter<std::vector<double>>("vx");
0116     int ny = pset.getParameter<uint32_t>("ny");
0117     double y0 = pset.getParameter<double>("y0");
0118     double y1 = pset.getParameter<double>("y1");
0119
0120     offsetPlots[name] = PlotProfile(name, title, dir, 0, 0, 0, vx, ny, y0, y1);
0121   }
0122
0123   //initialize th1d
0124   const auto& th1d_psets = iConfig.getParameter<std::vector<edm::ParameterSet>>("th1dPlots");
0125   for (auto& pset : th1d_psets) {
0126     std::string name = pset.getParameter<std::string>("name");
0127     std::string title = pset.getParameter<std::string>("title");
0128     std::string dir = pset.getParameter<std::string>("dir");
0129     int nx = pset.getParameter<uint32_t>("nx");
0130     double x0 = pset.getParameter<double>("x0");
0131     double x1 = pset.getParameter<double>("x1");
0132
0133     th1dPlots[name] = Plot1D(name, title, dir, nx, x0, x1);
0134   }
0135
0136   //create pdg map
0137   std::vector<uint32_t> pdgKeys = iConfig.getParameter<std::vector<uint32_t>>("pdgKeys");
0138   std::vector<std::string> pdgStrs = iConfig.getParameter<std::vector<std::string>>("pdgStrs");
0139   for (int i = 0, n = pdgKeys.size(); i < n; i++)
0140     pdgMap[pdgKeys[i]] = pdgStrs[i];
0141 }
0142
0143 void OffsetAnalyzerDQM::bookHistograms(DQMStore::IBooker& booker, edm::Run const&, edm::EventSetup const&) {
0144   //std::cout << "OffsetAnalyzerDQM booking offset histograms" << std::endl;
0145   for (auto& pair : offsetPlots) {
0146     pair.second.book(booker);
0147   }
0148   for (auto& pair : th1dPlots) {
0149     pair.second.book(booker);
0150   }
0151 }
0152
0153 void OffsetAnalyzerDQM::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0154   //npv//
0155   edm::Handle<edm::View<reco::Vertex>> vertexHandle;
0156   iEvent.getByToken(pvToken, vertexHandle);
0157
0158   unsigned int nPVall = vertexHandle->size();
0159   bool isGoodPV[nPVall];
0160   for (size_t i = 0; i < nPVall; ++i)
0161     isGoodPV[i] = false;
0162
0163   int npv = 0;
0164   for (unsigned int i = 0; i < nPVall; i++) {
0165     const auto& pv = vertexHandle->at(i);
0166     th1dPlots["pv_z"].fill(pv.z());
0167
0168     if (!pv.isFake() && pv.ndof() >= 4 && fabs(pv.z()) <= 24.0 && fabs(pv.position().rho()) <= 2.0) {
0169       npv++;
0170       isGoodPV[i] = true;
0171     }
0172   }
0173   th1dPlots["npv"].fill(npv);
0174   int npv_in_range = npv;
0175   if (npv_in_range < 0)
0176     npv_in_range = 0;
0177   else if (npv_in_range >= npvHigh)
0178     npv_in_range = npvHigh - 1;  // make sure int_mu won't lead to non-existing ME
0179
0180   //mu//
0181   int int_mu = -1;
0182   edm::Handle<edm::View<PileupSummaryInfo>> muHandle;
0183   if (iEvent.getByToken(muToken, muHandle)) {
0184     const auto& summary = *muHandle;
0185     auto it = std::find_if(summary.begin(), summary.end(), [](const auto& s) { return s.getBunchCrossing() == 0; });
0186
0187     if (it->getBunchCrossing() != 0) {
0188       edm::LogError("OffsetAnalyzerDQM") << "Cannot find the in-time pileup info " << it->getBunchCrossing();
0189     } else {
0190       float mu = it->getTrueNumInteractions();
0191       th1dPlots["mu"].fill(mu);
0192       int_mu = mu + 0.5;
0193     }
0194   }
0195   if (int_mu >= muHigh)
0196     int_mu = muHigh - 1;  // make sure int_mu won't lead to non-existing ME
0197
0198   //create map of pftypes vs total energy / eta
0199   std::map<std::string, std::vector<double>> m_pftype_etaE;
0200   int nEta = etabins.size() - 1;
0201   for (const auto& pftype : pftypes)
0202     m_pftype_etaE[pftype].assign(nEta, 0.0);
0203
0204   //pf particles//
0205   edm::Handle<edm::View<pat::PackedCandidate>> pfHandle;
0206   iEvent.getByToken(pfToken, pfHandle);
0207
0208   for (unsigned int i = 0, n = pfHandle->size(); i < n; i++) {
0209     const auto& cand = pfHandle->at(i);
0210
0211     int etaIndex = getEtaIndex(cand.eta());
0212     std::string pftype = pdgMap[abs(cand.pdgId())];
0213     if (etaIndex == -1 || pftype.empty())
0214       continue;
0215
0216     if (pftype == "chm") {  //check charged hadrons ONLY
0217       bool attached = false;
0218
0219       for (unsigned int ipv = 0; ipv < nPVall && !attached; ipv++) {
0220         if (isGoodPV[ipv] && cand.fromPV(ipv) == 3)
0221           attached = true;  //pv used in fit
0222       }
0223       if (!attached)
0224         pftype = "chu";  //unmatched charged hadron
0225     }
0226     ////AOD////
0227     /*
0228         reco::TrackRef candTrkRef( cand.trackRef() );
0229         if ( pftype == "chm" && !candTrkRef.isNull() ) { //check charged hadrons ONLY
0230             bool attached = false;
0231
0232             for (auto ipv=vertexHandle->begin(), endpv=vertexHandle->end(); ipv != endpv && !attached; ++ipv) {
0233                 if ( !ipv->isFake() && ipv->ndof() >= 4 && fabs(ipv->z()) < 24 ) { //must be attached to a good pv
0234
0235                     for(auto ivtrk=ipv->tracks_begin(), endvtrk=ipv->tracks_end(); ivtrk != endvtrk && !attached; ++ivtrk) {
0236                         reco::TrackRef pvTrkRef(ivtrk->castTo<reco::TrackRef>());
0237                         if (pvTrkRef == candTrkRef) attached = true;
0238                     }
0239                 }
0240             }
0241             if (!attached) pftype = "chu"; //unmatched charged hadron
0242         }
0243 */
0244     ///////////
0245     m_pftype_etaE[pftype][etaIndex] += cand.et();
0246   }
0247
0248   for (const auto& pair : m_pftype_etaE) {
0249     std::string pftype = pair.first;
0250     std::vector<double> etaE = pair.second;
0251
0252     std::string offset_name_npv = offsetPlotBaseName + "_npv" + std::to_string(npv_in_range) + "_" + pftype;
0253     if (offsetPlots.find(offset_name_npv) == offsetPlots.end())
0254       return;  //npv is out of range ()
0255
0256     for (int i = 0; i < nEta; i++) {
0257       double eta = 0.5 * (etabins[i] + etabins[i + 1]);
0258       offsetPlots[offset_name_npv].fill2D(eta, etaE[i]);
0259     }
0260
0261     if (int_mu != -1) {
0262       std::string offset_name_mu = offsetPlotBaseName + "_mu" + std::to_string(int_mu) + "_" + pftype;
0263       if (offsetPlots.find(offset_name_mu) == offsetPlots.end())
0264         return;  //mu is out of range
0265
0266       for (int i = 0; i < nEta; i++) {
0267         double eta = 0.5 * (etabins[i] + etabins[i + 1]);
0268         offsetPlots[offset_name_mu].fill2D(eta, etaE[i]);
0269       }
0270     }
0271   }
0272 }
0273
0274 int OffsetAnalyzerDQM::getEtaIndex(float eta) {
0275   int nEta = etabins.size() - 1;
0276
0277   for (int i = 0; i < nEta; i++) {
0278     if (etabins[i] <= eta && eta < etabins[i + 1])
0279       return i;
0280   }
0281   if (eta == etabins[nEta])
0282     return nEta - 1;
0283   else
0284     return -1;
0285 }
0286
0287 #include "FWCore/Framework/interface/MakerMacros.h"
0288 DEFINE_FWK_MODULE(OffsetAnalyzerDQM);