Normalisation Factor¶
The MULASSIS normalisation factor, which is set using the
/analysis/normalisecommand, is the scale factor that converts the simulated incident particle current density to the “real world” current density such that the analysis outputs can be provided in terms of real world values.The normalisation factor to be used in MULASSIS is composed of two components:
a geometric factor, converting incident fluence to a particle current density through a surface and
the number of particles in the energy range under consideration.
The normalisation factor is the product of the two values.
Note
In this section, we use the term “particle current density” somewhat loosely. Strictly, this quantity refers to the number of particles flowing through a surface per unit area per unit time, whereas here we use it to refer to the time-integrated quantity.
Furthermore, current density is conventionally a measure of the net flow across a surface, so that equal numbers of particles per unit area crossing a surface in opposite directions results in zero current. Instead, in the context of this discussion on MULASSIS normalisation, the two flows are additive.
Omnidirectional fluence to current density:¶
Consider a spatially uniform fluence \(\phi(T, \theta, \omega)\) of radiation in space, differential in particle kinetic energy T, obliquity \(\theta\) and azimuth \(\omega\) (e.g. \(\phi\) in units of MeV-1 cm-2 sr-1). Assume this fluence to be isotropic, i.e. independent of \(\theta\) and \(\omega\), so that \(\phi(T, \theta, \omega) = \phi(T)\). The total \(4\pi\)-fluence is simply:
(1)¶\[\phi_{4\pi}(T) = \int_0^{2\pi} { d\omega \int_{-1}^1 {\phi(T)d(cos\theta) = 4\pi\phi(T) } }\]The number of particles crossing a surface boundary per unit area per unit solid angle, at an angle \(\theta\) with respect to the normal to the boundary is \(\phi(T) cos(\theta)\), hence the use of the term “cosine-law” current density or planar fluence in describing an isotropic fluence incident on a plane area. Defining \(u = cos \theta\), the total current density crossing the boundary is then:
(2)¶\[j(T) = \int_0^{2\pi} { d\omega \int_0^1 {\phi(T)u du = \pi\phi(T) } }\]The total \(4\pi\)-fluence and the total current density have the same units (e.g. MeV-1cm-2) but from Eqs (1) and (2) are related according to the expression:
(3)¶\[\phi_{4\pi}(T) = 4 j(T)\]In general, the outputs of Monte-Carlo calculations are inherently normalized to unit incident current density, i.e. \(j(T)\), as in the case of MULASSIS where the incident particles are simulated from \(\theta=0\) to \(\pi/2\).
Fluence in energy range of the simulation:¶
The fluence of particles in the energy range under consideration in the simulation is given by integrating the differential particle spectrum over the energy limits of the simulation, or:
(4)¶\[\Phi ] _{T min}^{T max} = \int_{T min}^{T max} {\phi(T) dT}\]Where \(\Phi\) is the fluence of particles in the range Tmin to Tmax.
In the case where an integral spectrum is provided, then it is a simple matter of subtracting the integral fluence at the maximum energy in the simulation from the fluence at the minimum energy in the simulation.
Normalisation factor:¶
The normalization factor, or current density in the energy range, for a differential flux/fluence spectrum is then a combination of equation 3 and equation 4, or:
(5)¶\[n = J ] _{T min}^{T max} = = \int_{T min}^{T max} {j(T) dT} = = \int_{Tmin}^{Tmax} {\frac{\phi(T)}{4} dT} = \frac{1}{4} \Phi ]_{Tmin}^{Tmax}\]
Application to various spectral forms:¶
Linear in energy:¶
Power law:¶
Exponential:¶
SPENVIS environment spectra:¶
The integral particle spectra, i.e. the fluence from \(T \to \infty\), are provided, so it is a simple matter of subtracting the fluence at the maximum energy in the spectrum from the fluence of the minimum energy in the spectrum, and dividing by 4, e.g.
(12)¶\[n=\frac{1}{4} \left ( \Phi ]_{Tmin}^{\infty} - \Phi ]_{Tmax}^{\infty} \right )\]
Slab geometry considerations¶
The fluence output provided by a slab geometry simulation is an area output, i.e. it will have the units #/cm-2.
The dose outputs are adjusted for the area of the slab layers and are reported as rads. This is consistent with SHIELDOSE [R16].
Spherical shell considerations¶
The geometric considerations of a spherical shell normalisation factor are already accounted for in the MULASSIS code and the above relationships are directly applicable with no modification.
Source particle generation from \(0 \to \theta_{max}\)¶
The above relationships hold when the source particles are generated in a zenith range from 0 to \(\pi/2\). When the solid angle of the source particles is limited to a range 0 to \(\theta_{max}\), then the integration over the angle gives an additional factor of \(sin^2 \theta_{max}\):
(13)¶\[N_{\theta_{max}} = n_\frac{\pi}{2} sin^2\theta_{max}\]
Examples:¶
In the following examples, a slab shielding geometry has been used with 1 mm of Aluminium shielding followed by 100 \(\mu m\) of Silicon. The incident particles are assumed to be isotropic, i.e. a cosine law angular distribution has been used and 100,000 incident protons (primaries) generated in the simulation.
Power law spectrum example¶
From equation 8, the definition of the differential power law spectrum, where A=1.09E09 p+/cm2/MeV, B=-3.0 and we are interested in particles in the range from 10 MeV to 250 MeV, it follows from equation 9 that:
(14)¶\[ \begin{align}\begin{aligned}n=\frac{1}{4} \left [ \frac{A}{B+1}T^{B+1} \right ]_{Tmin}^{Tmax}\\n=\frac{1}{4} \left [ \frac{1.09E09}{-3.+1}T^{-3.0+1} \right ]_{10 MeV}^{250 MeV}\\n=\frac{1}{4} \frac{1.0E09}{-2} \left [ 250^{-2} - 10^{-2} \right ]\\n=1.36E06\ cm^{-2}\end{aligned}\end{align} \]MULASSIS macro to run the power law simulation.¶/geometry/layer/delete 0 /geometry/layer/shape slab /geometry/layer/add 0 Aluminium 7 1.000E+00 mm /geometry/layer/add 1 Silicon 1 1.000E+02 um /analysis/normalise 1.360E+06 cm2 /analysis/fluence/particle/add proton /analysis/file PowerLawExample /analysis/fluence/unit cm2 /analysis/fluence/type OMNI /analysis/fluence/add 2 /analysis/fluence/energy/default /analysis/fluence/angle/default /analysis/dose/add 2 /analysis/dose/unit rad /geometry/update /phys/scenario em /gps/particle proton /gps/ene/type Pow /gps/ene/alpha -3.000E+00 /gps/ene/min 1.000E+01 MeV /gps/ene/max 2.500E+02 MeV /gps/ang/type cos /gps/ang/mintheta 0.000E+00 deg /gps/ang/maxtheta 9.000E+01 deg /control/execute display5.mac /event/printModulo 10000 /run/beamOn 100000
Exponential in energy spectrum example¶
We define differential energy spectrum which is an exponential in energy, as in equation 10, where A=1.99E09 p+/cm2/MeV and B=-17.4 MeV. If we are interested in particles in the range from 10 MeV to 250 MeV, it follows from equation 11 that:
(15)¶\[ \begin{align}\begin{aligned}n=\frac{1}{4}\left[ABe^\frac{T}{B}\right]_{T_{min}}^{T_{max}}\\n=\frac{1}{4}\left[1.99E09\times(-17.4)e^\frac{T}{-17.4}\right]_{10\ MeV}^{250\ MeV}\\n=\frac{1}{4}\times\left(-3.46E10\right)\times\left(e^{-\frac{250}{17.4}}-e^{-\frac{10}{17.4}}\right)\\n=4.87E09 {cm}^{-2}\end{aligned}\end{align} \]MULASSIS macro to run the exponential in energy simulation.¶# some comments /geometry/layer/delete 0 /geometry/layer/shape slab /geometry/layer/add 0 Aluminium 7 1.000E+00 mm /geometry/layer/add 1 Silicon 1 1.000E+02 um /phys/scenario em /analysis/fluence/particle/add proton /analysis/file ExponentialInEnergy /analysis/dose/add 2 /analysis/dose/unit rad /geometry/update /gps/particle proton # Power law in Energy for Integral flux spectrum # from 2003/10/26 event # F = A'*exp(E/B) # A' = 3.46E10 # B = 17.4 # df/dE = (A'/B)*exp(E/B) # A = (A'/B) = 1.99E09 /gps/ene/type Exp /analysis/normalise 4.87E09 cm2 /gps/ene/ezero 17.4 /gps/ene/min 1.000E+01 MeV /gps/ene/max 2.500E+02 MeV /gps/ang/type cos /gps/ang/mintheta 0.000E+00 deg /gps/ang/maxtheta 9.000E+01 deg /control/execute display5.mac /event/printModulo 10000 /run/beamOn 100000
SPENVIS [R17] environment spectrum example¶
Having run the SPENVIS models for a GTO-type orbit (300 km by 36000 km, 0 deg inclination and 10 year mission) to produce an integral trapped proton spectrum, as seen in the table below, the normalisation factor is calculated from equation 12 as:
(16)¶\[ \begin{align}\begin{aligned}n=\frac{1}{4}\left(\left.\mathrm{\Phi}\right]_{T_{min}}^\infty-\left.\mathrm{\Phi}\right]_{T_{max}}^\infty\right)\\n=\frac{1}{4}\left(\left.\mathrm{\Phi}\right]_{0.1\ MeV}^\infty-\left.\mathrm{\Phi}\right]_{400\ MeV}^\infty\right)\\n=\frac{1}{4}\left(1.9986E16-5.1447E9\right)\end{aligned}\end{align} \](17)¶\[n=4.9965E15\ {cm}^{-2}\]The MULASSIS-calculated TID in the 100 mu m planar of silicon layer shielded by 1 mm aluminium is 444 krad(Si) (+/-54 krad(Si)), while SHIELDOSE [R16] calculates this to be 482 krad(Si).
SPENVIS Trapped Proton Integral Spectrum¶ Energy
Total Mission Average Flux
Total Mission Fluence
[MeV]
[/cm2/s]
[/cm2]
0.10
6.3376E+07
1.9986E+16
0.15
4.9998E+07
1.5767E+16
0.20
4.0210E+07
1.2681E+16
0.30
2.8164E+07
8.8819E+15
0.40
2.0229E+07
6.3794E+15
0.50
1.5121E+07
4.7685E+15
0.60
1.1489E+07
3.6232E+15
0.70
8.9398E+06
2.8192E+15
1.00
4.5197E+06
1.4253E+15
1.50
1.7644E+06
5.5643E+14
2.00
8.1013E+05
2.5548E+14
3.00
3.2925E+05
1.0383E+14
4.00
1.6395E+05
5.1704E+13
5.00
1.0148E+05
3.2002E+13
6.00
6.5922E+04
2.0789E+13
7.00
4.6104E+04
1.4539E+13
10.00
1.8483E+04
5.8287E+12
15.00
5.8347E+03
1.8400E+12
20.00
2.8677E+03
9.0434E+11
30.00
1.4041E+03
4.4280E+11
40.00
1.0496E+03
3.3101E+11
50.00
8.0483E+02
2.5381E+11
60.00
6.8183E+02
2.1502E+11
70.00
5.9238E+02
1.8681E+11
100.00
4.0035E+02
1.2626E+11
150.00
2.2042E+02
6.9512E+10
200.00
1.2262E+02
3.8668E+10
300.00
4.4557E+01
1.4051E+10
400.00
1.6314E+01
5.1447E+09
Note
In the following MULASSIS macros, the point-wise spectrum is the differential spectrum calculated by SPENVIS. Alternatively, the user can use the SPENVIS integral spectrum from the above table, and include the macro command
/gps/ene/diffspec Fto treat the input as an integral spectrum in energy.MULASSIS macros to run the SPENVIS GTO example.¶/geometry/layer/delete 0 /geometry/layer/shape slab /geometry/layer/add 0 Aluminium 7 1.000E+00 mm /geometry/layer/add 1 Silicon 1 1.000E+02 um /phys/scenario hadron-leem-ln /analysis/normalise 4.9965E15 cm2 /analysis/fluence/particle/add proton /analysis/file Mulassis_Example1 /analysis/fluence/unit cm2 /analysis/fluence/type OMNI /analysis/dose/add 2 /analysis/dose/unit rad /geometry/update /gps/particle proton /gps/ene/type Arb /gps/ene/min 1.000E-01 MeV /gps/ene/max 4.000E+02 MeV /gps/hist/type arb /gps/hist/point 1.000000E-01 3.034500E+08 /gps/hist/point 1.500000E-01 2.316600E+08 /gps/hist/point 2.000000E-01 1.706600E+08 /gps/hist/point 3.000000E-01 9.990400E+07 /gps/hist/point 4.000000E-01 6.521600E+07 /gps/hist/point 5.000000E-01 4.369900E+07 /gps/hist/point 6.000000E-01 3.090600E+07 /gps/hist/point 7.000000E-01 2.280400E+07 /gps/hist/point 1.000000E+00 1.127500E+07 /gps/hist/point 1.500000E+00 3.709500E+06 /gps/hist/point 2.000000E+00 1.432700E+06 /gps/hist/point 3.000000E+00 3.230900E+05 /gps/hist/point 4.000000E+00 1.138900E+05 /gps/hist/point 5.000000E+00 4.901600E+04 /gps/hist/point 6.000000E+00 2.768600E+04 /gps/hist/point 7.000000E+00 1.716600E+04 /gps/hist/point 1.000000E+01 6.703000E+03 /gps/hist/point 1.500000E+01 1.561500E+03 /gps/hist/point 2.000000E+01 4.443900E+02 /gps/hist/point 3.000000E+01 9.090100E+01 /gps/hist/point 4.000000E+01 2.996400E+01 /gps/hist/point 5.000000E+01 1.839000E+01 /gps/hist/point 6.000000E+01 1.062200E+01 /gps/hist/point 7.000000E+01 8.309300E+00 /gps/hist/point 1.000000E+02 5.350000E+00 /gps/hist/point 1.500000E+02 2.777400E+00 /gps/hist/point 2.000000E+02 1.564300E+00 /gps/hist/point 3.000000E+02 5.315200E-01 /gps/hist/point 4.000000E+02 3.334700E-02 /gps/hist/inter Lin /gps/ang/type cos /gps/ang/mintheta 0.000E+00 deg /gps/ang/maxtheta 9.000E+01 deg /event/printModulo 100000 /run/beamOn 1000000
SPENVIS trapped environment spectrum with angular biasing¶
Using the same particle spectrum from the previous example, we want to calculate the dose on a 100 \(\mu m\) silicon sphere inside of an aluminium spherical shell 1 mm thick. Because the target area is greatly reduced, most simulation particles will miss the target completely. By limiting the direction in which the particles are fired to, e.g. \(\theta\) = 0\(^\circ\) to 10\(^\circ\), instead of the default range of 0\(^\circ\) to 90\(^\circ\), the efficiency of the simulation can be greatly improved.
The additional components to the normalisation factor to be included are:
the angular difference, a factor of \(sin^2 \theta_{max}\), and
the additional \(4\pi\) of particle generation space (this is included automatically by MULASSIS for a spherical geometry).
This yields a normalisation factor for a MULASSIS spherical geometry of:
(18)¶\[ \begin{align}\begin{aligned}N= sin^2 (10 deg) n_{spenvis}\\N= 0.03015 \times 4.9965E15 cm^{-2}\\N=1.506E14 cm^{-2}\end{aligned}\end{align} \]And the previous macro file example becomes:
MULASSIS macros to run the SPENVIS angular biasing example.¶/geometry/layer/delete 0 /geometry/layer/shape sphere /geometry/layer/add 0 Aluminium 7 1.000E+00 mm /geometry/layer/add 1 Silicon 1 1.000E+02 um /phys/scenario hadron-leem-ln /analysis/normalise 1.506E+14 cm2 /analysis/fluence/particle/add proton /analysis/file Mulassis_Example2 /analysis/fluence/unit cm2 /analysis/fluence/type OMNI /analysis/dose/add 2 /analysis/dose/unit rad /geometry/update /gps/particle proton /gps/ene/type Arb /gps/hist/type arb /gps/ene/min 1.000E-01 MeV /gps/ene/max 4.000E+02 MeV /gps/hist/point 1.000000E-01 3.034500E+08 /gps/hist/point 1.500000E-01 2.316600E+08 /gps/hist/point 2.000000E-01 1.706600E+08 /gps/hist/point 3.000000E-01 9.990400E+07 /gps/hist/point 4.000000E-01 6.521600E+07 /gps/hist/point 5.000000E-01 4.369900E+07 /gps/hist/point 6.000000E-01 3.090600E+07 /gps/hist/point 7.000000E-01 2.280400E+07 /gps/hist/point 1.000000E+00 1.127500E+07 /gps/hist/point 1.500000E+00 3.709500E+06 /gps/hist/point 2.000000E+00 1.432700E+06 /gps/hist/point 3.000000E+00 3.230900E+05 /gps/hist/point 4.000000E+00 1.138900E+05 /gps/hist/point 5.000000E+00 4.901600E+04 /gps/hist/point 6.000000E+00 2.768600E+04 /gps/hist/point 7.000000E+00 1.716600E+04 /gps/hist/point 1.000000E+01 6.703000E+03 /gps/hist/point 1.500000E+01 1.561500E+03 /gps/hist/point 2.000000E+01 4.443900E+02 /gps/hist/point 3.000000E+01 9.090100E+01 /gps/hist/point 4.000000E+01 2.996400E+01 /gps/hist/point 5.000000E+01 1.839000E+01 /gps/hist/point 6.000000E+01 1.062200E+01 /gps/hist/point 7.000000E+01 8.309300E+00 /gps/hist/point 1.000000E+02 5.350000E+00 /gps/hist/point 1.500000E+02 2.777400E+00 /gps/hist/point 2.000000E+02 1.564300E+00 /gps/hist/point 3.000000E+02 5.315200E-01 /gps/hist/point 4.000000E+02 3.334700E-02 /gps/hist/inter Lin /gps/ang/type cos /gps/ang/mintheta 0.000E+00 deg /gps/ang/maxtheta 1.000E+01 deg /event/printModulo 100000 /run/beamOn 1000000The final dose in a spherical shell of silicon 100 mu m thick behind a 1 mm thick spherical shell of aluminium is 3.01+/-0.48 Mrad(Si), while SHIELDOSE [R16] calculates the value to be 3.06 Mrad(Si). Without the angular restriction, at least 33 times more primaries are required for a similar level of accuracy (~15%).
Note
In the above discussion, we refer to the time-integrated quantities particle fluence and dose. The example presented (SPENVIS [R17] environment spectrum example) analyzes the 10-year mission cumulative proton fluence for a GTO orbit and the resulting TID in silicon shielded by aluminium. The applied normalisation factor is the particle current density at the MULASSIS geometry surface corresponding to the 10 year mission, defined using the command:
/analysis/normalise 4.9965E15 cm2
However, we can also calculate shielded particle fluxes and dose rates by replacing the normalisation factor with the calculated number of particles incident upon the MULASSIS geometry surface per unit area per unit time. For the GTO orbit example, the mean particle current becomes:
The MULASSIS command to define the normalisation should not mention units of time:
/analysis/normalise 1.5844E+7 cm2
but the MULASSIS outputs can be interpreted by the user as average shielded particle fluxes and dose rates per second.
Mulassis/ml-v02-00/r342
