Variance reduction techniques

Introduction

Variance reduction are the techniques that serve to get the same precision of the estimates of a Monte Carlo simulation while reducing the CPU time, and without introducing a bias in the result. The usual criteria to measure the gain of a variance reduction is the efficiency, defined as the CPU time divided by the square of the error or the variable of interest:

\epsilon = \frac{T}{\sigma}

Importance sampling

Importance sampling is a variance reduction techniques that consists in estimating the properties of a particular distribution, while using a different sampling than the distribution of interest. It may help you saving CPU time if you choose to sample more times when the contribution to the distribution is bigger and viceversa. An example can be the case when you want to estimate the dose deposited by particles that reach a volume; in this case you may sample more particles when the distance to reach that volume is smaller. This can be done by multiplicating a particle N times (of course reducing its weight 1/N) when the particle comes closer to the volume, where N is inversely proportional to the distance to the volume.

In GAMOS you may do importance sampling using many different criteria: any of the GAMOS data can be used, as the importance sampling process uses a GAMOS distribution, taking the output value of the distribution as the splitting value that corresponds to each input value. If the splitting value (SPLIT_VALUE) is bigger than 1 it will duplicate the current particle SPLIT_VALUE-1 times (producing SPLIT_VALUE equal particles, with a weight reduced by 1/(SPLIT_VALUE). If the splitting value is smaller than 1, Russian roulette will be played with the particle with a surviving probability SPLIT_VALUE; if it survives its weight will be incremented by 1/(SPLIT_VALUE).

You may activate the importance sampling for a list of particles with the command:

/gamos/physics/VR/importanceSampling NAME DISTRIBUTION_NAME PARTICLE_NAME_1 PARTICLE_NAME_2 …

where the NAME you give to the importance sampling will be later used to set the parameters.

As the number of particle duplication will be set by the classifier index, you may want to set a different index than the default of the classifier. See section on Classifiers to learn how to do this.

Several options can be controlled with GAMOS parameters. First you may decide to deactivate the splitting, leaving only Russian roulette, or viceversa, with the parameters:

/gamos/setParam IMPORTANCE_SAMPLING_NAME:ApplySplitting 0/1

/gamos/setParam IMPORTANCE_SAMPLING_NAME:ApplyRussianRoulette 0/1

If you do not want that a particle already split is split again, you may control how many times a particle is split with the parameter:

/gamos/setParam IMPORTANCE_SAMPLING_NAME:MaxSplitTimes VALUE

You may also set that the same particle is not split several times, by setting the parameter:

/gamos/setParam IMPORTANCE_SAMPLING_NAME:SplitAtSeveralSteps

Finally you may add one of more filters, so that particles have to pass them before being split:

/gamos/setParam IMPORTANCE_SAMPLING_NAME:Filters FILTER_1 FILTER_2 …

Geometrical biasing

Geometrical biasing is a technique where the splitting is done as some function of the geometrical position. Many types of geometrical biasing can be done with GAMOS using the importance sampling by assigning different distributions with different data. For example you may use a different weight for each physical volume, logical volume or touchable or you can divide your space in bins along X, Y or Z or combinations of these variables. You may also consider the special distribution GmGeometricalBiasingDistribution (see chapter on Distributions), that calculates the weight as the division of the weights of the entering and exiting volumes.

Biasing operations

GAMOS provides commands to use the Geant4 biasing techniques in an easy way. Moreover several extra biasing techniques are implemented in GAMOS non existing in Geant4. We describe here those techniques and how to use them.

Following the Geant4 notation, the first thing you have to create a biasing operator, with the command:

/gamos/physics/biasing/createOperator OPERATOR_NAME OPERATOR_TYPE

where OPERATOR_NAME is the name you choose for your operator, and OPERATOR_TYPE has to be one of the following ones:

  • ChangeCrossSection:
  • ForceCollision
  • BremsSplitting
  • DirBremsSplitting
  • EWBS

These biasing operators will only valid for the particle processes selected by the user. Therefore you have to select the particles and processes. This can be done with three different commands:

Select a certain process of a certain particle:

/gamos/physics/biasing/addParticleProcesses2Oper PARTICLE_NAME PROCESS_NAME OPERATOR_NAME

Select all the processes of a certain particle:

/gamos/physics/biasing/addProcesses2Oper PROCESS_NAME OPERATOR_NAME

Select a certain process of all particles :

/gamos/physics/biasing/addParticle2Oper PARTICLE_NAME OPERATOR_NAME

Finally after the /run/initialize command (when the geometry is already built) you have to attach the biasing operator to one or several logical volumes, so that biasing will only occur when a particle is being tracked in one of those volumes:

/gamos/physics/biasing/associateOper2LogVol LOGICAL_VOLUMES OPERATOR_NAME

Example:

/gamos/physics/biasing/createOperator boper ForceCollision
/gamos/physics/biasing/addParticleProcesses2Oper gamma compt boper
/run/initialize
/gamos/physics/biasing/associateOper2LogVol crystal boper

Cross section biasing

The cross section of any process of any particle can be changed by using the biasing operator ChangeCrossSection. The cross section will be multiplied by a factor that has to be chosen by the user with the GAMOS parameter

/gamos/setParam OPERATOR_NAME:XSfactor VALUE

An example of its use in GAMOS is the following:

*/gamos/setParam myOperator:XSFactor 10.

*/gamos/physics/biasing/createOperator myOperator ChangeCrossSection

*/gamos/physics/biasing/addParticles2Oper gamma myOperator

*/run/initialize

*/gamos/physics/biasing/associateOper2LogVol target myOperator

Force collision

You can force a process to happen in a volume by using the biasing operator ForceCollision. We copy here the description of this biasing from the Geant4 user’s manual.

G4BOptrForceCollision is a biasing operator that implements a force collision scheme quite close to the one provided by MCNP. It handles the scheme though the following sequence:

  1. The operator starts by using a G4BOptnCloning cloning operation, making a copy of the primary entering the volume. The primary is given a zero weight.

  2. The primary is then transported through to the volume, without interactions. This is done with the operator requesting forced free flight G4BOptnForceFreeFlight operations to all physics processes. The weight is zero to prevent the primary to contribute to scores. This flight purpose is to accumulate the probability to fly through the volume without interaction. When the primary reaches the volume boundary, the first free flight operation restores the primary weight to its initial weight and all operations multiply this weight by their weight for non-interaction flight. The operator then abandons here the primary track, letting it back to normal tracking.

  3. The copy of the primary track starts and the track is forced to interact in the volume, using the G4BOptnForceCommonTruncatedExp operation, itself using the total cross-section to compute the forced interaction law (exponential law limited to path length in the volume). One of the physics processes is randomly selected (on the basis of cross-section values) for the interaction.

  4. Other processes are receiving a forced free flight operation, from the operator.

  5. The copy of the primary is transported up to its interaction point. With these operations configured, the G4BiasingProcessInterface instances have all needed information to automatically compute the weight of the primary track and of its interaction products. As this operation starts on the volume boundary, a single force interaction occurs: if the track survives the interaction (e.g Compton process), as it moved apart the boundary, the operator does not consider it further.

    An example of its use in GAMOS is the following:

    */gamos/physics/biasing/createOperator myOperator ForceCollision

    */gamos/physics/biasing/addParticles2Oper gamma myOperator

    */run/initialize

    */gamos/physics/biasing/associateOper2LogVol target myOperator

Uniform bremsstrahlung splitting

Bremsstrahlung splitting is a technique aimed to increase the statistics of gammas produced by bremsstrahlung in the region of interest, while reducing the time spent in tracking the electrons. It is a technique very frequently used to simulate radiotherapy linear accelerators. To use this technique you have to select the biasing operator BremsSplitting. The behaviour of this operator is controlled with a few parameters:

/gamos/setParam OPERATOR_NAME:NSplit VALUE: number of gammas that will be produced at each bremsstrahlung interaction. The weight of each gamma is set to 1./VALUE.

/gamos/setParam OPERATOR_NAME:BiasOnlyOnce 0/1: if set to 1 (default is 1) when the particle that suffers a bremsstrahlung collision is the descendant of a particle that already was split (e.g. electron(bremsstrahlung) -> gamma(compton) -> electron) is is not split.

/gamos/setParam OPERATOR_NAME:BiasPrimaryOnly 0/1: if set to 1 (default is 0) only primary particles will be split.

An example of its use in GAMOS is the following:

*/gamos/setParam boper:NSplit 10

*/gamos/setParam boper:BiasPrimaryOnly 0

*/gamos/setParam boper:BiasOnlyOnce 0

*/gamos/physics/biasing/createOperator boper BremsSplitting

*/gamos/physics/biasing/addParticleProcesses2Oper e- eBrem boper

*/run/initialize

/gamos/physics/biasing/associateOper2LogVol *target boper

Directional bremsstrahlung splitting

This technique first defines a rectangle in the XY plane situated at a selected Z coordinate. Then for each of the gammas generated by bremsstrahlung splitting it checks if its directions points towards this rectangle and if not plays Russian roulette with it. The result is that you end with many gammas in the region of interest and only a few gammas out of it. To use this technique you have to select the biasing operator DirBremsSplitting. The behaviour of this operator is controlled with a few parameters:

/gamos/setParam OPERATOR_NAME:NSplit VALUE: number of gammas that will be produced at each bremsstrahlung interaction. The weight of each gamma is set to 1./VALUE.

/gamos/setParam OPERATOR_NAME:XDim VALUE: dimension of the rectangle of intereset in the X coordinate (rectangel extends from -X to +X)

/gamos/setParam OPERATOR_NAME:YDim VALUE: dimension of the rectangle of intereset in the Y coordinate (rectangel extends from -Y to +Y)

/gamos/setParam OPERATOR_NAME:ZPos VALUE: Z position of the XY plane where the rectangle of interest is defined

/gamos/setParam OPERATOR_NAME:NKill VALUE: number used to play Russian roulette (probability is 1./VALUE). By default it is the same as NSplit

/gamos/setParam OPERATOR_NAME:BiasOnlyOnce 0/1: if set to 1 (default is 1) when the particle that suffers a bremsstrahlung collision is the descendant of a particle that already was split (e.g. electron(bremsstrahlung) -> gamma(compton) -> electron) is is not split.

/gamos/setParam OPERATOR_NAME:BiasPrimaryOnly 0/1: if set to 1 (default is 0) only primary particles will be split.

An example of its use in GAMOS is the following:

*/gamos/setParam boper:NSplit 10

*/gamos/setParam boper:BiasPrimaryOnly 0

*/gamos/setParam boper:BiasOnlyOnce 0

*/gamos/setParam boper:XDim 400

*/gamos/setParam boper:YDim 355

*/gamos/setParam boper:ZPos -1000

*/gamos/physics/biasing/createOperator boper BremsSplitting

*/gamos/physics/biasing/addParticleProcesses2Oper e- eBrem boper

*/run/initialize

/gamos/physics/biasing/associateOper2LogVol *target boper

Equal weight bremsstrahlung splitting

When bremsstrahlung splitting is used it often happens that some of the particles that contribute to the magnitude of interest (e.g. dose in a radiotherapy planning) have been split and have an small weight, while others have survived Russian roulette and their weight has been set back to 1. The fact of having particles of different weights contributing to a magnitude increases substantially the total error and therefore should be avoided. This technique is based in the one described at [biblio.dbs], and it handles the weights of all secondary particles created so that those that reach the region of interest have the same weight.

To use this technique you have to select the biasing operator EWPS. The behaviour of this operator is controlled by the same parameters than the direction bremsstrahlung technique described above.

To assure that all particle weights are treated correctly, it is necessary to include the following processes:

/gamos/physics/biasing/EWPS/addParticleProcesses2Oper e- eBrem OPERATOR_NAME

/gamos/physics/biasing/EWPS/addParticleProcesses2Oper e+ annihil OPERATOR_NAME

/gamos/physics/biasing/EWPS/addParticleProcesses2Oper gamma compt OPERATOR_NAME

/gamos/physics/biasing/EWPS/addParticleProcesses2Oper gamma conv OPERATOR_NAME

/gamos/physics/biasing/EWPS/addParticleProcesses2Oper gamma phot OPERATOR_NAME

An example of its use in GAMOS is the following:

*/gamos/setParam boper:NSplit 10

*/gamos/setParam boper:BiasPrimaryOnly 0

*/gamos/setParam boper:BiasOnlyOnce 0

*/gamos/setParam boper:XDim 400

*/gamos/setParam boper:YDim 355

*/gamos/physics/biasing/createOperator boper EWBS

*/gamos/physics/biasing/EWBS/addParticleProcesses2Oper e- eBrem boper

*/gamos/physics/biasing/EWBS/addParticleProcesses2Oper e+ annihil boper

*/gamos/physics/biasing/EWBS/addParticleProcesses2Oper gamma compt boper

*/gamos/physics/biasing/EWBS/addParticleProcesses2Oper gamma conv boper

*/gamos/physics/biasing/EWBS/addParticleProcesses2Oper gamma phot boper

*/run/initialize

/gamos/physics/biasing/associateOper2LogVol *target boper