PatExamples/test/patBJetTracks_efficiencies.C

0001 // name of analyzer
0002 static const char *directory = "analyzeBJetTracks";
0003
0004 static const char *flavours[] = { "b", "c", "udsg", 0 };
0005
0006 void patBJetTracks_efficiencies()
0007 {
0008     // define proper canvas style
0009     setNiceStyle();
0010     gStyle->SetOptStat(0);
0011
0012     // open file
0013     TFile* file = new TFile("analyzePatBJetTracks.root");
0014
0015     TLegend *legend[3] = { 0, 0, 0 };
0016
0017     // draw canvas with efficiencies
0018
0019     TCanvas *canv;
0020     canv = new TCanvas("canv0", "hand-crafted track counting efficiencies", 800, 300);
0021     canv->Divide(3, 1);
0022
0023     TH1 *total = (TH1*)file->Get(Form("%s/flavours", directory));
0024     TH1 *effVsCutB = 0;
0025     unsigned int i = 0;
0026     for(const char **flavour = flavours; *flavour; flavour++, i++) {
0027         TH1 *h = (TH1*)file->Get(Form("%s/trackIPSig_%s", directory, *flavour));
0028         TH1 *discrShape = (TH1*)h->Clone(Form("%s_discrShape", h->GetName()));
0029         discrShape->Scale(1.0 / discrShape->Integral());
0030         discrShape->SetMaximum(discrShape->GetMaximum() * 5);
0031         TH1 *effVsCut = computeEffVsCut(h, total->GetBinContent(4 - i));
0032         TH1 *effVsBEff = 0;
0033
0034         if (flavour == flavours)    // b-jets
0035             effVsCutB = effVsCut;
0036         else
0037             effVsBEff = computeEffVsBEff(effVsCut, effVsCutB);
0038
0039         discrShape->SetTitle("discriminator shape");
0040         effVsCut->SetTitle("efficiency versus discriminator cut");
0041         if (effVsBEff)
0042             effVsBEff->SetTitle("mistag versus b efficiency");
0043
0044         setHistStyle(discrShape);
0045         setHistStyle(effVsCut);
0046         setHistStyle(effVsBEff);
0047
0048         canv->cd(1);
0049         gPad->SetLogy(1);
0050         gPad->SetGridy(1);
0051         discrShape->SetLineColor(i + 1);
0052         discrShape->SetMarkerColor(i + 1);
0053         discrShape->Draw(i > 0 ? "same" : "");
0054         if (!legend[0])
0055             legend[0] = new TLegend(0.5, 0.7, 0.78, 0.88);
0056         legend[0]->AddEntry(discrShape, *flavour);
0057
0058         canv->cd(2);
0059         gPad->SetLogy(1);
0060         gPad->SetGridy(1);
0061         effVsCut->SetLineColor(i + 1);
0062         effVsCut->SetMarkerColor(i + 1);
0063         effVsCut->Draw(i > 0 ? "same" : "");
0064         if (!legend[1])
0065             legend[1] = new TLegend(0.12, 0.12, 0.40, 0.30);
0066         legend[1]->AddEntry(effVsCut, *flavour);
0067
0068         if (!effVsBEff)
0069             continue;
0070         canv->cd(3);
0071         gPad->SetLogy(1);
0072         gPad->SetGridx(1);
0073         gPad->SetGridy(1);
0074         effVsBEff->SetLineColor(i + 1);
0075         effVsBEff->SetMarkerColor(i + 1);
0076         effVsBEff->Draw(i > 1 ? "same" : "");
0077         if (!legend[2])
0078             legend[2] = new TLegend(0.12, 0.7, 0.40, 0.88);
0079         legend[2]->AddEntry(effVsBEff, *flavour);
0080     }
0081
0082     canv->cd(1);
0083     legend[0]->Draw();
0084
0085     canv->cd(2);
0086     legend[1]->Draw();
0087
0088     canv->cd(3);
0089     legend[2]->Draw();
0090
0091     ////////////////////////////////////////////
0092
0093     // canvas to compare negative tagger with light flavour mistag
0094
0095     TCanvas *canv;
0096     canv = new TCanvas("canv1", "comparing light flavour mistag with negative tagger", 530, 300);
0097     canv->Divide(2, 1);
0098
0099     TH1 *h1 = (TH1*)file->Get(Form("%s/trackIPSig_udsg", directory));
0100     TH1 *h2 = (TH1*)file->Get(Form("%s/negativeIPSig_all", directory));
0101     h2 = invertHisto(h2);   // invert x-axis
0102
0103     TH1 *discrShape1 = (TH1*)h1->Clone("discrShape1");
0104     TH1 *discrShape2 = (TH1*)h2->Clone("discrShape2");
0105
0106     discrShape1->Scale(1.0 / discrShape1->Integral());
0107     discrShape1->SetMaximum(discrShape1->GetMaximum() * 5);
0108     discrShape2->Scale(1.0 / discrShape2->Integral());
0109
0110     TH1 *effVsCut1 = computeEffVsCut(h1, total->GetBinContent(2));
0111     TH1 *effVsCut2 = computeEffVsCut(h2, total->GetBinContent(1));
0112
0113     discrShape1->SetTitle("discriminator shape");
0114     effVsCut1->SetTitle("efficiency versus discriminator cut");
0115
0116     setHistStyle(discrShape1);
0117     setHistStyle(discrShape2);
0118     setHistStyle(effVsCut1);
0119     setHistStyle(effVsCut2);
0120
0121     canv->cd(1);
0122     gPad->SetLogy(1);
0123     gPad->SetGridy(1);
0124     discrShape1->SetLineColor(1);
0125     discrShape1->SetMarkerColor(1);
0126     discrShape2->SetLineColor(2);
0127     discrShape2->SetMarkerColor(2);
0128
0129     discrShape1->Draw();
0130     discrShape2->Draw("same");
0131
0132     TLegend *l = new TLegend(0.5, 0.7, 0.78, 0.88);
0133     l->AddEntry(discrShape1, "udsg");
0134     l->AddEntry(discrShape2, "inv. neg");
0135     l->Draw();
0136
0137     canv->cd(2);
0138     gPad->SetLogy(1);
0139     gPad->SetGridy(1);
0140     effVsCut1->SetLineColor(1);
0141     effVsCut1->SetMarkerColor(1);
0142     effVsCut2->SetLineColor(2);
0143     effVsCut2->SetMarkerColor(2);
0144
0145     effVsCut1->Draw();
0146     effVsCut2->Draw("same");
0147
0148     l = new TLegend(0.5, 0.7, 0.78, 0.88);
0149     l->AddEntry(effVsCut1, "udsg");
0150     l->AddEntry(effVsCut2, "inv. neg");
0151     l->Draw();
0152 }
0153
0154 TH1 *computeEffVsCut(TH1 *discrShape, double total)
0155 {
0156     TH1 *h = discrShape->Clone(Form("%s_effVsCut", discrShape->GetName()));
0157     h->Sumw2();
0158     h->SetMaximum(1.5);
0159     h->SetMinimum(1e-3);
0160
0161     unsigned int n = h->GetNbinsX();
0162     for(unsigned int bin = 1; bin <= n; bin++) {
0163         double efficiency = h->Integral(bin, n + 1) / total;
0164         double error = sqrt(efficiency * (1 - efficiency) / total);
0165         h->SetBinContent(bin, efficiency);
0166         h->SetBinError(bin, error);
0167     }
0168
0169     return h;
0170 }
0171
0172 TH1 *computeEffVsBEff(TH1 *effVsCut, TH1 *effVsCutB)
0173 {
0174     TH1 *h = new TH1F(Form("%s_effVsBEff", effVsCut->GetName()), "effVsBEff",
0175                       100, 0, 1);
0176     h->SetMaximum(1.5);
0177     h->SetMinimum(1e-3);
0178
0179     unsigned int n = effVsCut->GetNbinsX();
0180     for(unsigned int bin = 1; bin <= n; bin++) {
0181         double eff = effVsCut->GetBinContent(bin);
0182         double error = effVsCut->GetBinError(bin);
0183         double effB = effVsCutB->GetBinContent(bin);
0184
0185         h->SetBinContent(h->FindBin(effB), eff);
0186         h->SetBinError(h->FindBin(effB), error);
0187         // FIXME: The error in effB is not propagated
0188     }
0189
0190     return h;
0191 }
0192
0193 TH1 *invertHisto(TH1 *h)
0194 {
0195     unsigned int n = h->GetNbinsX();
0196     TH1 *inv = new TH1F(Form("%s_inverted", h->GetName()), "inverted",
0197                         n, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
0198     for(unsigned int i = 0; i <= n + 1; i++)
0199         inv->SetBinContent(n + 1 - i, h->GetBinContent(i));
0200
0201     return inv;
0202 }
0203
0204 void setAxisStyle(TH1 *hist) {
0205     // --------------------------------------------------
0206     // define proper axsis style for a given histogram
0207     // --------------------------------------------------
0208     hist->GetXaxis()->SetTitleSize( 0.06);
0209     hist->GetXaxis()->SetTitleColor( 1);
0210     hist->GetXaxis()->SetTitleOffset( 0.8);
0211     hist->GetXaxis()->SetTitleFont( 62);
0212     hist->GetXaxis()->SetLabelSize( 0.05);
0213     hist->GetXaxis()->SetLabelFont( 62);
0214     hist->GetXaxis()->CenterTitle();
0215     hist->GetXaxis()->SetNdivisions( 505);
0216
0217     hist->GetYaxis()->SetTitleSize( 0.07);
0218     hist->GetYaxis()->SetTitleColor( 1);
0219     hist->GetYaxis()->SetTitleOffset( 0.5);
0220     hist->GetYaxis()->SetTitleFont( 62);
0221     hist->GetYaxis()->SetLabelSize( 0.05);
0222     hist->GetYaxis()->SetLabelFont( 62);
0223 }
0224
0225 void setHistStyle(TH1 *hist)
0226 {
0227     if (!hist)
0228         return;
0229
0230     // --------------------------------------------------
0231     // define proper histogram style
0232     // --------------------------------------------------
0233     setAxisStyle(hist);
0234     hist->GetXaxis()->SetTitle(hist->GetTitle());
0235     hist->SetTitle();
0236     hist->SetLineColor(4.);
0237     hist->SetLineWidth(2.);
0238     hist->SetMarkerSize(0.75);
0239     hist->SetMarkerColor(4.);
0240     hist->SetMarkerStyle(20.);
0241 }
0242
0243 void setNiceStyle()
0244 {
0245     gROOT->SetStyle("Plain");
0246
0247     // --------------------------------------------------
0248     // define proper canvas style
0249     // --------------------------------------------------
0250     TStyle *MyStyle = new TStyle ("MyStyle", "My style for nicer plots");
0251
0252     Float_t xoff = MyStyle->GetLabelOffset("X"),
0253             yoff = MyStyle->GetLabelOffset("Y"),
0254             zoff = MyStyle->GetLabelOffset("Z");
0255
0256     MyStyle->SetCanvasBorderMode ( 0 );
0257     MyStyle->SetPadBorderMode    ( 0 );
0258     MyStyle->SetPadColor         ( 0 );
0259     MyStyle->SetCanvasColor      ( 0 );
0260     MyStyle->SetTitleColor       ( 0 );
0261     MyStyle->SetStatColor        ( 0 );
0262     MyStyle->SetTitleBorderSize  ( 0 );
0263     MyStyle->SetTitleFillColor   ( 0 );
0264     MyStyle->SetTitleH        ( 0.07 );
0265     MyStyle->SetTitleW        ( 1.00 );
0266     MyStyle->SetTitleFont     (  132 );
0267
0268     MyStyle->SetLabelOffset (1.5*xoff, "X");
0269     MyStyle->SetLabelOffset (1.5*yoff, "Y");
0270     MyStyle->SetLabelOffset (1.5*zoff, "Z");
0271
0272     MyStyle->SetTitleOffset (0.9,      "X");
0273     MyStyle->SetTitleOffset (0.9,      "Y");
0274     MyStyle->SetTitleOffset (0.9,      "Z");
0275
0276     MyStyle->SetTitleSize   (0.045,    "X");
0277     MyStyle->SetTitleSize   (0.045,    "Y");
0278     MyStyle->SetTitleSize   (0.045,    "Z");
0279
0280     MyStyle->SetLabelFont   (132,      "X");
0281     MyStyle->SetLabelFont   (132,      "Y");
0282     MyStyle->SetLabelFont   (132,      "Z");
0283
0284     MyStyle->SetPalette(1);
0285
0286     MyStyle->cd();
0287 }