EMCal/macros/DrawPrototype2ShowerCalib_Sum.C

0001 // $Id: $
0002
0003 /*!
0004  * \file Draw.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include <cmath>
0012 #include <TFile.h>
0013 #include <TString.h>
0014 #include <TLine.h>
0015 #include <TTree.h>
0016 #include <TLatex.h>
0017 #include <TGraphErrors.h>
0018 #include <cassert>
0019 #include "SaveCanvas.C"
0020 #include "SetOKStyle.C"
0021 using namespace std;
0022
0023 //#include "Prototype2_DSTReader.h"
0024
0025 TCut event_sel;
0026 TString cuts;
0027 TFile * _file0 = NULL;
0028 TTree * T = NULL;
0029
0030 class lin_res
0031 {
0032 public:
0033   TString name;
0034   TGraphErrors * linearity;
0035   TGraphErrors * resolution;
0036   TF1 * f_res;
0037 };
0038
0039 void
0040 DrawPrototype2ShowerCalib_Sum( //
0041 )
0042 {
0043
0044   SetOKStyle();
0045   gStyle->SetOptStat(0);
0046   gStyle->SetOptFit(1111);
0047   gStyle->SetPadGridX(0);
0048   gStyle->SetPadGridY(0);
0049   TVirtualFitter::SetDefaultFitter("Minuit2");
0050
0051   gSystem->Load("libPrototype2.so");
0052   gSystem->Load("libProto2ShowCalib.so");
0053
0054   RejectionCompare(4);
0055   RejectionCompare(8);
0056   RejectionCompare(12);
0057
0058 //  LineShapeCompare(4, "QGSP_BERT_HP");
0059 //  LineShapeCompare(8, "QGSP_BERT_HP");
0060 ////  LineShapeCompare(12, "QGSP_BERT_HP");
0061 //  LineShapeCompare(8, "QGSP_BERT");
0062 //  LineShapeCompare(8, "FTFP_BERT_HP");
0063 //  LineShapeCompare(8, "FTFP_BERT");
0064
0065 //  LineShapeCompare_Electron(4, "QGSP_BERT_HP");
0066 //  LineShapeCompare_Electron(8, "QGSP_BERT_HP");
0067 //  LineShapeCompare_Electron(12, "QGSP_BERT_HP");
0068
0069 //  LineShapeCompare(8, "BirkConst0.151");
0070 //  LineShapeCompare(8, "BirkConst0.18");
0071 //    LineShapeCompare_Electron(8, "BirkConst0.151");
0072
0073 }
0074
0075 double
0076 GetMu2Pi(const double E)
0077 {
0078   double mu2pi = 0;
0079   if (E == 4)
0080     {
0081       mu2pi = 0.02 / 0.105;
0082     }
0083
0084   else if (E == 8)
0085     {
0086       mu2pi = 0.02 / 0.28;
0087     }
0088
0089   else if (E == 12)
0090     {
0091       mu2pi = 0.04 / 0.48;
0092     }
0093   else
0094     {
0095       cout << "GetMu2Pi - missing mu2pi data!!" << endl;
0096       assert(0);
0097     }
0098
0099   return mu2pi;
0100 }
0101
0102 void
0103 RejectionCompare(const double E = 4)
0104 {
0105
0106   double mu2pi = GetMu2Pi(E);
0107
0108   TString energy(Form("%.0f", E));
0109
0110   TString data_file(
0111       "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeData_Neg"
0112           + energy + "GeV_quality_h12345_v34567_col2_row2.root");
0113
0114   TFile * fdata = new TFile(data_file);
0115   assert(fdata->IsOpen());
0116
0117   TFile * fmu =
0118       new TFile(
0119           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll/Prototype_mu-_"
0120               + energy
0121               + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h5_v5_col2_row2.root");
0122   assert(fmu->IsOpen());
0123   cout << "Processing " << fmu->GetName() << endl;
0124
0125   TH1F * h_5x5sum_c2_sum_mu = fmu->Get("h_5x5sum_c2_sum")->Clone(
0126       "h_5x5sum_c2_sum_mu");
0127
0128   int idx = 1;
0129   TPad * p;
0130
0131   TH1F * h_5x5sum_c2_h = fdata->Get("h_5x5sum_c2_h");
0132
0133   TH1F * h_5x5sum_c2_sum_data_sub_mu = h_5x5sum_c2_h->Clone(
0134       "h_5x5sum_c2_sum_data_sub_mu");
0135
0136   cout << "h_5x5sum_c2_sum_data->GetNbinsX() = "
0137       << h_5x5sum_c2_sum_data_sub_mu->GetNbinsX() << endl;
0138   cout << "h_5x5sum_c2_sum_mu->GetNbinsX() = "
0139       << h_5x5sum_c2_sum_mu->GetNbinsX() << endl;
0140
0141   assert(
0142       h_5x5sum_c2_sum_data_sub_mu->GetNbinsX()
0143           == h_5x5sum_c2_sum_mu->GetNbinsX());
0144
0145   const double muon_count = h_5x5sum_c2_sum_data_sub_mu->GetSum()
0146       * (mu2pi / (1 + mu2pi));
0147   h_5x5sum_c2_sum_data_sub_mu->Add(h_5x5sum_c2_sum_mu,
0148       -muon_count / h_5x5sum_c2_sum_mu->GetSum());
0149
0150   TGraphErrors * ge_data = Distribution2Efficiency(h_5x5sum_c2_sum_data_sub_mu);
0151   ge_data->SetFillColor(kGray);
0152   ge_data->SetLineWidth(3);
0153
0154   TGraphErrors * ge_pi_QGSP_BERT_HP = GetSimRejCurve("QGSP_BERT_HP", energy,
0155       "pi-");
0156   TGraphErrors * ge_kaon_QGSP_BERT_HP = GetSimRejCurve("QGSP_BERT_HP", energy,
0157       "kaon-");
0158
0159   ge_pi_QGSP_BERT_HP->SetLineColor(kBlue);
0160   ge_kaon_QGSP_BERT_HP->SetLineColor(kCyan + 3);
0161
0162   TGraphErrors * ge_pi_FTFP_BERT_HP = GetSimRejCurve("FTFP_BERT_HP", energy,
0163       "pi-");
0164   TGraphErrors * ge_kaon_FTFP_BERT_HP = GetSimRejCurve("FTFP_BERT_HP", energy,
0165       "kaon-");
0166
0167   ge_pi_FTFP_BERT_HP->SetLineColor(kBlue);
0168   ge_kaon_FTFP_BERT_HP->SetLineColor(kCyan + 3);
0169   ge_pi_FTFP_BERT_HP->SetLineStyle(kDashed);
0170   ge_kaon_FTFP_BERT_HP->SetLineStyle(kDashed);
0171
0172   TGraphErrors * ge_pi_BirkConst18 = GetSimRejCurve("BirkConst0.18", energy,
0173       "pi-");
0174   TGraphErrors * ge_kaon_BirkConst18 = GetSimRejCurve("BirkConst0.18", energy,
0175       "kaon-");
0176
0177   ge_pi_BirkConst18->SetLineColor(kBlue);
0178   ge_kaon_BirkConst18->SetLineColor(kCyan + 3);
0179   ge_pi_BirkConst18->SetLineStyle(9);
0180   ge_kaon_BirkConst18->SetLineStyle(9);
0181
0182 //  new TCanvas;
0183 //  h_5x5sum_c2_sum_data_sub_mu->Draw();
0184 //  h_5x5sum_c2_sum_mu->Draw("same");
0185
0186   TGraphErrors * ge_ref = ge_pi_QGSP_BERT_HP->Clone("ge_ref");
0187
0188   for (int i = 0; i < ge_ref->GetN(); ++i)
0189     {
0190       (ge_ref->GetY())[i] = ((ge_pi_QGSP_BERT_HP->GetY())[i]
0191           + (ge_kaon_QGSP_BERT_HP->GetY())[i]) / 2.;
0192     }
0193
0194   TText * t;
0195   TCanvas *c1 = new TCanvas("RejectionCompare_" + energy + "GeV_",
0196       "RejectionCompare_" + energy + "GeV_", 800, 1000);
0197
0198   c1->Divide(1, 2);
0199
0200   p = (TPad *) c1->cd(idx++);
0201   c1->Update();
0202   p->SetLogy();
0203
0204   p->SetBottomMargin(0.03);
0205
0206   TH1 * hframe = p->DrawFrame(0, 5e-4, E * 1, 1);
0207   hframe->SetTitle(
0208       Form(
0209           ";;1/(Hadron Rejection)",
0210           E));
0211   hframe->GetXaxis()->SetLabelSize(.06);
0212   hframe->GetXaxis()->SetTitleSize(.06);
0213   hframe->GetYaxis()->SetLabelSize(.06);
0214   hframe->GetYaxis()->SetTitleSize(.06);
0215   hframe->GetYaxis()->SetTitleOffset(0.8);
0216
0217   ge_data->Draw("p3");
0218   ge_data->Draw("lX");
0219
0220   ge_pi_QGSP_BERT_HP->Draw("lX");
0221   ge_kaon_QGSP_BERT_HP->Draw("lX");
0222   ge_pi_FTFP_BERT_HP->Draw("lX");
0223   ge_kaon_FTFP_BERT_HP->Draw("lX");
0224   ge_pi_BirkConst18->Draw("lX");
0225   ge_kaon_BirkConst18->Draw("lX");
0226
0227   TLegend * leg = new TLegend(.12, .05, .72, .5);
0228   leg->SetHeader(Form("Beam momentum = -%.0f GeV/c",E));
0229   leg->AddEntry(ge_data, "T-1044 non-e^{-}  Data - Muon Sim.", "fl");
0230   leg->AddEntry(ge_pi_QGSP_BERT_HP, "\pi^{-}, QGSP_BERT_HP (Default)", "l");
0231   leg->AddEntry(ge_kaon_QGSP_BERT_HP, "K^{-}, QGSP_BERT_HP (Default)", "l");
0232   leg->AddEntry(ge_pi_FTFP_BERT_HP, "\pi^{-}, FTFP_BERT_HP", "l");
0233   leg->AddEntry(ge_kaon_FTFP_BERT_HP, "K^{-}, FTFP_BERT_HP", "l");
0234   leg->AddEntry(ge_pi_BirkConst18, "\pi^{-}, k_{B} = 0.18 mm/MeV", "l");
0235   leg->AddEntry(ge_kaon_BirkConst18, "K^{-}, k_{B} = 0.18 mm/MeV", "l");
0236
0237   leg->Draw();
0238
0239   p = (TPad *) c1->cd(idx++);
0240   c1->Update();
0241 //  p->SetLogy();
0242   p->SetTopMargin(0);
0243
0244   TH1 * hframe =
0245   p->DrawFrame(0, 0, E * 1, 2);
0246   hframe->SetTitle(
0247       Form(
0248           ";Minimal cut on 5x5 Cluster Energy (GeV/c);Ratio of 1/(Hadron Rejection) and Ref.",
0249           E));
0250   hframe->GetXaxis()->SetLabelSize(.06);
0251   hframe->GetXaxis()->SetTitleSize(.06);
0252   hframe->GetYaxis()->SetLabelSize(.06);
0253   hframe->GetYaxis()->SetTitleSize(.06);
0254   hframe->GetYaxis()->SetTitleOffset(0.8);
0255
0256   TGraphErrors * ge_data_ratio = GetTGraphRatio(ge_data, ge_ref);
0257
0258   ge_data_ratio->Draw("p3");
0259   ge_data_ratio->Draw("lX");
0260
0261   GetTGraphRatio(ge_pi_QGSP_BERT_HP, ge_ref)->Draw("lX");
0262   GetTGraphRatio(ge_kaon_QGSP_BERT_HP, ge_ref)->Draw("lX");
0263   GetTGraphRatio(ge_pi_FTFP_BERT_HP, ge_ref)->Draw("lX");
0264   GetTGraphRatio(ge_kaon_FTFP_BERT_HP, ge_ref)->Draw("lX");
0265   GetTGraphRatio(ge_pi_BirkConst18, ge_ref)->Draw("lX");
0266   GetTGraphRatio(ge_kaon_BirkConst18, ge_ref)->Draw("lX");
0267
0268   SaveCanvas(c1,
0269       data_file + "_DrawPrototype2ShowerCalib_Sum" + TString(c1->GetName()),
0270       kTRUE);
0271 }
0272
0273 TGraphErrors *
0274 GetSimRejCurve(const TString physics_lst, const TString energy,
0275     const TString particle)
0276 {
0277   const TString filename =
0278       "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_"
0279           + physics_lst + "/Prototype_" + particle + "_" + energy
0280           + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h5_v5_col2_row2.root";
0281   cout << "GetSimRejCurve - processing " << filename << endl;
0282
0283   TFile * f = new TFile(filename);
0284   assert(f->IsOpen());
0285
0286   TGraphErrors * ge = Distribution2Efficiency(
0287       (TH1F *) f->Get("h_5x5sum_c2_sum"));
0288 //  ge->SetLineColor(kBlue);
0289   ge->SetLineWidth(3);
0290
0291   return ge;
0292 }
0293
0294 void
0295 LineShapeCompare(const double E = 4, //
0296     const TString physics_lst = "QGSP_BERT_HP")
0297 {
0298
0299   double mu2pi = GetMu2Pi(E);
0300
0301   TString energy(Form("%.0f", E));
0302
0303   TText * t;
0304   TCanvas *c1 = new TCanvas("LineShapeCompare_" + energy + "GeV_" + physics_lst,
0305       "LineShapeCompare_" + energy + "GeV_" + physics_lst, 1000, 650);
0306
0307   int idx = 1;
0308   TPad * p;
0309
0310   p = (TPad *) c1->cd(idx++);
0311   c1->Update();
0312   p->SetLogy();
0313
0314   TString data_file(
0315       "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeData_Neg"
0316           + energy + "GeV_quality_h12345_v34567_col2_row2.root");
0317
0318   TFile * fdata = new TFile(data_file);
0319   assert(fdata->IsOpen());
0320
0321   TH1F * h_5x5sum_c2_sum_data = fdata->Get("h_5x5sum_c2_h3")->Clone(
0322       "h_5x5sum_c2_sum_data");
0323   assert(h_5x5sum_c2_sum_data);
0324
0325   h_5x5sum_c2_sum_data->Scale(1. / h_5x5sum_c2_sum_data->GetSum());
0326   h_5x5sum_c2_sum_data->SetLineColor(kBlack);
0327 //  h_5x5sum_c2_sum_data->SetLineStyle(kDashed);
0328
0329   h_5x5sum_c2_sum_data->Draw();
0330
0331   TFile * fpi =
0332       new TFile(
0333           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_"
0334               + physics_lst + "/Prototype_pi-_" + energy
0335               + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h5_v5_col2_row2.root");
0336   assert(fpi->IsOpen());
0337
0338   TH1F * h_5x5sum_c2_sum_pi = fpi->Get("h_5x5sum_c2_sum")->Clone(
0339       "h_5x5sum_c2_sum_pi");
0340   assert(h_5x5sum_c2_sum_pi);
0341
0342   h_5x5sum_c2_sum_pi->Scale((1. / (1 + mu2pi)) / h_5x5sum_c2_sum_pi->GetSum());
0343   h_5x5sum_c2_sum_pi->SetLineColor(kBlue);
0344 //  h_5x5sum_c2_sum_pi->SetLineStyle(kDashed);
0345
0346   h_5x5sum_c2_sum_pi->Draw("same");
0347
0348   TFile * fkaon =
0349       new TFile(
0350           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_"
0351               + physics_lst + "/Prototype_kaon-_" + energy
0352               + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h5_v5_col2_row2.root");
0353   assert(fkaon->IsOpen());
0354
0355   TH1F * h_5x5sum_c2_sum_kaon = fkaon->Get("h_5x5sum_c2_sum")->Clone(
0356       "h_5x5sum_c2_sum_kaon");
0357   assert(h_5x5sum_c2_sum_kaon);
0358
0359   h_5x5sum_c2_sum_kaon->Scale(
0360       (1. / (1 + mu2pi)) / h_5x5sum_c2_sum_kaon->GetSum());
0361   h_5x5sum_c2_sum_kaon->SetLineColor(kCyan + 3);
0362 //  h_5x5sum_c2_sum_pi->SetLineStyle(kDashed);
0363
0364   h_5x5sum_c2_sum_kaon->Draw("same");
0365
0366   TFile * fmu =
0367       new TFile(
0368           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_"
0369               + physics_lst + "/Prototype_mu-_" + energy
0370               + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h5_v5_col2_row2.root");
0371   assert(fmu->IsOpen());
0372   cout << "Processing " << fmu->GetName() << endl;
0373
0374   TH1F * h_5x5sum_c2_sum_mu = fmu->Get("h_5x5sum_c2_sum")->Clone(
0375       "h_5x5sum_c2_sum_mu");
0376   assert(h_5x5sum_c2_sum_mu);
0377
0378   h_5x5sum_c2_sum_mu->Scale(
0379       (mu2pi / (1 + mu2pi)) / h_5x5sum_c2_sum_mu->GetSum());
0380   h_5x5sum_c2_sum_mu->SetLineColor(kBlack);
0381 //  h_5x5sum_c2_sum_pi->SetLineStyle(kDashed);
0382
0383   h_5x5sum_c2_sum_mu->Draw("same");
0384
0385   c1->Update();
0386
0387   SaveCanvas(c1,
0388       data_file + "_DrawPrototype2ShowerCalib_Sum" + TString(c1->GetName()),
0389       kTRUE);
0390
0391 }
0392
0393 void
0394 LineShapeCompare_Electron(const double E = 4, //
0395     const TString physics_lst = "QGSP_BERT_HP")
0396 {
0397
0398   TString energy(Form("%.0f", E));
0399
0400   TText * t;
0401   TCanvas *c1 = new TCanvas(
0402       "LineShapeCompare_Electron_" + energy + "GeV_" + physics_lst,
0403       "LineShapeCompare_Electron_" + energy + "GeV_" + physics_lst, 1000, 650);
0404
0405   int idx = 1;
0406   TPad * p;
0407
0408   p = (TPad *) c1->cd(idx++);
0409   c1->Update();
0410   p->SetLogy();
0411
0412   TString data_file(
0413       "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeData_Neg"
0414           + energy + "GeV_quality_h3_v5_col2_row2.root");
0415
0416   TFile * fdata = new TFile(data_file);
0417   assert(fdata->IsOpen());
0418
0419   TH1F * h_5x5sum_c2_sum_data = fdata->Get("h_5x5sum_c2_e");
0420   assert(h_5x5sum_c2_sum_data);
0421
0422   h_5x5sum_c2_sum_data->Scale(1. / h_5x5sum_c2_sum_data->GetSum());
0423   h_5x5sum_c2_sum_data->SetLineColor(kBlack);
0424 //  h_5x5sum_c2_sum_data->SetLineStyle(kDashed);
0425
0426   h_5x5sum_c2_sum_data->Draw();
0427
0428   TFile * fe =
0429       new TFile(
0430           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_"
0431               + physics_lst + "/Prototype_e-_" + energy
0432               + "_SegALL_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeSim_0DegreeRot_h1_v1_col2_row2.root");
0433   assert(fe->IsOpen());
0434
0435   TH1F * h_5x5sum_c2_sum_e = fe->Get("h_5x5sum_c2_sum");
0436   assert(h_5x5sum_c2_sum_e);
0437
0438   h_5x5sum_c2_sum_e->Scale((0.27 / 0.30) / h_5x5sum_c2_sum_e->GetSum());
0439   h_5x5sum_c2_sum_e->SetLineColor(kBlue + 3);
0440 //  h_5x5sum_c2_sum_e->SetLineStyle(kDashed);
0441
0442   h_5x5sum_c2_sum_e->Draw("same");
0443
0444   c1->Update();
0445
0446   SaveCanvas(c1,
0447       data_file + "_DrawPrototype2ShowerCalib_Sum" + TString(c1->GetName()),
0448       kTRUE);
0449 }
0450
0451 TGraphErrors *
0452 FitProfile(const TH2 * h2)
0453 {
0454
0455   TProfile * p2 = h2->ProfileX();
0456
0457   int n = 0;
0458   double x[1000];
0459   double ex[1000];
0460   double y[1000];
0461   double ey[1000];
0462
0463   for (int i = 1; i <= h2->GetNbinsX(); i++)
0464     {
0465       TH1D * h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0466
0467       if (h1->GetSum() < 30)
0468         {
0469           cout << "FitProfile - ignore bin " << i << endl;
0470           continue;
0471         }
0472       else
0473         {
0474           cout << "FitProfile - fit bin " << i << endl;
0475         }
0476
0477       TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0478           p2->GetBinError(i) * 4);
0479
0480       TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0481           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0482
0483       fgaus.SetParameter(1, p2->GetBinContent(i));
0484       fgaus.SetParameter(2, p2->GetBinError(i));
0485
0486       h1->Fit(&fgaus, "MQ");
0487
0488       f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
0489           fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
0490           fgaus.GetParameter(2) / 4, 0);
0491
0492       h1->Fit(&f2, "MQ");
0493
0494 //      new TCanvas;
0495 //      h1->Draw();
0496 //      fgaus.Draw("same");
0497 //      break;
0498
0499       x[n] = p2->GetBinCenter(i);
0500       ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0501       y[n] = fgaus.GetParameter(1);
0502       ey[n] = fgaus.GetParError(1);
0503
0504 //      p2->SetBinContent(i, fgaus.GetParameter(1));
0505 //      p2->SetBinError(i, fgaus.GetParameter(2));
0506
0507       n++;
0508       delete h1;
0509     }
0510
0511   return new TGraphErrors(n, x, y, ex, ey);
0512 }
0513
0514 TGraphErrors *
0515 Distribution2Efficiency(TH1F * hCEMC3_Max)
0516 {
0517   double threshold[10000] =
0518     { 0 };
0519   double eff[10000] =
0520     { 0 };
0521   double eff_err[10000] =
0522     { 0 };
0523
0524   assert(hCEMC3_Max->GetNbinsX() < 10000);
0525
0526   const double n = hCEMC3_Max->GetSum();
0527   double pass = 0;
0528   int cnt = 0;
0529   for (int i = hCEMC3_Max->GetNbinsX(); i >= 1; i--)
0530     {
0531       pass += hCEMC3_Max->GetBinContent(i);
0532
0533       const double pp = pass / n;
0534 //      const double z = 1.96;
0535       const double z = 1.;
0536
0537       const double A = z * sqrt(1. / n * pp * (1 - pp) + z * z / 4 / n / n);
0538       const double B = 1 / (1 + z * z / n);
0539
0540       threshold[cnt] = hCEMC3_Max->GetBinCenter(i);
0541       eff[cnt] = (pp + z * z / 2 / n) * B;
0542       eff_err[cnt] = A * B;
0543
0544 //      cout << threshold[cnt] << ": " << "CL " << eff[cnt] << "+/-"
0545 //          << eff_err[cnt] << endl;
0546       cnt++;
0547     }
0548   TGraphErrors * ge = new TGraphErrors(cnt, threshold, eff, NULL, eff_err);
0549   ge->SetName(TString("ge_") + hCEMC3_Max->GetName());
0550   return ge;
0551 }
0552
0553 TGraphErrors *
0554 GetTGraphRatio(TGraphErrors *input, TGraphErrors * baseline)
0555 {
0556   assert(input);
0557   assert(baseline);
0558
0559   assert(input->GetN() == baseline->GetN());
0560
0561   TGraphErrors * ge = input->Clone(TString(input->GetName()) + "_Ratio");
0562
0563   for (int i = 0; i < input->GetN(); ++i)
0564     {
0565       (ge->GetY())[i] /= (baseline->GetY())[i];
0566       (ge->GetEY())[i] /= (baseline->GetY())[i];
0567     }
0568
0569   return ge;
0570 }
0571