Skip to content

Commit 6c41b07

Browse files
dobrienphddobrienphd
committed
Added examples directory and included drop-in code for using ParticleZoo in Geant4 to read phase space files.
1 parent c5a3adf commit 6c41b07

File tree

2 files changed

+311
-0
lines changed

2 files changed

+311
-0
lines changed

‎examples/geant4/ParticleZooPrimaryGeneratorAction.cc‎

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
2+
#include "ParticleZooPrimaryGeneratorAction.hh"
3+
4+
#include "G4GlobalConfig.hh"
5+
#include "G4Threading.hh"
6+
#include "G4Event.hh"
7+
#include "G4RunManager.hh"
8+
#ifdef G4MULTITHREADED
9+
#include "G4MTRunManager.hh"
10+
#endif
11+
12+
#include "particlezoo/utilities/formats.h"
13+
14+
namespace ParticleZoo
15+
{
16+
17+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
18+
19+
G4PrimaryGeneratorAction::G4PrimaryGeneratorAction(const std::string & phaseSpaceFilePath, const std::uint32_t partitionId, const std::uint32_t numberOfPartitions)
20+
: globalTranslation(G4ThreeVector(0,0,0)), partitionId(partitionId), numberOfPartitions(numberOfPartitions)
21+
{
22+
// Initialize the phase space reader
23+
UserOptions userOptions;
24+
phaseSpaceReader = FormatRegistry::CreateReader(phaseSpaceFilePath, userOptions);
25+
26+
// Determine total number of particles and partitioning
27+
std::uint64_t totalNumberOfParticles = phaseSpaceReader->getNumberOfParticles();
28+
std::uint64_t particlesPerPartition = totalNumberOfParticles / numberOfPartitions;
29+
30+
// Determine full range for this partition
31+
std::uint64_t fullRange_startIndex = partitionId * particlesPerPartition;
32+
std::uint64_t fullRange_endIndex = (partitionId == numberOfPartitions - 1) ? totalNumberOfParticles : fullRange_startIndex + particlesPerPartition;
33+
std::uint64_t particlesInThisPartition = fullRange_endIndex - fullRange_startIndex;
34+
35+
// Default to full range
36+
startIndex = fullRange_startIndex;
37+
endIndex = fullRange_endIndex;
38+
39+
// Recycling parameters
40+
recycleNumber = 0;
41+
recycleWeightFactor = 1.0;
42+
43+
// Incremental histories handling
44+
historiesToWait = 0;
45+
historiesWaited = 0;
46+
47+
// Get thread ID
48+
G4int threadId = G4Threading::G4GetThreadId();
49+
50+
#ifdef G4MULTITHREADED
51+
if (threadId >= 0) {
52+
// Configure for multithreading
53+
const G4RunManager* runManager = G4RunManager::GetRunManager();
54+
const G4MTRunManager* mtRunManager = dynamic_cast<const G4MTRunManager*>(runManager);
55+
56+
if (mtRunManager) {
57+
// Divide phase space particles among threads
58+
G4int nThreads = mtRunManager->GetNumberOfThreads();
59+
std::uint64_t particlesPerThread = particlesInThisPartition / nThreads;
60+
startIndex = fullRange_startIndex + threadId * particlesPerThread;
61+
endIndex = (threadId == nThreads - 1) ? fullRange_endIndex : startIndex + particlesPerThread;
62+
63+
if (threadId > 0) { // No need to move for thread 0
64+
// Move the reader to the start index for this thread
65+
phaseSpaceReader->moveToParticle(startIndex);
66+
// Make sure we are at the start of a new history
67+
while (phaseSpaceReader->hasMoreParticles() && !phaseSpaceReader->peekNextParticle().isNewHistory()) {
68+
phaseSpaceReader->getNextParticle();
69+
}
70+
}
71+
}
72+
73+
G4cout << "ParticleZoo::G4PrimaryGeneratorAction: Configured for multithreading. "
74+
<< "Thread ID: " << threadId
75+
<< ", Particle range: [" << startIndex << ", " << endIndex << ")" << G4endl;
76+
} else {
77+
#endif
78+
// Single-threaded configuration
79+
if (startIndex > 0) { // no need to move if starting at 0
80+
// Move to start index
81+
phaseSpaceReader->moveToParticle(startIndex);
82+
// Make sure we are at the start of a new history
83+
while (phaseSpaceReader->hasMoreParticles() && !phaseSpaceReader->peekNextParticle().isNewHistory()) {
84+
phaseSpaceReader->getNextParticle();
85+
}
86+
}
87+
G4cout << "ParticleZoo::G4PrimaryGeneratorAction: Running in single-threaded mode." << G4endl;
88+
#ifdef G4MULTITHREADED
89+
}
90+
#endif
91+
}
92+
93+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94+
95+
G4PrimaryGeneratorAction::~G4PrimaryGeneratorAction()
96+
{
97+
if (phaseSpaceReader) phaseSpaceReader->close();
98+
}
99+
100+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101+
102+
void G4PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
103+
{
104+
// Ensure the reader is valid
105+
if (!phaseSpaceReader) return;
106+
107+
// Handle incremental histories
108+
if (historiesToWait > 0 && historiesWaited < historiesToWait) {
109+
historiesWaited++;
110+
return;
111+
}
112+
113+
// Check if there are more particles available
114+
bool hasMoreParticles = phaseSpaceReader->hasMoreParticles();
115+
116+
// Peek at the next particle to determine incremental histories
117+
if (hasMoreParticles) {
118+
ParticleZoo::Particle particle = phaseSpaceReader->peekNextParticle();
119+
std::int32_t incrementalHistories = particle.getIncrementalHistories();
120+
if (incrementalHistories > historiesToWait + 1) {
121+
historiesToWait = incrementalHistories - 1;
122+
historiesWaited = 0;
123+
return;
124+
} else {
125+
historiesToWait = 0;
126+
historiesWaited = 0;
127+
}
128+
}
129+
130+
// Current global particle index
131+
std::uint64_t currentIndex = startIndex + phaseSpaceReader->getParticlesRead();
132+
133+
// Determine if we can accept the next particle
134+
bool acceptNextParticle = hasMoreParticles && currentIndex < endIndex;
135+
136+
// Generate primary vertices if accepted
137+
if (acceptNextParticle) {
138+
139+
// Lambda to read the next particle and create a G4PrimaryVertex
140+
auto ReadNextVertex = [this](const ParticleZoo::Particle & particle, float weightFactor = 1.0) -> G4PrimaryVertex* {
141+
// Define some unit conversions
142+
constexpr G4double energyUnit = CLHEP::MeV / ParticleZoo::MeV;
143+
constexpr G4double lengthUnit = CLHEP::cm / ParticleZoo::cm;
144+
145+
// Get particle properties
146+
G4int pdgCode = static_cast<G4int>(particle.getPDGCode());
147+
G4double kineticEnergy = static_cast<G4double>(particle.getKineticEnergy()) * energyUnit;
148+
G4double weight = static_cast<G4double>(particle.getWeight()) * weightFactor;
149+
150+
G4double particleX = static_cast<G4double>(particle.getX()) * lengthUnit;
151+
G4double particleY = static_cast<G4double>(particle.getY()) * lengthUnit;
152+
G4double particleZ = static_cast<G4double>(particle.getZ()) * lengthUnit;
153+
154+
G4double directionX = static_cast<G4double>(particle.getDirectionalCosineX());
155+
G4double directionY = static_cast<G4double>(particle.getDirectionalCosineY());
156+
G4double directionZ = static_cast<G4double>(particle.getDirectionalCosineZ());
157+
158+
// Apply global translation
159+
G4ThreeVector position(particleX, particleY, particleZ);
160+
position += this->globalTranslation;
161+
162+
// Create primary vertex and particle
163+
G4PrimaryVertex* vertex = new G4PrimaryVertex(position, 0.0);
164+
G4PrimaryParticle* primary = new G4PrimaryParticle(pdgCode,
165+
directionX,
166+
directionY,
167+
directionZ);
168+
169+
// Set kinetic energy and weight
170+
primary->SetKineticEnergy(kineticEnergy);
171+
primary->SetWeight(weight);
172+
173+
// Set the primary particle to the vertex
174+
vertex->SetPrimary(primary);
175+
176+
// Return the created vertex
177+
return vertex;
178+
};
179+
180+
do {
181+
// Read the next particle from the phase space file
182+
Particle particle = phaseSpaceReader->getNextParticle();
183+
// If recycling is requested, create multiple vertices
184+
for (std::uint32_t r = 0; r <= recycleNumber; r++) {
185+
// Create primary vertex and particle
186+
G4PrimaryVertex* vertex = ReadNextVertex(particle, recycleWeightFactor);
187+
// Add the primary vertex to the event
188+
anEvent->AddPrimaryVertex(vertex);
189+
}
190+
} while (phaseSpaceReader->hasMoreParticles() && !phaseSpaceReader->peekNextParticle().isNewHistory());
191+
192+
} else {
193+
194+
// No more particles available - warn only once
195+
static bool warned = false;
196+
if (!warned) {
197+
G4cout << "No more particles available in phase space file for thread " << G4Threading::G4GetThreadId();
198+
if (numberOfPartitions > 1) {
199+
G4cout << ", partition " << partitionId;
200+
}
201+
G4cout << G4endl;
202+
warned = true;
203+
}
204+
205+
}
206+
}
207+
208+
} // namespace ParticleZoo

‎examples/geant4/ParticleZooPrimaryGeneratorAction.hh‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
#pragma once
2+
3+
#include "G4VUserPrimaryGeneratorAction.hh"
4+
#include "G4Event.hh"
5+
#include "G4ThreeVector.hh"
6+
7+
#include <memory>
8+
#include <string>
9+
#include <particlezoo/PhaseSpaceFileReader.h>
10+
11+
namespace ParticleZoo
12+
{
13+
14+
/**
15+
* @class G4PrimaryGeneratorAction
16+
* @brief Primary generator action for ParticleZoo phase space files.
17+
*
18+
* This class reads particles from a ParticleZoo phase space file and
19+
* generates primary vertices for Geant4 events.
20+
*
21+
* Features:
22+
* - Handles incremental histories
23+
* - Supports multithreading by dividing the particle set among threads (i.e. can be used in multithreaded Geant4 applications)
24+
* - Supports partitioning for applications that run multiple instances of Geant4 to split the workload
25+
* - Particles can be recycled multiple times with adjusted weights
26+
*/
27+
class G4PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
28+
{
29+
public:
30+
/**
31+
* @brief Constructor.
32+
* @param phaseSpaceFilePath Path to the ParticleZoo phase space file.
33+
* @param translation Global translation to apply to all particle positions (default is no translation).
34+
* @param partitionId The partition ID for this instance (used for splitting workload across multiple application instances, default is 0).
35+
* @param numberOfPartitions The total number of partitions (used for splitting workload across multiple application instances, default is 1).
36+
*/
37+
G4PrimaryGeneratorAction(const std::string & phaseSpaceFilePath, const std::uint32_t partitionId = 0, const std::uint32_t numberOfPartitions = 1);
38+
~G4PrimaryGeneratorAction() override;
39+
40+
/**
41+
* @brief Generate primary vertices for the given event.
42+
* @param anEvent Pointer to the Geant4 event.
43+
*/
44+
void GeneratePrimaries(G4Event* anEvent) override;
45+
46+
/**
47+
* @brief Set the global translation to apply to all particle positions.
48+
* @param translation The translation vector.
49+
*/
50+
void SetTranslation(const G4ThreeVector & translation);
51+
52+
/**
53+
* @brief Set the number of times to recycle each particle.
54+
* @param n Number of times to recycle each particle.
55+
*/
56+
void SetRecycleNumber(std::uint32_t n);
57+
58+
private:
59+
// Phase space file reader
60+
std::unique_ptr<ParticleZoo::PhaseSpaceFileReader> phaseSpaceReader;
61+
62+
// Global translation to apply to all particle positions
63+
G4ThreeVector globalTranslation;
64+
65+
// Particle range for this generator action
66+
std::uint64_t startIndex;
67+
std::uint64_t endIndex;
68+
69+
// Phase-space partitioning information for multiple instances
70+
// for applications that don't use Geant4 multithreading but
71+
// want to split the workload across multiple G4 runs.
72+
std::uint32_t partitionId;
73+
std::uint32_t numberOfPartitions;
74+
75+
// Recycling parameters
76+
std::uint32_t recycleNumber;
77+
G4double recycleWeightFactor;
78+
79+
// Incremental histories handling
80+
std::int32_t historiesToWait;
81+
std::int32_t historiesWaited;
82+
};
83+
84+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85+
86+
inline void G4PrimaryGeneratorAction::SetTranslation(const G4ThreeVector & translation)
87+
{
88+
globalTranslation = translation;
89+
G4cout << "ParticleZoo::G4PrimaryGeneratorAction: Set global translation to "
90+
<< globalTranslation << G4endl;
91+
}
92+
93+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94+
95+
inline void G4PrimaryGeneratorAction::SetRecycleNumber(std::uint32_t n)
96+
{
97+
recycleNumber = n;
98+
recycleWeightFactor = 1.0 / static_cast<G4double>(recycleNumber + 1);
99+
G4cout << "ParticleZoo::G4PrimaryGeneratorAction: Set recycle number to "
100+
<< recycleNumber << ", weight factor: " << recycleWeightFactor << G4endl;
101+
}
102+
103+
}; // namespace ParticleZoo

0 commit comments

Comments
 (0)