Skip to content

Re write of variance reduction section #45

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions getting-started/install.rst
Original file line number Diff line number Diff line change
@@ -1,2 +1,5 @@
Installation
============

Installation instructions for the latest TOPAS version are provided at the top of the `Code Repository <https://sites.google.com/a/topasmc.org/home/code-repository-authorized-users-only>`_ page of the TOPAS web site.
This page is accessible only to licensed TOPAS users.
23 changes: 23 additions & 0 deletions parameters/variance/directional_filter.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
.. _vr_directional_filter:

Directional Filter
==================

Often one is only interested in having low noise data in some particular ROI, however a lot of simulation time can be wasted simulating particles that are very unlikely to reach this ROI. Therefore, in certain situations one may wish to simply kill such particles. To do this, use the settings below (analogous to :ref:`vr_geometrical_splitting`)::


# Directional particle killing
# ------------------------------------------------------------
s:Vr/DirectionalFilter/Type = "DirectionalRussianRoulette"
s:Vr/DirectionalFilter/ReferenceComponent = "tungsten_Target"
dv:Vr/DirectionalFilter/ForRegion/VarianceReduction/DirectionalSplitLimits = 2 -1.0 -1.0 m
sv:vr/DirectionalFilter/ForRegion/VarianceReduction/processesNamed = 2 "eBrem" "compt"
dv:Vr/DirectionalFilter/ForRegion/VarianceReduction/DirectionalSplitRadius = 2 5.0 5.0 cm

Note that this technique is similar to UseDirectionalSplitting technqiue described in :ref:`_vr_secondary_biasing`. There are two differences:

1. The directional filter technique does not need to be applied in conjunction with secondary biasing.
2. The directional filter technique kills particles, whereas the directional splitting technique applies russian roulette. (we know that the way these techniques are currently named is confusing; we are working on it!)

.. image:: secondary_biasing.png
Figure: Biasing particle of secondary photons after a bremsstrahlung process without (left) and with (right) a directional filter.
29 changes: 0 additions & 29 deletions parameters/variance/geometrical_particle_split.rst

This file was deleted.

17 changes: 0 additions & 17 deletions parameters/variance/importance_sampling.rst

This file was deleted.

5 changes: 2 additions & 3 deletions parameters/variance/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,9 @@ Variance Reduction
:maxdepth: 2

intro
geometrical_particle_split
importance_sampling
weight_window
topas_split_geometry_vr
specific_particles
secondary_biasing
directional_filter
forced_interaction
crosssection_enhancement
29 changes: 0 additions & 29 deletions parameters/variance/intro.rst
Original file line number Diff line number Diff line change
Expand Up @@ -20,32 +20,3 @@ All variance reduction techniques in topas can be switched on/off with a master
b:Vr/UseVarianceReduction = "true"

Specific variance reduction techniques are described next.

Specify the Split Geometry
~~~~~~~~~~~~~~~~~~~~~~~~~~

The geometry for variance reduction must be in a :ref:`parallel world <geometry_parallel>`. The type of component can be any standard solid (:ref:`geometry_dividable` or :ref:`geometry_generic`). The geometry must consist of a geometry component with a set of geometry sub-components as daughters. The sub-components must be located in such a way that the boundaries coincide. The split process or Russian roulette will occur at these boundaries. In the next figure a simple scheme is shown.

.. image:: split_geometry.png

:ref:`time_feature` can be used to move or rotate the component or sub-components. But there is a restriction: the implementation of VR does not allow you to change the dimensions of the component and sub-components.

To set the geometry for VR::

s:Vr/ParticleSplit/Component = "MyComponent"
sv:Vr/ParticleSplit/SubComponents = n "MySubComp_1" ... "MySubComp_n"



Define the Splitting Technique
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

There are three variance reduction techniques available:

* :ref:`GeometricalParticleSplit <vr_geometrical_splitting>`
* :ref:`ImportanceSampling <vr_importance_sampling>`
* :ref:`WeightWindow <vr_weight_window>`

To chose a technique, use for example::

s:Vr/ParticleSplit/Type = "GeometricalParticleSplit"
28 changes: 11 additions & 17 deletions parameters/variance/secondary_biasing.rst
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
.. _vr_secondary_biasing:

Secondary Biasing
-----------------
=================

The split of secondary particles created after an electromagnetic interaction is also supported. A common example is the split of secondary photons created in the bremsstrahlung process for conventional radiotherapy simulations.

Expand All @@ -11,26 +13,18 @@ The region ``MyRegion`` is automatically created if it does not exist. The next

s:Vr/ParticleSplit/Type = "SecondaryBiasing"

Then three vectors must be defined. One with the name of the electromagnetic processes, one with the split number for each process and one with the maximum energies for each processes. The biased particles with energies larger than these values are subject to Russian roulette::
Then three vectors must be defined. One with the name of the electromagnetic processes, one with the split number for each process and one with the maximum energies for each processes. Particles with energy less than this are split according to SplitNumber and their statistical weight reduced by a factor of SplitNumber. Particles with energies larger than these values are subject to Russian roulette::

sv:Vr/ParticleSplit/ForRegion/MyRegion/ProcessesNamed = 2 "eBrem" "compt"
uv:Vr/ParticleSplit/ForRegion/MyRegion/SplitNumber = 2 100 10
dv:Vr/ParticleSplit/ForRegion/MyRegion/MaximumEnergy = 2 6.0 0.511 MeV

If suitable, further CPU time can be saved with a directional Russian roulette for secondary particles created with split (analogous to :ref:`vr_geometrical_splitting`). For that, a new VRT must be set::

s:Vr/DirectionalFilter/Type = "DirectionalRussianRoulette"

Later, a reference component must to be chosen::

s:Vr/DirectionalFilter/ReferenceComponent = "Target"

And the directional filter is applied::

dv:Vr/DirectionalFilter/ForRegion/MyRegion/DirectionalSplitLimits = 2 1.0 1.0 m
dv:Vr/DirectionalFilter/ForRegion/MyRegion/DirectionalSplitRadius = 2 5.0 5.0 cm
sv:vr/DirectionalFilter/ForRegion/MyRegion/processesNamed = 2 "eBrem" "compt"
If suitable, further CPU time can be saved with a directional Russian roulette for secondary particles created::

.. image:: secondary_biasing.png
b:Vr/ParticleSplit/UseDirectionalSplitting = "True"
d:Vr/ParticleSplit/TransX = 0.0 mm
d:Vr/ParticleSplit/TransY = 0.0 mm
d:Vr/ParticleSplit/TransZ = 0.0 mm
d:Vr/ParticleSplit/Rmax = 21 mm

Figure: Biasing particle of secondary photons after a bremsstrahlung process. On the left, no directional Russian roulette is applied. On the right, a directional Russian roulette is applied.
This creates an ROI at the centre of the world. If the created secondaries are projected to land in this ROI, they are kept. If not, they have a 1/N chance of surviving, where N is the same as SplitNumber. Particles which do survive this roulette process have their statistical weight increased by a factor of N.
113 changes: 113 additions & 0 deletions parameters/variance/topas_split_geometry_vr.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
.. _vr_topas_split_geometry:

Topas split geometry variance reduction techniques
==================================================

One means to reduce variance is to split a particle and reduce it's statistical weight whenever it crosses some boundary. There are a number of ways to implement such an approach in topas.

Specify the Split Geometry
--------------------------

The geometry for variance reduction must be in a :ref:`parallel world <geometry_parallel>`. The type of component can be any standard solid (:ref:`geometry_dividable` or :ref:`geometry_generic`). The geometry must consist of a geometry component with a set of geometry sub-components as daughters. The sub-components must be located in such a way that the boundaries coincide. The split process or Russian roulette will occur at these boundaries. In the next figure a simple scheme is shown.

.. image:: split_geometry.png

:ref:`time_feature` can be used to move or rotate the component or sub-components. But there is a restriction: the implementation of VR does not allow you to change the dimensions of the component and sub-components.

To set the geometry for VR::

s:Vr/ParticleSplit/Component = "MyComponent"
sv:Vr/ParticleSplit/SubComponents = n "MySubComp_1" ... "MySubComp_n"



Define the Splitting Technique
------------------------------

There are three variance reduction techniques available:

* :ref:`GeometricalParticleSplit <vr_geometrical_splitting>`
* :ref:`ImportanceSampling <vr_importance_sampling>`
* :ref:`WeightWindow <vr_weight_window>`

To chose a technique, use for example::

s:Vr/ParticleSplit/Type = "GeometricalParticleSplit"

Each technique described in detail below.

.. _vr_geometrical_splitting:

Geometrical Particle Splitting
------------------------------

TOPAS variance reduction is further described in:

Ramos-Mendez et al, “Geometrical splitting technique to improve the computational efficiency in Monte Carlo calculations for proton therapy,” Med. Phys. 40, 041718 (2013)

This technique was designed for heavy charged particles. In this implementation, you must specify whether the beam entering into the sub-component has cylindrical symmetry or not. This is because the particles may or may not be randomly redistributed around the ``SplitAxis``.

The Russian roulette is applied in a particular direction. That is, at the split plane and prior to being split, the particle is subject to the Russian roulette if its direction does not point towards a Region of Interest (ROI). Then the radius of the ROI and its position on the ``SplitAxis`` must to be defined too. Further, the Russian roulette can be turned on/off at specific surfaces between sub-components.

.. code::

s:Vr/ParticleSplit/Type = "GeometricalParticleSplit"
s:Vr/ParticleSplit/SplitAxis = "zaxis"
d:Vr/ParticleSplit/RussianRoulette/ROIRadius = 7.8 cm
d:Vr/ParticleSplit/RussianRoulette/ROITrans = 10 cm
bv:Vr/ParticleSplit/RussianRoulette = 2 "false" "true"

To set whether the region at each sub-component is symmetric or not and to define the corresponding split number::

bv:Vr/ParticleSplit/Symmetric = 2 "false" "true"
uv:Vr/ParticleSplit/SplitNumber = 2 8 8

.. image:: geometrical_splitting.png

In addition for this technique, geometrical Russian roulette will be played if a particle leaves the component or the world in a scheme similar to the :ref:`vr_importance_sampling` technique.

.. _vr_importance_sampling:

Importance Sampling
-------------------

In this technique, an importance value is assigned to each sub-component. If a particle is transported into a sub-component with a higher importance, then the particle is split. If it is transported into a sub-component with a lower importance, then Russian roulette is played. By default an importance value of 1 is automatically assigned to the parent component and to the world.

.. warning::

It is desirable for the thickness of the sub-components to be similar to the mean free path of the physical process to be biased.

The sub-component importance values are defined by hand. For example, to split the particles by a factor of 2 between subsequent sub-components, one must to define::

s:Vr/ParticleSplit/Type = "ImportanceSampling"
uv:Vr/ParticleSplit/ImportanceValues = 5 1 2 4 8 16

.. image:: importance_sampling.png

.. _vr_weight_window:

Weight Window
-------------

In this technique, the split process or Russian roulette will be applied depending on the statistical weight of the particle. Every time that a particle crosses from a sub-component to the next one, the statistical weight is evaluated.

* Particles with weights greater than a lower bound and smaller than an upper bound will be tracked normally.
* Particles with weights smaller than a lower bound will be subject to Russian roulette. If it survives, the particle is tracked normally but its weight is increased.
* Particles with weights greater than an upper bound will be split, and the new particles will be assigned lower weights.

The split number and Russian roulette criteria are internally calculated from an energy map, a weight map, an upper limit factor and a survival factor. In simple geometries the maps can be input by hand.

The user must provide a double vector with upper energy bounds and a unitless vector with lower weight bounds for every sub-component: ``WeightMap`` and ``EnergyMap``. The inverse of a parameter named ``MaximumSplitNumber`` (100 by default) is used to specify the minimum survival probability to be used in Russian roulette. The parameter ``PlaceOfAction`` states whether the split process (or Russian roulette) will occur at the sub-component boundaries, at physics interactions or at both.

The follow configuration is equivalent to importance sampling with importance generator of 2::

s:Vr/ParticleSplit/Type = "WeightWindow"
uv:Vr/ParticleSplit/WeightMap = 4 1. 1. 0.125 0.0615
dv:Vr/ParticleSplit/EnergyMap = 4 1. 1. 1. 1. GeV
u:Vr/ParticleSplit/UpperLimitFactor = 1
u:Vr/ParticleSplit/SurvivalFactor = 1
i:Vr/ParticleSplit/MaximumSplitNumber = 100
s:Vr/ParticleSplit/PlaceOfAction = "onBoundary"
#Others options of PlaceOfAction: "OnCollision" and "OnBoundaryAndCollision"

.. image:: weight_window.png
27 changes: 0 additions & 27 deletions parameters/variance/weight_window.rst

This file was deleted.