CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Validation/RecoEgamma/plugins/PhotonPostprocessing.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 //
00003 
00004 #include "Validation/RecoEgamma/plugins/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   isRunCentrally_=   pset.getParameter<bool>("isRunCentrally");
00040   fastSim_ =   pset.getParameter<bool>("fastSim");
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   std::string photonPathName  = "EgammaV/PhotonValidator/Photons/";
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   if ( ! isRunCentrally_ ) {
00123     histname = "convVsEt";
00124     convVsEt_[0] =  dbe_->book1D(histname+"Barrel","Fraction of good conversions in R9<0.93 vs Et ",etBin,etMin, etMax) ;
00125     convVsEt_[1] =  dbe_->book1D(histname+"Endcap","Fraction of good conversions in R9<0.93 vs Et ",etBin,etMin, etMax) ;
00126   }
00127 
00128 
00129 
00130   // Conversion reconstruction efficiency
00131   histname = "convEffVsEtaTwoTracks";
00132   convEffEtaTwoTracks_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00133 
00134   histname = "convEffVsPhiTwoTracks";
00135   convEffPhiTwoTracks_ =  dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
00136 
00137   histname = "convEffVsRTwoTracks";
00138   convEffRTwoTracks_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00139 
00140   histname = "convEffVsZTwoTracks";
00141   convEffZTwoTracks_ =  dbe_->book1D(histname,histname,zBin,zMin,zMax);
00142 
00143   histname = "convEffVsEtTwoTracks";
00144   convEffEtTwoTracks_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00145   //
00146   histname = "convEffVsEtaTwoTracksAndVtxProbGT0";
00147   convEffEtaTwoTracksAndVtxProbGT0_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00148   histname = "convEffVsEtaTwoTracksAndVtxProbGT0005";
00149   convEffEtaTwoTracksAndVtxProbGT0005_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00150   histname = "convEffVsRTwoTracksAndVtxProbGT0";
00151   convEffRTwoTracksAndVtxProbGT0_ =  dbe_->book1D(histname,histname,rBin,rMin,rMax);
00152   histname = "convEffVsRTwoTracksAndVtxProbGT0005";
00153   convEffRTwoTracksAndVtxProbGT0005_ =  dbe_->book1D(histname,histname,rBin,rMin,rMax);
00154   //
00155   histname = "convEffVsEtaOneTrack";
00156   convEffEtaOneTrack_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00157   histname = "convEffVsROneTrack";
00158   convEffROneTrack_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00159   histname = "convEffVsEtOneTrack";
00160   convEffEtOneTrack_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00161   // Fake rate
00162   histname = "convFakeRateVsEtaTwoTracks";
00163   convFakeRateEtaTwoTracks_ =  dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
00164   histname = "convFakeRateVsPhiTwoTracks";
00165   convFakeRatePhiTwoTracks_ =  dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
00166   histname = "convFakeRateVsRTwoTracks";
00167   convFakeRateRTwoTracks_ =  dbe_->book1D(histname,histname,rBin,rMin, rMax);
00168   histname = "convFakeRateVsZTwoTracks";
00169   convFakeRateZTwoTracks_ =  dbe_->book1D(histname,histname,zBin,zMin,zMax);
00170   histname = "convFakeRateVsEtTwoTracks";
00171   convFakeRateEtTwoTracks_ =  dbe_->book1D(histname,histname,etBin,etMin, etMax);
00172 
00173   histname = "bkgEffVsEta";
00174   bkgRecoEffEta_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #eta",etaBin,etaMin, etaMax);
00175   histname = "bkgEffVsPhi";
00176   bkgRecoEffPhi_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #phi",phiBin,phiMin, phiMax);
00177   histname = "bkgEffVsEt";
00178   bkgRecoEffEt_ =  dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated Et",etBin,etMin, etMax) ;
00179   // Fraction of photons with at least one dead channel
00180   histname = "deadChVsEtaBkg";
00181   bkgDeadChEta_ =  dbe_->book1D(histname,"Fraction of bkg  with >=1 dead Xtal vs simulated #eta",etaBin,etaMin, etaMax);
00182   histname = "deadChVsPhiBkg";
00183   bkgDeadChPhi_ =  dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated #phi",phiBin,phiMin, phiMax);
00184   histname = "deadChVsEtBkg";
00185   bkgDeadChEt_ =  dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated Et",etBin,etMin, etMax) ;
00186 
00187 
00188 
00189   // efficiencies
00190   if ( ! isRunCentrally_ ) {
00191     dividePlots(dbe_->get(effPathName+"convVsEtBarrel"),dbe_->get(photonPathName+"EtR9Less093ConvBarrel"),dbe_->get(photonPathName+"EtR9Less093Barrel"), "effic");
00192     dividePlots(dbe_->get(effPathName+"convVsEtEndcap"),dbe_->get(photonPathName+"EtR9Less093ConvEndcap"),dbe_->get(photonPathName+"EtR9Less093Endcap"), "effic");
00193   }
00194 
00195   dividePlots(dbe_->get(effPathName+"recoEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"),dbe_->get(simInfoPathName+"h_SimPhoEta"), "effic");
00196   dividePlots(dbe_->get(effPathName+"recoEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),dbe_->get(simInfoPathName+"h_SimPhoPhi"),"effic");
00197   dividePlots(dbe_->get(effPathName+"recoEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),dbe_->get(simInfoPathName+"h_SimPhoEt"),"effic");
00198   // fraction of photons with at least one dead channel
00199   dividePlots(dbe_->get(effPathName+"deadChVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"), "effic");
00200   dividePlots(dbe_->get(effPathName+"deadChVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),"effic");
00201   dividePlots(dbe_->get(effPathName+"deadChVsEt"), dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),"effic");
00202   //
00203   if ( ! fastSim_ ) {
00204     dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
00205     dividePlots(dbe_->get(effPathName+"convEffVsPhiTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksPhi"),dbe_->get(simInfoPathName+"h_VisSimConvPhi"),"effic");
00206     dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
00207     dividePlots(dbe_->get(effPathName+"convEffVsZTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksZ"),dbe_->get(simInfoPathName+"h_VisSimConvZ"),"effic");
00208     dividePlots(dbe_->get(effPathName+"convEffVsEtTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
00209     dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
00210     dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
00211     dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
00212     dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
00213     //
00214     dividePlots(dbe_->get(effPathName+"convEffVsEtaOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
00215     dividePlots(dbe_->get(effPathName+"convEffVsROneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
00216     dividePlots(dbe_->get(effPathName+"convEffVsEtOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
00217     // fake rate
00218     dividePlots(dbe_->get(effPathName+"convFakeRateVsEtaTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEta"),dbe_->get(convPathName+"h_RecoConvTwoTracksEta"),"fakerate");
00219     dividePlots(dbe_->get(effPathName+"convFakeRateVsPhiTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksPhi"),dbe_->get(convPathName+"h_RecoConvTwoTracksPhi"),"fakerate");
00220     dividePlots(dbe_->get(effPathName+"convFakeRateVsRTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksR"),dbe_->get(convPathName+"h_RecoConvTwoTracksR"),"fakerate");
00221     dividePlots(dbe_->get(effPathName+"convFakeRateVsZTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksZ"),dbe_->get(convPathName+"h_RecoConvTwoTracksZ"),"fakerate");
00222     dividePlots(dbe_->get(effPathName+"convFakeRateVsEtTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEt"),dbe_->get(convPathName+"h_RecoConvTwoTracksEt"),"fakerate");
00223   }
00224   // Background efficiency
00225   dividePlots(dbe_->get(effPathName+"bkgEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"),dbe_->get(simInfoPathName+"h_SimJetEta"), "effic");
00226   dividePlots(dbe_->get(effPathName+"bkgEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),dbe_->get(simInfoPathName+"h_SimJetPhi"),"effic");
00227   dividePlots(dbe_->get(effPathName+"bkgEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),dbe_->get(simInfoPathName+"h_SimJetEt"),"effic");
00228   // fraction of photons with at least one dead channel
00229   dividePlots(dbe_->get(effPathName+"deadChVsEtaBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"), "effic");
00230   dividePlots(dbe_->get(effPathName+"deadChVsPhiBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),"effic");
00231   dividePlots(dbe_->get(effPathName+"deadChVsEtBkg"), dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),"effic");
00232 
00233 
00234 
00235   if(standAlone_) dbe_->save(outputFileName_);
00236   else if(batch_) dbe_->save(inputFileName_);
00237 
00238 
00239 
00240 }
00241 
00242 
00243 void PhotonPostprocessing::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup)
00244 {
00245 
00246 
00247 }
00248 
00249 
00250 
00251 void  PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, MonitorElement* denominator, std::string type ){
00252   double value,err;
00253   for (int j=1; j<=numerator->getNbinsX(); j++){
00254 
00255     if (denominator->getBinContent(j)!=0){
00256       if (type=="effic")
00257         value = ((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
00258       else if (type=="fakerate")
00259         value = 1-((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
00260       else return;
00261       err = sqrt( value*(1-value) / ((double) denominator->getBinContent(j)) );
00262       dividend->setBinContent(j, value);
00263       if ( err !=0 ) dividend->setBinError(j,err);
00264     }
00265     else {
00266       dividend->setBinContent(j, 0);
00267       dividend->setBinError(j,0);
00268     }
00269 
00270   }
00271 
00272 
00273 }
00274 
00275 
00276 void  PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, double denominator){
00277   double value,err;
00278 
00279   for (int j=1; j<=numerator->getNbinsX(); j++){
00280     if (denominator!=0){
00281       value = ((double) numerator->getBinContent(j))/denominator;
00282       err = sqrt( value*(1-value) / denominator);
00283       dividend->setBinContent(j, value);
00284       dividend->setBinError(j,err);
00285     }
00286     else {
00287       dividend->setBinContent(j, 0);
00288     }
00289   }
00290 
00291 }
00292