Transcript a/2 a/2 - Bose Institute

Monte Carlo Simulation Techniques
s
Pravata K Mohanty
Tata Institute of Fundamental Research, Mumbai
Winter School on Astroparticle Physics, Bose Institute, Darjeeling,
14 - 22 December 2009
Deterministic Process
For a given input, the outcome can be exactly determined.
s
Example:
X1
X2
X3
f(x) =
x2
Y1
Y2
Y3
Stochastic Process
Outcome can’t be exactly determined
Example: Interaction of primary cosmic rays in the atmosphere
s and development of EAS
E1
E2
E3
E1 = E2 = E3 = E
Atmosphere
N1
N2
N3
N1 = N2 = N3
Stochastic Process
P+P
s
0


-
+
+

-

Q. what are the life time of the muons?
Probability and Random numbers
Number of muons survive after time t
N (t) = N(0) e –t/
s
Number of decays after time t
Ndecay(t)= N(0) - N(t)= N(0) – N(0) e
-t/
Decay probability P = N decay(t) / N(0)
= 1 – e -t/
Or t = - ln(1-P),
0 < P <1
Replace P with R, where you call R as a random
number and 0<R<1
t = - ln(1-R)
Muon decay
s
Monte Carlo simulation
•
s
The average behaviour of the process
obtained from the measurements i.e. the
current knowledge about the process
Ex.
•
N (t) = N(0) e –t/
Convert it to a probability distribution
and use random numbers for
probability to generate the variates.
Random Numbers

True random numbers:
Obtained from natural processes

Pseudo random numbers:
generated by computers using some algorithm.
Example: Linear Congruential Generator
X n+1 = a X n + c mod m
m is the period, X0 is the SEED
m and c should be relatively prime
In C++, m = 232, a=214013, c = 2531011
How to generate Random Numbers
In C or C++, you can generate random numbers like this
for (int i=0; i<10000; i++) {
r = rand();
}
//Here rand() is the random
//number generator
Value of  using random numbers
Area of Circle/ Area of Square =  (a/2)2 /a2 = /4
Manual Method:
1. Randomly throw pebbles inside the square
2. Count the number of pebbles inside the circle
3. Take ratio of the number of pebbles inside
circle to the total.
-a/2
a/2
(0,0)
Using Computer
1. Generate points with x and y coordinates
uniformly inside the square of side a
x = -a/2 + a/2*R1 R1 and R2 are Random Numbers
y = -a/2 + a/2**R2
2. count the number of points inside the circle
r = (x2 + y2 ) < a/2
-a/2
a/2
History of Monte Carlo Simulation
The name "Monte Carlo" was popularized by Stanislaw Ulam, Enrico
Fermi, John von Neumann, and Nicholas Metropolis, among others; the
name is a reference to the Monte Carlo Casino in Monaco where Ulam's
uncle would borrow money to gamble. The use of randomness and the
repetitive nature of the process are analogous to the activities conducted
at a casino.
Perhaps the most famous early use was by Enrico Fermi in 1930, when he used a
random method to calculate the properties of the newly-discovered neutron. Monte
Carlo methods were central to the simulations required for the Manhattan Project,
though were severely limited by the computational tools at the time. Therefore, it
was only after electronic computers were first built (from 1945 on) that Monte Carlo
methods began to be studied in depth. In the 1950s they were used at Los Alamos
for early work relating to the development of the hydrogen bomb, and became
popularized in the fields of physics, physical chemistry, and operations research.
Applications of Monte Carlo
Monte Carlo method is used in almost every field of science, mathematics to
economics
Monte Carlo methods are very important in computational physics, physical
chemistry, and related applied fields, and have diverse applications from
complicated quantum chromodynamics calculations to designing heat shields
and aerodynamic forms.
The Monte Carlo method is widely used in statistical physics
In experimental particle physics, these methods are used for designing
detectors, understanding their behavior and comparing experimental data to
theory.
Designing a plastic scintillator
detector using Monte Carlo
Design Goals

High photon yield

Good spatial uniformity

Good timing

Low Cost

Ease for fabrication
Plastic scintillator detector
s
For total internal reflection,
sin  > sin c
c = 38.7o
Meridional and Skew Ray Mode
Along the fiber
1
1, 3
3
(a) A meridional
ray always
crosses the fiber
axis.
Meridio nal ray
Fiber axis
2
2
1
2
1
Skew ray
Fiber axis
5
3
5
4
Ray path along the fiber
4
2
3
(b) A skew ray
does not have
to cross the
fiber axis. It
zigzags around
the fiber axis.
Ray path projected
on to a plane normal
to fiber axis
Illustration of the difference between a meridional ray and a skew ray.
Numbers represent reflections of the ray.
© 1999 S.O. Kasap, Optoelectronics (Prentice Hall)
In skew ray mode
the incident angle
changes at every
reflection
Important steps of MC


Production and propagation of muons
Generation and propagation of photons inside scintillator
considering the various losses

Propagation of photons in WLS fiber

Collection of photon and convert then to photo electrons
Generating Cosmic Ray Muons
Angular distribution of muons
dN/d  cos2
 can be randomly chosen from this distribution
The probability of getting a muon between 0 to  is
P =  2 cos2 cos sin d d
602cos2 cos sin d d
= (1 – cos4)/0.9375
Hence  = cos-1(1-0.9375P)1/4
Continued
As we know the value of probability between is between 0 to 1, this can be
selected using a uniform random number between 0 to 1
P  R, R is the random number
As  distribution is uniform between 0 to 2 
Probability P =  / 2 
or  = 2P = 2R
Then distribute the muons uniformly over the surface of the scintillator.
Energy loss calculation
Though mean energy loss remains fairly constant over a large energy range
above minimum ionizing energy (Bethe-Block formula), however there is
fluctuation in the energy loss around the mean. The fluctuation of energy
loss is described by Landau distribution.
Landau distribution gives distribution of a universal
parameter called , which is independent of the
material and particle velocity.
The relation between energy loss ΔE and  is
ΔE =  [ + ln(5.596707 x109 2 )/(1- 2 )Z2+1 - 2 - E ]
Where  = (0.1536/ 2) (Z /A) S
S is mass density of the material
The probability distribution of  is used from ROOT
Mathematical function LandauI
Photon Production
Number of photons produced
N = ΔE Δl / ,  -> Energy required to produce a single photon
Δl -> Incremental path length
 = 100 eV, Navg = 20,000 /cm for vertical muons
Generate  and  for each photon
randomly from an isotropic
distribution
Track the photon and for every
reflection check for the critical angle
condition
sinI > sinc
c = 39o

Attenuation loss of photons
Scintillator is not fully transparent to the blue wave length photons because
of self-absorption in POPOP
The attenuation formula is
I = Io exp(-x/),
Here  = Mean attenuation length
I/Io = exp(-x/ )
P = 1 - I/Io = 1 - exp(-x/ ),
x = -  *ln(1-P) = -  *ln(1-R)
Determine the path length of each photon a priori and compare with total
path length traversed at each reflection.
Diffuse Reflection
Lambert’s cosine law
dI/d  cos
Hence probability of photons reflected between 0 to 
P =  cos d
/2 cos d
Hence P = sin,  = sin -1 R and  = 2R
Photon propagation in WLS fiber
Outer clad n2
Inner clad n1
Core n0
Conversion of photons to Photo-electrons

PMT converts the photons to photo-electrons. The conversion efficiency
depends on the quantum efficiency of the PMT.
Simulation Inputs
Photon statistics
No of photons produced
46,000
Fraction of photons escaped
25%
Fraction of photons lost due to
attenuation in scintillator
60%
Fraction of photons captured by
WLS fiber
15%
Trapping Efficiency in Fiber = 14% of the captured photons
Number of photon collected at PMT = 208
Collection efficiency = 0.45%
Photo electron yield and timing comparison
parallel
parallel
matrix
matrix


Photo–electron yield with Number of fibers
Ne  Nfib
Number of Fibers
The simulation code

This is a single C++ program of ~ 1000 lines of cde

The code can be compiled by g++ or C++ command

Easy to modify the inputs

Any one interested to use this can contact
[email protected]
Summary

The good agreement of simulation with
measurement would allow us to design and
optimize detector in future by doing simulation
prior to the actual construction which would save
lot of time and cost
THANKS
Plastic scintillator detector
s
Photon statistics
With Tyvek reflector
No of photons produced
46,000
No reflector
46,000
Fraction of photons escaped
25%
75%
Fraction of photons lost due to
attenuation in scintillator
60%
20%
Fraction of photons captured by
WLS fiber
15%
5%
Trapping Efficiency in Fiber:
For skew rays = 14 % (Real case)
For meridional rays = 11%
Number of photon collected at PMT = 208
Collection efficiency = 0.45%
Photo electron yield and timing comparison
parallel
parallel
matrix


Photo electron yield and timing comparison
parallel
parallel
matrix
matrix

