EMCal/macros/DrawPrototype2EMCalTower_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 DrawPrototype2EMCalTower_Resolution(
0025     //
0026     const TString base =
0027         "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/Prototype_2016/Production_0417_CEMC_MIP_set2_v3/",
0028     TString cut = "col1_row2_5x5_Valid_HODO_center_col1_row2")
0029 {
0030   SetOKStyle();
0031   gStyle->SetOptStat(0);
0032   gStyle->SetOptFit(0);
0033   TVirtualFitter::SetDefaultFitter("Minuit2");
0034
0035   const int N = 6;
0036   const double Es[] =
0037     { 2, 3, 4, 8, 12, 16 };
0038   const double runs[] =
0039     { 2042, 2040, 2039, 2038, 2067, 2063 };
0040
0041   double mean[1000];
0042   double dmean[1000];
0043   double res[1000];
0044   double dres[1000];
0045
0046   for (int i = 0; i < N; ++i)
0047     {
0048       const double E = Es[i];
0049
0050       vector<double> v(4);
0051
0052       v = GetOnePlot(base, runs[i], cut);
0053
0054       mean[i] = v[0];
0055       dmean[i] = v[1];
0056       res[i] = v[2] / v[0];
0057       dres[i] = v[3] / v[0];
0058
0059       cout << "mean[i] = " << mean[i] << ", " //
0060           << "dmean[i] = " << dmean[i] << ", " //
0061           << "res[i] = " << res[i] << ", " //
0062           << "dres[i] = " << dres[i] << endl;
0063     }
0064
0065   TGraphErrors * ge_linear = new TGraphErrors(N, Es, mean, 0, dmean);
0066   TGraphErrors * ge_res = new TGraphErrors(N, Es, res, 0, dres);
0067
0068   ge_linear->SetLineColor(kBlue + 3);
0069   ge_linear->SetMarkerColor(kBlue + 3);
0070   ge_linear->SetLineWidth(2);
0071   ge_linear->SetMarkerStyle(kFullCircle);
0072   ge_linear->SetMarkerSize(1.5);
0073
0074   ge_res->SetLineColor(kBlue + 3);
0075   ge_res->SetMarkerColor(kBlue + 3);
0076   ge_res->SetLineWidth(2);
0077   ge_res->SetMarkerStyle(kFullCircle);
0078   ge_res->SetMarkerSize(1.5);
0079
0080   TF1 * f_calo_r = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 25);
0081   TF1 * f_calo_l = new TF1("f_calo_l", "pol2", 0.5, 25);
0082
0083   TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
0084       25);
0085   f_calo_sim->SetParameters(2.4, 11.8);
0086   f_calo_sim->SetLineWidth(1);
0087   f_calo_sim->SetLineColor(kGreen + 2);
0088
0089   TF1 * f_calo_l_sim = new TF1("f_calo_l_sim", "pol2", 0.5, 25);
0090   f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
0091   f_calo_l_sim->SetLineWidth(1);
0092   f_calo_l_sim->SetLineColor(kGreen + 2);
0093
0094   TText * t;
0095   TCanvas *c1 = new TCanvas(
0096       Form("DrawPrototype2EMCalTower_Resolution_Run%.0f_", runs[0]) + cut,
0097       Form("DrawPrototype2EMCalTower_Resolution_Run%.0f_", runs[0]) + cut, 1300,
0098       600);
0099   c1->Divide(2, 1);
0100   int idx = 1;
0101   TPad * p;
0102
0103   p = (TPad *) c1->cd(idx++);
0104   c1->Update();
0105 //  p->SetLogy();
0106   p->SetGridx(0);
0107   p->SetGridy(0);
0108
0109   p->DrawFrame(0, 0, 25, 30,
0110       Form("Linearity;Input energy (GeV);Measured Energy (GeV)", E));
0111   TLine * l = new TLine(0, 0, 25, 25);
0112   l->SetLineColor(kGray);
0113 //  l->Draw();
0114
0115   ge_linear->Draw("p");
0116   ge_linear->Fit(f_calo_l, "RM0");
0117   f_calo_l->Draw("same");
0118   f_calo_l_sim->Draw("same");
0119
0120   p = (TPad *) c1->cd(idx++);
0121   c1->Update();
0122 //  p->SetLogy();
0123   p->SetGridx(0);
0124   p->SetGridy(0);
0125
0126   TH1 * hframe = p->DrawFrame(0, 0, 25, 0.2,
0127       Form("Resolution;Input energy (GeV);#DeltaE/<E>", E));
0128
0129   ge_res->Draw("p");
0130   ge_res->Fit(f_calo_r, "RM0");
0131   f_calo_r->Draw("same");
0132   f_calo_sim->Draw("same");
0133
0134   hframe->SetTitle(
0135       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0136           f_calo_r->GetParameter(0), f_calo_r->GetParameter(1)));
0137
0138   SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
0139 }
0140
0141 vector<double>
0142 GetOnePlot(const TString base, const int run, TString cut)
0143 {
0144 //  beam_00002063-0000_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_Energy_Sum_col1_row2_5x5_tower_center_col1_row2.png
0145 //  TString fname = base + TString("Prototype_") + particle + "_" + sE
0146 //      + "_SegALL_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_all_event.root";
0147   TString fname =
0148       base
0149           + Form(
0150               "/beam_0000%d-0000_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_Energy_Sum_",
0151               run) + cut + ".root";
0152 //  Prototype_e-_2_Seg0_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUMall_event.root
0153
0154   cout << "Process " << fname << endl;
0155   TFile * f = new TFile(fname);
0156   assert(f->IsOpen());
0157
0158   TH1 * hEnergySum = (TH1 *) f->GetObjectChecked("EnergySum_LG", "TH1");
0159   assert(hEnergySum);
0160   new TCanvas();
0161   TH1 * h = hEnergySum->DrawClone();
0162
0163   hEnergySum->Scale(1. / hEnergySum->Integral(1, -1));
0164 //  TF1 * fgaus = hEnergySum->GetFunction("fgaus_LG");
0165 //  assert(fgaus);
0166
0167   TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
0168       h->GetMean() + 4 * h->GetRMS());
0169   fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0170       h->GetMean() + 2 * h->GetRMS());
0171   h->Fit(fgaus_g, "MR0N");
0172
0173   TF1 * fgaus = new TF1("fgaus_LG", "gaus",
0174       fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
0175       fgaus_g->GetParameter(1) + 2 * fgaus_g->GetParameter(2));
0176   fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0177       fgaus_g->GetParameter(2));
0178   fgaus->SetLineColor(kRed);
0179   fgaus->SetLineWidth(4);
0180   h->Fit(fgaus, "MR");
0181
0182   vector<double> v(4);
0183   v[0] = fgaus->GetParameter(1);
0184   v[1] = fgaus->GetParError(1);
0185   v[2] = fgaus->GetParameter(2);
0186   v[3] = fgaus->GetParError(2);
0187
0188   return v;
0189 }
0190