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