0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – gwr Aug 15 '18 at 11:51
  • $\begingroup$ @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. $\endgroup$ – mattiav27 Aug 15 '18 at 12:48
  • $\begingroup$ 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. :) $\endgroup$ – gwr Aug 15 '18 at 14:05
  • $\begingroup$ 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. $\endgroup$ – gwr Aug 15 '18 at 14:22
  • $\begingroup$ Thans that is a mistake!! $\endgroup$ – mattiav27 Aug 15 '18 at 14:23
2
$\begingroup$

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

fig1

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.