EMCal/macros/DrawEMCalTower.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 TFile * _file0 = NULL;
0024 TTree * T = NULL;
0025 TString cuts = "";
0026
0027 void
0028 DrawEMCalTower( //
0029     const TString infile = // "data/Prototype2_CEMC.root_DSTReader.root" //
0030 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/15Degree_1Col_LightCollection//Prototype_e-_16_SegALL_DSTReader.root"//
0031 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_DSTReader.root"//
0032 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./10DegreeRot_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_DSTReader.root"//
0033 "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/sim/MuonCalib//1k_Muon_DSTReader.root"
0034 )
0035 {
0036   SetOKStyle();
0037   gStyle->SetOptStat(0);
0038   gStyle->SetOptFit(1111);
0039   TVirtualFitter::SetDefaultFitter("Minuit2");
0040   gSystem->Load("libg4eval.so");
0041   gSystem->Load("libqa_modules.so");
0042
0043   if (!_file0)
0044     {
0045       TString chian_str = infile;
0046       chian_str.ReplaceAll("ALL", "*");
0047
0048       TChain * t = new TChain("T");
0049       const int n = t->Add(chian_str);
0050
0051       cout << "Loaded " << n << " root files with " << chian_str << endl;
0052       assert(n > 0);
0053
0054       T = t;
0055
0056       _file0 = new TFile;
0057       _file0->SetName(infile);
0058     }
0059
0060   assert(_file0);
0061
0062   T->SetAlias("ActiveTower_LG",
0063       "TOWER_CALIB_LG_CEMC[].get_binphi()<8 && TOWER_CALIB_LG_CEMC[].get_bineta()<8");
0064   T->SetAlias("EnergySum_LG",
0065       "1*Sum$(TOWER_CALIB_LG_CEMC[].get_energy() * ActiveTower_LG)");
0066
0067   T->SetAlias("ActiveTower_HG",
0068       "TOWER_CALIB_HG_CEMC[].get_binphi()<8 && TOWER_CALIB_HG_CEMC[].get_bineta()<8");
0069   T->SetAlias("EnergySum_HG",
0070       "1*Sum$(TOWER_CALIB_HG_CEMC[].get_energy() * ActiveTower_HG)");
0071
0072 //
0073   const TCut event_sel = "1*1";
0074   cuts = "_all_event";
0075 //  const TCut event_sel = "Entry$==3";
0076 //  cuts = "ev3";
0077 //  const TCut event_sel = "Entry$<2";
0078 //  cuts = "2ev";
0079 //  const TCut event_sel = "Entry$<1000";
0080 //  cuts = "1000ev";
0081
0082   cout << "Build event selection of " << (const char *) event_sel << endl;
0083
0084   T->Draw(">>EventList", event_sel);
0085   TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0086   cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
0087       << endl;
0088
0089   T->SetEventList(elist);
0090
0091   EMCDistribution_Fast("HG");
0092 //  EMCDistribution_Fast("LG", true);
0093 //  EMCDistributionVSBeam_SUM();
0094 //  EMCDistributionVSBeam_SUM("TOWER_CEMC", -0, 5); // 0 degree tilted
0095 //  EMCDistributionVSBeam_SUM("TOWER_CEMC", -15); // 10 degree tilted
0096 //  EMCDistributionVSBeam_SUM("TOWER_CEMC", -15, 5); // 10 degree tilted
0097 //  EMCDistribution_SUM();
0098 }
0099
0100 void
0101 EMCDistribution_Fast(TString gain = "LG", bool log_scale = false)
0102 {
0103   TString hname = "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "")
0104       + cuts;
0105
0106   TH2 * h2 = NULL;
0107   if (log_scale)
0108     {
0109       h2 = new TH2F(hname,
0110           Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300, 10e-3,
0111           100, 64, -.5, 63.5);
0112       QAHistManagerDef::useLogBins(h2->GetXaxis());
0113     }
0114   else
0115     {
0116       h2 = new TH2F(hname,
0117           Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100, -.050,
0118           .5, 64, -.5, 63.5);
0119     }
0120
0121   T->Draw(
0122       "TOWER_CALIB_" + gain + "_CEMC[].get_bineta() + 8* TOWER_CALIB_" + gain
0123           + "_CEMC[].get_binphi():TOWER_CALIB_" + gain
0124           + "_CEMC[].get_energy()>>" + hname, "", "goff");
0125   TText * t;
0126   TCanvas *c1 = new TCanvas(
0127       "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
0128       "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
0129       950);
0130   c1->Divide(8, 8, 0., 0.01);
0131   int idx = 1;
0132   TPad * p;
0133
0134   for (int iphi = 8 - 1; iphi >= 0; iphi--)
0135     {
0136       for (int ieta = 0; ieta < 8; ieta++)
0137         {
0138           p = (TPad *) c1->cd(idx++);
0139           c1->Update();
0140
0141           if (log_scale)
0142             {
0143               p->SetLogy();
0144               p->SetLogx();
0145             }
0146           p->SetGridx(0);
0147           p->SetGridy(0);
0148
0149           TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
0150               + TString(log_scale ? "_Log" : "");
0151
0152           TH1 * h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
0153               ieta + 8 * iphi + 1); // axis bin number is encoded as ieta+8*iphi+1
0154
0155           h->SetLineWidth(0);
0156           h->SetLineColor(kBlue + 3);
0157           h->SetFillColor(kBlue + 3);
0158           h->GetXaxis()->SetTitleSize(.09);
0159           h->GetXaxis()->SetLabelSize(.08);
0160           h->GetYaxis()->SetLabelSize(.08);
0161
0162           h->Draw();
0163
0164           TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
0165           t->SetTextAlign(33);
0166           t->SetTextSize(.15);
0167           t->SetNDC();
0168           t->Draw();
0169         }
0170     }
0171
0172   SaveCanvas(c1,
0173       TString(_file0->GetName()) + TString("_DrawEMCalTower_")
0174           + TString(c1->GetName()), kTRUE);
0175
0176 }
0177
0178 void
0179 EMCDistribution(TString gain = "LG", bool log_scale = false)
0180 {
0181
0182   TText * t;
0183   TCanvas *c1 = new TCanvas(
0184       "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
0185       "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
0186       1000);
0187   c1->Divide(8, 8, 0., 0.01);
0188   int idx = 1;
0189   TPad * p;
0190
0191   for (int iphi = 8 - 1; iphi >= 0; iphi--)
0192     {
0193       for (int ieta = 0; ieta < 8; ieta++)
0194         {
0195           p = (TPad *) c1->cd(idx++);
0196           c1->Update();
0197
0198           if (log_scale)
0199             {
0200               p->SetLogy();
0201               p->SetLogx();
0202             }
0203           p->SetGridx(0);
0204           p->SetGridy(0);
0205
0206           TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
0207               + TString(log_scale ? "_Log" : "");
0208
0209           TH1 * h = NULL;
0210
0211           if (log_scale)
0212             h = new TH1F(hname,
0213                 Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300,
0214                 5e-3, 100);
0215           else
0216             h = new TH1F(hname,
0217                 Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100,
0218                 -.050, .5);
0219
0220           h->SetLineWidth(0);
0221           h->SetLineColor(kBlue + 3);
0222           h->SetFillColor(kBlue + 3);
0223           h->GetXaxis()->SetTitleSize(.09);
0224           h->GetXaxis()->SetLabelSize(.08);
0225           h->GetYaxis()->SetLabelSize(.08);
0226
0227           if (log_scale)
0228             QAHistManagerDef::useLogBins(h->GetXaxis());
0229
0230           T->Draw("TOWER_CALIB_" + gain + "_CEMC[].get_energy()>>" + hname,
0231               Form(
0232                   "TOWER_CALIB_%s_CEMC[].get_bineta()==%d && TOWER_CALIB_%s_CEMC[].get_binphi()==%d",
0233                   gain.Data(), ieta, gain.Data(), iphi), "");
0234
0235           TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
0236           t->SetTextAlign(33);
0237           t->SetTextSize(.15);
0238           t->SetNDC();
0239           t->Draw();
0240         }
0241     }
0242
0243   SaveCanvas(c1,
0244       TString(_file0->GetName()) + TString("_DrawEMCalTower_")
0245           + TString(c1->GetName()), kTRUE);
0246
0247 }
0248 void
0249 EMCDistribution_SUM(TString sTOWER = "TOWER_CEMC")
0250 {
0251   TH1 * EnergySum_LG = new TH1F("EnergySum_LG",
0252       ";Low-gain Tower Energy Sum (GeV);Count / bin", 1500, 0, 100);
0253   TH1 * EnergySum_HG = new TH1F("EnergySum_HG",
0254       ";High-gain Tower Energy Sum (GeV);Count / bin", 300, 0, 3);
0255
0256   T->Draw("EnergySum_LG>>EnergySum_LG", "", "goff");
0257   T->Draw("EnergySum_HG>>EnergySum_HG", "", "goff");
0258
0259   TText * t;
0260   TCanvas *c1 = new TCanvas("EMCDistribution_SUM" + cuts,
0261       "EMCDistribution_SUM" + cuts, 1000, 960);
0262   c1->Divide(2, 2);
0263   int idx = 1;
0264   TPad * p;
0265
0266   p = (TPad *) c1->cd(idx++);
0267   c1->Update();
0268   p->SetLogy();
0269   p->SetGridx(0);
0270   p->SetGridy(0);
0271
0272   TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
0273   h->SetLineWidth(2);
0274   h->SetLineColor(kBlue + 3);
0275 //  h->Sumw2();
0276   h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
0277
0278   p = (TPad *) c1->cd(idx++);
0279   c1->Update();
0280 //  p->SetLogy();
0281   p->SetGridx(0);
0282   p->SetGridy(0);
0283
0284   TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
0285
0286   TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
0287   fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0288       h->GetMean() + 2 * h->GetRMS());
0289   h->Fit(fgaus, "M");
0290
0291   h->Sumw2();
0292   h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0293       h->GetMean() + 4 * h->GetRMS());
0294
0295   h->SetLineWidth(2);
0296   h->SetMarkerStyle(kFullCircle);
0297
0298   h->SetTitle(
0299       Form("#DeltaE/<E> = %.1f%%",
0300           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0301
0302   p = (TPad *) c1->cd(idx++);
0303   c1->Update();
0304   p->SetLogy();
0305   p->SetGridx(0);
0306   p->SetGridy(0);
0307
0308   TH1 * h = (TH1 *) EnergySum_HG->DrawClone();
0309   h->SetLineWidth(2);
0310   h->SetLineColor(kBlue + 3);
0311 //  h->Sumw2();
0312   h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
0313
0314   p = (TPad *) c1->cd(idx++);
0315   c1->Update();
0316 //  p->SetLogy();
0317   p->SetGridx(0);
0318   p->SetGridy(0);
0319
0320   TH1 * h = (TH1 *) EnergySum_HG->DrawClone();
0321
0322   TF1 * fgaus = new TF1("fgaus_HG", "gaus", 0, 100);
0323   fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0324       h->GetMean() + 2 * h->GetRMS());
0325   h->Fit(fgaus, "M");
0326
0327   h->Sumw2();
0328   h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0329       h->GetMean() + 4 * h->GetRMS());
0330
0331   h->SetLineWidth(2);
0332   h->SetMarkerStyle(kFullCircle);
0333
0334   h->SetTitle(
0335       Form("#DeltaE/<E> = %.1f%%",
0336           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0337
0338   SaveCanvas(c1,
0339       TString(_file0->GetName()) + TString("_DrawEMCalTower_")
0340           + TString(c1->GetName()), kTRUE);
0341
0342 }
0343
0344 void
0345 EMCDistributionVSBeam_SUM(TString sTOWER = "TOWER_CEMC", const double z_shift =
0346     0, const int n_div = 1)
0347 {
0348   TH3F * EnergySum_LG3 =
0349       new TH3F("EnergySum_LG3",
0350           ";Beam Horizontal Pos (cm);Beam Horizontal Vertical (cm);Low-gain Tower Energy Sum (GeV)", //
0351           20 * n_div, z_shift - 5, z_shift + 5, //
0352           20 * n_div, -5, 5, //
0353           200, 0, 20);
0354
0355   T->Draw("EnergySum_LG:PHG4VtxPoint.vy:PHG4VtxPoint.vz>>EnergySum_LG3", "",
0356       "goff");
0357
0358   TProfile2D * EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
0359   TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
0360   TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
0361
0362   TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
0363   TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
0364
0365   TText * t;
0366   TCanvas *c1 = new TCanvas(
0367       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d", n_div)) + cuts,
0368       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d", n_div)) + cuts, 1000,
0369       960);
0370   c1->Divide(2, 2);
0371   int idx = 1;
0372   TPad * p;
0373
0374   p = (TPad *) c1->cd(idx++);
0375   c1->Update();
0376 //  p->SetLogy();
0377   p->SetGridx(0);
0378   p->SetGridy(0);
0379
0380   EnergySum_LG3_xy->Draw("colz");
0381   EnergySum_LG3_xy->SetTitle(
0382       "Position scan;Beam Horizontal Pos (cm);Beam Horizontal Vertical (cm)");
0383
0384   p = (TPad *) c1->cd(idx++);
0385   c1->Update();
0386 //  p->SetLogy();
0387   p->SetGridx(0);
0388   p->SetGridy(0);
0389
0390   EnergySum_LG3_zx->Draw("colz");
0391   EnergySum_LG3_zx->SetTitle(
0392       "Position scan;Beam Horizontal Pos (cm);Energy Sum (GeV)");
0393
0394   ge_EnergySum_LG3_zx->SetLineWidth(2);
0395   ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
0396   ge_EnergySum_LG3_zx->Draw("pe");
0397
0398   p = (TPad *) c1->cd(idx++);
0399   c1->Update();
0400 //  p->SetLogy();
0401   p->SetGridx(0);
0402   p->SetGridy(0);
0403
0404   EnergySum_LG3_zy->Draw("colz");
0405   EnergySum_LG3_zy->SetTitle(
0406       "Position scan;Beam Vertical Pos (cm);Energy Sum (GeV)");
0407
0408   ge_EnergySum_LG3_zy->SetLineWidth(2);
0409   ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
0410   ge_EnergySum_LG3_zy->Draw("pe");
0411
0412   p = (TPad *) c1->cd(idx++);
0413   c1->Update();
0414 //  p->SetLogy();
0415   p->SetGridx(0);
0416   p->SetGridy(0);
0417
0418   TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
0419
0420   TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
0421   fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0422       h->GetMean() + 2 * h->GetRMS());
0423   h->Fit(fgaus, "M");
0424
0425   h->Sumw2();
0426   h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0427       h->GetMean() + 4 * h->GetRMS());
0428   EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0429       h->GetMean() + 4 * h->GetRMS());
0430   EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0431       h->GetMean() + 4 * h->GetRMS());
0432
0433   h->SetLineWidth(2);
0434   h->SetMarkerStyle(kFullCircle);
0435
0436   h->SetTitle(
0437       Form("#DeltaE/<E> = %.1f%%",
0438           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0439
0440   SaveCanvas(c1,
0441       TString(_file0->GetName()) + TString("_DrawEMCalTower_")
0442           + TString(c1->GetName()), kTRUE);
0443
0444 }
0445
0446 TGraphErrors *
0447 FitProfile(const TH2 * h2)
0448 {
0449
0450   TProfile * p2 = h2->ProfileX();
0451
0452   int n = 0;
0453   double x[1000];
0454   double ex[1000];
0455   double y[1000];
0456   double ey[1000];
0457
0458   for (int i = 1; i <= h2->GetNbinsX(); i++)
0459     {
0460       TH1D * h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0461
0462       if (h1->GetSum() < 30)
0463         {
0464           cout << "FitProfile - ignore bin " << i << endl;
0465           continue;
0466         }
0467       else
0468         {
0469           cout << "FitProfile - fit bin " << i << endl;
0470         }
0471
0472       TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0473           p2->GetBinError(i) * 4);
0474
0475       TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0476           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0477
0478       fgaus.SetParameter(1, p2->GetBinContent(i));
0479       fgaus.SetParameter(2, p2->GetBinError(i));
0480
0481       h1->Fit(&fgaus, "MQ");
0482
0483       f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
0484           fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
0485           fgaus.GetParameter(2) / 4, 0);
0486
0487       h1->Fit(&f2, "MQ");
0488
0489 //      new TCanvas;
0490 //      h1->Draw();
0491 //      fgaus.Draw("same");
0492 //      break;
0493
0494       x[n] = p2->GetBinCenter(i);
0495       ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0496       y[n] = fgaus.GetParameter(1);
0497       ey[n] = fgaus.GetParError(1);
0498
0499 //      p2->SetBinContent(i, fgaus.GetParameter(1));
0500 //      p2->SetBinError(i, fgaus.GetParameter(2));
0501
0502       n++;
0503       delete h1;
0504     }
0505
0506   return new TGraphErrors(n, x, y, ex, ey);
0507 }
0508