CMS 3D CMS Logo

PFMETBenchmark.cc
Go to the documentation of this file.
2 
3 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
4 #define BOOK1D(name, title, nbinsx, lowx, highx) \
5  h##name = \
6  dbe_ ? dbe_->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)
7 
8 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
9 #define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy) \
10  h##name = dbe_ ? dbe_->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
11  : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
12 
13 // all versions OK
14 // preprocesor macro for setting axis titles
15 #define SETAXES(name, xtitle, ytitle) \
16  h##name->GetXaxis()->SetTitle(xtitle); \
17  h##name->GetYaxis()->SetTitle(ytitle)
18 
19 /*#define SET2AXES(name,xtitle,ytitle) \
20  hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle); hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
21 */
22 
23 #define PT (plotAgainstReco_) ? "reconstructed P_{T}" : "generated P_{T}"
24 
25 using namespace reco;
26 using namespace std;
27 
28 PFMETBenchmark::PFMETBenchmark() : file_(nullptr) {}
29 
31  if (file_)
32  file_->Close();
33 }
34 
36  // Store the DAQ Histograms
37  if (!outputFile_.empty()) {
38  if (dbe_)
40  // use bare Root if no DQM (FWLite applications)
41  else if (file_) {
42  file_->Write(outputFile_.c_str());
43  cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
44  file_->Close();
45  }
46  } else
47  cout << "No output file specified (" << outputFile_ << "). Results will not be saved!" << endl;
48 }
49 
51  string Filename, bool debug, bool plotAgainstReco, string benchmarkLabel_, DQMStore* dbe_store) {
52  debug_ = debug;
54  outputFile_ = Filename;
55  file_ = nullptr;
56  dbe_ = dbe_store;
57  // print parameters
58  //cout<< "PFMETBenchmark Setup parameters =============================================="<<endl;
59  cout << "Filename to write histograms " << Filename << endl;
60  cout << "PFMETBenchmark debug " << debug_ << endl;
61  cout << "plotAgainstReco " << plotAgainstReco_ << endl;
62  cout << "benchmarkLabel " << benchmarkLabel_ << endl;
63 
64  // Book histogram
65 
66  // Establish DQM Store
67  string path = "PFTask/Benchmarks/" + benchmarkLabel_ + "/";
68  if (plotAgainstReco)
69  path += "Reco";
70  else
71  path += "Gen";
72  if (dbe_) {
74  } else {
75  file_ = new TFile(outputFile_.c_str(), "recreate");
76  // TTree * tr = new TTree("PFTast");
77  // tr->Branch("Benchmarks/ParticleFlow")
78  cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. " << endl;
79  }
80 
81  // delta Pt or E quantities for Barrel
82  BOOK1D(MEX, "Particle Flow", 400, -200, 200);
83  BOOK1D(DeltaMEX, "Particle Flow", 400, -200, 200);
84  BOOK1D(DeltaMET, "Particle Flow", 400, -200, 200);
85  BOOK1D(DeltaPhi, "Particle Flow", 1000, -3.2, 3.2);
86  BOOK1D(DeltaSET, "Particle Flow", 400, -200, 200);
87  BOOK2D(SETvsDeltaMET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
88  BOOK2D(SETvsDeltaSET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
89  profileSETvsSETresp = new TProfile("#DeltaPFSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
90  profileMETvsMETresp = new TProfile("#DeltaPFMET / trueMET vs trueMET", "", 50, 0.0, 200.0, -1.0, 1.0);
91 
92  BOOK1D(CaloMEX, "Calorimeter", 400, -200, 200);
93  BOOK1D(DeltaCaloMEX, "Calorimeter", 400, -200, 200);
94  BOOK1D(DeltaCaloMET, "Calorimeter", 400, -200, 200);
95  BOOK1D(DeltaCaloPhi, "Calorimeter", 1000, -3.2, 3.2);
96  BOOK1D(DeltaCaloSET, "Calorimeter", 400, -200, 200);
97  BOOK2D(CaloSETvsDeltaCaloMET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
98  BOOK2D(CaloSETvsDeltaCaloSET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
99  profileCaloSETvsCaloSETresp = new TProfile("#DeltaCaloSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
100  profileCaloMETvsCaloMETresp = new TProfile("#DeltaCaloMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
101 
102  BOOK2D(DeltaPFMETvstrueMET, "Particle Flow", 500, 0.0, 1000.0, 400, -200.0, 200.0);
103  BOOK2D(DeltaCaloMETvstrueMET, "Calorimeter", 500, 0.0, 1000.0, 400, -200.0, 200.0);
104  BOOK2D(DeltaPFPhivstrueMET, "Particle Flow", 500, 0.0, 1000.0, 400, -3.2, 3.2);
105  BOOK2D(DeltaCaloPhivstrueMET, "Calorimeter", 500, 0.0, 1000.0, 400, -3.2, 3.2);
106  BOOK2D(CaloMETvstrueMET, "Calorimeter", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
107  BOOK2D(PFMETvstrueMET, "Particle Flow", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
108 
109  BOOK2D(DeltaCaloMEXvstrueSET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
110  BOOK2D(DeltaPFMEXvstrueSET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
111 
112  BOOK1D(TrueMET, "MC truth", 500, 0.0, 1000.0);
113  BOOK1D(CaloMET, "Calorimeter", 500, 0.0, 1000.0);
114  BOOK1D(PFMET, "Particle Flow", 500, 0.0, 1000.0);
115 
116  BOOK1D(TCMEX, "Track Corrected", 400, -200, 200);
117  BOOK1D(DeltaTCMEX, "Track Corrected", 400, -200, 200);
118  BOOK1D(DeltaTCMET, "Track Corrected", 400, -200, 200);
119  BOOK1D(DeltaTCPhi, "Track Corrected", 1000, -3.2, 3.2);
120  BOOK1D(DeltaTCSET, "Track Corrected", 400, -200, 200);
121  BOOK2D(TCSETvsDeltaTCMET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
122  BOOK2D(TCSETvsDeltaTCSET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
123  profileTCSETvsTCSETresp = new TProfile("#DeltaTCSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
124  profileTCMETvsTCMETresp = new TProfile("#DeltaTCMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
125 
126  BOOK1D(TCMET, "Track Corrected", 500, 0, 1000);
127  BOOK2D(DeltaTCMETvstrueMET, "Track Corrected", 500, 0.0, 1000.0, 400, -200.0, 200.0);
128  BOOK2D(DeltaTCPhivstrueMET, "Track Corrected", 500, 0.0, 1000.0, 400, -3.2, 3.2);
129 
130  BOOK2D(DeltaTCMEXvstrueSET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
131  BOOK2D(TCMETvstrueMET, "Track Corrected", 200, 0.0, 1000.0, 500, 0.0, 1000.0);
132 
133  //BOOK1D(meanPF, "Mean PFMEX", 100, 0.0, 1600.0);
134  //BOOK1D(meanCalo, "Mean CaloMEX", 100, 0.0, 1600.0);
135  //BOOK1D(sigmaPF, "#sigma(PFMEX)", 100, 0.0, 1600.0);
136  //BOOK1D(sigmaCalo, "#sigma(CaloMEX)", 100, 0.0, 1600.0);
137  //BOOK1D(rmsPF, "RMS(PFMEX)", 100, 0.0, 1600.0);
138  //BOOK1D(rmsCalo, "RMS(CaloMEX)", 100, 0.0, 1600.0);
139 
140  // Set Axis Titles
141  // delta Pt or E quantities for Barrel and Endcap
142  SETAXES(MEX, "MEX", "Events");
143  SETAXES(DeltaMEX, "#DeltaMEX", "Events");
144  SETAXES(DeltaMET, "#DeltaMET", "Events");
145  SETAXES(DeltaPhi, "#Delta#phi", "Events");
146  SETAXES(DeltaSET, "#DeltaSET", "Events");
147  SETAXES(SETvsDeltaMET, "SET", "#DeltaMET");
148  SETAXES(SETvsDeltaSET, "SET", "#DeltaSET");
149 
150  SETAXES(DeltaPFMETvstrueMET, "trueMET", "#DeltaMET");
151  SETAXES(DeltaCaloMETvstrueMET, "trueMET", "#DeltaCaloMET");
152  SETAXES(DeltaPFPhivstrueMET, "trueMET", "#DeltaPFPhi");
153  SETAXES(DeltaCaloPhivstrueMET, "trueMET", "#DeltaCaloPhi");
154  SETAXES(CaloMETvstrueMET, "trueMET", "CaloMET");
155  SETAXES(PFMETvstrueMET, "trueMET", "PFMET");
156  SETAXES(DeltaCaloMEXvstrueSET, "trueSET", "#DeltaCaloMEX");
157  SETAXES(DeltaPFMEXvstrueSET, "trueSET", "#DeltaPFMEX");
158  SETAXES(TrueMET, "trueMET", "Events");
159  SETAXES(CaloMET, "CaloMET", "Events");
160  SETAXES(PFMET, "PFMET", "Events");
161  SETAXES(TCMET, "TCMET", "Events");
162 
163  SETAXES(CaloMEX, "MEX", "Events");
164  SETAXES(DeltaCaloMEX, "#DeltaMEX", "Events");
165  SETAXES(DeltaCaloMET, "#DeltaMET", "Events");
166  SETAXES(DeltaCaloPhi, "#Delta#phi", "Events");
167  SETAXES(DeltaCaloSET, "#DeltaSET", "Events");
168  SETAXES(CaloSETvsDeltaCaloMET, "SET", "#DeltaMET");
169  SETAXES(CaloSETvsDeltaCaloSET, "SET", "#DeltaSET");
170 
171  SETAXES(TCMEX, "MEX", "Events");
172  SETAXES(DeltaTCMEX, "#DeltaMEX", "Events");
173  SETAXES(DeltaTCMET, "#DeltaMET", "Events");
174  SETAXES(DeltaTCPhi, "#Delta#phi", "Events");
175  SETAXES(DeltaTCSET, "#DeltaSET", "Events");
176  SETAXES(TCSETvsDeltaTCMET, "SET", "#DeltaMET");
177  SETAXES(TCSETvsDeltaTCSET, "SET", "#DeltaSET");
178 
179  SETAXES(DeltaTCMETvstrueMET, "trueMET", "#DeltaTCMET");
180  SETAXES(DeltaTCPhivstrueMET, "trueMET", "#DeltaTCPhi");
181 
182  SETAXES(DeltaTCMEXvstrueSET, "trueSET", "#DeltaTCMEX");
183  SETAXES(TCMETvstrueMET, "trueMET", "TCMET");
184 }
185 
186 //void PFMETBenchmark::process(const reco::PFMETCollection& pfMets, const reco::GenMETCollection& genMets) {
188  const reco::GenParticleCollection& genParticleList,
189  const reco::CaloMETCollection& caloMets,
190  const reco::METCollection& tcMets) {
191  calculateQuantities(pfMets, genParticleList, caloMets, tcMets);
192  if (debug_) {
193  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
194  cout << " =========GenMET " << true_met << ", " << true_phi << endl;
195  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
196  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
197  }
198  // fill histograms
199  hTrueMET->Fill(true_met);
200  // delta Pt or E quantities
201  // PF
202  hDeltaMET->Fill(rec_met - true_met);
203  hMEX->Fill(rec_mex);
204  hMEX->Fill(rec_mey);
205  hDeltaMEX->Fill(rec_mex - true_mex);
206  hDeltaMEX->Fill(rec_mey - true_mey);
207  hDeltaPhi->Fill(rec_phi - true_phi);
208  hDeltaSET->Fill(rec_set - true_set);
209  if (true_met > 5.0)
211  else
214  if (true_met > 5.0)
217 
221  hPFMET->Fill(rec_met);
222 
225 
226  // Calo
228  hCaloMEX->Fill(calo_mex);
229  hCaloMEX->Fill(calo_mey);
234  if (true_met > 5.0)
236  else
239  if (true_met > 5.0)
242 
246  hCaloMET->Fill(calo_met);
247 
250 
251  // TC
252  hDeltaTCMET->Fill(tc_met - true_met);
253  hTCMET->Fill(tc_met);
254  hTCMEX->Fill(tc_mex);
255  hTCMEX->Fill(tc_mey);
256  hDeltaTCMEX->Fill(tc_mex - true_mex);
257  hDeltaTCMEX->Fill(tc_mey - true_mey);
258  hDeltaTCPhi->Fill(tc_phi - true_phi);
259  hDeltaTCSET->Fill(tc_set - true_set);
260  if (true_met > 5.0)
262  else
265  if (true_met > 5.0)
268 
274 }
275 
277  const reco::GenParticleCollection& genParticleList,
278  const reco::CaloMETCollection& caloMets,
279  const reco::METCollection& tcMets) {
280  const reco::PFMET& pfm = pfMets[0];
281  const reco::CaloMET& cm = caloMets[0];
282  const reco::MET& tcm = tcMets[0];
283 
284  double trueMEY = 0.0;
285  double trueMEX = 0.0;
286  ;
287  true_set = 0.0;
288  ;
289 
290  // for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
291  for (unsigned i = 0; i < genParticleList.size(); i++) {
292  if (genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0) {
293  if (std::abs(genParticleList[i].pdgId()) == 12 || std::abs(genParticleList[i].pdgId()) == 14 ||
294  std::abs(genParticleList[i].pdgId()) == 16 || std::abs(genParticleList[i].pdgId()) < 7 ||
295  std::abs(genParticleList[i].pdgId()) == 21) {
296  trueMEX += genParticleList[i].px();
297  trueMEY += genParticleList[i].py();
298  } else {
299  true_set += genParticleList[i].pt();
300  }
301  }
302  }
303  true_mex = -trueMEX;
304  true_mey = -trueMEY;
305  true_met = sqrt(trueMEX * trueMEX + trueMEY * trueMEY);
306  true_phi = atan2(trueMEY, trueMEX);
307  rec_met = pfm.pt();
308  rec_mex = pfm.px();
309  rec_mex = pfm.py();
310  rec_phi = pfm.phi();
311  rec_set = pfm.sumEt();
312  calo_met = cm.pt();
313  calo_mex = cm.px();
314  calo_mey = cm.py();
315  calo_phi = cm.phi();
316  calo_set = cm.sumEt();
317  tc_met = tcm.pt();
318  tc_mex = tcm.px();
319  tc_mey = tcm.py();
320  tc_phi = tcm.phi();
321  tc_set = tcm.sumEt();
322 
323  if (debug_) {
324  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
325  cout << " =========trueMET " << true_met << ", " << true_phi << endl;
326  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
327  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
328  }
329 }
330 
332  const reco::GenParticleCollection& genParticleList,
333  const reco::CaloMETCollection& caloMets,
334  const reco::METCollection& tcMets,
335  const std::vector<reco::CaloJet>& caloJets,
336  const std::vector<reco::CaloJet>& corr_caloJets) {
337  const reco::PFMET& pfm = pfMets[0];
338  const reco::CaloMET& cm = caloMets[0];
339  const reco::MET& tcm = tcMets[0];
340 
341  double trueMEY = 0.0;
342  double trueMEX = 0.0;
343  ;
344  true_set = 0.0;
345  ;
346 
347  // for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
348  for (unsigned i = 0; i < genParticleList.size(); i++) {
349  if (genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0) {
350  if (std::abs(genParticleList[i].pdgId()) == 12 || std::abs(genParticleList[i].pdgId()) == 14 ||
351  std::abs(genParticleList[i].pdgId()) == 16 || std::abs(genParticleList[i].pdgId()) < 7 ||
352  std::abs(genParticleList[i].pdgId()) == 21) {
353  trueMEX += genParticleList[i].px();
354  trueMEY += genParticleList[i].py();
355  } else {
356  true_set += genParticleList[i].pt();
357  }
358  }
359  }
360  true_mex = -trueMEX;
361  true_mey = -trueMEY;
362  true_met = sqrt(trueMEX * trueMEX + trueMEY * trueMEY);
363  true_phi = atan2(trueMEY, trueMEX);
364  rec_met = pfm.pt();
365  rec_mex = pfm.px();
366  rec_mex = pfm.py();
367  rec_phi = pfm.phi();
368  rec_set = pfm.sumEt();
369 
370  // propagation of the JEC to the caloMET:
371  double caloJetCorPX = 0.0;
372  double caloJetCorPY = 0.0;
373 
374  for (unsigned int caloJetc = 0; caloJetc < caloJets.size(); ++caloJetc) {
375  //std::cout << "caloJets[" << caloJetc << "].pt() = " << caloJets[caloJetc].pt() << std::endl;
376  //std::cout << "caloJets[" << caloJetc << "].phi() = " << caloJets[caloJetc].phi() << std::endl;
377  //std::cout << "caloJets[" << caloJetc << "].eta() = " << caloJets[caloJetc].eta() << std::endl;
378  //}
379  for (unsigned int corr_caloJetc = 0; corr_caloJetc < corr_caloJets.size(); ++corr_caloJetc) {
380  //std::cout << "corr_caloJets[" << corr_caloJetc << "].pt() = " << corr_caloJets[corr_caloJetc].pt() << std::endl;
381  //std::cout << "corr_caloJets[" << corr_caloJetc << "].phi() = " << corr_caloJets[corr_caloJetc].phi() << std::endl;
382  //std::cout << "corr_caloJets[" << corr_caloJetc << "].eta() = " << corr_caloJets[corr_caloJetc].eta() << std::endl;
383  //}
384  Float_t DeltaPhi = corr_caloJets[corr_caloJetc].phi() - caloJets[caloJetc].phi();
385  Float_t DeltaEta = corr_caloJets[corr_caloJetc].eta() - caloJets[caloJetc].eta();
386  Float_t DeltaR2 = DeltaPhi * DeltaPhi + DeltaEta * DeltaEta;
387  if (DeltaR2 < 0.0001 && caloJets[caloJetc].pt() > 20.0) {
388  caloJetCorPX += (corr_caloJets[corr_caloJetc].px() - caloJets[caloJetc].px());
389  caloJetCorPY += (corr_caloJets[corr_caloJetc].py() - caloJets[caloJetc].py());
390  }
391  }
392  }
393  double corr_calomet =
394  sqrt((cm.px() - caloJetCorPX) * (cm.px() - caloJetCorPX) + (cm.py() - caloJetCorPY) * (cm.py() - caloJetCorPY));
395  calo_met = corr_calomet;
396  calo_mex = cm.px() - caloJetCorPX;
397  calo_mey = cm.py() - caloJetCorPY;
398  calo_phi = atan2((cm.py() - caloJetCorPY), (cm.px() - caloJetCorPX));
399  //calo_met = cm.pt();
400  //calo_mex = cm.px();
401  //calo_mey = cm.py();
402  //calo_phi = cm.phi();
403 
404  calo_set = cm.sumEt();
405  tc_met = tcm.pt();
406  tc_mex = tcm.px();
407  tc_mey = tcm.py();
408  tc_phi = tcm.phi();
409  tc_set = tcm.sumEt();
410 
411  if (debug_) {
412  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
413  cout << " =========trueMET " << true_met << ", " << true_phi << endl;
414  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
415  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
416  }
417 }
418 
419 //double fitf(double *x, double *par)
420 //{
421 // double fitval = sqrt( par[0]*par[0] +
422 // par[1]*par[1]*(x[0]-par[3]) +
423 // par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]) );
424 // return fitval;
425 //}
426 
428  //Define fit functions and histograms
429  //TF1* func1 = new TF1("fit1", fitf, 0, 40, 4);
430  //TF1* func2 = new TF1("fit2", fitf, 0, 40, 4);
431  //TF1* func3 = new TF1("fit3", fitf, 0, 40, 4);
432  //TF1* func4 = new TF1("fit4", fitf, 0, 40, 4);
433 
434  //fit gaussian to Delta MET corresponding to different slices in MET, store fit values (mean,sigma) in histos
439 
446 
447  // Make the MET resolution versus SET plot
448  /*
449  TCanvas* canvas_MetResVsRecoSet = new TCanvas("MetResVsRecoSet", "MET Sigma vs Reco SET", 500,500);
450  hsigmaPF->SetStats(0);
451  func1->SetLineColor(1);
452  func1->SetParNames("Noise", "Stochastic", "Constant", "Offset");
453  func1->SetParameters(10.0, 0.8, 0.1, 100.0);
454  hsigmaPF->Fit("fit1", "", "", 100.0, 900.0);
455  func2->SetLineColor(2);
456  func2->SetParNames("Noise", "Stochastic", "Constant", "Offset");
457  func2->SetParameters(10.0, 0.8, 0.1, 100.0);
458  hsigmaCalo->Fit("fit2", "", "", 100.0, 900.0);
459  func3->SetLineColor(4);
460  func3->SetParNames("Noise", "Stochastic", "Constant", "Offset");
461  func3->SetParameters(10.0, 0.8, 0.1, 100.0);
462  hrmsPF->Fit("fit3", "", "", 100.0, 900.0);
463  func4->SetLineColor(6);
464  func4->SetParNames("Noise", "Stochastic", "Constant", "Offset");
465  func4->SetParameters(10.0, 0.8, 0.1, 100.0);
466  hrmsCalo->Fit("fit4", "", "", 100.0, 900.0);
467  (hsigmaPF->GetYaxis())->SetRangeUser( 0.0, 50.0);
468  hsigmaPF->SetLineWidth(2);
469  hsigmaPF->SetLineColor(1);
470  hsigmaPF->Draw();
471  hsigmaCalo->SetLineWidth(2);
472  hsigmaCalo->SetLineColor(2);
473  hsigmaCalo->Draw("SAME");
474  hrmsPF->SetLineWidth(2);
475  hrmsPF->SetLineColor(4);
476  hrmsPF->Draw("SAME");
477  hrmsCalo->SetLineWidth(2);
478  hrmsCalo->SetLineColor(6);
479  hrmsCalo->Draw("SAME");
480  */
481 
482  // Make the SET response versus SET plot
483  /*
484  TCanvas* canvas_SetRespVsTrueSet = new TCanvas("SetRespVsTrueSet", "SET Response vs True SET", 500,500);
485  profileSETvsSETresp->SetStats(0);
486  profileSETvsSETresp->SetStats(0);
487  (profileSETvsSETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
488  profileSETvsSETresp->SetLineWidth(2);
489  profileSETvsSETresp->SetLineColor(4);
490  profileSETvsSETresp->Draw();
491  profileCaloSETvsCaloSETresp->SetLineWidth(2);
492  profileCaloSETvsCaloSETresp->SetLineColor(2);
493  profileCaloSETvsCaloSETresp->Draw("SAME");
494  */
495 
496  // Make the MET response versus MET plot
497  /*
498  TCanvas* canvas_MetRespVsTrueMet = new TCanvas("MetRespVsTrueMet", "MET Response vs True MET", 500,500);
499  profileMETvsMETresp->SetStats(0);
500  profileMETvsMETresp->SetStats(0);
501  (profileMETvsMETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
502  profileMETvsMETresp->SetLineWidth(2);
503  profileMETvsMETresp->SetLineColor(4);
504  profileMETvsMETresp->Draw();
505  profileCaloMETvsCaloMETresp->SetLineWidth(2);
506  profileCaloMETvsCaloMETresp->SetLineColor(2);
507  profileCaloMETvsCaloMETresp->Draw("SAME");
508  */
509 
510  //print the resulting plots to file
511  /*
512  canvas_MetResVsRecoSet->Print("MetResVsRecoSet.ps");
513  canvas_SetRespVsTrueSet->Print("SetRespVsTrueSet.ps");
514  canvas_MetRespVsTrueMet->Print("MetRespVsTrueMet.ps");
515  */
516 }
517 
518 //void PFMETBenchmark::FitSlicesInY(TH2F* h, TH1F* mean, TH1F* sigma, bool doGausFit, int type )
519 //{
520 // TAxis *fXaxis = h->GetXaxis();
521 // TAxis *fYaxis = h->GetYaxis();
522 // Int_t nbins = fXaxis->GetNbins();
523 //
524 // //std::cout << "nbins = " << nbins << std::endl;
525 //
526 // Int_t binmin = 1;
527 // Int_t binmax = nbins;
528 // TString option = "QNR";
529 // TString opt = option;
530 // opt.ToLower();
531 // Float_t ngroup = 1;
532 // ngroup = 1;
533 //
534 // //std::cout << "fYaxis->GetXmin() = " << fYaxis->GetXmin() << std::endl;
535 // //std::cout << "fYaxis->GetXmax() = " << fYaxis->GetXmax() << std::endl;
536 //
537 // //default is to fit with a gaussian
538 // TF1 *f1 = 0;
539 // if (f1 == 0)
540 // {
541 // //f1 = (TF1*)gROOT->GetFunction("gaus");
542 // if (f1 == 0) f1 = new TF1("gaus","gaus", fYaxis->GetXmin(), fYaxis->GetXmax());
543 // else f1->SetRange( fYaxis->GetXmin(), fYaxis->GetXmax());
544 // }
545 // Int_t npar = f1->GetNpar();
546 //
547 // //std::cout << "npar = " << npar << std::endl;
548 //
549 // if (npar <= 0) return;
550 //
551 // Double_t *parsave = new Double_t[npar];
552 // f1->GetParameters(parsave);
553 //
555 // Int_t ipar;
556 // TH1F **hlist = new TH1F*[npar];
557 // char *name = new char[2000];
558 // char *title = new char[2000];
559 // const TArrayD *bins = fXaxis->GetXbins();
560 // for( ipar=0; ipar < npar; ipar++ )
561 // {
562 // if( ipar == 1 )
563 // if( type == 1 ) sprintf(name,"meanPF");
564 // else sprintf(name,"meanCalo");
565 // else
566 // if( doGausFit )
567 // if( type == 1 ) sprintf(name,"sigmaPF");
568 // else sprintf(name,"sigmaCalo");
569 // else
570 // if( type == 1 ) sprintf(name,"rmsPF");
571 // else sprintf(name,"rmsCalo");
572 // if( type == 1 ) sprintf(title,"Particle Flow");
573 // else sprintf(title,"Calorimeter");
574 // delete gDirectory->FindObject(name);
575 // if (bins->fN == 0)
576 // hlist[ipar] = new TH1F(name,title, nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
577 // else
578 // hlist[ipar] = new TH1F(name,title, nbins,bins->fArray);
579 // hlist[ipar]->GetXaxis()->SetTitle(fXaxis->GetTitle());
580 // }
581 // sprintf(name,"test_chi2");
582 // delete gDirectory->FindObject(name);
583 // TH1F *hchi2 = new TH1F(name,"chisquare", nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
584 // hchi2->GetXaxis()->SetTitle(fXaxis->GetTitle());
585 //
586 // //Loop on all bins in X, generate a projection along Y
587 // Int_t bin;
588 // Int_t nentries;
589 // for( bin = (Int_t) binmin; bin <= (Int_t) binmax; bin += (Int_t)ngroup )
590 // {
591 // TH1F *hpy = (TH1F*) h->ProjectionY("_temp", (Int_t) bin, (Int_t) (bin + ngroup - 1), "e");
592 // if(hpy == 0) continue;
593 // nentries = Int_t( hpy->GetEntries() );
594 // if(nentries == 0 ) {delete hpy; continue;}
595 // f1->SetParameters(parsave);
596 // hpy->Fit( f1, opt.Data() );
597 // Int_t npfits = f1->GetNumberFitPoints();
598 // //cout << "bin = " << bin << "; Npfits = " << npfits << "; npar = " << npar << endl;
599 // if( npfits > npar )
600 // {
601 // Int_t biny = bin + (Int_t)ngroup/2;
602 // for( ipar=0; ipar < npar; ipar++ )
603 // {
604 // if( doGausFit ) hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), f1->GetParameter(ipar) );
605 // else hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), hpy->GetRMS() );
606 // //cout << "bin[" << bin << "]: RMS = " << hpy->GetRMS() << "; sigma = " << f1->GetParameter(ipar) << endl;
607 // hlist[ipar]->SetBinError( biny, f1->GetParError(ipar) );
608 // }
609 // hchi2->Fill( fXaxis->GetBinCenter(biny), f1->GetChisquare()/(npfits-npar) );
610 // }
611 // delete hpy;
612 // ngroup += ngroup*0.2;//0.1 //used for non-uniform binning
613 // }
614 // *mean = *hlist[1];
615 // *sigma = *hlist[2];
616 // cout << "Entries = " << hlist[0]->GetEntries() << endl;
617 //}
618 
620  const double pi = 3.14159265358979323;
621  const double pi2 = pi * 2.;
622  while (angle > pi)
623  angle -= pi2;
624  while (angle < -pi)
625  angle += pi2;
626  return angle;
627 }
TProfile * profileTCSETvsTCSETresp
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
TProfile * profileTCMETvsTCMETresp
TH2F * hCaloSETvsDeltaCaloSET
double pt() const final
transverse momentum
TH2F * hDeltaPFPhivstrueMET
TH2F * hDeltaTCPhivstrueMET
TH2F * hDeltaTCMETvstrueMET
TProfile * profileMETvsMETresp
TH2F * hDeltaCaloMETvstrueMET
virtual ~PFMETBenchmark()
TH2F * hCaloMETvstrueMET
#define BOOK1D(name, title, nbinsx, lowx, highx)
void calculateQuantities(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
TH2F * hDeltaCaloPhivstrueMET
void setCurrentFolder(std::string const &fullpath) override
Definition: DQMStore.h:646
const double pi2
Definition: Thrust.cc:4
TProfile * profileSETvsSETresp
#define SETAXES(name, xtitle, ytitle)
double sumEt() const
Definition: MET.h:56
TH2F * hCaloSETvsDeltaCaloMET
TH2F * hSETvsDeltaMET
std::string outputFile_
std::vector< reco::MET > METCollection
collection of MET objects
Definition: METCollection.h:22
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
const Double_t pi
double mpi_pi(double angle)
double px() const final
x coordinate of momentum vector
void process(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
TH2F * hSETvsDeltaSET
Definition: MET.h:41
TH2F * hDeltaPFMETvstrueMET
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TH2F * hDeltaTCMEXvstrueSET
string benchmarkLabel_
TH2F * hDeltaCaloMEXvstrueSET
double py() const final
y coordinate of momentum vector
TH2F * hDeltaPFMEXvstrueSET
TProfile * profileCaloSETvsCaloSETresp
bool plotAgainstReco
#define debug
Definition: HDRShower.cc:19
std::vector< reco::CaloMET > CaloMETCollection
collection of CaloMET objects
void setup(std::string Filename, bool debug, bool plotAgainstReco=false, std::string benchmarkLabel_="ParticleFlow", DQMStore *dbe_store=nullptr)
DQM_DEPRECATED void save(std::string const &filename, std::string const &path="")
Definition: DQMStore.cc:791
fixed size matrix
std::vector< reco::PFMET > PFMETCollection
collection of PFMET objects
TH2F * hTCSETvsDeltaTCSET
TH2F * hTCMETvstrueMET
DQMStore * dbe_
TProfile * profileCaloMETvsCaloMETresp
double phi() const final
momentum azimuthal angle
TH2F * hTCSETvsDeltaTCMET
TH2F * hPFMETvstrueMET
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11