Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:31

0001 #include "DQMOffline/PFTau/plugins/PFClient.h"
0002 
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 
0010 //
0011 // -- Constructor
0012 //
0013 PFClient::PFClient(const edm::ParameterSet &parameterSet) {
0014   folderNames_ = parameterSet.getParameter<std::vector<std::string>>("FolderNames");
0015   histogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNames");
0016   efficiencyFlag_ = parameterSet.getParameter<bool>("CreateEfficiencyPlots");
0017   effHistogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNamesForEfficiencyPlots");
0018   projectionHistogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNamesForProjectionPlots");
0019   profileFlag_ = parameterSet.getParameter<bool>("CreateProfilePlots");
0020   profileHistogramNames_ = parameterSet.getParameter<std::vector<std::string>>("HistogramNamesForProfilePlots");
0021 }
0022 
0023 //
0024 // -- EndJobBegin Run
0025 //
0026 void PFClient::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0027   doSummaries(ibooker, igetter);
0028   doProjection(ibooker, igetter);
0029   if (efficiencyFlag_)
0030     doEfficiency(ibooker, igetter);
0031   if (profileFlag_)
0032     doProfiles(ibooker, igetter);
0033 }
0034 
0035 //
0036 // -- Create Summaries
0037 //
0038 void PFClient::doSummaries(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0039   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0040        ifolder++) {
0041     std::string path = "ParticleFlow/" + (*ifolder);
0042 
0043     for (std::vector<std::string>::const_iterator ihist = histogramNames_.begin(); ihist != histogramNames_.end();
0044          ihist++) {
0045       std::string hname = (*ihist);
0046       createResolutionPlots(ibooker, igetter, path, hname);
0047     }
0048   }
0049 }
0050 
0051 //
0052 // -- Create Projection
0053 //
0054 void PFClient::doProjection(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0055   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0056        ifolder++) {
0057     std::string path = "ParticleFlow/" + (*ifolder);
0058 
0059     for (std::vector<std::string>::const_iterator ihist = projectionHistogramNames_.begin();
0060          ihist != projectionHistogramNames_.end();
0061          ihist++) {
0062       std::string hname = (*ihist);
0063       createProjectionPlots(ibooker, igetter, path, hname);
0064     }
0065   }
0066 }
0067 
0068 //
0069 // -- Create Profile
0070 //
0071 void PFClient::doProfiles(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0072   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0073        ifolder++) {
0074     std::string path = "ParticleFlow/" + (*ifolder);
0075 
0076     for (std::vector<std::string>::const_iterator ihist = profileHistogramNames_.begin();
0077          ihist != profileHistogramNames_.end();
0078          ihist++) {
0079       std::string hname = (*ihist);
0080       createProfilePlots(ibooker, igetter, path, hname);
0081     }
0082   }
0083 }
0084 
0085 //
0086 // -- Create Efficiency
0087 //
0088 void PFClient::doEfficiency(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0089   for (std::vector<std::string>::const_iterator ifolder = folderNames_.begin(); ifolder != folderNames_.end();
0090        ifolder++) {
0091     std::string path = "ParticleFlow/" + (*ifolder);
0092 
0093     for (std::vector<std::string>::const_iterator ihist = effHistogramNames_.begin(); ihist != effHistogramNames_.end();
0094          ihist++) {
0095       std::string hname = (*ihist);
0096       createEfficiencyPlots(ibooker, igetter, path, hname);
0097     }
0098   }
0099 }
0100 
0101 //
0102 // -- Create Resolution Plots
0103 //
0104 void PFClient::createResolutionPlots(DQMStore::IBooker &ibooker,
0105                                      DQMStore::IGetter &igetter,
0106                                      std::string &folder,
0107                                      std::string &name) {
0108   MonitorElement *me = igetter.get(folder + "/" + name);
0109 
0110   if (!me)
0111     return;
0112 
0113   MonitorElement *me_average;
0114   MonitorElement *me_rms;
0115   MonitorElement *me_mean;
0116   MonitorElement *me_sigma;
0117 
0118   if ((me->kind() == MonitorElement::Kind::TH2F) || (me->kind() == MonitorElement::Kind::TH2S) ||
0119       (me->kind() == MonitorElement::Kind::TH2D)) {
0120     TH2 *th = me->getTH2F();
0121     size_t nbinx = me->getNbinsX();
0122     size_t nbiny = me->getNbinsY();
0123 
0124     float ymin = th->GetYaxis()->GetXmin();
0125     float ymax = th->GetYaxis()->GetXmax();
0126     std::string xtit = th->GetXaxis()->GetTitle();
0127     std::string ytit = th->GetYaxis()->GetTitle();
0128     float *xbins = new float[nbinx + 1];
0129     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0130       xbins[ix - 1] = th->GetXaxis()->GetBinLowEdge(ix);
0131       if (ix == nbinx)
0132         xbins[ix] = th->GetXaxis()->GetBinUpEdge(ix);
0133     }
0134     std::string tit_new = ";" + xtit + ";" + ytit;
0135     ibooker.setCurrentFolder(folder);
0136     MonitorElement *me_slice = ibooker.book1D("PFlowSlice", "PFlowSlice", nbiny, ymin, ymax);
0137 
0138     tit_new = ";" + xtit + ";Average_" + ytit;
0139     me_average = ibooker.book1D("average_" + name, tit_new, nbinx, xbins);
0140     me_average->setEfficiencyFlag();
0141     tit_new = ";" + xtit + ";RMS_" + ytit;
0142     me_rms = ibooker.book1D("rms_" + name, tit_new, nbinx, xbins);
0143     me_rms->setEfficiencyFlag();
0144     tit_new = ";" + xtit + ";Mean_" + ytit;
0145     me_mean = ibooker.book1D("mean_" + name, tit_new, nbinx, xbins);
0146     me_mean->setEfficiencyFlag();
0147     tit_new = ";" + xtit + ";Sigma_" + ytit;
0148     me_sigma = ibooker.book1D("sigma_" + name, tit_new, nbinx, xbins);
0149     me_sigma->setEfficiencyFlag();
0150 
0151     double average, rms, mean, sigma;
0152 
0153     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0154       me_slice->Reset();
0155       for (size_t iy = 1; iy < nbiny + 1; ++iy) {
0156         me_slice->setBinContent(iy, th->GetBinContent(ix, iy));
0157       }
0158       getHistogramParameters(me_slice, average, rms, mean, sigma);
0159       me_average->setBinContent(ix, average);
0160       me_rms->setBinContent(ix, rms);
0161       me_mean->setBinContent(ix, mean);
0162       me_sigma->setBinContent(ix, sigma);
0163     }
0164     delete[] xbins;
0165   }
0166 }
0167 
0168 //
0169 // -- Create Projection Plots
0170 //
0171 void PFClient::createProjectionPlots(DQMStore::IBooker &ibooker,
0172                                      DQMStore::IGetter &igetter,
0173                                      std::string &folder,
0174                                      std::string &name) {
0175   MonitorElement *me = igetter.get(folder + "/" + name);
0176   if (!me)
0177     return;
0178 
0179   MonitorElement *projection = nullptr;
0180 
0181   if ((me->kind() == MonitorElement::Kind::TH2F) || (me->kind() == MonitorElement::Kind::TH2S) ||
0182       (me->kind() == MonitorElement::Kind::TH2D)) {
0183     TH2 *th = me->getTH2F();
0184     size_t nbinx = me->getNbinsX();
0185     size_t nbiny = me->getNbinsY();
0186 
0187     float ymin = th->GetYaxis()->GetXmin();
0188     float ymax = th->GetYaxis()->GetXmax();
0189     std::string xtit = th->GetXaxis()->GetTitle();
0190     std::string ytit = th->GetYaxis()->GetTitle();
0191     float *xbins = new float[nbinx + 1];
0192     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0193       xbins[ix - 1] = th->GetXaxis()->GetBinLowEdge(ix);
0194       if (ix == nbinx)
0195         xbins[ix] = th->GetXaxis()->GetBinUpEdge(ix);
0196     }
0197 
0198     std::string tit_new;
0199     ibooker.setCurrentFolder(folder);
0200 
0201     if (folder == "ParticleFlow/PFElectronValidation/CompWithGenElectron") {
0202       if (name == "delta_et_Over_et_VS_et_")
0203         projection = ibooker.book1D("delta_et_Over_et", "E_{T} resolution;#DeltaE_{T}/E_{T}", nbiny, ymin, ymax);
0204       if (name == "delta_et_VS_et_")
0205         projection = ibooker.book1D("delta_et_", "#DeltaE_{T};#DeltaE_{T}", nbiny, ymin, ymax);
0206       if (name == "delta_eta_VS_et_")
0207         projection = ibooker.book1D("delta_eta_", "#Delta#eta;#Delta#eta", nbiny, ymin, ymax);
0208       if (name == "delta_phi_VS_et_")
0209         projection = ibooker.book1D("delta_phi_", "#Delta#phi;#Delta#phi", nbiny, ymin, ymax);
0210     }
0211 
0212     if (projection) {
0213       for (size_t iy = 1; iy < nbiny + 1; ++iy) {
0214         projection->setBinContent(iy, th->ProjectionY("e")->GetBinContent(iy));
0215       }
0216       projection->setEntries(me->getEntries());
0217     }
0218 
0219     delete[] xbins;
0220   }
0221 }
0222 
0223 //
0224 // -- Create Profile Plots
0225 //
0226 void PFClient::createProfilePlots(DQMStore::IBooker &ibooker,
0227                                   DQMStore::IGetter &igetter,
0228                                   std::string &folder,
0229                                   std::string &name) {
0230   MonitorElement *me = igetter.get(folder + "/" + name);
0231   if (!me)
0232     return;
0233 
0234   if ((me->kind() == MonitorElement::Kind::TH2F) || (me->kind() == MonitorElement::Kind::TH2S) ||
0235       (me->kind() == MonitorElement::Kind::TH2D)) {
0236     TH2 *th = me->getTH2F();
0237     size_t nbinx = me->getNbinsX();
0238 
0239     float ymin = th->GetYaxis()->GetXmin();
0240     float ymax = th->GetYaxis()->GetXmax();
0241     std::string xtit = th->GetXaxis()->GetTitle();
0242     std::string ytit = th->GetYaxis()->GetTitle();
0243     double *xbins = new double[nbinx + 1];
0244     for (size_t ix = 1; ix < nbinx + 1; ++ix) {
0245       xbins[ix - 1] = th->GetXaxis()->GetBinLowEdge(ix);
0246       if (ix == nbinx)
0247         xbins[ix] = th->GetXaxis()->GetBinUpEdge(ix);
0248     }
0249 
0250     std::string tit_new;
0251     ibooker.setCurrentFolder(folder);
0252     // TProfiles
0253     MonitorElement *me_profile[2];
0254     me_profile[0] = ibooker.bookProfile("profile_" + name, tit_new, nbinx, xbins, ymin, ymax, "");
0255     me_profile[1] = ibooker.bookProfile("profileRMS_" + name, tit_new, nbinx, xbins, ymin, ymax, "s");
0256     TProfile *profileX = th->ProfileX();
0257     // size_t nbiny = me->getNbinsY();
0258     // TProfile* profileX = th->ProfileX("",0,nbiny+1); add underflow and
0259     // overflow
0260     static const Int_t NUM_STAT = 7;
0261     Double_t stats[NUM_STAT] = {0};
0262     th->GetStats(stats);
0263 
0264     for (Int_t i = 0; i < 2; i++) {
0265       if (me_profile[i]) {
0266         for (size_t ix = 0; ix <= nbinx + 1; ++ix) {
0267           me_profile[i]->setBinContent(ix, profileX->GetBinContent(ix) * profileX->GetBinEntries(ix));
0268           // me_profile[i]->Fill( profileX->GetBinCenter(ix),
0269           // profileX->GetBinContent(ix)*profileX->GetBinEntries(ix) ) ;
0270           me_profile[i]->setBinEntries(ix, profileX->GetBinEntries(ix));
0271           me_profile[i]->getTProfile()->GetSumw2()->fArray[ix] = profileX->GetSumw2()->fArray[ix];
0272           // me_profile[i]->getTProfile()->GetBinSumw2()->fArray[ix] =
0273           // profileX->GetBinSumw2()->fArray[ix]; // segmentation violation
0274         }
0275       }
0276       me_profile[i]->getTProfile()->PutStats(stats);
0277       me_profile[i]->setEntries(profileX->GetEntries());
0278     }
0279 
0280     delete[] xbins;
0281   }
0282 }
0283 
0284 //
0285 // -- Get Histogram Parameters
0286 //
0287 void PFClient::getHistogramParameters(
0288     MonitorElement *me_slice, double &average, double &rms, double &mean, double &sigma) {
0289   average = 0.0;
0290   rms = 0.0;
0291   mean = 0.0;
0292   sigma = 0.0;
0293 
0294   if (!me_slice)
0295     return;
0296   if (me_slice->kind() == MonitorElement::Kind::TH1F) {
0297     average = me_slice->getMean();
0298     rms = me_slice->getRMS();
0299     TH1F *th_slice = me_slice->getTH1F();
0300     if (th_slice && th_slice->GetEntries() > 0) {
0301       // need our own copy for thread safety
0302       TF1 gaus("mygaus", "gaus");
0303       th_slice->Fit(&gaus, "Q0 SERIAL");
0304       sigma = gaus.GetParameter(2);
0305       mean = gaus.GetParameter(1);
0306     }
0307   }
0308 }
0309 
0310 //
0311 // -- Create Efficiency Plots
0312 //
0313 void PFClient::createEfficiencyPlots(DQMStore::IBooker &ibooker,
0314                                      DQMStore::IGetter &igetter,
0315                                      std::string &folder,
0316                                      std::string &name) {
0317   MonitorElement *me1 = igetter.get(folder + "/" + name + "ref_");
0318   MonitorElement *me2 = igetter.get(folder + "/" + name + "gen_");
0319   if (!me1 || !me2)
0320     return;
0321 
0322   TH1F *me1_forEff = (TH1F *)me1->getTH1F()->Rebin(2, "me1_forEff");
0323   TH1F *me2_forEff = (TH1F *)me2->getTH1F()->Rebin(2, "me2_forEff");
0324 
0325   MonitorElement *me_eff;
0326   if ((me1->kind() == MonitorElement::Kind::TH1F) && (me1->kind() == MonitorElement::Kind::TH1F)) {
0327     // TH1* th1 = me1->getTH1F();
0328     // size_t nbinx = me1->getNbinsX();
0329     size_t nbinx = me1_forEff->GetNbinsX();
0330 
0331     // float xmin = th1->GetXaxis()->GetXmin();
0332     // float xmax = th1->GetXaxis()->GetXmax();
0333     float xmin = me1_forEff->GetXaxis()->GetXmin();
0334     float xmax = me1_forEff->GetXaxis()->GetXmax();
0335     std::string xtit = me1->getAxisTitle(1);
0336     std::string tit_new;
0337     tit_new = ";" + xtit + ";" + xtit + " efficiency";
0338 
0339     ibooker.setCurrentFolder(folder);
0340     me_eff = ibooker.book1D("efficiency_" + name, tit_new, nbinx, xmin, xmax);
0341 
0342     me_eff->Reset();
0343     me_eff->setEfficiencyFlag();
0344     /*
0345     double  efficiency;
0346     for (size_t ix = 1; ix < nbinx+1; ++ix) {
0347       float val1 = me1->getBinContent(ix);
0348       float val2 = me2->getBinContent(ix);
0349       if (val2 > 0.0) efficiency = val1/val2;
0350       else efficiency = 0;
0351       me_eff->setBinContent(ix,efficiency);
0352     }
0353     */
0354     // Binomial errors "B" asked by Florian
0355     /*me1_forEff->Sumw2(); me2_forEff->Sumw2();*/ me_eff->enableSumw2();
0356     me_eff->getTH1F()->Divide(me1_forEff, me2_forEff, 1, 1, "B");
0357   }
0358 }
0359 
0360 #include "FWCore/Framework/interface/MakerMacros.h"
0361 DEFINE_FWK_MODULE(PFClient);