Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:14:34 

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