Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "QAG4SimulationCalorimeterSum.h"
0002 
0003 #include <qautils/QAHistManagerDef.h>
0004 
0005 #include <g4eval/CaloEvalStack.h>
0006 #include <g4eval/CaloRawClusterEval.h>
0007 #include <g4eval/SvtxEvalStack.h>
0008 
0009 #include <g4main/PHG4Particle.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011 
0012 #include <calobase/RawCluster.h>
0013 #include <calobase/RawTower.h>
0014 #include <calobase/RawTowerContainer.h>
0015 #include <calobase/RawTowerGeomContainer.h>
0016 
0017 #include <trackbase_historic/SvtxTrack.h>
0018 
0019 #include <g4eval/SvtxTrackEval.h>  // for SvtxTrackEval
0020 
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/SubsysReco.h>
0024 
0025 #include <phool/getClass.h>
0026 
0027 #include <TAxis.h>
0028 #include <TH1.h>
0029 #include <TH2.h>
0030 #include <TString.h>
0031 
0032 #include <cassert>
0033 #include <cmath>
0034 #include <iostream>
0035 #include <iterator>  // for reverse_iterator
0036 #include <limits>
0037 #include <utility>  // for pair
0038 #include <vector>
0039 
0040 QAG4SimulationCalorimeterSum::QAG4SimulationCalorimeterSum(
0041     QAG4SimulationCalorimeterSum::enu_flags flags)
0042   : SubsysReco("QAG4SimulationCalorimeterSum")
0043   , _flags(flags)
0044   , _calo_name_cemc("CEMC")
0045   , _calo_name_hcalin("HCALIN")
0046   , _calo_name_hcalout("HCALOUT")
0047   , _truth_container(nullptr)
0048   , _magField(+1.4)
0049 {
0050 }
0051 
0052 int QAG4SimulationCalorimeterSum::InitRun(PHCompositeNode *topNode)
0053 {
0054   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0055                                                                 "G4TruthInfo");
0056   if (!_truth_container)
0057   {
0058     std::cout << "QAG4SimulationCalorimeterSum::InitRun - Fatal Error - "
0059               << "unable to find DST node "
0060               << "G4TruthInfo" << std::endl;
0061     assert(_truth_container);
0062   }
0063 
0064   if (flag(kProcessCluster))
0065   {
0066     if (!_caloevalstack_cemc)
0067     {
0068       _caloevalstack_cemc.reset(
0069           new CaloEvalStack(topNode, _calo_name_cemc));
0070       _caloevalstack_cemc->set_strict(true);
0071       _caloevalstack_cemc->set_verbosity(Verbosity() + 1);
0072     }
0073     if (!_caloevalstack_hcalin)
0074     {
0075       _caloevalstack_hcalin.reset(
0076           new CaloEvalStack(topNode, _calo_name_hcalin));
0077       _caloevalstack_hcalin->set_strict(true);
0078       _caloevalstack_hcalin->set_verbosity(Verbosity() + 1);
0079     }
0080     if (!_caloevalstack_hcalout)
0081     {
0082       _caloevalstack_hcalout.reset(
0083           new CaloEvalStack(topNode, _calo_name_hcalout));
0084       _caloevalstack_hcalout->set_strict(true);
0085       _caloevalstack_hcalout->set_verbosity(Verbosity() + 1);
0086     }
0087   }
0088 
0089   if (flag(kProcessTrackProj))
0090   {
0091     if (!_svtxevalstack)
0092     {
0093       _svtxevalstack.reset(new SvtxEvalStack(topNode));
0094       _svtxevalstack->set_strict(true);
0095       _svtxevalstack->set_verbosity(Verbosity() + 1);
0096     }
0097   }
0098   return Fun4AllReturnCodes::EVENT_OK;
0099 }
0100 
0101 int QAG4SimulationCalorimeterSum::Init(PHCompositeNode *topNode)
0102 {
0103   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0104   assert(hm);
0105   TH1D *h = new TH1D(TString(get_histo_prefix()) + "Normalization",  //
0106                      TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
0107   int i = 1;
0108   h->GetXaxis()->SetBinLabel(i++, "Event");
0109   h->GetXaxis()->SetBinLabel(i++, (_calo_name_cemc + " Tower").c_str());
0110   h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalin + " Tower").c_str());
0111   h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalout + " Tower").c_str());
0112   h->GetXaxis()->SetBinLabel(i++, (_calo_name_cemc + " Cluster").c_str());
0113   h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalin + " Cluster").c_str());
0114   h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalout + " Cluster").c_str());
0115   h->GetXaxis()->SetBinLabel(i++, "Track");
0116   h->GetXaxis()->LabelsOption("v");
0117   hm->registerHisto(h);
0118 
0119   //  if (flag(kProcessTower))
0120   //    {
0121   //      if (Verbosity() >= 1)
0122   //        std::cout << "QAG4SimulationCalorimeterSum::Init - Process tower occupancies"
0123   //            << std::endl;
0124   //      Init_Tower(topNode);
0125   //    }
0126   if (flag(kProcessCluster))
0127   {
0128     if (Verbosity() >= 1)
0129     {
0130       std::cout << "QAG4SimulationCalorimeterSum::Init - Process tower occupancies"
0131                 << std::endl;
0132     }
0133     Init_Cluster(topNode);
0134   }
0135 
0136   if (flag(kProcessTrackProj))
0137   {
0138     if (Verbosity() >= 1)
0139     {
0140       std::cout << "QAG4SimulationCalorimeterSum::Init - Process sampling fraction"
0141                 << std::endl;
0142     }
0143     Init_TrackProj(topNode);
0144   }
0145   return Fun4AllReturnCodes::EVENT_OK;
0146 }
0147 
0148 int QAG4SimulationCalorimeterSum::process_event(PHCompositeNode *topNode)
0149 {
0150   if (Verbosity() > 2)
0151   {
0152     std::cout << "QAG4SimulationCalorimeterSum::process_event() entered" << std::endl;
0153   }
0154 
0155   if (_caloevalstack_cemc)
0156   {
0157     _caloevalstack_cemc->next_event(topNode);
0158   }
0159   if (_caloevalstack_hcalin)
0160   {
0161     _caloevalstack_hcalin->next_event(topNode);
0162   }
0163   if (_caloevalstack_hcalout)
0164   {
0165     _caloevalstack_hcalout->next_event(topNode);
0166   }
0167   if (_svtxevalstack)
0168   {
0169     _svtxevalstack->next_event(topNode);
0170   }
0171 
0172   //  if (flag(kProcessTower))
0173   //    {
0174   //      int ret = process_event_Tower(topNode);
0175   //
0176   //      if (ret != Fun4AllReturnCodes::EVENT_OK)
0177   //        return ret;
0178   //    }
0179 
0180   if (flag(kProcessCluster))
0181   {
0182     int const ret = process_event_Cluster(topNode);
0183 
0184     if (ret != Fun4AllReturnCodes::EVENT_OK)
0185     {
0186       return ret;
0187     }
0188   }
0189 
0190   if (flag(kProcessTrackProj))
0191   {
0192     int const ret = process_event_TrackProj(topNode);
0193 
0194     if (ret != Fun4AllReturnCodes::EVENT_OK)
0195     {
0196       return ret;
0197     }
0198   }
0199 
0200   // at the end, count success events
0201   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0202   assert(hm);
0203   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0204       get_histo_prefix() + "Normalization"));
0205   assert(h_norm);
0206   h_norm->Fill("Event", 1);
0207 
0208   return Fun4AllReturnCodes::EVENT_OK;
0209 }
0210 
0211 std::string
0212 QAG4SimulationCalorimeterSum::get_histo_prefix()
0213 {
0214   return "h_QAG4Sim_CalorimeterSum_";
0215 }
0216 
0217 PHG4Particle *
0218 QAG4SimulationCalorimeterSum::get_truth_particle()
0219 {
0220   // get last primary
0221   assert(_truth_container);
0222   assert(!_truth_container->GetMap().empty());
0223   PHG4Particle *last_primary = _truth_container->GetMap().rbegin()->second;
0224   assert(last_primary);
0225 
0226   return last_primary;
0227 }
0228 
0229 int QAG4SimulationCalorimeterSum::Init_TrackProj(PHCompositeNode * /*topNode*/)
0230 {
0231   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0232   assert(hm);
0233 
0234   hm->registerHisto(
0235       new TH2F(
0236           TString(get_histo_prefix()) + TString(_calo_name_cemc.c_str()) + "_TrackProj",  //
0237           TString(_calo_name_cemc.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
0238           // NOLINTNEXTLINE(bugprone-integer-division)
0239           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
0240           // NOLINTNEXTLINE(bugprone-integer-division)
0241           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
0242 
0243   hm->registerHisto(
0244       new TH2F(
0245           TString(get_histo_prefix()) + TString(_calo_name_hcalin.c_str()) + "_TrackProj",  //
0246           TString(_calo_name_hcalin.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
0247           // NOLINTNEXTLINE(bugprone-integer-division)
0248           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
0249           // NOLINTNEXTLINE(bugprone-integer-division)
0250           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
0251 
0252   hm->registerHisto(
0253       new TH2F(
0254           TString(get_histo_prefix()) + TString(_calo_name_hcalout.c_str()) + "_TrackProj",  //
0255           TString(_calo_name_hcalout.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
0256           // NOLINTNEXTLINE(bugprone-integer-division)
0257           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
0258           // NOLINTNEXTLINE(bugprone-integer-division)
0259           (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
0260 
0261   hm->registerHisto(
0262       new TH1F(TString(get_histo_prefix()) + "TrackProj_3x3Tower_EP",  //
0263                "Tower 3x3 sum /E_{Truth};#Sigma_{3x3}[E_{Tower}] / total truth energy",
0264                150, 0, 1.5));
0265 
0266   hm->registerHisto(
0267       new TH1F(TString(get_histo_prefix()) + "TrackProj_5x5Tower_EP",  //
0268                "Tower 5x5 sum /E_{Truth};#Sigma_{5x5}[E_{Tower}] / total truth energy",
0269                150, 0, 1.5));
0270 
0271   return Fun4AllReturnCodes::EVENT_OK;
0272 }
0273 
0274 int QAG4SimulationCalorimeterSum::process_event_TrackProj(PHCompositeNode *topNode)
0275 {
0276   if (Verbosity() > 2)
0277   {
0278     std::cout << "QAG4SimulationCalorimeterSum::process_event_TrackProj() entered"
0279               << std::endl;
0280   }
0281 
0282   PHG4Particle *primary = get_truth_particle();
0283   if (!primary)
0284   {
0285     return Fun4AllReturnCodes::DISCARDEVENT;
0286   }
0287 
0288   SvtxTrackEval *trackeval = _svtxevalstack->get_track_eval();
0289   assert(trackeval);
0290   SvtxTrack *track = trackeval->best_track_from(primary);
0291   if (!track)
0292   {
0293     return Fun4AllReturnCodes::EVENT_OK;  // not through the whole event for missing track.
0294   }
0295 
0296   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0297   assert(hm);
0298   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0299       get_histo_prefix() + "Normalization"));
0300   assert(h_norm);
0301   h_norm->Fill("Track", 1);
0302 
0303   {
0304     TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
0305         (get_histo_prefix()) + "TrackProj_3x3Tower_EP"));
0306     assert(hsum);
0307 
0308     hsum->Fill(
0309         (track->get_cal_energy_3x3(SvtxTrack::CEMC) + track->get_cal_energy_3x3(SvtxTrack::HCALIN) + track->get_cal_energy_3x3(SvtxTrack::HCALOUT)) / (primary->get_e() + 1e-9));
0310   }
0311   {
0312     TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
0313         (get_histo_prefix()) + "TrackProj_5x5Tower_EP"));
0314     assert(hsum);
0315 
0316     hsum->Fill(
0317         (track->get_cal_energy_5x5(SvtxTrack::CEMC) + track->get_cal_energy_5x5(SvtxTrack::HCALIN) + track->get_cal_energy_5x5(SvtxTrack::HCALOUT)) / (primary->get_e() + 1e-9));
0318   }
0319 
0320   eval_trk_proj(_calo_name_cemc, track, topNode);
0321   eval_trk_proj(_calo_name_hcalin, track, topNode);
0322   eval_trk_proj(_calo_name_hcalout, track, topNode);
0323 
0324   return Fun4AllReturnCodes::EVENT_OK;
0325 }
0326 
0327 bool QAG4SimulationCalorimeterSum::eval_trk_proj(const std::string &detector, SvtxTrack *track,
0328                                                  PHCompositeNode *topNode)
0329 // Track projections
0330 {
0331   assert(track);
0332   assert(topNode);
0333 
0334   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0335   assert(hm);
0336 
0337   TH2F *h2_proj = dynamic_cast<TH2F *>(hm->getHisto(
0338       (get_histo_prefix()) + detector + "_TrackProj"));
0339   assert(h2_proj);
0340 
0341   // pull the tower geometry
0342   std::string const towergeonodename = "TOWERGEOM_" + detector;
0343   RawTowerGeomContainer *towergeo = findNode::getClass<RawTowerGeomContainer>(
0344       topNode, towergeonodename);
0345   assert(towergeo);
0346 
0347   if (Verbosity() > 2)
0348   {
0349     towergeo->identify();
0350   }
0351   // pull the towers
0352   std::string const towernodename = "TOWER_CALIB_" + detector;
0353   RawTowerContainer *towerList = findNode::getClass<RawTowerContainer>(topNode,
0354                                                                        towernodename);
0355   assert(towerList);
0356 
0357   if (Verbosity() > 3)
0358   {
0359     std::cout << __PRETTY_FUNCTION__ << " - info - handling track track p = ("
0360               << track->get_px() << ", " << track->get_py() << ", "
0361               << track->get_pz() << "), v = (" << track->get_x() << ", "
0362               << track->get_z() << ", " << track->get_z() << ")"
0363               << " size_states "
0364               << track->size_states() << std::endl;
0365   }
0366 
0367   // curved tracks inside mag field
0368   // straight projections thereafter
0369 
0370   std::vector<double> point;
0371   point.assign(3, std::numeric_limits<double>::quiet_NaN());
0372 
0373   //  const double radius = towergeo->get_radius() + towergeo->get_thickness() * 0.5;
0374 
0375   //  PHG4HoughTransform::projectToRadius(track, _magField, radius, point);
0376 
0377   if (std::isnan(point[0]) || std::isnan(point[1]) || std::isnan(point[2]))
0378   {
0379     // std::cout << __PRETTY_FUNCTION__ << "::" << Name()
0380     //      << " - Error - track extrapolation failure:";
0381     // track->identify();
0382     return false;
0383   }
0384 
0385   assert(!std::isnan(point[0]));
0386   assert(!std::isnan(point[1]));
0387   assert(!std::isnan(point[2]));
0388 
0389   double const x = point[0];
0390   double const y = point[1];
0391   double const z = point[2];
0392 
0393   double const phi = atan2(y, x);
0394   double const eta = asinh(z / sqrt(x * x + y * y));
0395 
0396   // calculate 3x3 tower energy
0397   int const binphi = towergeo->get_phibin(phi);
0398   int const bineta = towergeo->get_etabin(eta);
0399 
0400   double etabin_width = towergeo->get_etabounds(bineta).second - towergeo->get_etabounds(bineta).first;
0401   if (bineta > 1 && bineta < towergeo->get_etabins() - 1)
0402   {
0403     etabin_width = (towergeo->get_etacenter(bineta + 1) - towergeo->get_etacenter(bineta - 1)) / 2.;
0404   }
0405 
0406   double const phibin_width = towergeo->get_phibounds(binphi).second - towergeo->get_phibounds(binphi).first;
0407 
0408   assert(etabin_width > 0);
0409   assert(phibin_width > 0);
0410 
0411   const double etabin_shift = (towergeo->get_etacenter(bineta) - eta) / etabin_width;
0412   const double phibin_shift = (fmod(
0413                                    towergeo->get_phicenter(binphi) - phi + 5 * M_PI, 2 * M_PI) -
0414                                M_PI) /
0415                               phibin_width;
0416 
0417   if (Verbosity() > 3)
0418   {
0419     std::cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
0420               << ", " << y << ", " << z << ")"
0421               << ", eta " << eta << ", phi" << phi
0422               << ", bineta " << bineta << " - ["
0423               << towergeo->get_etabounds(bineta).first << ", "
0424               << towergeo->get_etabounds(bineta).second << "], etabin_shift = "
0425               << etabin_shift << ", binphi " << binphi << " - ["
0426               << towergeo->get_phibounds(binphi).first << ", "
0427               << towergeo->get_phibounds(binphi).second << "], phibin_shift = "
0428               << phibin_shift << std::endl;
0429   }
0430 
0431   const int bin_search_range = (Max_N_Tower - 1) / 2;
0432   for (int iphi = binphi - bin_search_range; iphi <= binphi + bin_search_range;
0433        ++iphi)
0434   {
0435     for (int ieta = bineta - bin_search_range;
0436          ieta <= bineta + bin_search_range; ++ieta)
0437     {
0438       // wrap around
0439       int wrapphi = iphi;
0440       if (wrapphi < 0)
0441       {
0442         wrapphi = towergeo->get_phibins() + wrapphi;
0443       }
0444       if (wrapphi >= towergeo->get_phibins())
0445       {
0446         wrapphi = wrapphi - towergeo->get_phibins();
0447       }
0448 
0449       // edges
0450       if (ieta < 0)
0451       {
0452         continue;
0453       }
0454       if (ieta >= towergeo->get_etabins())
0455       {
0456         continue;
0457       }
0458 
0459       if (Verbosity() > 3)
0460       {
0461         std::cout << __PRETTY_FUNCTION__ << " - info - processing tower"
0462                   << ", bineta " << ieta << " - ["
0463                   << towergeo->get_etabounds(ieta).first << ", "
0464                   << towergeo->get_etabounds(ieta).second << "]"
0465                   << ", binphi "
0466                   << wrapphi << " - [" << towergeo->get_phibounds(wrapphi).first
0467                   << ", " << towergeo->get_phibounds(wrapphi).second << "]" << std::endl;
0468       }
0469 
0470       RawTower *tower = towerList->getTower(ieta, wrapphi);
0471       double energy = 0;
0472       if (tower)
0473       {
0474         if (Verbosity() > 1)
0475         {
0476           std::cout << __PRETTY_FUNCTION__ << " - info - tower " << ieta << " "
0477                     << wrapphi << " energy = " << tower->get_energy() << std::endl;
0478         }
0479 
0480         energy = tower->get_energy();
0481       }
0482 
0483       h2_proj->Fill(ieta - bineta + etabin_shift,
0484                     iphi - binphi + phibin_shift, energy);
0485 
0486     }  //            for (int ieta = bineta-1; ieta < bineta+2; ++ieta) {
0487   }
0488 
0489   return true;
0490 }  //       // Track projections
0491 
0492 // int
0493 // QAG4SimulationCalorimeterSum::Init_Tower(PHCompositeNode *topNode)
0494 //{
0495 //
0496 //   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0497 //   assert(hm);
0498 //
0499 ////  TH1F * h = nullptr;
0500 //
0501 //  return Fun4AllReturnCodes::EVENT_OK;
0502 //}
0503 //
0504 // int
0505 // QAG4SimulationCalorimeterSum::process_event_Tower(PHCompositeNode *topNode)
0506 //{
0507 //  if (Verbosity() > 2)
0508 //    std::cout << "QAG4SimulationCalorimeterSum::process_event_Tower() entered"
0509 //        << std::endl;
0510 //
0511 //  return Fun4AllReturnCodes::EVENT_OK;
0512 //}
0513 
0514 int QAG4SimulationCalorimeterSum::Init_Cluster(PHCompositeNode * /*topNode*/)
0515 {
0516   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0517   assert(hm);
0518 
0519   hm->registerHisto(
0520       new TH2F(
0521           TString(get_histo_prefix()) + "Cluster_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str(),  //
0522           TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_cemc.c_str()) + ": best cluster energy;" + TString(_calo_name_cemc.c_str()) + " cluster energy (GeV);" + TString(_calo_name_hcalin.c_str()) + " cluster energy (GeV)",
0523           70, 0, 70, 70, 0, 70));
0524 
0525   hm->registerHisto(
0526       new TH2F(
0527           TString(get_histo_prefix()) + "Cluster_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str() + "_" + _calo_name_hcalout.c_str(),  //
0528           TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_hcalout.c_str()) + ": best cluster energy;" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " cluster energy (GeV);" + TString(_calo_name_hcalout.c_str()) + " cluster energy (GeV)",
0529           70, 0, 70, 70, 0, 70));
0530 
0531   hm->registerHisto(
0532       new TH1F(TString(get_histo_prefix()) + "Cluster_EP",  //
0533                "Total Cluster E_{Reco}/E_{Truth};Reco cluster energy sum / total truth energy",
0534                150, 0, 1.5));
0535 
0536   hm->registerHisto(
0537       new TH1F(
0538           TString(get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str(),  //
0539           "Energy ratio " + TString(_calo_name_cemc.c_str()) + " VS " + TString(_calo_name_hcalin.c_str()) + ";Best cluster " + TString(_calo_name_cemc.c_str()) + " / (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + ")", 110, 0, 1.1));
0540 
0541   hm->registerHisto(
0542       new TH1F(
0543           TString(get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str() + "_" + TString(_calo_name_hcalout.c_str()),  //
0544           "Energy ratio " + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_hcalout.c_str()) + ";Best cluster (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + ") / (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " + " + TString(_calo_name_hcalout.c_str()) + ")", 110, 0, 1.1));
0545 
0546   return Fun4AllReturnCodes::EVENT_OK;
0547 }
0548 
0549 int QAG4SimulationCalorimeterSum::process_event_Cluster(PHCompositeNode * /*topNode*/)
0550 {
0551   if (Verbosity() > 2)
0552   {
0553     std::cout << "QAG4SimulationCalorimeterSum::process_event_Cluster() entered"
0554               << std::endl;
0555   }
0556 
0557   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0558   assert(hm);
0559 
0560   PHG4Particle *primary = get_truth_particle();
0561   if (!primary)
0562   {
0563     return Fun4AllReturnCodes::DISCARDEVENT;
0564   }
0565 
0566   RawCluster *cluster_cemc =
0567       _caloevalstack_cemc->get_rawcluster_eval()->best_cluster_from(primary);
0568   RawCluster *cluster_hcalin =
0569       _caloevalstack_hcalin->get_rawcluster_eval()->best_cluster_from(primary);
0570   RawCluster *cluster_hcalout =
0571       _caloevalstack_hcalout->get_rawcluster_eval()->best_cluster_from(primary);
0572 
0573   const double cluster_cemc_e = cluster_cemc ? cluster_cemc->get_energy() : 0;
0574   const double cluster_hcalin_e =
0575       cluster_hcalin ? cluster_hcalin->get_energy() : 0;
0576   const double cluster_hcalout_e =
0577       cluster_hcalout ? cluster_hcalout->get_energy() : 0;
0578 
0579   if (cluster_cemc_e + cluster_hcalin_e > 0)
0580   {
0581     TH2F *h2 = dynamic_cast<TH2F *>(hm->getHisto(
0582         (get_histo_prefix()) + "Cluster_" + _calo_name_cemc + "_" + _calo_name_hcalin));
0583     assert(h2);
0584 
0585     h2->Fill(cluster_cemc_e, cluster_hcalin_e);
0586 
0587     TH1F *hr = dynamic_cast<TH1F *>(hm->getHisto(
0588         (get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc + "_" + _calo_name_hcalin));
0589     assert(hr);
0590 
0591     hr->Fill(cluster_cemc_e / (cluster_cemc_e + cluster_hcalin_e));
0592   }
0593 
0594   if (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e > 0)
0595   {
0596     if (Verbosity())
0597     {
0598       std::cout << "QAG4SimulationCalorimeterSum::process_event_Cluster - "
0599                 << " cluster_cemc_e = " << cluster_cemc_e
0600                 << " cluster_hcalin_e = " << cluster_hcalin_e
0601                 << " cluster_hcalout_e = " << cluster_hcalout_e << " hr = "
0602                 << (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e)
0603                 << std::endl;
0604     }
0605 
0606     TH2F *h2 = dynamic_cast<TH2F *>(hm->getHisto(
0607         (get_histo_prefix()) + "Cluster_" + _calo_name_cemc + "_" + _calo_name_hcalin + "_" + _calo_name_hcalout));
0608     assert(h2);
0609 
0610     h2->Fill((cluster_cemc_e + cluster_hcalin_e), cluster_hcalout_e);
0611 
0612     TH1F *hr = dynamic_cast<TH1F *>(hm->getHisto(
0613         (get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc + "_" + _calo_name_hcalin + "_" + _calo_name_hcalout));
0614     assert(hr);
0615 
0616     hr->Fill(
0617         (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e));
0618 
0619     TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
0620         (get_histo_prefix()) + "Cluster_EP"));
0621     assert(hsum);
0622 
0623     hsum->Fill(
0624         (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e) / (primary->get_e() + 1e-9));
0625   }
0626 
0627   return Fun4AllReturnCodes::EVENT_OK;
0628 }
[ 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!