Particle beams flattening simulation
Monte Carlo simulation of particle beam flattening
Using a dual-scattering-foil technique
E.Detsi
Summary
Using the Monte Carlo Method, the dual scattering technique to flatten a charged
particle beam carrying an energy of 160 MeV is simulated. An optimized shape of the
second scatterer is found. A flat particle field having a diameter of 8 cm is obtained.
The accuracy and the efficiency of the model are found to be approximately ± 1.35%
and 23% respectively.
1
Particle beams flattening simulation
Introduction
Ion beam produced by a particle accelerator can be used in radiation therapy to destroy
tumors cells by radiation damage. Unfortunately, in addition to destroying tumor cells,
healthy tissue cells are damaged as well [6].
In general, the particle beam delivered by a particle accelerator is very narrow and nonuniform. Such a narrow and non-uniform shape is difficult to handle during the
radiation therapy. This increases the probability of damaging healthy tissue cells [2].
The aim of this research simulation is to investigate how a broad and uniform particle
beam can be obtained from a narrow and non-uniform one. Nowadays, one of the most
used techniques [2] to obtain a flat particle beam is the dual scattering of the incoming
narrow particle beam into two foils placed at a certain distance apart from each other.
This scattering technique will be investigated by simulation.
Theory
Monte Carlo method [7]
The Monte Carlo method is a statistical simulation method that makes use of random
numbers to study the behavior of a mathematical or physical system. Because of the
randomness of the generated numbers, the Monte Carlo method appears to be a
nondeterministic method. This means that repeated simulations with the same sequence
of parameters never predict exactly the same behavior of the system of interest.
However in the case of a physical system, repeating simulations with larger sequences
of random numbers yields a stable behavior of the studied system. Similarly, in the case
of a mathematical system, repeating simulations with lager sequences of random
numbers gives results which converge to a specific value. It is proved that this specific
value to which the Monte Carlo simulations results converge is nothing else than the
solution of the considered mathematical problem.
Scattering theory
The notion of scattering is used to describe the collision of particles. Scattering can be
classically analyzed if particles are considered as ‘corpuscles’. However, based on the
particle-wave dual nature of particles [5], an alternative approach to study the scattering
involves the quantum mechanic, which widely describes the scattering in three
dimensions. In this section, only the classical approach namely the Coulomb scattering
is discussed.
The Rutherford scattering [5]
Quick description
It is hard to talk about scattering without mentioning the pioneers in this sector, Ernest
Rutherford and his two colleagues, Hans Geiger and Ernest Marsden who carried out
different scattering experiments and developed one of the first scattering theories at the
University of Manchester at the beginning of the twentieth century.
In their experiments, it was surprisingly observed that when alpha particles from a
radioactive source strike a thin gold foil, they are found to be deflected from their initial
trajectory. Some of those alpha particles were even found to be back-scattered [5] as
shown in figure1 below.
The deflection of the incoming alpha particles from their initial trajectory was explained
by the presence in the gold foil of positively charged nuclei generating a central
Coulomb repulsive force.
2
Particle beams flattening simulation
Figure 1 Coulomb scattering
Application of the Rutherford scattering
From the above described Rutherford experiment, the fraction of the scattered particles
with respect to the incident particles is found to be[10] :
Rscattered N avogadro d"#
(0)
=
Rincident
10 !3 Amass
With:
Rincident is the rate of incident particles per unit area
Rscattered is the rate of scattered particle per unit area
d is the target thickness
! is the target density
Amass is the mass number of the nuclei of the target atom
! is the scattering cross section
From relation (0) , it is clear that for a given is the scattering cross section ! , the
number of ‘collisions’ is proportional to the thickness (d) of the target. In other words,
the thicker the scattering foil, the weaker the energy of the outgoing particles and the
smaller the average scattered angles. That principle is used in this research project to
find out the optimal thickness of the scattering foils that generates a flat outgoing beam.
Coulomb scattering
The Coulomb scattering occurs when a charge particle traversing in a medium is
deflected by elastic collision in the Coulomb field of the nuclei. The Rutherford
scattering defined in section 3.2.1. is an example of Coulomb scattering.
Scattering angle distribution
The angular distribution of the scattered particle at large angle deflection is given by the
Rutherford scattering formula [5]:
N i ntZ 2 e 4
(1)
N (! ) =
2 2 2
4 !
(8"# 0 ) r E sin ( )
2
N (! ) : Number of alpha particles per unit area that reach the screen at a scattering angle !
N i : Total number of alpha particles that reach the screen
3
Particle beams flattening simulation
n : Number of atoms per unit volume in the foil
Z : Atomic number of the foil atoms
r : Distance of the screen from the foil
E : Kinetic energy of alpha particles
t : Foil thickness
Charged particles traversing a thick foil are submitted to repeated elastic Coulomb
scattering. The cumulative effect of these multiple scatterings results in an average
small angle deflection from the primitive direction as shown in figure 2 below.
Figure 2 Average angle after multiple scattering in a thick target
For small those angles deflection, the multiple scattering angular distribution is good
approximated by a Gaussian distribution [7]:
1
$" 2
P(" ) =
exp(
)
(2)
2#! 2
2! 2
Where the standard deviation ! is given by the Highland formula [8]:
" =
180 Mev
z d [1 + 0.038 ln(d )]
! .cp
(3)
d : Thickness of the foil
p : Momentum of the incident particle
! .c : Velocity of the incident particle
z : Charge number of the incident particle
Application of the scattering angle distribution
In this experiment, the relation " # c i d ! is applied to determine the scattering angle
distribution for both lead and tungsten. As shown on the curves below, when the
dimension (d) is measured in mm and the deviation ( ! ) is measured in mrad, the
values of the constant ! are 0.5591 and 0.5596 for lead and tungsten respectively.
Similarly, the values of the constants ci are 15.759 and 20.52 for lead and wolfram
respectively.
4
Particle beams flattening simulation
scatter angle (mrad)
Scatter
anlge
vs.
thickness
(Pb)
Scatter
30
30
25
25
20
anlge
vs.
thickness
(W)
20
15
15
10
10
5
0,5591
y = 15,759x
5
0,5596
y = 20,52x
0
0
0,5
1
1,5
2
thickness (m m )
0
0,05
0,55
1,05
thickness
1,55
(mm)
Figure 3 Scattering angle as function of the foil thickness for both lead and tungsten [9]
Modeling
Overview
Figure 3 below shows an overview of the model which is used for the simulation. This
model is simply formed by two scattering foils and a detector placed at a certain
distance apart from each other.
Figure 4 Scattering system formed by two foils and a screen [3]
The incident particle beam is scattered in foil 1 and foil 2 respectively. The intensity of
the scattered beam is recorded by the detector placed beyond the second scattering foil.
Incident particle beam
The incoming particle beam represented in blue in figure 3 an carries energy of 160
MeV.
First scattering foil
The first scattering foil is an isotropic material with a uniform thickness. For this
simulation, the chosen material is lead (Pb) and the uniform thickness is denoted by d1 .
Second scattering foil
The second scattering foil is also an isotropic material but has a non-uniform thickness.
In the present case, the maximum thickness is denoted by d 2 , at the center of the foil.
The uniform thickness is denoted by d flat , at the extremities of the foil. The chosen
material is tungsten (W). (see figure 5).
5
Particle beams flattening simulation
rbol
krom
d flat
Figure 5 Cross section of the second foil
d2
This second scattering foil is placed at distances l1 from the first one and l 2 from the
screen.
Motivation for the shape of the second scatterer
Before the optimal shape shown in figure 5a above is obtained, different models for the
shape of the second scatterer are tried without satisfactory results. Below follows a brief
presentation of these models:
A foil with a rectangular cross section
Figure 6. Left: rectangular shape of the second scatterer.
Right: correspond scattered field
When the second scattering foil has a flat shape with a rectangular cross section, no flat
field is obtained very thin scatterer. However, increasing the thickness the this
rectangular scatterer yieds a flat particle field. As one can see on figure 6 above (right),
the first curve is obtained for a thickness of 0.5 mm. It still has a Gaussian. The second
curve is obtained for a thickness of 5 mm. One sees this second particles field becomes
qualitatively flat. Unfortunately this is not a good result because the shape can not be
optimized: the thicker the scatterer, the flatter the field.
A foil with a half ciruclar cross section
Next a foil with a half circular cross section as shown in figure 7 is tested. The
motivation for this shape with a non constant thickness is directly related to the
scattering theory as presented in section 3.2.1. (“The thicker the scattering foil, the
6
Particle beams flattening simulation
weaker the energy of the outgoing particles and the smaller the average scattered
angles”). The peak of the Gaussian incident beam goes trough the centre of the
scattering foil. Since the centre of the foil is thicker than the extremities, the peak is
expected to be more flatten.
Figure 7. Left: foil with a half circular cross section.
Right: correspondent fields for as function of the curvature of the scattering foil.
By varying the curvature of the scatterer, the obtained particle field is about to be flat,
but it is not really flat. Depending on the curvature of the scatterer, this field presents a
hole or a bole at the center.
Correction of the hole.
At this stadium to get a flat field from the previous results shown in figure 7, one may
chose to correct the bole or the central hole. When chosing to correct the hole, the
extremities of the field are then flatten again to yield the same high as the central hole.
To do this, the shape of the scattering foil is adjusted as shown in figure 8 below.
Figure 8. Evolution in the shape of the scattering foil during the scattering shape
investigation and correspondent flat field.
7
Particle beams flattening simulation
Final design of the scatterer shape:
The previous shape presented in figure 8 above yields a flat field as expected, but this
shape has inconvenient, the variation of the curvature affect considerably all the other
parameters of that foil. To solve this problem, the alternative shape presented in figure 9
below is used. With this new shape, one can independently vary the curvature of the
scattering foil. The result is a larger flat field.
Figure 9. Final choice for the shape of the scatterer and correspondent flat field
Mathematical description of the model
Cylindrical coordinates
The above described model is mathematically schematized in figure 6 below. The
reference frame is chosen so that the direction of the incoming particle beam coincides
with the z-direction. The origin of the reference frame is taken at the impact of the
incoming particle with the first scattering foil.
8
Particle beams flattening simulation
eˆx
eˆz
r
p2
eˆy
r2
r
p1
" 2 ,! 2
r1
" 1 , !1
Particle
Beam
l2
l1
Foil 1
Foil 2
Figure 10. Description of the path follows by a particle from the source to the screen in
Cylindrical coordinates
The parameters used in this model are defined as follow:
d1 : Thickness first foil
d flat : Thickness at the extremities of the second foil
d 2 : Maximal thickness of the second foil (see figure 5 )
l1 : Distance between the source and the first foil
l 2 : Distance between the second foil and the detector
r
p1 : Position of the particle on foil 2 with respect to the origin of the reference frame
r
p2 : Position of the particle on the detector with respect to the impact point on foil 2.
r
r
!1 : Angle between p1 and e z
r
r
! 2 : Angle between p2 and e z
r
r
!1 : Rotation of p1 about e z
r
r
! 2 : Rotation of p2 about e z
In cylindrical coordinates the position of the particle after scattering in the first foil is
r
given by the norm of p1 and the angles ("1 , !1 ) as shown on figure 6 above. Similarly,
the position of the particle after scattering in the second foil is described by the norm of
r
p2 and the angles (" 2 , ! 2 ) .
Cartesian coordinates
For simplicity reasons the model is discussed in Cartesian coordinates rather than in
cylindrical coordinates. In Cartesian coordinates, the position of the particle after
9
Particle beams flattening simulation
r
scattering in foil 1 is completely defined by the norm of p1 and the angles (!1x , !1 y ) as
shown in figure 7 below. Similarly, the position of the particle after scattering on the
r
second foil is described by the norm of p2 and the angles (! 2 x ,! 2 y ) .
Figure 11. Position and angular distribution in Cartesian coordinates after scattering
in the first foil [3]
Scattering in the first foil
r
The components of p1 with respect to the origin of the reference frame are given by:
& x1 # &l1 tan (1x #
r $ ! $
!
p1 ' $ y1 ! = $l1 tan (1 y !
(4)
!
$%l1 !" $%l1
"
From relation (4) the radius on foil 2 of the disc where the particle falls is deduced:
r1 = x12 + y12 = (l1 tan !1x ) 2 + (l1 tan !1 y ) 2
(5)
Scattering in the second foil
r
Similarly, the components of p2 with respect to the origin of the reference frame are
given by:
& x2 # & x1 + l2 tan ( 2 x # &l1 tan (1x + l2 tan ( 2 x #
r
$
! $
!
p2 ' $$ y 2 !! = $ y1 + l2 tan ( 2 y ! = $l1 tan (1 y + l2 tan ( 2 y !
(6)
! $l + l
!
$%l1 + l2 !" $%l1 + l2
" %1 2
"
From relation (6) the radius on the screen of the disc where the particle falls is deduced:
10
Particle beams flattening simulation
r2 = x22 + y 22 = (l1 tan !1x + l 2 tan ! 2 x ) 2 + (l1 tan !1 y + l 2 tan ! 2 y ) 2
(7)
Scattering angles
As discussed in the theory (section 3.2.3), if c1 defines the characteristic constant of
lead (foil 1) then the outgoing angle after scattering in that foil has an average standard
deviation in the x-and y-direction given by:
! 1x = c1 d10.5591
(8)
! 1 y = c1 d10.5591
(9)
Similarly, if c 2 defines the characteristic constant of tungsten (foil 2) then the outgoing
angle after scattering in the second foil has an average standard deviation in the x-and
y-direction given by
! 2 x = c 2 d z0.5596
(10)
! 2 y = c 2 d z0.5596
(11)
Where d z denotes the non-constant thickness of the second foil. d z varies along the zdirection as follow:
d z = d flat
dz = d2 1!
(12)
At the extremities of the foil ( r1 ! rbol )
r12
2
rparameter
(13)
At the center of the foil ( r1 < rbol )
Implementation
Principal Algorithm
The principal algorithm for the simulation is now described. The full program can be
found in the appendix 1.
function hxy = scatterfoil(N,L1,hmax,krom,rbol,d2,dflat)
This principal algorithm has been implemented in a matlab, as a function which
describes the movement of each particle from the source where it is emitted to the
screen. This function is declared by scatterfoil( ).The inputs and output arguments are
defined as follow:
Input arguments
Depending on the parameters of interest, the function scatterfoil( ) expects seven input
arguments:
-
N: specifies the number of particle in the incident beam
L1: also called l1 in figure 6. With this parameter the distance between the
source and the first foil can be varied.
11
Particle beams flattening simulation
-
hmax: this is a measure for the number of bins in the output parameter hxy.
krom: this parameter is used to vary the curvature of the second foil.(see fig. 5)
rbol: a part of the N incident particles reach the second foil in a disc with radius
rbol. The parameter rbol is simply used to localize these particles. (see fig. 5)
dflat: another part of the N incident particles do not fall in the disc with radius
rbol. The correspondent foil thickness for this group of particle is dflat.
Output arguments
Also depending on the parameters of interest, many outputs argument can be associated
to the principal function scatterfoil( ). In the present case, one is mainly interested in the
position and density of the particles on the screen after two scatterings. One output
argument is sufficient to collect the data containing this information. This argument is
denoted by hxy, which is a two dimensional square histogram (20 x 20).
Generation of N incident particles
N=N;
n=4;
z=randn(n, N);
% number of particles --> N columns
% number of angle vectors --> 4 rows
% Normally distributed random numbers
Angular distribution after scattering in foil 1
Here the relation " # c i d ! described in section 3.2.3 is used to determine the angular
distribution. In this relation, c i ! lead , d ! d 1 and " ! 0.5591 .
sd1=lead*d1^0.5591;
sd1x=sd1/sqrt(2);
sd1y=sd1x;
% standard angular deviation after scattering in folie-1
% x-component of the standard angular deviation
% y-component of the standard angular deviation
theta1x=z(1,:)*sd1x;
theta1y=z(2,:)*sd1y;
% x-component of the angular distribution between
% y-component of the angular distribution between
Position of the scattered particles on the second foil
x1=l1*tan(theta1x); % x-coordinate on foil 2
y1=l1*tan(theta1y); % y-coordinate on foil 2
r1=sqrt(x1.^2+y1.^2); % radius on foil 2
Angular distribution after scattering in foil 2
Here again the relation " # c i d ! described in section 3.2.6 is used to determine the
angular distribution. In this relation tungsten, d ! d z and " ! 0.5596 .
12
Particle beams flattening simulation
dz ;
sd2=tungsten*dz.^0.5596;
sd2x=sd2/sqrt(2);
sd2y=sd2x;
% see formula (12) and (13) also figure 5 for dz
% standard angular deviation after scattering in folie-2
% x-component of the standard angular deviation
% y-component of the standard angular deviation
theta2x=z(3,:).*sd2x; % x-component of the angular distribution between
theta2y=z(4,:).*sd2y; % y-component of the angular distribution between
Position of the scattered particles on the screen
x2=x1+l2*tan(theta1x+theta2x);
y2=y1+l2*tan(theta1y+theta2y);
r2=sqrt(x2.^2+y2.^2);
% x-coordinate on the screen
% y-coordinate on the screen
% radius on the screen
Storing data in the two dimensional histogram hxy
nx=1+(nxybins/2)+floor((nxybins/2.0)*x2/xmax); % position of particles in de x-direction
ny=1+(nxybins/2)+floor((nxybins/2.0)*y2/ymax); % position of particles in de y-direction
for i=1:N
hxy(nx(i),ny(i))=hxy(nx(i),ny(i))+1; % all the data are stored in the matrix hxy
end
Algorithm for the flatness
This algorithm is briefly described. The full program can be found in the appendix 2.
function [deviation,binAverage] = flatness(hxy,radius)
This algorithm has been implemented in a matlab, as a function which analyzes the data
stored in hxy
The function is declared by flatness( ).The inputs and output arguments are defined as
follow:
Input arguments
Two inputs arguments are required namely the matrix hxy containing the data and the
radius of the disc where one wants to check the flatness.
13
Particle beams flattening simulation
Outputs arguments
The principal output argument is deviation which is a measure for the flatness of the
field in a given radius. Other secondary outputs arguments can be obtained:
- binAverage which denotes the average number of particle in each cell.
- k which denotes the number of cells
- sum which denotes the total number of particles in a chosen radius. With the output
argument sum, the efficiency of the system can be deduced since the quotient (
sum
*100 ) gives the percentage of incident particles in the flat field.
N
nbins=20;
rbin=radius/binLength;
% length in ‘bins’ of the radius ‘rbin’
for i=1:nbins
for j=1:nbins
if sqrt((i-ibins+0.5)^2+(j-ibins+0.5)^2)<rbin; % bin in the desired radius
k=k+1;
sum=sum+hxy(i,j);
sumSquare= sumSquare+hxy(i,j)^2;
end
end
end
k; % count, number of cells
sum; % sum of particle in the disc with radius rbin
sumSquare; % sum^2 of particle in the disc with radius rbin
binAverage=sum/k;
flatness=k/(k-1)*(((sumSquare/k)-(sum/k)^2)/binAverage);
deviation=flatness;
Finding the best parameters
Two matlab functions denoted by myfucntion( ) and investigate( ) are implemented and
used to find out the parameters for which flatness is minimal
Example
If one is interested in finding out the optimal distance l1 between to source and the first
foil and also the curvature parameter krom for the second foil, one can operate as
follow:
First in the function scatterfoil() all the ‘known’ parameters are replaced by their values
behalf of the parameters which one wants to investigate
function f = myfunction(parameter)
L1=parameter(1);
krom=parameter(2);
14
Particle beams flattening simulation
hxy=scatterfoil(1000000, L1, 0.09, krom, 0.0056, 0.000712, 0.00030 );
deviation=flatness(hxy, 0.04);
Secondly using the function investigate( ) one can find out the best value of L1 and
krom for which deviation is minimal.
function value=investigate(parameter)
for i = 1:1:10
for j=1:1:10
parameter(1)=i*0.0325;
parameter(2)=j*0.00065;
% L1 varies between 3.25 cm and 32.5 cm
% krom varies between 0.65 mn and 6.5 mn
value(i,j)=myfunction(parameter);
end
end
It is important to mention that finding out those parameters in two dimensions as
describes above is mainly possible when one already has an idea about the value of
these parameters. Otherwise it is not trivial to found the minima for deviation with a
given parameter. The main reason for which the minimum is found using a two
dimensional raster instead of making use of the matlab optimization Toolbox (fminunc
or fminsearch ) is that, since the function deviation is non-deterministic, its derivative
always yields different results.
Critical values
Once the optimal parameters are found, the correspondent field can only be considered
as flat if the related deviation is smaller than the theoretical critical value given by:
function criticalValue (numberCells)
zp95=1.65;
zp99=2.33;
% in 95% confidence interval
% in 99% confidence interval
k=numberCells ;
n=k-1;
% number of cells
% degree of freedom
criticalVaule_95percent = (0.5*(zp95 + sqrt(2*n-1))^2)/n
criticalValue_99percent = (0.5*(zp99 + sqrt(2*n-1))^2)/n
15
Particle beams flattening simulation
Results and discussions
Set of optimal parameters
De optimal parameters found for the function scatterfoil(N,L1,hmax,krom,rbol,d2,dflat)
are given by:
hxy=scatterfoil(1000000, 0.325, 0.09, 0.005624, 0.0056, 0.000712, 0.00030)
This means that:
-
L1 is found to be 32.5 cm (distance between the source and the first foil)
hmax is found to be 18.0 cm (size of the output matrix hxy)
krom is found to be 5.624 mm (curvature of the second foil)
rbol is found to be 5.6 mm with radius rbol.
d2 is found to be 0.7 mm ( maximal thickness of the second foil)
dflat is found to be 0.30 mm (foil thickness at the extremities)
Figure 12. Scattering foil 2 with the correspondent parameters
16
Particle beams flattening simulation
Qualitative images
Figure 13. left-3D Incident beam, Right-3D outgoing beam
Figure 14 Outgoing beam in 2D
Scaling: from bins to centimeters
The data obtained from the simulation are stored in 400 cells formed by the square
matrix hxy. Since the size of hxy is found to be 18 cm, converting the dimension of the
18cm
= 0,9cm for the size of each cell.
cells from bins to centimeters gives
20
Let consider various concentric discs with r = 1,2,3,4,...cm centered in the square as
shown in figure 10 below. The results of the simulation is given and analyzed with
respect to each of those concentric discs.
17
Particle beams flattening simulation
Figure 15 Scaling: from bins to centimeters
Quantitative results
For each of the above mentioned radii, the function flatness( ) is ran with the set of
found parameters given in section 6.1. The correspondent deviation is determined.
The obtained value of deviation which is a measure for the flatness is compared
with the theoretical value obtained from the function criticalValue( ).The flatness
criterion is satisfied when the simulation value is smaller than the theoretical one.
The results are summarized in table 1 below.
18
Particle beams flattening simulation
Table 1 Summarized results
Number
of particles
500,000
1,000,000
1,500,000
2,000,000
4,000,000
Radius
(in cm)
Number
of cells
(k)
Degree of
freedom
(k-1)
Flatness values
from
Simulation
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
4
16
32
60
96
140
4
16
32
60
96
140
4
16
32
60
96
140
4
16
32
60
96
140
4
16
32
60
96
140
3
15
31
59
95
139
3
15
31
59
95
139
3
15
31
59
95
139
3
15
31
59
95
139
3
15
31
59
95
139
0.3438
1.0861
0.9534
1.2270
2.5172
24.0942
1.4160
1.5587
1.3092
1.1708
3.7139
47.3368
1.6294
0.8287
1.2592
1.1641
5.1122
67.7539
0.9298
1.2198
0.9569
1.1233
7.4448
92.0938
0.5144
0.7385
1.1717
1.1515
12.8406
179.1374
Critical values
from the theory
95%
99%
confidence
confidence
interval
interval
2.5169
3.4748
1.6498
1.9841
1.4435
1.6585
1.3171
1.4647
1.2478
1.3605
1.2038
1.2949
2.5169
3.4748
1.6498
1.9841
1.4435
1.6585
1.3171
1.4647
1.2478
1.3605
1.2038
1.2949
2.5169
3.4748
1.6498
1.9841
1.4435
1.6585
1.3171
1.4647
1.2478
1.3605
1.2038
1.2949
2.5169
3.4748
1.6498
1.9841
1.4435
1.6585
1.3171
1.4647
1.2478
1.3605
1.2038
1.2949
2.5169
3.4748
1.6498
1.9841
1.4435
1.6585
1.3171
1.4647
1.2478
1.3605
1.2038
1.2949
From table 1, one sees that the particle field is flat for radii up to 4 cm as shown in
green, since the flatness criteria are fulfilled for those radii.
For a radius of 5 cm, the simulation value of the flatness is slightly larger than the
theoretical one. That radius can be defined as a transition radius between the flat and the
non-flat field.
For radii beyond 5 cm, the simulation value is largely bigger than the theoretical one, as
shown in red in table 1.
Based on these results, one is motivated to conclude that the particles field is flat for
radii up to 4 cm. Unfortunately, this single sequence of simulations is not statistically
representative. Therefore the results can not yet be generalized. The reason why one
19
Particle beams flattening simulation
must be careful with these results is related to the nature of the generated events, which
are non-deterministic. This means that repeated simulations with the same set of
parameters given in section 6.1 will give different results. In other words simulating a
flat field for radii up to 4 centimeters with this set of parameters does not necessary
imply that all other simulations with the same set of parameters will also yield a flat
field. However it is know from the Monte Carlo method as indicated in section 3.1 that
repeated simulations with larger sequences of random numbers yields a stable behavior
of the studied system. In the present case, the obtained results can only be generalized if
the Monte Carlo method confirms these results.
Confirmation with the Monte Carlo method
When applying the Monte Carlo method, each simulation is repeated 20 times with the
same set of parameters given in section 6.1. To run the simulation 20 times, a matlab
function is made. This function used the primary functions scatteroil( ) and flatness( ).
It is chosen to run the simulation 20 times and not more because of the speed of the
computer. Running a simulation 20 times with a large number of particles say 2 till 4
millions as shown in table 2 slows the speed of the computer down.
Table 2 Average critical value over 20 simulations
Average flatness values
Number of
particles
500,000
1,000,000
1,500,000
Radius
(in cm)
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
1
2
1 20
=
! [critical _ values]i
20 i =1
1.0229
1.0987
1.0895
1.1122
2.4997
23.1448
1.1696
1.1523
1.0594
1.1321
4.1880
45.3443
0.7422
1.3713
1.2130
1.2705
6.5370
72.5292
1.0182
1.2727
Critical values
from the theory
95%
confidence
interval
2.5169
1.6498
1.4435
1.3171
1.2478
1.2038
2.5169
1.6498
1.4435
1.3171
1.2478
1.2038
2.5169
1.6498
1.4435
1.3171
1.2478
1.2038
2.5169
1.6498
99%
confidence
interval
3.4748
1.9841
1.6585
1.4647
1.3605
1.2949
3.4748
1.9841
1.6585
1.4647
1.3605
1.2949
3.4748
1.9841
1.6585
1.4647
1.3605
1.2949
3.4748
1.9841
20
Particle beams flattening simulation
2,000,000
4,000,000
3
4
5
6
1
2
3
4
5
6
1.1539
1.1701
8.004
90.353
1.0499
1.2031
1.1652
1.1825
13.047
172.4713
1.4435
1.3171
1.2478
1.2038
2.5169
1.6498
1.4435
1.3171
1.2478
1.2038
1.6585
1.4647
1.3605
1.2949
3.4748
1.9841
1.6585
1.4647
1.3605
1.2949
Here again, a flat filed is obtained for radii up to 4 cm. These last results can now be
generalized since they are not the results of a single event, but the average results for
many events.
Efficiency of the system
The efficiency of the system is approximately 23% for a flat field with diameter 8 cm.
A few efficiency results are summarized on the table 3
Table 3
Number of incidents
particles
Radius of
the flat field
Number of particles
in the flat field
Efficiency
500,000
3 cm
4 cm
60911
114424
12.18 %
22,89 %
1,000,000
3 cm
4 cm
121198
227344
12.12 %
22.73 %
1,500,000
3 cm
4 cm
184110
342568
12.27 %
22.83%
21
Particle beams flattening simulation
Error analysis
Figure 16 Blue line: simulation results
Red stars: plane
When Fitting the obtained flat field with a plane, the average value of the high of that
flat field is calculated. An horizontal plane is assigned to this high. The lateral error
corresponds to the vertical deviation from that plane. See figure 11. The blue line
represents the simulation results and the red crosses represent the plane.
A more accurate scale (zoom) shows clearly the lateral fluctuation on see in Figure 12
below, which is a zoom of the previous cross section.
Figure 17 lateral error less than 1.35 %
After many simulations, the average value of the lateral error is found to less than
1.35 %. The way one single error is found is as follow:
If one denotes by H the average horizontal high of the plane formed by the flat field
and by ! the maximum deviation from H, then the error is given by:
!
Error = 100.
(14)
H
22
Particle beams flattening simulation
The single error for the cross section given in figure 12, this error is found to be:
90
Error = 100.
= 1.184%
7600
As mentioned above, the final error ( ± 1.35% ) is the average value of many single
errors. This is because the error is also non-deterministic.
Conclusion
The flattening of a charge particle beam using the dual scattering technique is
successfully simulated. The simulation results show that this technique really works. A
narrow incident particle beam carrying an energy of 160 MeV is successfully flattened
after double scattering in a lead foil and in a tugsten foil respectively. The obtained
particles field is flat in a diameter of 8 cm. These results are in agreement with other
experimental and theoretical results, namely these published [2] by the group of Erik
Grussel, from the university of Uppsala, Sweden. The proposed model for the
simulation is optimized. More attention is paid to the design of the shape of the second
scatterer. The accuracy of the simulation results ( ± 1.35% ) is in the same order as the
one proposed by Erik Grussel and al. ( ± 2.0% ). However, the efficiency of the system
(23%) is lower than the efficiency presented by Erik Grussel and al. (38%). A possible
explanation of this difference in efficiency can be related to the simplicity of model
used for the simulation. In the proposed model, only one interaction is taken in
consideration namely the Coulomb scattering. Furthermore, the scattered particles are
considered as moving in vacuum. Also the incoming particles beam is assimilated to a
punt source (without dimension). In reality, the physical situation is more complex. One
is satisfied with these results because if such results are obtained with a simplified
model, improvements in the results are then expected if one takes in account the other
neglected physical parameters.
Appendix
Appendix 1 Algorithm for the data matrix ‘hxy’
function hxy = scatterfoil(N,L1,hmax,krom,rbol,d2,dflat)
%--------------------------------------------------------------------------------------------------%
%
%
% DETERMINE THE POSITION AND DENSITY OF THE PARTICLE AFTER DUAL %
%
%
% SCATTERING AND STORE THESE INFORMSATION IN A SQUARE MATRIX HXY %
%
%
%= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = %
nxybins=20;
hxy(nxybins,nxybins)=0.0;
xmax=hmax;
ymax=hmax;
N=N;
n=4;
z=randn(n, N);
% number of bins in hxy --> 20 x 20 = 400
% initialization of the matrix hxy;
% half length of the matrix hxy in the x-direction
% half length of the matrix hxy in the y-direction
% number of particles --> N columns
% number of angle vectors --> 4 rows
% Normally distributed random numbers
theta1x(N)=0.0;
% initialization of the angle vector theta1x
23
Particle beams flattening simulation
theta1y(N)=0.0;
theta2x(N)=0.0;
theta2y(N)=0.0;
x1(N)=0.0;
y1(N)=0.0;
x2(N)=0.0;
y2(N)=0.0;
% initialization of the angle vector theta1y
% initialization of the angle vector theta2x
% initialization of the angle vector theta2y
% initialization of the position vector x1;
% initialization of the position vector y1;
% initialization of the position vector x2;
% initialization of the position vector y2;
lead=0.015759*1000^0.5591;
tungsten=0.02052*1000^0.5596;
d1=0.0012;
d2=d2;
dz(1:N)=dflat;
rbol=rbol;
d2par=krom;
l1=L1;
l2=3.0;
% characteristic constant --> lead
% characteristic constant --> wolfram
% thickness in meters --> foil 1
% maximal thickness --> foil 2
% minimal thickness --> foil 2
% bol radius --> foil 2
% curvature --> foil 2
% distance between the source and the first foil
% distance between the second foil and the detector in meters
sd1=lead*d1^0.5591;
sd1x=sd1/sqrt(2);
sd1y=sd1x;
% standard angular deviation after scattering in folie-1
% x-component of the standard angular deviation
% y-component of the standard angular deviation
theta1x=z(1,:)*sd1x; % x-component of the angular distribution between
theta1y=z(2,:)*sd1y; % y-component of the angular distribution between
x1=l1*tan(theta1x); % x-coordinate on foil 2
y1=l1*tan(theta1y); % y-coordinate on foil 2
r1=sqrt(x1.^2+y1.^2); % radius on foil 2
for i=1:N
if r1(i)<rbol
dz(i)=d2*(1-(r1(i)^2/d2par^2));
end
if dz(i)<0
dz(i)=0;
end
end
sd2=tungsten*dz.^0.5596;
sd2x=sd2/sqrt(2);
sd2y=sd2x;
% standard angular deviation after scattering in folie-2
% x-component of the standard angular deviation
% y-component of the standard angular deviation
theta2x=z(3,:).*sd2x;
theta2y=z(4,:).*sd2y;
% x-component of the angular distribution
% y-component of the angular distribution
x2=x1+l2*tan(theta1x+theta2x);
y2=y1+l2*tan(theta1y+theta2y);
r2=sqrt(x2.^2+y2.^2);
% x-coordinate on the screen
% y-coordinate on the screen
% radius on the screen
nx=1+(nxybins/2)+floor((nxybins/2.0)*x2/xmax); % position of particles in the direction
24
Particle beams flattening simulation
ny=1+(nxybins/2)+floor((nxybins/2.0)*y2/ymax); % position of particles in the y-direction
for i=1:N
if nx(i)>nxybins
nx(i)=nxybins;
end
if ny(i)>nxybins
ny(i)=nxybins;
end
end
for i=1:N
if nx(i)<1
nx(i)=1;
end
if ny(i)<1
ny(i)=1;
end
end
for i=1:N
hxy(nx(i),ny(i))=hxy(nx(i),ny(i))+1; % all the data are stored in the matrix hxy
end
end function
===========================================
Appendix 2 Algorithm for the flatness of the data stored in the matrix ‘hxy’
function [deviation,binAverage] = flatness(hxy,radius)
%----------------------------------------------------------------------------------------------------%
% DETERMINE THE FLATNESS OF THE PARTICLE FIELD FOR A GIVEN RADIUS
%
%
%
%
AND THE AVERAGE NUMBER OF PARTICLE IN EACH CELL
%
%= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = %
nbins=20;
ibins=nbins/2;
% number of bins in hxy --> 20 x 20 = 400
% centre of the matrix
hmax =0.09;
binLength=hmax/ibins;
% length in meters of ‘ibins’
% length of one bin := 0.09/10 = 9mm
rbin=radius/binLength;
% length in bins of the radius rbin
k=0;
sum=0;
sumSquare =0;
% count, number of cells
% sum of particle in the disc with radius rbin
% sum^2 of particle in the disc with radius rbin
25
Particle beams flattening simulation
for i=1:nbins
for j=1:nbins
if sqrt((i-ibins+0.5)^2+(j-ibins+0.5)^2)<rbin; % bin in the desired radius
k=k+1;
sum=sum+hxy(i,j);
sumSquare= sumSquare+hxy(i,j)^2;
end
end
end
k;
sum;
sumSquare;
binAverage=sum/k;
flatness=k/(k-1)*(((sumSquare/k)-(sum/k)^2)/binAverage);
deviation=flatness;
end function
===========================================
26
Particle beams flattening simulation
References
[1] Bron KVI
[2] Erik Grussel, Anders Montelius, Anders Brahme, Goran Rikner and Kellie
Russel, A general solution to charged particle beam flattening using an optimized
dual-scattering-foil technique, with application to proton therapy beams. Phys.
Med. Biol. 39 (1994) 2201-2216.
[3] Amles D.
[4] Richard L. Liboff, Introductory Quantum Mechanics, Fourth Edition
[5] Arthur Beiser, Concept of Modern Physics, Fifth Edition.
[6] rug/BMT
[7] wikipedi
[8] V. Highland, NIM 129, 497 (1975)
[9] R. Anne, J. Herault, R. Bimbot, H. Gauvin, G. Basten, F. Hubert, Multiple
angular scattering of heavy ions At intermediate energies.
[10] http://hyperphysics
27