Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ScoringManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ScoringManager.cc
1 // $Id: $
2 
11 #include "PHG4ScoringManager.h"
12 
13 #include "PHG4InEvent.h"
14 #include "PHG4Particle.h"
15 
16 #include <phhepmc/PHHepMCGenEvent.h>
17 #include <phhepmc/PHHepMCGenEventMap.h>
18 
19 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBO...
20 #include <fun4all/Fun4AllHistoManager.h>
21 #include <fun4all/Fun4AllReturnCodes.h>
22 #include <fun4all/Fun4AllServer.h>
23 #include <fun4all/PHTFileServer.h>
24 
25 #include <phool/getClass.h>
26 
27 #include <HepMC/SimpleVector.h> // for FourVector
28 
29 #include <TAxis.h> // for TAxis
30 #include <TDatabasePDG.h>
31 #include <TH1.h>
32 #include <TH3.h> // for TH3, TH3D
33 #include <TNamed.h> // for TNamed
34 #include <TParticlePDG.h> // for TParticlePDG
35 #include <TVector3.h>
36 
37 #include <Geant4/G4RunManager.hh>
38 #include <Geant4/G4ScoringManager.hh>
39 #include <Geant4/G4String.hh> // for G4String
40 #include <Geant4/G4SystemOfUnits.hh>
41 #include <Geant4/G4THitsMap.hh> // for G4THitsMap
42 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
43 #include <Geant4/G4Types.hh> // for G4int, G4double
44 #include <Geant4/G4UImanager.hh>
45 #include <Geant4/G4VScoringMesh.hh> // for G4VScoringMesh
46 #include <Geant4/G4Version.hh>
47 
48 #pragma GCC diagnostic push
49 #pragma GCC diagnostic ignored "-Wshadow"
50 #include <boost/format.hpp>
51 #pragma GCC diagnostic pop
52 
53 #include <cassert>
54 #include <cmath> // for fabs, M_PI
55 #include <iostream>
56 #include <limits> // for numeric_limits
57 #include <map> // for _Rb_tree_const_ite...
58 #include <utility> // for pair
59 
60 using namespace std;
61 
62 PHG4ScoringManager::PHG4ScoringManager()
63  : SubsysReco("PHG4ScoringManager")
64 {
65 }
66 
67 int PHG4ScoringManager::InitRun(PHCompositeNode */*topNode*/)
68 {
69  //1. check G4RunManager
70  G4RunManager *runManager = G4RunManager::GetRunManager();
71  if (runManager == nullptr)
72  {
73  cout << "PHG4ScoringManager::InitRun - fatal error: G4RunManager was not initialized yet. Please do include the Geant4 simulation in this Fun4All run." << endl;
74  return Fun4AllReturnCodes::ABORTRUN;
75  }
76 
77  //2. Init scoring manager
78  G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
79  assert(scoringManager);
80 
81  //3. run scoring commands
82  G4UImanager *UImanager = G4UImanager::GetUIpointer();
83  assert(UImanager);
84 
85  for (const string &cmd : m_commands)
86  {
87  if (Verbosity() >= VERBOSITY_SOME)
88  {
89  cout << "PHG4ScoringManager::InitRun - execute Geatn4 command: " << cmd << endl;
90  }
91  UImanager->ApplyCommand(cmd.c_str());
92  }
93 
94  //4 init IOs
95  if (Verbosity() >= VERBOSITY_SOME)
96  cout << "PHG4ScoringManager::InitRun - Making PHTFileServer " << m_outputFileName
97  << endl;
98  PHTFileServer::get().open(m_outputFileName, "RECREATE");
99 
100  Fun4AllHistoManager *hm = getHistoManager();
101  assert(hm);
102  TH1D *h = new TH1D("hNormalization", //
103  "Normalization;Items;Summed quantity", 10, .5, 10.5);
104  int i = 1;
105  h->GetXaxis()->SetBinLabel(i++, "Event count");
106  h->GetXaxis()->SetBinLabel(i++, "Collision count");
107  h->GetXaxis()->SetBinLabel(i++, "Event count accepted");
108  h->GetXaxis()->SetBinLabel(i++, "Collision count accepted");
109  // h->GetXaxis()->SetBinLabel(i++, "G4Hit count");
110  h->GetXaxis()->LabelsOption("v");
111  hm->registerHisto(h);
112 
113  hm->registerHisto(new TH1D("hNChEta", //
114  "Charged particle #eta distribution;#eta;Count",
115  1000, -5, 5));
116 
117  hm->registerHisto(new TH1D("hVertexZ", //
118  "Vertex z distribution;z [cm];Count",
119  10000, m_vertexHistRange.first, m_vertexHistRange.second));
120 
121  return Fun4AllReturnCodes::EVENT_OK;
122 }
123 
124 //_________________________________________________________________
125 void PHG4ScoringManager::G4Command(const string &cmd)
126 {
127  m_commands.push_back(cmd);
128  return;
129 }
130 //_________________________________________________________________
131 
132 //_________________________________________________________________
133 int PHG4ScoringManager::process_event(PHCompositeNode *topNode)
134 {
135  Fun4AllHistoManager *hm = getHistoManager();
136  assert(hm);
137 
138  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
139  assert(h_norm);
140 
141  h_norm->Fill("Event count", 1);
142 
143  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
144  if (!geneventmap)
145  {
146  static bool once = true;
147  if (once)
148  {
149  once = false;
150  cout << "PHG4ScoringManager::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
151  }
152  }
153  else
154  {
155  h_norm->Fill("Collision count", geneventmap->size());
156 
157  TH1D *hVertexZ = dynamic_cast<TH1D *>(hm->getHisto("hVertexZ"));
158  assert(hVertexZ);
159 
160  for (const auto genevntpair : geneventmap->get_map())
161  {
162  const PHHepMCGenEvent *genevnt = genevntpair.second;
163  assert(genevnt);
164 
165  if (genevnt->get_collision_vertex().z() < m_vertexAcceptanceRange.first
166  or genevnt->get_collision_vertex().z() > m_vertexAcceptanceRange.second)
167  {
168  if (Verbosity() >= 2)
169  {
170  cout <<__PRETTY_FUNCTION__<<": get vertex "<<genevnt->get_collision_vertex().z()
171  <<" which is outside range "<<m_vertexAcceptanceRange.first <<" to "<<m_vertexAcceptanceRange.second<<" cm:";
172  genevnt->identify();
173  }
174 
175  return Fun4AllReturnCodes::ABORTEVENT;
176  }
177 
178  hVertexZ->Fill(genevnt->get_collision_vertex().z());
179  }
180 
181  h_norm->Fill("Collision count accepted", geneventmap->size());
182  }
183 
184  TH1D *hNChEta = dynamic_cast<TH1D *>(hm->getHisto("hNChEta"));
185  assert(hNChEta);
186 
187  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
188  if (!ineve)
189  {
190  cout << "PHG4ScoringManager::process_event - Error - "
191  << "unable to find DST node "
192  << "PHG4INEVENT" << endl;
193  }
194  else
195  {
196  const auto primary_range = ineve->GetParticles();
197  for (auto particle_iter = primary_range.first;
198  particle_iter != primary_range.second;
199  ++particle_iter)
200  {
201  const PHG4Particle *p = particle_iter->second;
202  assert(p);
203  TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(p->get_pid());
204  assert(pdg_p);
205  if (fabs(pdg_p->Charge()) > 0)
206  {
207  TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
208  if (pvec.Perp2() > 0)
209  {
210  assert(hNChEta);
211  hNChEta->Fill(pvec.PseudoRapidity());
212  }
213  }
214  } // if (_load_all_particle) else
215  }
216 
217  h_norm->Fill("Event count accepted", 1);
218 
219  return Fun4AllReturnCodes::EVENT_OK;
220 }
221 
222 //_________________________________________________________________
223 int PHG4ScoringManager::End(PHCompositeNode *)
224 {
225  if (not m_outputFileName.empty())
226  {
227  PHTFileServer::get().cd(m_outputFileName);
228  cout << "PHG4ScoringManager::End - save results to " << m_outputFileName << endl;
229 
230  makeScoringHistograms();
231 
232  Fun4AllHistoManager *hm = getHistoManager();
233  assert(hm);
234  for (unsigned int i = 0; i < hm->nHistos(); i++)
235  hm->getHisto(i)->Write();
236 
237  } // if (not m_outputFileName.empty())
238 
239  return Fun4AllReturnCodes::EVENT_OK;
240 }
241 
243 void PHG4ScoringManager::makeScoringHistograms()
244 {
245  Fun4AllHistoManager *hm = getHistoManager();
246  assert(hm);
247 
248  G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
249  assert(scoringManager);
250 
251  for (unsigned int imesh = 0; imesh < scoringManager->GetNumberOfMesh(); ++imesh)
252  {
253  G4VScoringMesh *g4mesh = scoringManager->GetMesh(imesh);
254  assert(g4mesh);
255 
256  const string meshName(g4mesh->GetWorldName().data());
257  if (Verbosity())
258  {
259  cout << "PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName << ": " << endl;
260  g4mesh->List();
261  }
262 
263  // descriptors
264  G4int nMeshSegments[3]; // number of segments of the mesh
265  g4mesh->GetNumberOfSegments(nMeshSegments);
266 
267  G4String divisionAxisNames[3];
268  g4mesh->GetDivisionAxisNames(divisionAxisNames);
269 
270  // process shape
271  const G4ThreeVector meshSize = g4mesh->GetSize();
272  const G4ThreeVector meshTranslate = g4mesh->GetTranslation();
273 #if G4VERSION_NUMBER >= 1060
274  const G4VScoringMesh::MeshShape meshShape = g4mesh->GetShape();
275 #else
276  const MeshShape meshShape = g4mesh->GetShape();
277 #endif
278  //PHENIX units
279  vector<double> meshBoundMin = {std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN()};
280  //PHENIX units
281  vector<double> meshBoundMax = {std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN()};
282 #if G4VERSION_NUMBER >= 1060
283  if (meshShape == G4VScoringMesh::MeshShape::box)
284 #else
285  if (meshShape == boxMesh)
286 #endif
287  {
288  meshBoundMin[0] = (-meshSize[0] + meshTranslate[0]) / cm;
289  meshBoundMax[0] = (meshSize[0] + meshTranslate[0]) / cm;
290  meshBoundMin[1] = (-meshSize[1] + meshTranslate[1]) / cm;
291  meshBoundMax[1] = (meshSize[1] + meshTranslate[1]) / cm;
292  meshBoundMin[2] = (-meshSize[2] + meshTranslate[2]) / cm;
293  meshBoundMax[2] = (meshSize[2] + meshTranslate[2]) / cm;
294 
295  divisionAxisNames[0] += " [cm]";
296  divisionAxisNames[1] += " [cm]";
297  divisionAxisNames[2] += " [cm]";
298  }
299 #if G4VERSION_NUMBER >= 1060
300  else if (meshShape == G4VScoringMesh::MeshShape::cylinder)
301 #else
302  else if (meshShape == cylinderMesh)
303 #endif
304  {
305  // fDivisionAxisNames[0] = "Z";
306  // fDivisionAxisNames[1] = "PHI";
307  // fDivisionAxisNames[2] = "R";
308  // G4VSolid * tubsSolid = new G4Tubs(tubsName+"0", // name
309  // 0., // R min
310  // fSize[0], // R max
311  // fSize[1], // Dz
312  // 0., // starting phi
313  // twopi*rad); // segment phi
314  meshBoundMin[0] = (-meshSize[1] + meshTranslate[0]) / cm;
315  meshBoundMax[0] = (meshSize[1] + meshTranslate[0]) / cm;
316  meshBoundMin[1] = 0;
317  meshBoundMax[1] = 2 * M_PI;
318  meshBoundMin[2] = 0;
319  meshBoundMax[2] = meshSize[0] / cm;
320 
321  divisionAxisNames[0] += " [cm]";
322  divisionAxisNames[1] += " [rad]";
323  divisionAxisNames[2] += " [cm]";
324  }
325  else
326  {
327  cout << "PHG4ScoringManager::makeScoringHistograms - Error - unsupported mesh shape " << (int) meshShape << ". Skipping this mesh!" << endl;
328  g4mesh->List();
329  continue;
330  }
331 
332 #if G4VERSION_NUMBER >= 1060
333  G4VScoringMesh::MeshScoreMap fSMap = g4mesh->GetScoreMap();
334  G4VScoringMesh::MeshScoreMap::const_iterator msMapItr = fSMap.begin();
335 #else
336  MeshScoreMap fSMap = g4mesh->GetScoreMap();
337  MeshScoreMap::const_iterator msMapItr = fSMap.begin();
338 #endif
339  for (; msMapItr != fSMap.end(); ++msMapItr)
340  {
341  G4String psname = msMapItr->first;
342 #if G4VERSION_NUMBER >= 1040
343  std::map<G4int, G4StatDouble *> &score = *(msMapItr->second->GetMap());
344 #else
345  std::map<G4int, G4double *> &score = *(msMapItr->second->GetMap());
346 #endif
347  G4double unitValue = g4mesh->GetPSUnitValue(psname);
348  G4String unit = g4mesh->GetPSUnit(psname);
349 
350  const string hname = boost::str(boost::format("hScore_%1%_%2%") % meshName.data() % psname.data());
351  const string htitle = boost::str(boost::format("Mesh %1%, Primitive scorer %2%: score [%3%]") % meshName.c_str() % psname.data() % unit.data());
352 
353  if (Verbosity())
354  {
355  cout << "PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName
356  << " scorer " << psname
357  << " with axis: "
358  << "# i" << divisionAxisNames[0]
359  << ", i" << divisionAxisNames[1]
360  << ", i" << divisionAxisNames[2]
361  << ", value "
362  << "[unit: " << unit << "]."
363  << " Saving to histogram " << hname << " : " << htitle
364  << endl;
365  }
366  //book histogram
367  TH3 *h = new TH3D(hname.c_str(), //
368  htitle.c_str(), //
369  nMeshSegments[0],
370  meshBoundMin[0], meshBoundMax[0],
371  nMeshSegments[1],
372  meshBoundMin[1], meshBoundMax[1],
373  nMeshSegments[2],
374  meshBoundMin[2], meshBoundMax[2]);
375  hm->registerHisto(h);
376 
377  h->GetXaxis()->SetTitle(divisionAxisNames[0].data());
378  h->GetYaxis()->SetTitle(divisionAxisNames[1].data());
379  h->GetZaxis()->SetTitle(divisionAxisNames[2].data());
380 
381  // write quantity
382  for (int x = 0; x < nMeshSegments[0]; x++)
383  {
384  for (int y = 0; y < nMeshSegments[1]; y++)
385  {
386  for (int z = 0; z < nMeshSegments[2]; z++)
387  {
388  const int idx = x * nMeshSegments[1] * nMeshSegments[2] + y * nMeshSegments[2] + z;
389 
390 #if G4VERSION_NUMBER >= 1040
391  std::map<G4int, G4StatDouble *>::iterator value = score.find(idx);
392 #else
393  std::map<G4int, G4double *>::iterator value = score.find(idx);
394 #endif
395  if (value != score.end())
396  {
397 #if G4VERSION_NUMBER >= 1040
398  h->SetBinContent(x + 1, y + 1, z + 1, (value->second->mean()) / unitValue);
399 #else
400  h->SetBinContent(x + 1, y + 1, z + 1, *(value->second) / unitValue);
401 #endif
402  }
403 
404  } // for (int z = 0; z < fNMeshSegments[2]; z++)
405  }
406  } // for (int x = 0; x < nMeshSegments[0]; x++)
407 
408  } // for (; msMapItr != fSMap.end(); msMapItr++)
409 
410  } // for (int imesh = 0; imesh < scoringManager->GetNumberOfMesh(); ++imesh)
411 }
412 
413 Fun4AllHistoManager *
414 PHG4ScoringManager::getHistoManager()
415 {
416  static string histname("PHG4ScoringManager_HISTOS");
417  Fun4AllServer *se = Fun4AllServer::instance();
418  Fun4AllHistoManager *hm = se->getHistoManager(histname);
419  if (!hm)
420  {
421  if (Verbosity())
422  cout
423  << "PHG4ScoringManager::get_HistoManager - Making Fun4AllHistoManager " << histname
424  << endl;
425  hm = new Fun4AllHistoManager(histname);
426  se->registerHistoManager(hm);
427  }
428  assert(hm);
429  return hm;
430 }