EMCal/macros/DrawPrototype2EMCalTower_ResolutionRecalib.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 bool mask = true;
0024
0025 void
0026 DrawPrototype2EMCalTower_ResolutionRecalib(
0027 //
0028     const TString base =
0029         "/phenix/u/jinhuang/tmp/miliped_work/Production_0414/calirbated_")
0030 {
0031   SetOKStyle();
0032   gStyle->SetOptStat(0);
0033   gStyle->SetGridStyle(0);
0034   if (mask)
0035     gStyle->SetOptFit(0);
0036   else
0037     gStyle->SetOptFit(0);
0038   TVirtualFitter::SetDefaultFitter("Minuit2");
0039
0040   const int N = 6;
0041   const double Es[] =
0042     { 2, 3, 4, 8, 12, 16 };
0043   const double runs[] =
0044     { 2042, 2040, 2039, 2038, 2067, 2063 };
0045
0046   double mean[1000];
0047   double dmean[1000];
0048   double res[1000];
0049   double dres[1000];
0050
0051   TText * t;
0052   TCanvas *c1 = new TCanvas(
0053       "DrawPrototype2EMCalTower_ResolutionRecalib_RunSummary"
0054           + TString(mask ? "Mask" : ""),
0055       "DrawPrototype2EMCalTower_ResolutionRecalib_RunSummary", 1800, 900);
0056   c1->Divide(3, 2);
0057   int idx = 1;
0058   TPad * p;
0059   for (int i = 0; i < N; ++i)
0060     {
0061
0062       p = (TPad *) c1->cd(idx++);
0063       c1->Update();
0064
0065       p->SetGridx(0);
0066       p->SetGridy(0);
0067
0068       const double E = Es[i];
0069
0070       vector<double> v(4);
0071
0072 //      v = GetOnePlot(base, runs[i], cut);
0073       v = RecalibratorMakeOnePlot(base, runs[i], E);
0074
0075       mean[i] = v[0];
0076       dmean[i] = v[1];
0077       res[i] = v[2] / v[0];
0078       dres[i] = v[3] / v[0];
0079
0080       cout << "mean[i] = " << mean[i] << ", " //
0081           << "dmean[i] = " << dmean[i] << ", " //
0082           << "res[i] = " << res[i] << ", " //
0083           << "dres[i] = " << dres[i] << endl;
0084     }
0085
0086   SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
0087
0088   TGraphErrors * ge_linear = new TGraphErrors(N, Es, mean, 0, dmean);
0089   TGraphErrors * ge_res = new TGraphErrors(N, Es, res, 0, dres);
0090
0091   ge_linear->SetLineColor(kBlue + 3);
0092   ge_linear->SetMarkerColor(kBlue + 3);
0093   ge_linear->SetLineWidth(2);
0094   ge_linear->SetMarkerStyle(kFullCircle);
0095   ge_linear->SetMarkerSize(2);
0096
0097   ge_res->SetLineColor(kBlue + 3);
0098   ge_res->SetMarkerColor(kBlue + 3);
0099   ge_res->SetLineWidth(2);
0100   ge_res->SetMarkerStyle(kFullCircle);
0101   ge_res->SetMarkerSize(2);
0102
0103   TF1 * f_calo_r = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 25);
0104   TF1 * f_calo_l = new TF1("f_calo_l", "pol2", 0.5, 25);
0105
0106   f_calo_r->SetLineColor(kBlue + 3);
0107   f_calo_r->SetLineWidth(3);
0108
0109   TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
0110       25);
0111   f_calo_sim->SetParameters(2.4, 11.8);
0112   f_calo_sim->SetLineWidth(3);
0113   f_calo_sim->SetLineColor(kGreen + 2);
0114
0115   TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
0116 //  f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
0117   f_calo_l_sim->SetParameters(-0., 1, -0.);
0118 //  f_calo_l_sim->SetLineWidth(1);
0119   f_calo_l_sim->SetLineColor(kGreen + 2);
0120   f_calo_l_sim->SetLineWidth(3);
0121
0122   TFile * f_ref =
0123       new TFile(
0124           "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/Prototype_2016/Production_0417_CEMC_MIP_set2_v3/DrawPrototype2EMCalTower_Resolution_Run2042_col1_row2_5x5_Valid_HODO_center_col1_row2.root");
0125   assert(f_ref->IsOpen());
0126
0127   TGraphErrors * ge_res_ref = (TGraphErrors *) f_ref->GetObjectChecked("Graph",
0128       "TGraphErrors");
0129   assert(ge_res_ref);
0130   ge_res_ref->SetObjectStat(0);
0131
0132   TF1 * f_calo_r_ref = (TF1 *) f_ref->GetObjectChecked("f_calo_r", "TF1");
0133   assert(f_calo_r_ref);
0134
0135   ge_res_ref->SetLineColor(kBlue);
0136   ge_res_ref->SetMarkerColor(kBlue);
0137   ge_res_ref->SetLineWidth(2);
0138   ge_res_ref->SetMarkerStyle(kOpenCircle);
0139   ge_res_ref->SetMarkerSize(2);
0140
0141   f_calo_r_ref->SetLineColor(kBlue);
0142   f_calo_r_ref->SetMarkerColor(kBlue);
0143   f_calo_r_ref->SetLineWidth(2);
0144
0145   TText * t;
0146   TCanvas *c1 = new TCanvas(
0147       Form("DrawPrototype2EMCalTower_ResolutionRecalib%.0f_", runs[0]),
0148       Form("DrawPrototype2EMCalTower_ResolutionRecalib%.0f_", runs[0]), 1300,
0149       600);
0150   c1->Divide(2, 1);
0151   int idx = 1;
0152   TPad * p;
0153
0154   p = (TPad *) c1->cd(idx++);
0155   c1->Update();
0156 //  p->SetLogy();
0157   p->SetGridx(0);
0158   p->SetGridy(0);
0159
0160   TLegend* leg = new TLegend(.15, .7, .6, .85);
0161
0162   p->DrawFrame(0, 0, 25, 25,
0163       Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
0164   TLine * l = new TLine(0, 0, 25, 25);
0165   l->SetLineColor(kGray);
0166 //  l->Draw();
0167
0168   f_calo_l_sim->Draw("same");
0169   ge_linear->Draw("p");
0170 //  ge_linear->Fit(f_calo_l, "RM0");
0171 //  f_calo_l->Draw("same");
0172
0173   leg->AddEntry(ge_linear, "T-1044 data", "ep");
0174   leg->AddEntry(f_calo_l_sim, "Unity", "l");
0175   leg->Draw();
0176
0177   p = (TPad *) c1->cd(idx++);
0178   c1->Update();
0179 //  p->SetLogy();
0180   p->SetGridx(0);
0181   p->SetGridy(0);
0182
0183   TH1 * hframe = p->DrawFrame(0, 0, 25, 0.2,
0184       Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
0185
0186   ge_res->Fit(f_calo_r, "RM0Q");
0187
0188   f_calo_r_ref->Draw("same");
0189   ge_res_ref->Draw("p");
0190   f_calo_r->Draw("same");
0191   f_calo_sim->Draw("same");
0192   ge_res->Draw("p");
0193
0194   TLegend* leg = new TLegend(.2, .6, .85, .9);
0195   leg->AddEntry(ge_res_ref, "T-1044, MIP calib.", "ep");
0196   leg->AddEntry(f_calo_r_ref,
0197       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0198           f_calo_r_ref->GetParameter(0), f_calo_r_ref->GetParameter(1)), "l");
0199 //  leg->AddEntry(new TH1(), "", "l");
0200   leg->AddEntry((TObject*)0, " ", "");
0201   leg->AddEntry(ge_res, "T-1044, e-shower calib.", "ep");
0202   leg->AddEntry(f_calo_r,
0203       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0204           f_calo_r->GetParameter(0), f_calo_r->GetParameter(1)), "l");
0205   leg->Draw();
0206
0207   TLegend* leg = new TLegend(.2, .15, .85, .3);
0208
0209   leg->AddEntry(f_calo_sim,
0210       Form("Simulation, #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0211           f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
0212   leg->Draw();
0213
0214   hframe->SetTitle("Electron Resolution");
0215
0216   SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
0217 }
0218
0219 vector<double>
0220 RecalibratorMakeOnePlot(const TString base, const int run, const double E)
0221 {
0222 //  /phenix/u/jinhuang/tmp/miliped_work/Production_0414/calirbated_2040.dat
0223
0224   TString fname = base + Form("%d.dat", run);
0225 //  Prototype_e-_2_Seg0_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUMall_event.root
0226
0227   cout << "Process " << fname << endl;
0228
0229   fstream f(fname, ios_base::in);
0230   assert(f.is_open());
0231
0232   TH1 * EnergySum_LG_production = new TH1F(
0233       Form("EnergySum_LG_production_%d", run),
0234       Form(";Run %d: 5x5 EMCal Energy Sum (GeV);Count / bin", run), 400, 0, 30);
0235   TH1 * EnergySum_LG_recalib = new TH1F(Form("EnergySum_LG_recalib_%d", run),
0236       Form(";Run %d: 5x5 EMCal Energy Sum (GeV);Count / bin", run), 400, 0, 30);
0237
0238   int count = 0;
0239   while (!f.eof())
0240     {
0241       string line;
0242       getline(f, line);
0243
0244       if (line.length() == 0)
0245         continue;
0246
0247       double old_E = -1, new_E = -1;
0248       char tab;
0249
0250       stringstream sline(line);
0251       sline >> old_E >> tab >> new_E;
0252
0253 //      cout <<"line -> old_E = "<<old_E<<", new_E = "<<new_E<<endl;
0254       assert(old_E > 0);
0255       assert(new_E > 0);
0256
0257       EnergySum_LG_production->Fill(old_E);
0258       EnergySum_LG_recalib->Fill(new_E);
0259
0260     }
0261
0262   EnergySum_LG_production->SetLineColor(kBlue + 3);
0263   EnergySum_LG_production->SetLineWidth(2);
0264   EnergySum_LG_production->SetMarkerStyle(kFullCircle);
0265   EnergySum_LG_production->SetMarkerColor(kBlue + 3);
0266   EnergySum_LG_recalib->SetLineColor(kRed + 3);
0267   EnergySum_LG_recalib->SetLineWidth(2);
0268   EnergySum_LG_recalib->SetMarkerStyle(kFullCircle);
0269   EnergySum_LG_recalib->SetMarkerColor(kRed + 3);
0270
0271   EnergySum_LG_recalib->Sumw2();
0272
0273   TH1 * h = (TH1 *) EnergySum_LG_recalib->DrawClone();
0274   if (!mask)
0275     EnergySum_LG_production->DrawClone("same");
0276
0277   TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
0278       h->GetMean() + 4 * h->GetRMS());
0279   fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0280       h->GetMean() + 2 * h->GetRMS());
0281   h->Fit(fgaus_g, "MR0N");
0282
0283   TF1 * fgaus = new TF1("fgaus_LG", "gaus",
0284       fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
0285       fgaus_g->GetParameter(1) + 2 * fgaus_g->GetParameter(2));
0286   fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0287       fgaus_g->GetParameter(2));
0288   fgaus->SetLineColor(kRed);
0289   fgaus->SetLineWidth(4);
0290   h->Fit(fgaus, "MR");
0291
0292 //  h->Sumw2();
0293   h->GetXaxis()->SetRangeUser(h->GetMean() - 5 * h->GetRMS(),
0294       h->GetMean() + 5 * h->GetRMS());
0295   if (!mask)
0296     h->SetTitle(
0297         Form("#DeltaE/<E> = %.1f%% @ Beam = %.0f GeV",
0298             100 * fgaus->GetParameter(2) / fgaus->GetParameter(1), E));
0299   else
0300     h->SetTitle(Form("Beam = %.0f GeV", E));
0301
0302   assert(fgaus);
0303
0304   vector<double> v(4);
0305   v[0] = fgaus->GetParameter(1);
0306   v[1] = fgaus->GetParError(1);
0307   v[2] = fgaus->GetParameter(2);
0308   v[3] = fgaus->GetParError(2);
0309
0310   return v;
0311 }