CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GenericBenchmark.cc
Go to the documentation of this file.
4 
5 // CMSSW_2_X_X
7 
8 #include <TROOT.h>
9 #include <TFile.h>
10 
11 //Colin: it seems a bit strange to use the preprocessor for that kind of
12 //thing. Looks like all these macros could be replaced by plain functions.
13 
14 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
15 #define BOOK1D(name,title,nbinsx,lowx,highx) \
16  h##name = DQM ? DQM->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
17  : new TH1F(#name,title,nbinsx,lowx,highx)
18 
19 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
20 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
21  h##name = DQM ? DQM->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
22  : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
23 
24 
25 // all versions OK
26 // preprocesor macro for setting axis titles
27 
28 #define SETAXES(name,xtitle,ytitle) \
29  h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
30 
31 #define ET (PlotAgainstReco_)?"reconstructed E_{T} [GeV]":"generated E_{T} [GeV]"
32 #define ETA (PlotAgainstReco_)?"reconstructed #eta":"generated #eta"
33 #define PHI (PlotAgainstReco_)?"reconstructed #phi (rad)":"generated #phi (rad)"
34 
35 using namespace std;
36 
38 
40 
41 void GenericBenchmark::setup(DQMStore *DQM, bool PlotAgainstReco_, float minDeltaEt, float maxDeltaEt,
42  float minDeltaPhi, float maxDeltaPhi, bool doMetPlots) {
43 
44  //std::cout << "minDeltaPhi = " << minDeltaPhi << std::endl;
45 
46  // CMSSW_2_X_X
47  // use bare Root if no DQM (FWLite applications)
48  //if (!DQM)
49  //{
50  // file_ = new TFile("pfmetBenchmark.root", "recreate");
51  // cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
52  //}
53  // Book Histograms
54 
55  //std::cout << "FL : pwd = ";
56  //gDirectory->pwd();
57  //std::cout << std::endl;
58 
59  int nbinsEt = 1000;
60  float minEt = 0;
61  float maxEt = 1000;
62 
63  //float minDeltaEt = -100;
64  //float maxDeltaEt = 50;
65 
66  int nbinsEta = 200;
67  float minEta = -5;
68  float maxEta = 5;
69 
70  int nbinsDeltaEta = 1000;
71  float minDeltaEta = -0.5;
72  float maxDeltaEta = 0.5;
73 
74  int nbinsDeltaPhi = 1000;
75  //float minDeltaPhi = -0.5;
76  //float maxDeltaPhi = 0.5;
77 
78 
79  // delta et quantities
80  BOOK1D(DeltaEt,"#DeltaE_{T}", nbinsEt, minDeltaEt, maxDeltaEt);
81  BOOK1D(DeltaEx,"#DeltaE_{X}", nbinsEt, minDeltaEt, maxDeltaEt);
82  BOOK1D(DeltaEy,"#DeltaE_{Y}", nbinsEt, minDeltaEt, maxDeltaEt);
83  BOOK2D(DeltaEtvsEt,"#DeltaE_{T} vs E_{T}",
84  nbinsEt, minEt, maxEt,
85  nbinsEt,minDeltaEt, maxDeltaEt);
86  BOOK2D(DeltaEtOverEtvsEt,"#DeltaE_{T}/E_{T} vsE_{T}",
87  nbinsEt, minEt, maxEt,
88  nbinsEt,-1,1);
89  BOOK2D(DeltaEtvsEta,"#DeltaE_{T} vs #eta",
90  nbinsEta, minEta, maxEta,
91  nbinsEt,minDeltaEt, maxDeltaEt);
92  BOOK2D(DeltaEtOverEtvsEta,"#DeltaE_{T}/E_{T} vs #eta",
93  nbinsEta, minEta, maxEta,
94  100,-1,1);
95  BOOK2D(DeltaEtvsPhi,"#DeltaE_{T} vs #phi",
96  200,-M_PI,M_PI,
97  nbinsEt,minDeltaEt, maxDeltaEt);
98  BOOK2D(DeltaEtOverEtvsPhi,"#DeltaE_{T}/E_{T} vs #Phi",
99  200,-M_PI,M_PI,
100  100,-1,1);
101  BOOK2D(DeltaEtvsDeltaR,"#DeltaE_{T} vs #DeltaR",
102  100,0,1,
103  nbinsEt,minDeltaEt, maxDeltaEt);
104  BOOK2D(DeltaEtOverEtvsDeltaR,"#DeltaE_{T}/E_{T} vs #DeltaR",
105  100,0,1,
106  100,-1,1);
107 
108  // delta eta quantities
109  BOOK1D(DeltaEta,"#Delta#eta",nbinsDeltaEta,minDeltaEta,maxDeltaEta);
110  BOOK2D(DeltaEtavsEt,"#Delta#eta vs E_{T}",
111  nbinsEt, minEt, maxEt,
112  nbinsDeltaEta,minDeltaEta,maxDeltaEta);
113  BOOK2D(DeltaEtavsEta,"#Delta#eta vs #eta",
114  nbinsEta, minEta, maxEta,
115  nbinsDeltaEta,minDeltaEta,maxDeltaEta);
116 
117  // delta phi quantities
118  BOOK1D(DeltaPhi,"#Delta#phi",nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
119  BOOK2D(DeltaPhivsEt,"#Delta#phi vs E_{T}",
120  nbinsEt, minEt, maxEt,
121  nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
122  BOOK2D(DeltaPhivsEta,"#Delta#phi vs #eta",
123  nbinsEta, minEta, maxEta,
124  nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
125 
126  // delta R quantities
127  BOOK1D(DeltaR,"#DeltaR",100,0,1);
128  BOOK2D(DeltaRvsEt,"#DeltaR vs E_{T}",
129  nbinsEt, minEt, maxEt,
130  100,0,1);
131  BOOK2D(DeltaRvsEta,"#DeltaR vs #eta",
132  nbinsEta, minEta, maxEta,
133  100,0,1);
134 
135  BOOK1D(NRec,"Number of reconstructed objects",20,0,20);
136 
137  // seen and gen distributions, for efficiency computation
138  BOOK1D(EtaSeen,"seen #eta",100,-5,5);
139  BOOK1D(PhiSeen,"seen #phi",100,-3.5,3.5);
140  BOOK1D(EtSeen,"seen E_{T}",nbinsEt, minEt, maxEt);
141  BOOK2D(EtvsEtaSeen,"seen E_{T} vs eta",100,-5,5,200, 0, 200);
142  BOOK2D(EtvsPhiSeen,"seen E_{T} vs seen #phi",100,-3.5,3.5,200, 0, 200);
143 
144  BOOK1D(PhiRec,"Rec #phi",100,-3.5,3.5);
145  BOOK1D(EtRec,"Rec E_{T}",nbinsEt, minEt, maxEt);
146  BOOK1D(ExRec,"Rec E_{X}",nbinsEt, -maxEt, maxEt);
147  BOOK1D(EyRec,"Rec E_{Y}",nbinsEt, -maxEt, maxEt);
148 
149  BOOK2D(EtRecvsEt,"Rec E_{T} vs E_{T}",
150  nbinsEt, minEt, maxEt,
151  nbinsEt, minEt, maxEt);
152  BOOK2D(EtRecOverTrueEtvsTrueEt,"Rec E_{T} / E_{T} vs E_{T}",
153  nbinsEt, minEt, maxEt,
154  1000, 0., 100.);
155 
156  BOOK1D(EtaGen,"generated #eta",100,-5,5);
157  BOOK1D(PhiGen,"generated #phi",100,-3.5,3.5);
158  BOOK1D(EtGen,"generated E_{T}",nbinsEt, minEt, maxEt);
159  BOOK2D(EtvsEtaGen,"generated E_{T} vs generated #eta",100,-5,5,200, 0, 200);
160  BOOK2D(EtvsPhiGen,"generated E_{T} vs generated #phi",100,-3.5,3.5, 200, 0, 200);
161 
162  BOOK1D(NGen,"Number of generated objects",20,0,20);
163 
164  if (doMetPlots)
165  {
166  BOOK1D(SumEt,"SumEt", 1000, 0., 3000.);
167  BOOK1D(TrueSumEt,"TrueSumEt", 1000, 0., 3000.);
168  BOOK2D(DeltaSetvsSet,"#DeltaSEt vs trueSEt",
169  3000, 0., 3000.,
170  1000,-1000., 1000.);
171  BOOK2D(DeltaMexvsSet,"#DeltaMEX vs trueSEt",
172  3000, 0., 3000.,
173  1000,-400., 400.);
174  BOOK2D(DeltaSetOverSetvsSet,"#DeltaSetOverSet vs trueSet",
175  3000, 0., 3000.,
176  1000,-1., 1.);
177  BOOK2D(RecSetvsTrueSet,"Set vs trueSet",
178  3000, 0., 3000.,
179  3000,0., 3000.);
180  BOOK2D(RecSetOverTrueSetvsTrueSet,"Set/trueSet vs trueSet",
181  3000, 0., 3000.,
182  500,0., 2.);
183  BOOK2D(TrueMexvsTrueSet,"trueMex vs trueSet",
184  3000,0., 3000.,
185  nbinsEt, -maxEt, maxEt);
186  }
187  // number of truth particles found within given cone radius of reco
188  //BOOK2D(NumInConeVsConeSize,NumInConeVsConeSize,100,0,1,25,0,25);
189 
190  // Set Axis Titles
191 
192  // delta et quantities
193  SETAXES(DeltaEt,"#DeltaE_{T} [GeV]","");
194  SETAXES(DeltaEx,"#DeltaE_{X} [GeV]","");
195  SETAXES(DeltaEy,"#DeltaE_{Y} [GeV]","");
196  SETAXES(DeltaEtvsEt,ET,"#DeltaE_{T} [GeV]");
197  SETAXES(DeltaEtOverEtvsEt,ET,"#DeltaE_{T}/E_{T}");
198  SETAXES(DeltaEtvsEta,ETA,"#DeltaE_{T} [GeV]");
199  SETAXES(DeltaEtOverEtvsEta,ETA,"#DeltaE_{T}/E_{T}");
200  SETAXES(DeltaEtvsPhi,PHI,"#DeltaE_{T} [GeV]");
201  SETAXES(DeltaEtOverEtvsPhi,PHI,"#DeltaE_{T}/E_{T}");
202  SETAXES(DeltaEtvsDeltaR,"#DeltaR","#DeltaE_{T} [GeV]");
203  SETAXES(DeltaEtOverEtvsDeltaR,"#DeltaR","#DeltaE_{T}/E_{T}");
204 
205  // delta eta quantities
206  SETAXES(DeltaEta,"#Delta#eta","Events");
207  SETAXES(DeltaEtavsEt,ET,"#Delta#eta");
208  SETAXES(DeltaEtavsEta,ETA,"#Delta#eta");
209 
210  // delta phi quantities
211  SETAXES(DeltaPhi,"#Delta#phi [rad]","Events");
212  SETAXES(DeltaPhivsEt,ET,"#Delta#phi [rad]");
213  SETAXES(DeltaPhivsEta,ETA,"#Delta#phi [rad]");
214 
215  // delta R quantities
216  SETAXES(DeltaR,"#DeltaR","Events");
217  SETAXES(DeltaRvsEt,ET,"#DeltaR");
218  SETAXES(DeltaRvsEta,ETA,"#DeltaR");
219 
220  SETAXES(NRec,"Number of Rec Objects","");
221 
222  SETAXES(EtaSeen,"seen #eta","");
223  SETAXES(PhiSeen,"seen #phi [rad]","");
224  SETAXES(EtSeen,"seen E_{T} [GeV]","");
225  SETAXES(EtvsEtaSeen,"seen #eta","seen E_{T}");
226  SETAXES(EtvsPhiSeen,"seen #phi [rad]","seen E_{T}");
227 
228  SETAXES(PhiRec,"#phi [rad]","");
229  SETAXES(EtRec,"E_{T} [GeV]","");
230  SETAXES(ExRec,"E_{X} [GeV]","");
231  SETAXES(EyRec,"E_{Y} [GeV]","");
232  SETAXES(EtRecvsEt,ET,"Rec E_{T} [GeV]");
233  SETAXES(EtRecOverTrueEtvsTrueEt,ET,"Rec E_{T} / E_{T} [GeV]");
234 
235  SETAXES(EtaGen,"generated #eta","");
236  SETAXES(PhiGen,"generated #phi [rad]","");
237  SETAXES(EtGen,"generated E_{T} [GeV]","");
238  SETAXES(EtvsPhiGen,"generated #phi [rad]","generated E_{T} [GeV]");
239  SETAXES(EtvsEtaGen,"generated #eta","generated E_{T} [GeV]");
240 
241  SETAXES(NGen,"Number of Gen Objects","");
242 
243  if (doMetPlots)
244  {
245  SETAXES(SumEt,"SumEt [GeV]","");
246  SETAXES(TrueSumEt,"TrueSumEt [GeV]","");
247  SETAXES(DeltaSetvsSet,"TrueSumEt","#DeltaSumEt [GeV]");
248  SETAXES(DeltaMexvsSet,"TrueSumEt","#DeltaMEX [GeV]");
249  SETAXES(DeltaSetOverSetvsSet,"TrueSumEt","#DeltaSumEt/trueSumEt");
250  SETAXES(RecSetvsTrueSet,"TrueSumEt","SumEt");
251  SETAXES(RecSetOverTrueSetvsTrueSet,"TrueSumEt","SumEt/trueSumEt");
252  SETAXES(TrueMexvsTrueSet,"TrueSumEt","TrueMEX");
253  }
254 
255  TDirectory* oldpwd = gDirectory;
256 
257 
258  TIter next( gROOT->GetListOfFiles() );
259 
260  const bool debug=false;
261 
262  while ( TFile *file = (TFile *)next() )
263  {
264  if (debug) cout<<"file "<<file->GetName()<<endl;
265  }
266  if (DQM)
267  {
268  cout<<"DQM subdir"<<endl;
269  cout<< DQM->pwd().c_str()<<endl;
270 
271  DQM->cd( DQM->pwd() );
272  }
273 
274  if (debug)
275  {
276  cout<<"current dir"<<endl;
277  gDirectory->pwd();
278  }
279 
280  oldpwd->cd();
281  //gDirectory->pwd();
282 
283  doMetPlots_=doMetPlots;
284 
285  // tree_ = new BenchmarkTree("genericBenchmark", "Generic Benchmark TTree");
286 }
287 
289  double ptCut,
290  double minEtaCut,
291  double maxEtaCut ) const {
292 
293  //skip reconstructed PFJets with p_t < recPt_cut
294  if (particle->pt() < ptCut and ptCut != -1.)
295  return false;
296 
297  if (fabs(particle->eta())>maxEtaCut and maxEtaCut > 0)
298  return false;
299  if (fabs(particle->eta())<minEtaCut and minEtaCut > 0)
300  return false;
301 
302  //accepted!
303  return true;
304 
305 }
306 
307 
309  const reco::Candidate* recParticle,
310  double deltaR_cut,
311  bool plotAgainstReco ) {
312 
313  // get the quantities to place on the denominator and/or divide by
314  double et = genParticle->et();
315  double eta = genParticle->eta();
316  double phi = genParticle->phi();
317  //std::cout << "FL : et = " << et << std::endl;
318  //std::cout << "FL : eta = " << eta << std::endl;
319  //std::cout << "FL : phi = " << phi << std::endl;
320  //std::cout << "FL : rec et = " << recParticle->et() << std::endl;
321  //std::cout << "FL : rec eta = " << recParticle->eta() << std::endl;
322  //std::cout << "FL : rec phi = " <<recParticle-> phi() << std::endl;
323 
324  if (plotAgainstReco) {
325  et = recParticle->et();
326  eta = recParticle->eta();
327  phi = recParticle->phi();
328  }
329 
330 
331  // get the delta quantities
332  double deltaEt = algo_->deltaEt(recParticle,genParticle);
333  double deltaR = algo_->deltaR(recParticle,genParticle);
334  double deltaEta = algo_->deltaEta(recParticle,genParticle);
335  double deltaPhi = algo_->deltaPhi(recParticle,genParticle);
336 
337  //std::cout << "FL :deltaR_cut = " << deltaR_cut << std::endl;
338  //std::cout << "FL :deltaR = " << deltaR << std::endl;
339 
340  if (deltaR>deltaR_cut && deltaR_cut != -1.)
341  return;
342 
343  hDeltaEt->Fill(deltaEt);
344  hDeltaEx->Fill(recParticle->px()-genParticle->px());
345  hDeltaEy->Fill(recParticle->py()-genParticle->py());
346  hDeltaEtvsEt->Fill(et,deltaEt);
347  hDeltaEtOverEtvsEt->Fill(et,deltaEt/et);
348  hDeltaEtvsEta->Fill(eta,deltaEt);
349  hDeltaEtOverEtvsEta->Fill(eta,deltaEt/et);
350  hDeltaEtvsPhi->Fill(phi,deltaEt);
351  hDeltaEtOverEtvsPhi->Fill(phi,deltaEt/et);
352  hDeltaEtvsDeltaR->Fill(deltaR,deltaEt);
353  hDeltaEtOverEtvsDeltaR->Fill(deltaR,deltaEt/et);
354 
355  hDeltaEta->Fill(deltaEta);
356  hDeltaEtavsEt->Fill(et,deltaEta);
357  hDeltaEtavsEta->Fill(eta,deltaEta);
358 
359  hDeltaPhi->Fill(deltaPhi);
360  hDeltaPhivsEt->Fill(et,deltaPhi);
361  hDeltaPhivsEta->Fill(eta,deltaPhi);
362 
363  hDeltaR->Fill(deltaR);
364  hDeltaRvsEt->Fill(et,deltaR);
365  hDeltaRvsEta->Fill(eta,deltaR);
366 
368  entry.deltaEt = deltaEt;
369  entry.deltaEta = deltaEta;
370  entry.et = et;
371  entry.eta = eta;
372 
373  if (doMetPlots_)
374  {
375  const reco::MET* met1=static_cast<const reco::MET*>(genParticle);
376  const reco::MET* met2=static_cast<const reco::MET*>(recParticle);
377  if (met1!=NULL && met2!=NULL)
378  {
379  //std::cout << "FL: met1.sumEt() = " << (*met1).sumEt() << std::endl;
380  hTrueSumEt->Fill((*met1).sumEt());
381  hSumEt->Fill((*met2).sumEt());
382  hDeltaSetvsSet->Fill((*met1).sumEt(),(*met2).sumEt()-(*met1).sumEt());
383  hDeltaMexvsSet->Fill((*met1).sumEt(),recParticle->px()-genParticle->px());
384  hDeltaMexvsSet->Fill((*met1).sumEt(),recParticle->py()-genParticle->py());
385  if ((*met1).sumEt()>0.01) hDeltaSetOverSetvsSet->Fill((*met1).sumEt(),((*met2).sumEt()-(*met1).sumEt())/(*met1).sumEt());
386  hRecSetvsTrueSet->Fill((*met1).sumEt(),(*met2).sumEt());
387  hRecSetOverTrueSetvsTrueSet->Fill((*met1).sumEt(),(*met2).sumEt()/((*met1).sumEt()));
388  hTrueMexvsTrueSet->Fill((*met1).sumEt(),(*met1).px());
389  hTrueMexvsTrueSet->Fill((*met1).sumEt(),(*met1).py());
390  }
391  else
392  {
393  std::cout << "Warning : reco::MET* == NULL" << std::endl;
394  }
395  }
396 
397  // tree_->Fill(entry);
398 
399 }
400 
402 
403  if (Filename.size() != 0 && file_)
404  file_->Write(Filename.c_str());
405 
406 }
407 
409 {
410  file_=file;
411 }
void fillHistos(const reco::Candidate *genParticle, const reco::Candidate *recParticle, double deltaR_cut, bool plotAgainstReco)
virtual ~GenericBenchmark()
Definition: deltaR.h:79
virtual double et() const =0
transverse energy
virtual float eta() const =0
momentum pseudorapidity
void cd(void)
go to top directory (ie. root)
Definition: DQMStore.cc:411
#define BOOK1D(name, title, nbinsx, lowx, highx)
void setfile(TFile *file)
virtual float phi() const =0
momentum azimuthal angle
#define NULL
Definition: scimark2.h:8
double maxEta
T eta() const
bool accepted(const reco::Candidate *particle, double ptCut, double minEtaCut, double maxEtaCut) const
#define ETA
virtual float pt() const =0
transverse momentum
tuple DQM
Definition: DQM_cfg.py:6
Definition: MET.h:32
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
std::pair< std::string, MonitorElement * > entry
Definition: ME_MAP.h:8
virtual double py() const =0
y coordinate of momentum vector
virtual double px() const =0
x coordinate of momentum vector
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
bool plotAgainstReco
#define M_PI
Definition: BFit3D.cc:3
#define PHI
void write(std::string Filename)
tuple cout
Definition: gather_cfg.py:121
#define ET
#define debug
Definition: MEtoEDMFormat.h:34
#define SETAXES(name, xtitle, ytitle)
void setup(DQMStore *DQM=NULL, bool PlotAgainstReco_=true, float minDeltaEt=-100., float maxDeltaEt=50., float minDeltaPhi=-0.5, float maxDeltaPhi=0.5, bool doMetPlots=false)
const std::string & pwd(void) const
Definition: DQMStore.cc:406
Definition: DDAxes.h:10