EMCal/macros/DrawEMCalTower_Resolution.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 void
0024 DrawEMCalTower_Resolution( //
0025     const TString base =
0026         "../../sPHENIX_work/Prototype_2016/EMCal_sim/15Degree_1Col_LightCollection/")
0027 {
0028   SetOKStyle();
0029   gStyle->SetOptStat(0);
0030   gStyle->SetOptFit(1111);
0031   TVirtualFitter::SetDefaultFitter("Minuit2");
0032
0033   const int N = 7;
0034   const double Es[] =
0035     { 1, 2, 4, 8, 12, 16, 24 };
0036   double mean[1000];
0037   double dmean[1000];
0038   double res[1000];
0039   double dres[1000];
0040
0041   for (int i = 0; i < N; ++i)
0042     {
0043       const double E = Es[i];
0044
0045       vector<double> v(4);
0046
0047       const TString sE = Form("%.0f", E);
0048       v = GetOnePlot(base, "e-", sE);
0049
0050       mean[i] = v[0];
0051       dmean[i] = v[1];
0052       res[i] = v[2] / v[0];
0053       dres[i] = v[3] / v[0];
0054
0055       cout << "mean[i] = " << mean[i] << ", " //
0056           << "dmean[i] = " << dmean[i] << ", " //
0057           << "res[i] = " << res[i] << ", " //
0058           << "dres[i] = " << dres[i] << endl;
0059     }
0060
0061   TGraphErrors * ge_linear = new TGraphErrors(N, Es, mean, 0, dmean);
0062   TGraphErrors * ge_res = new TGraphErrors(N, Es, res, 0, dres);
0063
0064   ge_linear->SetLineColor(kBlue + 3);
0065   ge_linear->SetMarkerColor(kBlue + 3);
0066   ge_linear->SetLineWidth(2);
0067   ge_linear->SetMarkerStyle(kFullCircle);
0068   ge_linear->SetMarkerSize(1.5);
0069
0070   ge_res->SetLineColor(kBlue + 3);
0071   ge_res->SetMarkerColor(kBlue + 3);
0072   ge_res->SetLineWidth(2);
0073   ge_res->SetMarkerStyle(kFullCircle);
0074   ge_res->SetMarkerSize(1.5);
0075
0076   TF1 * f_calo_r = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 60);
0077   TF1 * f_calo_l = new TF1("f_calo_l", "pol2", 0.5, 60);
0078
0079   TText * t;
0080   TCanvas *c1 = new TCanvas("DrawEMCalTower_Resolution",
0081       "DrawEMCalTower_Resolution", 1300, 600);
0082   c1->Divide(2, 1);
0083   int idx = 1;
0084   TPad * p;
0085
0086   p = (TPad *) c1->cd(idx++);
0087   c1->Update();
0088 //  p->SetLogy();
0089   p->SetGridx(0);
0090   p->SetGridy(0);
0091
0092   p->DrawFrame(0, 0, 30, 30,
0093       Form("Linearity;Input energy (GeV);Measured Energy (GeV)", E));
0094   TLine * l = new TLine(0, 0, 30, 30);
0095   l->SetLineColor(kGray);
0096   l->Draw();
0097
0098   ge_linear->Draw("p");
0099   ge_linear->Fit(f_calo_l, "RM0");
0100   f_calo_l->Draw("same");
0101
0102   p = (TPad *) c1->cd(idx++);
0103   c1->Update();
0104 //  p->SetLogy();
0105   p->SetGridx(0);
0106   p->SetGridy(0);
0107
0108   TH1 * hframe = p->DrawFrame(0, 0, 30, 0.2,
0109       Form("Resolution;Input energy (GeV);#DeltaE/<E>", E));
0110
0111   ge_res->Draw("p");
0112   ge_res->Fit(f_calo_r, "RM0");
0113   f_calo_r->Draw("same");
0114
0115   hframe->SetTitle(
0116       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0117           f_calo_r->GetParameter(0), f_calo_r->GetParameter(1)));
0118
0119   SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
0120 }
0121
0122 vector<double>
0123 GetOnePlot(const TString base, const TString particle, const TString sE)
0124 {
0125
0126   TString fname =
0127       base + TString("Prototype_") + particle + "_" + sE
0128           + "_SegALL_DSTReader.root_DrawEMCalTower_EMCDistribution_SUM_all_event.root";
0129 //  Prototype_e-_2_Seg0_DSTReader.root_DrawEMCalTower_EMCDistribution_SUMall_event.root
0130
0131   cout << "Process " << fname << endl;
0132   TFile * f = new TFile(fname);
0133   assert(f->IsOpen());
0134
0135   TH1 * hEnergySum = (TH1 *) f->GetObjectChecked("EnergySum_LG", "TH1");
0136   assert(hEnergySum);
0137   new TCanvas();
0138   hEnergySum->DrawClone();
0139
0140   hEnergySum->Scale(1. / hEnergySum->Integral(1, -1));
0141   TF1 * fgaus = hEnergySum->GetFunction("fgaus_LG");
0142   assert(fgaus);
0143
0144   vector<double> v(4);
0145   v[0] = fgaus->GetParameter(1);
0146   v[1] = fgaus->GetParError(1);
0147   v[2] = fgaus->GetParameter(2);
0148   v[3] = fgaus->GetParError(2);
0149
0150   return v;
0151 }
0152
0153 void
0154 DrawEMCalTower_Resolution_Compare()
0155 {
0156   gStyle->SetOptStat(0);
0157
0158   TFile * f1 =
0159       new TFile(
0160           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/15Degree_1Col/DrawEMCalTower_Resolution.root");
0161
0162   TGraphErrors * g1 = f1->GetObjectChecked("Graph", "TGraphErrors");
0163   TF1 * f_calo_r1 = f1->GetObjectChecked("f_calo_r", "TF1");
0164
0165   g1->SetLineColor(kBlue + 3);
0166   g1->SetMarkerColor(kBlue + 3);
0167   f_calo_r1->SetLineColor(kBlue + 3);
0168 //  g1->GetHistogram()->SetStats(0);
0169   g1->FindObject("stats")->Delete();
0170
0171   TFile * f2 =
0172       new TFile(
0173           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/15Degree_1Col_LightCollection/DrawEMCalTower_Resolution.root");
0174
0175   TGraphErrors * g2 = f2->GetObjectChecked("Graph", "TGraphErrors");
0176   TF1 * f_calo_r2 = f2->GetObjectChecked("f_calo_r", "TF1");
0177
0178   g2->SetLineColor(kRed + 3);
0179   g2->SetMarkerColor(kRed + 3);
0180   f_calo_r2->SetLineColor(kRed + 3);
0181 //  g2->GetHistogram()->SetStats(false);
0182   g2->FindObject("stats")->Delete();
0183
0184   TCanvas *c1 = new TCanvas("DrawEMCalTower_Resolution_Compare",
0185       "DrawEMCalTower_Resolution_Compare", 700, 600);
0186   c1->Divide(1, 1);
0187   int idx = 1;
0188   TPad * p;
0189
0190   TLegend * leg = new TLegend(0.3, 0.5, 0.9, 0.9);
0191
0192   p = (TPad *) c1->cd(idx++);
0193   c1->Update();
0194 //  p->SetLogy();
0195   p->SetGridx(0);
0196   p->SetGridy(0);
0197
0198   TH1 * hframe = p->DrawFrame(0, 0, 30, 0.2,
0199       Form("Resolution Comparison;Input energy (GeV);#DeltaE/<E>"));
0200
0201   g1->Draw("ep");
0202   f_calo_r1->Draw("same");
0203
0204   leg->AddEntry(g1, "Pre-test beam simulation", "p");
0205   leg->AddEntry(f_calo_r1,
0206       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0207           f_calo_r1->GetParameter(0), f_calo_r1->GetParameter(1)), "l");
0208
0209   g2->Draw("ep");
0210   f_calo_r2->Draw("same");
0211
0212   leg->AddEntry(g2, "Pre-test beam simulation", "p");
0213   leg->AddEntry(f_calo_r2,
0214       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0215           f_calo_r2->GetParameter(0), f_calo_r2->GetParameter(1)), "l");
0216
0217   leg->Draw();
0218
0219   SaveCanvas(c1, TString(c1->GetName()), kTRUE);
0220 }
0221