File indexing completed on 2024-10-08 05:12:15
0001 #include <TFile.h>
0002 #include <TH1.h>
0003 #include <TH2.h>
0004 #include <TPaveStats.h>
0005 #include <TStyle.h>
0006
0007 #include <cassert>
0008 #include <cstdlib>
0009 #include <sstream>
0010
0011 #include "Validation/RecoParticleFlow/interface/Comparator.h"
0012 #include "Validation/RecoParticleFlow/interface/NicePlot.h"
0013 #include "Validation/RecoParticleFlow/interface/TH2Analyzer.h"
0014
0015 using namespace std;
0016
0017 void Comparator::SetDirs(const char *file0, const char *dir0, const char *file1, const char *dir1) {
0018 file0_ = new TFile(file0);
0019 if (file0_->IsZombie())
0020 exit(1);
0021 dir0_ = file0_->GetDirectory(dir0);
0022 if (!dir0_)
0023 exit(1);
0024
0025 file1_ = new TFile(file1);
0026 if (file1_->IsZombie())
0027 exit(1);
0028 dir1_ = file1_->GetDirectory(dir1);
0029 if (!dir1_)
0030 exit(1);
0031 }
0032
0033 void Comparator::SetStyles(Style *s0, Style *s1, const char *leg0, const char *leg1) {
0034 s0_ = s0;
0035 s1_ = s1;
0036
0037 legend_.Clear();
0038 legend_.AddEntry(s0_, leg0, "mlf");
0039 legend_.AddEntry(s1_, leg1, "mlf");
0040 }
0041
0042 void Comparator::DrawSlice(const char *key, int binxmin, int binxmax, Mode mode) {
0043 static int num = 0;
0044
0045 ostringstream out0;
0046 out0 << "h0_2d_" << num;
0047 ostringstream out1;
0048 out1 << "h1_2d_" << num;
0049 num++;
0050
0051 string name0 = out0.str();
0052 string name1 = out1.str();
0053
0054 TH1 *h0 = Histo(key, 0);
0055 TH1 *h1 = Histo(key, 1);
0056
0057 TH2 *h0_2d = dynamic_cast<TH2 *>(h0);
0058 TH2 *h1_2d = dynamic_cast<TH2 *>(h1);
0059
0060 if (h0_2d->GetNbinsY() == 1 || h1_2d->GetNbinsY() == 1) {
0061 cerr << key << " is not 2D" << endl;
0062 return;
0063 }
0064
0065 TH1::AddDirectory(false);
0066
0067 TH1D *h0_slice = h0_2d->ProjectionY(name0.c_str(), binxmin, binxmax, "");
0068 TH1D *h1_slice = h1_2d->ProjectionY(name1.c_str(), binxmin, binxmax, "");
0069 TH1::AddDirectory(true);
0070 Draw(h0_slice, h1_slice, mode);
0071 }
0072
0073 void Comparator::DrawMeanSlice(const char *key, const int rebinFactor, Mode mode) {
0074 TDirectory *dir = dir1_;
0075 dir->cd();
0076 TH2D *h2 = (TH2D *)dir->Get(key);
0077 TH2Analyzer TH2Ana(h2, rebinFactor);
0078 TH1D *ha = TH2Ana.Average();
0079
0080 dir = dir0_;
0081 dir->cd();
0082 TH2D *h2b = (TH2D *)dir->Get(key);
0083 TH2Analyzer TH2Anab(h2b, rebinFactor);
0084 TH1D *hb = TH2Anab.Average();
0085
0086 Draw(hb, ha, mode);
0087 }
0088
0089 void Comparator::DrawSigmaSlice(const char *key, const int rebinFactor, Mode mode) {
0090 TDirectory *dir = dir1_;
0091 dir->cd();
0092 TH2D *h2 = (TH2D *)dir->Get(key);
0093 TH2Analyzer TH2Ana(h2, rebinFactor);
0094 TH1D *ha = TH2Ana.RMS();
0095
0096 dir = dir0_;
0097 dir->cd();
0098 TH2D *h2b = (TH2D *)dir->Get(key);
0099 TH2Analyzer TH2Anab(h2b, rebinFactor);
0100 TH1D *hb = TH2Anab.RMS();
0101
0102 Draw(hb, ha, mode);
0103 }
0104
0105 void Comparator::DrawGaussSigmaSlice(const char *key, const int rebinFactor, Mode mode) {
0106 TDirectory *dir = dir1_;
0107 dir->cd();
0108 TH2D *h2 = (TH2D *)dir->Get(key);
0109 TH2Analyzer TH2Ana(h2, rebinFactor);
0110 TH1D *ha = TH2Ana.SigmaGauss();
0111
0112 dir = dir0_;
0113 dir->cd();
0114 TH2D *h2b = (TH2D *)dir->Get(key);
0115 TH2Analyzer TH2Anab(h2b, rebinFactor);
0116 TH1D *hb = TH2Anab.SigmaGauss();
0117
0118 Draw(hb, ha, mode);
0119 }
0120
0121 Double_t fitFunction_f(Double_t *x, Double_t *par) {
0122 const Double_t value =
0123 sqrt(par[0] * par[0] + par[1] * par[1] * (x[0] - par[3]) + par[2] * par[2] * (x[0] - par[3]) * (x[0] - par[3])) /
0124 x[0];
0125 return value;
0126 }
0127
0128 void Comparator::DrawGaussSigmaSlice(
0129 const char *key, const int rebinFactor, const int binxmin, const int binxmax, const bool cst_binning, Mode mode) {
0130 TDirectory *dir = dir1_;
0131 dir->cd();
0132 TH2D *h2 = (TH2D *)dir->Get(key);
0133 TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
0134 TH1D *hrms = TH2Ana.RMS();
0135 TF1 *fitfcndgssrms3 = new TF1("fitfcndgssrms3",
0136 fitFunction_f,
0137 hrms->GetXaxis()->GetBinLowEdge(1),
0138 hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),
0139 4);
0140 fitfcndgssrms3->SetNpx(500);
0141 fitfcndgssrms3->SetLineWidth(3);
0142 fitfcndgssrms3->SetLineStyle(2);
0143 fitfcndgssrms3->SetLineColor(4);
0144 hrms->Fit(fitfcndgssrms3, "0R");
0145
0146 TH1D *ha = TH2Ana.SigmaGauss();
0147 TF1 *fitfcndgsse3 = new TF1("fitfcndgsse3",
0148 fitFunction_f,
0149 hrms->GetXaxis()->GetBinLowEdge(1),
0150 hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),
0151 4);
0152 fitfcndgsse3->SetNpx(500);
0153 fitfcndgsse3->SetLineWidth(3);
0154 fitfcndgsse3->SetLineStyle(1);
0155 fitfcndgsse3->SetLineColor(4);
0156 ha->Fit(fitfcndgsse3, "0R");
0157
0158 dir = dir0_;
0159 dir->cd();
0160 TH2D *h2b = (TH2D *)dir->Get(key);
0161 TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
0162 TH1D *hrmsb = TH2Anab.RMS();
0163 TF1 *fitfcndgssrmsb3 = new TF1("fitfcndgssrmsb3",
0164 fitFunction_f,
0165 hrms->GetXaxis()->GetBinLowEdge(1),
0166 hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),
0167 4);
0168 fitfcndgssrmsb3->SetNpx(500);
0169 fitfcndgssrmsb3->SetLineWidth(3);
0170 fitfcndgssrmsb3->SetLineStyle(2);
0171 fitfcndgssrmsb3->SetLineColor(2);
0172 hrmsb->Fit(fitfcndgssrmsb3, "0R");
0173
0174 TH1D *hb = TH2Anab.SigmaGauss();
0175 TF1 *fitfcndgsseb3 = new TF1("fitfcndgsseb3",
0176 fitFunction_f,
0177 hrms->GetXaxis()->GetBinLowEdge(1),
0178 hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),
0179 4);
0180 fitfcndgsseb3->SetNpx(500);
0181 fitfcndgsseb3->SetLineWidth(3);
0182 fitfcndgsseb3->SetLineStyle(1);
0183 fitfcndgsseb3->SetLineColor(2);
0184 hb->Fit(fitfcndgsseb3, "0R");
0185
0186 Draw(hb, ha, mode);
0187
0188
0189 fitfcndgssrms3->Draw("same");
0190 fitfcndgsse3->Draw("same");
0191 fitfcndgssrmsb3->Draw("same");
0192 fitfcndgsseb3->Draw("same");
0193 }
0194
0195 void Comparator::DrawGaussSigmaOverMeanXSlice(
0196 const char *key, const int rebinFactor, const int binxmin, const int binxmax, const bool cst_binning, Mode mode) {
0197 TDirectory *dir = dir1_;
0198 dir->cd();
0199 TH2D *h2 = (TH2D *)dir->Get(key);
0200 TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
0201 TH1D *hrms = TH2Ana.RMS();
0202
0203 TH1D *meanXslice = TH2Ana.MeanX();
0204
0205
0206
0207
0208
0209
0210
0211 hrms->Divide(meanXslice);
0212
0213 TF1 *fitXfcndgssrms3 = new TF1("fitXfcndgssrms3",
0214 fitFunction_f,
0215 hrms->GetXaxis()->GetBinLowEdge(1),
0216 hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),
0217 4);
0218 fitXfcndgssrms3->SetNpx(500);
0219 fitXfcndgssrms3->SetLineWidth(3);
0220 fitXfcndgssrms3->SetLineStyle(2);
0221 fitXfcndgssrms3->SetLineColor(4);
0222 hrms->Fit(fitXfcndgssrms3, "0R");
0223
0224 TH1D *ha = TH2Ana.SigmaGauss();
0225 ha->Divide(meanXslice);
0226 TF1 *fitXfcndgsse3 = new TF1("fitXfcndgsse3",
0227 fitFunction_f,
0228 ha->GetXaxis()->GetBinLowEdge(1),
0229 ha->GetXaxis()->GetBinUpEdge(ha->GetNbinsX()),
0230 4);
0231 fitXfcndgsse3->SetNpx(500);
0232 fitXfcndgsse3->SetLineWidth(3);
0233 fitXfcndgsse3->SetLineStyle(1);
0234 fitXfcndgsse3->SetLineColor(4);
0235 ha->Fit(fitXfcndgsse3, "0R");
0236
0237 dir = dir0_;
0238 dir->cd();
0239 TH2D *h2b = (TH2D *)dir->Get(key);
0240 TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
0241 TH1D *hrmsb = TH2Anab.RMS();
0242 hrmsb->Divide(meanXslice);
0243 TF1 *fitXfcndgssrmsb3 = new TF1("fitXfcndgssrmsb3",
0244 fitFunction_f,
0245 hrmsb->GetXaxis()->GetBinLowEdge(1),
0246 hrmsb->GetXaxis()->GetBinUpEdge(hrmsb->GetNbinsX()),
0247 4);
0248 fitXfcndgssrmsb3->SetNpx(500);
0249 fitXfcndgssrmsb3->SetLineWidth(3);
0250 fitXfcndgssrmsb3->SetLineStyle(2);
0251 fitXfcndgssrmsb3->SetLineColor(2);
0252 hrmsb->Fit(fitXfcndgssrmsb3, "0R");
0253
0254 TH1D *hb = TH2Anab.SigmaGauss();
0255 hb->Divide(meanXslice);
0256 TF1 *fitXfcndgsseb3 = new TF1("fitXfcndgsseb3",
0257 fitFunction_f,
0258 hb->GetXaxis()->GetBinLowEdge(1),
0259 hb->GetXaxis()->GetBinUpEdge(hb->GetNbinsX()),
0260 4);
0261 fitXfcndgsseb3->SetNpx(500);
0262 fitXfcndgsseb3->SetLineWidth(3);
0263 fitXfcndgsseb3->SetLineStyle(1);
0264 fitXfcndgsseb3->SetLineColor(2);
0265 hb->Fit(fitXfcndgsseb3, "0R");
0266
0267 Draw(hb, ha, mode);
0268
0269
0270 fitXfcndgssrms3->Draw("same");
0271 fitXfcndgsse3->Draw("same");
0272 fitXfcndgssrmsb3->Draw("same");
0273 fitXfcndgsseb3->Draw("same");
0274 }
0275
0276 void Comparator::DrawGaussSigmaOverMeanSlice(const char *key, const char *key2, const int rebinFactor, Mode mode) {
0277 TDirectory *dir = dir1_;
0278 dir->cd();
0279
0280 TH2D *h2_b = (TH2D *)dir->Get(key2);
0281 TH2Analyzer TH2Ana_b(h2_b, rebinFactor);
0282 TH1D *meanslice = TH2Ana_b.Average();
0283
0284 TH2D *h2 = (TH2D *)dir->Get(key);
0285 TH2Analyzer TH2Ana(h2, rebinFactor);
0286
0287 TH1D *ha = TH2Ana.SigmaGauss();
0288 ha->Divide(meanslice);
0289
0290 dir = dir0_;
0291 dir->cd();
0292 TH2D *h2b = (TH2D *)dir->Get(key);
0293 TH2Analyzer TH2Anab(h2b, rebinFactor);
0294
0295 TH2D *h2b_b = (TH2D *)dir->Get(key2);
0296 TH2Analyzer TH2Anab_b(h2b_b, rebinFactor);
0297 TH1D *meansliceb = TH2Anab_b.Average();
0298
0299 TH1D *hb = TH2Anab.SigmaGauss();
0300 hb->Divide(meansliceb);
0301
0302 Draw(hb, ha, mode);
0303
0304 }
0305
0306 void Comparator::Draw(const char *key, Mode mode) {
0307 TH1::AddDirectory(false);
0308 TH1 *h0 = Histo(key, 0);
0309 TH1 *h1 = (TH1 *)Histo(key, 1)->Clone("h1");
0310
0311 TH1::AddDirectory(true);
0312 Draw(h0, h1, mode);
0313 }
0314
0315 void Comparator::Draw(const char *key0, const char *key1, Mode mode) {
0316 TH1 *h0 = nullptr;
0317 TH1 *h1 = nullptr;
0318 if (mode != EFF) {
0319 h0 = Histo(key0, 0);
0320 h1 = Histo(key1, 1);
0321 } else {
0322 h0 = Histo(key0, 0);
0323 TH1 *h0b = Histo(key1, 0);
0324 h1 = Histo(key0, 1);
0325 TH1 *h1b = Histo(key1, 1);
0326 if (rebin_ > 1) {
0327 h0->Rebin(rebin_);
0328 h1->Rebin(rebin_);
0329 h0b->Rebin(rebin_);
0330 h1b->Rebin(rebin_);
0331 }
0332 if (resetAxis_) {
0333 h0->GetXaxis()->SetRangeUser(xMin_, xMax_);
0334 h1->GetXaxis()->SetRangeUser(xMin_, xMax_);
0335 h0b->GetXaxis()->SetRangeUser(xMin_, xMax_);
0336 h1b->GetXaxis()->SetRangeUser(xMin_, xMax_);
0337 }
0338
0339 h0b->Sumw2();
0340 h0->Sumw2();
0341 h0->Divide(h0, h0b, 1., 1., "B");
0342 h1b->Sumw2();
0343 h1->Sumw2();
0344 h1->Divide(h1, h1b, 1., 1., "B");
0345 }
0346 Draw(h0, h1, mode);
0347 }
0348
0349 TH1 *Comparator::Histo(const char *key, unsigned dirIndex) {
0350 if (dirIndex > 1U) {
0351 cerr << "bad dir index: " << dirIndex << endl;
0352 return nullptr;
0353 }
0354 TDirectory *dir = nullptr;
0355 if (dirIndex == 0)
0356 dir = dir0_;
0357 if (dirIndex == 1)
0358 dir = dir1_;
0359 assert(dir);
0360
0361 dir->cd();
0362
0363 TH1 *h = (TH1 *)dir->Get(key);
0364 if (!h)
0365 cerr << "no key " << key << " in directory " << dir->GetName() << endl;
0366 return h;
0367 }
0368
0369 void Comparator::Draw(TH1 *h0, TH1 *h1, Mode mode) {
0370 if (!(h0 && h1)) {
0371 cerr << "invalid histo" << endl;
0372 return;
0373 }
0374
0375 TH1::AddDirectory(false);
0376 h0_ = (TH1 *)h0->Clone("h0_");
0377 h1_ = (TH1 *)h1->Clone("h1_");
0378 TH1::AddDirectory(true);
0379
0380
0381
0382
0383 h0_->SetTitle("");
0384 h1_->SetTitle("");
0385
0386
0387
0388
0389 if (mode != EFF) {
0390 if (rebin_ > 1) {
0391 h0_->Rebin(rebin_);
0392 h1_->Rebin(rebin_);
0393 }
0394 if (resetAxis_) {
0395 h0_->GetXaxis()->SetRangeUser(xMin_, xMax_);
0396 h1_->GetXaxis()->SetRangeUser(xMin_, xMax_);
0397 }
0398 }
0399
0400 if (mode != GRAPH) {
0401 TPaveStats *ptstats = new TPaveStats(0.7385057, 0.720339, 0.9396552, 0.8792373, "brNDC");
0402 ptstats->SetName("stats");
0403 ptstats->SetBorderSize(1);
0404 ptstats->SetLineColor(2);
0405 ptstats->SetFillColor(10);
0406 ptstats->SetTextAlign(12);
0407 ptstats->SetTextColor(2);
0408 ptstats->SetOptStat(1111);
0409 ptstats->SetOptFit(0);
0410 ptstats->Draw();
0411 h0_->GetListOfFunctions()->Add(ptstats);
0412 ptstats->SetParent(h0_->GetListOfFunctions());
0413
0414
0415
0416
0417
0418
0419 TPaveStats *ptstats2 = new TPaveStats(0.7399425, 0.529661, 0.941092, 0.6885593, "brNDC");
0420 ptstats2->SetName("stats");
0421 ptstats2->SetBorderSize(1);
0422 ptstats2->SetLineColor(4);
0423 ptstats2->SetFillColor(10);
0424 ptstats2->SetTextAlign(12);
0425 ptstats2->SetTextColor(4);
0426 TText *text = ptstats2->AddText("h1_");
0427 text->SetTextSize(0.03654661);
0428
0429 std::ostringstream oss3;
0430 oss3 << h1_->GetEntries();
0431 const std::string txt_entries = "Entries = " + oss3.str();
0432 ptstats2->AddText(txt_entries.c_str());
0433 std::ostringstream oss;
0434 oss << h1_->GetMean();
0435 const std::string txt_mean = "Mean = " + oss.str();
0436 ptstats2->AddText(txt_mean.c_str());
0437 std::ostringstream oss2;
0438 oss2 << h1_->GetRMS();
0439 const std::string txt_rms = "RMS = " + oss2.str();
0440 ptstats2->AddText(txt_rms.c_str());
0441 ptstats2->SetOptStat(1111);
0442 ptstats2->SetOptFit(0);
0443 ptstats2->Draw();
0444 h1_->GetListOfFunctions()->Add(ptstats2);
0445 ptstats2->SetParent(h1_->GetListOfFunctions());
0446 } else {
0447 TPaveStats *ptstats = new TPaveStats(0.0, 0.0, 0.0, 0.0, "brNDC");
0448 ptstats->Draw();
0449 h0_->GetListOfFunctions()->Add(ptstats);
0450 ptstats->SetParent(h0_->GetListOfFunctions());
0451 }
0452
0453 float min = -999.;
0454 float max = +999.;
0455 switch (mode) {
0456 case SCALE:
0457 h1_->Scale(h0_->GetEntries() / h1_->GetEntries());
0458 break;
0459 case NORMAL:
0460 if (s0_)
0461 Styles::FormatHisto(h0_, s0_);
0462 if (s1_)
0463 Styles::FormatHisto(h1_, s1_);
0464
0465 if (h1_->GetMaximum() > h0_->GetMaximum()) {
0466 h0_->SetMaximum(h1_->GetMaximum() * 1.15);
0467 }
0468
0469 h0_->Draw();
0470 h1_->Draw("same");
0471
0472 break;
0473 case EFF:
0474 if (s0_)
0475 Styles::FormatHisto(h0_, s0_);
0476 if (s1_)
0477 Styles::FormatHisto(h1_, s1_);
0478
0479
0480 max = h0_->GetMaximum();
0481 if (h1_->GetMaximum() > max)
0482 max = h1_->GetMaximum();
0483 if (max > 0.8)
0484 max = 1;
0485
0486 max *= 1.1;
0487
0488 min = h0_->GetMinimum();
0489 if (h1_->GetMinimum() < min)
0490 min = h1_->GetMinimum();
0491 if (min > 0.2)
0492 min = 0.;
0493
0494 min *= 0.8;
0495 h0_->SetMaximum(max);
0496 h0_->SetMinimum(min);
0497
0498 h0_->Draw("E");
0499 h1_->Draw("Esame");
0500 break;
0501 case GRAPH:
0502 if (s0_)
0503 Styles::FormatHisto(h0_, s0_);
0504 if (s1_)
0505 Styles::FormatHisto(h1_, s1_);
0506
0507 if (h1_->GetMaximum() > h0_->GetMaximum()) {
0508 h0_->SetMaximum(h1_->GetMaximum() * 1.15);
0509 }
0510 if (h1_->GetMinimum() < h0_->GetMinimum()) {
0511 h0_->SetMinimum(h1_->GetMinimum() * 1.15);
0512 }
0513 h0_->SetMarkerStyle(21);
0514 h0_->SetMarkerColor(2);
0515 h1_->SetMarkerStyle(21);
0516 h1_->SetMarkerColor(4);
0517
0518 h0_->Draw("E1");
0519 h1_->Draw("E1same");
0520 break;
0521 case RATIO:
0522 h0_->Sumw2();
0523 h1_->Sumw2();
0524 h0_->Divide(h1_);
0525 if (s0_)
0526 Styles::FormatHisto(h0_, s0_);
0527 h0_->Draw();
0528 break;
0529 default:
0530 break;
0531 }
0532 }