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.

closed as unclear what you're asking by gwr, m_goldberg, AccidentalFourierTransform, MarcoB, Henrik SchumacherAug 16 '18 at 21:26

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• 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