CMS 3D CMS Logo

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