Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:46

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 
0167     if (!pv.isFake() && pv.ndof() >= 4 && fabs(pv.z()) <= 24.0 && fabs(pv.position().rho()) <= 2.0) {
0168       npv++;
0169       isGoodPV[i] = true;
0170     }
0171   }
0172   th1dPlots["npv"].fill(npv);
0173   int npv_in_range = npv;
0174   if (npv_in_range < 0)
0175     npv_in_range = 0;
0176   else if (npv_in_range >= npvHigh)
0177     npv_in_range = npvHigh - 1;  // make sure int_mu won't lead to non-existing ME
0178 
0179   //mu//
0180   int int_mu = -1;
0181   edm::Handle<edm::View<PileupSummaryInfo>> muHandle;
0182   if (iEvent.getByToken(muToken, muHandle)) {
0183     const auto& summary = *muHandle;
0184     auto it = std::find_if(summary.begin(), summary.end(), [](const auto& s) { return s.getBunchCrossing() == 0; });
0185 
0186     if (it->getBunchCrossing() != 0) {
0187       edm::LogError("OffsetAnalyzerDQM") << "Cannot find the in-time pileup info " << it->getBunchCrossing();
0188     } else {
0189       float mu = it->getTrueNumInteractions();
0190       th1dPlots["mu"].fill(mu);
0191       int_mu = mu + 0.5;
0192     }
0193   }
0194   if (int_mu >= muHigh)
0195     int_mu = muHigh - 1;  // make sure int_mu won't lead to non-existing ME
0196 
0197   //create map of pftypes vs total energy / eta
0198   std::map<std::string, std::vector<double>> m_pftype_etaE;
0199   int nEta = etabins.size() - 1;
0200   for (const auto& pftype : pftypes)
0201     m_pftype_etaE[pftype].assign(nEta, 0.0);
0202 
0203   //pf particles//
0204   edm::Handle<edm::View<pat::PackedCandidate>> pfHandle;
0205   iEvent.getByToken(pfToken, pfHandle);
0206 
0207   for (unsigned int i = 0, n = pfHandle->size(); i < n; i++) {
0208     const auto& cand = pfHandle->at(i);
0209 
0210     int etaIndex = getEtaIndex(cand.eta());
0211     std::string pftype = pdgMap[abs(cand.pdgId())];
0212     if (etaIndex == -1 || pftype.empty())
0213       continue;
0214 
0215     if (pftype == "chm") {  //check charged hadrons ONLY
0216       bool attached = false;
0217 
0218       for (unsigned int ipv = 0; ipv < nPVall && !attached; ipv++) {
0219         if (isGoodPV[ipv] && cand.fromPV(ipv) == 3)
0220           attached = true;  //pv used in fit
0221       }
0222       if (!attached)
0223         pftype = "chu";  //unmatched charged hadron
0224     }
0225     ////AOD////
0226     /*
0227         reco::TrackRef candTrkRef( cand.trackRef() );
0228         if ( pftype == "chm" && !candTrkRef.isNull() ) { //check charged hadrons ONLY
0229             bool attached = false;
0230 
0231             for (auto ipv=vertexHandle->begin(), endpv=vertexHandle->end(); ipv != endpv && !attached; ++ipv) {
0232                 if ( !ipv->isFake() && ipv->ndof() >= 4 && fabs(ipv->z()) < 24 ) { //must be attached to a good pv
0233 
0234                     for(auto ivtrk=ipv->tracks_begin(), endvtrk=ipv->tracks_end(); ivtrk != endvtrk && !attached; ++ivtrk) {
0235                         reco::TrackRef pvTrkRef(ivtrk->castTo<reco::TrackRef>());
0236                         if (pvTrkRef == candTrkRef) attached = true;
0237                     }
0238                 }
0239             }
0240             if (!attached) pftype = "chu"; //unmatched charged hadron
0241         }
0242 */
0243     ///////////
0244     m_pftype_etaE[pftype][etaIndex] += cand.et();
0245   }
0246 
0247   for (const auto& pair : m_pftype_etaE) {
0248     std::string pftype = pair.first;
0249     std::vector<double> etaE = pair.second;
0250 
0251     std::string offset_name_npv = offsetPlotBaseName + "_npv" + std::to_string(npv_in_range) + "_" + pftype;
0252     if (offsetPlots.find(offset_name_npv) == offsetPlots.end())
0253       return;  //npv is out of range ()
0254 
0255     for (int i = 0; i < nEta; i++) {
0256       double eta = 0.5 * (etabins[i] + etabins[i + 1]);
0257       offsetPlots[offset_name_npv].fill2D(eta, etaE[i]);
0258     }
0259 
0260     if (int_mu != -1) {
0261       std::string offset_name_mu = offsetPlotBaseName + "_mu" + std::to_string(int_mu) + "_" + pftype;
0262       if (offsetPlots.find(offset_name_mu) == offsetPlots.end())
0263         return;  //mu is out of range
0264 
0265       for (int i = 0; i < nEta; i++) {
0266         double eta = 0.5 * (etabins[i] + etabins[i + 1]);
0267         offsetPlots[offset_name_mu].fill2D(eta, etaE[i]);
0268       }
0269     }
0270   }
0271 }
0272 
0273 int OffsetAnalyzerDQM::getEtaIndex(float eta) {
0274   int nEta = etabins.size() - 1;
0275 
0276   for (int i = 0; i < nEta; i++) {
0277     if (etabins[i] <= eta && eta < etabins[i + 1])
0278       return i;
0279   }
0280   if (eta == etabins[nEta])
0281     return nEta - 1;
0282   else
0283     return -1;
0284 }
0285 
0286 #include "FWCore/Framework/interface/MakerMacros.h"
0287 DEFINE_FWK_MODULE(OffsetAnalyzerDQM);