2 #include <qautils/QAHistManagerDef.h>
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawTower.h>
13 #include <calobase/RawTowerContainer.h>
14 #include <calobase/RawTowerGeomContainer.h>
43 , _calo_name_cemc(
"CEMC")
44 , _calo_name_hcalin(
"HCALIN")
45 , _calo_name_hcalout(
"HCALOUT")
46 , _truth_container(nullptr)
57 std::cout <<
"QAG4SimulationCalorimeterSum::InitRun - Fatal Error - "
58 <<
"unable to find DST node "
59 <<
"G4TruthInfo" << std::endl;
107 h->GetXaxis()->SetBinLabel(i++,
"Event");
111 h->GetXaxis()->SetBinLabel(i++, (
_calo_name_cemc +
" Cluster").c_str());
114 h->GetXaxis()->SetBinLabel(i++,
"Track");
115 h->GetXaxis()->LabelsOption(
"v");
128 std::cout <<
"QAG4SimulationCalorimeterSum::Init - Process tower occupancies"
136 std::cout <<
"QAG4SimulationCalorimeterSum::Init - Process sampling fraction"
146 std::cout <<
"QAG4SimulationCalorimeterSum::process_event() entered" << std::endl;
184 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
187 h_norm->Fill(
"Event", 1);
195 return "h_QAG4Sim_CalorimeterSum_";
218 TString(
_calo_name_cemc.c_str()) +
" Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
225 TString(
_calo_name_hcalin.c_str()) +
" Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
232 TString(
_calo_name_hcalout.c_str()) +
" Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
238 "Tower 3x3 sum /E_{Truth};#Sigma_{3x3}[E_{Tower}] / total truth energy",
243 "Tower 5x5 sum /E_{Truth};#Sigma_{5x5}[E_{Tower}] / total truth energy",
252 std::cout <<
"QAG4SimulationCalorimeterSum::process_event_TrackProj() entered"
267 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
270 h_norm->Fill(
"Track", 1);
273 TH1F *hsum =
dynamic_cast<TH1F *
>(hm->
getHisto(
281 TH1F *hsum =
dynamic_cast<TH1F *
>(hm->
getHisto(
306 TH2F *h2_proj =
dynamic_cast<TH2F *
>(hm->
getHisto(
313 topNode, towergeonodename);
328 std::cout << __PRETTY_FUNCTION__ <<
" - info - handling track track p = ("
330 << track->
get_pz() <<
"), v = (" << track->
get_x() <<
", "
331 << track->
get_z() <<
", " << track->
get_z() <<
")"
339 std::vector<double>
point;
340 point.assign(3, NAN);
346 if (std::isnan(point[0]) or std::isnan(point[1]) or std::isnan(point[2]))
354 assert(not std::isnan(point[0]));
355 assert(not std::isnan(point[1]));
356 assert(not std::isnan(point[2]));
362 double phi = atan2(y, x);
363 double eta = asinh(z / sqrt(x * x + y * y));
370 if (bineta > 1 and bineta < towergeo->get_etabins() - 1)
378 const double etabin_shift = (towergeo->
get_etacenter(bineta) -
eta) / etabin_width;
379 const double phibin_shift = (fmod(
380 towergeo->
get_phicenter(binphi) - phi + 5 * M_PI, 2 * M_PI) -
385 std::cout << __PRETTY_FUNCTION__ <<
" - info - handling track proj. (" << x
386 <<
", " << y <<
", " << z <<
")"
387 <<
", eta " << eta <<
", phi" << phi
388 <<
", bineta " << bineta <<
" - ["
390 << towergeo->
get_etabounds(bineta).second <<
"], etabin_shift = "
391 << etabin_shift <<
", binphi " << binphi <<
" - ["
393 << towergeo->
get_phibounds(binphi).second <<
"], phibin_shift = "
394 << phibin_shift << std::endl;
396 const int bin_search_range = (
Max_N_Tower - 1) / 2;
397 for (
int iphi = binphi - bin_search_range; iphi <= binphi + bin_search_range;
400 for (
int ieta = bineta - bin_search_range;
401 ieta <= bineta + bin_search_range; ++ieta)
421 std::cout << __PRETTY_FUNCTION__ <<
" - info - processing tower"
422 <<
", bineta " << ieta <<
" - ["
426 << wrapphi <<
" - [" << towergeo->
get_phibounds(wrapphi).first
427 <<
", " << towergeo->
get_phibounds(wrapphi).second <<
"]" << std::endl;
434 std::cout << __PRETTY_FUNCTION__ <<
" - info - tower " << ieta <<
" "
435 << wrapphi <<
" energy = " << tower->
get_energy() << std::endl;
440 h2_proj->Fill(ieta - bineta + etabin_shift,
441 iphi - binphi + phibin_shift, energy);
479 70, 0, 70, 70, 0, 70));
485 70, 0, 70, 70, 0, 70));
489 "Total Cluster E_{Reco}/E_{Truth};Reco cluster energy sum / total truth energy",
508 std::cout <<
"QAG4SimulationCalorimeterSum::process_event_Cluster() entered"
525 const double cluster_cemc_e = cluster_cemc ? cluster_cemc->
get_energy() : 0;
526 const double cluster_hcalin_e =
527 cluster_hcalin ? cluster_hcalin->
get_energy() : 0;
528 const double cluster_hcalout_e =
529 cluster_hcalout ? cluster_hcalout->
get_energy() : 0;
531 if (cluster_cemc_e + cluster_hcalin_e > 0)
537 h2->Fill(cluster_cemc_e, cluster_hcalin_e);
539 TH1F *hr =
dynamic_cast<TH1F *
>(hm->
getHisto(
543 hr->Fill(cluster_cemc_e / (cluster_cemc_e + cluster_hcalin_e));
546 if (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e > 0)
550 std::cout <<
"QAG4SimulationCalorimeterSum::process_event_Cluster - "
551 <<
" cluster_cemc_e = " << cluster_cemc_e
552 <<
" cluster_hcalin_e = " << cluster_hcalin_e
553 <<
" cluster_hcalout_e = " << cluster_hcalout_e <<
" hr = "
554 << (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e)
562 h2->Fill((cluster_cemc_e + cluster_hcalin_e), cluster_hcalout_e);
564 TH1F *hr =
dynamic_cast<TH1F *
>(hm->
getHisto(
569 (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e));
571 TH1F *hsum =
dynamic_cast<TH1F *
>(hm->
getHisto(
576 (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e) / (primary->
get_e() + 1
e-9));