PFTau/plugins/PFClient_JetRes.cc

0001 #include "DQMOffline/PFTau/plugins/PFClient_JetRes.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/Utilities/interface/InputTag.h"
0006
0007 #include "DQMServices/Core/interface/DQMStore.h"
0008
0009 #include "TCanvas.h"
0010 #include "TGraph.h"
0011
0012 //
0013 // -- Constructor
0014 //
0015 PFClient_JetRes::PFClient_JetRes(const edm::ParameterSet &parameterSet) {
0016   folderNames_ = parameterSet.getParameter<std::vector<std::string>>("FolderNames");
0017   histogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNames");
0018   efficiencyFlag_ = parameterSet.getParameter<bool>("CreateEfficiencyPlots");
0019   effHistogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNamesForEfficiencyPlots");
0020   PtBins_ = parameterSet.getParameter<std::vector<int>>("VariablePtBins");
0021 }
0022
0023 //
0024 // -- EndJobBegin Run
0025 //
0026 void PFClient_JetRes::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0027   doSummaries(ibooker, igetter);
0028   if (efficiencyFlag_)
0029     doEfficiency(ibooker, igetter);
0030 }
0031
0032 //
0033 // -- Create Summaries
0034 //
0035 void PFClient_JetRes::doSummaries(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0036   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0037        ifolder++) {
0038     std::string path = "ParticleFlow/" + (*ifolder);
0039
0040     for (std::vector<std::string>::const_iterator ihist = histogramNames_.begin(); ihist != histogramNames_.end();
0041          ihist++) {
0042       std::string hname = (*ihist);
0043       createResolutionPlots(ibooker, igetter, path, hname);
0044     }
0045   }
0046 }
0047
0048 //
0049 // -- Create Efficiency
0050 //
0051 void PFClient_JetRes::doEfficiency(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0052   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0053        ifolder++) {
0054     std::string path = "ParticleFlow/" + (*ifolder);
0055
0056     for (std::vector<std::string>::const_iterator ihist = effHistogramNames_.begin(); ihist != effHistogramNames_.end();
0057          ihist++) {
0058       std::string hname = (*ihist);
0059       createEfficiencyPlots(ibooker, igetter, path, hname);
0060     }
0061   }
0062 }
0063
0064 //
0065 // -- Create Resolution Plots
0066 //
0067 void PFClient_JetRes::createResolutionPlots(DQMStore::IBooker &ibooker,
0068                                             DQMStore::IGetter &igetter,
0069                                             std::string &folder,
0070                                             std::string &name) {
0071   MonitorElement *me = igetter.get(folder + "/" + name);
0072   if (!me)
0073     return;
0074
0075   MonitorElement *pT[PtBins_.size() - 1];
0076   std::vector<double> pTEntries(PtBins_.size() - 1, 0);
0077
0078   // std::vector<std::string> pTRange (PtBins_.size() -1) ;
0079   // char* pTRange[PtBins_.size() -1] ;
0080   std::vector<TString> pTRange(PtBins_.size() - 1);
0081   // float pTCenter[PtBins_.size() -1] ;
0082
0083   MonitorElement *me_average;
0084   MonitorElement *me_rms;
0085   MonitorElement *me_mean;
0086   MonitorElement *me_sigma;
0087
0088   if ((me->kind() == MonitorElement::Kind::TH2F) || (me->kind() == MonitorElement::Kind::TH2S) ||
0089       (me->kind() == MonitorElement::Kind::TH2D)) {
0090     TH2 *th = me->getTH2F();
0091     // size_t nbinx = me->getNbinsX();
0092     size_t nbinx = PtBins_.size() - 1;
0093     size_t nbiny = me->getNbinsY();
0094
0095     float ymin = th->GetYaxis()->GetXmin();
0096     float ymax = th->GetYaxis()->GetXmax();
0097     std::string xtit = th->GetXaxis()->GetTitle();
0098     // std::string ytit = th->GetYaxis()->GetTitle();
0099     std::string ytit = "#Deltap_{T}/p_{T}";
0100
0101     float *xbins = new float[nbinx + 1];
0102     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0103       // xbins[ix-1] = th->GetBinLowEdge(ix);
0104       xbins[ix - 1] = PtBins_[ix - 1];
0105       // snprintf(pTRange[ix-1].data(), 15, "Pt%d_%d", PtBins_[ix-1],
0106       // PtBins_[ix]);
0107       pTRange[ix - 1] = TString::Format("Pt%d_%d", PtBins_[ix - 1], PtBins_[ix]);
0108       if (name == "BRdelta_et_Over_et_VS_et_")
0109         pTRange[ix - 1] = TString::Format("BRPt%d_%d", PtBins_[ix - 1], PtBins_[ix]);
0110       else if (name == "ERdelta_et_Over_et_VS_et_")
0111         pTRange[ix - 1] = TString::Format("ERPt%d_%d", PtBins_[ix - 1], PtBins_[ix]);
0112
0113       // pTCenter[ix-1] = (PtBins_[ix] - PtBins_[ix-1]) / 2. ;
0114       if (ix == nbinx) {
0115         // xbins[ix] = th->GetXaxis()->GetBinUpEdge(ix);
0116         xbins[ix] = PtBins_[ix];
0117       }
0118     }
0119
0120     std::string tit_new;
0121     ibooker.setCurrentFolder(folder);
0122     // MonitorElement* me_slice = ibooker.book1D("PFlowSlice", "PFlowSlice",
0123     // nbiny, ymin, ymax);
0124
0125     tit_new = "Average " + ytit + ";" + xtit + ";Average_" + ytit;
0126     me_average = ibooker.book1D("average_" + name, tit_new, nbinx, xbins);
0127     me_average->setEfficiencyFlag();
0128     tit_new = "RMS " + ytit + ";" + xtit + ";RMS_" + ytit;
0129     me_rms = ibooker.book1D("rms_" + name, tit_new, nbinx, xbins);
0130     me_rms->setEfficiencyFlag();
0131     tit_new = ";" + xtit + ";Mean_" + ytit;
0132     me_mean = ibooker.book1D("mean_" + name, tit_new, nbinx, xbins);
0133     me_mean->setEfficiencyFlag();
0134     tit_new = ";" + xtit + ";Sigma_" + ytit;
0135     me_sigma = ibooker.book1D("sigma_" + name, tit_new, nbinx, xbins);
0136     me_sigma->setEfficiencyFlag();
0137
0138     double average, rms, mean, sigma;
0139     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0140       // me_slice->Reset();
0141       if (name == "delta_et_Over_et_VS_et_")
0142         pT[ix - 1] = ibooker.book1D(
0143             pTRange[ix - 1], TString::Format("Total %s;%s;Events", ytit.data(), ytit.data()), nbiny, ymin, ymax);
0144       if (name == "BRdelta_et_Over_et_VS_et_")
0145         pT[ix - 1] = ibooker.book1D(
0146             pTRange[ix - 1], TString::Format("Barrel %s;%s;Events", ytit.data(), ytit.data()), nbiny, ymin, ymax);
0147       else if (name == "ERdelta_et_Over_et_VS_et_")
0148         pT[ix - 1] = ibooker.book1D(
0149             pTRange[ix - 1], TString::Format("Endcap %s;%s;Events", ytit.data(), ytit.data()), nbiny, ymin, ymax);
0150
0151       for (size_t iy = 0; iy <= nbiny + 1; ++iy)  // add under and overflow
0152         if (th->GetBinContent(ix, iy)) {
0153           // me_slice->setBinContent(iy,th->GetBinContent(ix,iy));
0154           pT[ix - 1]->setBinContent(iy, th->GetBinContent(ix, iy));
0155           pT[ix - 1]->setBinError(iy, th->GetBinError(ix, iy));
0156           pTEntries[ix - 1] += th->GetBinContent(ix, iy);
0157         }
0158
0159       pT[ix - 1]->setEntries(pTEntries[ix - 1]);
0160
0161       // getHistogramParameters(me_slice, average, rms, mean, sigma);
0162       getHistogramParameters(pT[ix - 1], average, rms, mean, sigma);
0163       me_average->setBinContent(ix, average);
0164       me_rms->setBinContent(ix, rms);
0165       me_mean->setBinContent(ix, mean);
0166       me_sigma->setBinContent(ix, sigma);
0167     }
0168     delete[] xbins;
0169   }
0170 }
0171
0172 //
0173 // -- Get Histogram Parameters
0174 //
0175 void PFClient_JetRes::getHistogramParameters(
0176     MonitorElement *me_slice, double &average, double &rms, double &mean, double &sigma) {
0177   average = 0.0;
0178   rms = 0.0;
0179   mean = 0.0;
0180   sigma = 0.0;
0181
0182   if (!me_slice)
0183     return;
0184   if (me_slice->kind() == MonitorElement::Kind::TH1F) {
0185     average = me_slice->getMean();
0186     rms = me_slice->getRMS();
0187     TH1F *th_slice = me_slice->getTH1F();
0188     if (th_slice && th_slice->GetEntries() > 0) {
0189       // need our own copy for thread safety
0190       TF1 gaus("mygaus", "gaus");
0191       th_slice->Fit(&gaus, "Q0 SERIAL");
0192       sigma = gaus.GetParameter(2);
0193       mean = gaus.GetParameter(1);
0194     }
0195   }
0196 }
0197
0198 //
0199 // -- Create Resolution Plots
0200 //
0201 void PFClient_JetRes::createEfficiencyPlots(DQMStore::IBooker &ibooker,
0202                                             DQMStore::IGetter &igetter,
0203                                             std::string &folder,
0204                                             std::string &name) {
0205   MonitorElement *me1 = igetter.get(folder + "/" + name);
0206   MonitorElement *me2 = igetter.get(folder + "/" + name + "ref_");
0207   if (!me1 || !me2)
0208     return;
0209   MonitorElement *me_eff;
0210   if ((me1->kind() == MonitorElement::Kind::TH1F) && (me1->kind() == MonitorElement::Kind::TH1F)) {
0211     TH1 *th1 = me1->getTH1F();
0212     size_t nbinx = me1->getNbinsX();
0213
0214     float xmin = th1->GetXaxis()->GetXmin();
0215     float xmax = th1->GetXaxis()->GetXmax();
0216     std::string xtit = me1->getAxisTitle(1);
0217     std::string tit_new;
0218     tit_new = ";" + xtit + ";Efficiency";
0219
0220     ibooker.setCurrentFolder(folder);
0221     me_eff = ibooker.book1D("efficiency_" + name, tit_new, nbinx, xmin, xmax);
0222
0223     double efficiency;
0224     me_eff->Reset();
0225     me_eff->setEfficiencyFlag();
0226     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0227       float val1 = me1->getBinContent(ix);
0228       float val2 = me2->getBinContent(ix);
0229       if (val2 > 0.0)
0230         efficiency = val1 / val2;
0231       else
0232         efficiency = 0;
0233       me_eff->setBinContent(ix, efficiency);
0234     }
0235   }
0236 }
0237
0238 #include "FWCore/Framework/interface/MakerMacros.h"
0239 DEFINE_FWK_MODULE(PFClient_JetRes);