g4simulation/g4main/PHG4ScoringManager.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PHG4ScoringManager.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PHG4ScoringManager.h"
0012
0013 #include "PHG4InEvent.h"
0014 #include "PHG4Particle.h"
0015
0016 #include <phhepmc/PHHepMCGenEvent.h>
0017 #include <phhepmc/PHHepMCGenEventMap.h>
0018
0019 #include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERBO...
0020 #include <fun4all/Fun4AllHistoManager.h>
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/Fun4AllServer.h>
0023 #include <fun4all/PHTFileServer.h>
0024
0025 #include <phool/getClass.h>
0026
0027 #include <HepMC/SimpleVector.h>  // for FourVector
0028
0029 #include <TAxis.h>  // for TAxis
0030 #include <TDatabasePDG.h>
0031 #include <TH1.h>
0032 #include <TH3.h>           // for TH3, TH3D
0033 #include <TNamed.h>        // for TNamed
0034 #include <TParticlePDG.h>  // for TParticlePDG
0035 #include <TVector3.h>
0036
0037 #include <Geant4/G4RunManager.hh>
0038 #include <Geant4/G4ScoringManager.hh>
0039 #include <Geant4/G4String.hh>  // for G4String
0040 #include <Geant4/G4SystemOfUnits.hh>
0041 #include <Geant4/G4THitsMap.hh>     // for G4THitsMap
0042 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0043 #include <Geant4/G4Types.hh>        // for G4int, G4double
0044 #include <Geant4/G4UImanager.hh>
0045 #include <Geant4/G4VScoringMesh.hh>  // for G4VScoringMesh
0046 #include <Geant4/G4Version.hh>
0047
0048 #include <cassert>
0049 #include <cmath>  // for fabs, M_PI
0050 #include <format>
0051 #include <iostream>
0052 #include <limits>   // for numeric_limits
0053 #include <map>      // for _Rb_tree_const_ite...
0054 #include <utility>  // for pair
0055
0056 PHG4ScoringManager::PHG4ScoringManager()
0057   : SubsysReco("PHG4ScoringManager")
0058 {
0059 }
0060
0061 int PHG4ScoringManager::InitRun(PHCompositeNode * /*topNode*/)
0062 {
0063   // 1. check G4RunManager
0064   G4RunManager *runManager = G4RunManager::GetRunManager();
0065   if (runManager == nullptr)
0066   {
0067     std::cout << "PHG4ScoringManager::InitRun - fatal error: G4RunManager was not initialized yet. Please do include the Geant4 simulation in this Fun4All run." << std::endl;
0068     return Fun4AllReturnCodes::ABORTRUN;
0069   }
0070
0071   // 2. Init scoring manager
0072   G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
0073   assert(scoringManager);
0074
0075   // 3. run scoring commands
0076   G4UImanager *UImanager = G4UImanager::GetUIpointer();
0077   assert(UImanager);
0078
0079   for (const std::string &cmd : m_commands)
0080   {
0081     if (Verbosity() >= VERBOSITY_SOME)
0082     {
0083       std::cout << "PHG4ScoringManager::InitRun - execute Geatn4 command: " << cmd << std::endl;
0084     }
0085     UImanager->ApplyCommand(cmd.c_str());
0086   }
0087
0088   // 4 init IOs
0089   if (Verbosity() >= VERBOSITY_SOME)
0090   {
0091     std::cout << "PHG4ScoringManager::InitRun - Making PHTFileServer " << m_outputFileName
0092               << std::endl;
0093   }
0094   PHTFileServer::open(m_outputFileName, "RECREATE");
0095
0096   Fun4AllHistoManager *hm = getHistoManager();
0097   assert(hm);
0098   TH1D *h = new TH1D("hNormalization",  //
0099                      "Normalization;Items;Summed quantity", 10, .5, 10.5);
0100   int i = 1;
0101   h->GetXaxis()->SetBinLabel(i++, "Event count");
0102   h->GetXaxis()->SetBinLabel(i++, "Collision count");
0103   h->GetXaxis()->SetBinLabel(i++, "Event count accepted");
0104   h->GetXaxis()->SetBinLabel(i++, "Collision count accepted");
0105   //  h->GetXaxis()->SetBinLabel(i++, "G4Hit count");
0106   h->GetXaxis()->LabelsOption("v");
0107   hm->registerHisto(h);
0108
0109   hm->registerHisto(new TH1D("hNChEta",  //
0110                              "Charged particle #eta distribution;#eta;Count",
0111                              1000, -5, 5));
0112
0113   hm->registerHisto(new TH1D("hVertexZ",  //
0114                              "Vertex z distribution;z [cm];Count",
0115                              10000, m_vertexHistRange.first, m_vertexHistRange.second));
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 //_________________________________________________________________
0121 void PHG4ScoringManager::G4Command(const std::string &cmd)
0122 {
0123   m_commands.push_back(cmd);
0124   return;
0125 }
0126 //_________________________________________________________________
0127
0128 //_________________________________________________________________
0129 int PHG4ScoringManager::process_event(PHCompositeNode *topNode)
0130 {
0131   Fun4AllHistoManager *hm = getHistoManager();
0132   assert(hm);
0133
0134   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
0135   assert(h_norm);
0136
0137   h_norm->Fill("Event count", 1);
0138
0139   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0140   if (!geneventmap)
0141   {
0142     static bool once = true;
0143     if (once)
0144     {
0145       once = false;
0146       std::cout << "PHG4ScoringManager::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
0147     }
0148   }
0149   else
0150   {
0151     h_norm->Fill("Collision count", geneventmap->size());
0152
0153     TH1D *hVertexZ = dynamic_cast<TH1D *>(hm->getHisto("hVertexZ"));
0154     assert(hVertexZ);
0155
0156     for (const auto genevntpair : geneventmap->get_map())
0157     {
0158       const PHHepMCGenEvent *genevnt = genevntpair.second;
0159       assert(genevnt);
0160
0161       if (genevnt->get_collision_vertex().z() < m_vertexAcceptanceRange.first || genevnt->get_collision_vertex().z() > m_vertexAcceptanceRange.second)
0162       {
0163         if (Verbosity() >= 2)
0164         {
0165           std::cout << __PRETTY_FUNCTION__ << ": get vertex " << genevnt->get_collision_vertex().z()
0166                     << " which is outside range " << m_vertexAcceptanceRange.first << " to " << m_vertexAcceptanceRange.second << " cm:";
0167           genevnt->identify();
0168         }
0169
0170         return Fun4AllReturnCodes::ABORTEVENT;
0171       }
0172
0173       hVertexZ->Fill(genevnt->get_collision_vertex().z());
0174     }
0175
0176     h_norm->Fill("Collision count accepted", geneventmap->size());
0177   }
0178
0179   TH1D *hNChEta = dynamic_cast<TH1D *>(hm->getHisto("hNChEta"));
0180   assert(hNChEta);
0181
0182   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0183   if (!ineve)
0184   {
0185     std::cout << "PHG4ScoringManager::process_event - Error - "
0186               << "unable to find DST node "
0187               << "PHG4INEVENT" << std::endl;
0188   }
0189   else
0190   {
0191     const auto primary_range = ineve->GetParticles();
0192     for (auto particle_iter = primary_range.first;
0193          particle_iter != primary_range.second;
0194          ++particle_iter)
0195     {
0196       const PHG4Particle *p = particle_iter->second;
0197       assert(p);
0198       TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(p->get_pid());
0199       assert(pdg_p);
0200       if (fabs(pdg_p->Charge()) > 0)
0201       {
0202         TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
0203         if (pvec.Perp2() > 0)
0204         {
0205           assert(hNChEta);
0206           hNChEta->Fill(pvec.PseudoRapidity());
0207         }
0208       }
0209     }  //          if (_load_all_particle) else
0210   }
0211
0212   h_norm->Fill("Event count accepted", 1);
0213
0214   return Fun4AllReturnCodes::EVENT_OK;
0215 }
0216
0217 //_________________________________________________________________
0218 int PHG4ScoringManager::End(PHCompositeNode * /*unused*/)
0219 {
0220   if (!m_outputFileName.empty())
0221   {
0222     PHTFileServer::cd(m_outputFileName);
0223     std::cout << "PHG4ScoringManager::End - save results to " << m_outputFileName << std::endl;
0224
0225     makeScoringHistograms();
0226
0227     Fun4AllHistoManager *hm = getHistoManager();
0228     assert(hm);
0229     for (unsigned int i = 0; i < hm->nHistos(); i++)
0230     {
0231       hm->getHisto(i)->Write();
0232     }
0233
0234   }  //   if (not m_outputFileName.empty())
0235
0236   return Fun4AllReturnCodes::EVENT_OK;
0237 }
0238
0239 //! based on G4VScoreWriter::DumpAllQuantitiesToFile()
0240 void PHG4ScoringManager::makeScoringHistograms()
0241 {
0242   Fun4AllHistoManager *hm = getHistoManager();
0243   assert(hm);
0244
0245   G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
0246   assert(scoringManager);
0247
0248   for (unsigned int imesh = 0; imesh < scoringManager->GetNumberOfMesh(); ++imesh)
0249   {
0250     G4VScoringMesh *g4mesh = scoringManager->GetMesh(imesh);
0251     assert(g4mesh);
0252
0253     const std::string meshName(g4mesh->GetWorldName().data());
0254     if (Verbosity())
0255     {
0256       std::cout << "PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName << ": " << std::endl;
0257       g4mesh->List();
0258     }
0259
0260     // descriptors
0261     G4int nMeshSegments[3];  // number of segments of the mesh
0262     g4mesh->GetNumberOfSegments(nMeshSegments);
0263
0264     G4String divisionAxisNames[3];
0265     g4mesh->GetDivisionAxisNames(divisionAxisNames);
0266
0267     // process shape
0268     const G4ThreeVector meshSize = g4mesh->GetSize();
0269     const G4ThreeVector meshTranslate = g4mesh->GetTranslation();
0270 #if G4VERSION_NUMBER >= 1060
0271     const G4VScoringMesh::MeshShape meshShape = g4mesh->GetShape();
0272 #else
0273     const MeshShape meshShape = g4mesh->GetShape();
0274 #endif
0275     // PHENIX units
0276     std::vector<double> meshBoundMin = {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
0277     // PHENIX units
0278     std::vector<double> meshBoundMax = {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
0279 #if G4VERSION_NUMBER >= 1060
0280     if (meshShape == G4VScoringMesh::MeshShape::box)
0281 #else
0282     if (meshShape == boxMesh)
0283 #endif
0284     {
0285       meshBoundMin[0] = (-meshSize[0] + meshTranslate[0]) / cm;
0286       meshBoundMax[0] = (meshSize[0] + meshTranslate[0]) / cm;
0287       meshBoundMin[1] = (-meshSize[1] + meshTranslate[1]) / cm;
0288       meshBoundMax[1] = (meshSize[1] + meshTranslate[1]) / cm;
0289       meshBoundMin[2] = (-meshSize[2] + meshTranslate[2]) / cm;
0290       meshBoundMax[2] = (meshSize[2] + meshTranslate[2]) / cm;
0291
0292       divisionAxisNames[0] += " [cm]";
0293       divisionAxisNames[1] += " [cm]";
0294       divisionAxisNames[2] += " [cm]";
0295     }
0296 #if G4VERSION_NUMBER >= 1060
0297     else if (meshShape == G4VScoringMesh::MeshShape::cylinder)
0298 #else
0299     else if (meshShape == cylinderMesh)
0300 #endif
0301     {
0302       //      fDivisionAxisNames[0] = "Z";
0303       //      fDivisionAxisNames[1] = "PHI";
0304       //      fDivisionAxisNames[2] = "R";
0305       //      G4VSolid * tubsSolid = new G4Tubs(tubsName+"0", // name
0306       //                0.,           // R min
0307       //                fSize[0],     // R max
0308       //                fSize[1],     // Dz
0309       //                0.,           // starting phi
0310       //                                        twopi*rad);   // segment phi
0311       meshBoundMin[0] = (-meshSize[1] + meshTranslate[0]) / cm;
0312       meshBoundMax[0] = (meshSize[1] + meshTranslate[0]) / cm;
0313       meshBoundMin[1] = 0;
0314       meshBoundMax[1] = 2 * M_PI;
0315       meshBoundMin[2] = 0;
0316       meshBoundMax[2] = meshSize[0] / cm;
0317
0318       divisionAxisNames[0] += " [cm]";
0319       divisionAxisNames[1] += " [rad]";
0320       divisionAxisNames[2] += " [cm]";
0321     }
0322     else
0323     {
0324       std::cout << "PHG4ScoringManager::makeScoringHistograms - Error - unsupported mesh shape " << (int) meshShape << ". Skipping this mesh!" << std::endl;
0325       g4mesh->List();
0326       continue;
0327     }
0328
0329 #if G4VERSION_NUMBER >= 1060
0330     G4VScoringMesh::MeshScoreMap fSMap = g4mesh->GetScoreMap();
0331     G4VScoringMesh::MeshScoreMap::const_iterator msMapItr = fSMap.begin();
0332 #else
0333     MeshScoreMap fSMap = g4mesh->GetScoreMap();
0334     MeshScoreMap::const_iterator msMapItr = fSMap.begin();
0335 #endif
0336     for (; msMapItr != fSMap.end(); ++msMapItr)
0337     {
0338       G4String psname = msMapItr->first;
0339 #if G4VERSION_NUMBER >= 1040
0340       std::map<G4int, G4StatDouble *> &score = *(msMapItr->second->GetMap());
0341 #else
0342       std::map<G4int, G4double *> &score = *(msMapItr->second->GetMap());
0343 #endif
0344       G4double unitValue = g4mesh->GetPSUnitValue(psname);
0345       G4String unit = g4mesh->GetPSUnit(psname);
0346
0347       const std::string hname = std::format("hScore_{}_{}", meshName, std::string(psname));
0348       const std::string htitle = std::format("Mesh {}, Primitive scorer {}: score [{}]", meshName, std::string(psname), std::string(unit));
0349
0350       if (Verbosity())
0351       {
0352         std::cout << "PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName
0353                   << "  scorer " << psname
0354                   << "  with axis: "
0355                   << "# i" << divisionAxisNames[0]
0356                   << ", i" << divisionAxisNames[1]
0357                   << ", i" << divisionAxisNames[2]
0358                   << ", value "
0359                   << "[unit: " << unit << "]."
0360                   << " Saving to histogram " << hname << " : " << htitle
0361                   << std::endl;
0362       }
0363       // book histogram
0364       TH3 *h = new TH3D(hname.c_str(),   //
0365                         htitle.c_str(),  //
0366                         nMeshSegments[0],
0367                         meshBoundMin[0], meshBoundMax[0],
0368                         nMeshSegments[1],
0369                         meshBoundMin[1], meshBoundMax[1],
0370                         nMeshSegments[2],
0371                         meshBoundMin[2], meshBoundMax[2]);
0372       hm->registerHisto(h);
0373
0374       h->GetXaxis()->SetTitle(divisionAxisNames[0].data());
0375       h->GetYaxis()->SetTitle(divisionAxisNames[1].data());
0376       h->GetZaxis()->SetTitle(divisionAxisNames[2].data());
0377
0378       // write quantity
0379       for (int x = 0; x < nMeshSegments[0]; x++)
0380       {
0381         for (int y = 0; y < nMeshSegments[1]; y++)
0382         {
0383           for (int z = 0; z < nMeshSegments[2]; z++)
0384           {
0385             const int idx = x * nMeshSegments[1] * nMeshSegments[2] + y * nMeshSegments[2] + z;
0386
0387 #if G4VERSION_NUMBER >= 1040
0388             std::map<G4int, G4StatDouble *>::iterator value = score.find(idx);
0389 #else
0390             std::map<G4int, G4double *>::iterator value = score.find(idx);
0391 #endif
0392             if (value != score.end())
0393             {
0394 #if G4VERSION_NUMBER >= 1040
0395               h->SetBinContent(x + 1, y + 1, z + 1, (value->second->mean()) / unitValue);
0396 #else
0397               h->SetBinContent(x + 1, y + 1, z + 1, *(value->second) / unitValue);
0398 #endif
0399             }
0400
0401           }  //           for (int z = 0; z < fNMeshSegments[2]; z++)
0402         }
0403       }  //       for (int x = 0; x < nMeshSegments[0]; x++)
0404
0405     }  //     for (; msMapItr != fSMap.end(); msMapItr++)
0406
0407   }  //   for (int imesh = 0; imesh < scoringManager->GetNumberOfMesh(); ++imesh)
0408 }
0409
0410 Fun4AllHistoManager *
0411 PHG4ScoringManager::getHistoManager()
0412 {
0413   static std::string histname("PHG4ScoringManager_HISTOS");
0414   Fun4AllServer *se = Fun4AllServer::instance();
0415   Fun4AllHistoManager *hm = se->getHistoManager(histname);
0416   if (!hm)
0417   {
0418     if (Verbosity())
0419     {
0420       std::cout
0421           << "PHG4ScoringManager::get_HistoManager - Making Fun4AllHistoManager " << histname
0422           << std::endl;
0423     }
0424     hm = new Fun4AllHistoManager(histname);
0425     se->registerHistoManager(hm);
0426   }
0427   assert(hm);
0428   return hm;
0429 }