Tracking/BenchmarkingTools/CommonTools.h

0001 void myText(Double_t x,Double_t y,Color_t color,const char *text, Double_t tsize = 0.05, double angle = -1);
0002 void DivideCanvas(TVirtualPad *p, int npads);
0003 TGraphErrors* FitProfile(const TH2 *h2);
0004
0005 void myText(Double_t x,Double_t y,Color_t color,
0006         const char *text, Double_t tsize, double angle) {
0007
0008   TLatex l; //l.SetTextAlign(12);
0009   l.SetTextSize(tsize);
0010   l.SetNDC();
0011   l.SetTextColor(color);
0012   if (angle > 0) l.SetTextAngle(angle);
0013   l.DrawLatex(x,y,text);
0014 }
0015
0016 //! Draw a horizontal line in a pad at given x coordinate
0017 TLine *HorizontalLine(TVirtualPad *p, Double_t y)
0018 {
0019   p->Update();
0020
0021   Double_t xMin = p->GetUxmin();
0022   Double_t xMax = p->GetUxmax();
0023
0024   if (p->GetLogx())
0025   {
0026     xMin = std::pow(10, xMin);
0027     xMax = std::pow(10, xMax);
0028   }
0029
0030   TLine *line = new TLine(xMin, y, xMax, y);
0031   line->SetLineStyle(2);
0032   line->SetLineWidth(1);
0033   line->SetLineColor(1);
0034   return line;
0035 }
0036
0037 //! Fit for resolution of TH2F
0038 TGraphErrors *
0039 FitResolution(const TH2F *h2, const bool normalize_mean = true, const int param = 2)
0040 {
0041   TProfile *p2 = h2->ProfileX();
0042
0043   int n = 0;
0044   double x[1000];
0045   double ex[1000];
0046   double y[1000];
0047   double ey[1000];
0048
0049   for (int i = 1; i <= h2->GetNbinsX(); i++)
0050   {
0051     TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0052
0053     if (h1->GetSum() < 10)
0054       continue;
0055
0056     TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0057               p2->GetBinError(i) * 4);
0058
0059     TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0060            -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0061
0062     fgaus.SetParameter(1, p2->GetBinContent(i));
0063     fgaus.SetParameter(2, p2->GetBinError(i));
0064
0065     h1->Fit(&fgaus, "MQ0");
0066
0067     x[n] = p2->GetBinCenter(i);
0068     ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0069
0070     const double norm = normalize_mean ? fgaus.GetParameter(1) : 1;
0071
0072     y[n] = fgaus.GetParameter(param) / norm;
0073     ey[n] = fgaus.GetParError(param) / norm;
0074
0075     n++;
0076     delete h1;
0077   }
0078
0079   TGraphErrors *ge = new TGraphErrors(n, x, y, 0, ey);
0080   ge->SetName(TString(h2->GetName()) + "_FitResolution");
0081
0082   ge->SetLineColor(kBlue + 3);
0083   ge->SetMarkerColor(kBlue + 3);
0084   ge->SetLineWidth(2);
0085   ge->SetMarkerStyle(kFullCircle);
0086   ge->SetMarkerSize(1);
0087   return ge;
0088 }
0089
0090 //! Fit for profile along the Y direction of TH2F
0091 TGraphErrors *
0092 FitProfile(const TH2 *h2)
0093 {
0094   TProfile *p2 = h2->ProfileX();
0095
0096   int n = 0;
0097   double x[1000];
0098   double ex[1000];
0099   double y[1000];
0100   double ey[1000];
0101
0102   for (int i = 1; i <= h2->GetNbinsX(); i++)
0103   {
0104     TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0105
0106     if (h1->GetSum() < 10)
0107       continue;
0108
0109     TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0110               p2->GetBinError(i) * 4);
0111
0112     TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0113            -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0114
0115     fgaus.SetParameter(1, p2->GetBinContent(i));
0116     fgaus.SetParameter(2, p2->GetBinError(i));
0117
0118     h1->Fit(&fgaus, "MQ0");
0119
0120     //      f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
0121     //          fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
0122     //          fgaus.GetParameter(2) / 4, 0);
0123     //
0124     //      h1->Fit(&f2, "MQ0");
0125
0126     //      new TCanvas;
0127     //      h1->Draw();
0128     //      fgaus.Draw("same");
0129     //      break;
0130
0131     x[n] = p2->GetBinCenter(i);
0132     ex[n] = p2->GetBinWidth(i) / 2;
0133     y[n] = fgaus.GetParameter(1);
0134     ey[n] = fgaus.GetParameter(2);
0135
0136     //      p2->SetBinContent(i, fgaus.GetParameter(1));
0137     //      p2->SetBinError(i, fgaus.GetParameter(2));
0138
0139     n++;
0140     delete h1;
0141   }
0142
0143   TGraphErrors *ge = new TGraphErrors(n, x, y, ex, ey);
0144   ge->SetName(TString(h2->GetName()) + "_FitProfile");
0145   ge->SetLineColor(kBlue + 3);
0146   ge->SetMarkerColor(kBlue + 3);
0147   ge->SetLineWidth(2);
0148   ge->SetMarkerStyle(kFullCircle);
0149   ge->SetMarkerSize(1);
0150   return ge;
0151 }
0152
0153 //!ratio between two histograms with binominal error based on Wilson score interval. Assuming each histogram is count.
0154 TH1 *GetBinominalRatio(TH1 *h_pass, TH1 *h_n_trial, bool process_zero_bins = false)
0155 {
0156   assert(h_pass);
0157   assert(h_n_trial);
0158
0159   assert(h_pass->GetNbinsX() == h_n_trial->GetNbinsX());
0160   assert(h_pass->GetNbinsY() == h_n_trial->GetNbinsY());
0161   assert(h_pass->GetNbinsZ() == h_n_trial->GetNbinsZ());
0162
0163   TH1 *h_ratio = (TH1 *) h_pass->Clone(TString(h_pass->GetName()) + "_Ratio");
0164   assert(h_ratio);
0165   h_ratio->Divide(h_n_trial);  // a rough estimation first, also taking care of the overflow bins and zero bins
0166
0167   for (int x = 1; x <= h_n_trial->GetNbinsX(); ++x)
0168     for (int y = 1; y <= h_n_trial->GetNbinsY(); ++y)
0169       for (int z = 1; z <= h_n_trial->GetNbinsZ(); ++z)
0170       {
0171         const double n_trial = h_n_trial->GetBinContent(x, y, z);
0172
0173         if (n_trial > 0)
0174         {
0175           const double p = h_pass->GetBinContent(x, y, z) / n_trial;
0176
0177           // Wilson score interval
0178           h_ratio->SetBinContent(x, y, z,  //
0179                                  (p + 1 / (2 * n_trial)) / (1 + 1 / n_trial));
0180           h_ratio->SetBinError(x, y,
0181                                z,  //
0182                                std::sqrt(
0183                                    1. / n_trial * p * (1 - p) + 1. / (4 * n_trial * n_trial)) /
0184                                    (1 + 1 / n_trial));
0185         }
0186         else if (process_zero_bins)
0187         {
0188           h_ratio->SetBinContent(x, y, z, 0.5);
0189           h_ratio->SetBinError(x, y, z, 0.5);
0190         }
0191       }
0192
0193   return h_ratio;
0194 }
0195
0196
0197 //! Draw a vertical line in a pad at given x coordinate
0198 TLine *VerticalLine(TVirtualPad *p, Double_t x)
0199 {
0200   p->Update();
0201
0202   Double_t yMin = p->GetUymin();
0203   Double_t yMax = p->GetUymax();
0204
0205   if (p->GetLogy())
0206   {
0207     yMin = std::pow(10, yMin);
0208     yMax = std::pow(10, yMax);
0209   }
0210
0211   TLine *line = new TLine(x, yMin, x, yMax);
0212   line->SetLineStyle(2);
0213   line->SetLineWidth(1);
0214   line->SetLineColor(1);
0215   return line;
0216 }
0217 double DrawReference(TH1 *hnew, TH1 *href, bool draw_href_error = false, bool do_kstest = true)
0218 {
0219   hnew->SetLineColor(kBlue + 3);
0220   hnew->SetMarkerColor(kBlue + 3);
0221   //  hnew->SetLineWidth(2);
0222   hnew->SetMarkerStyle(kFullCircle);
0223   //  hnew->SetMarkerSize(1);
0224
0225   if (href)
0226   {
0227     if (draw_href_error)
0228     {
0229       href->SetLineColor(kGreen + 1);
0230       href->SetFillColor(kGreen + 1);
0231       href->SetLineStyle(kSolid);
0232       href->SetMarkerColor(kGreen + 1);
0233       //      href->SetLineWidth(2);
0234       href->SetMarkerStyle(kDot);
0235       href->SetMarkerSize(0);
0236     }
0237     else
0238     {
0239       href->SetLineColor(kGreen + 1);
0240       href->SetFillColor(kGreen + 1);
0241       href->SetLineStyle(0);
0242       href->SetMarkerColor(kGreen + 1);
0243       href->SetLineWidth(0);
0244       href->SetMarkerStyle(kDot);
0245       href->SetMarkerSize(0);
0246     }
0247   }
0248
0249   hnew->Draw();  // set scale
0250
0251   double ks_test = numeric_limits<double>::signaling_NaN();
0252
0253   if (href)
0254   {
0255     if (draw_href_error)
0256     {
0257       href->DrawClone("E2 same");
0258       href->SetFillStyle(0);
0259       href->SetLineWidth(8);
0260       href->DrawClone("HIST same ][");
0261     }
0262     else
0263       href->Draw("HIST same");
0264     hnew->Draw("same");  // over lay data points
0265
0266     if (do_kstest)
0267       ks_test = hnew->KolmogorovTest(href, "");
0268   }
0269
0270   // ---------------------------------
0271   // now, make summary header
0272   // ---------------------------------
0273   if (href)
0274   {
0275     gPad->SetTopMargin(.14);
0276     TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
0277     legend->Draw();
0278     legend = new TLegend(0, .86, .3, .93, NULL, "NB NDC");
0279     legend->AddEntry(href, Form("Reference"), "f");
0280     legend->Draw();
0281     legend = new TLegend(0.3, .86, 1, .93, NULL, "NB NDC");
0282
0283     if (do_kstest)
0284     {
0285       TLegendEntry *le = legend->AddEntry(hnew, Form("New: KS-Test P=%.3f", ks_test), "lpe");
0286       if (ks_test >= 1)
0287         le->SetTextColor(kBlue + 1);
0288       else if (ks_test >= .2)
0289         le->SetTextColor(kGreen + 2);
0290       else if (ks_test >= .05)
0291         le->SetTextColor(kYellow + 1);
0292       else
0293         le->SetTextColor(kRed + 1);
0294       legend->Draw();
0295     }
0296     else
0297     {
0298       TLegendEntry *le = legend->AddEntry(hnew, Form("New Result"), "lpe");
0299       legend->Draw();
0300     }
0301   }
0302   else
0303   {
0304     gPad->SetTopMargin(.07);
0305     TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
0306     legend->Draw();
0307   }
0308
0309
0310   return ks_test;
0311 }
0312
0313 //! Draw 1D TGraph along with its reference as shade
0314 //! @param[in] draw_href_error whether to draw error band for reference plot. Otherwise, it is a filled histogram (default)
0315 void DrawReference(TGraph *hnew, TGraph *href, bool draw_href_error = true)
0316 {
0317   hnew->SetLineColor(kBlue + 3);
0318   hnew->SetMarkerColor(kBlue + 3);
0319   hnew->SetLineWidth(2);
0320   hnew->SetMarkerStyle(kFullCircle);
0321   hnew->SetMarkerSize(1);
0322
0323   if (href)
0324   {
0325     if (draw_href_error)
0326     {
0327       href->SetLineColor(kGreen + 1);
0328       href->SetFillColor(kGreen + 1);
0329       href->SetFillStyle(0);
0330       href->SetLineStyle(kSolid);
0331       href->SetMarkerColor(kGreen + 1);
0332       href->SetLineWidth(4);
0333       href->SetMarkerStyle(kDot);
0334       href->SetMarkerSize(0);
0335     }
0336     else
0337     {
0338       href->SetLineColor(kGreen + 1);
0339       href->SetFillColor(kGreen + 1);
0340       href->SetLineStyle(0);
0341       href->SetMarkerColor(kGreen + 1);
0342       href->SetLineWidth(0);
0343       href->SetMarkerStyle(kDot);
0344       href->SetMarkerSize(0);
0345     }
0346   }
0347
0348   if (href)
0349   {
0350     if (draw_href_error)
0351     {
0352       href->DrawClone("E2");
0353     }
0354     else
0355       href->Draw("HIST same");
0356     hnew->Draw("p e");  // over lay data points
0357   }
0358
0359   // ---------------------------------
0360   // now, make summary header
0361   // ---------------------------------
0362   if (href)
0363   {
0364     gPad->SetTopMargin(.14);
0365     TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
0366     legend->Draw();
0367     legend = new TLegend(0, .86, .3, .93, NULL, "NB NDC");
0368     legend->AddEntry(href, Form("Reference"), "f");
0369     legend->Draw();
0370     legend = new TLegend(0.3, .86, 1, .93, NULL, "NB NDC");
0371     TLegendEntry *le = legend->AddEntry(hnew, Form("New"), "lpe");
0372     legend->Draw();
0373   }
0374   else
0375   {
0376     gPad->SetTopMargin(.07);
0377     TLegend *legend = new TLegend(0, .93, 0, 1, hnew->GetTitle(), "NB NDC");
0378     legend->Draw();
0379   }
0380 }
0381
0382
0383 //! Divide canvas in a given number of pads, with a number of pads in both directions that match the width/height ratio of the canvas
0384 void DivideCanvas(TVirtualPad *p, int npads)
0385 {
0386   if (!p) return;
0387   if (!npads) return;
0388   const double ratio = double(p->GetWw()) / p->GetWh();
0389   Int_t columns = std::max(1, int(std::sqrt(npads * ratio)));
0390   Int_t rows = std::max(1, int(npads / columns));
0391   while (rows * columns < npads)
0392   {
0393     columns++;
0394     if (rows * columns < npads) rows++;
0395   }
0396   p->Divide(rows, columns);
0397 }