17 #include <TGraphErrors.h>
33 "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/Production_1101_EMCalSet2_HCalPR12/beam_00002187-0000_DSTReader.root",
34 bool plot_all =
false,
const double momentum = 120)
37 gROOT->SetStyle(
"Modern");
41 gStyle->SetOptStat(0);
42 gStyle->SetOptFit(1111);
43 TVirtualFitter::SetDefaultFitter(
"Minuit2");
44 gSystem->Load(
"libg4eval.so");
45 gSystem->Load(
"libqa_modules.so");
46 gSystem->Load(
"libPrototype2.so");
52 TString chian_str =
infile;
53 chian_str.ReplaceAll(
"ALL",
"*");
55 TChain *
t =
new TChain(
"T");
56 const int n = t->Add(chian_str);
58 cout <<
"Loaded " << n <<
" root files with " << chian_str << endl;
69 T->SetAlias(
"ActiveTower_LG",
70 "TOWER_LG_CEMC[].get_binphi()<8 && TOWER_LG_CEMC[].get_bineta()<8");
71 T->SetAlias(
"EnergySum_LG",
72 "1*Sum$(TOWER_LG_CEMC[].get_energy() * ActiveTower_LG)");
74 T->SetAlias(
"ActiveTower_HG",
75 "TOWER_HG_CEMC[].get_binphi()<8 && TOWER_HG_CEMC[].get_bineta()<8");
76 T->SetAlias(
"EnergySum_HG",
77 "1*Sum$(TOWER_HG_CEMC[].get_energy() * ActiveTower_HG)");
80 T->SetAlias(
"C2_Inner_e",
"1*abs(TOWER_RAW_C2[2].energy)");
81 T->SetAlias(
"C2_Outer_e",
"1*abs(TOWER_RAW_C2[3].energy)");
82 T->SetAlias(
"C2_Sum_e",
"C2_Inner_e + C2_Outer_e");
86 T->SetAlias(
"Average_column",
87 "Sum$(TOWER_CALIB_CEMC.get_column() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
88 T->SetAlias(
"Average_row",
89 "Sum$(TOWER_CALIB_CEMC.get_row() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
91 T->SetAlias(
"Average_HODO_VERTICAL",
92 "Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
96 T->SetAlias(
"Valid_HODO_VERTICAL",
97 "Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30 && TOWER_CALIB_HODO_VERTICAL.towerid==3 ) ");
100 T->SetAlias(
"Average_HODO_HORIZONTAL",
101 "Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
104 T->SetAlias(
"Valid_HODO_HORIZONTAL",
105 "Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30)==1");
107 T->SetAlias(
"No_Triger_VETO",
108 "Sum$(abs(TOWER_RAW_TRIGGER_VETO.energy)>15)==0");
110 T->SetAlias(
"Energy_Sum_col1_row2_3x3",
111 "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-1)<=1 && abs(TOWER_CALIB_CEMC.get_row()-2)<=1 ) * TOWER_CALIB_CEMC.get_energy())");
112 T->SetAlias(
"Energy_Sum_col1_row2_5x5",
113 "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-1)<=2 && abs(TOWER_CALIB_CEMC.get_row()-2)<=2 ) * TOWER_CALIB_CEMC.get_energy())");
114 T->SetAlias(
"Energy_Sum_col2_row2_5x5",
115 "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-2)<=2 && abs(TOWER_CALIB_CEMC.get_row()-2)<=2 ) * TOWER_CALIB_CEMC.get_energy())");
116 T->SetAlias(
"Energy_Sum_CEMC",
"1*Sum$(TOWER_CALIB_CEMC.get_energy())");
127 T->SetAlias(
"Energy_Sum_Hadron_CEMC",
128 "1.14*12./8.71776e+00*Sum$(TOWER_CALIB_CEMC.get_energy())");
131 T->SetAlias(
"CEMC_MIP",
"Energy_Sum_Hadron_CEMC<0.7");
146 T->SetAlias(
"Energy_Sum_Hadron_HCALIN",
147 "12./6.99727e+00*Sum$(TOWER_CALIB_LG_HCALIN.get_energy())");
148 T->SetAlias(
"HCALIN_MIP",
"Energy_Sum_Hadron_HCALIN<0.5");
149 T->SetAlias(
"Energy_Sum_Hadron_HCALOUT",
150 "12./9.50430e+00*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");
152 T->SetAlias(
"MIP_Count_Col2",
153 "Sum$( abs( TOWER_RAW_CEMC.get_energy() )>20 && abs( TOWER_RAW_CEMC.get_energy() )<400 && TOWER_CALIB_CEMC.get_column() == 2)");
155 T->SetAlias(
"Pedestal_Count_AllCEMC",
156 "Sum$( abs( TOWER_RAW_CEMC.get_energy() )<20 && TOWER_CALIB_CEMC.get_column() != 2)");
157 T->SetAlias(
"TowerTimingGood_Count_AllCEMC",
158 "Sum$( abs( TOWER_RAW_CEMC.get_time() )>6 && abs( TOWER_RAW_CEMC.get_time() )<13 && TOWER_CALIB_CEMC.get_column() == 2)");
170 "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO";
171 cuts =
"_Valid_HODO_Trigger_VETO";
189 "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && (Pedestal_Count_AllCEMC >= 64 -8) && (TowerTimingGood_Count_AllCEMC == 8)";
190 cuts =
"_Valid_HODO_MIP_Col2_PedestalOther_TimingGood";
220 cout <<
"Build event selection of " << (
const char *) event_sel << endl;
222 T->Draw(
">>EventList", event_sel);
223 TEventList * elist = gDirectory->GetObjectChecked(
"EventList",
"TEventList");
224 cout << elist->GetN() <<
" / " <<
T->GetEntriesFast() <<
" events selected"
227 T->SetEventList(elist);
236 gDirectory->mkdir(Form(
"dir_%d", rnd));
237 gDirectory->cd(Form(
"dir_%d", rnd));
242 gDirectory->mkdir(Form(
"dir_%d", rnd));
243 gDirectory->cd(Form(
"dir_%d", rnd));
248 gDirectory->mkdir(Form(
"dir_%d", rnd));
249 gDirectory->cd(Form(
"dir_%d", rnd));
254 gDirectory->mkdir(Form(
"dir_%d", rnd));
255 gDirectory->cd(Form(
"dir_%d", rnd));
269 TString hname =
"EMCDistribution_" +
gain
270 + TString(full_gain ?
"_FullGain" :
"") +
cuts;
273 if (
gain.Contains(
"CALIB"))
278 Form(
";Calibrated Tower Energy Sum (GeV);Count / bin"), 100, .05,
285 Form(
";Calibrated Tower Energy Sum (GeV);Count / bin"), 260, -.2,
289 "TOWER_" +
gain +
"_CEMC[].get_bineta() + 8* TOWER_" +
gain
290 +
"_CEMC[].get_binphi():TOWER_" +
gain +
"_CEMC[].get_energy()>>"
291 + hname,
"",
"goff");
293 else if (
gain.Contains(
"RAW"))
298 Form(
";Calibrated Tower Energy Sum (ADC);Count / bin"), 100,
299 .05 * 100, 25 * 100, 64, -.5, 63.5);
305 Form(
";Calibrated Tower Energy Sum (ADC);Count / bin"), 260,
306 -.2 * 100, 5 * 100, 64, -.5, 63.5);
309 "TOWER_" +
gain +
"_CEMC[].get_bineta() + 8* TOWER_" +
gain
310 +
"_CEMC[].get_binphi():TOWER_" +
gain
311 +
"_CEMC[].get_energy()*(-1)>>" + hname,
"",
"goff");
315 TCanvas *c1 =
new TCanvas(
316 "EMCDistribution_" +
gain + TString(full_gain ?
"_FullGain" :
"") + cuts,
317 "EMCDistribution_" +
gain + TString(full_gain ?
"_FullGain" :
"") + cuts,
319 c1->Divide(8, 8, 0., 0.01);
323 for (
int iphi = 8 - 1; iphi >= 0; iphi--)
325 for (
int ieta = 0; ieta < 8; ieta++)
327 p = (TPad *) c1->cd(idx++);
338 TString hname = Form(
"hEnergy_ieta%d_iphi%d", ieta, iphi)
339 + TString(full_gain ?
"_FullGain" :
"");
341 TH1 *
h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
342 ieta + 8 * iphi + 1);
345 h->SetLineColor(kBlue + 3);
346 h->SetFillColor(kBlue + 3);
348 h->GetXaxis()->SetTitleSize(.09);
349 h->GetXaxis()->SetLabelSize(.08);
350 h->GetYaxis()->SetLabelSize(.08);
355 h->Fit(
"x*gaus",
"M");
357 h->Fit(
"landau",
"M");
361 TF1 *
fit = ((TF1 *) (h->GetListOfFunctions()->At(0)));
365 fit->SetLineColor(kRed);
366 peak = fit->GetParameter(1);
370 cout << Form(
"Finished <Col%d Row%d> = %.1f", ieta, iphi, peak)
373 TText *t =
new TText(.9, .9,
374 Form(
"<Col%d Row%d> = %.1f", ieta, iphi, peak));
383 TString(
_file0->GetName()) + TString(
"_DrawPrototype2MIPAnalysis_")
384 + TString(c1->GetName()),
false);
392 const TString
gain =
"RAW";
394 TString hname =
"EMCDistribution_" + gain
395 + TString(full_gain ?
"_FullGain" :
"") +
cuts;
402 Form(
";Calibrated Tower Energy Sum (ADC);Count / bin"), 100,
403 .05 * 100, 25 * 100, 64, -.5, 63.5);
409 Form(
";Calibrated Tower Energy Sum (ADC);Count / bin"), 260,
410 -.2 * 100, 5 * 100, 64, -.5, 63.5);
413 "TOWER_" + gain +
"_CEMC[].get_bineta() + 8* TOWER_" + gain
414 +
"_CEMC[].get_binphi():(TOWER_RAW_CEMC[].signal_samples[10] - TOWER_RAW_CEMC[].signal_samples[0])*(-1)>>"
415 + hname,
"",
"goff");
419 TCanvas *c1 =
new TCanvas(
420 "EMCDistribution_PeakSample_Fast_" + TString(full_gain ?
"_FullGain" :
"")
422 "EMCDistribution_PeakSample_Fast_" + TString(full_gain ?
"_FullGain" :
"")
424 c1->Divide(8, 8, 0., 0.01);
428 for (
int iphi = 8 - 1; iphi >= 0; iphi--)
430 for (
int ieta = 0; ieta < 8; ieta++)
432 p = (TPad *) c1->cd(idx++);
443 TString hname = Form(
"hEnergy_ieta%d_iphi%d", ieta, iphi)
444 + TString(full_gain ?
"_FullGain" :
"");
446 TH1 *
h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
447 ieta + 8 * iphi + 1);
450 h->SetLineColor(kBlue + 3);
451 h->SetFillColor(kBlue + 3);
453 h->GetXaxis()->SetTitleSize(.09);
454 h->GetXaxis()->SetLabelSize(.08);
455 h->GetYaxis()->SetLabelSize(.08);
460 h->Fit(
"x*gaus",
"M");
462 h->Fit(
"landau",
"M");
466 TF1 *
fit = ((TF1 *) (h->GetListOfFunctions()->At(0)));
470 fit->SetLineColor(kRed);
471 peak = fit->GetParameter(1);
475 cout << Form(
"Finished <Col%d Row%d> = %.1f", ieta, iphi, peak)
478 TText *t =
new TText(.9, .9,
479 Form(
"<Col%d Row%d> = %.1f", ieta, iphi, peak));
488 TString(
_file0->GetName()) + TString(
"_DrawPrototype2MIPAnalysis_")
489 + TString(c1->GetName()),
false);
496 TString
gain =
"RAW";
499 TCanvas *c1 =
new TCanvas(
500 "EMCDistribution_ADC_" + gain + TString(log_scale ?
"_Log" :
"") +
cuts,
501 "EMCDistribution_ADC_" + gain + TString(log_scale ?
"_Log" :
"") +
cuts,
503 c1->Divide(8, 8, 0., 0.01);
507 for (
int iphi = 8 - 1; iphi >= 0; iphi--)
509 for (
int ieta = 0; ieta < 8; ieta++)
511 p = (TPad *) c1->cd(idx++);
521 TString hname = Form(
"hEnergy_ieta%d_iphi%d", ieta, iphi)
522 + TString(log_scale ?
"_Log" :
"");
528 Form(
";Calibrated Tower Energy Sum (GeV);Count / bin"), 24, -.5,
538 h->SetLineColor(kBlue + 3);
539 h->SetFillColor(kBlue + 3);
540 h->GetXaxis()->SetTitleSize(.09);
541 h->GetXaxis()->SetLabelSize(.08);
542 h->GetYaxis()->SetLabelSize(.08);
547 TString sdraw =
"TOWER_" + gain
548 +
"_CEMC[].signal_samples[]:fmod(Iteration$,24)>>" + hname;
551 "TOWER_%s_CEMC[].get_bineta()==%d && TOWER_%s_CEMC[].get_binphi()==%d",
552 gain.Data(), ieta, gain.Data(), iphi);
554 cout <<
"T->Draw(\"" << sdraw <<
"\",\"" << scut <<
"\");" << endl;
556 T->Draw(sdraw, scut,
"colz");
558 TText *t =
new TText(.9, .9, Form(
"Col%d Row%d", ieta, iphi));
569 TString(
_file0->GetName()) + TString(
"_DrawPrototype2MIPAnalysis_")
570 + TString(c1->GetName()),
false);