CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PhotonPostprocessing.cc
Go to the documentation of this file.
1 #include <iostream>
2 //
3 
4 #include "Validation/RecoEgamma/plugins/PhotonPostprocessing.h"
5 
6 
7 //#define TWOPI 6.283185308
8 //
9 
23 using namespace std;
24 
25 
26 PhotonPostprocessing::PhotonPostprocessing(const edm::ParameterSet& pset)
27 {
28 
29  dbe_ = 0;
30  dbe_ = edm::Service<DQMStore>().operator->();
31  dbe_->setVerbose(0);
32  parameters_ = pset;
33 
34 
35  analyzerName_ = pset.getParameter<std::string>("analyzerName");
36  standAlone_ = pset.getParameter<bool>("standAlone");
37  batch_ = pset.getParameter<bool>("batch");
38  outputFileName_ = pset.getParameter<string>("OutputFileName");
39  inputFileName_ = pset.getParameter<std::string>("InputFileName");
40  isRunCentrally_= pset.getParameter<bool>("isRunCentrally");
41  fastSim_ = pset.getParameter<bool>("fastSim");
42 
43  etMin = parameters_.getParameter<double>("etMin");
44  etMax = parameters_.getParameter<double>("etMax");
45  etBin = parameters_.getParameter<int>("etBin");
46 
47 
48  etaMin = parameters_.getParameter<double>("etaMin");
49  etaMax = parameters_.getParameter<double>("etaMax");
50  etaBin = parameters_.getParameter<int>("etaBin");
51  etaBin2 = parameters_.getParameter<int>("etaBin2");
52 
53  phiMin = parameters_.getParameter<double>("phiMin");
54  phiMax = parameters_.getParameter<double>("phiMax");
55  phiBin = parameters_.getParameter<int>("phiBin");
56 
57  rMin = parameters_.getParameter<double>("rMin");
58  rMax = parameters_.getParameter<double>("rMax");
59  rBin = parameters_.getParameter<int>("rBin");
60 
61  zMin = parameters_.getParameter<double>("zMin");
62  zMax = parameters_.getParameter<double>("zMax");
63  zBin = parameters_.getParameter<int>("zBin");
64 
65 
66 
67 }
68 
69 
70 
71 PhotonPostprocessing::~PhotonPostprocessing()
72 {}
73 
74 void PhotonPostprocessing::beginJob()
75 {
76 
77 }
78 
79 void PhotonPostprocessing::analyze(const edm::Event& e, const edm::EventSetup& esup)
80 {}
81 
82 
83 void PhotonPostprocessing::endJob() {
84 
85 if(standAlone_) runPostprocessing();
86 
87 }
88 
89 void PhotonPostprocessing::endRun(const edm::Run& run, const edm::EventSetup& setup) {
90 
91  if(!standAlone_) runPostprocessing();
92 
93 }
94 
95 
96 void PhotonPostprocessing::runPostprocessing()
97 {
98 
99 
100  std::string simInfoPathName = "EgammaV/"+ analyzerName_+"/SimulationInfo/";
101  std::string convPathName = "EgammaV/"+ analyzerName_+"/ConversionInfo/";
102  std::string effPathName = "EgammaV/"+ analyzerName_+"/Efficiencies/";
103  std::string photonPathName = "EgammaV/"+ analyzerName_+"/Photons/";
104 
105  if(batch_) dbe_->open(inputFileName_);
106 
107  dbe_->setCurrentFolder(effPathName);
108  // Photon reconstruction efficiencies
109  string histname = "recoEffVsEta";
110  phoRecoEffEta_ = dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated #eta",etaBin,etaMin, etaMax);
111 
112  histname = "recoEffVsPhi";
113  phoRecoEffPhi_ = dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated #phi",phiBin,phiMin, phiMax);
114  histname = "recoEffVsEt";
115  phoRecoEffEt_ = dbe_->book1D(histname,"Photon reconstruction efficiency vs simulated Et",etBin,etMin, etMax) ;
116  // Fraction of photons with at least one dead channel
117  histname = "deadChVsEta";
118  phoDeadChEta_ = dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated #eta",etaBin,etaMin, etaMax);
119  histname = "deadChVsPhi";
120  phoDeadChPhi_ = dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated #phi",phiBin,phiMin, phiMax);
121  histname = "deadChVsEt";
122  phoDeadChEt_ = dbe_->book1D(histname,"Fraction of photons with >=1 dead Xtal vs simulated Et",etBin,etMin, etMax) ;
123 
124  if ( ! isRunCentrally_ ) {
125  histname = "convVsEt";
126  convVsEt_[0] = dbe_->book1D(histname+"Barrel","Fraction of good conversions in R9<0.93 vs Et ",etBin,etMin, etMax) ;
127  convVsEt_[1] = dbe_->book1D(histname+"Endcap","Fraction of good conversions in R9<0.93 vs Et ",etBin,etMin, etMax) ;
128  }
129 
130 
131 
132  // Conversion reconstruction efficiency
133  histname = "convEffVsEtaTwoTracks";
134  convEffEtaTwoTracks_ = dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
135 
136  histname = "convEffVsPhiTwoTracks";
137  convEffPhiTwoTracks_ = dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
138 
139  histname = "convEffVsRTwoTracks";
140  convEffRTwoTracks_ = dbe_->book1D(histname,histname,rBin,rMin, rMax);
141 
142  histname = "convEffVsZTwoTracks";
143  convEffZTwoTracks_ = dbe_->book1D(histname,histname,zBin,zMin,zMax);
144 
145  histname = "convEffVsEtTwoTracks";
146  convEffEtTwoTracks_ = dbe_->book1D(histname,histname,etBin,etMin, etMax);
147  //
148  histname = "convEffVsEtaTwoTracksAndVtxProbGT0";
149  convEffEtaTwoTracksAndVtxProbGT0_ = dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
150  histname = "convEffVsEtaTwoTracksAndVtxProbGT0005";
151  convEffEtaTwoTracksAndVtxProbGT0005_ = dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
152  histname = "convEffVsRTwoTracksAndVtxProbGT0";
153  convEffRTwoTracksAndVtxProbGT0_ = dbe_->book1D(histname,histname,rBin,rMin,rMax);
154  histname = "convEffVsRTwoTracksAndVtxProbGT0005";
155  convEffRTwoTracksAndVtxProbGT0005_ = dbe_->book1D(histname,histname,rBin,rMin,rMax);
156  //
157  histname = "convEffVsEtaOneTrack";
158  convEffEtaOneTrack_ = dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
159  histname = "convEffVsROneTrack";
160  convEffROneTrack_ = dbe_->book1D(histname,histname,rBin,rMin, rMax);
161  histname = "convEffVsEtOneTrack";
162  convEffEtOneTrack_ = dbe_->book1D(histname,histname,etBin,etMin, etMax);
163  // Fake rate
164  histname = "convFakeRateVsEtaTwoTracks";
165  convFakeRateEtaTwoTracks_ = dbe_->book1D(histname,histname,etaBin2,etaMin, etaMax);
166  histname = "convFakeRateVsPhiTwoTracks";
167  convFakeRatePhiTwoTracks_ = dbe_->book1D(histname,histname,phiBin,phiMin,phiMax);
168  histname = "convFakeRateVsRTwoTracks";
169  convFakeRateRTwoTracks_ = dbe_->book1D(histname,histname,rBin,rMin, rMax);
170  histname = "convFakeRateVsZTwoTracks";
171  convFakeRateZTwoTracks_ = dbe_->book1D(histname,histname,zBin,zMin,zMax);
172  histname = "convFakeRateVsEtTwoTracks";
173  convFakeRateEtTwoTracks_ = dbe_->book1D(histname,histname,etBin,etMin, etMax);
174 
175  histname = "bkgEffVsEta";
176  bkgRecoEffEta_ = dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #eta",etaBin,etaMin, etaMax);
177  histname = "bkgEffVsPhi";
178  bkgRecoEffPhi_ = dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated #phi",phiBin,phiMin, phiMax);
179  histname = "bkgEffVsEt";
180  bkgRecoEffEt_ = dbe_->book1D(histname,"Bkg reconstruction efficiency vs simulated Et",etBin,etMin, etMax) ;
181  // Fraction of photons with at least one dead channel
182  histname = "deadChVsEtaBkg";
183  bkgDeadChEta_ = dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated #eta",etaBin,etaMin, etaMax);
184  histname = "deadChVsPhiBkg";
185  bkgDeadChPhi_ = dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated #phi",phiBin,phiMin, phiMax);
186  histname = "deadChVsEtBkg";
187  bkgDeadChEt_ = dbe_->book1D(histname,"Fraction of bkg with >=1 dead Xtal vs simulated Et",etBin,etMin, etMax) ;
188 
189 
190 
191  // efficiencies
192  if ( ! isRunCentrally_ ) {
193  dividePlots(dbe_->get(effPathName+"convVsEtBarrel"),dbe_->get(photonPathName+"EtR9Less093ConvBarrel"),dbe_->get(photonPathName+"EtR9Less093Barrel"), "effic");
194  dividePlots(dbe_->get(effPathName+"convVsEtEndcap"),dbe_->get(photonPathName+"EtR9Less093ConvEndcap"),dbe_->get(photonPathName+"EtR9Less093Endcap"), "effic");
195  }
196 
197  dividePlots(dbe_->get(effPathName+"recoEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"),dbe_->get(simInfoPathName+"h_SimPhoEta"), "effic");
198  dividePlots(dbe_->get(effPathName+"recoEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),dbe_->get(simInfoPathName+"h_SimPhoPhi"),"effic");
199  dividePlots(dbe_->get(effPathName+"recoEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),dbe_->get(simInfoPathName+"h_SimPhoEt"),"effic");
200  // fraction of photons with at least one dead channel
201  dividePlots(dbe_->get(effPathName+"deadChVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEta"), "effic");
202  dividePlots(dbe_->get(effPathName+"deadChVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimPhoPhi"),"effic");
203  dividePlots(dbe_->get(effPathName+"deadChVsEt"), dbe_->get(simInfoPathName+"h_MatchedSimPhoBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimPhoEt"),"effic");
204  //
205  if ( ! fastSim_ ) {
206  dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
207  dividePlots(dbe_->get(effPathName+"convEffVsPhiTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksPhi"),dbe_->get(simInfoPathName+"h_VisSimConvPhi"),"effic");
208  dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
209  dividePlots(dbe_->get(effPathName+"convEffVsZTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksZ"),dbe_->get(simInfoPathName+"h_VisSimConvZ"),"effic");
210  dividePlots(dbe_->get(effPathName+"convEffVsEtTwoTracks"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
211  dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
212  dividePlots(dbe_->get(effPathName+"convEffVsEtaTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEtaAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksEta"),"effic");
213  dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
214  dividePlots(dbe_->get(effPathName+"convEffVsRTwoTracksAndVtxProbGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksRAndVtxPGT0005"),dbe_->get(simInfoPathName+"h_SimConvTwoMTracksR"),"effic");
215  //
216  dividePlots(dbe_->get(effPathName+"convEffVsEtaOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEta"),dbe_->get(simInfoPathName+"h_VisSimConvEta"),"effic");
217  dividePlots(dbe_->get(effPathName+"convEffVsROneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksR"),dbe_->get(simInfoPathName+"h_VisSimConvR"),"effic");
218  dividePlots(dbe_->get(effPathName+"convEffVsEtOneTrack"),dbe_->get(simInfoPathName+"h_SimConvOneMTracksEt"),dbe_->get(simInfoPathName+"h_VisSimConvEt"),"effic");
219  // fake rate
220  dividePlots(dbe_->get(effPathName+"convFakeRateVsEtaTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEta"),dbe_->get(convPathName+"h_RecoConvTwoTracksEta"),"fakerate");
221  dividePlots(dbe_->get(effPathName+"convFakeRateVsPhiTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksPhi"),dbe_->get(convPathName+"h_RecoConvTwoTracksPhi"),"fakerate");
222  dividePlots(dbe_->get(effPathName+"convFakeRateVsRTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksR"),dbe_->get(convPathName+"h_RecoConvTwoTracksR"),"fakerate");
223  dividePlots(dbe_->get(effPathName+"convFakeRateVsZTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksZ"),dbe_->get(convPathName+"h_RecoConvTwoTracksZ"),"fakerate");
224  dividePlots(dbe_->get(effPathName+"convFakeRateVsEtTwoTracks"),dbe_->get(convPathName+"h_RecoConvTwoMTracksEt"),dbe_->get(convPathName+"h_RecoConvTwoTracksEt"),"fakerate");
225  }
226  // Background efficiency
227  dividePlots(dbe_->get(effPathName+"bkgEffVsEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"),dbe_->get(simInfoPathName+"h_SimJetEta"), "effic");
228  dividePlots(dbe_->get(effPathName+"bkgEffVsPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),dbe_->get(simInfoPathName+"h_SimJetPhi"),"effic");
229  dividePlots(dbe_->get(effPathName+"bkgEffVsEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),dbe_->get(simInfoPathName+"h_SimJetEt"),"effic");
230  // fraction of photons with at least one dead channel
231  dividePlots(dbe_->get(effPathName+"deadChVsEtaBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEta"),dbe_->get(simInfoPathName+"h_MatchedSimJetEta"), "effic");
232  dividePlots(dbe_->get(effPathName+"deadChVsPhiBkg"),dbe_->get(simInfoPathName+"h_MatchedSimJetBadChPhi"),dbe_->get(simInfoPathName+"h_MatchedSimJetPhi"),"effic");
233  dividePlots(dbe_->get(effPathName+"deadChVsEtBkg"), dbe_->get(simInfoPathName+"h_MatchedSimJetBadChEt"),dbe_->get(simInfoPathName+"h_MatchedSimJetEt"),"effic");
234 
235 
236 
237  if(standAlone_) dbe_->save(outputFileName_);
238  else if(batch_) dbe_->save(inputFileName_);
239 
240 
241 
242 }
243 
244 
245 void PhotonPostprocessing::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup)
246 {
247 
248 
249 }
250 
251 
252 
253 void PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, MonitorElement* denominator, std::string type ){
254  double value,err;
255 
256  for (int j=1; j<=numerator->getNbinsX(); j++){
257  dividend->setEfficiencyFlag();
258  if (denominator->getBinContent(j)!=0){
259  if (type=="effic")
260  value = ((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
261  else if (type=="fakerate")
262  value = 1-((double) numerator->getBinContent(j))/((double) denominator->getBinContent(j));
263  else return;
264  err = sqrt( value*(1-value) / ((double) denominator->getBinContent(j)) );
265  dividend->setBinContent(j, value);
266  if ( err !=0 ) dividend->setBinError(j,err);
267  }
268  else {
269  dividend->setBinContent(j, 0);
270  dividend->setBinError(j,0);
271  }
272 
273  }
274 
275 
276 }
277 
278 
279 void PhotonPostprocessing::dividePlots(MonitorElement* dividend, MonitorElement* numerator, double denominator){
280  double value,err;
281 
282  for (int j=1; j<=numerator->getNbinsX(); j++){
283  if (denominator!=0){
284  value = ((double) numerator->getBinContent(j))/denominator;
285  err = sqrt( value*(1-value) / denominator);
286  dividend->setBinContent(j, value);
287  dividend->setBinError(j,err);
288  }
289  else {
290  dividend->setBinContent(j, 0);
291  }
292  }
293 
294 }
295 
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
MonitorElement::setBinContent
void setBinContent(int binx, double content)
set content of bin (1-D)
Definition: MonitorElement.cc:813
cuy.numerator
list numerator
Definition: cuy.py:483
DQMStore::book1D
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:722
fjr2json.lumi
tuple lumi
Definition: fjr2json.py:35
DQMStore::save
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:2118
jptDQMConfig_cff.etMax
int etMax
Definition: jptDQMConfig_cff.py:19
jptDQMConfig_cff.phiMax
float phiMax
Definition: jptDQMConfig_cff.py:7
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
cuy.denominator
list denominator
Definition: cuy.py:484
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
j
int j
Definition: DBlmapReader.cc:9
MonitorElement::setBinError
void setBinError(int binx, double error)
set uncertainty on content of bin (1-D)
Definition: MonitorElement.cc:838
DQMStore::setVerbose
void setVerbose(unsigned level)
Definition: DQMStore.cc:398
DQMStore::get
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1473
PhotonPostprocessing::PhotonPostprocessing
PhotonPostprocessing(const edm::ParameterSet &pset)
Definition: PhotonPostprocessing.cc:26
dbe_
DQMStore * dbe_
Definition: PFJetBenchmarkAnalyzer.cc:78
reco::tau::qcuts::etMin
bool etMin(const PFCandidate &cand, double cut)
Definition: RecoTauQualityCuts.cc:36
jptDQMConfig_cff.etaMin
int etaMin
Definition: jptDQMConfig_cff.py:52
PhotonPostprocessing::runPostprocessing
virtual void runPostprocessing()
Definition: PhotonPostprocessing.cc:96
MonitorElement::getBinContent
double getBinContent(int binx) const
get content of bin (1-D)
Definition: MonitorElement.cc:712
jptDQMConfig_cff.etaMax
float etaMax
Definition: jptDQMConfig_cff.py:9
PhotonPostprocessing::dividePlots
void dividePlots(MonitorElement *dividend, MonitorElement *numerator, MonitorElement *denominator, std::string type)
Definition: PhotonPostprocessing.cc:253
MonitorElement::getNbinsX
int getNbinsX(void) const
get # of bins in X-axis
Definition: MonitorElement.cc:694
jptDQMConfig_cff.phiMin
int phiMin
Definition: jptDQMConfig_cff.py:50
DQMStore::open
bool open(const std::string &filename, bool overwrite=false, const std::string &path="", const std::string &prepend="", OpenRunDirs stripdirs=KeepRunDirs, bool fileMustExist=true)
Definition: DQMStore.cc:2437
PhotonPostprocessing::endRun
virtual void endRun(const edm::Run &, const edm::EventSetup &)
Definition: PhotonPostprocessing.cc:89
MonitorElement::setEfficiencyFlag
void setEfficiencyFlag(void)
Definition: MonitorElement.h:141
PhotonPostprocessing::endLuminosityBlock
virtual void endLuminosityBlock(const edm::LuminosityBlock &, const edm::EventSetup &)
Definition: PhotonPostprocessing.cc:245
HcalObjRepresent::setup
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Definition: HcalObjRepresent.h:628
DQMStore::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:434
PhotonPostprocessing::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PhotonPostprocessing.cc:79
edm::Run
Definition: Run.h:36
Generated for CMSSW Reference Manual by doxygen 1.8.5