RecoEgamma/plugins/ConversionPostprocessing.cc

0001 #include <iostream>
0002 //
0003
0004 #include "Validation/RecoEgamma/plugins/ConversionPostprocessing.h"
0005
0006 //
0007
0008 /** \class ConversionPostprocessing
0009  **
0010  **
0011  **  $Id: ConversionPostprocessing
0012  **  author:
0013  **   Nancy Marinelli, U. of Notre Dame, US
0014  **
0015  **
0016  ***/
0017
0018 using namespace std;
0019
0020 ConversionPostprocessing::ConversionPostprocessing(const edm::ParameterSet& pset) {
0021   dbe_ = nullptr;
0022   dbe_ = edm::Service<DQMStore>().operator->();
0023   parameters_ = pset;
0024
0025   standAlone_ = pset.getParameter<bool>("standAlone");
0026   batch_ = pset.getParameter<bool>("batch");
0027   outputFileName_ = pset.getParameter<string>("OutputFileName");
0028   inputFileName_ = pset.getParameter<std::string>("InputFileName");
0029
0030   etMin = parameters_.getParameter<double>("etMin");
0031   etMax = parameters_.getParameter<double>("etMax");
0032   etBin = parameters_.getParameter<int>("etBin");
0033
0034   etaMin = parameters_.getParameter<double>("etaMin");
0035   etaMax = parameters_.getParameter<double>("etaMax");
0036   etaBin = parameters_.getParameter<int>("etaBin");
0037   etaBin2 = parameters_.getParameter<int>("etaBin2");
0038
0039   phiMin = parameters_.getParameter<double>("phiMin");
0040   phiMax = parameters_.getParameter<double>("phiMax");
0041   phiBin = parameters_.getParameter<int>("phiBin");
0042
0043   rMin = parameters_.getParameter<double>("rMin");
0044   rMax = parameters_.getParameter<double>("rMax");
0045   rBin = parameters_.getParameter<int>("rBin");
0046
0047   zMin = parameters_.getParameter<double>("zMin");
0048   zMax = parameters_.getParameter<double>("zMax");
0049   zBin = parameters_.getParameter<int>("zBin");
0050 }
0051
0052 ConversionPostprocessing::~ConversionPostprocessing() {}
0053
0054 void ConversionPostprocessing::beginJob() {}
0055
0056 void ConversionPostprocessing::analyze(const edm::Event& e, const edm::EventSetup& esup) {}
0057
0058 void ConversionPostprocessing::endJob() {
0059   if (standAlone_)
0060     runPostprocessing();
0061 }
0062
0063 void ConversionPostprocessing::endRun(const edm::Run& run, const edm::EventSetup& setup) {
0064   if (!standAlone_)
0065     runPostprocessing();
0066 }
0067
0068 void ConversionPostprocessing::runPostprocessing() {
0069   std::string simInfoPathName = "EgammaV/ConversionValidator/SimulationInfo/";
0070   std::string convPathName = "EgammaV/ConversionValidator/ConversionInfo/";
0071   std::string effPathName = "EgammaV/ConversionValidator/EfficienciesAndFakeRate/";
0072
0073   if (batch_)
0074     dbe_->open(inputFileName_);
0075
0076   dbe_->setCurrentFolder(effPathName);
0077   // Conversion reconstruction efficiency
0078   std::string histname = "convEffVsEtaTwoTracks";
0079   convEffEtaTwoTracks_ = dbe_->book1D(histname, histname, etaBin2, etaMin, etaMax);
0080
0081   histname = "convEffVsPhiTwoTracks";
0082   convEffPhiTwoTracks_ = dbe_->book1D(histname, histname, phiBin, phiMin, phiMax);
0083
0084   histname = "convEffVsRTwoTracks";
0085   convEffRTwoTracks_ = dbe_->book1D(histname, histname, rBin, rMin, rMax);
0086
0087   histname = "convEffVsZTwoTracks";
0088   convEffZTwoTracks_ = dbe_->book1D(histname, histname, zBin, zMin, zMax);
0089
0090   histname = "convEffVsEtTwoTracks";
0091   convEffEtTwoTracks_ = dbe_->book1D(histname, histname, etBin, etMin, etMax);
0092   //
0093   histname = "convEffVsEtaTwoTracksAndVtxProbGT0";
0094   convEffEtaTwoTracksAndVtxProbGT0_ = dbe_->book1D(histname, histname, etaBin2, etaMin, etaMax);
0095   histname = "convEffVsEtaTwoTracksAndVtxProbGT0005";
0096   convEffEtaTwoTracksAndVtxProbGT0005_ = dbe_->book1D(histname, histname, etaBin2, etaMin, etaMax);
0097   histname = "convEffVsRTwoTracksAndVtxProbGT0";
0098   convEffRTwoTracksAndVtxProbGT0_ = dbe_->book1D(histname, histname, rBin, rMin, rMax);
0099   histname = "convEffVsRTwoTracksAndVtxProbGT0005";
0100   convEffRTwoTracksAndVtxProbGT0005_ = dbe_->book1D(histname, histname, rBin, rMin, rMax);
0101   //
0102   // Fake rate
0103   histname = "convFakeRateVsEtaTwoTracks";
0104   convFakeRateEtaTwoTracks_ = dbe_->book1D(histname, histname, etaBin2, etaMin, etaMax);
0105   histname = "convFakeRateVsPhiTwoTracks";
0106   convFakeRatePhiTwoTracks_ = dbe_->book1D(histname, histname, phiBin, phiMin, phiMax);
0107   histname = "convFakeRateVsRTwoTracks";
0108   convFakeRateRTwoTracks_ = dbe_->book1D(histname, histname, rBin, rMin, rMax);
0109   histname = "convFakeRateVsZTwoTracks";
0110   convFakeRateZTwoTracks_ = dbe_->book1D(histname, histname, zBin, zMin, zMax);
0111   histname = "convFakeRateVsEtTwoTracks";
0112   convFakeRateEtTwoTracks_ = dbe_->book1D(histname, histname, etBin, etMin, etMax);
0113
0114   // efficiencies
0115   dividePlots(dbe_->get(effPathName + "convEffVsEtaTwoTracks"),
0116               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEta"),
0117               dbe_->get(simInfoPathName + "h_VisSimConvEta"),
0118               "effic");
0119   dividePlots(dbe_->get(effPathName + "convEffVsPhiTwoTracks"),
0120               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksPhi"),
0121               dbe_->get(simInfoPathName + "h_VisSimConvPhi"),
0122               "effic");
0123   dividePlots(dbe_->get(effPathName + "convEffVsRTwoTracks"),
0124               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksR"),
0125               dbe_->get(simInfoPathName + "h_VisSimConvR"),
0126               "effic");
0127   dividePlots(dbe_->get(effPathName + "convEffVsZTwoTracks"),
0128               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksZ"),
0129               dbe_->get(simInfoPathName + "h_VisSimConvZ"),
0130               "effic");
0131   dividePlots(dbe_->get(effPathName + "convEffVsEtTwoTracks"),
0132               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEt"),
0133               dbe_->get(simInfoPathName + "h_VisSimConvEt"),
0134               "effic");
0135   dividePlots(dbe_->get(effPathName + "convEffVsEtaTwoTracksAndVtxProbGT0"),
0136               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEtaAndVtxPGT0"),
0137               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEta"),
0138               "effic");
0139   dividePlots(dbe_->get(effPathName + "convEffVsEtaTwoTracksAndVtxProbGT0005"),
0140               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEtaAndVtxPGT0005"),
0141               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksEta"),
0142               "effic");
0143   dividePlots(dbe_->get(effPathName + "convEffVsRTwoTracksAndVtxProbGT0"),
0144               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksRAndVtxPGT0"),
0145               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksR"),
0146               "effic");
0147   dividePlots(dbe_->get(effPathName + "convEffVsRTwoTracksAndVtxProbGT0005"),
0148               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksRAndVtxPGT0005"),
0149               dbe_->get(simInfoPathName + "h_SimConvTwoMTracksR"),
0150               "effic");
0151
0152   // fake rate
0153   dividePlots(dbe_->get(effPathName + "convFakeRateVsEtaTwoTracks"),
0154               dbe_->get(convPathName + "convEtaAss2"),
0155               dbe_->get(convPathName + "convEta2"),
0156               "fakerate");
0157   dividePlots(dbe_->get(effPathName + "convFakeRateVsPhiTwoTracks"),
0158               dbe_->get(convPathName + "convPhiAss"),
0159               dbe_->get(convPathName + "convPhi"),
0160               "fakerate");
0161   dividePlots(dbe_->get(effPathName + "convFakeRateVsRTwoTracks"),
0162               dbe_->get(convPathName + "convRAss"),
0163               dbe_->get(convPathName + "convR"),
0164               "fakerate");
0165   dividePlots(dbe_->get(effPathName + "convFakeRateVsZTwoTracks"),
0166               dbe_->get(convPathName + "convZAss"),
0167               dbe_->get(convPathName + "convZ"),
0168               "fakerate");
0169   dividePlots(dbe_->get(effPathName + "convFakeRateVsEtTwoTracks"),
0170               dbe_->get(convPathName + "convPtAss"),
0171               dbe_->get(convPathName + "convPt"),
0172               "fakerate");
0173
0174   if (standAlone_)
0175     dbe_->save(outputFileName_);
0176   else if (batch_)
0177     dbe_->save(inputFileName_);
0178 }
0179
0180 void ConversionPostprocessing::dividePlots(MonitorElement* dividend,
0181                                            MonitorElement* numerator,
0182                                            MonitorElement* denominator,
0183                                            std::string type) {
0184   double value, err;
0185
0186   //quick fix to avoid seg. faults due to null pointers.
0187   if (dividend == nullptr || numerator == nullptr || denominator == nullptr)
0188     return;
0189
0190   for (int j = 1; j <= numerator->getNbinsX(); j++) {
0191     if (denominator->getBinContent(j) != 0) {
0192       if (type == "effic")
0193         value = ((double)numerator->getBinContent(j)) / ((double)denominator->getBinContent(j));
0194       else if (type == "fakerate")
0195         value = 1 - ((double)numerator->getBinContent(j)) / ((double)denominator->getBinContent(j));
0196       else
0197         return;
0198       err = sqrt(value * (1 - value) / ((double)denominator->getBinContent(j)));
0199       dividend->setBinContent(j, value);
0200       if (err != 0)
0201         dividend->setBinError(j, err);
0202     } else {
0203       dividend->setBinContent(j, 0);
0204       dividend->setBinError(j, 0);
0205     }
0206   }
0207 }
0208
0209 void ConversionPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, double denominator) {
0210   double value, err;
0211
0212   //quick fix to avoid seg. faults due to null pointers.
0213   if (dividend == nullptr || numerator == nullptr)
0214     return;
0215
0216   for (int j = 1; j <= numerator->getNbinsX(); j++) {
0217     if (denominator != 0) {
0218       value = ((double)numerator->getBinContent(j)) / denominator;
0219       err = sqrt(value * (1 - value) / denominator);
0220       dividend->setBinContent(j, value);
0221       dividend->setBinError(j, err);
0222     } else {
0223       dividend->setBinContent(j, 0);
0224     }
0225   }
0226 }