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