RecoParticleFlow/src/Comparator.cc

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   // Draw(hrms,ha,mode);
0188   // Draw(ha,ha,mode);
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   // for( int i=1; i<=meanXslice->GetNbinsX(); ++i) {
0205   //  std::cout << "meanXslice->GetBinContent(" << i << ") = "
0206   //          << meanXslice->GetBinContent(i) << std::endl;
0207   //  std::cout << "meanXslice->GetBinError(" << i << ") = "
0208   //          << meanXslice->GetBinError(i) << std::endl;
0209   //}
0210   // Draw(meanXslice,meanXslice,mode);
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   // Draw(hrms,ha,mode);
0269   // Draw(ha,ha,mode);
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   // Draw(meansliceb,meanslice,mode);
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) {  // dirIndex >= 0, since dirIndex is unsigned
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   // unsetting the title, since the title of projections
0381   // is still the title of the 2d histo
0382   // and this is better anyway
0383   h0_->SetTitle("");
0384   h1_->SetTitle("");
0385
0386   // h0_->SetStats(1);
0387   // h1_->SetStats(1);
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     // std::cout << "FL: h0_->GetMean() = " << h0_->GetMean() << std::endl;
0415     // std::cout << "FL: h0_->GetRMS() = " << h0_->GetRMS() << std::endl;
0416     // std::cout << "FL: h1_->GetMean() = " << h1_->GetMean() << std::endl;
0417     // std::cout << "FL: h1_->GetRMS() = " << h1_->GetRMS() << std::endl;
0418     // std::cout << "FL: test2" << std::endl;
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     text = ptstats2->AddText(txt_entries.c_str());
0433     std::ostringstream oss;
0434     oss << h1_->GetMean();
0435     const std::string txt_mean = "Mean  = " + oss.str();
0436     text = ptstats2->AddText(txt_mean.c_str());
0437     std::ostringstream oss2;
0438     oss2 << h1_->GetRMS();
0439     const std::string txt_rms = "RMS  = " + oss2.str();
0440     text = 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       // rather arbitrary but useful
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 }