CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Validation/RecoEgamma/src/PhotonPostprocessing.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 //
00003 
00004 #include "Validation/RecoEgamma/interface/PhotonPostprocessing.h"
00005 
00006 
00007 //#define TWOPI 6.283185308
00008 // 
00009 
00023 using namespace std;
00024 
00025  
00026 PhotonPostprocessing::PhotonPostprocessing(const edm::ParameterSet& pset) 
00027 {
00028 
00029   dbe_ = 0;
00030   dbe_ = edm::Service<DQMStore>().operator->();
00031   dbe_->setVerbose(0);
00032   parameters_ = pset;
00033 
00034 
00035   standAlone_ = pset.getParameter<bool>("standAlone");
00036   batch_ = pset.getParameter<bool>("batch");
00037   outputFileName_ = pset.getParameter<string>("OutputFileName");
00038   inputFileName_  = pset.getParameter<std::string>("InputFileName");
00039 
00040 
00041 
00042   etMin = parameters_.getParameter<double>("etMin");
00043   etMax = parameters_.getParameter<double>("etMax");
00044   etBin = parameters_.getParameter<int>("etBin");
00045 
00046 
00047   etaMin = parameters_.getParameter<double>("etaMin");
00048   etaMax = parameters_.getParameter<double>("etaMax");
00049   etaBin = parameters_.getParameter<int>("etaBin");
00050   etaBin2 = parameters_.getParameter<int>("etaBin2");
00051 
00052   phiMin = parameters_.getParameter<double>("phiMin");
00053   phiMax = parameters_.getParameter<double>("phiMax");
00054   phiBin = parameters_.getParameter<int>("phiBin");
00055 
00056   rMin = parameters_.getParameter<double>("rMin");
00057   rMax = parameters_.getParameter<double>("rMax");
00058   rBin = parameters_.getParameter<int>("rBin");
00059 
00060   zMin = parameters_.getParameter<double>("zMin");
00061   zMax = parameters_.getParameter<double>("zMax");
00062   zBin = parameters_.getParameter<int>("zBin");
00063 
00064  
00065 
00066 }
00067 
00068 
00069 
00070 PhotonPostprocessing::~PhotonPostprocessing()
00071 {}
00072 
00073 void PhotonPostprocessing::beginJob()
00074 {
00075 
00076 }
00077 
00078 void PhotonPostprocessing::analyze(const edm::Event& e, const edm::EventSetup& esup)
00079 {}
00080 
00081 
00082 void PhotonPostprocessing::endJob() {
00083 
00084 if(standAlone_) runPostprocessing();
00085 
00086 }
00087 
00088 void PhotonPostprocessing::endRun(const edm::Run& run, const edm::EventSetup& setup) {
00089 
00090   if(!standAlone_) runPostprocessing();
00091 
00092 }
00093 
00094 
00095 void PhotonPostprocessing::runPostprocessing()
00096 {
00097 
00098 
00099   std::string simInfoPathName = "EgammaV/PhotonValidator/SimulationInfo/";
00100   std::string convPathName    = "EgammaV/PhotonValidator/ConversionInfo/";
00101   std::string effPathName     = "EgammaV/PhotonValidator/Efficiencies/";
00102 
00103 
00104   if(batch_)  dbe_->open(inputFileName_);
00105 
00106   dbe_->setCurrentFolder(effPathName);
00107   //  Photon reconstruction efficiencies
00108   string histname = "recoEffVsEta";
00109   phoRecoEffEta_ =  dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated #eta",etaBin,etaMin, etaMax);
00110   histname = "recoEffVsPhi";
00111   phoRecoEffPhi_ =  dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated #phi",phiBin,phiMin, phiMax);
00112   histname = "recoEffVsEt";
00113   phoRecoEffEt_ =  dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated Et",etBin,etMin, etMax) ;
00114   // Fraction of photons with at least one dead channel
00115   histname = "deadChVsEta";
00116   phoDeadChEta_ =  dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated #eta",etaBin,etaMin, etaMax);
00117   histname = "deadChVsPhi";
00118   phoDeadChPhi_ =  dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated #phi",phiBin,phiMin, phiMax);
00119   histname = "deadChVsEt";
00120   phoDeadChEt_ =  dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated Et",etBin,etMin, etMax) ;
00121 
00122   // Conversion reconstruction efficiency
00123   histname = "convEffVsEtaTwoTracks";
00124   convEffEtaTwoTracks_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00125 
00126   histname = "convEffVsPhiTwoTracks";
00127   convEffPhiTwoTracks_ =  dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
00128 
00129   histname = "convEffVsRTwoTracks";
00130   convEffRTwoTracks_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00131 
00132   histname = "convEffVsZTwoTracks";
00133   convEffZTwoTracks_ =  dbe_->book1D(histname,histname,zBin,zMin,zMax);
00134 
00135   histname = "convEffVsEtTwoTracks";
00136   convEffEtTwoTracks_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00137   // 
00138   histname = "convEffVsEtaTwoTracksAndVtxProbGT0";
00139   convEffEtaTwoTracksAndVtxProbGT0_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00140   histname = "convEffVsEtaTwoTracksAndVtxProbGT0005";
00141   convEffEtaTwoTracksAndVtxProbGT0005_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00142   histname = "convEffVsRTwoTracksAndVtxProbGT0";
00143   convEffRTwoTracksAndVtxProbGT0_ =  dbe_->book1D(histname,histname,rBin,rMin,rMax);
00144   histname = "convEffVsRTwoTracksAndVtxProbGT0005";
00145   convEffRTwoTracksAndVtxProbGT0005_ =  dbe_->book1D(histname,histname,rBin,rMin,rMax);
00146   //
00147   histname = "convEffVsEtaOneTrack";
00148   convEffEtaOneTrack_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00149   histname = "convEffVsROneTrack";
00150   convEffROneTrack_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00151   histname = "convEffVsEtOneTrack";
00152   convEffEtOneTrack_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00153   // Fake rate
00154   histname = "convFakeRateVsEtaTwoTracks";
00155   convFakeRateEtaTwoTracks_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00156   histname = "convFakeRateVsPhiTwoTracks";
00157   convFakeRatePhiTwoTracks_ =  dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
00158   histname = "convFakeRateVsRTwoTracks";
00159   convFakeRateRTwoTracks_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00160   histname = "convFakeRateVsZTwoTracks";
00161   convFakeRateZTwoTracks_ =  dbe_->book1D(histname,histname,zBin,zMin,zMax);
00162   histname = "convFakeRateVsEtTwoTracks";
00163   convFakeRateEtTwoTracks_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00164 
00165   histname = "bkgEffVsEta";
00166   bkgRecoEffEta_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #eta",etaBin,etaMin, etaMax);
00167   histname = "bkgEffVsPhi";
00168   bkgRecoEffPhi_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #phi",phiBin,phiMin, phiMax);
00169   histname = "bkgEffVsEt";
00170   bkgRecoEffEt_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated Et",etBin,etMin, etMax) ;
00171   // Fraction of photons with at least one dead channel
00172   histname = "deadChVsEtaBkg";
00173   bkgDeadChEta_ =  dbe_->book1D(histname,"Fraction of bkg  with >=1 dead Xtal vs simulated #eta",etaBin,etaMin, etaMax);
00174   histname = "deadChVsPhiBkg";
00175   bkgDeadChPhi_ =  dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated #phi",phiBin,phiMin, phiMax);
00176   histname = "deadChVsEtBkg";
00177   bkgDeadChEt_ =  dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated Et",etBin,etMin, etMax) ;
00178 
00179   // efficiencies
00180   dividePlots(dbe_->get(effPathName+"recoEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"),dbe_->get(simInfoPathName+"h_SimPhoEta"), "effic");
00181   dividePlots(dbe_->get(effPathName+"recoEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),dbe_->get(simInfoPathName+"h_SimPhoPhi"),"effic");
00182   dividePlots(dbe_->get(effPathName+"recoEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),dbe_->get(simInfoPathName+"h_SimPhoEt"),"effic");
00183   // fraction of photons with at least one dead channel
00184   dividePlots(dbe_->get(effPathName+"deadChVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"), "effic");
00185   dividePlots(dbe_->get(effPathName+"deadChVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),"effic");
00186   dividePlots(dbe_->get(effPathName+"deadChVsEt"), dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),"effic");
00187   //
00188   dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
00189   dividePlots(dbe_->get(effPathName+"convEffVsPhiTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksPhi"),dbe_->get(simInfoPathName+"h_VisSimConvPhi"),"effic");
00190   dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
00191   dividePlots(dbe_->get(effPathName+"convEffVsZTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksZ"),dbe_->get(simInfoPathName+"h_VisSimConvZ"),"effic");
00192   dividePlots(dbe_->get(effPathName+"convEffVsEtTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
00193   dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
00194   dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
00195   dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
00196   dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
00197   //
00198   dividePlots(dbe_->get(effPathName+"convEffVsEtaOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
00199   dividePlots(dbe_->get(effPathName+"convEffVsROneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
00200   dividePlots(dbe_->get(effPathName+"convEffVsEtOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
00201   // fake rate
00202   dividePlots(dbe_->get(effPathName+"convFakeRateVsEtaTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEta"),dbe_->get(convPathName+"h_RecoConvTwoTracksEta"),"fakerate");
00203   dividePlots(dbe_->get(effPathName+"convFakeRateVsPhiTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksPhi"),dbe_->get(convPathName+"h_RecoConvTwoTracksPhi"),"fakerate");
00204   dividePlots(dbe_->get(effPathName+"convFakeRateVsRTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksR"),dbe_->get(convPathName+"h_RecoConvTwoTracksR"),"fakerate");
00205   dividePlots(dbe_->get(effPathName+"convFakeRateVsZTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksZ"),dbe_->get(convPathName+"h_RecoConvTwoTracksZ"),"fakerate");
00206   dividePlots(dbe_->get(effPathName+"convFakeRateVsEtTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEt"),dbe_->get(convPathName+"h_RecoConvTwoTracksEt"),"fakerate");
00207   // Background efficiency
00208   dividePlots(dbe_->get(effPathName+"bkgEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"),dbe_->get(simInfoPathName+"h_SimJetEta"), "effic");
00209   dividePlots(dbe_->get(effPathName+"bkgEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),dbe_->get(simInfoPathName+"h_SimJetPhi"),"effic");
00210   dividePlots(dbe_->get(effPathName+"bkgEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),dbe_->get(simInfoPathName+"h_SimJetEt"),"effic");
00211   // fraction of photons with at least one dead channel
00212   dividePlots(dbe_->get(effPathName+"deadChVsEtaBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"), "effic");
00213   dividePlots(dbe_->get(effPathName+"deadChVsPhiBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),"effic");
00214   dividePlots(dbe_->get(effPathName+"deadChVsEtBkg"), dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),"effic");
00215   
00216 
00217 
00218   if(standAlone_) dbe_->save(outputFileName_);
00219   else if(batch_) dbe_->save(inputFileName_);
00220  
00221   
00222 
00223 }
00224 
00225 
00226 void PhotonPostprocessing::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup)
00227 {
00228 
00229  
00230 }
00231 
00232 
00233 
00234 void  PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, MonitorElement* denominator, std::string type ){
00235   double value,err;
00236   for (int j=1; j<=numerator->getNbinsX(); j++){
00237 
00238     if (denominator->getBinContent(j)!=0){
00239       if (type=="effic")
00240         value = ((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
00241       else if (type=="fakerate")
00242         value = 1-((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
00243       else return; 
00244       err = sqrt( value*(1-value) / ((double) denominator->getBinContent(j)) );
00245       dividend->setBinContent(j, value);
00246       if ( err !=0 ) dividend->setBinError(j,err);
00247     }
00248     else {
00249       dividend->setBinContent(j, 0);
00250       dividend->setBinError(j,0);
00251     }
00252 
00253   }
00254 
00255 
00256 }
00257 
00258 
00259 void  PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, double denominator){
00260   double value,err;
00261 
00262   for (int j=1; j<=numerator->getNbinsX(); j++){
00263     if (denominator!=0){
00264       value = ((double) numerator->getBinContent(j))/denominator;
00265       err = sqrt( value*(1-value) / denominator);
00266       dividend->setBinContent(j, value);
00267       dividend->setBinError(j,err);
00268     }
00269     else {
00270       dividend->setBinContent(j, 0);
00271     }
00272   }
00273 
00274 }
00275