ECCE @ 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 
18 
19 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBO...
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 
60  : SubsysReco("PHG4ScoringManager")
61 {
62 }
63 
65 {
66  //1. check G4RunManager
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;
72  }
73 
74  //2. Init scoring manager
76  assert(scoringManager);
77 
78  //3. run scoring commands
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;
96 
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 
119 }
120 
121 //_________________________________________________________________
122 void PHG4ScoringManager::G4Command(const string &cmd)
123 {
124  m_commands.push_back(cmd);
125  return;
126 }
127 //_________________________________________________________________
128 
129 //_________________________________________________________________
131 {
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 
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 
217 }
218 
219 //_________________________________________________________________
221 {
222  if (not m_outputFileName.empty())
223  {
225  cout << "PHG4ScoringManager::End - save results to " << m_outputFileName << endl;
226 
228 
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 
237 }
238 
241 {
243  assert(hm);
244 
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 
412 {
413  static string histname("PHG4ScoringManager_HISTOS");
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 }