Egamma/src/EmDQMPostProcessor.cc

0001 #include "HLTriggerOffline/Egamma/interface/EmDQMPostProcessor.h"
0002
0003 #include "DQMServices/Core/interface/DQMStore.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006
0007 #include <cmath>
0008 #include <cstring>
0009 #include <fstream>
0010 #include <iomanip>
0011 #include <iostream>
0012
0013 #include <TEfficiency.h>
0014
0015 //----------------------------------------------------------------------
0016
0017 EmDQMPostProcessor::EmDQMPostProcessor(const edm::ParameterSet &pset) {
0018   subDir_ = pset.getUntrackedParameter<std::string>("subDir");
0019
0020   dataSet_ = pset.getUntrackedParameter<std::string>("dataSet", "unknown");
0021
0022   noPhiPlots = pset.getUntrackedParameter<bool>("noPhiPlots", true);
0023
0024   normalizeToReco = pset.getUntrackedParameter<bool>("normalizeToReco", false);
0025
0026   ignoreEmpty = pset.getUntrackedParameter<bool>("ignoreEmpty", true);
0027 }
0028
0029 //----------------------------------------------------------------------
0030
0031 void EmDQMPostProcessor::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0032   // go to the directory to be processed
0033   if (igetter.dirExists(subDir_))
0034     ibooker.cd(subDir_);
0035   else {
0036     edm::LogWarning("EmDQMPostProcessor") << "cannot find directory: " << subDir_ << " , skipping";
0037     return;
0038   }
0039
0040   //--------------------
0041   // with respect to what are (some) efficiencies calculated ?
0042   std::string shortReferenceName;
0043   if (normalizeToReco)
0044     shortReferenceName = "RECO";
0045   else
0046     shortReferenceName = "gen";
0047   //--------------------
0048
0049   //////////////////////////////////
0050   // loop over all triggers/samples//
0051   //////////////////////////////////
0052
0053   // store dataset name (if defined) in output file
0054   // DQMStore:bookString seems to put a key in the file which is
0055   // of the form <dataSet>s=....</dataSet> which is not very convenient
0056   // (it points to a null pointer, one would have to loop over
0057   // all keys of the corresponding directory in the ROOT file
0058   // and check whether it is of the desired form and then parse
0059   // it from this string...).
0060   //
0061   // So we store the name of the dataset as the title of a histogram,
0062   // which is much easier to access...
0063   // TH1D *dataSetNameHisto =
0064   ibooker.book1D("DataSetNameHistogram", dataSet_, 1, 0, 1);
0065
0066   std::vector<std::string> subdirectories = igetter.getSubdirs();
0067   ////////////////////////////////////////////////////////
0068   // Do everything twice: once for mc-matched histos,   //
0069   // once for unmatched histos                          //
0070   ////////////////////////////////////////////////////////
0071
0072   std::vector<std::string> postfixes;
0073   postfixes.push_back("");               // unmatched histograms
0074   postfixes.push_back("_RECO_matched");  // for data
0075   // we put this on the list even when we're running on
0076   // data (where there is no generator information).
0077   // The first test in the loop will then fail and
0078   // the iteration is skipped.
0079   postfixes.push_back("_MC_matched");
0080
0081   std::vector<TProfile *> allElePaths;
0082   int nEle = 0;
0083   int nPhoton = 0;
0084
0085   // find the number of electron and photon paths
0086   for (std::vector<std::string>::iterator dir = subdirectories.begin(); dir != subdirectories.end(); ++dir) {
0087     if (dir->find("Ele") != std::string::npos || dir->find("_SC") != std::string::npos)
0088       ++nEle;
0089     else if (dir->find("Photon") != std::string::npos)
0090       ++nPhoton;
0091   }
0092
0093   std::vector<TProfile *> allPhotonPaths;
0094   for (std::vector<std::string>::iterator postfix = postfixes.begin(); postfix != postfixes.end(); postfix++) {
0095     bool pop = false;
0096     int elePos = 1;
0097     int photonPos = 1;
0098
0099     /////////////////////////////////////
0100     // computer per-event efficiencies //
0101     /////////////////////////////////////
0102
0103     std::string histoName = "efficiency_by_step" + *postfix;
0104     std::string baseName = "total_eff" + *postfix;
0105
0106     std::string allEleHistoName = "EfficiencyByPath_Ele" + *postfix;
0107     std::string allEleHistoLabel = "Efficiency_for_each_validated_electron_path" + *postfix;
0108     allElePaths.push_back(
0109         new TProfile(allEleHistoName.c_str(), allEleHistoLabel.c_str(), nEle, 0., (double)nEle, 0., 1.2));
0110     std::string allPhotonHistoName = "EfficiencyByPath_Photon" + *postfix;
0111     std::string allPhotonHistoLabel = "Efficiency_for_each_validated_photon_path" + *postfix;
0112     allPhotonPaths.push_back(
0113         new TProfile(allPhotonHistoName.c_str(), allPhotonHistoLabel.c_str(), nPhoton, 0., (double)nPhoton, 0., 1.2));
0114
0115     for (std::vector<std::string>::iterator dir = subdirectories.begin(); dir != subdirectories.end(); dir++) {
0116       ibooker.cd(*dir);
0117
0118       // get the current trigger name
0119       std::string trigName = dir->substr(dir->rfind("/") + 1);
0120       trigName = trigName.replace(trigName.rfind("_DQM"), 4, "");
0121
0122       // Get the gen-level (or reco, for data) plots
0123       std::string genName;
0124
0125       // only generate efficiency plots if there are generated/recoed events
0126       // selected but label the bin in the overview plot, even if the bin is
0127       // empty
0128       if (ignoreEmpty) {
0129         if (normalizeToReco)
0130           genName = ibooker.pwd() + "/reco_et";
0131         else
0132           genName = ibooker.pwd() + "/gen_et";
0133         TH1F *genHist = getHistogram(ibooker, igetter, genName);
0134         if (genHist->GetEntries() == 0) {
0135           edm::LogInfo("EmDQMPostProcessor")
0136               << "Zero events were selected for path '" << trigName << "'. No efficiency plots will be generated.";
0137           if (trigName.find("Ele") != std::string::npos || trigName.find("_SC") != std::string::npos) {
0138             allElePaths.back()->GetXaxis()->SetBinLabel(elePos, trigName.c_str());
0139             ++elePos;
0140           } else if (trigName.find("Photon") != std::string::npos) {
0141             allPhotonPaths.back()->GetXaxis()->SetBinLabel(photonPos, trigName.c_str());
0142             ++photonPos;
0143           }
0144           ibooker.goUp();
0145           continue;
0146         }
0147       }
0148
0149       TH1F *basehist = getHistogram(ibooker, igetter, ibooker.pwd() + "/" + baseName);
0150       if (basehist == nullptr) {
0151         // edm::LogWarning("EmDQMPostProcessor") << "histogram " <<
0152         // (ibooker.pwd() + "/" + baseName) << " does not exist, skipping
0153         // postfix
0154         // '" << *postfix << "'";
0155         pop = true;
0156         ibooker.goUp();
0157         continue;
0158       }
0159       // at least one histogram with postfix was found
0160       pop = false;
0161
0162       ibooker.goUp();
0163       MonitorElement *meTotal = ibooker.bookProfile(trigName + "__" + histoName,
0164                                                     trigName + "__" + histoName,
0165                                                     basehist->GetXaxis()->GetNbins(),
0166                                                     basehist->GetXaxis()->GetXmin(),
0167                                                     basehist->GetXaxis()->GetXmax(),
0168                                                     0.,
0169                                                     1.2);
0170       meTotal->setEfficiencyFlag();
0171       TProfile *total = meTotal->getTProfile();
0172       ibooker.cd(*dir);
0173       total->GetXaxis()->SetBinLabel(1, basehist->GetXaxis()->GetBinLabel(1));
0174
0175       //       std::vector<std::string> mes = igetter.getMEs();
0176       //       for(std::vector<std::string>::iterator me = mes.begin() ;me!=
0177       //       mes.end(); me++ )
0178       //    std::cout <<*me <<std::endl;
0179       //       std::cout <<std::endl;
0180
0181       double value = 0;
0182       double errorh = 0, errorl = 0, error = 0;
0183       // compute stepwise total efficiencies
0184       for (int bin = total->GetNbinsX() - 2; bin > 1; bin--) {
0185         value = 0;
0186         errorl = 0;
0187         errorh = 0;
0188         error = 0;
0189         if (basehist->GetBinContent(bin - 1) != 0) {
0190           Efficiency(
0191               (int)basehist->GetBinContent(bin), (int)basehist->GetBinContent(bin - 1), 0.683, value, errorl, errorh);
0192           error = value - errorl > errorh - value ? value - errorl : errorh - value;
0193         }
0194         total->SetBinContent(bin, value);
0195         total->SetBinEntries(bin, 1);
0196         total->SetBinError(bin, sqrt(value * value + error * error));
0197         total->GetXaxis()->SetBinLabel(bin, basehist->GetXaxis()->GetBinLabel(bin));
0198       }
0199
0200       // set first bin to L1 efficiency
0201       if (basehist->GetBinContent(basehist->GetNbinsX()) != 0) {
0202         Efficiency((int)basehist->GetBinContent(1),
0203                    (int)basehist->GetBinContent(basehist->GetNbinsX()),
0204                    0.683,
0205                    value,
0206                    errorl,
0207                    errorh);
0208         error = value - errorl > errorh - value ? value - errorl : errorh - value;
0209       } else {
0210         value = 0;
0211         error = 0;
0212       }
0213       total->SetBinContent(1, value);
0214       total->SetBinEntries(1, 1);
0215       total->SetBinError(1, sqrt(value * value + error * error));
0216
0217       // total efficiency relative to gen or reco
0218       if (basehist->GetBinContent(basehist->GetNbinsX()) != 0) {
0219         Efficiency((int)basehist->GetBinContent(basehist->GetNbinsX() - 2),
0220                    (int)basehist->GetBinContent(basehist->GetNbinsX()),
0221                    0.683,
0222                    value,
0223                    errorl,
0224                    errorh);
0225         error = value - errorl > errorh - value ? value - errorl : errorh - value;
0226       } else {
0227         value = 0;
0228         error = 0;
0229       }
0230       total->SetBinContent(total->GetNbinsX(), value);
0231       total->SetBinEntries(total->GetNbinsX(), 1);
0232       total->SetBinError(total->GetNbinsX(), sqrt(value * value + error * error));
0233       total->GetXaxis()->SetBinLabel(total->GetNbinsX(), ("total efficiency rel. " + shortReferenceName).c_str());
0234
0235       // total efficiency relative to L1
0236       if (basehist->GetBinContent(1) != 0) {
0237         Efficiency((int)basehist->GetBinContent(basehist->GetNbinsX() - 2),
0238                    (int)basehist->GetBinContent(1),
0239                    0.683,
0240                    value,
0241                    errorl,
0242                    errorh);
0243         error = value - errorl > errorh - value ? value - errorl : errorh - value;
0244       } else {
0245         value = 0;
0246         error = 0;
0247       }
0248       total->SetBinContent(total->GetNbinsX() - 1, value);
0249       total->SetBinError(total->GetNbinsX() - 1, sqrt(value * value + error * error));
0250       total->SetBinEntries(total->GetNbinsX() - 1, 1);
0251       total->GetXaxis()->SetBinLabel(total->GetNbinsX() - 1, "total efficiency rel. L1");
0252
0253       //----------------------------------------
0254
0255       ///////////////////////////////////////////
0256       // compute per-object efficiencies       //
0257       ///////////////////////////////////////////
0258       // MonitorElement *eff, *num, *denom, *genPlot, *effVsGen, *effL1VsGen;
0259       std::vector<std::string> varNames;
0260       varNames.push_back("et");
0261       varNames.push_back("eta");
0262
0263       if (!noPhiPlots) {
0264         varNames.push_back("phi");
0265       }
0266       varNames.push_back("etaphi");
0267
0268       std::string filterName;
0269       std::string filterName2;
0270       std::string denomName;
0271       std::string numName;
0272
0273       // Get the L1 over gen filter first
0274       filterName2 = total->GetXaxis()->GetBinLabel(1);
0275       // loop over variables (eta/phi/et)
0276       for (std::vector<std::string>::iterator var = varNames.begin(); var != varNames.end(); var++) {
0277         numName = ibooker.pwd() + "/" + filterName2 + *var + *postfix;
0278
0279         if (normalizeToReco)
0280           genName = ibooker.pwd() + "/reco_" + *var;
0281         else
0282           genName = ibooker.pwd() + "/gen_" + *var;
0283
0284         if ((*var).find("etaphi") != std::string::npos) {
0285           if (!dividehistos2D(ibooker,
0286                               igetter,
0287                               numName,
0288                               genName,
0289                               "efficiency_" + filterName2 + "_vs_" + *var + *postfix,
0290                               *var,
0291                               "eff. of" + filterName2 + " vs " + *var + *postfix))
0292             break;
0293         } else if (!dividehistos(ibooker,
0294                                  igetter,
0295                                  numName,
0296                                  genName,
0297                                  "efficiency_" + filterName2 + "_vs_" + *var + *postfix,
0298                                  *var,
0299                                  "eff. of" + filterName2 + " vs " + *var + *postfix))
0300           break;
0301       }  // loop over variables
0302
0303       // get the filter names from the bin-labels of the master-histogram
0304       for (int filter = 1; filter < total->GetNbinsX() - 2; filter++) {
0305         filterName = total->GetXaxis()->GetBinLabel(filter);
0306         filterName2 = total->GetXaxis()->GetBinLabel(filter + 1);
0307
0308         // loop over variables (eta/et/phi)
0309         for (std::vector<std::string>::iterator var = varNames.begin(); var != varNames.end(); var++) {
0310           numName = ibooker.pwd() + "/" + filterName2 + *var + *postfix;
0311           denomName = ibooker.pwd() + "/" + filterName + *var + *postfix;
0312
0313           // Is this the last filter? Book efficiency vs gen (or reco, for data)
0314           // level
0315           std::string temp = *postfix;
0316           if (filter == total->GetNbinsX() - 3 && temp.find("matched") != std::string::npos) {
0317             if (normalizeToReco)
0318               genName = ibooker.pwd() + "/reco_" + *var;
0319             else
0320               genName = ibooker.pwd() + "/gen_" + *var;
0321
0322             if ((*var).find("etaphi") != std::string::npos) {
0323               if (!dividehistos2D(ibooker,
0324                                   igetter,
0325                                   numName,
0326                                   genName,
0327                                   "final_eff_vs_" + *var,
0328                                   *var,
0329                                   "Efficiency Compared to " + shortReferenceName + " vs " + *var))
0330                 break;
0331             } else if (!dividehistos(ibooker,
0332                                      igetter,
0333                                      numName,
0334                                      genName,
0335                                      "final_eff_vs_" + *var,
0336                                      *var,
0337                                      "Efficiency Compared to " + shortReferenceName + " vs " + *var))
0338               break;
0339           }
0340
0341           if ((*var).find("etaphi") != std::string::npos) {
0342             if (!dividehistos2D(ibooker,
0343                                 igetter,
0344                                 numName,
0345                                 denomName,
0346                                 "efficiency_" + filterName2 + "_vs_" + *var + *postfix,
0347                                 *var,
0348                                 "efficiency_" + filterName2 + "_vs_" + *var + *postfix))
0349               break;
0350           } else if (!dividehistos(ibooker,
0351                                    igetter,
0352                                    numName,
0353                                    denomName,
0354                                    "efficiency_" + filterName2 + "_vs_" + *var + *postfix,
0355                                    *var,
0356                                    "efficiency_" + filterName2 + "_vs_" + *var + *postfix))
0357             break;
0358
0359         }  // loop over variables
0360       }    // loop over monitoring modules within path
0361
0362       ibooker.goUp();
0363
0364       // fill overall efficiency histograms
0365       double totCont = total->GetBinContent(total->GetNbinsX());
0366       double totErr = total->GetBinError(total->GetNbinsX());
0367       if (trigName.find("Ele") != std::string::npos || trigName.find("_SC") != std::string::npos) {
0368         allElePaths.back()->SetBinContent(elePos, totCont);
0369         allElePaths.back()->SetBinEntries(elePos, 1);
0370         allElePaths.back()->SetBinError(elePos, sqrt(totCont * totCont + totErr * totErr));
0371         allElePaths.back()->GetXaxis()->SetBinLabel(elePos, trigName.c_str());
0372         ++elePos;
0373       } else if (trigName.find("Photon") != std::string::npos) {
0374         allPhotonPaths.back()->SetBinContent(photonPos, totCont);
0375         allPhotonPaths.back()->SetBinEntries(photonPos, 1);
0376         allPhotonPaths.back()->SetBinError(photonPos, sqrt(totCont * totCont + totErr * totErr));
0377         allPhotonPaths.back()->GetXaxis()->SetBinLabel(photonPos, trigName.c_str());
0378         ++photonPos;
0379       }
0380
0381     }  // loop over dirs
0382     if (pop) {
0383       allElePaths.pop_back();
0384       allPhotonPaths.pop_back();
0385     } else {
0386       allElePaths.back()->GetXaxis()->SetLabelSize(0.03);
0387       allPhotonPaths.back()->GetXaxis()->SetLabelSize(0.03);
0388       ibooker.bookProfile(allEleHistoName, allElePaths.back())->setEfficiencyFlag();
0389       ibooker.bookProfile(allPhotonHistoName, allPhotonPaths.back())->setEfficiencyFlag();
0390     }
0391   }  // loop over postfixes
0392 }
0393
0394 //----------------------------------------------------------------------
0395
0396 TProfile *EmDQMPostProcessor::dividehistos(DQMStore::IBooker &ibooker,
0397                                            DQMStore::IGetter &igetter,
0398                                            const std::string &numName,
0399                                            const std::string &denomName,
0400                                            const std::string &outName,
0401                                            const std::string &label,
0402                                            const std::string &titel) {
0403   // std::cout << numName <<std::endl;
0404
0405   TH1F *num = getHistogram(ibooker, igetter, numName);
0406
0407   // std::cout << denomName << std::endl;
0408   TH1F *denom = getHistogram(ibooker, igetter, denomName);
0409
0410   if (num == nullptr)
0411     edm::LogWarning("EmDQMPostProcessor") << "numerator histogram " << numName << " does not exist";
0412
0413   if (denom == nullptr)
0414     edm::LogWarning("EmDQMPostProcessor") << "denominator histogram " << denomName << " does not exist";
0415
0416   // Check if histograms actually exist
0417
0418   if (num == nullptr || denom == nullptr)
0419     return nullptr;
0420
0421   MonitorElement *meOut = ibooker.bookProfile(
0422       outName, titel, num->GetXaxis()->GetNbins(), num->GetXaxis()->GetXmin(), num->GetXaxis()->GetXmax(), 0., 1.2);
0423   meOut->setEfficiencyFlag();
0424   TProfile *out = meOut->getTProfile();
0425   out->GetXaxis()->SetTitle(label.c_str());
0426   out->SetYTitle("Efficiency");
0427   out->SetOption("PE");
0428   out->SetLineColor(2);
0429   out->SetLineWidth(2);
0430   out->SetMarkerStyle(20);
0431   out->SetMarkerSize(0.8);
0432   out->SetStats(kFALSE);
0433   for (int i = 1; i <= num->GetNbinsX(); i++) {
0434     double e, low, high;
0435     Efficiency((int)num->GetBinContent(i), (int)denom->GetBinContent(i), 0.683, e, low, high);
0436     double err = e - low > high - e ? e - low : high - e;
0437     // here is the trick to store info in TProfile:
0438     out->SetBinContent(i, e);
0439     out->SetBinEntries(i, 1);
0440     out->SetBinError(i, sqrt(e * e + err * err));
0441   }
0442
0443   return out;
0444 }
0445 //----------------------------------------------------------------------
0446 TH2F *EmDQMPostProcessor::dividehistos2D(DQMStore::IBooker &ibooker,
0447                                          DQMStore::IGetter &igetter,
0448                                          const std::string &numName,
0449                                          const std::string &denomName,
0450                                          const std::string &outName,
0451                                          const std::string &label,
0452                                          const std::string &titel) {
0453   // std::cout << numName <<std::endl;
0454   TH2F *num = get2DHistogram(ibooker, igetter, numName);
0455   // std::cout << denomName << std::endl;
0456   TH2F *denom = get2DHistogram(ibooker, igetter, denomName);
0457   if (num == nullptr)
0458     edm::LogWarning("EmDQMPostProcessor") << "2D numerator histogram " << numName << " does not exist";
0459
0460   if (denom == nullptr)
0461     edm::LogWarning("EmDQMPostProcessor") << "2D denominator histogram " << denomName << " does not exist";
0462
0463   // Check if histograms actually exist
0464   if (num == nullptr || denom == nullptr)
0465     return nullptr;
0466
0467   MonitorElement *meOut = ibooker.book2D(outName,
0468                                          titel,
0469                                          num->GetXaxis()->GetNbins(),
0470                                          num->GetXaxis()->GetXmin(),
0471                                          num->GetXaxis()->GetXmax(),
0472                                          num->GetYaxis()->GetNbins(),
0473                                          num->GetYaxis()->GetXmin(),
0474                                          num->GetYaxis()->GetXmax());
0475   TH2F *out = meOut->getTH2F();
0476   out->Add(num);
0477   out->Divide(denom);
0478   out->GetXaxis()->SetTitle(label.c_str());
0479   out->SetYTitle("#phi");
0480   out->SetXTitle("#eta");
0481   out->SetOption("COLZ");
0482   out->SetStats(kFALSE);
0483
0484   return out;
0485 }
0486 //--------------------
0487 TH1F *EmDQMPostProcessor::getHistogram(DQMStore::IBooker &ibooker,
0488                                        DQMStore::IGetter &igetter,
0489                                        const std::string &histoPath) {
0490   MonitorElement *monElement = igetter.get(histoPath);
0491   if (monElement != nullptr)
0492     return monElement->getTH1F();
0493   else
0494     return nullptr;
0495 }
0496 //----------------------------------------------------------------------
0497 TH2F *EmDQMPostProcessor::get2DHistogram(DQMStore::IBooker &ibooker,
0498                                          DQMStore::IGetter &igetter,
0499                                          const std::string &histoPath) {
0500   MonitorElement *monElement = igetter.get(histoPath);
0501   if (monElement != nullptr)
0502     return monElement->getTH2F();
0503   else
0504     return nullptr;
0505 }
0506
0507 //----------------------------------------------------------------------
0508
0509 void EmDQMPostProcessor::Efficiency(
0510     int passing, int total, double level, double &mode, double &lowerBound, double &upperBound) {
0511   // protection (see also TGraphAsymmErrors::Efficiency(..), mimick the old
0512   // behaviour )
0513   if (total == 0) {
0514     mode = 0.5;
0515     lowerBound = 0;
0516     upperBound = 1;
0517     return;
0518   }
0519
0520   mode = passing / ((double)total);
0521
0522   // note that the order of the two first arguments ('total' and 'passed') is
0523   // the opposited with respect to the earlier TGraphAsymmErrors::Efficiency(..)
0524   // method
0525   //
0526   // see http://root.cern.ch/root/html/TEfficiency.html#compare for
0527   // an illustration of the coverage of the different methods provided by
0528   // TEfficiency
0529
0530   lowerBound = TEfficiency::Wilson(total, passing, level, false);
0531   upperBound = TEfficiency::Wilson(total, passing, level, true);
0532 }
0533
0534 //----------------------------------------------------------------------
0535
0536 DEFINE_FWK_MODULE(EmDQMPostProcessor);
0537
0538 //----------------------------------------------------------------------