# What is wrong with this MonteCarlo? [closed]

I am making a trial for a MonteCarlo. In this code I simulate 10 protons hitting on a slab of gold:

(*Geometry of the system*)
ClearAll["Global*"]
slab = Cuboid[{0, 0, 0}, {3, 3, 0.3}];
fslab = RegionMember[slab];

(*Definition of total cross section*)
me =(*511000/29979200*)511000;
M = 938000000;
(*c=29979200;*)
z = 1;
ρ = 19.32;
τ[T_?NumericQ] := T/M ;
β[T_?NumericQ] := Sqrt[1 - (1/(τ[T] + 1))^2];
wm[T_?NumericQ] := (2 me β[T]^2)/(1 - β[T]^2);
Zs[T_?NumericQ] = z (1 - Exp[-((125 β[T])/z^(2/3))]);
(*Numero delta prodotti da protone da 100 MeV*)
Σ[en_] :=
Integrate[
0.307075*14/28.0855 10^6 ρ (
10^-4 Zs[en]^2(*Per micrometri*))/(β[
en]^2 (w + 790)^2) (1 - (β[en]^2 (w + 790))/
wm[en] + (Pi β[en] Zs[en]^2)/
137 Sqrt[(w + 790)/wm[en]] (1 - (w + 790)/wm[en])), {w, 1,
wm[en]}, GenerateConditions -> False]

(*Function to select energy of produced electron*)
pr[s_] :=
ProbabilityDistribution[(0.307075*14/28.0855 10^6 ρ (
10^-4 Zs[w]^2(*Per micrometri*))/(β[
w]^2 (w + 790)^2) (1 - (β[w]^2 (w + 790))/
wm[w] + (Pi β[w] Zs[w]^2)/
137 Sqrt[(w + 790)/wm[w]] (1 - (w + 790)/wm[w])))/
NIntegrate[
0.307075*79/19.96655 10^6 ρ (
10^-4 Zs[w]^2(*Per micrometri*))/(β[
w]^2 (w + 790)^2) (1 - (β[w]^2 (w + 790))/
wm[w] + (Pi β[w] Zs[w]^2)/
137 Sqrt[(w + 790)/wm[w]] (1 - (w + 790)/wm[w])), {w, 1,
s}], {w, 1, s}]

(*proper MonteCarlo*)
e = 10^8; (*proton initial energy*)
λ[e_] :=
1/Σ[
e];(*mean free path of proton*)
ne = 0;

elettrone = {};
dx = 0.001; (*Step in μm*)

Do[{x, y, s} = {0, 0, 0.15};
enp = e;
While[fslab[{x, y, s}], nl = -Log[RandomReal[]]; (**)
While[fslab[{x, y, s}] && nl > 0, nl -= dx/λ[e]; s -= dx];
ne++;

energy = RandomVariate[pr[wm[enp]]];
θ = ArcCos[(wm[enp] + me)/wm[enp]*energy/(energy + me)];
ϕ = RandomReal[{-Pi, Pi}];
elettrone =
AppendTo[elettrone, {{x, y, s}, θ, ϕ, energy}];
enp -= energy;], {10}]
ne
elettrone // TableForm


this produces an output of this kind:

{{{0, 0, -0.001}, 1.52521, -0.32064, 7321.34}, {{0, 0,
0.065}, 1.5549, 0.836828, 2530.24}, {{0,
0, -0.001}, 1.56153, 0.227091, 1472.38}}


The problem are the electrons with position {0,0,-0.001}: these are all spurious and I do not understand where they come from.

• What makes the positions "spurious"? If nl remains positive, the second While test will reject {0, 0, -0.001} having stepped s downwards from 0.15 in steps of -0.001 - in other words any condition that has nl remain positive makes {0,0, 0.001} a legitimite position according to your code. – gwr Aug 15 '18 at 11:51
• @gwr They are spurios because if I run the code with 10 proton, I have ten of these electrons and if I run the code with 100 proton, I have 100 of these electrons: this can't be right. – mattiav27 Aug 15 '18 at 12:48
• That is not true for me: ne will vary and quite often be larger than 10 for me. I am not a physicist, so you rather have to explain why s should not have a value of -0.001 as the logic of your algorithm will quite allow it - which is what I was trying to say. :) – gwr Aug 15 '18 at 14:05
• I note the following: $\lambda[e]$ is a constant as e will not change as opposed to enp. ne is solely determined by s and nl in your algorithm. I find your question rather hard to answer with the information given so far. – gwr Aug 15 '18 at 14:22
• Thans that is a mistake!! – mattiav27 Aug 15 '18 at 14:23

I simplified the code as much as I could. A simple example is given with {x,y,s} ={0,0,.1} and dx = .1. It is clear from the data that the boundary does not intersect.

(*Geometry of the system*)
ClearAll["Global*"]
slab = Cuboid[{0, 0, 0}, {3, 3, 0.3}];
(*Definition of total cross section*)
me = 511000;
M = 938000000;
(*c=29979200*)
z = 1;
ρ = 19.32;
τ[T_] := T/M;
β[T_] := Sqrt[1 - (1/(τ[T] + 1))^2];
wm[T_] := (2 me β[T]^2)/(1 - β[T]^2);
Zs[T_] := z (1 - Exp[-((125 β[T])/z^(2/3))]);
(*Numero delta prodotti da protone da 100 MeV*)
Σ[en_] :=
0.307075*14/28.0855*10^6 *ρ*
Integrate[(10^-4 Zs[en]^2)/(β[
en]^2 (w + 790)^2)*(1 - (β[en]^2 (w + 790))/
wm[en] + (Pi β[en] Zs[en]^2)/
137 Sqrt[(w + 790)/wm[en]] (1 - (w + 790)/wm[en])), {w, 1,
wm[en]}, GenerateConditions -> False]

(*Function to select energy of produced electron*)
pr[s_] :=
ProbabilityDistribution[(0.307075*14/
28.0855 10^6 ρ (10^-4 Zs[w]^2)/(β[
w]^2 (w + 790)^2) (1 - (β[w]^2 (w + 790))/
wm[w] + (Pi β[w] Zs[w]^2)/
137 Sqrt[(w + 790)/wm[w]] (1 - (w + 790)/wm[w])))/
NIntegrate[
0.307075*79/
19.96655 10^6 ρ (10^-4 Zs[w1]^2)/(β[
w1]^2 (w1 + 790)^2) (1 - (β[w1]^2 (w1 + 790))/
wm[w1] + (Pi β[w1] Zs[w1]^2)/
137 Sqrt[(w1 + 790)/wm[w1]] (1 - (w1 + 790)/wm[w1])), {w1, 1,
s}], {w, 1, s}]

(*proper MonteCarlo*)
e = 10^8; λ[e_] := 1/Σ[e]; ne = 0;

elettrone = {};
dx = 0.1; Do[{x, y, s} = {0, 0, 0.1};
enp = e;
While[RegionMember[slab, {x, y, s - dx}], nl = -Log[RandomReal[]];
While[RegionMember[slab, {x, y, s - dx}] && nl > 0,
nl -= dx/λ[e]; s -= dx];
ne++;
energy = RandomVariate[pr[wm[enp]]];
θ = ArcCos[(wm[enp] + me)/wm[enp]*energy/(energy + me)];
ϕ = RandomReal[{-Pi, Pi}];
elettrone =
AppendTo[elettrone, {{x, y, s}, θ, ϕ, energy}];
enp -= energy;], {10}]
ne
elettrone // TableForm


• This only changes the spurious electron from {0,0,-0.001} to {0,0,-7.63278*10^-17} – mattiav27 Aug 15 '18 at 13:14
• This is less than the machine number 0, i.e. this is 0. – Alex Trounev Aug 15 '18 at 13:38
• You calculate with automatic accuracy, therefore, the number -7.633278*10^-17 is 0. – Alex Trounev Aug 15 '18 at 13:55