9   Scoring

9.1   Creating a scorer

Geant4 provides several classes to score different quantities in the selected volumes. GAMOS provides all the Geant4 functionality through user commands and also some extra functionality that we describe in this section.

The first thing you should do to use GAMOS scoring is to create a multi-functional detector [biblio.g4doc-scoring] and associate it with a list of logical volumes, with the user command

/gamos/scoring/createMFDetector MFD_NAME LOGICAL_VOLUME_NAME(s)

where MFD_NAME is the detector name that will be used later. and LOGICAL_VOLUME_NAME(s) is a list of logical volumes that you associate to this detector.

Then you should add to the detector one of the GAMOS scorers, or your own ones, with the user command

/gamos/scoring/addScorer2MFD SCORER_NAME SCORER_CLASS MFD_NAME

SCORER_NAME is a name you give to the scorer to be used later, SCORER_CLASS is one of the available scorer classes, MFD_NAME is one of the multi-functional detectors created above.

If the scorer class needs some extra parameters (see scorer classes below), you cannot use it directly, but you have first to create a scorer:

/gamos/scoring/createScorer NEW_SCORER_NAME SCORER_CLASS SCORER_PARAMETERS

and then attach the scorer to the detector

/gamos/scoring/addScorer2MFD NEW_SCORER_NAME SCORER_CLASS MFD_NAME

You may repeat these commands to associate several scorers to the same detector.

To each of the defined scorers you can add a filter, to select for which track conditions the scoring will be done:

/gamos/scoring/addFilter2Scorer FILTER_NAME/CLASS SCORER_NAME

FILTER_NAME/CLASS is the name of a GmVFilter class or the name you gave to a filter built from a filter class by using the command /gamos/filter (see section on Filters) and SCORER_NAME is one of the scorers defined above.

You may repeat this command to associate several filters to the same scorer. See section on Filters for a description of the available filters and how to create your own one.

By default a different count is scored for each of the copies of the selected volumes with different copy number. This is managed by the scorer classifier GmScorerClassifierBy1Ancestor. The user can attach different classifiers to the different scorers so that the counts are done in different ways:

/gamos/scoring/assignClassifier2Scorer CLASSIFIER_NAME/CLASS SCORER_NAME

CLASSIFIER_NAME/CLASS is the name of a GmVClassifier class or the name you gave to a classifier built from a classifier class by using the command /gamos/classifier (see section on Classifiers) and SCORER_NAME is one of the scorers defined above.

Finally you may select the format of the scoring results by associating one of the available scorer printers for each scorer,

/gamos/scoring/addPrinter2Scorer PRINTER_NAME SCORER_NAME

PRINTER_NAME is the name of a GmVPSPrinter class or the name you gave to a printer built from a printer class by using the command /gamos/printer (see section on Scorer printers below) and SCORER_NAME is one of the scorers defined above.

If no printer is attached to a scorer, it will use the printer type GmG4PSPrinterCout.

The scoring is done by default taking into account the track weight, except for the scorers when it is explicitly mentioned (see scorers description). If you do not want to take weights into account you can switch them off with the command

/gamos/scoring/useTrackWeight SCORER_NAME FALSE

The quantities scored are given per event by default. If you want the score without dividing by the number of events, use the command

/gamos/scoring/printByEvent PRINTER_NAME FALSE

The error in the scored quantity per voxel is also calculated by default, using the following formula:

sqrt{ (SumV2*nEvents - SumV*SumV) / (nEvents-1)} / nEvents

where nEvents is the total number of events in the run, SumV is the sum of the scored quantity value times its weight (i.e. the scored quantity itself) and SumV2 is the sum of squares of scored quantity value times its weight.

When the scoring is done per event, this sum of squares is done summing first all the values times its weight belonging to all the particles of the same event and then squaring this quantity. In this way the correlations between particles of the same event is properly taken into account. If the scoring is not by event, the error calculation uses the same formula, but the sum of squares is not done summing the contribution of the particles of the same event, but squaring each contribution individually.

Calculating the errors makes it necessary to store the square of the weights, increasing substantially the memory usage and CPU time. If you want to deactivate this option for a scorer, use the command

/gamos/scoring/scoreErrors SCORER_NAME FALSE

SCORER_NAME is one of the scorers defined above.

As mentioned above the errors that are calculated taking into account the number of events. You have to be careful then if you set to off the option of scoring by event and keep on the option of calculating the errors. In the default GAMOS scorer printers, the errors are printed are relative, i.e. the error divided by the value, so no caution is necessary, but be careful if you define a printer yourself.

Each scorer has a default unit (see section below), but it can be overridden with the command:

/gamos/scoring/scorerUnit SCORER_NAME UNIT_NAME UNIT_VALUE

9.2   Scorer classes

All the available scorers in Geant4 are also available in GAMOS. The classes have been slightly changed to provide the extra functionality.

The scorers can be classified in the following types:

  • Track length scorers - GmG4PSTrackLength: The track length is defined by the sum of step lengths of the particles inside the cell (i.e., the volume where the scoring happens). A particle weight is not applied by default. There are two extra parameters, that are FALSE by default and can be set TRUE or FALSE: to multiply by the kinetic energy and to divide by the velocity.

    If the energy track flux is required then you should set them to TRUE FALSE. Alternatively to measure the flux per unit velocity then you should set them to FALSE TRUE. Finally to measure the flux energy per unit velocity then you should set them to TRUE TRUE.

    • GmG4PSPassageTrackLength: The passage track length is the same as the track length in GmG4PSTrackLength, except that only tracks which pass through the volume are taken into account. This means that newly-generated or stopped tracks inside the cell are excluded from the calculation. A particle weight is not applied by default.
  • Deposited energy scorers - GmG4PSEnergyDeposit: This scorer stores a sum of particles’ energy deposits at each step in the cell. - GmG4PSDoseDeposit: In some cases, dose is a more convenient way to evaluate the effect of energy deposit in a cell than simple deposited energy. The dose deposit is defined by the sum of energy deposits at each step in a cell divided by the mass of the cell. The mass is calculated from the density and volume of the cell taken from the methods of G4VSolid and G4LogicalVolume. - GmG4PSKerma: This scorer stores a sum of particles’ kerma each step in the cell, that is the energies of the charged secondary particles produced by non-charged particles, divided by the mass of the cell. The following scorers serve to calculate the deposited dose in a volume dividing it virtually in smaller pieces. This means that the energy deposited is not divided by the volume the track is in, but by the division volume as defined by the scorer. As you may find other type of use for these scorers, you can choose also to divided by the total volume defined with the scorer, by setting the parameter: /gamos/setParam SCORER_NAME:UseTotalVolume 1 - GmPSSphericalDoseDeposit: This scorer is designed to compute the dose deposit and assign a different index to each spherical radius interval, despite no spherical division is defined in the geometry. It creates intervals as a function of spherical radius with the parameters given at score definition: NSLICES_R MIN_R MAX_R and calculates the volume of each spherical shell so that dose can be computed. Three extra parameters may be given: CENTRE_X CENTRE_Y CENTRE_Z, if the you do not want the sphere to be centred at the origin of coordinates, (0,0,0). - GmPSCylindricalRZDoseDeposit: This scorer is designed to compute the dose deposit and assign a different index to each interval in cylindrical coordinates (in 2D radius and Z), despite no cylindrical division is defined in the geometry. It creates intervals as a function of the cylindrical radius and the Z coordinate with the parameters given at score definition: NSLICES_R MIN_R MAX_R NSLICES_Z MIN_Z MAX_Z and calculates the volume of each cylindrical section so that dose can be computed. Three extra parameters may be given: CENTRE_X CENTRE_Y CENTRE_Z, if you do not want the cylinder to be centred at the origin of coordinates, (0,0,0), and another three: ZAXIS_X ZAXIS_Y ZAXIS_Z, if the you do not want the cylinder axis to be along the Z global axis, (0,0,1).

    The index of the scorer is build as NSLICES_R*indexZ + indexR.

    • GmPSCylindricalRPhiDoseDeposit: This scorer is designed to compute the dose deposit and assign a different index to each interval in cylindrical coordinates (in 2D radius and phi angle), despite no cylindrical division is defined in the geometry. It creates intervals as a function of the cylindrical radius and the phi angle with the parameters given at score definition: NSLICES_R MIN_R MAX_R NSLICES_PHI MIN_PHI MAX_PHI and calculates the volume of each cylindrical section so that dose can be computed. Tow more parameters are mandatory, to define the dimension in Z: MIN_Z MAX_Z. Three extra parameters are optional: CENTRE_X CENTRE_Y CENTRE_Z, if you do not want the cylinder to be centred at the origin of coordinates, (0,0,0), and another three: ZAXIS_X ZAXIS_Y ZAXIS_Z, if the you do not want the cylinder axis to be along the Z global axis, (0,0,1). The index of the scorer is build as NSLICES_R*indexPHI + indexR.
    • GmPSCylindricalZPhiDoseDeposit: This scorer is designed to compute the dose deposit and assign a different index to each interval in cylindrical coordinates (in 2D radius and phi angle), despite no cylindrical division is defined in the geometry. It creates intervals as a function of the cylindrical radius and the phi angle with the parameters given at score definition: NSLICES_Z MIN_Z MAX_Z NSLICES_PHI MIN_PHI MAX_PHI and calculates the volume of each cylindrical section so that dose can be computed. Tow more parameters are mandatory, to define the dimension in radius: MIN_R MAX_R. Three extra parameters are optional: CENTRE_X CENTRE_Y CENTRE_Z, if you do not want the cylinder to be centred at the origin of coordinates, (0,0,0), and another three: ZAXIS_X ZAXIS_Y ZAXIS_Z, if the you do not want the cylinder axis to be along the Z global axis, (0,0,1). The index of the scorer is build as NSLICES_Z*indexPHI + indexZ.
  • Current and flux scorers - There are two different definitions of a particles flow for a given geometry. One is a current and the other is a flux. In our scorers, the current is simply defined as the number of particle’s (with the particle’s weight) passing through a certain surface or volume, while the flux takes the particle’s injection angle to the geometry into account (dividing the area of the surface that is traversed and the cosine of the angle between the track direction and the surface normal). The current and flux are usually defined at a surface, but volume current and volume flux are also provided. - GmPSSurfaceFlux: This is a surface based scorer.

    The quantity is the flux, i.e. the number of tracks that reach the surface, divided by the surface area and the cosine of the angle between the track and the surface normal. Instead of the flux, the user may choose not to divide by the area, by setting the parameter: /gamos/setParam SCORER_NAME:DivideByArea 0 And also not to divide by the cosine of the angle between the track direction and the surface normal, by setting the parameter): /gamos/setParam SCORER_NAME:DivideByAngle 0 It serves to score the flux in the surfaces of a box, tube or sphere. The user can choose the direction of the particle to be scored with the parameter /gamos/setParam SCORER_NAME:Direction DIRECTION The choices are In, Out or InOut. In scores incoming particles to the volume, while Out scores only outgoing particles from the volume. InOut scores both directions. The user should also define the list of volume surfaces in which the scoring will be done: /gamos/setParam SCORER_NAME:Surfaces SURFACE_1 SURFACE_2 … These surfaces depend on the volume solid: Box - X+: YZ plane at positive X - X-: YZ plane at positive X - Y+: XZ plane at positive Y - Y-: XZ plane at positive Y - Z+: XY plane at positive Z - Z-: XY plane at positive Z Tube - INNER: inner radius surface - OUTER: outer radius surface - PHI: if the full phi (360 degrees) is not defined, it is the two surfaces composed of the points of same phi that limit the tube volume - UP: the upper surface - DOWN: the lower surface Sphere - INNER: inner radius surface - OUTER: outer radius surface - THETA: if the full theta (180 degrees) is not defined, it is the two surfaces composed of the points of same theta that limit the sphere volume - PHI: if the full phi (360 degrees) is not defined, it is the two surfaces composed of the points of same phi that limit the sphere volume

    • GmPSSurfaceFlux: This is a surface based scorer. The quantity is the flux, i.e. the number of tracks that reach the surface, divided by the surface area and the cosine of the angle between the track and the surface normal.
    • GmG4PSPassageCellCurrent: Passage current is a volume-based scorer. The current is defined by the number of tracks that pass through the volume.
    • GmG4PSCellFlux: Cell flux is a volume based flux scorer. The cell flux is defined by a track length (L) of the particle inside a volume divided by the volume (V) of this cell. The track length is calculated by a sum of the step lengths in the volume. The expression for cell flux is given by the sum of (W*L)/V, where W is a particle weight, and is multiplied by the track length at each step. This scorer is also called GmPSVolumeFlux.
    • GmG4PSPassageCellFlux: Passage cell flux is a volume based scorer similar to G4PSCellFlux. The only difference is that tracks which pass through a cell are taken into account. It means that tracks generated or stopped inside the volume are excluded from the calculation.
  • In/Out behaviour For the scorers GmPSSurfaceFlux and GmG4PSTrackCounter you can make the scoring only for tracks that are entering, only for tracks that are exiting or both for tracks that are entering or exiting (default behaviour). To select among these three options you can add an extra parameter when defining the scorer that can be In, Out or InOut

  • Other scorers - GmG4PSMinKinEAtGeneration: This scorer records the minimum kinetic energy of secondary particles at their production point in the volume in an event. This primitive scorer does not integrate the quantity, but records the minimum quantity. - GmG4PSNofSecondary: scores the number of secondary particles generated in the volume. A particle weight is not applied by default.

    The user can choose if the scoring is done for all types of particles (default) or only for a set of particles, by naming them as extra parameters.

    • GmG4PSNofStep: scores the number of steps in the cell. A particle weight is not applied by default. If an extra parameter is set to TRUE those steps with step length zero will not be taken into account.
    • GmG4PSCellCharge: scores the total charge of particles which have stopped or have been created in the volume, i.e. the tracks that enter count as +1 and the tracks that exit count as -1.
    • GmG4PSTrackCounter: scores the number of tracks in a cell.
  • For Event Biasing Scoring for event biasing is a very specific use case whereby particle weights and fluxes through importance cells are required. The goals of the scoring technique are to: - appraise particle quantities related to special regions or surfaces, - be applicable to all “cells” (physical volumes or replicas) of a given geometry, - be customizable. .. COMMENT: %Standard scoring is provided for quantities such as tracks entering a cell, average weight of entering tracks, energy of entering tracks, and collisions inside the cell. A number of scorers have been created for this specific application: - GmG4PSNofCollision: This scorer records the number of collisions that occur within a scored volume/cell. - GmG4PSPopulation: This scores the number of tracks within in a given cell per event. A particle weight is not applied by default. - GmG4PSTermination: This scores the number of tracks that are terminated in a given cell per event. A particle weight is not applied by default.

  • From data - GmG4PSData: This score uses a GAMOS data as quantity to score. See section on GAMOS data.

All scorers that must divide by the volume to obtain the relevant quantity use by default the cubic volume of the physical volume in which the particle is being propagated. If you want to use instead the total volume of all copies of the logical volume with the same name as the volume the particle is in, you have to set the parameter:

/gamos/setParam SCORER_NAME:IntegrateVolumes 1

9.2.1   Scoring in voxelised phantoms

When you are navigating in a voxelised phantom geometry and you are using the regular navigation you have the option of skipping voxel frontiers when the material does not change (see section on Reading DICOM files). In this case the scorer GmG4PSDoseDeposit will automatically distribute in the voxels traversed the dose corresponding to the energy deposited along each step. A correction is nevertheless needed: when the particle travels it loses energy, and therefore it is not right to distribute the dose proportionally to the length traversed inside each voxel. The needed corrections to the multiple scattering and energy loss when the particle loses energy are done in an iterative way: the increase of multiple scattering increases the path length, what increases the energy lost, what makes the multiple scattering even bigger, what increases the path length and therefore the energy lost. By default the number of iterations is two, but you may change it with the parameter:

/gamos/setParam GmEnergySplitter:NIterations NITER

Other scorers need the same corrections, namely GmG4PSEnergyDeposit, GmG4PSEnergyLost, GmG4PSKerma and GmG4PSTrackLength.

In this case the index used is the copy number of the voxels traversed, what means that the classifier index is not used at all. You may nevertheless force the use of the classifier index with the parameter:

/gamos/scoring/printer SCORER_NAME:UseClassifierIndex 1

but you have to be aware that if you do so the score will not be distributed in the voxels traversed, even if the classifier is GmClassifierBy1Ancestor:

Other scorers are also affected by the feature of skipping voxel boundaries:

  • GmG4PSNofSecondary: the score is done at the last voxel traversed, while using the default classifier GmClassifierBy1Ancestor would assign it to the one corresponding to the origin of the step, that is the first voxel.
  • GmG4PSNofCollision: the same as GmG4PSNofSecondary.
  • GmG4PSNofStep: a score is assigned to each one of the voxels traversed.
  • GmG4PSTrackCounter: for the first step it is always that the track exited; for intermediate voxels it is always that the track entered and exited; for the last voxel it is always that the track entered.

The flux and current scorers are not corrected, as it likely has no sense to calculate the flux entering and exiting voxels. Please contact the GAMOS team if you think you need this feature.

9.3   Compound scorer

The compound scorer allows to score in any arithmetical expression of other scorers. For example

/gamos/scoring/createScorer myScorer GmCompoundScorer GmG4PSEnergyDeposit/GmG4PSTrackLength

calculates for each classifier index (for example for each voxel if you are scoring in a phantom with classifier GmClassifierBy1Ancestor) the sum of the energy deposited of all the steps of all the tracks in an event and divides it by the sum of track lenghts.

This is different than using:

/gamos/scoring/createScorer myScorer GmPSData GmG4PSEnergyDeposit/GmG4PSTrackLength

In this case for each step in the scoring volume the division of the energy deposited by the track length, and therefore it will only give the same result as the previous scorer if there is never two steps in an event contributing to the scoring that are classified by the same index.

One example of the utility of the compound scorer are the Linear Energy Transfer scorers.

9.3.1   Linear Energy Transfer (LET) scorer classes

Apart from defining the LET scorer with the command above :

/gamos/scoring/createScorer myScorer GmPSData GmG4PSEnergyDeposit/GmG4PSTrackLength

there are other ways to define it used in dosimetry:

  • GmPSLETByEDepPhi:

LET = {\sum_{i}^{N}dx_i \times LET_i \over \sum_{i}^{N}dx_i}

where dx_i is the step length

  • GmPSLETByEDepD:

LET = {\sum_{i}^{N}E_{dep,i} \times LET_i \over \sum_{i}^{N}E_{dep,i}}

where E_{dep,i} is the energy deposited in the step

  • GmPSLETBydEdxPhi:

LET = {\sum_{i}^{N}dEdx_i \times LET_i \over \sum_{i}^{N}dEdx_i}

where dEdx_i is the electronic variation of energy per unit length, computed using G4EmCalculator, at the initial point of the step

  • GmPSLETBydEdxD:

LET = {\sum_{i}^{N}dEdx_i \times LET_i \times LET_i \over \sum_{i}^{N}dEdx_i \times LET_i}

  • GmPSLETESpectPhi:

LET = {\sum_{i}^{N}\Phi_i \times LET_i \over \sum_{i}^{N}dx_i}

where \Phi_i is the particle fluence as a function of kinetic energy

  • GmPSLETESpectD:

LET = {\sum_{i}^{N}\Phi_i \times LET_i \times LET_i \over \sum_{i}^{N}\Phi_i \times LET_i }

9.4   Filters

See section on Filters.

9.5   Scorer printers

A scorer printers serves to select the format of the output of a scorer. These classes are unique to GAMOS, as Geant4 does not provide this functionality. As mentioned above, several printers can be associated to the same scorer.

To associate a printer to an scorer, you can use the command mentioned above:

/gamos/scoring/addPrinter2Scorer PRINTER_NAME SCORER_NAME

The name of the printer is the name of the printer class (see list below). If you want to change the printer name, or you want to add some parameters to a printer, you can use the command

/gamos/scoring/printer PRINTER_NAME PRINTER_CLASS PARAMETERS

where PRINTER_NAME is the new name you want to give to the printer and PRINTER_CLASS is the name of the printer class (see list below).

The general-use scorer printers currently in GAMOS are the following:

  • GmPSPrinterCout: Prints in the standard output the summary of scoring:

    • “MultiFunctionalDet: ” Detector name
    • “PrimitiveScorer: ” Scorer name
    • “Number of entries: ” Number of scores
    • For each per score one line with: “index: ” IndexName ” = ” Score value ” +-(REL) ” Relative score error (=error/value) ” sumV2= ” Square of score values. This last value is needed instead of the error for properly taking into account the correlations when the errors of several files are summed. If you do not want it to be printed you have to set the parameter: /gamos/setParam PRINTER_NAME:PrintSumV2 1
    • ScorerName ” SUM ALL: ” Sum of all scores ” +-(REL) ” Relative error ” ” Unit name

    An example is the following:

    MultiFunctionalDet: fluxDet
    PrimitiveScorer: fluxScorer
    Number of entries= 2
    index: 1  = 0.0422424 +-(REL) 0.02251 cm-2 sumV2= 108.251
    index: 2  = 0.010207 +-(REL) 0.0710814 cm-2 sumV2= 53.6754
    fluxScorer SUM ALL: 0.0524494 +-(REL) 0.00119606 cm-2
    

    where MFD_NAME is the name of multi-functional detector and SCORER_NAME is the name of the scorer. The columns after index have the following meaning: Index, scorer value, scorer error (relative, i.e. error/value), unit name. At the end it is printer the sum of all scores, with error (non relative).

  • GmPSPrinterTextFile: Prints in a text file the same information than the GmPSPrinterCout printer. The name of the file is PRINTER-NAME.out. It can be changed with the parameter /gamos/setParam PRINTER-NAME_SCORER-NAME:FileName FILE-NAME

  • GmPSPrinterBinFile: Prints in a binary file almost the same information than the GmPSPrinterCout printer. The format of the data is the following:

    • Detector name (char[30])
    • Scorer name (char[20])
    • Unit name (char[10])
    • Unit value (float)
    • Number of entries (unsigned int)

    For each classifier index: - Index name (char[20]) - Score value (float) - Square of score value (float). This value is needed instead of the error for properly taking into account the correlations when the errors of several files are summed. The name of the file is PRINTER_NAME.out It can be changed with the parameter /gamos/setParam PRINTER-NAME_SCORER-NAME:FileName FILE-NAME

  • GmPSPrinterHistos: Prints scoring dat in histograms. The name of the file is PRINTER_NAME.root. It can be changed with the parameter /gamos/setParam PRINTER_NAME:FileName FILE_NAME The user must define the histogram name, number of bins and axis minimum and maximum, giving them as arguments of the /gamos/scoring/printer command (this means that the class GmPSPrinterHistos cannot be used directly). The scoring data corresponding to index IDX will be used to set the content of the histogram bin IDX. It should be noted that the first histogram bin is number 1, while it may happen that the first index is not 1. If this is the case, a number (positive or negative) can be added to the index by defining an offset as a parameter:

    /gamos/setParam PRINTER_NAME:OffsetX OFFSET

    If and only if you are using a GmCompoundClassifier built from two classifiers, you may want to produce a 2-dimensional histogram. In this case you should give eight instead of four parameters to the /gamos/scoring/printer command: the first four should correspond to the data of the first classifier and the second four to the second one. Then three histogram will be filled: the compound index will be split in the two indices (as defined with the NShift of the GmCompoundClassifier); then an histogram will be filled for the first index, another for the second one and a 2-dimensional histogram combining the two indices. A second offset may be defined for the second index: /gamos/setParam PRINTER_NAME:OffsetY OFFSET

Other scorer printers are provided for specific applications. See corresponding sections in this guide.

9.6   Classifiers

See section on Classifiers.

9.7   Multiplying by data

You can multiply the quantity you are scoring by any of the GAMOS data (see section on GAMOS data), so that before scoring the quantity value it will be multiplied by the value of that the GAMOS data has in that step. To do it you have to use the parameter:

/gamos/setParam SCORER_NAME:MultiplyByData DATA_NAME

where DATA_NAME is the name of the data you want to use.

9.8   Multiplying by a distribution

You can multiply the quantity you are scoring by a GAMOS distribution (see section on GAMOS distributions), so that before scoring the quantity value it will be multiplied by the value of that the GAMOS distribution has in that step. To do it you have to use the parameter:

/gamos/setParam SCORER_NAME:MultiplyByDistribution DISTRIBUTION_NAME

where DISTRIBUTION_NAME is the name of a distribution you have previously created.

9.9   Convergence testing

The fact of having a small relative error in a score does not always guarantee that it has converged to a good result. This may be specially true when there are contributions of very different weights to the scorer, and the high weight scores have not been properly sampling with your statistics. If you suspect of a misbehaviour you can activate the convergence test with the parameter (as these tests consume some CPU time and memory, by default they are deactivated):

/gamos/setParam SCORER_NAME:ConvergenceTester CONVERGENCE_NAME

where SCORER_NAME is the name of the scorer and CONVERGENCE_NAME is the name you want to give to the convergence tester, which will be used in the report. If you are using a point detector scorer, you should substitute SCORER_NAME by GmPDS.

The convergence tests are taken from the Geant4 code, and are inspired in the MCNP tests. The tests are based in the analysis of the sum of scores at the end of each event. The following variables are printed about the sum of scores :

  • EFFICIENCY = Proportion of events that have a non zero value.
  • MEAN = Average value*
  • VAR = variance of values
  • SD = standard deviation
  • R = The estimated relative error = SD / MEAN / sqrt(number of values)
  • SHIFT = Shift in the midpoint of the confidence interval to a higher value = ( SUM(Xi-MEAN)*(Xi-MEAN)*(Xi-MEAN) - (N-Number_of_non_zero_values)*MEAN*MEAN*MEAN ) / (2*VAR*N)
  • VOV = Variance of variance = ( SUM(Xi-MEAN)*(Xi-MEAN)*(Xi-MEAN)*(Xi-MEAN) + (N-Number_of_non_zero_values)*MEAN*MEAN*MEAN*MEAN ) / (VAR*VAR) -1./N
  • FOM = Figure of merit = 1 / (R*R) / CPU_time_of_last_event
  • THE LARGEST VALUE and where in which event it happened

To get a feeling of how big are the fluctuations the same variables, i.e. mean, variance, R, shift and FOM, are printed again but adding to the values a new one equal to the largest value (so that this value is counted twice). And also the ratio of this affects to the original ones.

Then the results of eight convergence tests are shown:

  • MEAN distribution is RANDOM
  • r follows 1/std::sqrt(N)
  • r is monotonically decrease
  • r is less than 0.1. This value can be changed with the parameter: /gamos/setParam GmConvergenceTester:Rmin VALUE
  • VOV follows 1/std::sqrt(N)
  • VOV is monotonically decrease
  • FOM distribution is RANDOM
  • SLOPE is large enough

Finally it prints the evolution of several variables, i.e. the change of these variables when more events are added. The variables are printed each N/16 events, where N is the total number of events. The variables printed are the following:

  • mean
  • var
  • sd
  • r
  • vov
  • shift
  • e
  • r2eff
  • r2int

9.10   Point detector scorer

The point detector scorer covers the problems where the quantity to be calculated is the flux or the dose of neutral particles (neutrons or gammas) in a very small detector (small with respect to the setup dimensions) situated far from the primary particles source, of behind a thick shielding layer. In this kind of problems the fraction of particles that reach the detector is very small, and therefore if one wanted to calculate the flux or dose by conventional methods, the statistics needed would be prohibitive.

The technique implemented in GAMOS is similar to the F5 tally implemented in MCNP. It is based on the following idea: normal tracks are propagated and for each neutron or gamma interaction is calculated the probability that it would be deviated in the direction of the point detector (instead of the real direction towards which it is deviated) multiplied by the probability that it reaches the point detector without any further interaction. The first probability is based on a precalculated table of angle probabilities for each interaction type, each material and each energy, multiplied by the solid angle covered by the point detector. The second probability is based of the number of mean free paths that the neutron or gamma would have to traverse along the path to reach the point detector. A similar calculation is done for each neutron or gamma initial step, with the difference that the angle probability is simplified as a constant distribution in the local reference system of the parent particle (i.e. probability equal to 0.5). Summing up these probabilities for a number of events, usually several orders of magnitude less than with conventional methods, one can get the flux and the energy spectrum at the point detector, and with these magnitude the equivalent doses (H*, Hp(0,15,…)) can be obtained.

In this chapter we describe first the theoretical basis of the method. Then we explain how it is implemented in GAMOS, how to select the different available options, and which are the possible outputs (tables or histograms). In the third part we describe the several variance reduction techniques that can be used to reduce the CPU time, and give some recommendations on how to use the output table and histograms to help in determining the best variance reduction techniques for a given problem.

9.10.1   Theoretical basis

The point detector scorer calculates the flux at a point based on the probability that particles reach it. Suppose we call p(\mu,\phi) the probability that at a source or neutron/gamma interaction the the particle produced or scattered goes into the solid angle d\Omega about the direction (\mu,\phi), where \mu = cos(\theta). if R is the distance to the detector from the source or interaction point

p(\mu,\phi)d\Omega e^{ \int_{o}^{R} \sum_{t}(s)ds }

yields the probability that the particle reaches the detector point with no further collisions, where

e^{ \int_{o}^{R} \sum_{t}(s)ds }

is the attenuation of a beam of mono-energetic particles passing through a material medium, where s is measured along the direction from the source or interaction point to the detector and \Sigma(s) is the macroscopic total cross section at s. If dA is an element of area normal to the line of flight to the detector, d\Omega = dA / R^2. As the flux is the number of particles passing through a unit area normal to the line of flight direction, we can write the flux as

\Phi = p(\mu,\phi) e^{\lambda} / R^{2}

If there is azimuthal geometry, then p(\mu,\phi) = p(\mu) / 2\phi, and we can finally write

The attenuation term, d\Omega = dA / R^2. As the flux is the number of particles passing through a unit area normal to the line of flight direction, we can write the flux as

\Phi = p(\mu,\phi) e^{\lambda} / R^{2}

If there is azimuthal geometry, then p(\mu,\phi) = p(\mu) / 2\phi, and we can finally write

The attenaution term, e{\lambda} is calculated summing up the number of interaction lengths from the source or interaction point until the detector; this calculation depends on the neutron/gamma energy and the materials traversed, and it is computed in GAMOS through the use of geantino

A geantino is a pseudo-particle that has only transportation, but no physical process. .

The R2 term in the denominator causes a singularity: when a source or interaction event occurs near the detector point, R approaches zero and the flux approaches infinity. The technique is still valid and unbiased, but convergence is slower and often impractical. To avoid this singularity, a fictitious sphere of radius R0 surrounding the detector point is defined (which we call the exclusion sphere). If we can assume that the flux is uniformly distributed inside this sphere, then it can be demonstrated that the average flux in the sphere is

\Phi(R<R{0}) = \frac{p(\mu)(1- e^{-\sum_{t}R_{0}})}{\frac{2}{3}\pi R_{0}^{3}\sum_{t}}

The exclusion sphere radius has a default value of 1 mm, which can be changed with the parameter

/gamos/setParam GmPDS:ExclusionRadius RADIUS

9.10.2   GAMOS implementation

The point detector scorer does not really behave as the other scorers, although the output format is similar to then. To activate you hasve to select the user action:

/gamos/userAction GmPDSUA

But before this command you have to set the parameters that drive the behaviour of the point detector scoring. You may first decide if you want to score the flux for neutrons, gammas or for both with the two parameters

/gamos/setParam GmPDS:ScoreNeutrons 1/0

/gamos/setParam GmPDS:ScoreGammas 1/0

There are several parameters that you may change distinctly for neutrons or gammas, or that you may want to set the same value for both. All the parameters start with GmPDS and then they are followed by :neutron or :gamma; if the second word is not one of these two, then the parameter serves for both. We will then omit this word in the explanation of the parameters in this section, but remember you can always add it, except in the cases where it is explicitly mentioned that there is a unique parameter for both particles.

9.10.2.1   Creating angle deviation probability densities

The first thing you should do before running your point detector scoring job, is to create a file that contains the probability density functions of the emission angles for each energy, process and material. Else GAMOS will assume flat probabilities for all angles (what may be enough in your case).

To do this you can run a job that writes in a ROOT file the probability of emission angle for all the materials that constitute your setup. You may follow the example found in tutorials/ShieldingTutorial/exercise3/exercise3.angle.gg.in, that we explain here. The class that makes the work is called GmPDSCreateAngleTablesUA, but remember that the parameters have to go before the class that uses them, so you should add this line at the end of your script

/gamos/userAction GmPDSCreateAngleTablesUA

The histogram will only be saved for the secondary particle type defined with the parameter

/gamos/setParam USER_ACTION_NAME:SecondaryParticleName NAME

The histogram file name is by default angleDeviation, but it can be changed by the parameter

/gamos/setParam USER_ACTION_NAME:HistoFileName NAME

To create the tables for different materials and tables you may use the generator GmGeneratorChangeEnergyAndMaterials (see section on Event generator changing energy and material).

If you have in your setup several particles that produce neutrons (or gammas) you should repeat the process with different primary particle types. Then you may add all histograms in a file by using the ROOT command

hadd MERGED_FILE.root FILE1.root FILE2.root …

The file created in this way should be passed to the point detector scorer.

/gamos/setParam GmPDS:AngleDeviationFileName FILE_NAME

the default name is angleDeviation.neutron.root for neutrons and angleDeviation.gamma.root for gammas.

When this file is read it will be checked if the number of events in each histogram is big enough, else the histogram will not be used, and a constant value of 0.5 will be used. You may control the minimum number of events in the histogram with the parameter:

/gamos/setParam GmPDS:InteractionAngleManager::MinimumNumberOfEntries VALUE

9.10.2.2   Point detector volume

The detector where the flux will be scored has to be a real volume in your geometry. It should be an small volume, as small as you want, although you may decide to make it the size of the real detector in your experiment. It is important to take into account that despite the flux is calculated as the probability that tracks reach the detector, if a track actually reaches it, it will also be counted. The name of the detector volume is given with the parameter

/gamos/setParam GmPDS:DetectorName VOLUME_NAME

If you want to score the flux in several points you may place the detector volume at several places.

9.10.2.3   Energy and dose equivalent bins

The flux will be scored in different energy bins. These energy bins should be defined in a file whose name is

/gamos/setParam GmPDS:EnergiesFileName FILE_NAME

The format of this file is a set of one-column lines with the energy values.

There are several dose equivalent magnitudes that can be calculated, namely the ambient dose and the personal dose at several angles: 0, 15, 30, 45, 60, 76. To select them you have to set the following parameters (not a different one for neutrons and gammas)

/gamos/setParam GmPDS:PrintHstar 1

/gamos/setParam GmPDS:PrintHp0 1

/gamos/setParam GmPDS:PrintHp15 1

/gamos/setParam GmPDS:PrintHp30 1

/gamos/setParam GmPDS:PrintHp45 1

/gamos/setParam GmPDS:PrintHp60 1

/gamos/setParam GmPDS:PrintHp75 1

For each magnitude you have to define the conversion factor for each energy bin. These are defined in a file whose name is given by the parameter

/gamos/setParam GmPDS:Flux2DoseFileName FILE_NAME

This file must contain in each line the energy and the conversion factors in pSv cm:sup:`2` for the magnitudes: H*(10) Hp(10,0) Hp(10,15) Hp(10,30) Hp(10,45) Hp(10,60) Hp(10,75). The default name of this file is Flux2Dose.neutron.ICRU57.lis for neutrons and Flux2Dose.gamma.ICRU57.lis for gammas. These two files contain the conversion factor given by the ICRU57 report, and can be found in the data directory of the GAMOS release. By default these files will also be used to read the list of energy bins to be used for the flux results.

9.10.2.4   Scoring with filters

As for any user action, you may add one or more filters so that only track steps that pass them will be scored. The filter is applied to the neutron or gamma that creates the geantino, not to the geantino that reaches the detector.

By default the filter will act on the G4Track, but it may act on the G4Step if the following parameter is set.

/gamos/setParam GmPDS:FiltersOnTrack 0

9.10.2.5   Scoring per category

The scoring of flux and equivalent dose can be split following different criteria.

The first possibility is to score separately the contributions from the real particles that reach the point detector and the indirect contributions calculated with the geantinos. The first one will have in the score name the word neutron or gamma, while the second one will have the word geantino. To activate this option the following parameter should be used

/gamos/setParam GmPDS:ScoreSeparatelyTrueAndGeantino 1

The second possibility is use a classifier and get a different score for each classifier index. You can add a classifier to GmPDSUA*as you do for any other user action. By default the classifier will act on the *G4Track, but it may act on the G4Step if the following parameter is set.

/gamos/setParam GmPDS:ClassifierOnTrack 0

9.10.2.6   Output format

At the end of the run the results are written in the standard output or in a file given by the parameter

/gamos/setParam GmPDS:ResultsFileName FILE_NAME

The following table is written for each score or neutrons or gammas (one per point detector copy, one per category as describe above). First a line with the filter names, classifier name, particle type, detector name and copy, score type (named ALL if it is not a sub category), the detector position and the number of events in the job. Then follows a line for each energy with the filter names, classifier name, particle type, detector name and copy, score type, energy value, flux value, relative flux error (error divided by value), number of particles contributing to the flux and then, for later statistical processing, sum of flux values squared, and to the third and fourth power. After all energies comes a line with the total sum of flux in all energies, which has the filter names, classifier name, particle type, detector name and copy, score type, the words FLUX_TOTAL/particle followed by the total flux value and relative error and the total number of particles contributing to the flux. If the dose equivalent magnitudes are required, they come at the end in a line starting also with the filter names, classifier name, particle type, detector name and copy and score type followed by the magnitude name (Hstar, Hp(10,0), Hp(10,15), …). An example output is the following

%%%%%% SCORE IN POINT DETECTOR FOR set ALL: NeutronInelastic: gamma:
   PD1: ALL: at (4800,0,0) NEVENTS= 1000
ALL: NeutronInelastic: gamma: PD1: ALL: ENERGY= 0.01 FLUX= 0 +-(REL) 0
   N 6 Fwei2 0 Fwei3 0 Fwei4 0
ALL: NeutronInelastic: gamma: PD1: ALL: ENERGY= 0.015 FLUX= 0 +-(REL) 0
   N 0 Fwei2 0 Fwei3 0 Fwei4 0
...
ALL: NeutronInelastic: gamma: PD1: ALL:  FLUX_TOTAL/particle=
   1.33272e-13 +-(REL) 0.97118  N 12
ALL: NeutronInelastic: gamma: PD1: ALL: Hstar= 1.49381e-13 pSv/particle

9.10.2.7   Control Histograms

There are several optional histograms that may help you in better understanding the behaviour of your point detector scoring. By default they are created, but you can switch them off with the parameter

/gamos/setParam GmPDS:ControlHistograms 0

The histograms are written in a file named pds, but this name can be changed with the parameter:

/gamos/setParam GmPDS:HistosFileName VALUE

All histograms start with the filter names, classifier name and particle type. There is an histogram about the particle interaction.

  • interaction log Energy (MeV): Kinetic energy of neutrons/gamma at each step initial position.

A set histograms are about the distance to each detector, so they also have the detector name and copy:

  • interaction dist to detector (mm): Distance from each neutron interaction or origin to the detector.
  • interaction dist to detector (mm) weighted by Hstar: Distance from each neutron interaction or origin to the detector, weighted by Hstar.
  • interaction dist to detector (mm) vs weight: Distance from each neutron interaction or origin to the detector versus weight when detector is reached.

The rest of histograms are for each detector copy and score type. The variables are calculated when the geantino or neutron/gamma reaches the detector, so they have the name At detector, and they also have the score type. There is a set for the geantinos and another set for the neutron/gamma particles. They are the following:

  • log10(energy) (MeV): log10 of track kinetic energy times the weight.
  • log10(energy) no weighted (MeV): log10 of track kinetic energy.
  • log10(energy) weighted by Hstar: log10 of track kinetic energy times the weight and Hstar value.
  • log10(weight): log10 of track weight.
  • log10(weight) weighted by Hstar: log10 of track weight times the weight and Hstar value.
  • log10(energy) vs log10(weight): log10 of track kinetic energy versus track weight.

9.10.3   Variance reduction techniques

There are three variance reduction techniques available, that will provide a better efficiency, i.e. less CPU time for the same error.

9.10.3.1   Kill contributions of low weight

The first technique consist on not taking into account the contributions of low weight. It works checking the weight of the geantinos at each step of their propagation (each time they change volume) and playing Russian roulette with them if their way is smaller than the value given by the parameter

/gamos/setParam GmPDS:UseMinimumGeantinoWeight MIN_WEIGHT

which by default takes a value of 1.E-30. The Russian roulette factor is given by the parameter

/gamos/setParam GmPDS:MinimumGeantinoWeightRRFactor VALUE

which by default takes a value of 100. To activate this technique you have to set the parameter

/gamos/setParam GmPDS:UseVRMinimumGeantinoWeight 1

You may use the control histograms, specially the histogram log10 of geantino kinetic energy versus geantino weight to determine which are the minimum weight and Russian roulette factor that best match your application.

9.10.3.2   Kill contributions too far from detector

A more efficient technique consists on not taking into account the contributions of interactions or sources that are too far from the detector, without the need to send a geantino to calculate the probability to reach it. But this technique assumes that the probability is smaller if the point is farther from the detector, what is not always the case, as it may be that a point closer to the detector traverses denser materials. Indeed the contributions are not eliminated, but, to avoid biasing Russian roulette is played. The maximum distance is set with the parameter

/gamos/setParam GmPDS:MaximumDistance DISTANCE

which by default takes a value of 1000 mm. The Russian roulette factor is given by the parameter

/gamos/setParam GmPDS:MaximumDistanceRRFactor VALUE

which by default takes a value of 100. To activate this technique you have to set the parameter

/gamos/setParam GmPDS:UseVRMaximumDistance 1

You may use the control histograms, specially the histogram Neutron: PD N: interaction dist to detector (mm) vs weight to determine which are the maximum distance and Russian roulette factor that best match your application.

9.10.3.3   Split particles of big weight

The idea of this technique is to sample more those cases which bigger probabilities. It may be very useful also to eliminate the seldom occurring big weight contributions that will make your error big. The user defines a minimum and maximum weight for the particles to be split:

/gamos/setParam GmPDS:MinimumWeightForSplitting MIN_VALUE

which by default takes a value of 1.E-30

/gamos/setParam GmPDS:MaximumWeightForSplitting MAX_VALUE

which by default takes a value of 1.E-5; and the maximum splitting number:

/gamos/setParam GmPDS:MaximumSplitting MAX_SPLIT

which by default takes a value of 100.

Each time a geantino reaches the detector, its weight will be checked. If it is between the minimum and maximum, the splitting number will be calculated making a logarithmic interpolation using 0 at the minimum weight and the maximum splitting number at the maximum weight. If this number, nSplit, is bigger than one, nSplit-1 particles will be created at the position, direction, energy and time of the neutron or gamma that created the geantino.

To activate this technique you have to set the parameter

/gamos/setParam GmPDS:UseVRSplitting 1

You may use the control histograms, specially the histogram SCORE_NAME At detector : geantino : log10(energy) vs log10(weight) to determine which are the minimum and maximum weight and splitting number that best match your application.

9.10.4   Sum scoring results and plot them

If you want to sum the scoring results from several jobs, for example, because you have split a long job in several ones, there is an utility to do it. You first have to build a file with the list of files you want to sum, one line per file. Then you have to run the utility sumScores. This utility will scan all files and will looks for all scorers with the same name (same MultiFunctionalDet: and PrimitiveScorer:). It merges this scorers calculating the average of each score. It also propagates the error and gives it relative to the value (divided by it); and it also gets the error as the deviation of the values, i.e. as the