Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:59 

0001 
0002 
0003 #include "SynRadAna.h"
0004 
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 
0007 #include <fun4all/Fun4AllHistoManager.h>
0008 #include <fun4all/Fun4AllServer.h>
0009 
0010 #include <phhepmc/PHHepMCGenEvent.h>
0011 #include <phhepmc/PHHepMCGenEventMap.h>
0012 
0013 #include <g4eval/CaloEvalStack.h>
0014 #include <g4eval/CaloRawClusterEval.h>
0015 #include <g4eval/SvtxEvalStack.h>
0016 
0017 #include <g4main/PHG4Hit.h>
0018 #include <g4main/PHG4HitContainer.h>
0019 #include <g4main/PHG4Particle.h>
0020 #include <g4main/PHG4TruthInfoContainer.h>
0021 #include <g4main/PHG4VtxPoint.h>
0022 
0023 #include <calobase/RawCluster.h>
0024 #include <calobase/RawTower.h>
0025 #include <calobase/RawTowerContainer.h>
0026 #include <calobase/RawTowerGeomContainer.h>
0027 
0028 #include <mvtx/CylinderGeom_Mvtx.h>
0029 #include <mvtx/MvtxDefs.h>
0030 //#include <mvtx/MvtxHit.h>
0031 #include <trackbase/TrkrDefs.h>
0032 #include <trackbase/TrkrHit.h>  // for TrkrHit
0033 #include <trackbase/TrkrHitSet.h>
0034 #include <trackbase/TrkrHitSetContainer.h>
0035 #include <trackbase/TrkrHitTruthAssoc.h>
0036 #include <trackbase_historic/SvtxTrack.h>
0037 
0038 #include <g4eval/SvtxTrackEval.h>  // for SvtxTrackEval
0039 #include <g4eval/SvtxTruthEval.h>  // for SvtxTruthEval
0040 
0041 #include <fun4all/Fun4AllHistoManager.h>
0042 #include <fun4all/Fun4AllReturnCodes.h>
0043 #include <fun4all/SubsysReco.h>
0044 
0045 #include <phool/PHCompositeNode.h>
0046 #include <phool/PHDataNode.h>
0047 #include <phool/PHNode.h>  // for PHNode
0048 #include <phool/PHNodeIterator.h>
0049 #include <phool/PHObject.h>
0050 #include <phool/PHRandomSeed.h>
0051 #include <phool/getClass.h>
0052 #include <phool/phool.h>
0053 #include <phool/recoConsts.h>
0054 
0055 #include <HepMC/GenEvent.h>
0056 #include <HepMC/GenParticle.h>
0057 #include <HepMC/GenVertex.h>
0058 #include <HepMC/IteratorRange.h>
0059 #include <HepMC/SimpleVector.h>
0060 #include <HepMC/Units.h>
0061 
0062 #include <gsl/gsl_const.h>
0063 #include <gsl/gsl_randist.h>
0064 #include <gsl/gsl_rng.h>
0065 
0066 #include <TAxis.h>
0067 #include <TDatabasePDG.h>
0068 #include <TH1.h>
0069 #include <TH2.h>
0070 #include <TNamed.h>
0071 #include <TString.h>
0072 #include <TVector3.h>
0073 
0074 #include <phool/PHCompositeNode.h>
0075 #include <cassert>
0076 #include <cmath>
0077 #include <cstdlib>
0078 #include <iostream>
0079 #include <iterator>
0080 #include <list>
0081 #include <utility>
0082 
0083 using namespace std;
0084 
0085 //____________________________________________________________________________..
0086 SynRadAna::SynRadAna(const std::string &fname)
0087   : SubsysReco("SynRadAna")
0088   , _embedding_id(1)
0089   , m_eventWeight(0)
0090   , do_photon(true)
0091   , do_MVTXHits(true)
0092   , do_TPCHits(true)
0093   , m_outputFIle(fname)
0094 {
0095   if (Verbosity())
0096     cout << "SynRadAna::SynRadAna(const std::string &name) Calling ctor" << endl;
0097 }
0098 
0099 //____________________________________________________________________________..
0100 SynRadAna::~SynRadAna()
0101 {
0102   if (Verbosity())
0103     cout << "SynRadAna::~SynRadAna() Calling dtor" << endl;
0104 }
0105 
0106 //____________________________________________________________________________..
0107 int SynRadAna::Init(PHCompositeNode *topNode)
0108 {
0109   if (Verbosity())
0110     cout << "SynRadAna::Init(PHCompositeNode *topNode) Initializing" << endl;
0111 
0112   Fun4AllHistoManager *hm = getHistoManager();
0113   assert(hm);
0114 
0115   {
0116     TH1 *h(nullptr);
0117     // n events and n tracks histogram
0118     h = new TH1D(TString(get_histo_prefix()) + "Normalization",
0119                  TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
0120     int i = 1;
0121     h->GetXaxis()->SetBinLabel(i++, "Event");
0122     h->GetXaxis()->SetBinLabel(i++, "Photon");
0123     h->GetXaxis()->SetBinLabel(i++, "Flux");
0124     h->GetXaxis()->SetBinLabel(i++, "G4Particle");
0125     h->GetXaxis()->SetBinLabel(i++, "G4PrimaryParticle");
0126     h->GetXaxis()->SetBinLabel(i++, "G4Vertex");
0127     h->GetXaxis()->LabelsOption("v");
0128     hm->registerHisto(h);
0129   }
0130 
0131   for (const auto &hitnode : _node_postfix)
0132   {
0133     TH2D *h2(nullptr);
0134 
0135     h2 = new TH2D(TString(get_histo_prefix()) + hitnode + "_nHit",
0136                   "Hit sum;Sum number of hits", 2000, -.5, 2000 - .5, 2, .5, 2.5);
0137     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0138     h2->GetYaxis()->SetBinLabel(1, "Flux");
0139     h2->GetYaxis()->SetBinLabel(2, "Photon");
0140     hm->registerHisto(h2);
0141 
0142     h2 = new TH2D(TString(get_histo_prefix()) + hitnode + "_sumEdep",
0143                   "Hit sum energy distribution;Ionizing energy deposition [keV]", 2000, 0, 100, 2, .5, 2.5);
0144     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0145     h2->GetYaxis()->SetBinLabel(1, "Flux");
0146     h2->GetYaxis()->SetBinLabel(2, "Photon");
0147     hm->registerHisto(h2);
0148 
0149     h2 = new TH2D(TString(get_histo_prefix()) + hitnode + "_photonEnergy",
0150                   "Hit source photon;Photon energy [keV]", 2000, 0, 100, 2, .5, 2.5);
0151     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0152     h2->GetYaxis()->SetBinLabel(1, "Flux");
0153     h2->GetYaxis()->SetBinLabel(2, "Photon");
0154     hm->registerHisto(h2);
0155 
0156     h2 = new TH2D(TString(get_histo_prefix()) + hitnode + "_photonZ",
0157                   "#gamma z crossing interface facet [cm]", 800, -400, 400, 2, .5, 2.5);
0158     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0159     h2->GetYaxis()->SetBinLabel(1, "Flux");
0160     h2->GetYaxis()->SetBinLabel(2, "Photon");
0161     hm->registerHisto(h2);
0162   }
0163 
0164   for (vector<string>::const_iterator it = _tower_postfix.begin();
0165        it != _tower_postfix.end(); ++it)
0166   {
0167   }
0168 
0169   if (do_photon)
0170   {
0171     TH2D *h2 = new TH2D(TString(get_histo_prefix()) + "photonEnergy",
0172                         "Source photon;Photon energy [keV]", 2000, 0, 100, 2, .5, 2.5);
0173     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0174     h2->GetYaxis()->SetBinLabel(1, "Flux");
0175     h2->GetYaxis()->SetBinLabel(2, "Photon");
0176     hm->registerHisto(h2);
0177 
0178     h2 = new TH2D(TString(get_histo_prefix()) + "photonZ",
0179                   "#gamma z crossing interface facet [cm]", 800, -400, 400, 2, .5, 2.5);
0180     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0181     h2->GetYaxis()->SetBinLabel(1, "Flux");
0182     h2->GetYaxis()->SetBinLabel(2, "Photon");
0183     hm->registerHisto(h2);
0184   }
0185 
0186   if (do_MVTXHits)
0187   {
0188     TH2D *h2(nullptr);
0189 
0190     h2 = new TH2D(TString(get_histo_prefix()) + "MVTXHit" + "_nHit",
0191                   "Hit sum;Sum number of hits", 2000, -.5, 2000 - .5, 2, .5, 2.5);
0192     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0193     h2->GetYaxis()->SetBinLabel(1, "Flux");
0194     h2->GetYaxis()->SetBinLabel(2, "Photon");
0195     hm->registerHisto(h2);
0196 
0197     h2 = new TH2D(TString(get_histo_prefix()) + "MVTXHit" + "_nHit_Layer",
0198                   "Hit sum;Sum number of hits;layer", 2000, -.5, 2000 - .5, 3, -.5, 2.5);
0199     hm->registerHisto(h2);
0200   }
0201 
0202   if (do_TPCHits)
0203   {
0204     TH2D *h2(nullptr);
0205 
0206     h2 = new TH2D(TString(get_histo_prefix()) + "TPCHit" + "_nHit",
0207                   "Hit sum;Sum number of hits", 2000, -.5, 2000 - .5, 2, .5, 2.5);
0208     //  QAHistManagerDef::useLogBins(h->GetXaxis());
0209     h2->GetYaxis()->SetBinLabel(1, "Flux");
0210     h2->GetYaxis()->SetBinLabel(2, "Photon");
0211     hm->registerHisto(h2);
0212 
0213     h2 = new TH2D(TString(get_histo_prefix()) + "TPCHit" + "_nHit_Layer",
0214                   "Hit sum;Sum number of hits;layer", 2000, -.5, 2000 - .5, 61, -.5, 60.5);
0215     hm->registerHisto(h2);
0216   }
0217 
0218   return Fun4AllReturnCodes::EVENT_OK;
0219 }
0220 
0221 //____________________________________________________________________________..
0222 int SynRadAna::InitRun(PHCompositeNode *topNode)
0223 {
0224   if (Verbosity())
0225     cout << "SynRadAna::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << endl;
0226 
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229 
0230 //____________________________________________________________________________..
0231 int SynRadAna::process_event(PHCompositeNode *topNode)
0232 {
0233   if (Verbosity())
0234     cout << "SynRadAna::process_event(PHCompositeNode *topNode) Processing Event" << endl;
0235 
0236   static const double GeV2keV = 1e6;
0237 
0238   // histogram manager
0239   Fun4AllHistoManager *hm = getHistoManager();
0240   assert(hm);
0241   // n events and n tracks histogram
0242   TH1 *h_norm = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "Normalization"));
0243   assert(h_norm);
0244   h_norm->Fill("Event", 1);
0245 
0246   PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
0247       PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0248   assert(truthInfoList);
0249   h_norm->Fill("G4Particle", truthInfoList->size());
0250   h_norm->Fill("G4Vertex", truthInfoList->GetNumVertices());
0251   h_norm->Fill("G4PrimaryParticle", truthInfoList->GetNumPrimaryVertexParticles());
0252 
0253   // For pile-up simulation: define GenEventMap
0254   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0255   assert(genevtmap);
0256   PHHepMCGenEvent *genev = genevtmap->get(_embedding_id);
0257   assert(genev);
0258   HepMC::GenEvent *hepmc_evt = genev->getEvent();
0259   assert(hepmc_evt);
0260 
0261   if (Verbosity() >= 2)
0262   {
0263     hepmc_evt->print();
0264   }
0265 
0266   //weight calculation/normalization
0267   const auto &weightcontainer = hepmc_evt->weights();
0268   assert(weightcontainer.size() >= 2);
0269   h_norm->Fill("Flux", weightcontainer[0]);
0270   h_norm->Fill("Photon", weightcontainer[1]);
0271   if (weightcontainer[1] != 1)
0272   {
0273     cout << __PRETTY_FUNCTION__ << "- warning - weightcontainer[1] = " << weightcontainer[1]
0274          << " != 1. Can not precisely pin point each photon's flux!"
0275          << endl;
0276   }
0277   m_eventWeight = weightcontainer[0] / weightcontainer[1];
0278 
0279   // initial photons
0280   if (do_photon)
0281   {
0282     TH2 *h_photonEnergy = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "photonEnergy"));
0283     assert(h_photonEnergy);
0284     TH2 *h_photonZ = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "photonZ"));
0285     assert(h_photonZ);
0286 
0287     PHG4TruthInfoContainer::ConstRange primary_range =
0288         truthInfoList->GetPrimaryParticleRange();
0289 
0290     for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0291              primary_range.first;
0292          particle_iter != primary_range.second;
0293          ++particle_iter)
0294     {
0295       const auto particle = particle_iter->second;
0296 
0297       assert(particle);
0298       if (particle->get_pid() == 22)
0299       {
0300         const double photon_e_keV = particle->get_e() * GeV2keV;
0301         h_photonEnergy->Fill(photon_e_keV, "Flux", m_eventWeight);
0302         h_photonEnergy->Fill(photon_e_keV, "Photon", 1);
0303 
0304         PHG4VtxPoint *vtx = truthInfoList->GetVtx(particle->get_vtx_id());
0305         assert(vtx);
0306 
0307         const double photon_z = vtx->get_z();
0308         h_photonZ->Fill(photon_z, "Flux", m_eventWeight);
0309         h_photonZ->Fill(photon_z, "Photon", 1);
0310       }
0311       else
0312       {
0313         cout << __PRETTY_FUNCTION__ << "- warning - non-photon source!";
0314         particle->identify();
0315       }
0316     }  //          if (_load_all_particle) else
0317   }
0318 
0319   for (const auto &hitnode : _node_postfix)
0320   {
0321     PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode,
0322                                                                   string("G4HIT_" + hitnode));
0323     if (!hits)
0324     {
0325       cout << __PRETTY_FUNCTION__ << " - Fatal Error - missing " << string("G4HIT_" + hitnode) << endl;
0326     }
0327 
0328     assert(hits);
0329     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0330 
0331     if (Verbosity() >= 2)
0332       cout << "SynRadAna::process_event - processing "
0333            << hitnode << " and received " << hits->size() << " hits"
0334            << endl;
0335 
0336     int nhit(0);
0337     double sumEdep_keV(0);
0338 
0339     map<PHG4Particle *, pair<double, double> > primary_photon_map;
0340     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
0341          hit_iter != hit_range.second; hit_iter++)
0342     {
0343       PHG4Hit *hit = hit_iter->second;
0344       assert(hit);
0345 
0346       ++nhit;
0347 
0348       const double eion_keV = hit->get_eion() * GeV2keV;
0349       sumEdep_keV += eion_keV;
0350 
0351       PHG4Particle *hit_particle =
0352           truthInfoList->GetParticle(hit->get_trkid());
0353       assert(hit_particle);
0354 
0355       PHG4Particle *primary_particle = truthInfoList->GetParticle(hit_particle->get_primary_id());
0356       assert(primary_particle);
0357       if (primary_particle->get_pid() != 22)
0358       {
0359         cout << "SynRadAna::process_event - WARNING - unexpected primary particle that is not photon: ";
0360         primary_particle->identify();
0361 
0362         continue;
0363       }
0364       const double photon_e_keV = primary_particle->get_e() * GeV2keV;
0365 
0366       PHG4VtxPoint *vtx = truthInfoList->GetVtx(primary_particle->get_vtx_id());
0367       assert(vtx);
0368       const double photon_z = vtx->get_z();
0369 
0370       primary_photon_map.insert(
0371           make_pair(primary_particle,
0372                     make_pair(photon_e_keV, photon_z)));
0373     }
0374 
0375     TH2 *h_nHit = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + hitnode + "_nHit"));
0376     assert(h_nHit);
0377 
0378     TH2 *h_sumEdep = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + hitnode + "_sumEdep"));
0379     assert(h_sumEdep);
0380 
0381     TH2 *h_photonEnergy = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + hitnode + "_photonEnergy"));
0382     assert(h_photonEnergy);
0383 
0384     TH2 *h_photonZ = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + hitnode + "_photonZ"));
0385     assert(h_photonZ);
0386 
0387     h_nHit->Fill(nhit, "Flux", m_eventWeight);
0388     h_nHit->Fill(nhit, "Photon", 1);
0389 
0390     if (nhit > 0)
0391     {
0392       if (primary_photon_map.size() != 1)
0393       {
0394         cout << "SynRadAna::process_event - WARNING - primary_photon_map.size() = " << primary_photon_map.size() << endl;
0395 
0396         continue;
0397       }
0398 
0399       for (auto &pair : primary_photon_map)
0400       {
0401         const std::pair<double, double> & e_z = pair.second;
0402 
0403         h_photonEnergy->Fill(e_z.first, "Flux", m_eventWeight);
0404         h_photonEnergy->Fill(e_z.first, "Photon", 1);
0405 
0406         h_photonZ->Fill(e_z.second, "Flux", m_eventWeight);
0407         h_photonZ->Fill(e_z.second, "Photon", 1);
0408       }
0409 
0410       h_sumEdep->Fill(sumEdep_keV, "Flux", m_eventWeight);
0411       h_sumEdep->Fill(sumEdep_keV, "Photon", 1);
0412     }
0413 
0414   }  //  for (const auto &hitnode : _node_postfix)
0415 
0416   if (do_MVTXHits)
0417   {
0418     // Get the TrkrHitSetContainer node
0419     TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0420     if (!trkrhitsetcontainer)
0421     {
0422       cout << "Could not locate TRKR_HITSET node, quit! " << endl;
0423       exit(1);
0424     }
0425 
0426     int nhit(0);
0427     map<int, int> layer_nhit{{0, 0}, {0, 0}, {0, 0}};
0428     // We want all hitsets for the Mvtx
0429     TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::mvtxId);
0430     for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
0431          hitset_iter != hitset_range.second;
0432          ++hitset_iter)
0433     {
0434       // we have an itrator to one TrkrHitSet for the mvtx from the trkrHitSetContainer
0435       // get the hitset key so we can find the layer
0436       TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
0437       int layer = TrkrDefs::getLayer(hitsetkey);
0438       if (Verbosity() >= 2) cout << __PRETTY_FUNCTION__ << ": found MVTX hitset with key: " << hitsetkey << " in layer " << layer << endl;
0439 
0440       TrkrHitSet *hitset = hitset_iter->second;
0441       TrkrHitSet::ConstRange hit_range = hitset->getHits();
0442       for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
0443            hit_iter != hit_range.second;
0444            ++hit_iter)
0445       {
0446         TrkrHit *hit = hit_iter->second;
0447         assert(hit);
0448         if (Verbosity() >= 2)
0449         {
0450           cout << hit->getAdc() << "ADC hit: ";
0451           hit->identify();
0452         }
0453 
0454         if (hit->getAdc())
0455         {
0456           ++nhit;
0457 
0458           layer_nhit[layer] += 1;
0459         }
0460       }
0461     }
0462 
0463     TH2 *h_nHit = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "MVTXHit" + "_nHit"));
0464     assert(h_nHit);
0465 
0466     TH2 *h_nHit_Layer = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "MVTXHit" + "_nHit_Layer"));
0467     assert(h_nHit_Layer);
0468 
0469     h_nHit->Fill(nhit, "Flux", m_eventWeight);
0470     h_nHit->Fill(nhit, "Photon", 1);
0471 
0472     for (auto &pair : layer_nhit)
0473     {
0474       if (Verbosity() >= 2)
0475         cout << __PRETTY_FUNCTION__ << ": MVTX Summary: found " << pair.second << " hits in layer " << pair.first << endl;
0476 
0477       h_nHit_Layer->Fill(pair.second, pair.first, m_eventWeight);
0478     }
0479 
0480   }  //  if (do_MVTXHits)
0481 
0482   if (do_TPCHits)
0483   {
0484     // Get the TrkrHitSetContainer node
0485     TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0486     if (!trkrhitsetcontainer)
0487     {
0488       cout << "Could not locate TRKR_HITSET node, quit! " << endl;
0489       exit(1);
0490     }
0491 
0492     int nhit(0);
0493     map<int, int> layer_nhit{};
0494     // We want all hitsets for the TPC
0495     TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
0496     for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
0497          hitset_iter != hitset_range.second;
0498          ++hitset_iter)
0499     {
0500       // we have an itrator to one TrkrHitSet for the mvtx from the trkrHitSetContainer
0501       // get the hitset key so we can find the layer
0502       TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
0503       int layer = TrkrDefs::getLayer(hitsetkey);
0504       if (Verbosity() >= 2) cout << __PRETTY_FUNCTION__ << ": found TPC hitset with key: " << hitsetkey << " in layer " << layer << endl;
0505 
0506       TrkrHitSet *hitset = hitset_iter->second;
0507       TrkrHitSet::ConstRange hit_range = hitset->getHits();
0508       for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
0509            hit_iter != hit_range.second;
0510            ++hit_iter)
0511       {
0512         TrkrHit *hit = hit_iter->second;
0513         assert(hit);
0514         if (Verbosity() >= 2)
0515         {
0516           cout << hit->getAdc() << " ADC hit. " << endl;
0517         }
0518 
0519         if (hit->getAdc())
0520         {
0521           ++nhit;
0522 
0523           if (layer_nhit.find(layer) == layer_nhit.end())
0524             layer_nhit.insert(make_pair(layer, 0));
0525 
0526           layer_nhit[layer] += 1;
0527         }
0528       }
0529     }
0530 
0531     TH2 *h_nHit = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "TPCHit" + "_nHit"));
0532     assert(h_nHit);
0533 
0534     TH2 *h_nHit_Layer = dynamic_cast<TH2 *>(hm->getHisto(string(get_histo_prefix()) + "TPCHit" + "_nHit_Layer"));
0535     assert(h_nHit_Layer);
0536 
0537     h_nHit->Fill(nhit, "Flux", m_eventWeight);
0538     h_nHit->Fill(nhit, "Photon", 1);
0539 
0540     for (auto &pair : layer_nhit)
0541     {
0542       if (Verbosity() >= 2)
0543         cout << __PRETTY_FUNCTION__ << ": TPC summary:  found " << pair.second << " hits in layer " << pair.first << endl;
0544 
0545       h_nHit_Layer->Fill(pair.second, pair.first, m_eventWeight);
0546     }
0547 
0548   }  //  if (do_TPCHits)
0549 
0550   return Fun4AllReturnCodes::EVENT_OK;
0551 }
0552 
0553 //____________________________________________________________________________..
0554 int SynRadAna::ResetEvent(PHCompositeNode *topNode)
0555 {
0556   if (Verbosity())
0557     cout << "SynRadAna::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << endl;
0558   return Fun4AllReturnCodes::EVENT_OK;
0559 }
0560 
0561 //____________________________________________________________________________..
0562 int SynRadAna::EndRun(const int runnumber)
0563 {
0564   if (Verbosity())
0565     cout << "SynRadAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << endl;
0566   return Fun4AllReturnCodes::EVENT_OK;
0567 }
0568 
0569 //____________________________________________________________________________..
0570 int SynRadAna::End(PHCompositeNode *topNode)
0571 {
0572   if (Verbosity())
0573     cout << "SynRadAna::End(PHCompositeNode *topNode) This is the End..." << endl;
0574 
0575   getHistoManager()->dumpHistos(m_outputFIle);
0576 
0577   return Fun4AllReturnCodes::EVENT_OK;
0578 }
0579 
0580 //____________________________________________________________________________..
0581 int SynRadAna::Reset(PHCompositeNode *topNode)
0582 {
0583   if (Verbosity())
0584     cout << "SynRadAna::Reset(PHCompositeNode *topNode) being Reset" << endl;
0585   return Fun4AllReturnCodes::EVENT_OK;
0586 }
0587 
0588 //____________________________________________________________________________..
0589 void SynRadAna::Print(const std::string &what) const
0590 {
0591   cout << "SynRadAna::Print(const std::string &what) const Printing info for " << what << endl;
0592 }
0593 
0594 //! Get a pointer to the default hist manager for QA modules
0595 Fun4AllHistoManager *SynRadAna::
0596     getHistoManager()
0597 {
0598   static const string HistoManagerName("HistoManager_SynRadAna");
0599 
0600   Fun4AllServer *se = Fun4AllServer::instance();
0601   Fun4AllHistoManager *hm = se->getHistoManager(HistoManagerName);
0602 
0603   if (not hm)
0604   {
0605     //        cout
0606     //            << "QAHistManagerDef::get_HistoManager - Making Fun4AllHistoManager EMCalAna_HISTOS"
0607     //            << endl;
0608     hm = new Fun4AllHistoManager(HistoManagerName);
0609     se->registerHistoManager(hm);
0610   }
0611 
0612   assert(hm);
0613 
0614   return hm;
0615 }
0616 
0617 string
0618 SynRadAna::get_histo_prefix()
0619 {
0620   return string("h_") + Name() + string("_");
0621 }
[ 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!