CMS 3D CMS Logo

JetResolution_PayloadInspector.cc
Go to the documentation of this file.
2 
6 
7 // the data format of the condition to be inspected
11 
12 #include <memory>
13 #include <sstream>
14 #include <fstream>
15 #include <iostream>
16 
17 // include ROOT
18 #include "TH2F.h"
19 #include "TLegend.h"
20 #include "TCanvas.h"
21 #include "TLine.h"
22 #include "TStyle.h"
23 #include "TLatex.h"
24 #include "TPave.h"
25 #include "TPaveStats.h"
26 
27 #define MIN_ETA -5.05
28 #define MAX_ETA 5.05
29 #define NBIN_ETA 51
30 #define MIN_PT -5.
31 #define MAX_PT 3005.
32 #define NBIN_PT 301
33 
34 namespace JME {
35 
36  using namespace cond::payloadInspector;
37 
38  enum index { NORM = 0, DOWN = 1, UP = 2 };
39 
40  /*******************************************************
41  *
42  * 1d histogram of JetResolution of 1 IOV
43  *
44  *******************************************************/
45 
46  // inherit from one of the predefined plot class: Histogram1D
47 
48  class JetResolutionVsEta : public cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV> {
49  public:
51  : cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV>(
52  "Jet Resolution", "#eta", NBIN_ETA, MIN_ETA, MAX_ETA, "Resolution") {
55  }
56 
57  bool fill() override {
58  double par_Pt = 100.;
59  double par_Eta = 1.;
60  double par_Rho = 20.;
61 
62  // Default values will be used if no input parameters
64  auto ip = paramValues.find("Jet_Pt");
65  if (ip != paramValues.end()) {
66  par_Pt = std::stod(ip->second);
67  edm::LogPrint("JER_PI") << "Jet Pt: " << par_Pt;
68  }
69  ip = paramValues.find("Jet_Rho");
70  if (ip != paramValues.end()) {
71  par_Rho = std::stod(ip->second);
72  edm::LogPrint("JER_PI") << "Jet Rho: " << par_Rho;
73  }
74 
75  auto tag = PlotBase::getTag<0>();
76  for (auto const& iov : tag.iovs) {
77  std::shared_ptr<JetResolutionObject> payload = Base::fetchPayload(std::get<1>(iov));
78  if (payload.get()) {
79  if (!payload->getRecords().empty() && // No formula for SF
80  payload->getDefinition().getFormulaString().compare("") == 0)
81  return false;
82 
83  for (const auto& record : payload->getRecords()) {
84  // Check Pt & Rho
85  if (!record.getVariablesRange().empty() && payload->getDefinition().getVariableName(0) == "JetPt" &&
86  record.getVariablesRange()[0].is_inside(par_Pt)) {
87  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(1) == "Rho" &&
88  record.getBinsRange()[1].is_inside(par_Rho)) {
89  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta") {
90  reco::FormulaEvaluator f(payload->getDefinition().getFormulaString());
91 
92  for (size_t idx = 0; idx <= NBIN_ETA; idx++) {
93  par_Eta = (idx + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
94  if (record.getBinsRange()[0].is_inside(par_Eta)) {
95  std::vector<double> var = {par_Pt};
96  std::vector<double> param;
97  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
98  double par = record.getParametersValues()[i];
99  param.push_back(par);
100  }
101  float res = f.evaluate(var, param);
102  fillWithBinAndValue(idx + 1, res);
103  }
104  }
105  }
106  }
107  }
108  } // records
109  return true;
110  }
111  }
112  return false;
113  }
114  }; // class
115 
116  class JetResolutionVsPt : public cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV> {
117  public:
119  : cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV>(
120  "Jet Energy Resolution", "p_T", NBIN_PT, MIN_PT, MAX_PT, "Resolution") {
123  }
124 
125  bool fill() override {
126  double par_Pt = 100.;
127  double par_Eta = 1.;
128  double par_Rho = 20.;
129 
131  auto ip = paramValues.find("Jet_Eta");
132  if (ip != paramValues.end()) {
133  par_Eta = std::stod(ip->second);
134  edm::LogPrint("JER_PI") << "Jet Eta: " << par_Eta;
135  }
136  ip = paramValues.find("Jet_Rho");
137  if (ip != paramValues.end()) {
138  par_Rho = std::stod(ip->second);
139  edm::LogPrint("JER_PI") << "Jet Rho: " << par_Rho;
140  }
141 
142  auto tag = PlotBase::getTag<0>();
143  for (auto const& iov : tag.iovs) {
144  std::shared_ptr<JetResolutionObject> payload = Base::fetchPayload(std::get<1>(iov));
145  if (payload.get()) {
146  if (!payload->getRecords().empty() && // No formula for SF
147  payload->getDefinition().getFormulaString().compare("") == 0)
148  return false;
149 
150  for (const auto& record : payload->getRecords()) {
151  // Check Eta & Rho
152  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta" &&
153  record.getBinsRange()[0].is_inside(par_Eta)) {
154  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(1) == "Rho" &&
155  record.getBinsRange()[1].is_inside(par_Rho)) {
156  if (!record.getVariablesRange().empty() && payload->getDefinition().getVariableName(0) == "JetPt") {
157  reco::FormulaEvaluator f(payload->getDefinition().getFormulaString());
158 
159  for (size_t idx = 0; idx <= NBIN_PT; idx++) {
160  par_Pt = (idx + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
161  if (record.getVariablesRange()[0].is_inside(par_Pt)) {
162  std::vector<double> var = {par_Pt};
163  std::vector<double> param;
164  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
165  double par = record.getParametersValues()[i];
166  param.push_back(par);
167  }
168  float res = f.evaluate(var, param);
169  fillWithBinAndValue(idx + 1, res);
170  }
171  }
172  }
173  }
174  }
175  }
176  return true;
177  }
178  }
179  return false;
180  }
181  };
182 
183  class JetResolutionSummary : public cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV> {
184  public:
186  : cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV>("Jet Resolution Summary") {
190  }
191 
192  bool fill() override {
193  double par_Pt = 100.;
194  double par_Eta = 1.;
195  double par_Rho = 20.;
196 
198  auto ip = paramValues.find("Jet_Pt");
199  if (ip != paramValues.end()) {
200  par_Pt = std::stod(ip->second);
201  }
202  ip = paramValues.find("Jet_Eta");
203  if (ip != paramValues.end()) {
204  par_Eta = std::stod(ip->second);
205  }
206  ip = paramValues.find("Jet_Rho");
207  if (ip != paramValues.end()) {
208  par_Rho = std::stod(ip->second);
209  }
210 
211  TH1D* resol_eta = new TH1D("Jet Resolution vs #eta", "", NBIN_ETA, MIN_ETA, MAX_ETA);
212  TH1D* resol_pt = new TH1D("Jet Resolution vs p_T", "", NBIN_PT, MIN_PT, MAX_PT);
213  TLegend* leg_eta = new TLegend(0.26, 0.73, 0.935, 0.90);
214  TLegend* leg_pt = new TLegend(0.26, 0.73, 0.935, 0.90);
215 
216  leg_eta->SetBorderSize(0);
217  leg_eta->SetLineStyle(0);
218  leg_eta->SetFillStyle(0);
219 
220  leg_eta->SetTextFont(42);
221  leg_pt->SetBorderSize(0);
222  leg_pt->SetLineStyle(0);
223  leg_pt->SetFillStyle(0);
224 
225  auto tag = PlotBase::getTag<0>();
226  auto iov = tag.iovs.front();
227  std::shared_ptr<JetResolutionObject> payload = fetchPayload(std::get<1>(iov));
228  unsigned int run = std::get<0>(iov);
229  std::string tagname = tag.name;
230  std::stringstream ss_tagname(tag.name);
231  std::string stmp;
232 
233  std::string tag_ver;
234  std::string tag_res;
235  std::string tag_jet;
236 
237  getline(ss_tagname, stmp, '_'); // drop first
238  getline(ss_tagname, stmp, '_'); // year
239  tag_ver = stmp;
240  getline(ss_tagname, stmp, '_'); // ver
241  tag_ver += '_' + stmp;
242  getline(ss_tagname, stmp, '_'); // cmssw
243  tag_ver += '_' + stmp;
244  getline(ss_tagname, stmp, '_'); // data/mc
245  tag_ver += '_' + stmp;
246  getline(ss_tagname, stmp, '_'); // bin
247  tag_res = stmp;
248  getline(ss_tagname, stmp, '_'); // jet algorithm
249  tag_jet = stmp;
250 
251  if (payload.get()) {
252  if (!payload->getRecords().empty() && // No formula for SF
253  payload->getDefinition().getFormulaString().compare("") == 0)
254  return false;
255 
256  for (const auto& record : payload->getRecords()) {
257  // Check Pt & Rho
258  if (!record.getVariablesRange().empty() && payload->getDefinition().getVariableName(0) == "JetPt" &&
259  record.getVariablesRange()[0].is_inside(par_Pt)) {
260  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(1) == "Rho" &&
261  record.getBinsRange()[1].is_inside(par_Rho)) {
262  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta") {
263  reco::FormulaEvaluator f(payload->getDefinition().getFormulaString());
264 
265  for (size_t idx = 0; idx <= NBIN_ETA; idx++) {
266  double x_axis = (idx + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
267  if (record.getBinsRange()[0].is_inside(x_axis)) {
268  std::vector<double> var = {par_Pt};
269  std::vector<double> param;
270  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
271  double par = record.getParametersValues()[i];
272  param.push_back(par);
273  }
274  float res = f.evaluate(var, param);
275  resol_eta->SetBinContent(idx + 1, res);
276  }
277  }
278  }
279  }
280  }
281 
282  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta" &&
283  record.getBinsRange()[0].is_inside(par_Eta)) {
284  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(1) == "Rho" &&
285  record.getBinsRange()[1].is_inside(par_Rho)) {
286  if (!record.getVariablesRange().empty() && payload->getDefinition().getVariableName(0) == "JetPt") {
287  reco::FormulaEvaluator f(payload->getDefinition().getFormulaString());
288 
289  for (size_t idx = 0; idx <= NBIN_PT + 2; idx++) {
290  double x_axis = (idx + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
291  if (record.getVariablesRange()[0].is_inside(x_axis)) {
292  std::vector<double> var = {x_axis};
293  std::vector<double> param;
294  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
295  double par = record.getParametersValues()[i];
296  param.push_back(par);
297  }
298  float res = f.evaluate(var, param);
299  resol_pt->SetBinContent(idx + 1, res);
300  }
301  }
302  }
303  }
304  }
305  } // records
306 
307  gStyle->SetOptStat(0);
308  gStyle->SetLabelFont(42, "XYZ");
309  gStyle->SetLabelSize(0.05, "XYZ");
310  gStyle->SetFrameLineWidth(3);
311 
312  std::string title = Form("Summary Run %i", run);
313  TCanvas canvas("Jet Resolution Summary", title.c_str(), 800, 1200);
314  canvas.Divide(1, 2);
315 
316  canvas.cd(1);
317  resol_eta->SetTitle(tag_res.c_str());
318  resol_eta->SetXTitle("#eta");
319  resol_eta->SetYTitle("Resolution");
320  resol_eta->SetLineWidth(3);
321  resol_eta->SetMaximum(resol_eta->GetMaximum() * 1.25);
322  resol_eta->Draw("");
323 
324  leg_eta->AddEntry(resol_eta, (tag_ver + '_' + tag_jet).c_str(), "l");
325  leg_eta->AddEntry((TObject*)nullptr, Form("JetPt=%.2f; JetRho=%.2f", par_Pt, par_Rho), "");
326  leg_eta->Draw();
327 
328  canvas.cd(2);
329  resol_pt->SetXTitle("p_{T} [GeV]");
330  resol_pt->SetYTitle("Resolution");
331  resol_pt->SetLineWidth(3);
332  resol_pt->Draw("][");
333 
334  leg_pt->AddEntry(resol_pt, (tag_ver + '_' + tag_jet).c_str(), "l");
335  leg_pt->AddEntry((TObject*)nullptr, Form("JetEta=%.2f; JetRho=%.2f", par_Eta, par_Rho), "");
336  leg_pt->Draw();
337 
338  canvas.SaveAs(m_imageFileName.c_str());
339 
340  return true;
341  } else // no payload.get()
342  return false;
343  } // fill
344 
345  }; // class
346  class JetResolutionCompare : public cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV, 2> {
347  public:
349  : cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV, 2>("Jet Resolution Compare") {
353  }
354 
355  bool fill() override {
356  double par_Pt = 100.;
357  double par_Eta = 1.;
358  double par_Rho = 20.;
359 
361  auto ip = paramValues.find("Jet_Pt");
362  if (ip != paramValues.end()) {
363  par_Pt = std::stod(ip->second);
364  }
365  ip = paramValues.find("Jet_Eta");
366  if (ip != paramValues.end()) {
367  par_Eta = std::stod(ip->second);
368  }
369  ip = paramValues.find("Jet_Rho");
370  if (ip != paramValues.end()) {
371  par_Rho = std::stod(ip->second);
372  }
373 
374  TH1D* resol_eta1 = new TH1D("Jet Resolution vs #eta one", "", NBIN_ETA, MIN_ETA, MAX_ETA);
375  TH1D* resol_eta2 = new TH1D("Jet Resolution vs #eta two", "", NBIN_ETA, MIN_ETA, MAX_ETA);
376  TH1D* resol_pt1 = new TH1D("Jet Resolution vs p_T one", "", NBIN_PT, MIN_PT, MAX_PT);
377  TH1D* resol_pt2 = new TH1D("Jet Resolution vs p_T two", "", NBIN_PT, MIN_PT, MAX_PT);
378  TLegend* leg_eta = new TLegend(0.26, 0.73, 0.935, 0.90);
379  TLegend* leg_pt = new TLegend(0.26, 0.73, 0.935, 0.90);
380 
381  leg_eta->SetBorderSize(0);
382  leg_eta->SetLineStyle(0);
383  leg_eta->SetFillStyle(0);
384 
385  leg_eta->SetTextFont(42);
386  leg_pt->SetBorderSize(0);
387  leg_pt->SetLineStyle(0);
388  leg_pt->SetFillStyle(0);
389 
390  auto tag1 = PlotBase::getTag<0>();
391  auto iov1 = tag1.iovs.front();
392  std::shared_ptr<JetResolutionObject> payload1 = fetchPayload(std::get<1>(iov1));
393  auto tag2 = PlotBase::getTag<1>();
394  auto iov2 = tag2.iovs.front();
395  std::shared_ptr<JetResolutionObject> payload2 = fetchPayload(std::get<1>(iov2));
396 
397  std::stringstream ss_tagname1(tag1.name);
398  std::stringstream ss_tagname2(tag1.name);
399  std::string stmp;
400 
401  std::string tag_ver1;
402  std::string tag_ver2;
403 
404  getline(ss_tagname1, stmp, '_'); // drop first
405  getline(ss_tagname1, stmp); // year
406  tag_ver1 = stmp;
407 
408  getline(ss_tagname2, stmp, '_'); // drop first
409  getline(ss_tagname2, stmp); // year
410  tag_ver2 = stmp;
411 
412  if (payload1.get() && payload2.get()) {
413  if (!payload1->getRecords().empty() && // No formula for SF
414  payload1->getDefinition().getFormulaString().compare("") == 0)
415  return false;
416 
417  for (const auto& record : payload1->getRecords()) {
418  // Check Pt & Rho
419  if (!record.getVariablesRange().empty() && payload1->getDefinition().getVariableName(0) == "JetPt" &&
420  record.getVariablesRange()[0].is_inside(par_Pt)) {
421  if (record.getBinsRange().size() > 1 && payload1->getDefinition().getBinName(1) == "Rho" &&
422  record.getBinsRange()[1].is_inside(par_Rho)) {
423  if (!record.getBinsRange().empty() && payload1->getDefinition().getBinName(0) == "JetEta") {
424  reco::FormulaEvaluator f(payload1->getDefinition().getFormulaString());
425 
426  for (size_t idx = 0; idx <= NBIN_ETA; idx++) {
427  double x_axis = (idx + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
428  if (record.getBinsRange()[0].is_inside(x_axis)) {
429  std::vector<double> var = {par_Pt};
430  std::vector<double> param;
431  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
432  double par = record.getParametersValues()[i];
433  param.push_back(par);
434  }
435  float res = f.evaluate(var, param);
436  resol_eta1->SetBinContent(idx + 1, res);
437  }
438  }
439  }
440  }
441  }
442 
443  if (!record.getBinsRange().empty() && payload1->getDefinition().getBinName(0) == "JetEta" &&
444  record.getBinsRange()[0].is_inside(par_Eta)) {
445  if (record.getBinsRange().size() > 1 && payload1->getDefinition().getBinName(1) == "Rho" &&
446  record.getBinsRange()[1].is_inside(par_Rho)) {
447  if (!record.getVariablesRange().empty() && payload1->getDefinition().getVariableName(0) == "JetPt") {
448  reco::FormulaEvaluator f(payload1->getDefinition().getFormulaString());
449 
450  for (size_t idx = 0; idx <= NBIN_PT + 2; idx++) {
451  double x_axis = (idx + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
452  if (record.getVariablesRange()[0].is_inside(x_axis)) {
453  std::vector<double> var = {x_axis};
454  std::vector<double> param;
455  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
456  double par = record.getParametersValues()[i];
457  param.push_back(par);
458  }
459  float res = f.evaluate(var, param);
460  resol_pt1->SetBinContent(idx + 1, res);
461  }
462  }
463  }
464  }
465  }
466  } // records
467 
468  for (const auto& record : payload2->getRecords()) {
469  // Check Pt & Rho
470  if (!record.getVariablesRange().empty() && payload2->getDefinition().getVariableName(0) == "JetPt" &&
471  record.getVariablesRange()[0].is_inside(par_Pt)) {
472  if (record.getBinsRange().size() > 1 && payload2->getDefinition().getBinName(1) == "Rho" &&
473  record.getBinsRange()[1].is_inside(par_Rho)) {
474  if (!record.getBinsRange().empty() && payload2->getDefinition().getBinName(0) == "JetEta") {
475  reco::FormulaEvaluator f(payload2->getDefinition().getFormulaString());
476 
477  for (size_t idx = 0; idx <= NBIN_ETA; idx++) {
478  double x_axis = (idx + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
479  if (record.getBinsRange()[0].is_inside(x_axis)) {
480  std::vector<double> var = {par_Pt};
481  std::vector<double> param;
482  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
483  double par = record.getParametersValues()[i];
484  param.push_back(par);
485  }
486  float res = f.evaluate(var, param);
487  resol_eta2->SetBinContent(idx + 1, res);
488  }
489  }
490  }
491  }
492  }
493 
494  if (!record.getBinsRange().empty() && payload2->getDefinition().getBinName(0) == "JetEta" &&
495  record.getBinsRange()[0].is_inside(par_Eta)) {
496  if (record.getBinsRange().size() > 1 && payload2->getDefinition().getBinName(1) == "Rho" &&
497  record.getBinsRange()[1].is_inside(par_Rho)) {
498  if (!record.getVariablesRange().empty() && payload2->getDefinition().getVariableName(0) == "JetPt") {
499  reco::FormulaEvaluator f(payload2->getDefinition().getFormulaString());
500 
501  for (size_t idx = 0; idx <= NBIN_PT + 2; idx++) {
502  double x_axis = (idx + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
503  if (record.getVariablesRange()[0].is_inside(x_axis)) {
504  std::vector<double> var = {x_axis};
505  std::vector<double> param;
506  for (size_t i = 0; i < record.getParametersValues().size(); i++) {
507  double par = record.getParametersValues()[i];
508  param.push_back(par);
509  }
510  float res = f.evaluate(var, param);
511  resol_pt2->SetBinContent(idx + 1, res);
512  }
513  }
514  }
515  }
516  }
517  } // records
518 
519  gStyle->SetOptStat(0);
520  gStyle->SetLabelFont(42, "XYZ");
521  gStyle->SetLabelSize(0.05, "XYZ");
522  gStyle->SetFrameLineWidth(3);
523 
524  std::string title = Form("Comparison between %s and %s", tag1.name.c_str(), tag2.name.c_str());
525  TCanvas canvas("Jet Resolution Comparison", title.c_str(), 800, 1200);
526  canvas.Divide(1, 2);
527 
528  canvas.cd(1);
529  resol_eta1->SetTitle("Jet Resolution Comparison");
530  resol_eta1->SetXTitle("#eta");
531  resol_eta1->SetYTitle("Resolution");
532  resol_eta1->SetLineWidth(3);
533  resol_eta1->SetMaximum(resol_eta1->GetMaximum() * 1.25);
534  resol_eta1->Draw("][");
535 
536  resol_eta2->SetLineColor(2);
537  resol_eta2->SetLineWidth(3);
538  resol_eta2->SetLineStyle(2);
539  resol_eta2->Draw("][same");
540 
541  leg_eta->AddEntry(resol_eta1, tag_ver1.c_str(), "l");
542  leg_eta->AddEntry(resol_eta2, tag_ver2.c_str(), "l");
543  leg_eta->AddEntry((TObject*)nullptr, Form("JetPt=%.2f; JetRho=%.2f", par_Pt, par_Rho), "");
544  leg_eta->Draw();
545 
546  canvas.cd(2);
547  resol_pt1->SetXTitle("p_{T} [GeV]");
548  resol_pt1->SetYTitle("Resolution");
549  resol_pt1->SetLineWidth(3);
550  resol_pt1->Draw("][");
551 
552  resol_pt2->SetLineColor(2);
553  resol_pt2->SetLineWidth(3);
554  resol_pt2->SetLineStyle(2);
555  resol_pt2->Draw("][same");
556 
557  leg_pt->AddEntry(resol_pt1, tag_ver1.c_str(), "l");
558  leg_pt->AddEntry(resol_pt2, tag_ver2.c_str(), "l");
559  leg_pt->AddEntry((TObject*)nullptr, Form("JetEta=%.2f; JetRho=%.2f", par_Eta, par_Rho), "");
560  leg_pt->Draw();
561 
562  canvas.SaveAs(m_imageFileName.c_str());
563 
564  return true;
565  } else // no payload.get()
566  return false;
567  } // fill
568 
569  }; // class
570 
571  template <index ii>
572  class JetScaleFactorVsEta : public cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV> {
573  public:
575  : cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV>(
576  "Jet Energy Scale Factor", "#eta", NBIN_ETA, MIN_ETA, MAX_ETA, "Scale Factor") {
579  }
580 
581  bool fill() override {
582  double par_Pt = 100.;
583  double par_Eta = 1.;
584 
586  auto ip = paramValues.find("Jet_Pt");
587  if (ip != paramValues.end()) {
588  par_Pt = std::stod(ip->second);
589  edm::LogPrint("JER_PI") << "Jet Pt: " << par_Pt;
590  }
591  ip = paramValues.find("Jet_Eta");
592  if (ip != paramValues.end()) {
593  par_Eta = std::stod(ip->second);
594  edm::LogPrint("JER_PI") << "Jet Eta: " << par_Eta;
595  }
596 
597  auto tag = PlotBase::getTag<0>();
598  for (auto const& iov : tag.iovs) {
599  std::shared_ptr<JetResolutionObject> payload = Base::fetchPayload(std::get<1>(iov));
600  if (payload.get()) {
601  if (!payload->getRecords().empty() && // No formula for SF
602  payload->getDefinition().getFormulaString().compare("") != 0)
603  return false;
604 
605  for (const auto& record : payload->getRecords()) {
606  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta" &&
607  record.getParametersValues().size() == 3) { // norm, down, up
608 
609  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(1) == "JetPt" &&
610  !record.getBinsRange()[1].is_inside(par_Pt))
611  continue; // for 2-bin payload, take jetpt=500
612 
613  for (size_t it = 0; it <= NBIN_ETA; it++) {
614  par_Eta = (it + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
615  if (record.getBinsRange()[0].is_inside(par_Eta)) {
616  double sf = 0.;
617  sf = record.getParametersValues()[ii];
618  fillWithBinAndValue(it, sf);
619  }
620  }
621  }
622  } // records
623  return true;
624  } else
625  return false;
626  } // for
627  return false;
628  } // fill
629  }; // class
630 
634 
635  template <index ii>
636  class JetScaleFactorVsPt : public cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV> {
637  public:
639  : cond::payloadInspector::Histogram1D<JetResolutionObject, SINGLE_IOV>(
640  "Jet Energy Scale Factor", "p_T", NBIN_PT, MIN_PT, MAX_PT, "Scale Factor") {
643  }
644 
645  bool fill() override {
646  double par_Pt = 100.;
647  double par_Eta = 1.;
648 
650  auto ip = paramValues.find("Jet_Pt");
651  if (ip != paramValues.end()) {
652  par_Pt = std::stod(ip->second);
653  edm::LogPrint("JER_PI") << "Jet Pt: " << par_Pt;
654  }
655  ip = paramValues.find("Jet_Eta");
656  if (ip != paramValues.end()) {
657  par_Eta = std::stod(ip->second);
658  edm::LogPrint("JER_PI") << "Jet Eta: " << par_Eta;
659  }
660 
661  auto tag = PlotBase::getTag<0>();
662  for (auto const& iov : tag.iovs) {
663  std::shared_ptr<JetResolutionObject> payload = Base::fetchPayload(std::get<1>(iov));
664  if (payload.get()) {
665  if (!payload->getRecords().empty() && // No formula for SF
666  payload->getDefinition().getFormulaString().compare("") != 0)
667  return false;
668 
669  for (const auto& record : payload->getRecords()) {
670  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(0) == "JetEta" &&
671  record.getBinsRange()[0].is_inside(par_Eta) && // take jeteta=2.5
672  payload->getDefinition().getBinName(1) == "JetPt" && // 2-bin
673  record.getParametersValues().size() == 3) { // norm, down, up
674 
675  for (size_t it = 0; it <= NBIN_PT; it++) {
676  par_Pt = (it + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
677  if (record.getBinsRange()[1].is_inside(par_Pt)) {
678  double sf = 0.;
679  sf = record.getParametersValues()[ii];
680  fillWithBinAndValue(it, sf);
681  }
682  }
683  }
684  } // records
685  return true;
686  } else
687  return false;
688  } // for
689  return false;
690  } // fill
691  }; // class
692 
696 
697  class JetScaleFactorSummary : public cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV> {
698  public:
700  : cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV>("Jet ScaleFactor Summary") {
703  }
704 
705  bool fill() override {
706  double par_Pt = 100.;
707  double par_Eta = 1.;
708 
710  auto ip = paramValues.find("Jet_Pt");
711  if (ip != paramValues.end()) {
712  par_Pt = std::stod(ip->second);
713  }
714  ip = paramValues.find("Jet_Eta");
715  if (ip != paramValues.end()) {
716  par_Eta = std::stod(ip->second);
717  }
718 
719  TH1D* sf_eta_norm = new TH1D("Jet SF vs #eta NORM", "", NBIN_ETA, MIN_ETA, MAX_ETA);
720  TH1D* sf_eta_down = new TH1D("Jet SF vs #eta DOWN", "", NBIN_ETA, MIN_ETA, MAX_ETA);
721  TH1D* sf_eta_up = new TH1D("Jet SF vs #eta UP", "", NBIN_ETA, MIN_ETA, MAX_ETA);
722  TH1D* sf_pt_norm = new TH1D("Jet SF vs p_T NORM", "", NBIN_PT, MIN_PT, MAX_PT);
723  TH1D* sf_pt_down = new TH1D("Jet SF vs p_T DOWN", "", NBIN_PT, MIN_PT, MAX_PT);
724  TH1D* sf_pt_up = new TH1D("Jet SF vs p_T UP", "", NBIN_PT, MIN_PT, MAX_PT);
725 
726  TLegend* leg_eta = new TLegend(0.26, 0.73, 0.935, 0.90);
727  TLegend* leg_pt = new TLegend(0.26, 0.73, 0.935, 0.90);
728 
729  leg_eta->SetBorderSize(0);
730  leg_eta->SetLineStyle(0);
731  leg_eta->SetFillStyle(0);
732 
733  leg_eta->SetTextFont(42);
734  leg_pt->SetBorderSize(0);
735  leg_pt->SetLineStyle(0);
736  leg_pt->SetFillStyle(0);
737 
738  auto tag = PlotBase::getTag<0>();
739  auto iov = tag.iovs.front();
740  std::shared_ptr<JetResolutionObject> payload = fetchPayload(std::get<1>(iov));
741  unsigned int run = std::get<0>(iov);
742  std::string tagname = tag.name;
743  std::stringstream ss_tagname(tag.name);
744  std::string stmp;
745 
746  std::string tag_ver;
747  std::string tag_res;
748  std::string tag_jet;
749 
750  getline(ss_tagname, stmp, '_'); // drop first
751  getline(ss_tagname, stmp, '_'); // year
752  tag_ver = stmp;
753  getline(ss_tagname, stmp, '_'); // ver
754  tag_ver += '_' + stmp;
755  getline(ss_tagname, stmp, '_'); // cmssw
756  tag_ver += '_' + stmp;
757  getline(ss_tagname, stmp, '_'); // data/mc
758  tag_ver += '_' + stmp;
759  getline(ss_tagname, stmp, '_'); // bin
760  tag_res = stmp;
761  getline(ss_tagname, stmp, '_'); // jet algorithm
762  tag_jet = stmp;
763 
764  bool is_2bin = false;
765 
766  if (payload.get()) {
767  if (!payload->getRecords().empty() && // No formula for SF
768  payload->getDefinition().getFormulaString().compare("") != 0)
769  return false;
770 
771  is_2bin = false;
772  for (const auto& record : payload->getRecords()) {
773  if (record.getBinsRange().size() > 1)
774  is_2bin = true;
775 
776  if (!record.getBinsRange().empty() && payload->getDefinition().getBinName(0) == "JetEta" &&
777  record.getParametersValues().size() == 3) { // norm, down, up
778 
779  for (size_t it = 0; it <= NBIN_ETA; it++) {
780  double x_axis = (it + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
781  if (((is_2bin == false) || (is_2bin == true && record.getBinsRange()[1].is_inside(par_Pt))) &&
782  record.getBinsRange()[0].is_inside(x_axis)) {
783  sf_eta_norm->SetBinContent(it + 1, record.getParametersValues()[0]);
784  sf_eta_down->SetBinContent(it + 1, record.getParametersValues()[1]);
785  sf_eta_up->SetBinContent(it + 1, record.getParametersValues()[2]);
786  }
787  }
788  }
789 
790  if (record.getBinsRange().size() > 1 && payload->getDefinition().getBinName(0) == "JetEta" &&
791  record.getBinsRange()[0].is_inside(par_Eta) && // take jeteta=2.5
792  payload->getDefinition().getBinName(1) == "JetPt" &&
793  record.getParametersValues().size() == 3) { // norm, down, up
794 
795  is_2bin = true;
796 
797  for (size_t it = 0; it <= NBIN_PT; it++) {
798  double x_axis = (it + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
799  if (record.getBinsRange()[1].is_inside(x_axis)) {
800  sf_pt_norm->SetBinContent(it + 1, record.getParametersValues()[0]);
801  sf_pt_down->SetBinContent(it + 1, record.getParametersValues()[1]);
802  sf_pt_up->SetBinContent(it + 1, record.getParametersValues()[2]);
803  }
804  }
805  } // 2-bin
806  } // records
807 
808  gStyle->SetOptStat(0);
809  gStyle->SetLabelFont(42, "XYZ");
810  gStyle->SetLabelSize(0.05, "XYZ");
811  gStyle->SetFrameLineWidth(3);
812 
813  std::string title = Form("Summary Run %i", run);
814  TCanvas canvas("Jet ScaleFactor Summary", title.c_str(), 800, 1200);
815  canvas.Divide(1, 2);
816 
817  canvas.cd(1);
818  sf_eta_up->SetTitle("ScaleFactor vs. #eta");
819  sf_eta_up->SetXTitle("#eta");
820  sf_eta_up->SetYTitle("Scale Factor");
821  sf_eta_up->SetLineStyle(7);
822  sf_eta_up->SetLineWidth(3);
823  sf_eta_up->SetFillColorAlpha(kGray, 0.5);
824  sf_eta_up->SetMaximum(sf_eta_up->GetMaximum() * 1.25);
825  sf_eta_up->SetMinimum(0.);
826  sf_eta_up->Draw("][");
827 
828  sf_eta_down->SetLineStyle(7);
829  sf_eta_down->SetLineWidth(3);
830  sf_eta_down->SetFillColorAlpha(kWhite, 1);
831  sf_eta_down->Draw("][ same");
832 
833  sf_eta_norm->SetLineStyle(1);
834  sf_eta_norm->SetLineWidth(5);
835  sf_eta_norm->SetFillColor(0);
836  sf_eta_norm->Draw("][ same");
837  sf_eta_norm->Draw("axis same");
838 
839  leg_eta->AddEntry(sf_eta_norm, (tag_ver + '_' + tag_jet).c_str(), "l");
840  leg_eta->AddEntry((TObject*)nullptr, Form("JetPt=%.2f", par_Pt), "");
841  leg_eta->Draw();
842 
843  if (is_2bin == true) {
844  canvas.cd(2);
845  sf_pt_up->SetTitle("ScaleFactor vs. p_{T}");
846  sf_pt_up->SetXTitle("p_{T} [GeV]");
847  sf_pt_up->SetYTitle("Scale Factor");
848  sf_pt_up->SetLineStyle(7);
849  sf_pt_up->SetLineWidth(3);
850  sf_pt_up->SetFillColorAlpha(kGray, 0.5);
851  sf_pt_up->SetMaximum(sf_pt_up->GetMaximum() * 1.25);
852  sf_pt_up->SetMinimum(0.);
853  sf_pt_up->Draw("][");
854 
855  sf_pt_down->SetLineStyle(7);
856  sf_pt_down->SetLineWidth(3);
857  sf_pt_down->SetFillColorAlpha(kWhite, 1);
858  sf_pt_down->Draw("][ same");
859 
860  sf_pt_norm->SetLineStyle(1);
861  sf_pt_norm->SetLineWidth(5);
862  sf_pt_norm->SetFillColor(0);
863  sf_pt_norm->Draw("][ same");
864  sf_pt_norm->Draw("axis same");
865 
866  leg_pt->AddEntry(sf_pt_norm, (tag_ver + '_' + tag_jet).c_str(), "l");
867  leg_pt->AddEntry((TObject*)nullptr, Form("JetEta=%.2f", par_Eta), "");
868  leg_pt->Draw();
869  }
870 
871  canvas.SaveAs(m_imageFileName.c_str());
872 
873  return true;
874  } else // no payload.get()
875  return false;
876  } // fill
877 
878  }; // class
879 
880  class JetScaleFactorCompare : public cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV, 2> {
881  public:
883  : cond::payloadInspector::PlotImage<JetResolutionObject, SINGLE_IOV, 2>("Jet ScaleFactor Summary") {
886  }
887 
888  bool fill() override {
889  double par_Pt = 100.;
890  double par_Eta = 1.;
891 
893  auto ip = paramValues.find("Jet_Pt");
894  if (ip != paramValues.end()) {
895  par_Pt = std::stod(ip->second);
896  }
897  ip = paramValues.find("Jet_Eta");
898  if (ip != paramValues.end()) {
899  par_Eta = std::stod(ip->second);
900  }
901 
902  TH1D* sf_eta_one = new TH1D("Jet SF vs #eta one", "", NBIN_ETA, MIN_ETA, MAX_ETA);
903  TH1D* sf_eta_two = new TH1D("Jet SF vs #eta two", "", NBIN_ETA, MIN_ETA, MAX_ETA);
904  TH1D* sf_pt_one = new TH1D("Jet SF vs p_T one", "", NBIN_PT, MIN_PT, MAX_PT);
905  TH1D* sf_pt_two = new TH1D("Jet SF vs p_T two", "", NBIN_PT, MIN_PT, MAX_PT);
906 
907  TLegend* leg_eta = new TLegend(0.26, 0.73, 0.935, 0.90);
908  TLegend* leg_pt = new TLegend(0.26, 0.73, 0.935, 0.90);
909 
910  leg_eta->SetBorderSize(0);
911  leg_eta->SetLineStyle(0);
912  leg_eta->SetFillStyle(0);
913 
914  leg_eta->SetTextFont(42);
915  leg_pt->SetBorderSize(0);
916  leg_pt->SetLineStyle(0);
917  leg_pt->SetFillStyle(0);
918 
919  auto tag1 = PlotBase::getTag<0>();
920  auto iov1 = tag1.iovs.front();
921  std::shared_ptr<JetResolutionObject> payload1 = fetchPayload(std::get<1>(iov1));
922 
923  auto tag2 = PlotBase::getTag<1>();
924  auto iov2 = tag2.iovs.front();
925  std::shared_ptr<JetResolutionObject> payload2 = fetchPayload(std::get<1>(iov2));
926 
927  std::stringstream ss_tagname1(tag1.name);
928  std::stringstream ss_tagname2(tag2.name);
929  std::string stmp;
930 
931  std::string tag_ver1;
932  std::string tag_ver2;
933 
934  getline(ss_tagname1, stmp, '_'); // drop first
935  getline(ss_tagname1, stmp); // year
936  tag_ver1 = stmp;
937 
938  getline(ss_tagname2, stmp, '_'); // drop first
939  getline(ss_tagname2, stmp); // year
940  tag_ver2 = stmp;
941 
942  bool is_2bin = false;
943 
944  if (payload1.get() && payload2.get()) {
945  if (!payload1->getRecords().empty() && // No formula for SF
946  payload1->getDefinition().getFormulaString().compare("") != 0)
947  return false;
948 
949  is_2bin = false;
950  for (const auto& record : payload1->getRecords()) {
951  if (record.getBinsRange().size() > 1)
952  is_2bin = true;
953 
954  if (!record.getBinsRange().empty() && payload1->getDefinition().getBinName(0) == "JetEta" &&
955  record.getParametersValues().size() == 3) { // norm, down, up
956 
957  for (size_t it = 0; it <= NBIN_ETA; it++) {
958  double x_axis = (it + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
959  if (((is_2bin == false) || (is_2bin == true && record.getBinsRange()[1].is_inside(par_Pt))) &&
960  record.getBinsRange()[0].is_inside(x_axis)) {
961  sf_eta_one->SetBinContent(it + 1, record.getParametersValues()[0]);
962  }
963  }
964  }
965 
966  if (record.getBinsRange().size() > 1 && payload1->getDefinition().getBinName(0) == "JetEta" &&
967  record.getBinsRange()[0].is_inside(par_Eta) && // take jeteta=2.5
968  payload1->getDefinition().getBinName(1) == "JetPt" &&
969  record.getParametersValues().size() == 3) { // norm, down, up
970 
971  is_2bin = true;
972 
973  for (size_t it = 0; it <= NBIN_PT; it++) {
974  double x_axis = (it + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
975  if (record.getBinsRange()[1].is_inside(x_axis)) {
976  sf_pt_one->SetBinContent(it + 1, record.getParametersValues()[0]);
977  }
978  }
979  } // 2-bin
980  } // records
981 
982  is_2bin = false;
983  for (const auto& record : payload2->getRecords()) {
984  if (record.getBinsRange().size() > 1)
985  is_2bin = true;
986 
987  if (!record.getBinsRange().empty() && payload2->getDefinition().getBinName(0) == "JetEta" &&
988  record.getParametersValues().size() == 3) { // norm, down, up
989 
990  for (size_t it = 0; it <= NBIN_ETA; it++) {
991  double x_axis = (it + 0.5) * (MAX_ETA - MIN_ETA) / NBIN_ETA + MIN_ETA;
992  if (((is_2bin == false) || (is_2bin == true && record.getBinsRange()[1].is_inside(par_Pt))) &&
993  record.getBinsRange()[0].is_inside(x_axis)) {
994  sf_eta_two->SetBinContent(it + 1, record.getParametersValues()[0]);
995  }
996  }
997  }
998 
999  if (record.getBinsRange().size() > 1 && payload2->getDefinition().getBinName(0) == "JetEta" &&
1000  record.getBinsRange()[0].is_inside(par_Eta) && // take jeteta=2.5
1001  payload2->getDefinition().getBinName(1) == "JetPt" &&
1002  record.getParametersValues().size() == 3) { // norm, down, up
1003 
1004  is_2bin = true;
1005 
1006  for (size_t it = 0; it <= NBIN_PT; it++) {
1007  double x_axis = (it + 0.5) * (MAX_PT - MIN_PT) / NBIN_PT + MIN_PT;
1008  if (record.getBinsRange()[1].is_inside(x_axis)) {
1009  sf_pt_two->SetBinContent(it + 1, record.getParametersValues()[0]);
1010  }
1011  }
1012  } // 2-bin
1013  } // records
1014 
1015  gStyle->SetOptStat(0);
1016  gStyle->SetLabelFont(42, "XYZ");
1017  gStyle->SetLabelSize(0.05, "XYZ");
1018  gStyle->SetFrameLineWidth(3);
1019 
1020  std::string title = Form("Comparison");
1021  TCanvas canvas("Jet ScaleFactor", title.c_str(), 800, 1200);
1022  canvas.Divide(1, 2);
1023 
1024  canvas.cd(1);
1025  sf_eta_one->SetTitle("Jet ScaleFactor vs. #eta");
1026  sf_eta_one->SetXTitle("#eta");
1027  sf_eta_one->SetYTitle("Scale Factor");
1028  sf_eta_one->SetLineWidth(3);
1029  sf_eta_one->SetMaximum(sf_eta_one->GetMaximum() * 1.25);
1030  sf_eta_one->SetMinimum(0.);
1031  sf_eta_one->Draw("][");
1032 
1033  sf_eta_two->SetLineColor(2);
1034  sf_eta_two->SetLineWidth(3);
1035  sf_eta_two->SetLineStyle(2);
1036  sf_eta_two->Draw("][ same");
1037 
1038  leg_eta->AddEntry(sf_eta_one, tag_ver1.c_str(), "l");
1039  leg_eta->AddEntry(sf_eta_two, tag_ver2.c_str(), "l");
1040  leg_eta->AddEntry((TObject*)nullptr, Form("JetPt=%.2f", par_Pt), "");
1041  leg_eta->Draw();
1042 
1043  if (is_2bin == true) {
1044  canvas.cd(2);
1045  sf_pt_one->SetTitle("Jet ScaleFactor vs. p_{T}");
1046  sf_pt_one->SetXTitle("p_{T} [GeV]");
1047  sf_pt_one->SetYTitle("Scale Factor");
1048  sf_pt_one->SetLineWidth(3);
1049  sf_pt_one->SetMaximum(sf_pt_one->GetMaximum() * 1.25);
1050  sf_pt_one->SetMinimum(0.);
1051  sf_pt_one->Draw("][");
1052 
1053  sf_pt_two->SetLineColor(2);
1054  sf_pt_two->SetLineWidth(3);
1055  sf_pt_two->SetLineStyle(2);
1056  sf_pt_two->Draw("][ same");
1057 
1058  leg_pt->AddEntry(sf_pt_one, tag_ver1.c_str(), "l");
1059  leg_pt->AddEntry(sf_pt_two, tag_ver2.c_str(), "l");
1060  leg_pt->AddEntry((TObject*)nullptr, Form("JetEta=%.2f", par_Eta), "");
1061  leg_pt->Draw();
1062  }
1063 
1064  canvas.SaveAs(m_imageFileName.c_str());
1065 
1066  return true;
1067  } else // no payload.get()
1068  return false;
1069  } // fill
1070 
1071  }; // class
1072 
1073  // Register the classes as boost python plugin
1083 
1088  }
1089 
1090 } // namespace JME
JetScaleFactorVsEta< NORM > JetScaleFactorVsEtaNORM
JetScaleFactorVsEta< UP > JetScaleFactorVsEtaUP
JetScaleFactorVsPt< NORM > JetScaleFactorVsPtNORM
Definition: Electron.h:6
#define PAYLOAD_INSPECTOR_CLASS(CLASS_NAME)
JetScaleFactorVsPt< DOWN > JetScaleFactorVsPtDOWN
double f[11][100]
#define PAYLOAD_INSPECTOR_MODULE(PAYLOAD_TYPENAME)
Log< level::Warning, true > LogPrint
ii
Definition: cuy.py:589
const std::map< std::string, std::string > & inputParamValues() const
void addInputParam(const std::string &paramName)
Definition: plugin.cc:23
JetScaleFactorVsPt< UP > JetScaleFactorVsPtUP
def canvas(sub, attr)
Definition: svgfig.py:482
JetScaleFactorVsEta< DOWN > JetScaleFactorVsEtaDOWN