3 #include <qautils/QAHistManagerDef.h>
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawClusterContainer.h>
13 #include <calobase/RawTower.h>
14 #include <calobase/RawTowerContainer.h>
15 #include <calobase/RawTowerGeomContainer.h>
36 #include <CLHEP/Vector/ThreeVector.h>
48 :
SubsysReco(
"QAG4SimulationCalorimeter_" + calo_name)
49 , _calo_name(calo_name)
51 , _calo_hit_container(nullptr)
52 , _calo_abs_hit_container(nullptr)
53 , _truth_container(nullptr)
69 std::cout <<
"QAG4SimulationCalorimeter::InitRun - Fatal Error - "
70 <<
"unable to find DST node "
79 std::cout <<
"QAG4SimulationCalorimeter::InitRun - Fatal Error - "
80 <<
"unable to find DST node "
91 std::cout <<
"QAG4SimulationCalorimeter::InitRun - Fatal Error - "
92 <<
"unable to find DST node "
93 <<
"G4TruthInfo" << std::endl;
118 TString(
_calo_name) +
" Normalization;Items;Count", 10, .5, 10.5);
120 h->GetXaxis()->SetBinLabel(i++,
"Event");
121 h->GetXaxis()->SetBinLabel(i++,
"G4Hit Active");
122 h->GetXaxis()->SetBinLabel(i++,
"G4Hit Absor.");
123 h->GetXaxis()->SetBinLabel(i++,
"Tower");
124 h->GetXaxis()->SetBinLabel(i++,
"Tower Hit");
125 h->GetXaxis()->SetBinLabel(i++,
"Cluster");
126 h->GetXaxis()->LabelsOption(
"v");
132 std::cout <<
"QAG4SimulationCalorimeter::Init - Process sampling fraction"
139 std::cout <<
"QAG4SimulationCalorimeter::Init - Process tower occupancies"
146 std::cout <<
"QAG4SimulationCalorimeter::Init - Process tower occupancies"
157 std::cout <<
"QAG4SimulationCalorimeter::process_event() entered" << std::endl;
189 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
192 h_norm->Fill(
"Event", 1);
210 TString(
_calo_name) +
" RZ projection;G4 Hit Z (cm);G4 Hit R (cm)", 1200, -300, 300,
215 TString(
_calo_name) +
" XY projection;G4 Hit X (cm);G4 Hit Y (cm)", 1200, -300, 300,
220 TString(
_calo_name) +
" shower lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
221 200, -30, 30, 200, -30, 30));
224 TString(
_calo_name) +
" sampling fraction;Sampling fraction", 1000, 0, .2));
228 TString(
_calo_name) +
" visible sampling fraction;Visible sampling fraction", 1000, 0,
233 TString(
_calo_name) +
" hit time (edep weighting);Hit time - T0 (ns);Geant4 energy density",
240 TString(
_calo_name) +
" fraction truth energy ;G4 edep / particle energy",
245 TString(
_calo_name) +
" fraction visible energy from EM; visible energy from e^{#pm} / total visible energy",
254 std::cout <<
"QAG4SimulationCalorimeter::process_event_G4Hit() entered" << std::endl;
261 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
269 double total_primary_energy = 1
e-9;
271 particle_iter != primary_range.second; ++particle_iter)
275 total_primary_energy += particle->
get_e();
286 <<
"QAG4SimulationCalorimeter::process_event_G4Hit() handle this truth particle"
288 last_primary->identify();
296 <<
"QAG4SimulationCalorimeter::process_event_G4Hit() handle this vertex"
301 const double t0 = primary_vtx->
get_t();
303 primary_vtx->
get_z());
306 TVector3 axis_proj(last_primary->get_px(), last_primary->get_py(),
307 last_primary->get_pz());
308 if (axis_proj.Mag() == 0)
309 axis_proj.SetXYZ(0, 0, 1);
310 axis_proj = axis_proj.Unit();
313 TVector3 axis_azimuth = axis_proj.Cross(TVector3(0, 0, 1));
314 if (axis_azimuth.Mag() == 0)
315 axis_azimuth.SetXYZ(1, 0, 0);
316 axis_azimuth = axis_azimuth.Unit();
319 TVector3 axis_polar = axis_proj.Cross(axis_azimuth);
320 assert(axis_polar.Mag() > 0);
321 axis_polar = axis_polar.Unit();
324 double ev_calo = 0.0;
325 double ea_calo = 0.0;
326 double ev_calo_em = 0.0;
330 TH2F *hrz =
dynamic_cast<TH2F *
>(hm->
getHisto(
333 TH2F *hxy =
dynamic_cast<TH2F *
>(hm->
getHisto(
336 TH1F *ht =
dynamic_cast<TH1F *
>(hm->
getHisto(
339 TH2F *hlat =
dynamic_cast<TH2F *
>(hm->
getHisto(
347 hit_iter != calo_hit_range.second; hit_iter++)
349 PHG4Hit *this_hit = hit_iter->second;
360 std::cout << __PRETTY_FUNCTION__ <<
" - Error - this PHG4hit missing particle: ";
364 if (abs(particle->
get_pid()) == 11)
370 hrz->Fill(hit.Z(), hit.Perp(), this_hit->
get_edep());
371 hxy->Fill(hit.X(), hit.Y(), this_hit->
get_edep());
374 const double hit_azimuth = axis_azimuth.Dot(hit -
vertex);
375 const double hit_polar = axis_polar.Dot(hit -
vertex);
376 hlat->Fill(hit_polar, hit_azimuth, this_hit->
get_edep());
387 hit_iter != calo_abs_hit_range.second; hit_iter++)
389 PHG4Hit *this_hit = hit_iter->second;
397 std::cout <<
"QAG4SimulationCalorimeter::process_event_G4Hit::" <<
_calo_name
398 <<
" - SF = " << e_calo / (e_calo + ea_calo + 1
e-9) <<
", VSF = "
399 << ev_calo / (e_calo + ea_calo + 1
e-9) << std::endl;
401 if (e_calo + ea_calo > 0)
405 h->Fill(e_calo / (e_calo + ea_calo));
409 h->Fill(ev_calo / (e_calo + ea_calo));
412 h =
dynamic_cast<TH1F *
>(hm->
getHisto(
415 h->Fill((e_calo + ea_calo) / total_primary_energy);
419 h =
dynamic_cast<TH1F *
>(hm->
getHisto(
422 h->Fill(ev_calo_em / (ev_calo));
426 std::cout <<
"QAG4SimulationCalorimeter::process_event_G4Hit::" <<
_calo_name
427 <<
" - histogram " << h->GetName() <<
" Get Sum = " << h->GetSum()
439 TString(
_calo_name) +
" 1x1 tower;1x1 TOWER Energy (GeV)", 100, 9
e-4, 100);
445 TString(
_calo_name) +
" 1x1 tower max per event;1x1 tower max per event (GeV)", 5000,
449 TString(
_calo_name) +
" 2x2 tower;2x2 TOWER Energy (GeV)", 100, 9
e-4, 100);
454 TString(
_calo_name) +
" 2x2 tower max per event;2x2 tower max per event (GeV)", 5000,
458 TString(
_calo_name) +
" 3x3 tower;3x3 TOWER Energy (GeV)", 100, 9
e-4, 100);
463 TString(
_calo_name) +
" 3x3 tower max per event;3x3 tower max per event (GeV)", 5000,
467 TString(
_calo_name) +
" 4x4 tower;4x4 TOWER Energy (GeV)", 100, 9
e-4, 100);
472 TString(
_calo_name) +
" 4x4 tower max per event;4x4 tower max per event (GeV)", 5000,
476 TString(
_calo_name) +
" 5x5 tower;5x5 TOWER Energy (GeV)", 100, 9
e-4, 100);
481 TString(
_calo_name) +
" 5x5 tower max per event;5x5 tower max per event (GeV)", 5000,
492 std::cout <<
"QAG4SimulationCalorimeter::process_event_Tower() entered" << std::endl;
496 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
506 std::cout <<
PHWHERE <<
": Could not find node " << towernodename
512 topNode, towergeomnodename);
515 std::cout <<
PHWHERE <<
": Could not find node " << towergeomnodename
520 static const int max_size = 5;
521 std::map<int, std::string> size_label;
522 size_label[1] =
"1x1";
523 size_label[2] =
"2x2";
524 size_label[3] =
"3x3";
525 size_label[4] =
"4x4";
526 size_label[5] =
"5x5";
527 std::map<int, double> max_energy;
528 std::map<int, TH1F *> energy_hist_list;
529 std::map<int, TH1F *> max_energy_hist_list;
533 max_energy[
size] = 0;
535 TH1F *
h =
dynamic_cast<TH1F *
>(hm->
getHisto(
538 energy_hist_list[
size] =
h;
539 h =
dynamic_cast<TH1F *
>(hm->
getHisto(
542 max_energy_hist_list[
size] =
h;
545 h_norm->Fill(
"Tower", towergeom->
size());
546 h_norm->Fill(
"Tower Hit", towers->
size());
548 for (
int binphi = 0; binphi < towergeom->
get_phibins(); ++binphi)
550 for (
int bineta = 0; bineta < towergeom->
get_etabins(); ++bineta)
555 if ((
size == 2 or
size == 4) and ((binphi % 2 != 0) and (bineta % 2 != 0)))
561 for (
int iphi = binphi; iphi < binphi +
size; ++iphi)
563 for (
int ieta = bineta; ieta < bineta +
size; ++ieta)
587 energy_hist_list[
size]->Fill(energy == 0 ? 9.1
e-4 : energy);
589 if (energy > max_energy[size])
598 max_energy_hist_list[
size]->Fill(max_energy[
size]);
610 TString(
_calo_name) +
" best matched cluster E/E_{Truth};E_{Cluster}/E_{Truth}", 150,
615 TString(
_calo_name) +
" best cluster lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
616 200, -15, 15, 200, -15, 15));
625 std::cout <<
"QAG4SimulationCalorimeter::process_event_Cluster() entered"
633 topNode, towergeomnodename);
636 std::cout <<
PHWHERE <<
": Could not find node " << towergeomnodename << std::endl;
647 h_norm->Fill(
"Cluster", clusters->
size());
658 <<
"QAG4SimulationCalorimeter::process_event_Cluster() handle this truth particle"
660 last_primary->identify();
667 TH1F *
h =
dynamic_cast<TH1F *
>(hm->
getHisto(
677 std::cout <<
"QAG4SimulationCalorimeter::process_event_Cluster::"
678 <<
_calo_name <<
" - get cluster with energy "
679 << cluster->
get_energy() <<
" VS primary energy "
680 << last_primary->get_e() << std::endl;
682 h->Fill(cluster->
get_energy() / (last_primary->get_e() + 1
e-9));
693 <<
"QAG4SimulationCalorimeter::process_event_Cluster() handle this vertex"
695 primary_vtx->identify();
698 const CLHEP::Hep3Vector
vertex(primary_vtx->get_x(), primary_vtx->get_y(),
699 primary_vtx->get_z());
702 CLHEP::Hep3Vector axis_proj(last_primary->get_px(), last_primary->get_py(),
703 last_primary->get_pz());
704 if (axis_proj.mag() == 0)
705 axis_proj.set(0, 0, 1);
706 axis_proj = axis_proj.unit();
709 CLHEP::Hep3Vector axis_azimuth = axis_proj.cross(CLHEP::Hep3Vector(0, 0, 1));
710 if (axis_azimuth.mag() == 0)
711 axis_azimuth.set(1, 0, 0);
712 axis_azimuth = axis_azimuth.unit();
715 CLHEP::Hep3Vector axis_polar = axis_proj.cross(axis_azimuth);
716 assert(axis_polar.mag() > 0);
717 axis_polar = axis_polar.unit();
719 TH2F *hlat =
dynamic_cast<TH2F *
>(hm->
getHisto(
723 const double hit_azimuth = axis_azimuth.dot(hit -
vertex);
724 const double hit_polar = axis_polar.dot(hit -
vertex);
725 hlat->Fill(hit_polar, hit_azimuth);
730 std::cout <<
"QAG4SimulationCalorimeter::process_event_Cluster::"