CMS 3D CMS Logo
PFMETBenchmark.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/Benchmark/interface/PFMETBenchmark.h"
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 
30 PFMETBenchmark::~PFMETBenchmark() {
31  if (file_)
32  file_->Close();
33 }
34 
35 void PFMETBenchmark::write() {
36  // Store the DAQ Histograms
37  if (!outputFile_.empty()) {
38  if (dbe_)
39  dbe_->save(outputFile_);
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 
50 void PFMETBenchmark::setup(
51  string Filename, bool debug, bool plotAgainstReco, string benchmarkLabel_, DQMStore* dbe_store) {
52  debug_ = debug;
53  plotAgainstReco_ = plotAgainstReco;
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_) {
73  dbe_->setCurrentFolder(path);
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) {
187 void PFMETBenchmark::process(const reco::PFMETCollection& pfMets,
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)
210  hSETvsDeltaMET->Fill(rec_set, rec_met - true_met);
211  else
212  hSETvsDeltaMET->Fill(rec_set, rec_mex - true_mex);
213  hSETvsDeltaSET->Fill(rec_set, rec_set - true_set);
214  if (true_met > 5.0)
215  profileMETvsMETresp->Fill(true_met, (rec_met - true_met) / true_met);
216  profileSETvsSETresp->Fill(true_set, (rec_set - true_set) / true_set);
217 
218  hDeltaPFMETvstrueMET->Fill(true_met, rec_met - true_met);
219  hDeltaPFPhivstrueMET->Fill(true_met, rec_phi - true_phi);
220  hPFMETvstrueMET->Fill(true_met, rec_met);
221  hPFMET->Fill(rec_met);
222 
223  hDeltaPFMEXvstrueSET->Fill(true_set, rec_mex - true_mex);
224  hDeltaPFMEXvstrueSET->Fill(true_set, rec_mey - true_mey);
225 
226  // Calo
227  hDeltaCaloMET->Fill(calo_met - true_met);
228  hCaloMEX->Fill(calo_mex);
229  hCaloMEX->Fill(calo_mey);
230  hDeltaCaloMEX->Fill(calo_mex - true_mex);
231  hDeltaCaloMEX->Fill(calo_mey - true_mey);
232  hDeltaCaloPhi->Fill(calo_phi - true_phi);
233  hDeltaCaloSET->Fill(calo_set - true_set);
234  if (true_met > 5.0)
235  hCaloSETvsDeltaCaloMET->Fill(calo_set, calo_met - true_met);
236  else
237  hCaloSETvsDeltaCaloMET->Fill(calo_set, calo_mex - true_mex);
238  hCaloSETvsDeltaCaloSET->Fill(calo_set, calo_set - true_set);
239  if (true_met > 5.0)
240  profileCaloMETvsCaloMETresp->Fill(true_met, (calo_met - true_met) / true_met);
241  profileCaloSETvsCaloSETresp->Fill(true_set, (calo_set - true_set) / true_set);
242 
243  hDeltaCaloMETvstrueMET->Fill(true_met, calo_met - true_met);
244  hDeltaCaloPhivstrueMET->Fill(true_met, calo_phi - true_phi);
245  hCaloMETvstrueMET->Fill(true_met, calo_met);
246  hCaloMET->Fill(calo_met);
247 
248  hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mex - true_mex);
249  hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mey - true_mey);
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)
261  hTCSETvsDeltaTCMET->Fill(tc_set, tc_met - true_met);
262  else
263  hTCSETvsDeltaTCMET->Fill(tc_set, tc_mex - true_mex);
264  hTCSETvsDeltaTCSET->Fill(tc_set, tc_set - true_set);
265  if (true_met > 5.0)
266  profileTCMETvsTCMETresp->Fill(true_met, (tc_met - true_met) / true_met);
267  profileTCSETvsTCSETresp->Fill(true_set, (tc_set - true_set) / true_set);
268 
269  hDeltaTCMETvstrueMET->Fill(true_met, tc_met - true_met);
270  hDeltaTCPhivstrueMET->Fill(true_met, tc_phi - true_phi);
271  hDeltaTCMEXvstrueSET->Fill(true_set, tc_mex - true_mex);
272  hDeltaTCMEXvstrueSET->Fill(true_set, tc_mey - true_mey);
273  hTCMETvstrueMET->Fill(true_met, tc_met);
274 }
275 
276 void PFMETBenchmark::calculateQuantities(const reco::PFMETCollection& pfMets,
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 
331 void PFMETBenchmark::calculateQuantities(const reco::PFMETCollection& pfMets,
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 
427 void PFMETBenchmark::analyse() {
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 
619 double PFMETBenchmark::mpi_pi(double angle) {
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 }
PFMETBenchmark::profileTCMETvsTCMETresp
TProfile * profileTCMETvsTCMETresp
Definition: PFMETBenchmark.h:130
PFMETBenchmark::debug_
bool debug_
Definition: PFMETBenchmark.h:174
PFMETBenchmark::rec_mey
double rec_mey
Definition: PFMETBenchmark.h:158
reco::PFMETCollection
std::vector< reco::PFMET > PFMETCollection
collection of PFMET objects
Definition: PFMETCollection.h:20
PFMETBenchmark::hDeltaMET
TH1F * hDeltaMET
Definition: PFMETBenchmark.h:100
PFMETBenchmark::hTCMETvstrueMET
TH2F * hTCMETvstrueMET
Definition: PFMETBenchmark.h:140
mps_fire.i
i
Definition: mps_fire.py:428
PFMETBenchmark::hSETvsDeltaMET
TH2F * hSETvsDeltaMET
Definition: PFMETBenchmark.h:97
PFMETBenchmark::hDeltaMEX
TH1F * hDeltaMEX
Definition: PFMETBenchmark.h:99
PFMETBenchmark::calculateQuantities
void calculateQuantities(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
Definition: PFMETBenchmark.cc:276
PFMETBenchmark::hCaloMEX
TH1F * hCaloMEX
Definition: PFMETBenchmark.h:123
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
mps_update.status
status
Definition: mps_update.py:69
PFMETBenchmark::hDeltaCaloMET
TH1F * hDeltaCaloMET
Definition: PFMETBenchmark.h:125
PFMETBenchmark::hDeltaSET
TH1F * hDeltaSET
Definition: PFMETBenchmark.h:102
PFMETBenchmark::tc_set
double tc_set
Definition: PFMETBenchmark.h:170
PFMETBenchmark::hDeltaCaloPhi
TH1F * hDeltaCaloPhi
Definition: PFMETBenchmark.h:126
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PFMETBenchmark::dbe_
DQMStore * dbe_
Definition: PFMETBenchmark.h:176
BOOK2D
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
Definition: PFMETBenchmark.cc:9
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:13
PFMETBenchmark::hDeltaPFPhivstrueMET
TH2F * hDeltaPFPhivstrueMET
Definition: PFMETBenchmark.h:106
PFMETBenchmark::profileCaloSETvsCaloSETresp
TProfile * profileCaloSETvsCaloSETresp
Definition: PFMETBenchmark.h:119
SETAXES
#define SETAXES(name, xtitle, ytitle)
Definition: PFMETBenchmark.cc:15
reco::MET::sumEt
double sumEt() const
Definition: MET.h:56
PFMETBenchmark::true_set
double true_set
Definition: PFMETBenchmark.h:151
PFMETBenchmark::profileTCSETvsTCSETresp
TProfile * profileTCSETvsTCSETresp
Definition: PFMETBenchmark.h:129
DeltaPhi
Definition: deltaPhi.h:49
PFMETBenchmark::hCaloMET
TH1F * hCaloMET
Definition: PFMETBenchmark.h:113
PFMETBenchmark::hTCSETvsDeltaTCMET
TH2F * hTCSETvsDeltaTCMET
Definition: PFMETBenchmark.h:132
HLT_FULL_cff.DeltaEta
DeltaEta
Definition: HLT_FULL_cff.py:454
PFMETBenchmark::tc_phi
double tc_phi
Definition: PFMETBenchmark.h:169
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
PFMETBenchmark.h
PFMETBenchmark::hDeltaCaloPhivstrueMET
TH2F * hDeltaCaloPhivstrueMET
Definition: PFMETBenchmark.h:107
pi2
const double pi2
Definition: Thrust.cc:4
PFMETBenchmark::hDeltaCaloSET
TH1F * hDeltaCaloSET
Definition: PFMETBenchmark.h:127
PFMETBenchmark::true_met
double true_met
Definition: PFMETBenchmark.h:152
PFMETBenchmark::rec_met
double rec_met
Definition: PFMETBenchmark.h:156
reco::CaloMET
Definition: CaloMET.h:21
reco::METCollection
std::vector< reco::MET > METCollection
collection of MET objects
Definition: METCollection.h:22
PFMETBenchmark::hDeltaPFMEXvstrueSET
TH2F * hDeltaPFMEXvstrueSET
Definition: PFMETBenchmark.h:111
PFMETBenchmark::hCaloSETvsDeltaCaloSET
TH2F * hCaloSETvsDeltaCaloSET
Definition: PFMETBenchmark.h:121
PFMETBenchmark::hDeltaTCPhivstrueMET
TH2F * hDeltaTCPhivstrueMET
Definition: PFMETBenchmark.h:117
PFMETBenchmark::profileSETvsSETresp
TProfile * profileSETvsSETresp
Definition: PFMETBenchmark.h:94
reco::MET
Definition: MET.h:41
PFMETBenchmark::file_
TFile * file_
Definition: PFMETBenchmark.h:90
debug
#define debug
Definition: HDRShower.cc:19
PFMETBenchmark::true_mey
double true_mey
Definition: PFMETBenchmark.h:154
dqm::legacy::DQMStore
Definition: DQMStore.h:727
PFMETBenchmark::process
void process(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
Definition: PFMETBenchmark.cc:187
PFMETBenchmark::hCaloSETvsDeltaCaloMET
TH2F * hCaloSETvsDeltaCaloMET
Definition: PFMETBenchmark.h:122
PFMETBenchmark::hDeltaCaloMEX
TH1F * hDeltaCaloMEX
Definition: PFMETBenchmark.h:124
reco::LeafCandidate::py
double py() const final
y coordinate of momentum vector
Definition: LeafCandidate.h:142
PFMETBenchmark::hDeltaCaloMEXvstrueSET
TH2F * hDeltaCaloMEXvstrueSET
Definition: PFMETBenchmark.h:110
PFMETBenchmark::hDeltaTCSET
TH1F * hDeltaTCSET
Definition: PFMETBenchmark.h:137
PVValHelper::eta
Definition: PVValidationHelpers.h:70
PFMETBenchmark::outputFile_
std::string outputFile_
Definition: PFMETBenchmark.h:149
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
PFMETBenchmark::hTCSETvsDeltaTCSET
TH2F * hTCSETvsDeltaTCSET
Definition: PFMETBenchmark.h:131
PFMETBenchmark::analyse
void analyse()
Definition: PFMETBenchmark.cc:427
PFMETBenchmark::mpi_pi
double mpi_pi(double angle)
Definition: PFMETBenchmark.cc:619
dqm::implementation::DQMStore::setCurrentFolder
void setCurrentFolder(std::string const &fullpath) override
Definition: DQMStore.h:569
PFMETBenchmark::calo_mex
double calo_mex
Definition: PFMETBenchmark.h:162
PFMETBenchmark::setup
void setup(std::string Filename, bool debug, bool plotAgainstReco=false, std::string benchmarkLabel_="ParticleFlow", DQMStore *dbe_store=nullptr)
Definition: PFMETBenchmark.cc:50
isolatedTracks_cfi.caloJets
caloJets
Definition: isolatedTracks_cfi.py:33
dqm::implementation::DQMStore::save
DQM_DEPRECATED void save(std::string const &filename, std::string const &path="")
Definition: DQMStore.cc:761
PFMETBenchmark::hTCMEX
TH1F * hTCMEX
Definition: PFMETBenchmark.h:133
PFMETBenchmark::calo_mey
double calo_mey
Definition: PFMETBenchmark.h:163
PFMETBenchmark::hPFMET
TH1F * hPFMET
Definition: PFMETBenchmark.h:114
reco::PFMET
Definition: PFMET.h:18
PFMETBenchmark::rec_phi
double rec_phi
Definition: PFMETBenchmark.h:159
PFMETBenchmark::hCaloMETvstrueMET
TH2F * hCaloMETvstrueMET
Definition: PFMETBenchmark.h:108
PFMETBenchmark::hDeltaTCMET
TH1F * hDeltaTCMET
Definition: PFMETBenchmark.h:135
EgammaValidation_cff.pdgId
pdgId
Definition: EgammaValidation_cff.py:118
PFMETBenchmark::hSETvsDeltaSET
TH2F * hSETvsDeltaSET
Definition: PFMETBenchmark.h:96
PFMETFilter_cfi.TrueMET
TrueMET
Definition: PFMETFilter_cfi.py:23
plotAgainstReco
bool plotAgainstReco
Definition: PFJetBenchmarkAnalyzer.cc:71
PFMETBenchmark::rec_set
double rec_set
Definition: PFMETBenchmark.h:160
PFMETBenchmark::tc_mex
double tc_mex
Definition: PFMETBenchmark.h:167
PFMETBenchmark::true_phi
double true_phi
Definition: PFMETBenchmark.h:155
PFMETBenchmark::~PFMETBenchmark
virtual ~PFMETBenchmark()
Definition: PFMETBenchmark.cc:30
std
Definition: JetResolutionObject.h:76
PFMETBenchmark::calo_set
double calo_set
Definition: PFMETBenchmark.h:165
PFMETBenchmark::calo_met
double calo_met
Definition: PFMETBenchmark.h:161
PFMETBenchmark::true_mex
double true_mex
Definition: PFMETBenchmark.h:153
PFMETBenchmark::plotAgainstReco_
bool plotAgainstReco_
Definition: PFMETBenchmark.h:175
reco::LeafCandidate::phi
double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:148
PFMETBenchmark::hDeltaPhi
TH1F * hDeltaPhi
Definition: PFMETBenchmark.h:101
PFMETBenchmark::calo_phi
double calo_phi
Definition: PFMETBenchmark.h:164
PFMETBenchmark::hTrueMET
TH1F * hTrueMET
Definition: PFMETBenchmark.h:112
PFMETBenchmark::tc_met
double tc_met
Definition: PFMETBenchmark.h:166
PFMETBenchmark::hDeltaTCPhi
TH1F * hDeltaTCPhi
Definition: PFMETBenchmark.h:136
PFMETBenchmark::profileCaloMETvsCaloMETresp
TProfile * profileCaloMETvsCaloMETresp
Definition: PFMETBenchmark.h:120
PFMETBenchmark::hDeltaTCMEXvstrueSET
TH2F * hDeltaTCMEXvstrueSET
Definition: PFMETBenchmark.h:139
angle
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
PFMETBenchmark::rec_mex
double rec_mex
Definition: PFMETBenchmark.h:157
PFMETBenchmark::hMEX
TH1F * hMEX
Definition: PFMETBenchmark.h:98
PFMETBenchmark::write
void write()
Definition: PFMETBenchmark.cc:35
PFMETBenchmark::hDeltaTCMETvstrueMET
TH2F * hDeltaTCMETvstrueMET
Definition: PFMETBenchmark.h:116
PFMETBenchmark::hDeltaCaloMETvstrueMET
TH2F * hDeltaCaloMETvstrueMET
Definition: PFMETBenchmark.h:105
PFMETBenchmark::hDeltaPFMETvstrueMET
TH2F * hDeltaPFMETvstrueMET
Definition: PFMETBenchmark.h:104
pi
const Double_t pi
Definition: trackSplitPlot.h:36
castor_dqm_sourceclient_file_cfg.path
path
Definition: castor_dqm_sourceclient_file_cfg.py:37
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::CaloMETCollection
std::vector< reco::CaloMET > CaloMETCollection
collection of CaloMET objects
Definition: CaloMETCollection.h:20
BOOK1D
#define BOOK1D(name, title, nbinsx, lowx, highx)
Definition: PFMETBenchmark.cc:4
reco::LeafCandidate::px
double px() const final
x coordinate of momentum vector
Definition: LeafCandidate.h:140
PFMETBenchmark::hPFMETvstrueMET
TH2F * hPFMETvstrueMET
Definition: PFMETBenchmark.h:109
PFMETBenchmark::profileMETvsMETresp
TProfile * profileMETvsMETresp
Definition: PFMETBenchmark.h:95
PFMETBenchmark::tc_mey
double tc_mey
Definition: PFMETBenchmark.h:168
PFMETBenchmark::hDeltaTCMEX
TH1F * hDeltaTCMEX
Definition: PFMETBenchmark.h:134
PFMETBenchmark::PFMETBenchmark
PFMETBenchmark()
Definition: PFMETBenchmark.cc:28
PFMETBenchmark::hTCMET
TH1F * hTCMET
Definition: PFMETBenchmark.h:115
benchmarkLabel_
string benchmarkLabel_
Definition: PFJetBenchmarkAnalyzer.cc:74
Generated for CMSSW Reference Manual by doxygen 1.8.16