Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:35:00

0001 // call with
0002 // root -l pileupJetAnalysis_mkplots.C'("analysis.root", "pileupJetAnalyzer", 0.5)'
0003 
0004 #include <iostream>
0005 #include <stdlib.h>
0006 
0007 #include <RVersion.h>
0008 #include <TH1.h>
0009 #include <TFile.h>
0010 #include <TList.h>
0011 #include <TColor.h>
0012 #include <TGraph.h>
0013 #include <TString.h>
0014 #include <TSystem.h>
0015 #include <TStyle.h>
0016 #include <TCanvas.h>
0017 #include <TLegend.h>
0018 #include <TObject.h>
0019 #include <TObjString.h>
0020 #include <TIterator.h>
0021 #include <TControlBar.h>
0022 #include <TGWindow.h>
0023 #include <TGMenu.h>
0024 #include <TGClient.h>
0025 #include <TGResourcePool.h>
0026 #include <TVirtualX.h>
0027 #ifndef __CINT__
0028 #   include <GuiTypes.h>
0029 #endif
0030 
0031 static TFile *file = 0;
0032 static TStyle *style = 0;
0033 static const TGWindow *root = 0;
0034 
0035 static const char *base = "abs(eta) < 2.4 && et > 15";
0036 
0037 static Double_t Min(Double_t a, Double_t b)
0038 { return a < b ? a : b; }
0039 
0040 static Double_t Max(Double_t a, Double_t b)
0041 { return a > b ? a : b; }
0042 
0043 static void FindRange(TH1 *histo, Int_t &b1, Int_t &b2,
0044                       Double_t &x1, Double_t &x2)
0045 {
0046     for(Int_t i = 1; i <= histo->GetNbinsX(); i++) {
0047         if (histo->GetBinContent(i) > 0.0) {
0048             if (i < b1)
0049                 b1 = i;
0050             if (i > b2)
0051                 b2 = i;
0052 
0053             Double_t low = histo->GetBinLowEdge(i);
0054             Double_t high = histo->GetBinLowEdge(i + 1);
0055 
0056             if (low < x1)
0057                 x1 = low;
0058             if (high > x2)
0059                 x2 = high;
0060         }
0061     }
0062 }
0063 
0064 class PadService {
0065     public:
0066     PadService(TString name, TString title, Int_t nPlots);
0067     ~PadService();
0068 
0069     TVirtualPad *Next();
0070 
0071     private:
0072     TCanvas *last;
0073     Int_t   index, count;
0074 
0075     TString name, title;
0076 
0077     Int_t   width, height;
0078     Int_t   nPadsX, nPadsY;
0079 };
0080 
0081 PadService::PadService(TString name, TString title, Int_t nPlots) :
0082     last(0), index(0), count(0), name(name), title(title)
0083 {
0084     switch (nPlots) {
0085         case 1:
0086         nPadsX = 1; nPadsY = 1; width = 500; height = 1.00 * width;
0087         break;
0088         case 2:
0089         nPadsX = 2; nPadsY = 1; width = 600; height = 0.55 * width;
0090         break;
0091         case 3:
0092         nPadsX = 3; nPadsY = 1; width = 900; height = 0.40 * width;
0093         break;
0094         case 4:
0095         nPadsX = 2; nPadsY = 2; width = 600; height = 1.00 * width;
0096         break;
0097         default:
0098         nPadsX = 3; nPadsY = 2; width = 800; height = 0.55 * width;
0099         break;
0100     }
0101 }
0102 
0103 PadService::~PadService()
0104 {
0105 }
0106 
0107 TVirtualPad *PadService::Next()
0108 {
0109     if (!last) {
0110         last = new TCanvas(Form("%s_%d", (const char*)name,
0111                                 ++count),
0112                            title, count * 50 + 200, count * 20,
0113                            width, height);
0114         last->SetBorderSize(2);
0115         last->SetFrameFillColor(0);
0116         last->SetHighLightColor(0);
0117         last->Divide(nPadsX, nPadsY);
0118         index = 0;
0119         count++;
0120     }
0121 
0122     last->cd(++index);
0123     TVirtualPad *pad = last->GetPad(index);
0124 
0125     if (index == nPadsX * nPadsY)
0126         last = 0;
0127 
0128     return pad;
0129 }
0130 
0131 void pileupJetAnalysis_mkplots(TString fname, TString dir, double cut)
0132 {
0133     TFile *f = TFile::Open(fname);
0134     TString name = "discriminator";
0135 
0136     TH1 *sig = new TH1F("sig", "sig", 4000, 0, 5);
0137     TH1 *bkg = new TH1F("bkg", "bkg", 4000, 0, 5);
0138 
0139     TTree *t = dynamic_cast<TTree*>(f->Get(dir + "/jets"));
0140     TString sigCut = Form("%s && mc >= %f", base, cut);
0141     TString bkgCut = Form("%s && mc < %f", base, cut);
0142 
0143     t->Draw("min(9.999, tag) >> sig", sigCut, "goff");
0144     t->Draw("min(9.999, tag) >> bkg", bkgCut, "goff");
0145 
0146     Int_t bin1 = sig->GetNbinsX();
0147     Int_t bin2 = 1;
0148     Double_t x1 = sig->GetXaxis()->GetXmax();
0149     Double_t x2 = sig->GetXaxis()->GetXmin();
0150     FindRange(bkg, bin1, bin2, x1, x2);
0151     FindRange(sig, bin1, bin2, x1, x2);
0152 
0153     PadService pads("discriminator", "discriminator & performance", 3);
0154 
0155     TVirtualPad *pad = pads.Next();
0156 
0157     TH1 *bkg_ = (TH1*)bkg->Clone(name + "_tmpO1");
0158     bkg_->Rebin(8);
0159     double bkgN = bkg_->Integral() / bkg_->Integral("width");
0160     bkg_->SetNormFactor(bkg_->Integral() / bkg_->Integral("width"));
0161 
0162     TH1 *sig_ = (TH1*)sig->Clone(name + "_tmpO2");
0163     sig_->Rebin(8);
0164     double sigN = sig_->Integral() / sig_->Integral("width");
0165     sig_->SetNormFactor(sig_->Integral() / sig_->Integral("width"));
0166 
0167     Double_t y = Max(bkg_->GetMaximum() / bkgN, sig_->GetMaximum() / sigN);
0168     bkg_->SetStats(0);
0169     bkg_->GetXaxis()->SetRangeUser(x1, x2);
0170     bkg_->GetYaxis()->SetRangeUser(0.0, y * 1.275);
0171     bkg_->SetTitle("discriminator");
0172     bkg_->SetXTitle(name);
0173     bkg_->SetFillColor(2);
0174     bkg_->SetLineColor(2);
0175     bkg_->SetLineWidth(2);
0176     bkg_->SetFillStyle(3554);
0177     bkg_->Draw();
0178     sig_->GetXaxis()->SetRangeUser(x1, x2);
0179     sig_->SetFillColor(0);
0180     sig_->SetLineColor(4);
0181     sig_->SetLineWidth(2);
0182     sig_->SetFillStyle(1);
0183     sig_->Draw("same");
0184 
0185     TLegend *leg = new TLegend(0.6 - pad->GetRightMargin(),
0186                                1.0 - pad->GetTopMargin() - 0.15,
0187                                1.0 - pad->GetRightMargin(),
0188                                1.0 - pad->GetTopMargin());
0189     leg->SetFillStyle(1);
0190     leg->AddEntry(sig_, "Signal", "F");
0191     leg->AddEntry(bkg_, "Background", "F");
0192     leg->SetBorderSize(1);
0193     leg->SetMargin(0.3);
0194     leg->Draw("same");
0195 
0196     pad->RedrawAxis();
0197 
0198     TVirtualPad *pad = pads.Next();
0199 
0200     unsigned int n = bin2 - bin1 + 1;
0201     Double_t incr = (x2 - x1) * 0.5 / n;
0202     x1 -= incr;
0203     x2 += incr;
0204     TH1F *sum[4];
0205     for(unsigned int i = 0; i < 4; i++) {
0206         sum[i] = new TH1F(Form("%s_tmpO%d", (const char*)name, 3 + i),
0207                           "cut efficiencies",
0208                           n + 1, x1 - incr, x2 + incr);
0209         sum[i]->SetBit(kCanDelete);
0210     }
0211 
0212     Double_t sumBkg = 0.0, sumSig = 0.0;
0213     for(unsigned int i = 0; i < n; i++) {
0214         sum[0]->SetBinContent(i + 1, sumBkg);
0215         sum[1]->SetBinContent(i + 1, sumSig);
0216         sumBkg += bkg->GetBinContent(i + bin1);
0217         sumSig += sig->GetBinContent(i + bin1);
0218     }
0219     sum[0]->SetBinContent(n + 1, sumBkg);
0220     sum[1]->SetBinContent(n + 1, sumSig);
0221     for(unsigned int i = 0; i <= n; i++) {
0222         sum[2]->SetBinContent(i + 1, sumBkg);
0223         sum[3]->SetBinContent(i + 1, sumSig);
0224     }
0225     for(unsigned int i = 0; i <= n; i++) {
0226         sum[0]->SetBinContent(
0227             i + 1, sumBkg - sum[0]->GetBinContent(i + 1));
0228         sum[1]->SetBinContent(
0229             i + 1, sumSig - sum[1]->GetBinContent(i + 1));
0230     }
0231 
0232     for(unsigned int i = 0; i < 2; i++)
0233         sum[i]->Divide(sum[i], sum[2 + i], 1, 1, "b");
0234 
0235     sum[2]->Delete();
0236     sum[3]->Delete();
0237 
0238     pad->SetGridy();
0239     pad->SetTicky();
0240     sum[0]->SetLineColor(2);
0241     sum[0]->SetLineWidth(2);
0242     sum[0]->GetYaxis()->SetRangeUser(0.0, 1.005);
0243     sum[0]->SetXTitle(name);
0244     sum[0]->GetYaxis()->SetTitleOffset(1.25);
0245     sum[0]->SetYTitle("efficiency");
0246     sum[0]->SetStats(0);
0247     sum[0]->Draw("C");
0248     sum[1]->SetLineColor(4);
0249     sum[1]->SetLineWidth(2);
0250     sum[1]->Draw("C same");
0251 
0252     pad->RedrawAxis();
0253 
0254     gStyle->SetFillColor(kWhite);
0255     pad->SetFillColor(kWhite);
0256     TVirtualPad *pad = pads.Next();
0257 
0258     TH1F *tmp = new TH1F(name + "_tmpO6",
0259                          "efficiency vs. purity", n, 0, 1);
0260     tmp->SetBit(kCanDelete);
0261 
0262     pad->SetGridx();
0263     pad->SetGridy();
0264     pad->SetTickx();
0265     pad->SetTicky();
0266     tmp->GetXaxis()->SetRangeUser(0.0, 1.005);
0267     tmp->GetYaxis()->SetRangeUser(0.0, 1.005);
0268     tmp->SetXTitle("efficiency");
0269     tmp->SetYTitle("purity");
0270     tmp->SetStats(0);
0271     tmp->SetLineWidth(2);
0272     tmp->SetFillColor(kWhite);
0273     tmp->Draw();
0274 
0275     Float_t *values[2];
0276     values[0] = new Float_t[n + 1];
0277     values[1] = new Float_t[n + 1];
0278     for(unsigned int i = 0; i <= n; i++) {
0279         values[0][i] = sum[1]->GetBinContent(i + 1);
0280         values[1][i] = 1.0 - sum[0]->GetBinContent(i + 1);
0281     }
0282     TGraph *graph = new TGraph(n + 1, values[0], values[1]);
0283     delete[] values[0];
0284     delete[] values[1];
0285 
0286     graph->SetName(name + "_tmpO7");
0287     graph->SetBit(kCanDelete);
0288     graph->SetLineColor(4);
0289     graph->SetLineWidth(2);
0290     graph->Draw("C");
0291 
0292     pad->RedrawAxis();
0293 }