Radiotherapy application

The Radiotherapy application contains some utilities for the simulation of a teletherapy linear accelerator and dose calculations in the patient. Many of the utilities can also be useful in the simulation of brachytherapy treatments.

Geometrical modules

Two geometrical modules are defined related to Radiotherapy: JAWS and MLC (multi-leaf collimators). They can be defined in the usual geometry text file (see chapter on Geometry) with the format described below.

JAWS module

This module serves to describe the jaws of a radiotherapy accelerator. The jaws are described by their dimensions, and their rotations are calculated automatically using the parameters of the projection into a plane. It is important to remember that the beam is assumed to run along the negative Z axis, following the IEC 61217 convention (if you prefer to define the beam along the positive Z axis, you have to define the environmental variable BEAMZPOS and recompile GAMOS).

The endleaf straight implies a circular leaf shape, and the movement is made along a circle. The endleaf rounded implies an straight leaf shape and the movement is made along a straight line. The rounded type has a sub-type, when the center of the circle that shapes the rounded leaf tip is not at the centre of the leaf, so that the tip is not symmetrical with respect to the tip half plane. And the displaced rounded tip has in turn a sub-type, which allows to define the tip positions as a function of the field aperture reading the positions from a file, instead of calculating them automatically (see below).

The endleaf cross positions are described by their projection at the isocenter plane, as it is common in the DICOM RTPlan files. For the straight leafs, this position is calculated so that the line between the focus and the projection at the isocenter plane is parallel to the leaf end. For the rounded leafs, this position is calculated so that this line is tangent to the leaf end. If a half-value layer is defined this line traverses the leaf end by a distance equals to this value. Alternative, as a function of the field aperture reading the positions from a file, instead of calculating them automatically.

The following parameters have to be provided (between parenthesis the name of each parameters as they appear in fig. RT1):

:MODULE JAWS

Module name

Module orientation: direction of the jaws. It has to be X or Y.

End leaf type It must be STRAIGHT, ROUND, ROUND_DISP. ROUND_DISP means that the centre of the circle that defines the profile is not at the line in the middle of the leaf. It must be STRAIGHT, ROUND, ROUND_DISP, ROUND_DISP_POSFILE. STRAIGHT means that the leaf profile is a rectangle and it moves along a circle, while ROUND means that the leaf profile tip closest to the beam has circular shape and it moves along a straight line. ROUND_DISP means that the centre of the circle that defines the profile is not at the line (perpendicular to the beam axis) in the middle of the leaf. And ROUND_DISP_POSFILE means that the endleaf positions are a function of the field aperture and are read from a file, instead of calculating them automatically.

  • X half dimension
  • Y half dimension (LONGDIM)
  • Z half dimension (ZDIM)
  • If the end leaf type is ROUND or ROUND_DISP, end leaf radius (TIP_RADIUS)
  • If the end leaf type is ROUND_DISP, distance between leaf top and radius centre (DISP).

If the end leaf type is ROUND, the radius center is at the middle of the leaf, therefore this parameter is not needed.

  • If the end leaf type is ROUND or ROUND_DISP, half value layer (HVL).

This half value layer should correspond to the distance at which only half of the gammas from the accelerator traverse the material of the leaf. It is used to displace the leafs so that the line from the focus to the aperture position in the isocenter given below traverses a distance of the rounded leaf tip equal to this value.

  • Z position of the focus (Z_FOCUS)
  • If the end leaf type is STRAIGHT, radius of rotation of the aperture movement (usually the same as the distance between the focus and the volume centre)
  • Z of the position where the center of the volume would be if the field were 0 (no rotation) (Z_CENTER)
  • Z position of the isocenter, where the projections are calculated (Z_ISOCENTER)
  • Field negative projection (usually a negative value) (FIELD)
  • Field positive projection (usually a positive value)
  • Material
  • Parent volume name
../_images/jaws.jpg

Fig. RT1. XZ plane view of a cut along Y=40 mm of the jaws created from the file below

Module Jaws example

On the next lines we include a simplified geometrical description, which can be used as example. On figure RT1 it can we seen the corresponding representation.

JAWS_Y Module_name

Y ROUND_DISP Orientation, Tip type

:P JAWS_ZPOS -150.

:P JAWS_Z 100./2.

:P JAWS_TIP_CIRCLE_Z 40.

:P JAWS_TIP_CIRCLE_R 150.

:P JAWS_HVL 70.

:P JAWS_FIELD 5.*cm/2.

:MODULE JAWS

10.*cm/2 10.*cm/2 $JAWS_Z // Half-dimensions

$JAWS_TIP_CIRCLE_R $JAWS_TIP_CIRCLE_Z // Tip_circle_radius Tip_circle_centre_Z

$JAWS_HVL // Half value layer

  1. $JAWS_ZPOS-$JAWS_Z -30.*cm // Z focus Z_centre Z_isocentre(field)

-$JAWS_FIELD $JAWS_FIELD // Field apertures: RIGHT & LEFT

RTUW ACCEL // Material Mother_volume

Projection in the isocenter plane of the light field, XL(i), and leaf tip position, W(i)

MLC module

Introduction

A MLC (Multi-Leaf Collimator) is a collimator system composed of a set of leaf pairs usually a few millimeters thick which can be open or close to form a hole in the field of view of an irradiation treatment beam. By opening the different leaves a different distance an irregular treatment field can be shaped.

This module can describe a high variety of MLC models, whether the leaf profiles are focused to a point in the Z axis or to an off-axis point.

Two endleaf types are supported: straight or rounded. See Module Jaws for details.

Any leaf cross profile (i.e. the 2D profile in the direction orthogonal to the leaf movement) shape can be described by enumerating the list of 2-dimensional points that describe it. These points have to be given in clock-wise order. When the leaf is placed, the profile shape is deformed so the lines that are not perpendicular to the Z axis are focused towards the focus point; this has to be taken into account so that you do not define lines not perpendicular if you do not want that they are modified to point to the focus).

Geometrical module definition and variables list

The MLC module description has to start with the two words (words have to be separated by blank spaces, and follow the order described here, but there is no constraint on the number of lines they may occupy). The following parameters have to be provided (between parenthesis the name of each parameters as they appears in fig. RT2 and RT3):

:MODULE MLC

Then the following parameters have to be filled

- Module name

- Module type. The type of MLC has to be FOCUSED, what means that all leaves profiles in the cross place are focused to a point. Type UNFOCUSED will be supported in future releases.

- Leaves orientation. Orientation of the leaves, i.e. the axis of movement, can be X or Y

- Leaf tip type. It must be STRAIGHT, ROUND, ROUND_DISP, ROUND_DISP_POSFILE. As explained above STRAIGHT means that the leaf profile is a rectangle and it moves along a circle, while ROUND means that the leaf profile tip closest to the beam has circular shape and it moves along a straight line. ROUND_DISP means that the centre of the circle that defines the profile is not at the line in the middle of the leaf. And ROUND_DISP_POSFILE ,means that the tip positions are a function of the field aperture and are read from a file, instead of calculating them automatically.

- Leaf half dimension in the inline plane (LONGDIM). It is X half dimension if the orientation is X and Y half dimension if the orientation is Y. If the endleaf is curved this length is defined as the maximum extent. If the endleaf is straight, and therefore it has a circular shape, the length is defined as the arc length of the upper part of the leaf (the side closest to the source).

- If the leaf tip type is ROUND, ROUND_DISP or ROUND_DISP_POSFILE, leaf tip radius (TIP_RADIUS)

- If the leaf tip type is ROUND_DISP or ROUND_DISP_POSFILE, distance between leaf top (part closest to the beam origin) and radius centre (DIST). If the leaf tip type is ROUND, the radius center is at the middle of the leaf, so this parameter is not needed.

- If the leaf tip type is ROUND or ROUND_DISP, half value layer (HVL). This half value layer should correspond to the distance at which only half of the gammas with the spectrum generated by the accelerator traverse the material of the leaf. It is used to move the leafs so that the line from the focus to the aperture position in the isocenter traverses a distance of the rounded leaf tip equal to the half value layer.

- If the leaf tip type is ROUND_DISP_POSFILE, name of file with leaf aperture (the number that defines the field) vs. leaf positions (the intersection of the leaf tip at the center of curvature with the isocentre plane. In figure RT1 these two magnitudes correspond to the light field edge, XL(i), and leaf tip position, W(i).

- Z position of the focus

- Cross position of the focus. It the leaves orientation is X the cross position is in the Y axis and viceversa.

- Z of position of the isocenter. The XY plane at this Z is used to calculate the cross positions that define the leaf profiles as well as the positions that the define the leaf apertures. Following the IEC 61217 convention, the beam is assumed to run along the negative Z axis, therefore this position should be a value smaller than the Z position of the focus. F

- Z of position of the top of the leaves (the position closer to the beam source) (Z_TOP). This position will be used as the 0 reference to place the leaves. Indeed, it is used as the Z=0 position for defining the leaf profiles; therefore the Z positions that define the leaf profile should be negative (further from the beam source).

- Distance in Z between Z_TOP and the plane where leaf gap is defined (Z_GAP). Remember that if you define it at the bottom of the leaves, it should have a negative value.

- Gap between leafs in the cross direction (CROSS_GAP). Separation between leaves defined at Z_GAP.

- Cross leaf starting point. Cross position of the leaf situated at the most negative position. It is defined at the same Z plane as the leaf gap. (CROSS_START)

- Number of different types of leaf cross profiles.

For each leaf cross profile type the following data has to be defined:

  • - Leaf type. It will be a LEAF or BLOCK, but this options is not implemented yet, so it must be LEAF
  • - Number of 2-dimensional points that define the cross profile
  • For each point the Z and cross coordinates must be defined - Z coordinate - Cross coordinate

Number of leaf pairs

For each leaf pair the following data has to be defined:

  • Leaf type number 1, 2, 3, …; following the order of leaf types defined
  • - Cross projection of the left part (negative coordinate) of the leaf on isocenter plane The projection is taken as that of the leaf profile point that gives an smaller (more negative) coordinate.
  • - Opening distance of negative leaf, projected on the isocenter plane In case of rounded leaf tips, the point projected is not the end of the leaf, but a line is traced from the source towards the isocenter plane traversing the leaf, so that the distance traversed is the same as the half value layer given above.
  • - Opening distance of positive leaf, projected on the isocenter plane

Finally the material of the leaves have to be defined

material

Mother volume name

../_images/MLC.along.jpg

Fig. RT2. YZ plane view of a cut along X=30 mm of the MLC created from the file below

../_images/MLC.cross.jpg

Fig. RT3. XZ plane view of a cut along Y=9 mm of the MLC created from the file below

../_images/MLC.leafs.jpg

Fig. RT4. Leafs geometry of the MLC created from the file below

Module MLC example

On the next lines we include a simplified geometrical description, which can be used as example. On figures RT2, RT3 and RT4 it can we seen the corresponding representation.

:P MLC_LH 2.

:P MLC_LVD 2.

:P MLC_LHD 0.2

:P MLC_LX 2.

:P MLC_LVN1 20.

:P MLC_LVN2 50.

:P MLC_LVN3 60.

:P MLC_LVN4 100.

:P MLC_LVP1 10.

:P MLC_LVP2 $MLC_LVN1

:P MLC_LVP3 $MLC_LVN2

:P MLC_LVP4 75.

:P MLC_LVP5 $MLC_LVN4

:P MLC_TIP_CIRCLE_R 150

:P MLC_TIP_CIRCLE_Z 40.

:P MLC_ZPOS_TOP -100

:P MLC_CROSS_FOCUS 5.

:P MLC_IsocentrePositionZ -300.

:P MLC_GAP 0.2

:P MLC_NLEAF_PAIRS 7

:P MLC_HVL 70.

:P MLC_FIELD 30.

:MODULE MLC

MLC_X Module_name

FOCUSED X ROUND_DISP_POSFILE // Type, Orientation, LeafTip type

100./2. // Leaf long half-length

$MLC_TIP_CIRCLE_R $MLC_TIP_CIRCLE_Z // Tip_circle_radius Tip_circle_centre_Z

$MLC_HVL // Half value layer

  1. $MLC_CROSS_FOCUS $MLC_IsocentrePositionZ // Z focus, CrossLeaf focus, Z isocentre

$MLC_ZPOS_TOP // Z top

-$MLC_LVN4 $MLC_GAP -($MLC_LH*($MLC_NLEAF_PAIRS-2)+($MLC_LH*4)*2+$MLC_GAP*($MLC_NLEAF_PAIRS-1))/2. // Z gap, Cross Leaf Gap, Cross Leaf Start Point

3 // N Types of Leaf Cross Profiles

LEAF 27 // Leaf type 1: Leaf Mode, N Leaves Cross Points

  1. -$MLC_LH/2.
  1. -$MLC_LH/6.

$MLC_LX -$MLC_LH/6.

$MLC_LX $MLC_LH/6.

  1. $MLC_LH/6.
  1. $MLC_LH/2.

-$MLC_LVP1 $MLC_LH/2.

-($MLC_LVP1+$MLC_LVD/2.) $MLC_LH/2.-$MLC_LHD

-($MLC_LVP1+$MLC_LVD) $MLC_LH/2.

-$MLC_LVP2 $MLC_LH/2.

-($MLC_LVP2+$MLC_LVD) $MLC_LH/2.-$MLC_LHD

-$MLC_LVP3 $MLC_LH/2.-$MLC_LHD

-($MLC_LVP3+$MLC_LVD) $MLC_LH/2.

-$MLC_LVP4 $MLC_LH/2.

-($MLC_LVP4+$MLC_LVD/2.) $MLC_LH/2.-$MLC_LHD

-($MLC_LVP4+$MLC_LVD) $MLC_LH/2.

-$MLC_LVP5 $MLC_LH/2.

-$MLC_LVP5 0.

-($MLC_LVP5+$MLC_LX) 0.

-($MLC_LVN4+$MLC_LX) -$MLC_LH/2.

-($MLC_LVN3+$MLC_LVD) -$MLC_LH/2.

-($MLC_LVN3+$MLC_LVD/2.) -$MLC_LH/2.+$MLC_LHD

-$MLC_LVN3 -$MLC_LH/2.

-($MLC_LVN2+$MLC_LVD) -$MLC_LH/2.

-$MLC_LVN2 -$MLC_LH/2.-$MLC_LHD

-($MLC_LVN1+$MLC_LVD) -$MLC_LH/2.-$MLC_LHD

-$MLC_LVN1 -$MLC_LH/2.

LEAF 12 // Leaf type 2: Leaf Mode, N Leaves Cross Points

$MLC_LX -$MLC_LH*2.

$MLC_LX 0.

  1. $MLC_LH*2.

-$MLC_LVN1 $MLC_LH*2.

-($MLC_LVN1+$MLC_LVD) $MLC_LH*2-$MLC_LHD

-$MLC_LVN2 $MLC_LH*2-$MLC_LHD

-($MLC_LVN2+$MLC_LVD) $MLC_LH*2.

-$MLC_LVN4 $MLC_LH*2.

-$MLC_LVN4 0.

-($MLC_LVN4+$MLC_LX) 0.

-($MLC_LVN4+$MLC_LX) -$MLC_LH*2.

LEAF 12 // Leaf type 3: Leaf Mode, N Leaves Cross Points

  1. -$MLC_LH*2.

$MLC_LX 0.

$MLC_LX $MLC_LH*2.

-($MLC_LVN4+$MLC_LX) $MLC_LH*2.

-($MLC_LVN4+$MLC_LX) 0.

-$MLC_LVN4 0.

-$MLC_LVN4 -$MLC_LH*2.

-($MLC_LVN2+$MLC_LVD) -$MLC_LH*2.

-$MLC_LVN2 -$MLC_LH*2-$MLC_LHD

-($MLC_LVN1+$MLC_LVD) -$MLC_LH*2-$MLC_LHD

-$MLC_LVN1 -$MLC_LH*2.

$MLC_NLEAF_PAIRS // Number of leaves

2 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

1 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

1 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

1 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

1 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

1 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

3 -$MLC_FIELD $MLC_FIELD // Leaf type, negative and positive aperture at isocentre

RTUW ACCEL // Material mother_volume

Range modulator module

EJEMPLO:

//:MODULE RANGE_MODULATOR NAME RMIN RIN ROUT NBLADE STEPS n{ thigness weight } :MODULE RANGE_MODULATOR range_modulators_MODULE 85/2 85*2/2 300 4 4 6.5 0.1811111 11 0.12433333 17.1 0.09644444 22.3 0.0953889

EJEMPLO:

//:MODULE RANGE_MODULATOR NAME RMIN RIN ROUT NBLADE STEPS n{ thigness weight } :MODULE RANGE_MODULATOR range_modulators_MODULE 85/2 85*2/2 300 4 4 6.5 0.1811111 11 0.12433333 17.1 0.09644444 22.3 0.0953889

RMIN es el radio del agujero del soporte central

Rin es el radio a partir del cual comienza el modulador escalonado

Rout es el radio donde termina el modulador escalonado

NBLADE indica el n?mero de brazos del modulador

STEPS indica el numero de espesores que tiene el modulador

Luego se da la lista de espesores y el porcentaje en peso que representa ese espesor en el modulador, la suma de porcentajes es menor de 1 ya que la diferencia entre la suma de los pesos y 1 es la parte que queda libre.

Using phase spaces

A very common utility in teletherapy simulation is the writing of a phase space, i.e. the set of particles that reach a certain plane in Z, and the later starting of the next simulation from this set of particles. This technique can save a lot of time in cases several calculations share some accelerator parts and in the case where the accelerator simulation is very slow compared to the dose calculation.

Writing phase spaces

The writing of phase space can be done automatically in GAMOS by selecting the user action

/gamos/userAction RTPhaseSpaceUA

When a particle crosses any of the planes defined by the parameter

/gamos/setParam RTPhaseSpaceUA:ZStops Z_1 Z_2 Z_3 …

its information is stored in a file whose named is given by the parameter

/gamos/setParam RTPhaseSpaceUA:FileName MY_FILENAME

plus a suffix .IAEAphsp. The default value of this parameter is test. If there are several Z planes, the Z value of each plane will be added to the name, with a ``_’’ in front, so that each phase space will go to a different file. If there is only one Z defined, the Z value may be written in the file name if the following parameter is set

/gamos/setParam RTPhaseSpaceUA:ZStopInFileName TRUE

If the plane crossed is the one with maximum Z, the particle may be stopped, if the following parameter is set

/gamos/setParam RTPhaseSpaceUA:KillAfterLastZStop TRUE

The format of the phase space files is the one defined by the IAEA, generated using the official C files from IAEA. First there is a header file, that will have the same name, but with the suffix .IAEAheader. It is generated by the IAEA code and you can find the description of it at [biblio.iaeaphsp].

The variables stored for each particle are the following

  • X coordinate (cm)
  • Y coordinate (cm)
  • X direction cosine
  • Y direction cosine
  • Kinetic energy (MeV)
  • Particle statistical weight
  • Type of the particle (1 = gamma, 2 = electron, 3 = positron, 4 = neutron, 5 = proton)
  • Z direction cosine * particle charge
  • Is new history
  • Extra integers (0, 1 or 2 values)
  • Extra floats (0, 1 or 2 values)

The Z value may optionally be stored if the following parameter is set

/gamos/setParam RTPhaseSpaceUA:StoreZ TRUE

The header file may be written each N events, so that if the job ends abnormally the first N events may be recovered in the phase space. The number of events elapsed between header writing is controlled with the parameter

/gamos/setParam RTPhaseSpaceUA:NEventsToSave TRUE/FALSE

By default this parameter takes a value of -1 and no header is saved until the end of the run.

It is also possible to limit the number of particles in a file. This can be useful if you want to store a phase space at different Z values and you do not want too many particles at the smallest Z planes. To do this you can set the parameter:

/gamos/setParam RTPhaseSpaceUA:MaxNTracksInFile NTRACKS

The number of tracks stored in each file will be stopped at NTRACKS, but it may continue in other files.

As the Z plane is probably not a physical plane in the geometry, particle steps will start before the plane and end after the plane. Therefore, the position and energy are re-scaled as if the particle would have stopped at the Z plane (by a simple linear interpolation).

By default particles that travel towards the negative Z direction are not stored, as it is assumed that the original direction was towards positive Z and therefore this particles (or their ancestors) should already have been stored. You may nevertheless eliminate this check and store all particles by setting the parameter:

/gamos/setParam RTPhaseSpaceUA:NotStoreBackwards FALSE

You can limit the number of particles stored in each phase space (what may be useful when you are filling several phase spaces in a job and do not want that the first phase spaces end occupying too much disk space. To do this you have to set the parameter:

/gamos/setParam RTPhaseSpaceUA:MaxNTracksRecorded NTRACKS

Phase space text file

You can also write the phase space in a simple text format. The following variables will be written for each track that reaches the phase space plane: particle_code kinetic_energy(MeV) x(cm) y(cm) z(cm) direction_x direction_y direction_z weight extra_information_word_1 extra_information_word_2 … extra_information_word_n. The name of the file is set with the parameter:

/gamos/setParam RTPhaseSpaceUA:TextFileName FILE_NAME

Phase space histograms

When a phase space is generated, a set of histograms are produced to give you some information about the particles in the phase space. The name of all these histograms start with PhaseSpace: and they are all dumped into the file phaseSpace.root/csv. For each of the following particle types a different set of histograms are produced (containing the particle type in the histogram name:) “gamma”, “e-“, “e+”, plus a set of histograms for all particles, named “ALL”. Also histograms for “neutron” and “proton”, will be produce if the following parameter is set

/gamos/setParam RTPhaseSpaceHistos:Hadrons TRUE

A full set of histograms is produced for each Z plane, with the Z value on its name.

The following histograms are produced:

  • Position X (cm) (“X at Zstop”)
  • Position Y (cm) (“Y at Zstop”)
  • Position X and Y (cm) (“XY at Zstop”)
  • Radial position in the XY plane (cm) (“R at Zstop”)
  • Theta angle of direction (degrees) (“Direction Theta at Zstop”)
  • Phi angle of direction (degrees) “Direction Phi at Zstop”)
  • Energy (MeV) (“Energy at Zstop”)
  • X direction cosine (“Vx at Zstop”)
  • Y direction cosine (“Vy at Zstop”)
  • Z direction cosine (“Vz at Zstop”)
  • Radial position in the XY plane (cm) vs theta angle of direction (degrees) (“R vs Direction Theta at Zstop”);
  • Radial position in the XY plane (cm) vs energy (MeV) (“R vs Energy at Zstop”);
  • Theta angle of direction (degrees) vs energy (MeV) (“Direction Theta vs Energy at Zstop”);

The user can choose not to produce any of these histograms by setting the following parameter

/gamos/setParam RTPhaseSpaceUA:Histos FALSE

The name of the histogram file can be controlled with the parameter

/gamos/setParam RTPhaseSpaceHistos:FileName MY_FILENAME

that by default is “phaseSpace”.

use the commands in /gamos/analysis/ (see section on Histogramming) to control the minimum, maximum and number of bins of the histograms

Reading phase spaces

To use the generated phase space you have to define as your primary generator

/gamos/generator RTGeneratorPhaseSpace

You can change the file name (by default test) with the parameter

/gamos/setParam RTGeneratorPhaseSpace:FileName MY_FILENAME

The particles in the phase space can be translated or rotated by using the parameters

/gamos/setParam RTGeneratorPhaseSpace:InitialDisplacement POS_X POS_Y POS_Z

/gamos/setParam RTGeneratorPhaseSpace:InitialRotAngles ANG_X ANG_Y ANG_Z

the position of the particles is first changed by the initial displacement, then the momentum vector is rotated around the X axis by ANG_X, around the Y axis by ANG_Y and around the Z axis by ANG_Z; and finally the position vector is rotated around the X axis by ANG_X, around the Y axis by ANG_Y and around the Z axis by ANG_Z.

An alternative way to provide the initial rotation and displacement is by setting the transformation. The transformations are set with the parameter:

/gamos/setParam RTGeneratorPhaseSpace:Transformations TYPE_1 VAL1_1 VAL2_1 VAL3_1 TYPE_2 VAL1_2 VAL2_2 VAL3_2

Several transformations can be set at the same type with this parameter, and they will be executed in the order provided. For each transformation four parameters have to be supplied: transformation type, first value, second value, third value.

Three types of transformation are supported:

  • Displacement (D): the three values are the displacement in X, Y and Z
  • Rotations around XYZ (RXYZ): the three values are the angles of rotation around the X, Y and Z axis (always executed in this order)
  • Rotation around accelerator axis (RTPS): the three values correspond to the angles of rotation in theta, phi and around the accelerator axis itself (in other words, a rotation around the Z axis is done with angle VAL3, then a rotation around Y with angle VAL1 and finally a rotation around Z with angle VAL2

If the number of phase space particles is not enough to calculate the dose in the phantom with enough precision you may reuse the particles several times by setting the parameter

/gamos/setParam RTGeneratorPhaseSpace:MaxNReuse N

to a value bigger than 1. In fact, if this parameter is not explicitly set to 1, GAMOS calculates it automatically by dividing the number of events asked for by the number of events in the input phase space.

Alternatively you may recycle the full phase space several times (i.e. when all particles in the phase space are read, the file is closed and restarted again). The number of times a phase space is recycled is controlled by the parameter

/gamos/setParam RTGeneratorPhaseSpace:MaxNRecycle N

which must take a value bigger than 1. If no reusing is explicitly set, the number of reuses will be automatically calculated as mentioned above and the number of recyclings will be set to 1, so that no recycling is done.

Caution should be taken when recycling phase spaces, as the error correlations due to the reusing of the same particles is not taken account. Therefore you should first consider reusing instead of recycling.

If you think your phase space has X/Y symmetry you may reduce the correlations due to reusing the same particle several times by mirroring your particles each time they are reused, by setting the parameter

/gamos/setParam RTGeneratorPhaseSpace:MirrorWhenReuse OPT

where OPT can be X or Y or XY. If it is X it means that when a particle is used an even time the X original value of position and momentum direction will be changed sign, while when it is an uneven time it will not be changed. If it is Y the same but mirroring in Y. And if it XY the second time a particle is used the X values will be changed sign, the third time the Y values will be changed sign, the fourth time both X and Y values will be changed sign, the fifth time it will be no change, and the cycles restarts.

When reading a phase space file you may skip the first N events by using the parameter

/gamos/setParam RTGeneratorPhaseSpace:NEventsSkip N

You can produce the same histograms that were produced when writing the phase space file by setting the following parameter

/gamos/setParam RTGeneratorPhaseSpace:Histos TRUE

It is also possible to read a phase space in EGS format, by using as your primary generator

/gamos/generator RTGeneratorPhaseSpace_EGS

Adding extra information to a phase space

The IAEA format allows to store two long integers and two floats in the phase space file as extra information. In GAMOS we have extended this functionality by putting the 32+32 bits of the two integers in a continuous stream, that the user can divide in groups of bits of the desired length to store several data. The user can even add more sets of 32 bits by changing the the following line at source/RadioTherapy/include/iaeaRecord.hh:

#define NUM_EXTRA_LONG  2

**

and recompiling (cd source/RadioTherapy; make).

To store some information in this format the user has to instantiate one of the following user actions

  • RTExtraInfoProviderCrossings: fills a bit if the track has crossed the corresponding region before reaching the phase space
  • RTExtraInfoProviderInteractions: fills a bit if the track has interacted in the corresponding region before reaching the phase space
  • RTExtraInfoProviderOrigin: fills a bit if the track has been created in the corresponding region before reaching the phase space

The user may select how many bits each information must occupy by using the GAMOS parameter :

/gamos/setParam EXTRA_INFO_NAME:NBits NBITS

where EXTRA_INFO_NAME is the name of the extra information class (see above). GAMOS will check that the index to be filled by a class is not bigger than the number of bits reserved for it. And it will also check that the total number of bits is not bigger than the available quantity (32*NUM_EXTRA_LONG).

The order of declaration of the extra information user actions sets the order of bits occupancy. At the end of the job each of these extra information user actions fills a file explaining the information contained in each bit. By default this file is called RTExtraInfoProvider.summ, but the user may change the name of it with the parameter

/gamos/setParam RTExtraInfoProviderLong:FileName FILE_NAME

An example of a file is the following one:

RTExtraInfoProviderOrigin INDEX/REGION 0 DefaultRegionForTheWorld
RTExtraInfoProviderOrigin INDEX/REGION 1 targetReg
RTExtraInfoProviderOrigin INDEX/REGION 2 collimatorReg
RTExtraInfoProviderOrigin INDEX/REGION 3 filterReg
RTExtraInfoProviderOrigin INDEX/REGION 4 monitorReg
RTExtraInfoProviderOrigin INDEX/REGION 5 shieldingReg
RTExtraInfoProviderOrigin INDEX/REGION 6 jawsReg

The float extra information providers fill one of the two available words in the order they have been defined.

The user can add more float words by changing the the following line at source/RadioTherapy/include/iaeaRecord.hh:

#define NUM_EXTRA_FLOAT  2

and recompile the GAMOS code after that.

To store some information as float the user has to instantiate one of the following user actions

  • RTExtraInfoProviderZOrigin: stores the Z position of the origin of the track
  • RTExtraInfoProviderZLast: stores the Z position of the last interaction before reaching the phase space
  • RTExtraInfoProviderZLastNoMaterial: stores the Z position of the last interaction before reaching the phase space, but you may define a list of materials that will not be taken into account. by default this list is G4_AIR, but it can be changed with the parameter /gamos/setParam RTExtraInfoProviderZLastNoMaterial:Materials MATERIAL_1 MATERIAL_2 …

At the end of the job each of these extra information user actions fills a file explaining the information contained in each float. By default this file is called RTExtraInfoProvider.summ, but the user may change the name of it with the parameter

/gamos/setParam RTExtraInfoProviderFloat:FileName FILE_NAME

If you have written a phase space file in a 32-bit machine and try to read it in a 64-bit one, you will not be able to recover the numbers properly. You may avoid it with the parameter

/gamos/setParam IAEArecord:ExtraLongSize 32

Similarly if you wrote the file in a 64-bit machine and try to read it in a 32-bit one.

Filtering and classifying on extra information

It is possible to define a filter and a classifier (see sections on Filters and Classifiers) based on the long extra information.

The filter name is RTEILongFilter. When it is defined it must receive as extra argument the list of extra info values that are accepted:

/gamos/filter MYFILTER RTEILongFilter VALUE_1 VALUE_2 …

The classifier is named RTClassifierByEILong. It will read the names of the indices of the extra long information from a file named RTExtraInfoProvider.summ. The name of this file can be changed with the parameter:

/gamos/setParam CLASSIFIER_NAME:FileName VALUE

Reusing a particle at a phase space without writing the phase space file

It is common in radiotherapy treatment simulations that in a first job a phase space file is created and in a second job the particles therein are read and reused several times to calculate the dose in the patient. GAMOS provides the possibility of reusing particles in a phase space file without having to divide the jobs in two. To use this utility the user has to instantiate the user action x /gamos/userAction RTReuseAtZPlaneUA

and then use the command

/gamos/RT/ReuseAtZPlane

and before the user has to define the z value of the plane where particles will be replicated, with the parameter

/gamos/setParam RTReuseAtZPlane:ZReusePlane Z_POS

and the number of times particles will be reused (indeed the particles is copied NREUSE-1 times)

/gamos/setParam RTReuseAtZPlane:NReuse NREUSE

The user may independently create the phase space file or not at the same Z position or at others.

An important thing to take into account is that the particle is copied from the position of the PreStepPoint. A possible option would be to extrapolate it to the Z plane and copied it at that position, but if the extrapolation is done with a straight line, it would not take into account the straggling of charged particles.

If you think your phase space has X/Y symmetry you may reduce the correlations due to reusing the same particle several times by mirroring your particles each time they are reused, by setting the parameter

/gamos/setParam RTReuseAtZPlane:MirrorWhenReuse OPT

where OPT can be X or Y or XY. If it is X it means that when a particle is used an even time the X original value of position and momentum direction will be changed sign, while when it is an uneven time it will not be changed. If it is Y the same but mirroring in Y. And if it XY the second time a particle is used the X values will be changed sign, the third time the Y values will be changed sign, the fourth time both X and Y values will be changed sign, the fifth time it will be no change, and the cycles restarts.

Optimisation of a linac simulation

To optimise your simulation in Geant4 you may tune the physics parameters (production cuts, special cuts, multiple scattering options, …) as well as try some reduction variance techniques.

Cuts optimisation

Please read the sections on automatic optimisation of production cuts and user limits for accelerator and dose calculation simulations.

Electromagnetic parameters optimisation

There are many parameters that control the speed of the electromagnetic physics vs its precision. Please see the web page http://fismed.ciemat.es/GAMOS/RToptim to get a list of the best electromagnetic physics parameters that can optimise your simulation.

Particle splitting

Particle splitting is a non-biased variance reduction technique that may reduce the CPU time of your accelerator simulation by a big factor. It basically consists on splitting the secondary particles (in radiotherapy mainly gammas generated by bremsstrahlung) N times, so that each time a bremsstrahlung gamma is created, another N-1 gammas are created at the same position, with weight 1/N, re-sampling the energy and/or angle distribution. As most of the particles that reach a patient in a radiotherapy accelerator are bremsstrahlung gammas, by using this technique you can spare the time spent simulating the original electron and also the time spent tracking gammas that have small possibility of reaching the patient.

There are three splitting techniques implemented in GAMOS:

  • Uniform bremsstrahlung splitting
  • Directional bremsstrahlung splitting
  • Equal weight particle splitting

The three techniques are explained in the sections on Biasing operations in the chapter on Variance reduction techniques.

Killing particles at big X/Y

This utility serves to kill the particles that would probably not reach your detector because they have too big X and Y positions (it is assumed that that your detector is described along Z and that your initial particles move in this direction in the positive sense).

It makes a list of the volumes that are placed directly on the world volume and computes the minimum and maximum extension in X and Y (using the method G4VSolid::GetExtent()). This rectangular area extends from the minimum Z value of this volume until the minimum Z value of the next volume, so that all tracks that are outside it will be killed. This utility is activated by selecting the user action

/gamos/userAction RTZRLimitsAutoUA

You may argue that too many particles (or too few) are killed with this automatic definition of limits. You may simply changed the dimensions of the volumes placed at world (adding container volumes made of air will not change your physics). or you can define the limits your self by selecting the user action

/gamos/userAction RTZRLimitsUA

Then you have to write a file named xylimits.dat with the list of values you want to use. The format of that file should be the following: a set of lines with three numbers representing the Z value of the plane and then the X and Y limits. The planes will be extended in Z until the previously defined Z (for the minimum Z defined the world negative Z limit will be used).

If you want to change the name of the limits file, you can do it with the command (remember to define a parameter always before selecting the user action)

/gamos/setParam RTZRLimitsUA:FileName MY_FILENAME

The user has also the option that particles rejected are not killed, but Russian Roulette is played with them and only if they are rejected they are killed. If they are accepted their weight will be increased by the corresponding value. To use this option the following parameter must be set:

/gamos/setParam RTVZRLimitsUA:UseRussianRoulette TRUE

The value of the Russian roulette threshold is set with the parameter

/gamos/setParam RTVZRLimitsUA:RussianRouletteValue VALUE

which by default takes a value of 100.

Scoring dose in phantom

To use a phantom geometry you can use the GAMOS utilities to read phantom geometries in EGS or Geant4 formats or build simple phantoms with a few user commands. The scoring of the dose in the phantom volumes can be done using the scorer GmG4PSDoseDeposit and selecting as detector the voxels, that are named phantom. For example you can use the following commands:

/gamos/scoring/createMFDetector DoseDet patient

/gamos/scoring/addScorer2MFD DoseScorer GmG4PSDoseDeposit DoseDet

This is all you need to get on the standard output the dose by event deposited in each voxel. When you use a phase space as input the number of events is not the number of events run in your job (what would equal be the number of particles in the phase space multiplied by the number of times each particles is reused), but GAMOS uses as number of events the original number of events that were used to generate the phase space. In fact, GAMOS gets the ratio of particles written in the phase space per original event transported through the accelerator and multiplies this ratio by the number of phase space particles used. This allows to get the best approximation to the correct number of events when the phase space is not used fully, or when particles are reused (it is indeed not the exact number if for example you use only the 10 first particles in the phase space: the number of events that generated those 10 particles may be bigger or smaller than the number of events that generated the following 10 particles, due to statistical fluctuations).

If you are using the option of skipping voxel frontiers when the material does not change (see section on Reading DICOM files) you may take into account the possibility of tuning the distributions of dose deposited in the voxels traversed by one step (see section on Scoring in voxelised phantoms).

Saving scores and scores squared in binary file

The second available file is a binary file where it is written the dose and dose squared in each voxel. The binary format allows for a faster writing and reading and the fact of writing the dose squared instead of the error serves to calculate in a proper way the error correlations when the doses from several jobs are added. To obtain it you just have to add the scorer printer GmPSPrinterSqdose to your scorer, for example

/gamos/scoring/addPrinter2Scorer GmPSPrinterSqdose DoseScorer

The output file is named by default sqdose.out. You may change it with the parameter

/gamos/setParam PRINTERNAME_SCORERNAME:FileName MY_FILENAME

where*PRINTERNAME* is the name of the printer, by default GmPSPrinterSqdose, and SCORERNAME is the name you have given to the scorer.

All variables are of type float and the format is the following:

  • Number of events (original number of events used to generate a phase space file is phase space file is used as primary generator)
  • Number of voxels in the X, Y and Z directions (e.g. n_x, n_y, n_z)
  • Array of voxel boundaries (cm) in the X direction (n_x+1 values)
  • Array of voxel boundaries (cm) in the Y direction (n_y+1 values)
  • Array of voxel boundaries (cm) in the Z direction (n_z+1 values)
  • Array of dose values (Gy/cm^3) (n_x n_y n_z values)
  • Array of dose squared values ((Gy/cm^3)^2) (n_x n_y n_z values)

To optimise the disk space and the memory needed to read the file back, it is possible to only store the dose and dose squares of those voxels that are filled. This can be done by setting the parameter:

/gamos/setParam GmSqdose:FileType FILLED

The default value of this parameter is ALL.

Saving scores and score errors in text file

You can also store the dose in each voxel in a file, what allows to calculate the average dose calculated with several jobs (see dose analysis section). Nevertheless we recommend you to use the previous format instead of this one, as it is faster for writing and reading.

The first file is a text file where it is written the dose and dose error in each voxel. To obtain it you just have to add the scorer printer GmPSPrinter3ddose to your scorer, for example

/gamos/scoring/addPrinter2Scorer GmPSPrinter3ddose DoseScorer

The output file is named by default 3ddose.out. You may change it with the parameter

/gamos/setParam PRINTERNAME_SCORERNAME:FileName MY_FILENAME

where*PRINTERNAME* is the name of the printer, by default GmPSPrinter3ddose, and SCORERNAME is the name you have given to the scorer.

The format is the same as the one 3ddose format used in DOSXYZnrc, except that the first line contains the number of events:

  • Number of voxels in the X, Y and Z directions (e.g. n_x, n_y, n_z)
  • Array of voxel boundaries (cm) in the X direction (n_x+1 values)
  • Array of voxel boundaries (cm) in the Y direction (n_y+1 values)
  • Array of voxel boundaries (cm) in the Z direction (n_z+1 values)
  • Array of dose values (Gy/cm^3) (n_x n_y n_z values)
  • Array of dose relative error values (n_x n_y n_z values)

Saving scores in histograms

If you want to store the scores in the phantom in an histogram file, you can use the scorer printer RTPSPDoseHistos

/gamos/scoring/addPrinter2Scorer RTPSPDoseHistos DoseScorer

The name of the histogram file if by default dose_PRINTER_NAME, but it may be chanted with the parameter

/gamos/setParam PRINTER_NAME:FileName VALUE

The number of bins in the histograms will be equal to the number of voxels in the phantom (nVoxel) and the histogram limits are the -nVoxel/2 +to nVoxel/2.

The name of all these histograms start with RTPSPDoseHistos: and they are all dumped into the file dose.root/csv.

The following histograms are produced by default:

  • Dose in X direction (“Dose Profile X_merged”)
  • Dose in Y direction (“Dose Profile Y_merged”)
  • Dose in Z direction (“Dose Profile Z_merged”)
  • Dose in XY direction (“Dose XY_merged”)
  • Dose in XZ direction (“Dose XZ_merged”)
  • Dose in YZ direction (“Dose YZ_merged”)
  • Dose of each voxel (“Dose”)
  • Integrated dose of each voxel, i.e. all the voxels that have dose equal or greater that a given bin fill that bin (“Dose-volume”)

The word merged refers to the fact that this histograms count the dose in one direction (or two) integrating the dose in all the voxels in the other two (or one) direction.

If you want to obtain other histograms out of the dose in the voxels, you may list what you want in a file in the following way: first give the name of your file with the parameter:

/gamos/setParam RTPSPDoseHistos:HistosFileName FILE_NAME

In this file you have to fill a line for each histogram you want with the following information:

  • Histogram type: it must be 1X, 1Y, 1Z, 2XY, 2XZ or 2YZ
  • Histogram name
  • Minimum voxel in X
  • Maximum voxel in X
  • Minimum voxel in Y
  • Maximum voxel in Y
  • Minimum voxel in Z
  • Maximum voxel in Z

For the minimum and maximum voxel value you may use the total number of voxels by writing one of the words NVoxelX, NVoxelY or NVoxelZ. Also arithmetic expressions are comment (starting the line with // are allowed. One example file can be the following (for a phantom with number of voxels 101 101 100):

// recovering the 1D merged histograms
1X "RTPSPDoseHistos:Dose Profile X_merged" 0 NVoxelX-1 0 NVoxelY-1 0 NVoxelZ-1
1Y "RTPSPDoseHistos:Dose Profile Y_merged" 0 NVoxelX-1 0 NVoxelY-1 0 NVoxelZ-1
1Z "RTPSPDoseHistos:Dose Profile Z_merged" 0 NVoxelX-1 0 NVoxelY-1 0 NVoxelZ-1
// Y profiles at several depths
1Y "RTPSPDoseHistos:Dose Profile Y 1.400 cm" 50 50 0 NVoxelY-1 6 6
1Y "RTPSPDoseHistos:Dose Profile Y 10.000 cm" 50 50 0 NVoxelY-1 49 49
// PDDs merging voxels
1Z "RTPSPDoseHistos:Dose Profile Z 1x1" 50 50 50 50 0 NVoxelZ-1
1Z "RTPSPDoseHistos:Dose Profile Z 3x3" 49 51 49 51 0 NVoxelZ-1

It is also possible to produce in a simple way a full set of 2-dimensional plots covering a whole phantom in any of the three dimensions. To do one must use the following parameter:

/gamos/setParam RTPSPDoseHistos:AllHistos2D 2DIM_2 2DIM_2 ..

The 2DIM_i values can be XY, XZ or YZ (you can use one of several of these three in your command). To visualise them in a simple way you may use the printAll.C (see Appendix) ROOT utility to convert all the histograms in graphics files.

Analysis utilities

This is a set of utilities that may serve in the analysis of phase space and dose files, for example to sum phase space or dose files from different jobs, to get basic information of the file contents, to fill histograms out of the files or to compare files from two jobs.

All these utilities are under the directory analysis. They are compiled by default with the rest of GAMOS code and after that they are available as executables as mentioned in the following sub-sections.

Summing phase space files

This utility serves to sum phase space files corresponding to different jobs with the same setup. To use it you have to write a file containing the list of header phase space files, one file per line, for example

ps.20000.IAEAheader
ps.20001.IAEAheader
ps.20002.IAEAheader

Then you just have to run the executable

sumPS INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to add and OUTPUT_FILE_NAME is the name of the output file that will contain the sum of all the files (two files indeed as usual for IAEA phase space files: OUTPUT_FILE_NAME.IAEAheader and OUTPUT_FILE_NAME.IAEAphsp ).

When running you will see on the screen something similar to this:

.. COMMENT: \\tiny
Opening phase space contained in     ps.20000.IAEAheader of type IAEA
PARTICLES 225437 NPART_TOT 225437 NPARTORIG_TOT 5000000
RATIO 0.0450874 +- 0.000101661 RATIO_TOT 0.0450874
Opening phase space contained in     ps.20001.IAEAheader of type IAEA
PARTICLES 224635 NPART_TOT 450072 NPARTORIG_TOT 10000000
RATIO 0.044927 +- 0.000101455 RATIO_TOT 0.0450072
Opening phase space contained in     ps.20002.IAEAheader of type IAEA
PARTICLES 224813 NPART_TOT 674885 NPARTORIG_TOT 15000000
RATIO 0.0449626 +- 0.000101501 RATIO_TOT 0.0449923
* * * * N Particles   = 674885
* * * * N Photons   = 673804
* * * * N Electrons = 1057
* * * * N Positrons = 24
* * * * N Original Histories =  1.5e+07

For each phase space files after the name of the file comes a line with the file statistics: number of particles, accumulated number of particles of all files, accumulated number of original histories, ratio particles/original histories +- ratio error, accumulated ratio particles/original histories. And at the end the statistics on the total number of particles in total, photons, electron and positrons, and the total number of original histories in all summed files.

Making histograms out of a phase space file

You can make histograms out of the particles contained in a phase space file by running gamos with the input script that you can file in RadioTherapy/analysis/phaseSpace/analysePS/rt.analysePS.in. If you edit it and change the name of the input phase space file, at /gamos/setParam RTGeneratorPhaseSpace:FileName, and run

gamos rt.analysePS.in

You will get a file named phaseSpace.root, that contains the same histograms that you get when you run the job to write the phase space file (see section on Phase space histograms).

Alternatively you can directly analyse a phase space file with the executable analysePS. You just have to type in your shell:

analysePS FILE_NAME

where FILE_NAME is the name of the phase space header file (the one that ends with .IAEAheader; you may indeed omit the suffix .IAEAheader if you want).

When running you will see on the screen something similar to this:

.. COMMENT: \\tiny
READING FILE ps.ubs.10x10.6
ORIGINAL HISTORIES= 2000000
N PARTICLES= 5.13303e+07 PER ORIG_HIST= 25.6651 +- 0.0184982
N GAMMAS= 5.11882e+07 PER ORIG_HIST= 25.5941 +- 0.0184479
N ELECTRONS= 137257 PER ORIG_HIST 0.0686285 +- 0.000191492
N POSITRONS= 4839 PER ORIG_HIST= 0.0024195 +- 3.48235e-05
Reading PHSP track 0
...
Reading PHSP track 51000000
=== saving histograms in file === phaseSpace.root

The output contains the total number of original histories and after the number of particles, gammas, electrons and positrons in the phase space file, and these numbers divided by the number of original histories, with the error. Also an histogram file equal to the file that was created when producing the phase space file will be created.

Several arguments can be supplied to the executable in the standard Unix format:

  • -f

    phase space file name (in this case do not use the file name alone as first argument as before)
  • -NRead

    number of particles to be read from the phase space file
  • -fOut

    output file name
  • -EMax

    maximum limit of the energy histograms
  • -RMax

    maximum limit of the position histograms
  • -NBins

    number of bins of histograms
  • -verb

    verbosity: it sets the RTVerbosity. Default is warning, that will print the above lines; debug, that will print each particle read form the phase space file
  • -help

    prints the set of arguments

Merging ‘sqdose’ files

This utility serves to merge dose files in the format sqdose corresponding to different jobs with the same setup, and getting the average dose of them. To use it you have to write a file containing the list of header phase space files, one file per line, for example

sqddose.water.20000.out
sqdose.water.20001.out
sqdose.water.20002.out

Then you just have to run the executable

mergeSqdose INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to merge and OUTPUT_FILE_NAME is the name of the output file that will contain the merging of all the files.

When merging files it will be checked that they correspond to the same phantom by checking the number of voxels and voxel limits.

Merging ‘3ddose’ files

This utility serves to merge dose files in the format 3ddose corresponding to different jobs with the same setup, and getting the average dose of them. To use it you have to write a file containing the list of header phase space files, one file per line, for example

3ddose.water.20000.out
3ddose.water.20001.out
3ddose.water.20002.out

Then you just have to run the executable

merge3ddose INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to merge and OUTPUT_FILE_NAME is the name of the output file that will contain the merging of all the files.

When merging files it will be checked that they correspond to the same phantom by checking the number of voxels and voxel limits.

Making histograms out of a ‘sqdose’ file

You can make histograms out of dose information contained in a sqdose file by running

analyseSqdose SQDOSE_FILE_NAME

When running you will see on the screen something similar to this:

READING FILE sqdose.ubs.10x10.6.out
GmSqdoseHeader::Read NEvent 2.5e+08
GmSqdoseHeader::Read NVoxels 21 112 150
GmSqdose::Read type 1
USING std::map to store doses
Number of voxels= 352800
RTPSPDoseHistos nvoxel 21 112 150
RTPSPDoseHistos dim 1 1 1
RTPSPDoseHistos translation (0,0,-4)
RTPSPDoseHistos rotation
[ (           1             0             0)
(           0             1             0)
(           0             0             1) ]
MINIMUM DOSE 2.352e-15
MAXIMUM DOSE 4.76001e-13
RTPSPDoseHistos AVERAGE ERROR 20% = 0.0139839
RTPSPDoseHistos AVERAGE ERROR 50% = 0.0111726
RTPSPDoseHistos AVERAGE ERROR 90% = 0.00996564
=== saving histograms in file === dose_analyseSqdose.root
N EVENTS IN SOURCE 2.5e+08

First it is printer the sqdose file header information: number of events (original events run) and number of voxels in X, Y and Z. Then the type of the sqdose (see above). Then the STL container that will be used to store the doses and dose errors (see argument descriptions below). Then the total number of voxels to be read and the information from the class RTPSPDoseHistos: number of voxels in X, Y and Z, voxel dimensions in X, Y and Z, translation and rotation applied to phantom when dose file was written and the statistics, i.e. minimum and maximum voxel dose, and average dose error in 20%, 50% and 90% highest dose voxels.

You will also get a file named dose_analyseSqdose.root, that contains the same histograms that you get when you run the job to write the dose file using the scorer printer RTPSPDoseHistos.

Several arguments can be supplied to the executable in the standard Unix format:

  • -f

    phase space file name (in this case do not use the file name alone as first argument as before)
  • -NRead

    number of particles to be read from the phase space file
  • -fOut

    output file name
  • -fHistos

    name of file with list of histograms
  • -DoseMin

    minimum limit of dose histograms
  • -DoseMax

    maximum limit of the dose histograms
  • -cont

    type of the STL container that will be used to store the doses and dose errors. It can be MAP or map, that uses a std::map; it is the default one, the one used by the Geant4 scorers.

    It can also be VECTOR or vector, that uses a std::vector; it occupies about ten times less than the std::map. The std::map container occupies a big amount of space, about 500 Mb for 10 million voxels, so we recommend you that you use std::vector if your phantom is big

  • -NVoxels

    total number of voxels in phantom (argument needed if sqdose file is of type FILLED
  • -verb

    verbosity: it sets the RTVerbosity. Default is warning, that will print the above lines; debug, that will print each particle read form the phase space file
  • -help

    prints the set of arguments

Automatic determination of production cuts for an accelerator simulation

The method used in GAMOS to determine the best production cuts is based on what we can call an ‘inverse reasoning’. We count each particle that reaches a given Z plane (corresponding to the phantom surface) and we calculate first the range of the particle in the region where it is created. Then we can know that if we put a range cut in that region smaller than the calculated range, that particle would not reach our target plane. We also compute the range of the mother particle in the region where it was created and the same consecutively for all the ancestors. We know then that if we set in any of the regions where each of the ancestor particles are created a cut smaller than the corresponding range, we would stop the chain of particles and therefore we would have no particle in the target plane. After running a good number of tracks we can know for each particle type and for each region which is the biggest range we can put if we do not want to lose any particle. Indeed we may allow to lose a small amount of particles if this speeds up our simulation. To know easily which is the biggest cut you can use to lose less than a given percentage of particles, GAMOS provides a set of plots (one per each particle type and per each region) and a simple script to get automatically the cut values.

One warning is due here: as mentioned above when a track reaches the target, its range fills a histogram, but also the range of all the ancestors of this track. It may happen then that when you set a certain cut and the above-mentioned script gives you how many tracks would be killed, more than one killed track correspond to the same track reaching the target (i.e., with a cut you kill the track that reaches the target and the parent track). Therefore you might have an overcounting of the number of tracks killed by a cut. To avoid this the total number of tracks (the last lines of output) is not computed as the sum of tracks in the region. This number uses a histogram that contains only one entry per track reaching the target, the one corresponding to the track with the smallest range. If you want to set a different cut for each region and are worried for this double counting, you may have a look at the histogram named “trackInfos per Track in target”, that plots per each track reaching the target how many track informations are kept in the histograms. Another useful histogram for this case may be the 2D histogram “trackInfo Region vs trackInfo Region”, that plots all the region number of all the pairs of track informations that correspond to the same track reaching the target (you can get a list of which region number corresponds to which region at the end of the standard output file).

To use this utility in GAMOS it is only needed to add this command in your script:

/gamos/userAction GmProdCutsStudyUA RTCutsStudyFilter

or that will use as target condition that a track reaches a plane perpendicular to the Z axis defined with the parameters

/gamos/setParam RTCutsStudyFilter:PlaneZ ZPOS

/gamos/setParam RTCutsStudyFilter:PlaneXDim XDIM

/gamos/setParam RTCutsStudyFilter:PlaneYDim YDIM

This command will produce at the end of run a table and a histogram file with the needed information. The table will contain the minimum range that can be applied for each region/particle/process not to lose any track reaching the target, and it will look like this

%%%%% PRODUCTION CUTS STUDY RESULTS
GmProdCutsStudyUA: REGION= DefaultRegionForTheWorld PARTICLE=
gamma PROCESS= ALL MIN RANGE= 353161.38
GmProdCutsStudyUA: REGION= DefaultRegionForTheWorld PARTICLE=
gamma PROCESS= eBrem MIN RANGE= 353161.38

To get the cuts values for not losing a given percentage of particles in the target plane you can execute the ROOT script that can be found at GamosCore/GamosPhysics/GamosCuts/getProdCutsEffect.C :

*root -b -p -q .x getProdCutsEffect.C++\\(\\"prodcuts.root\\",percentage\\)*

and look at the last lines of output, those that contain the word ‘FINAL’, like the following ones

PARTICLE= e+ FINAL= 17 / 19
PARTICLE= e- FINAL= 72 / 34185
PARTICLE= gamma FINAL= 72 / 34184

Automatic determination of production cuts for a dose in a phantom simulation

The use of production cuts in a dose computation may introduce a bias when a particle is killed (and then its energy is deposited locally) and it has enough energy to reach the next phantom voxel, or enough energy to create a particle that reaches the next phantom voxel (this happens mainly for electrons creating gammas, which have a much higher range).

To calculate automatically the best production cut, that is the one that gives the smallest CPU while biasing the dose computation a minimal amount, GAMOS uses an inverse reasoning. For a given set of cuts for electron and gamma it does not apply them but tags the particles that would have been killed by them. It also tags the voxel in which the particle is produced and then it computes all the dose deposited by the tagged particle or any of its children in a voxel that is not the same as the tagged voxel.

To use this utility in GAMOS you just have to associate to your dose scorer a filter of type GmProdCutOutsideVoxelFilter, passing to it as arguments the gamma cut and the electron cut, like in the following example

/gamos/scoring/addFilter2Scorer ProdCutFilter GmProdCutOutsideVoxelFilter PDDscorerPC10.1. 10.*mm 1.*mm

You should add another scorer without filter to get the total dose. After running your job with as many scorer-filter combinations as you like, you can look at the total dose deposited by each scorer and compare it with the total dose. It may happen that the dose lost with certain cuts, despite being a small proportion of the total dose, is distributed in a different manner than the total dose, introducing some bias in some region that you consider not acceptable. To check in detail the dose produced with a certain filter you can add a scorer printer of type RTPSPDoseHistos, that will produce several histograms of the dose (PDD, X & Y profiles, dose, dose-volume):

/gamos/scoring/addPrinter2Scorer RTPSPDoseHistos DoseScorerPC10.1.

The name of the printer will be passed to the name of the file containing the histograms.

Automatic determination of user limits for an accelerator simulation

The method used in GAMOS to determine the minimum range user limits is similar to the one used to determine the best production cuts. The main difference is that when a track reaches the target we do not have to look at the range it had when created, but at the range it had in every step. This is because even if we want the minimum step, the track may have crossed several regions and the smallest range may not correspond to the last step. What we do nevertheless is only consider the last step when there are a set of contiguous steps in the same region. Also for the ancestor tracks we have to store the information of each step, starting of course with the one when the track that reached the target (or its n-th ancestor if we are looking at the (n+1)-th ancestor) was created.

The same warning as for the production cuts should be mentioned here, but in this case it is more than a mere warning: when a track reaches the target, we accumulate one-track information of the last step in each region, for each of the ancestor tracks. Therefore it is very likely that there are more than one track information per track reaching the target, and therefore there will be overcounting of the number of tracks killed by a cut. As for the production cuts you should keep an eye on this.

To use this utility in GAMOS it is only needed to add the command in your script:

/gamos/userAction GmMinRangeLimitsStudyUA RTCutsStudyFilter

what will produce at the end of run a table and a histogram file with the needed information.

root -b -p -q .x getMinRangeCutsEffect.C++\\(\\"prodcuts.root\\",percentage\\) \|& tee out

and look at the last lines of output, those that contain the word ‘FINAL’, like the following ones

PARTICLE= e+ FINAL= 17 / 19
PARTICLE= e- FINAL= 72 / 34185
PARTICLE= gamma FINAL= 72 / 34184

Automatic determination of user limits for a dose in phantom simulation

The method used in GAMOS to determine the minimum range user limits is similar to the one used to determine the best production cut.

To use this utility in GAMOS you just have to associate to your dose scorer a filter of type GmProdCutOutsideVoxelFilter, passing to it as arguments the gamma cut and the electron cut, like in the following example

/gamos/filter ProdCutFilter GmMinRangeCutOutsideVoxelFilter 10.*mm 1.*mm

/gamos/scoring/addFilter2Scorer ProdCutFilter DoseScorerPC10.1.

After running your job with as many scorer-filter combinations as desired, you can look at the total dose deposited and compare it with the dose with very small cuts. It may happen that the dose with certain cuts, despite being a small number, is distributed in a different manner than with very small cuts, introducing some bias that you consider not acceptable. To check in detail the dose produced with a certain filter you can add a scorer printer of type RTPSPDoseHistos, that will produce several histograms of the dose (PDD, X & Y profiles, dose, dose-volume):

/gamos/scoring/addPrinter2Scorer RTPSPDoseHistos DoseScorerPC10.1.

The name of the printer will be passed to the name of the file containing the histograms.