27 #include <HepMC/SimpleVector.h>
30 #include <TDatabasePDG.h>
34 #include <TParticlePDG.h>
37 #include <Geant4/G4RunManager.hh>
38 #include <Geant4/G4ScoringManager.hh>
39 #include <Geant4/G4String.hh>
40 #include <Geant4/G4SystemOfUnits.hh>
41 #include <Geant4/G4THitsMap.hh>
42 #include <Geant4/G4ThreeVector.hh>
43 #include <Geant4/G4Types.hh>
44 #include <Geant4/G4UImanager.hh>
45 #include <Geant4/G4VScoringMesh.hh>
46 #include <Geant4/G4Version.hh>
48 #include <boost/format.hpp>
68 if (runManager ==
nullptr)
70 cout <<
"PHG4ScoringManager::InitRun - fatal error: G4RunManager was not initialized yet. Please do include the Geant4 simulation in this Fun4All run." << endl;
76 assert(scoringManager);
86 cout <<
"PHG4ScoringManager::InitRun - execute Geatn4 command: " << cmd << endl;
93 cout <<
"PHG4ScoringManager::InitRun - Making PHTFileServer " <<
m_outputFileName
99 TH1D *
h =
new TH1D(
"hNormalization",
100 "Normalization;Items;Summed quantity", 10, .5, 10.5);
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");
107 h->GetXaxis()->LabelsOption(
"v");
111 "Charged particle #eta distribution;#eta;Count",
115 "Vertex z distribution;z [cm];Count",
135 TH1D *h_norm =
dynamic_cast<TH1D *
>(hm->
getHisto(
"hNormalization"));
138 h_norm->Fill(
"Event count", 1);
140 PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
143 static bool once =
true;
147 cout <<
"PHG4ScoringManager::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
152 h_norm->Fill(
"Collision count", geneventmap->
size());
154 TH1D *hVertexZ =
dynamic_cast<TH1D *
>(hm->
getHisto(
"hVertexZ"));
157 for (
const auto genevntpair : geneventmap->
get_map())
178 h_norm->Fill(
"Collision count accepted", geneventmap->
size());
181 TH1D *hNChEta =
dynamic_cast<TH1D *
>(hm->
getHisto(
"hNChEta"));
184 PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode,
"PHG4INEVENT");
187 cout <<
"PHG4ScoringManager::process_event - Error - "
188 <<
"unable to find DST node "
189 <<
"PHG4INEVENT" << endl;
194 for (
auto particle_iter = primary_range.first;
195 particle_iter != primary_range.second;
200 TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(p->
get_pid());
202 if (fabs(pdg_p->Charge()) > 0)
205 if (pvec.Perp2() > 0)
208 hNChEta->Fill(pvec.PseudoRapidity());
214 h_norm->Fill(
"Event count accepted", 1);
225 cout <<
"PHG4ScoringManager::End - save results to " <<
m_outputFileName << endl;
231 for (
unsigned int i = 0; i < hm->
nHistos(); i++)
246 assert(scoringManager);
248 for (
unsigned int imesh = 0; imesh < scoringManager->
GetNumberOfMesh(); ++imesh)
256 cout <<
"PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName <<
": " << endl;
261 G4int nMeshSegments[3];
270 #if G4VERSION_NUMBER >= 1060
273 const MeshShape meshShape = g4mesh->
GetShape();
276 vector<double> meshBoundMin = {std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN(), std::numeric_limits<double>::signaling_NaN()};
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
282 if (meshShape == boxMesh)
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;
292 divisionAxisNames[0] +=
" [cm]";
293 divisionAxisNames[1] +=
" [cm]";
294 divisionAxisNames[2] +=
" [cm]";
296 #if G4VERSION_NUMBER >= 1060
297 else if (meshShape == G4VScoringMesh::MeshShape::cylinder)
299 else if (meshShape == cylinderMesh)
311 meshBoundMin[0] = (-meshSize[1] + meshTranslate[0]) /
cm;
312 meshBoundMax[0] = (meshSize[1] + meshTranslate[0]) /
cm;
314 meshBoundMax[1] = 2 *
M_PI;
316 meshBoundMax[2] = meshSize[0] /
cm;
318 divisionAxisNames[0] +=
" [cm]";
319 divisionAxisNames[1] +=
" [rad]";
320 divisionAxisNames[2] +=
" [cm]";
324 cout <<
"PHG4ScoringManager::makeScoringHistograms - Error - unsupported mesh shape " << (
int) meshShape <<
". Skipping this mesh!" << endl;
329 #if G4VERSION_NUMBER >= 1060
331 G4VScoringMesh::MeshScoreMap::const_iterator msMapItr = fSMap.begin();
334 MeshScoreMap::const_iterator msMapItr = fSMap.begin();
336 for (; msMapItr != fSMap.end(); ++msMapItr)
339 #if G4VERSION_NUMBER >= 1040
340 std::map<G4int, G4StatDouble *> &score = *(msMapItr->second->GetMap());
342 std::map<G4int, G4double *> &score = *(msMapItr->second->GetMap());
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());
352 cout <<
"PHG4ScoringManager::makeScoringHistograms - processing mesh " << meshName
353 <<
" scorer " << psname
355 <<
"# i" << divisionAxisNames[0]
356 <<
", i" << divisionAxisNames[1]
357 <<
", i" << divisionAxisNames[2]
359 <<
"[unit: " << unit <<
"]."
360 <<
" Saving to histogram " << hname <<
" : " << htitle
364 TH3 *
h =
new TH3D(hname.c_str(),
367 meshBoundMin[0], meshBoundMax[0],
369 meshBoundMin[1], meshBoundMax[1],
371 meshBoundMin[2], meshBoundMax[2]);
374 h->GetXaxis()->SetTitle(divisionAxisNames[0].
data());
375 h->GetYaxis()->SetTitle(divisionAxisNames[1].
data());
376 h->GetZaxis()->SetTitle(divisionAxisNames[2].
data());
379 for (
int x = 0;
x < nMeshSegments[0];
x++)
381 for (
int y = 0;
y < nMeshSegments[1];
y++)
383 for (
int z = 0;
z < nMeshSegments[2];
z++)
385 const int idx =
x * nMeshSegments[1] * nMeshSegments[2] +
y * nMeshSegments[2] +
z;
387 #if G4VERSION_NUMBER >= 1040
388 std::map<G4int, G4StatDouble *>::iterator
value = score.find(idx);
390 std::map<G4int, G4double *>::iterator value = score.find(idx);
392 if (value != score.end())
394 #if G4VERSION_NUMBER >= 1040
395 h->SetBinContent(
x + 1,
y + 1, z + 1, (value->second->mean()) / unitValue);
397 h->SetBinContent(
x + 1,
y + 1, z + 1, *(value->second) / unitValue);
413 static string histname(
"PHG4ScoringManager_HISTOS");
420 <<
"PHG4ScoringManager::get_HistoManager - Making Fun4AllHistoManager " << histname