File indexing completed on 2025-08-05 08:12:15
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include <cmath>
0012 #include <TFile.h>
0013 #include <TString.h>
0014 #include <TLine.h>
0015 #include <TTree.h>
0016 #include <cassert>
0017 #include "SaveCanvas.C"
0018 #include "SetOKStyle.C"
0019 using namespace std;
0020
0021 TFile * _file0 = NULL;
0022 TTree * T = NULL;
0023 TString cuts = "";
0024
0025 void
0026 DrawEcal_Likelihood(
0027
0028 TString base_dir =
0029 "../..//sPHENIX_work/production_analysis_cemc2x2/emcstudies/pidstudies/spacal2d/fieldmap/",
0030
0031
0032
0033 TString pid = "e-",
0034 TString kine_config = "eta0_8GeV", int eval_mode = 1)
0035 {
0036
0037 const TString infile = base_dir + "G4Hits_sPHENIX_" + pid + "_" + kine_config
0038 + "-ALL.root_Ana.root.lst_EMCalLikelihood.root";
0039
0040
0041 SetOKStyle();
0042 gStyle->SetOptStat(0);
0043 gStyle->SetOptFit(1111);
0044 TVirtualFitter::SetDefaultFitter("Minuit2");
0045
0046 gSystem->Load("libg4eval.so");
0047 gSystem->Load("libemcal_ana.so");
0048 gSystem->Load("libg4vertex.so");
0049
0050 if (!_file0)
0051 {
0052 TString chian_str = infile;
0053 chian_str.ReplaceAll("ALL", "*");
0054 chian_str.ReplaceAll("+", "\\+");
0055
0056
0057 TChain * t = new TChain("T");
0058 const int n = t->Add(chian_str);
0059
0060 cout << "Loaded " << n << " root files with " << chian_str << endl;
0061 assert(n>0);
0062
0063 T = t;
0064
0065 _file0 = new TFile;
0066 _file0->SetName(infile);
0067
0068 fstream flst(infile + ".lst", ios_base::out);
0069
0070 TObjArray *fileElements = t->GetListOfFiles();
0071 TIter next(fileElements);
0072 TChainElement *chEl = 0;
0073 while ((chEl = (TChainElement*) next()))
0074 {
0075 flst << chEl->GetTitle() << endl;
0076 }
0077 flst.close();
0078
0079 cout << "Saved file list to " << infile + ".lst" << endl;
0080 }
0081
0082 assert(_file0);
0083
0084 T->SetAlias("UpsilonPair_trk_gpt",
0085 "1*sqrt(DST.UpsilonPair.trk.gpx**2 + DST.UpsilonPair.trk.gpy**2)");
0086 T->SetAlias("UpsilonPair_trk_pt",
0087 "1*sqrt(DST.UpsilonPair.trk.px**2 + DST.UpsilonPair.trk.py**2)");
0088
0089 T->SetAlias("EMCalTrk_pt", "1*sqrt(DST.EMCalTrk.px**2 + DST.EMCalTrk.py**2)");
0090 T->SetAlias("EMCalTrk_gpt",
0091 "1*sqrt(DST.EMCalTrk.gpx**2 + DST.EMCalTrk.gpy**2)");
0092
0093 TCut event_sel = "1*1";
0094 cuts = "_all_event";
0095 if (eval_mode == 0)
0096 {
0097 event_sel = "Entry$<50000 && fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05";
0098 cuts = "_ll_sample";
0099 }
0100 else if (eval_mode == 1)
0101 {
0102
0103 event_sel = " fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05";
0104 cuts = "_eval_sample";
0105 }
0106 cout << "Build event selection of " << (const char *) event_sel << endl;
0107
0108 T->Draw(">>EventList", event_sel);
0109 TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0110 cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
0111 << endl;
0112
0113 T->SetEventList(elist);
0114
0115 if (eval_mode == 0)
0116 {
0117 Edep_Distribution(infile);
0118 }
0119 else if (eval_mode == 1)
0120 {
0121 Edep_LL_Distribution(infile);
0122 EP_LL_Distribution(infile);
0123 }
0124
0125
0126
0127 Edep_Checks(infile, "1");
0128 Ep_Checks(infile, "1");
0129
0130 }
0131
0132 void
0133 EP_LL_Distribution(TString infile)
0134 {
0135
0136 TCanvas *c1 = new TCanvas("EP_LL_Distribution" + cuts,
0137 "EP_LL_Distribution" + cuts, 1900, 900);
0138 c1->Divide(2, 2);
0139 int idx = 1;
0140 TPad * p;
0141
0142 p = (TPad *) c1->cd(idx++);
0143 c1->Update();
0144 p->SetLogy();
0145 T->Draw("DST.EMCalTrk.ll_ep_e>>hll_ep_e(300,-30,30)", "", "");
0146 hll_ep_e->SetTitle(
0147 Form(
0148 "EMCal only: Electron likelihood distribution;log(Electron likelihood);A. U."));
0149
0150 p = (TPad *) c1->cd(idx++);
0151 c1->Update();
0152 p->SetLogy();
0153 T->Draw("DST.EMCalTrk.ll_ep_h>>hll_ep_h(300,-30,30)", "", "");
0154 hll_ep_h->SetTitle(
0155 Form(
0156 "EMCal only: Hadron likelihood distribution;log(Hadron likelihood);A. U."));
0157
0158 p = (TPad *) c1->cd(idx++);
0159 c1->Update();
0160 p->SetLogy();
0161 T->Draw(
0162 "DST.EMCalTrk.ll_ep_e - DST.EMCalTrk.ll_ep_h>>hll_ep_diff(300,-30,30)",
0163 "", "");
0164 TH1F *hll_ep_diff = (TH1F *) gDirectory->GetObjectChecked("hll_ep_diff",
0165 "TH1F");
0166
0167 hll_ep_diff->SetTitle(
0168 Form(
0169 "EMCal only: log likelihood cut;log(Electron likelihood) - log(Hadron likelihood);Count / bin"));
0170
0171
0172 p = (TPad *) c1->cd(idx++);
0173 c1->Update();
0174
0175
0176 p->DrawFrame(-30, 1e-4, 30, 1,
0177 "EMCal only: Cut Efficiency;Cut on log(Electron likelihood) - log(Hadron likelihood);Cut Efficiency");
0178 TGraphErrors * ge = Distribution2Efficiency(hll_ep_diff);
0179 ge->SetLineColor(kBlue + 2);
0180 ge->SetMarkerColor(kBlue + 21);
0181 ge->SetMarkerColor(kFullCircle);
0182 ge->SetLineWidth(3);
0183 ge->Draw("lp");
0184
0185 SaveCanvas(c1,
0186 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0187 + TString(c1->GetName()), kFALSE);
0188
0189 }
0190
0191 void
0192 Edep_LL_Distribution(TString infile)
0193 {
0194
0195 TCanvas *c1 = new TCanvas("Edep_LL_Distribution" + cuts,
0196 "Edep_LL_Distribution" + cuts, 1900, 900);
0197 c1->Divide(2, 2);
0198 int idx = 1;
0199 TPad * p;
0200
0201 p = (TPad *) c1->cd(idx++);
0202 c1->Update();
0203 p->SetLogy();
0204 T->Draw("DST.EMCalTrk.ll_edep_e>>hll_edep_e(300,-30,30)", "", "");
0205 hll_edep_e->SetTitle(
0206 Form(
0207 "EMCal + HCal_{in}: Electron likelihood distribution;log(Electron likelihood);A. U."));
0208
0209 p = (TPad *) c1->cd(idx++);
0210 c1->Update();
0211 p->SetLogy();
0212 T->Draw("DST.EMCalTrk.ll_edep_h>>hll_edep_h(300,-30,30)", "", "");
0213 hll_edep_h->SetTitle(
0214 Form(
0215 "EMCal + HCal_{in}: Hadron likelihood distribution;log(Hadron likelihood);A. U."));
0216
0217 p = (TPad *) c1->cd(idx++);
0218 c1->Update();
0219 p->SetLogy();
0220 T->Draw(
0221 "DST.EMCalTrk.ll_edep_e - DST.EMCalTrk.ll_edep_h>>hll_edep_diff(300,-30,30)",
0222 "", "");
0223 TH1F *hll_edep_diff = (TH1F *) gDirectory->GetObjectChecked("hll_edep_diff",
0224 "TH1F");
0225
0226 hll_edep_diff->SetTitle(
0227 Form(
0228 "EMCal + HCal_{in}: log likelihood cut;log(Electron likelihood) - log(Hadron likelihood);Count / bin"));
0229
0230
0231 p = (TPad *) c1->cd(idx++);
0232 c1->Update();
0233
0234
0235 p->DrawFrame(-30, 1e-4, 30, 1,
0236 "EMCal + HCal_{in}: Cut Efficiency;log(Electron likelihood) - log(Hadron likelihood);Cut Efficiency");
0237 TGraphErrors * ge = Distribution2Efficiency(hll_edep_diff);
0238 ge->SetLineColor(kBlue + 2);
0239 ge->SetMarkerColor(kBlue + 21);
0240 ge->SetMarkerColor(kFullCircle);
0241 ge->SetLineWidth(3);
0242 ge->Draw("lp");
0243
0244 SaveCanvas(c1,
0245 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0246 + TString(c1->GetName()), kFALSE);
0247
0248 }
0249
0250 void
0251 Edep_Distribution(TString infile)
0252 {
0253
0254 double N_Event = T->GetEntries();
0255
0256 TCanvas *c1 = new TCanvas("Edep_Distribution" + cuts,
0257 "Edep_Distribution" + cuts, 1900, 900);
0258 c1->Divide(2, 1);
0259 int idx = 1;
0260 TPad * p;
0261
0262 p = (TPad *) c1->cd(idx++);
0263 c1->Update();
0264 p->SetLogz();
0265 T->Draw(
0266 "DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.get_ep()>>h2_Edep_Distribution_raw(240,-.0,2, 240,-.0,12)",
0267 "", "colz");
0268 h2_Edep_Distribution_raw->SetTitle(
0269 Form(
0270 "Energy distribution;CEMC Cluster Energy in GeV;HCal_{in} Cluster Energy in GeV"));
0271 h2_Edep_Distribution_raw->Scale(1. / N_Event);
0272
0273
0274 p = (TPad *) c1->cd(idx++);
0275 c1->Update();
0276 p->SetLogz();
0277
0278 TH2F * h2_Edep_Distribution = (TH2F *) h2_Edep_Distribution_raw->Clone(
0279 "h2_Edep_Distribution");
0280
0281 h2_Edep_Distribution->Smooth(1, "k5b");
0282 h2_Edep_Distribution->Scale(1. / h2_Edep_Distribution->GetSum());
0283 h2_Edep_Distribution->Draw("colz");
0284
0285 SaveCanvas(c1,
0286 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0287 + TString(c1->GetName()), kFALSE);
0288
0289 }
0290
0291 void
0292 ShowerShape_Checks(TString infile, TCut good_track_cut)
0293 {
0294
0295 double N_Event = T->GetEntries();
0296
0297 TCanvas *c1 = new TCanvas("ShowerShape_Checks" + cuts,
0298 "ShowerShape_Checks" + cuts, 1900, 900);
0299 c1->Divide(2, 1);
0300 int idx = 1;
0301 TPad * p;
0302
0303 p = (TPad *) c1->cd(idx++);
0304 c1->Update();
0305 p->SetLogz();
0306
0307 T->Draw(
0308 "DST.EMCalTrk.cemc_energy/DST.EMCalTrk.cemc_sum_energy:DST.EMCalTrk.cemc_radius>>hEMCalTrk_cemc_shape(30,0,2,100,0,1)",
0309 good_track_cut, "goff");
0310
0311 TH2 * hEMCalTrk_cemc_shape = (TH2 *) gDirectory->GetObjectChecked(
0312 "hEMCalTrk_cemc_shape", "TH2");
0313 hEMCalTrk_cemc_shape->SetTitle(
0314 Form(
0315 "CEMC Shower Shape;Distance to track projection (Cluster width);Tower Energy/Cluster Energy"));
0316
0317
0318 TH1D * hEMCalTrk_cemc_shape_px = hEMCalTrk_cemc_shape->ProjectionX(
0319 "hEMCalTrk_cemc_shape_px");
0320
0321 for (int r = 1; r <= hEMCalTrk_cemc_shape->GetNbinsX(); r++)
0322 for (int e = 1; e <= hEMCalTrk_cemc_shape->GetNbinsY(); e++)
0323 {
0324 hEMCalTrk_cemc_shape->SetBinContent(r, e,
0325 hEMCalTrk_cemc_shape->GetBinContent(r, e)
0326 / hEMCalTrk_cemc_shape_px->GetBinContent(r));
0327 }
0328
0329 hEMCalTrk_cemc_shape->Draw("colz");
0330
0331
0332 p = (TPad *) c1->cd(idx++);
0333 c1->Update();
0334 p->SetLogz();
0335
0336 T->Draw(
0337 "DST.EMCalTrk.hcalin_energy/DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.hcalin_radius>>hEMCalTrk_hcalin_shape(30,0,2,100,0,1)",
0338 good_track_cut, "colz");
0339 TH2 * hEMCalTrk_hcalin_shape = (TH2 *) gDirectory->GetObjectChecked(
0340 "hEMCalTrk_hcalin_shape", "TH2");
0341
0342 hEMCalTrk_hcalin_shape->SetTitle(
0343 Form(
0344 "HCal_{in} Shower Shape;Distance to track projection (Cluster width);Tower Energy/Cluster Energy"));
0345
0346
0347 TH1D * hEMCalTrk_hcalin_shape_px = hEMCalTrk_hcalin_shape->ProjectionX(
0348 "hEMCalTrk_hcalin_shape_px");
0349
0350 for (int r = 1; r <= hEMCalTrk_hcalin_shape->GetNbinsX(); r++)
0351 for (int e = 1; e <= hEMCalTrk_hcalin_shape->GetNbinsY(); e++)
0352 {
0353 hEMCalTrk_hcalin_shape->SetBinContent(r, e,
0354 hEMCalTrk_hcalin_shape->GetBinContent(r, e)
0355 / hEMCalTrk_hcalin_shape_px->GetBinContent(r));
0356 }
0357
0358 hEMCalTrk_hcalin_shape->Draw("colz");
0359
0360 SaveCanvas(c1,
0361 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0362 + TString(c1->GetName()), kFALSE);
0363
0364 }
0365
0366 void
0367 Edep_Checks(TString infile, TCut good_track_cut)
0368 {
0369
0370 double N_Event = T->GetEntries();
0371
0372 TCanvas *c1 = new TCanvas("Edep_Checks" + cuts, "Edep_Checks" + cuts, 1900,
0373 950);
0374 c1->Divide(3, 2);
0375 int idx = 1;
0376 TPad * p;
0377
0378 p = (TPad *) c1->cd(idx++);
0379 c1->Update();
0380 p->SetLogy();
0381 T->Draw("DST.EMCalTrk.cemc_sum_n_tower>>hEMCalTrk_cemc_ntower(16,-.5,15.5)",
0382 good_track_cut);
0383 hEMCalTrk_cemc_ntower->SetTitle(
0384 Form("CEMC Cluster Size;Cluster Size (Towers);Probability"));
0385 hEMCalTrk_cemc_ntower->Scale(1. / N_Event);
0386
0387 p = (TPad *) c1->cd(idx++);
0388 c1->Update();
0389 p->SetLogy();
0390 T->Draw("DST.EMCalTrk.cemc_sum_energy>>hEMCalTrk_cemc_e(240,-.0,12)",
0391 good_track_cut);
0392 hEMCalTrk_cemc_e->SetTitle(
0393 Form("CEMC Cluster Energy;Cluster Energy (GeV);Count/bin"));
0394
0395
0396 p = (TPad *) c1->cd(idx++);
0397 c1->Update();
0398 p->SetLogy();
0399 p->DrawFrame(-.0, 1e-3, 12, 1,
0400 "CEMC Cut Eff;Cut on Cluster Energy (GeV);Efficiency");
0401
0402 TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_cemc_e);
0403 ge->SetLineColor(kBlue + 2);
0404 ge->SetMarkerColor(kBlue + 21);
0405 ge->SetMarkerColor(kFullCircle);
0406 ge->SetLineWidth(3);
0407 ge->Draw("lp");
0408
0409 p = (TPad *) c1->cd(idx++);
0410 c1->Update();
0411 p->SetLogy();
0412 T->Draw(
0413 "DST.EMCalTrk.hcalin_sum_n_tower>>hEMCalTrk_hcalin_ntower(16,-.5,15.5)",
0414 good_track_cut);
0415 hEMCalTrk_hcalin_ntower->SetTitle(
0416 Form("HCal_{in} Cluster Size;Cluster Size (Towers);Probability"));
0417 hEMCalTrk_hcalin_ntower->Scale(1. / N_Event);
0418
0419 p = (TPad *) c1->cd(idx++);
0420 c1->Update();
0421 p->SetLogy();
0422 T->Draw("DST.EMCalTrk.hcalin_sum_energy>>hEMCalTrk_hcalin_e(240,-.0,12)",
0423 good_track_cut);
0424 hEMCalTrk_hcalin_e->SetTitle(
0425 Form("HCal_{in} Cluster Energy;Cluster Energy (GeV);Count/bin"));
0426
0427
0428 p = (TPad *) c1->cd(idx++);
0429 c1->Update();
0430 p->SetLogy();
0431 p->DrawFrame(-.0, 1e-3, 12, 1,
0432 "HCal_{in} Cut Eff;Cut on Cluster Energy (GeV);Efficiency");
0433
0434 TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_hcalin_e);
0435 ge->SetLineColor(kBlue + 2);
0436 ge->SetMarkerColor(kBlue + 21);
0437 ge->SetMarkerColor(kFullCircle);
0438 ge->SetLineWidth(3);
0439 ge->Draw("lp");
0440
0441 SaveCanvas(c1,
0442 TString(_file0->GetName()) + TString("_DrawEcal_pDST_")
0443 + TString(c1->GetName()), kFALSE);
0444
0445 TCanvas *c1 = new TCanvas("Edep_Checks_2D" + cuts, "Edep_Checks_2D" + cuts,
0446 900, 900);
0447
0448
0449
0450
0451 p = (TPad *) c1->cd(idx++);
0452 c1->Update();
0453 p->SetLogz();
0454 T->Draw(
0455 "DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.cemc_sum_energy>>h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e(240,-.0,8, 240,-.0,8)",
0456 good_track_cut, "colz");
0457 h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->SetTitle(
0458 Form(
0459 "Energy distribution;CEMC Cluster Energy in GeV;HCal_{in} Cluster Energy in GeV"));
0460 h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->Scale(1. / N_Event);
0461 h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->GetZaxis()->SetRangeUser(1. / N_Event,
0462 1);
0463
0464 SaveCanvas(c1,
0465 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0466 + TString(c1->GetName()), kFALSE);
0467
0468 }
0469
0470 void
0471 Ep_Checks(TString infile, TCut good_track_cut)
0472 {
0473
0474 double N_Event = T->GetEntries();
0475
0476 TCanvas *c1 = new TCanvas("Ep_Checks" + cuts, "Ep_Checks" + cuts, 1900, 950);
0477 c1->Divide(2, 1);
0478 int idx = 1;
0479 TPad * p;
0480
0481 p = (TPad *) c1->cd(idx++);
0482 c1->Update();
0483 p->SetLogy();
0484 T->Draw("DST.EMCalTrk.get_ep()>>hEMCalTrk_get_ep(240,-.0,2)",
0485 good_track_cut);
0486 hEMCalTrk_get_ep->SetTitle(
0487 Form("CEMC Cluster Energy/Track Momentum;E/p;Count/bin"));
0488
0489
0490 p = (TPad *) c1->cd(idx++);
0491 c1->Update();
0492
0493 if (infile.Contains("e-") || infile.Contains("e+"))
0494 {
0495 p->DrawFrame(-.0, 0.8, 1.5, 1,
0496 "CEMC E/p Cut Eff;Cut on E/p;Signal Efficiency");
0497 }
0498 else
0499 {
0500 p->DrawFrame(-.0, 1e-3, 1.5, 1,
0501 "CEMC E/p Cut Eff;Cut on E/p;Background Efficiency or 1/Rejection");
0502 p->SetLogy();
0503 }
0504 TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_get_ep);
0505 ge->SetLineColor(kBlue + 2);
0506 ge->SetMarkerColor(kBlue + 21);
0507 ge->SetMarkerColor(kFullCircle);
0508 ge->SetLineWidth(3);
0509 ge->Draw("lp");
0510
0511 SaveCanvas(c1,
0512 TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
0513 + TString(c1->GetName()), kFALSE);
0514
0515 }
0516
0517 TGraphErrors *
0518 Distribution2Efficiency(TH1F * hCEMC3_Max)
0519 {
0520 double threshold[10000] =
0521 { 0 };
0522 double eff[10000] =
0523 { 0 };
0524 double eff_err[10000] =
0525 { 0 };
0526
0527 assert(hCEMC3_Max->GetNbinsX()<10000);
0528
0529 const double n = hCEMC3_Max->GetSum();
0530 double pass = 0;
0531 int cnt = 0;
0532 for (int i = hCEMC3_Max->GetNbinsX(); i >= 1; i--)
0533 {
0534 pass += hCEMC3_Max->GetBinContent(i);
0535
0536 const double pp = pass / n;
0537
0538 const double z = 1.;
0539
0540 const double A = z * sqrt(1. / n * pp * (1 - pp) + z * z / 4 / n / n);
0541 const double B = 1 / (1 + z * z / n);
0542
0543 threshold[cnt] = hCEMC3_Max->GetBinCenter(i);
0544 eff[cnt] = (pp + z * z / 2 / n) * B;
0545 eff_err[cnt] = A * B;
0546
0547
0548
0549 cnt++;
0550 }
0551 TGraphErrors * ge = new TGraphErrors(cnt, threshold, eff, NULL, eff_err);
0552
0553 return ge;
0554 }
