File indexing completed on 2024-04-06 12:25:29
0001
0002
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 }
