Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:15:42 

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 "Prototype3_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 DrawPrototype3ShowerCalib_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("libPrototype3.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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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 + "_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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 + "_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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_DrawPrototype3ShowerCalib_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 + "_DrawPrototype3ShowerCalib_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
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!