Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:22 

0001 #include "QAG4SimulationCalorimeter.h"
0002 
0003 #include <qautils/QAHistManagerDef.h>
0004 
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4TruthInfoContainer.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010 
0011 #include <calobase/RawCluster.h>
0012 #include <calobase/RawClusterContainer.h>
0013 #include <calobase/RawTower.h>
0014 #include <calobase/RawTowerContainer.h>
0015 #include <calobase/RawTowerGeomContainer.h>
0016 #include <calobase/TowerInfo.h>
0017 #include <calobase/TowerInfoContainer.h>
0018 #include <calobase/TowerInfoDefs.h>
0019 
0020 #include <g4eval/CaloEvalStack.h>
0021 #include <g4eval/CaloRawClusterEval.h>
0022 
0023 #include <fun4all/Fun4AllHistoManager.h>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/SubsysReco.h>
0026 
0027 #include <phool/PHCompositeNode.h>
0028 #include <phool/PHNodeIterator.h>
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>
0031 
0032 #include <TAxis.h>
0033 #include <TH1.h>
0034 #include <TH2.h>
0035 #include <TString.h>
0036 #include <TVector3.h>
0037 
0038 #include <CLHEP/Vector/ThreeVector.h>  // for Hep3Vector
0039 
0040 #include <algorithm>
0041 #include <cassert>
0042 #include <cmath>
0043 #include <cstdlib>
0044 #include <iostream>
0045 #include <iterator>  // for reverse_iterator
0046 #include <map>
0047 #include <string>
0048 #include <utility>
0049 
0050 QAG4SimulationCalorimeter::QAG4SimulationCalorimeter(const std::string &calo_name,
0051                                                      QAG4SimulationCalorimeter::enu_flags flags)
0052   : SubsysReco("QAG4SimulationCalorimeter_" + calo_name)
0053   , _calo_name(calo_name)
0054   , _flags(flags)
0055   , _calo_hit_container(nullptr)
0056   , _calo_abs_hit_container(nullptr)
0057   , _truth_container(nullptr)
0058 {
0059 }
0060 
0061 int QAG4SimulationCalorimeter::InitRun(PHCompositeNode *topNode)
0062 {
0063   PHNodeIterator iter(topNode);
0064   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0065   assert(dstNode);
0066 
0067   if (flag(kProcessG4Hit))
0068   {
0069     _calo_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0070                                                                "G4HIT_" + _calo_name);
0071     if (!_calo_hit_container)
0072     {
0073       std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
0074                 << "unable to find DST node "
0075                 << "G4HIT_" + _calo_name << std::endl;
0076       assert(_calo_hit_container);
0077     }
0078 
0079     _calo_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0080                                                                    "G4HIT_ABSORBER_" + _calo_name);
0081     if (!_calo_abs_hit_container)
0082     {
0083       std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
0084                 << "unable to find DST node "
0085                 << "G4HIT_ABSORBER_" + _calo_name
0086                 << std::endl;
0087       assert(_calo_abs_hit_container);
0088     }
0089   }
0090 
0091   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0092                                                                 "G4TruthInfo");
0093   if (!_truth_container)
0094   {
0095     std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
0096               << "unable to find DST node "
0097               << "G4TruthInfo" << std::endl;
0098     assert(_truth_container);
0099   }
0100 
0101   if (flag(kProcessCluster))
0102   {
0103     if (!_caloevalstack)
0104     {
0105       _caloevalstack.reset(new CaloEvalStack(topNode, _calo_name));
0106       _caloevalstack->set_strict(true);
0107       _caloevalstack->set_verbosity(Verbosity() + 1);
0108     }
0109     else
0110     {
0111       _caloevalstack->next_event(topNode);
0112     }
0113   }
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116 
0117 int QAG4SimulationCalorimeter::Init(PHCompositeNode *topNode)
0118 {
0119   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0120   assert(hm);
0121   TH1D *h = new TH1D(TString(get_histo_prefix()) + "_Normalization",  //
0122                      TString(_calo_name) + " Normalization;Items;Count", 10, .5, 10.5);
0123   int i = 1;
0124   h->GetXaxis()->SetBinLabel(i++, "Event");
0125   h->GetXaxis()->SetBinLabel(i++, "G4Hit Active");
0126   h->GetXaxis()->SetBinLabel(i++, "G4Hit Absor.");
0127   h->GetXaxis()->SetBinLabel(i++, "Tower");
0128   h->GetXaxis()->SetBinLabel(i++, "Tower Hit");
0129   h->GetXaxis()->SetBinLabel(i++, "Cluster");
0130   h->GetXaxis()->LabelsOption("v");
0131   hm->registerHisto(h);
0132 
0133   if (flag(kProcessG4Hit))
0134   {
0135     if (Verbosity() >= 1)
0136     {
0137       std::cout << "QAG4SimulationCalorimeter::Init - Process sampling fraction"
0138                 << std::endl;
0139     }
0140     Init_G4Hit(topNode);
0141   }
0142   if (flag(kProcessTower))
0143   {
0144     if (Verbosity() >= 1)
0145     {
0146       std::cout << "QAG4SimulationCalorimeter::Init - Process tower occupancies"
0147                 << std::endl;
0148     }
0149     Init_Tower(topNode);
0150   }
0151   if (flag(kProcessCluster))
0152   {
0153     if (Verbosity() >= 1)
0154     {
0155       std::cout << "QAG4SimulationCalorimeter::Init - Process tower occupancies"
0156                 << std::endl;
0157     }
0158     Init_Cluster(topNode);
0159   }
0160 
0161   return Fun4AllReturnCodes::EVENT_OK;
0162 }
0163 
0164 int QAG4SimulationCalorimeter::process_event(PHCompositeNode *topNode)
0165 {
0166   if (Verbosity() > 2)
0167   {
0168     std::cout << "QAG4SimulationCalorimeter::process_event() entered" << std::endl;
0169   }
0170 
0171   if (_caloevalstack)
0172   {
0173     _caloevalstack->next_event(topNode);
0174   }
0175 
0176   if (flag(kProcessG4Hit))
0177   {
0178     int const ret = process_event_G4Hit(topNode);
0179 
0180     if (ret != Fun4AllReturnCodes::EVENT_OK)
0181     {
0182       return ret;
0183     }
0184   }
0185 
0186   if (flag(kProcessTower))
0187   {
0188     int const ret = process_event_Tower(topNode);
0189 
0190     if (ret != Fun4AllReturnCodes::EVENT_OK)
0191     {
0192       return ret;
0193     }
0194   }
0195 
0196   if (flag(kProcessCluster))
0197   {
0198     int const ret = process_event_Cluster(topNode);
0199 
0200     if (ret != Fun4AllReturnCodes::EVENT_OK)
0201     {
0202       return ret;
0203     }
0204   }
0205 
0206   // at the end, count success events
0207   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0208   assert(hm);
0209   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0210       get_histo_prefix() + "_Normalization"));
0211   assert(h_norm);
0212   h_norm->Fill("Event", 1);
0213 
0214   return Fun4AllReturnCodes::EVENT_OK;
0215 }
0216 
0217 std::string
0218 QAG4SimulationCalorimeter::get_histo_prefix()
0219 {
0220   return "h_QAG4Sim_" + std::string(_calo_name);
0221 }
0222 
0223 int QAG4SimulationCalorimeter::Init_G4Hit(PHCompositeNode * /*topNode*/)
0224 {
0225   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0226   assert(hm);
0227 
0228   hm->registerHisto(
0229       new TH2F(TString(get_histo_prefix()) + "_G4Hit_RZ",  //
0230                TString(_calo_name) + " RZ projection;G4 Hit Z (cm);G4 Hit R (cm)", 1200, -300, 300,
0231                600, -000, 300));
0232 
0233   hm->registerHisto(
0234       new TH2F(TString(get_histo_prefix()) + "_G4Hit_XY",  //
0235                TString(_calo_name) + " XY projection;G4 Hit X (cm);G4 Hit Y (cm)", 1200, -300, 300,
0236                1200, -300, 300));
0237 
0238   hm->registerHisto(
0239       new TH2F(TString(get_histo_prefix()) + "_G4Hit_LateralTruthProjection",  //
0240                TString(_calo_name) + " shower lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
0241                200, -30, 30, 200, -30, 30));
0242 
0243   hm->registerHisto(new TH1F(TString(get_histo_prefix()) + "_G4Hit_SF",  //
0244                              TString(_calo_name) + " sampling fraction;Sampling fraction", 1000, 0, .2));
0245 
0246   hm->registerHisto(
0247       new TH1F(TString(get_histo_prefix()) + "_G4Hit_VSF",  //
0248                TString(_calo_name) + " visible sampling fraction;Visible sampling fraction", 1000, 0,
0249                .2));
0250 
0251   TH1F *h =
0252       new TH1F(TString(get_histo_prefix()) + "_G4Hit_HitTime",  //
0253                TString(_calo_name) + " hit time (edep weighting);Hit time - T0 (ns);Geant4 energy density",
0254                1000, 0.5, 10000);
0255   QAHistManagerDef::useLogBins(h->GetXaxis());
0256   hm->registerHisto(h);
0257 
0258   hm->registerHisto(
0259       new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionTruthEnergy",  //
0260                TString(_calo_name) + " fraction truth energy ;G4 edep / particle energy",
0261                1000, 0, 1));
0262 
0263   hm->registerHisto(
0264       new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionEMVisibleEnergy",  //
0265                TString(_calo_name) + " fraction visible energy from EM; visible energy from e^{#pm} / total visible energy",
0266                100, 0, 1));
0267 
0268   return Fun4AllReturnCodes::EVENT_OK;
0269 }
0270 
0271 int QAG4SimulationCalorimeter::process_event_G4Hit(PHCompositeNode * /*topNode*/)
0272 {
0273   if (Verbosity() > 2)
0274   {
0275     std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit() entered" << std::endl;
0276   }
0277 
0278   TH1F *h = nullptr;
0279 
0280   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0281   assert(hm);
0282 
0283   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0284       get_histo_prefix() + "_Normalization"));
0285   assert(h_norm);
0286 
0287   // get primary
0288   assert(_truth_container);
0289   PHG4TruthInfoContainer::ConstRange const primary_range =
0290       _truth_container->GetPrimaryParticleRange();
0291   double total_primary_energy = 1e-9;  // make it zero energy epsilon samll so it can be used for denominator
0292   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0293        particle_iter != primary_range.second; ++particle_iter)
0294   {
0295     const PHG4Particle *particle = particle_iter->second;
0296     assert(particle);
0297     total_primary_energy += particle->get_e();
0298   }
0299 
0300   assert(!_truth_container->GetMap().empty());
0301   const PHG4Particle *last_primary =
0302       _truth_container->GetMap().rbegin()->second;
0303   assert(last_primary);
0304 
0305   if (Verbosity() > 2)
0306   {
0307     std::cout
0308         << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this truth particle"
0309         << std::endl;
0310     last_primary->identify();
0311   }
0312   const PHG4VtxPoint *primary_vtx =  //
0313       _truth_container->GetPrimaryVtx(last_primary->get_vtx_id());
0314   assert(primary_vtx);
0315   if (Verbosity() > 2)
0316   {
0317     std::cout
0318         << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this vertex"
0319         << std::endl;
0320     primary_vtx->identify();
0321   }
0322 
0323   const double t0 = primary_vtx->get_t();
0324   const TVector3 vertex(primary_vtx->get_x(), primary_vtx->get_y(),
0325                         primary_vtx->get_z());
0326 
0327   // projection axis
0328   TVector3 axis_proj(last_primary->get_px(), last_primary->get_py(),
0329                      last_primary->get_pz());
0330   if (axis_proj.Mag() == 0)
0331   {
0332     axis_proj.SetXYZ(0, 0, 1);
0333   }
0334   axis_proj = axis_proj.Unit();
0335 
0336   // azimuthal direction axis
0337   TVector3 axis_azimuth = axis_proj.Cross(TVector3(0, 0, 1));
0338   if (axis_azimuth.Mag() == 0)
0339   {
0340     axis_azimuth.SetXYZ(1, 0, 0);
0341   }
0342   axis_azimuth = axis_azimuth.Unit();
0343 
0344   // polar direction axis
0345   TVector3 axis_polar = axis_proj.Cross(axis_azimuth);
0346   assert(axis_polar.Mag() > 0);
0347   axis_polar = axis_polar.Unit();
0348 
0349   double e_calo = 0.0;      // active energy deposition
0350   double ev_calo = 0.0;     // visible energy
0351   double ea_calo = 0.0;     // absorber energy
0352   double ev_calo_em = 0.0;  // EM visible energy
0353 
0354   if (_calo_hit_container)
0355   {
0356     TH2F *hrz = dynamic_cast<TH2F *>(hm->getHisto(
0357         get_histo_prefix() + "_G4Hit_RZ"));
0358     assert(hrz);
0359     TH2F *hxy = dynamic_cast<TH2F *>(hm->getHisto(
0360         get_histo_prefix() + "_G4Hit_XY"));
0361     assert(hxy);
0362     TH1F *ht = dynamic_cast<TH1F *>(hm->getHisto(
0363         get_histo_prefix() + "_G4Hit_HitTime"));
0364     assert(ht);
0365     TH2F *hlat = dynamic_cast<TH2F *>(hm->getHisto(
0366         get_histo_prefix() + "_G4Hit_LateralTruthProjection"));
0367     assert(hlat);
0368 
0369     h_norm->Fill("G4Hit Active", _calo_hit_container->size());
0370     PHG4HitContainer::ConstRange const calo_hit_range =
0371         _calo_hit_container->getHits();
0372     for (PHG4HitContainer::ConstIterator hit_iter = calo_hit_range.first;
0373          hit_iter != calo_hit_range.second; hit_iter++)
0374     {
0375       PHG4Hit *this_hit = hit_iter->second;
0376       assert(this_hit);
0377 
0378       e_calo += this_hit->get_edep();
0379       ev_calo += this_hit->get_light_yield();
0380 
0381       // EM visible energy that is only associated with electron energy deposition
0382       PHG4Particle *particle = _truth_container->GetParticle(
0383           this_hit->get_trkid());
0384       if (!particle)
0385       {
0386         std::cout << __PRETTY_FUNCTION__ << " - Error - this PHG4hit missing particle: ";
0387         this_hit->identify();
0388       }
0389       assert(particle);
0390       if (abs(particle->get_pid()) == 11)
0391       {
0392         ev_calo_em += this_hit->get_light_yield();
0393       }
0394 
0395       const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0396                          this_hit->get_avg_z());
0397 
0398       hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_edep());
0399       hxy->Fill(hit.X(), hit.Y(), this_hit->get_edep());
0400       ht->Fill(this_hit->get_avg_t() - t0, this_hit->get_edep());
0401 
0402       const double hit_azimuth = axis_azimuth.Dot(hit - vertex);
0403       const double hit_polar = axis_polar.Dot(hit - vertex);
0404       hlat->Fill(hit_polar, hit_azimuth, this_hit->get_edep());
0405     }
0406   }
0407 
0408   if (_calo_abs_hit_container)
0409   {
0410     h_norm->Fill("G4Hit Absor.", _calo_abs_hit_container->size());
0411 
0412     PHG4HitContainer::ConstRange const calo_abs_hit_range =
0413         _calo_abs_hit_container->getHits();
0414     for (PHG4HitContainer::ConstIterator hit_iter = calo_abs_hit_range.first;
0415          hit_iter != calo_abs_hit_range.second; hit_iter++)
0416     {
0417       PHG4Hit *this_hit = hit_iter->second;
0418       assert(this_hit);
0419 
0420       ea_calo += this_hit->get_edep();
0421     }
0422   }
0423 
0424   if (Verbosity() > 3)
0425   {
0426     std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
0427               << " - SF = " << e_calo / (e_calo + ea_calo + 1e-9) << ", VSF = "
0428               << ev_calo / (e_calo + ea_calo + 1e-9) << std::endl;
0429   }
0430 
0431   if (e_calo + ea_calo > 0)
0432   {
0433     h = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "_G4Hit_SF"));
0434     assert(h);
0435     h->Fill(e_calo / (e_calo + ea_calo));
0436 
0437     h = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "_G4Hit_VSF"));
0438     assert(h);
0439     h->Fill(ev_calo / (e_calo + ea_calo));
0440   }
0441 
0442   h = dynamic_cast<TH1F *>(hm->getHisto(
0443       get_histo_prefix() + "_G4Hit_FractionTruthEnergy"));
0444   assert(h);
0445   h->Fill((e_calo + ea_calo) / total_primary_energy);
0446 
0447   if (ev_calo > 0)
0448   {
0449     h = dynamic_cast<TH1F *>(hm->getHisto(
0450         get_histo_prefix() + "_G4Hit_FractionEMVisibleEnergy"));
0451     assert(h);
0452     h->Fill(ev_calo_em / (ev_calo));
0453   }
0454 
0455   if (Verbosity() > 3)
0456   {
0457     std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
0458               << " - histogram " << h->GetName() << " Get Sum = " << h->GetSum()
0459               << std::endl;
0460   }
0461 
0462   return Fun4AllReturnCodes::EVENT_OK;
0463 }
0464 
0465 int QAG4SimulationCalorimeter::Init_Tower(PHCompositeNode * /*topNode*/)
0466 {
0467   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0468   assert(hm);
0469 
0470   TH1F *h = new TH1F(TString(get_histo_prefix()) + "_Tower_1x1",  //
0471                      TString(_calo_name) + " 1x1 tower;1x1 TOWER Energy (GeV)", 100, 9e-4, 100);
0472   QAHistManagerDef::useLogBins(h->GetXaxis());
0473   hm->registerHisto(h);
0474 
0475   hm->registerHisto(
0476       new TH1F(TString(get_histo_prefix()) + "_Tower_1x1_max",  //
0477                TString(_calo_name) + " 1x1 tower max per event;1x1 tower max per event (GeV)", 5000,
0478                0, 50));
0479 
0480   h = new TH1F(TString(get_histo_prefix()) + "_Tower_2x2",  //
0481                TString(_calo_name) + " 2x2 tower;2x2 TOWER Energy (GeV)", 100, 9e-4, 100);
0482   QAHistManagerDef::useLogBins(h->GetXaxis());
0483   hm->registerHisto(h);
0484   hm->registerHisto(
0485       new TH1F(TString(get_histo_prefix()) + "_Tower_2x2_max",  //
0486                TString(_calo_name) + " 2x2 tower max per event;2x2 tower max per event (GeV)", 5000,
0487                0, 50));
0488 
0489   h = new TH1F(TString(get_histo_prefix()) + "_Tower_3x3",  //
0490                TString(_calo_name) + " 3x3 tower;3x3 TOWER Energy (GeV)", 100, 9e-4, 100);
0491   QAHistManagerDef::useLogBins(h->GetXaxis());
0492   hm->registerHisto(h);
0493   hm->registerHisto(
0494       new TH1F(TString(get_histo_prefix()) + "_Tower_3x3_max",  //
0495                TString(_calo_name) + " 3x3 tower max per event;3x3 tower max per event (GeV)", 5000,
0496                0, 50));
0497 
0498   h = new TH1F(TString(get_histo_prefix()) + "_Tower_4x4",  //
0499                TString(_calo_name) + " 4x4 tower;4x4 TOWER Energy (GeV)", 100, 9e-4, 100);
0500   QAHistManagerDef::useLogBins(h->GetXaxis());
0501   hm->registerHisto(h);
0502   hm->registerHisto(
0503       new TH1F(TString(get_histo_prefix()) + "_Tower_4x4_max",  //
0504                TString(_calo_name) + " 4x4 tower max per event;4x4 tower max per event (GeV)", 5000,
0505                0, 50));
0506 
0507   h = new TH1F(TString(get_histo_prefix()) + "_Tower_5x5",  //
0508                TString(_calo_name) + " 5x5 tower;5x5 TOWER Energy (GeV)", 100, 9e-4, 100);
0509   QAHistManagerDef::useLogBins(h->GetXaxis());
0510   hm->registerHisto(h);
0511   hm->registerHisto(
0512       new TH1F(TString(get_histo_prefix()) + "_Tower_5x5_max",  //
0513                TString(_calo_name) + " 5x5 tower max per event;5x5 tower max per event (GeV)", 5000,
0514                0, 50));
0515 
0516   return Fun4AllReturnCodes::EVENT_OK;
0517 }
0518 
0519 int QAG4SimulationCalorimeter::process_event_Tower(PHCompositeNode *topNode)
0520 {
0521   const std::string detector(_calo_name);
0522 
0523   if (Verbosity() > 2)
0524   {
0525     std::cout << "QAG4SimulationCalorimeter::process_event_Tower() entered" << std::endl;
0526   }
0527 
0528   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0529   assert(hm);
0530   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0531       get_histo_prefix() + "_Normalization"));
0532   assert(h_norm);
0533 
0534   std::string const towernodename = "TOWER_CALIB_" + detector;
0535   // Grab the towers
0536   RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,
0537                                                                     towernodename);
0538   std::string const towerinfonodename = "TOWERINFO_CALIB_" + detector;
0539   TowerInfoContainer *towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerinfonodename);
0540   if (!towers && !flag(kProcessTowerinfo))
0541   {
0542     std::cout << PHWHERE << ": Could not find node " << towernodename
0543               << std::endl;
0544     return Fun4AllReturnCodes::ABORTRUN;
0545   }
0546   if (!towerinfos && flag(kProcessTowerinfo))
0547   {
0548     std::cout << PHWHERE << ": Could not find node " << towerinfonodename
0549               << std::endl;
0550     return Fun4AllReturnCodes::ABORTRUN;
0551   }
0552   std::string const towergeomnodename = "TOWERGEOM_" + detector;
0553   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
0554       topNode, towergeomnodename);
0555   if (!towergeom)
0556   {
0557     std::cout << PHWHERE << ": Could not find node " << towergeomnodename
0558               << std::endl;
0559     return Fun4AllReturnCodes::ABORTRUN;
0560   }
0561 
0562   static const int max_size = 5;
0563   std::map<int, std::string> size_label;
0564   size_label[1] = "1x1";
0565   size_label[2] = "2x2";
0566   size_label[3] = "3x3";
0567   size_label[4] = "4x4";
0568   size_label[5] = "5x5";
0569   std::map<int, double> max_energy;
0570   std::map<int, TH1F *> energy_hist_list;
0571   std::map<int, TH1F *> max_energy_hist_list;
0572 
0573   for (int size = 1; size <= max_size; ++size)
0574   {
0575     max_energy[size] = 0;
0576 
0577     TH1F *h = dynamic_cast<TH1F *>(hm->getHisto(
0578         get_histo_prefix() + "_Tower_" + size_label[size]));
0579     assert(h);
0580     energy_hist_list[size] = h;
0581     h = dynamic_cast<TH1F *>(hm->getHisto(
0582         get_histo_prefix() + "_Tower_" + size_label[size] + "_max"));
0583     assert(h);
0584     max_energy_hist_list[size] = h;
0585   }
0586 
0587   h_norm->Fill("Tower", towergeom->size());  // total tower count
0588   h_norm->Fill("Tower Hit", towers->size());
0589 
0590   for (int binphi = 0; binphi < towergeom->get_phibins(); ++binphi)
0591   {
0592     for (int bineta = 0; bineta < towergeom->get_etabins(); ++bineta)
0593     {
0594       for (int size = 1; size <= max_size; ++size)
0595       {
0596         // for 2x2 and 4x4 use slide-2 window as implemented in DAQ
0597         if ((size == 2 || size == 4) && ((binphi % 2 != 0) && (bineta % 2 != 0)))
0598         {
0599           continue;
0600         }
0601 
0602         double energy = 0;
0603 
0604         // sliding window made from 2x2 sums
0605         for (int iphi = binphi; iphi < binphi + size; ++iphi)
0606         {
0607           for (int ieta = bineta; ieta < bineta + size; ++ieta)
0608           {
0609             if (ieta > towergeom->get_etabins())
0610             {
0611               continue;
0612             }
0613 
0614             // wrap around
0615             int wrapphi = iphi;
0616             assert(wrapphi >= 0);
0617             if (wrapphi >= towergeom->get_phibins())
0618             {
0619               wrapphi = wrapphi - towergeom->get_phibins();
0620             }
0621             if (!flag(kProcessTowerinfo))
0622             {
0623               RawTower *tower = towers->getTower(ieta, wrapphi);
0624 
0625               if (tower)
0626               {
0627                 const double e_intput = tower->get_energy();
0628 
0629                 energy += e_intput;
0630               }
0631             }
0632             else
0633             {
0634               unsigned int const towerkey = _calo_name == "CEMC" ? TowerInfoDefs::encode_emcal(ieta, wrapphi) : TowerInfoDefs::encode_hcal(ieta, wrapphi);
0635               TowerInfo *towerinfo = towerinfos->get_tower_at_key(towerkey);
0636               if (towerinfo)
0637               {
0638                 const double e_intput = towerinfo->get_energy();
0639                 energy += e_intput;
0640               }
0641             }
0642           }
0643         }
0644 
0645         energy_hist_list[size]->Fill(energy == 0 ? 9.1e-4 : energy);  // trick to fill 0 energy tower to the first bin
0646 
0647         max_energy[size] = std::max(energy, max_energy[size]);
0648 
0649       }  //          for (int size = 1; size <= 4; ++size)
0650     }
0651   }
0652 
0653   for (int size = 1; size <= max_size; ++size)
0654   {
0655     max_energy_hist_list[size]->Fill(max_energy[size]);
0656   }
0657   return Fun4AllReturnCodes::EVENT_OK;
0658 }
0659 
0660 int QAG4SimulationCalorimeter::Init_Cluster(PHCompositeNode * /*topNode*/)
0661 {
0662   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0663   assert(hm);
0664 
0665   hm->registerHisto(
0666       new TH1F(TString(get_histo_prefix()) + "_Cluster_BestMatchERatio",  //
0667                TString(_calo_name) + " best matched cluster E/E_{Truth};E_{Cluster}/E_{Truth}", 150,
0668                0, 1.5));
0669 
0670   hm->registerHisto(
0671       new TH2F(TString(get_histo_prefix()) + "_Cluster_LateralTruthProjection",  //
0672                TString(_calo_name) + " best cluster lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
0673                200, -15, 15, 200, -15, 15));
0674 
0675   return Fun4AllReturnCodes::EVENT_OK;
0676 }
0677 
0678 int QAG4SimulationCalorimeter::process_event_Cluster(PHCompositeNode *topNode)
0679 {
0680   if (Verbosity() > 2)
0681   {
0682     std::cout << "QAG4SimulationCalorimeter::process_event_Cluster() entered"
0683               << std::endl;
0684   }
0685 
0686   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0687   assert(hm);
0688   std::string const towergeomnodename = "TOWERGEOM_" + _calo_name;
0689   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
0690       topNode, towergeomnodename);
0691   if (!towergeom)
0692   {
0693     std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
0694     return Fun4AllReturnCodes::ABORTRUN;
0695   }
0696 
0697   // get a cluster count
0698   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(get_histo_prefix() + "_Normalization"));
0699   assert(h_norm);
0700 
0701   std::string const nodename = "CLUSTER_" + _calo_name;
0702   RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, nodename);
0703   assert(clusters);
0704   h_norm->Fill("Cluster", clusters->size());
0705 
0706   // get primary
0707   assert(_truth_container);
0708   assert(!_truth_container->GetMap().empty());
0709   PHG4Particle *last_primary = _truth_container->GetMap().rbegin()->second;
0710   assert(last_primary);
0711 
0712   if (Verbosity() > 2)
0713   {
0714     std::cout
0715         << "QAG4SimulationCalorimeter::process_event_Cluster() handle this truth particle"
0716         << std::endl;
0717     last_primary->identify();
0718   }
0719 
0720   assert(_caloevalstack);
0721   CaloRawClusterEval *clustereval = _caloevalstack->get_rawcluster_eval();
0722   clustereval->set_usetowerinfo(flag(kProcessTowerinfo));
0723   assert(clustereval);
0724 
0725   TH1F *h = dynamic_cast<TH1F *>(hm->getHisto(
0726       get_histo_prefix() + "_Cluster_BestMatchERatio"));
0727   assert(h);
0728 
0729   RawCluster *cluster = clustereval->best_cluster_from(last_primary);
0730   if (cluster)
0731   {
0732     // has a cluster matched and best cluster selected
0733 
0734     if (Verbosity() > 3)
0735     {
0736       std::cout << "QAG4SimulationCalorimeter::process_event_Cluster::"
0737                 << _calo_name << " - get cluster with energy "
0738                 << cluster->get_energy() << " VS primary energy "
0739                 << last_primary->get_e() << std::endl;
0740     }
0741 
0742     h->Fill(cluster->get_energy() / (last_primary->get_e() + 1e-9));  // avoids divide zero
0743 
0744     // now work on the projection:
0745     const CLHEP::Hep3Vector hit(cluster->get_position());
0746 
0747     const PHG4VtxPoint *primary_vtx =  //
0748         _truth_container->GetPrimaryVtx(last_primary->get_vtx_id());
0749     assert(primary_vtx);
0750     if (Verbosity() > 2)
0751     {
0752       std::cout
0753           << "QAG4SimulationCalorimeter::process_event_Cluster() handle this vertex"
0754           << std::endl;
0755       primary_vtx->identify();
0756     }
0757 
0758     const CLHEP::Hep3Vector vertex(primary_vtx->get_x(), primary_vtx->get_y(),
0759                                    primary_vtx->get_z());
0760 
0761     // projection axis
0762     CLHEP::Hep3Vector axis_proj(last_primary->get_px(), last_primary->get_py(),
0763                                 last_primary->get_pz());
0764     if (axis_proj.mag() == 0)
0765     {
0766       axis_proj.set(0, 0, 1);
0767     }
0768     axis_proj = axis_proj.unit();
0769 
0770     // azimuthal direction axis
0771     CLHEP::Hep3Vector axis_azimuth = axis_proj.cross(CLHEP::Hep3Vector(0, 0, 1));
0772     if (axis_azimuth.mag() == 0)
0773     {
0774       axis_azimuth.set(1, 0, 0);
0775     }
0776     axis_azimuth = axis_azimuth.unit();
0777 
0778     // polar direction axis
0779     CLHEP::Hep3Vector axis_polar = axis_proj.cross(axis_azimuth);
0780     assert(axis_polar.mag() > 0);
0781     axis_polar = axis_polar.unit();
0782 
0783     TH2F *hlat = dynamic_cast<TH2F *>(hm->getHisto(
0784         get_histo_prefix() + "_Cluster_LateralTruthProjection"));
0785     assert(hlat);
0786 
0787     const double hit_azimuth = axis_azimuth.dot(hit - vertex);
0788     const double hit_polar = axis_polar.dot(hit - vertex);
0789     hlat->Fill(hit_polar, hit_azimuth);
0790   }
0791   else
0792   {
0793     if (Verbosity() > 3)
0794     {
0795       std::cout << "QAG4SimulationCalorimeter::process_event_Cluster::"
0796                 << _calo_name << " - missing cluster !";
0797     }
0798     h->Fill(0);  // no cluster matched
0799   }
0800 
0801   return Fun4AllReturnCodes::EVENT_OK;
0802 }
[ 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!