Skip to content

Commit

Permalink
Phs source fixed segmentation fault (#328)
Browse files Browse the repository at this point in the history 
* Fixed segmentation faults due to wrong variable type when passing from python to cpp. Was double, should be float

* Enforce correct data type on python side (in case different data types in root file or delivered via uproot

* restructured data transfer from python to cpp

* minor adjustments

* added method to rotate the phase space source in addition to movement required by treatment plan for treatment plan phs source

* fixed wrongly initialized value for entry_start for multithreaded mode

* fixed wrongly initialized value for entry_start for multithreaded mode

* fixed wrongly initialized value for entry_start for multithreaded mode

* fixed initialization for multithreading of phs source

* Update phspsources.py

In multithreading mode
entry_start needs to be an array containing one entry for every thread.
If it is not created by the user,
create a entry_start array with the correct number of start entries
all entries are spaced by the ceil(number of particles/thread)+1

* added improved initialization of starting values for phase space files and sources in multi threading

* Improved phase space source documentation

---------

Co-authored-by: Thomas BAUDIER <thomas.baudier@creatis.insa-lyon.fr>
  • Loading branch information
@GitFuchs @tbaudier
GitFuchs and tbaudier committed May 15, 2024
1 parent 7f4ce0b commit c7174a0
Show file tree
Hide file tree
Showing 7 changed files with 184 additions and 63 deletions.
3 changes: 2 additions & 1 deletion core/opengate_core/g4_bindings/pyG4ParticleTable.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

namespace py = pybind11;

#include "G4IonTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "G4String.hh"
Expand Down Expand Up @@ -107,7 +108,7 @@ void init_G4ParticleTable(py::module &m) {

.def("DumpTable", &G4ParticleTable::DumpTable) //, f_DumpTable())

//.def("GetIonTable", &G4ParticleTable::GetIonTable,
.def("GetIonTable", &G4ParticleTable::GetIonTable)
//...)
//.def("GetShortLivedTable", &G4ParticleTable::GetShortLivedTable,
//...)
Expand Down
63 changes: 45 additions & 18 deletions core/opengate_core/opengate_lib/GatePhaseSpaceSource.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
-------------------------------------------------- */

#include "GatePhaseSpaceSource.h"
#include "G4IonTable.hh"
#include "G4ParticleTable.hh"
#include "G4UnitsTable.hh"
#include "GateHelpersDict.h"
Expand All @@ -17,6 +18,7 @@ GatePhaseSpaceSource::GatePhaseSpaceSource() : GateVSource() {
fCurrentBatchSize = 0;
// fMaxN = 0;
fGlobalFag = false;
fVerbose = false;
}

GatePhaseSpaceSource::~GatePhaseSpaceSource() {
Expand All @@ -37,6 +39,8 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) {
// global (world) or local (mother volume) coordinate system
fGlobalFag = DictGetBool(user_info, "global_flag");

fVerbose = DictGetInt(user_info, "verbose");

// This is done in GateSingleParticleSource, but we need charge/mass later
auto pname = DictGetStr(user_info, "particle");
fParticleTable = G4ParticleTable::GetParticleTable();
Expand All @@ -47,6 +51,9 @@ void GatePhaseSpaceSource::InitializeUserInfo(py::dict &user_info) {
else {
fUseParticleTypeFromFile = false;
fParticleDefinition = fParticleTable->FindParticle(pname);
if (fParticleDefinition == 0) {
Fatal("GatePhaseSpaceSource: PDGCode not found. Aborting.");
}
fCharge = fParticleDefinition->GetPDGCharge();
fMass = fParticleDefinition->GetPDGMass();
}
Expand Down Expand Up @@ -116,6 +123,13 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event,
// check if we should simulate until next primary
// in this case, generate until a second primary is in the list, excluding the
// second primary

// If batch is empty, we generate some particles, could happen at the first
// execution
if (l.fCurrentBatchSize == 0) {
GenerateBatchOfParticles();
}

if (l.fgenerate_until_next_primary) {
int num_primaries = 0;
while (num_primaries <= 2) {
Expand All @@ -129,14 +143,13 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event,

// don't generate the second primary
if (num_primaries < 2) {
// Go
// Create vertex and particle, if called more than once
// multiple vertices and particles exist in the event
// Geant4 will only print the first one, but generate all of them
GenerateOnePrimary(event, current_simulation_time);

// update the index;
// update the root file index;
l.fCurrentIndex++;

// // update the number of generated event
// fNumberOfGeneratedEvents++;
} else
break;
}
Expand All @@ -148,11 +161,10 @@ void GatePhaseSpaceSource::GeneratePrimaries(G4Event *event,
// If batch is empty, we generate some particles
if (l.fCurrentIndex >= l.fCurrentBatchSize)
GenerateBatchOfParticles();

// Go
GenerateOnePrimary(event, current_simulation_time);

// update the index;
// update the root file index;
l.fCurrentIndex++;

// update the number of generated event
Expand Down Expand Up @@ -182,7 +194,6 @@ void GatePhaseSpaceSource::GenerateOnePrimary(G4Event *event,
direction = direction / direction.mag();
direction = ls.fGlobalRotation * direction;
}

// Create the final vertex
AddOnePrimaryVertex(event, position, direction, energy,
current_simulation_time, weight);
Expand All @@ -197,8 +208,18 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event,
auto &l = fThreadLocalDataPhsp.Get();
if (fUseParticleTypeFromFile) {
if (l.fPDGCode[l.fCurrentIndex] != 0) {
// find if particle exists
fParticleDefinition =
fParticleTable->FindParticle(l.fPDGCode[l.fCurrentIndex]);
fParticleTable->FindParticle((int)l.fPDGCode[l.fCurrentIndex]);
// if not, find if it is an ion
if (fParticleDefinition == 0) {
G4IonTable *ionTable = fParticleTable->GetIonTable();
fParticleDefinition =
ionTable->GetIon((G4int)l.fPDGCode[l.fCurrentIndex]);
}
if (fParticleDefinition == 0) {
Fatal("GatePhaseSpaceSource: PDGCode not found. Aborting.");
}
particle->SetParticleDefinition(fParticleDefinition);
} else {
Fatal("GatePhaseSpaceSource: PDGCode not available. Aborting.");
Expand All @@ -215,9 +236,15 @@ void GatePhaseSpaceSource::AddOnePrimaryVertex(G4Event *event,
auto *vertex = new G4PrimaryVertex(position, time);
vertex->SetPrimary(particle);
event->AddPrimaryVertex(vertex);

// weights
event->GetPrimaryVertex(0)->SetWeight(w);
if (fVerbose) {
std::cout << "Particle PDGCode: " << fParticleDefinition->GetPDGEncoding()
<< " Energy: " << energy << " Weight: " << w
<< " Position: " << position << " Direction: " << direction
<< " Time: " << time << " EventID: " << event->GetEventID()
<< std::endl;
}
}

void GatePhaseSpaceSource::SetPDGCodeBatch(
Expand All @@ -227,49 +254,49 @@ void GatePhaseSpaceSource::SetPDGCodeBatch(
}

void GatePhaseSpaceSource::SetEnergyBatch(
const py::array_t<double> &fEnergy) const {
const py::array_t<std::float_t> &fEnergy) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fEnergy = PyBindGetVector(fEnergy);
}

void GatePhaseSpaceSource::SetWeightBatch(
const py::array_t<double> &fWeight) const {
const py::array_t<std::float_t> &fWeight) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fWeight = PyBindGetVector(fWeight);
}

void GatePhaseSpaceSource::SetPositionXBatch(
const py::array_t<double> &fPositionX) const {
const py::array_t<std::float_t> &fPositionX) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fPositionX = PyBindGetVector(fPositionX);
}

void GatePhaseSpaceSource::SetPositionYBatch(
const py::array_t<double> &fPositionY) const {
const py::array_t<std::float_t> &fPositionY) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fPositionY = PyBindGetVector(fPositionY);
}

void GatePhaseSpaceSource::SetPositionZBatch(
const py::array_t<double> &fPositionZ) const {
const py::array_t<std::float_t> &fPositionZ) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fPositionZ = PyBindGetVector(fPositionZ);
}

void GatePhaseSpaceSource::SetDirectionXBatch(
const py::array_t<double> &fDirectionX) const {
const py::array_t<std::float_t> &fDirectionX) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fDirectionX = PyBindGetVector(fDirectionX);
}

void GatePhaseSpaceSource::SetDirectionYBatch(
const py::array_t<double> &fDirectionY) const {
const py::array_t<std::float_t> &fDirectionY) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fDirectionY = PyBindGetVector(fDirectionY);
}

void GatePhaseSpaceSource::SetDirectionZBatch(
const py::array_t<double> &fDirectionZ) const {
const py::array_t<std::float_t> &fDirectionZ) const {
auto &l = fThreadLocalDataPhsp.Get();
l.fDirectionZ = PyBindGetVector(fDirectionZ);
}
Expand Down
40 changes: 21 additions & 19 deletions core/opengate_core/opengate_lib/GatePhaseSpaceSource.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,39 +53,41 @@ class GatePhaseSpaceSource : public GateVSource {

G4ParticleDefinition *fParticleDefinition;
G4ParticleTable *fParticleTable;
double fCharge;
double fMass;

float fCharge;
float fMass;
bool fGlobalFag;
bool fUseParticleTypeFromFile;
bool fVerbose;

// unsigned long fMaxN;
long fNumberOfGeneratedEvents;
size_t fCurrentBatchSize;

void SetPDGCodeBatch(const py::array_t<std::int32_t> &fPDGCode) const;

void SetEnergyBatch(const py::array_t<double> &fEnergy) const;
void SetEnergyBatch(const py::array_t<std::float_t> &fEnergy) const;

void SetWeightBatch(const py::array_t<double> &fWeight) const;
void SetWeightBatch(const py::array_t<std::float_t> &fWeight) const;

void SetPositionXBatch(const py::array_t<double> &fPositionX) const;
void SetPositionXBatch(const py::array_t<std::float_t> &fPositionX) const;

void SetPositionYBatch(const py::array_t<double> &fPositionY) const;
void SetPositionYBatch(const py::array_t<std::float_t> &fPositionY) const;

void SetPositionZBatch(const py::array_t<double> &fPositionZ) const;
void SetPositionZBatch(const py::array_t<std::float_t> &fPositionZ) const;

void SetDirectionXBatch(const py::array_t<double> &fDirectionX) const;
void SetDirectionXBatch(const py::array_t<std::float_t> &fDirectionX) const;

void SetDirectionYBatch(const py::array_t<double> &fDirectionY) const;
void SetDirectionYBatch(const py::array_t<std::float_t> &fDirectionY) const;

void SetDirectionZBatch(const py::array_t<double> &fDirectionZ) const;
void SetDirectionZBatch(const py::array_t<std::float_t> &fDirectionZ) const;

// For MT, all threads local variables are gathered here
struct threadLocalTPhsp {

bool fgenerate_until_next_primary;
int fprimary_PDGCode;
double fprimary_lower_energy_threshold;
float fprimary_lower_energy_threshold;

ParticleGeneratorType fGenerator;
unsigned long fNumberOfGeneratedEvents;
Expand All @@ -94,16 +96,16 @@ class GatePhaseSpaceSource : public GateVSource {

int *fPDGCode;

double *fPositionX;
double *fPositionY;
double *fPositionZ;
float *fPositionX;
float *fPositionY;
float *fPositionZ;

double *fDirectionX;
double *fDirectionY;
double *fDirectionZ;
float *fDirectionX;
float *fDirectionY;
float *fDirectionZ;

double *fEnergy;
double *fWeight;
float *fEnergy;
float *fWeight;
// double * fTime;
};
G4Cache<threadLocalTPhsp> fThreadLocalDataPhsp;
Expand Down
30 changes: 28 additions & 2 deletions docs/source/user_guide_2_2_sources.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ Like all objects, by default, the source is located according to the coordinate

### Phase-Space sources

A phase-space source reads particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example:
A phase-space source reads particles properties (position, direction, energy, etc.) from a root file and use them as events. Typically one particle read is counted as one particle. There is an option to change it, see Enhanced particle counting below. Here is an example to use a phase space source:

```python
source = sim.add_source("PhaseSpaceSource", "phsp_source")
Expand All @@ -152,6 +152,8 @@ In that case, the key "PrePositionLocal" in the root tree file will be used to d

Limitation: the particle timestamps is NOT read from the phsp and not considered (yet)

#### Particle type

The particle type can be set by ```source.particle = "proton"``` option. Using this option all generated particles will be for example protons, overriding the particle type specified in the phase space.

Alternatively, by setting ```source.particle = None``` the particle type is read from the phase space file using the PDGCode. ```source.PDGCode_key = PDGCode``` specifies the name of the entry in the phase space file.
Expand All @@ -171,11 +173,12 @@ gamma: 22
carbon ion C12: 1000060120
```

#### Naming of the phase space file keys
The naming of the phsp file entries generated by e.g. a GATE phase space actor changed over time, most notably from GATE v9 to GATE v10.
Setting ```source.position_key = "PrePositionLocal"``` will cause the phsp source to look for particle positions in ```PrePositionLocal_X, PrePositionLocal_Y, PrePositionLocal_Z```.
```source.direction_key = "PreDirectionLocal"``` will do the corresponding for the particle direction vector components in ```PreDirectionLocal_X, PreDirectionLocal_y, PreDirectionLocal_Z```.

Alternatively, it is possible to directly set the individual components:
It is possible to directly set the individual keys of the phase space file:
```
source.position_key = None"PrePositionLocal"
source.position_key_x = Position_X"
Expand All @@ -185,8 +188,29 @@ source.direction_key = None
source.direction_key_x = Direction_X
source.direction_key_y = Direction_X
source.direction_key_z = Direction_X
source.energy_key = "KineticEnergy"
source.weight_key = "Weight"
source.PDGCode_key = "PDGCode"
```

#### Source rotation and translation

The starting position and direction from eah particle is read from the phase space fiel. It is possible to shift the origin as well as rotate the source.
```
source.translate_position = False
source.rotate_direction = False
source.position.translation = [0, 0, 0]
source.position.rotation = Rotation.identity().as_matrix()
```
If translate_position is set to true, the source.position.translation is evaluated and translates the starting point of the particles by this vector. If rotate_direction is set to true, the source.position.rotation is evaluated to rotate the initial particle vectors. It can be set usign the code below, resulting in a rotation of 30 degrees around the x axis.

```
from scipy.spatial.transform import Rotation
rotation = Rotation.from_euler("x", [30], degrees=True)
source.position.rotation = rotation.as_matrix()
```

#### Enhanced particle counting - realistic particle mix
In case of simulating a realistic particle mix, for example the output after a linac, a phsp file could contain a mixture of particles.
Typically, one would be interested in simulating a given number of primary particles (e.g. protons), simulating, but not counting as secondary particles (e.g. secondary electrons) in the number of particles to simulate.
This can be acchieved by setting ```generate_until_next_primary = True```. Furthermore, the PDG code of the primary particle needs to be specified, as well as a
Expand All @@ -199,6 +223,8 @@ source.primary_lower_energy_threshold = 90.0 * MeV
source.primary_PDGCode = 2212
```

#### Multithreading - where to start reading in a phase space file

For multithread: you need to indicate the ```entry_start``` for all threads, as an array, so that each thread starts in the phsp file at a different position. This done for example as follows (see ```test019_linac_phsp_source_MT.py```). Warning, if the phsp reach its end, it will cycle and start back at the beginning.

```python
Expand Down
Loading

0 comments on commit c7174a0

Please sign in to comment.