# The correct assignment of probability distribution to generate random variables [closed]

I'm asking you for a help to solve a difficult (for myself) task. Now I will tell you what and why I'm trying to implement.

The task is to simulate the passage of elementary particles through a material consisting of nano-objects. To solve the task, you need to generate many random variables - the angles, on which an elementary particle will be scattered after the interaction with nano-objects. The random variables are distributed according to the law - the differential cross section of scattering, which has a rather complicated form.

In the beginning sets the scattering amplitude $fB[theta, v, npr]$, where $theta$ is the scattering angle, $v$ is the velocity of an elementary particle, $npr$ is the radius of the nanoparticle:

(* h - Planck constant *)
(* m - the mass of elementary particle *)
k[V_] := m V/h;
q[theta_, v_] := 2 k[v] Sin[theta/2];
fBintr[theta_, v_, npr_] := -((2 m)/h^2) Uvnyt (Sin[npr q[theta, v]] -
npr Cos[npr q[theta, v]] q[theta, v])/q[theta, v]^3;
fBextr[theta_, v_, npr_] := -((2 m)/h^2) Uvnesh (Sin[npr q[theta, v]] -
npr Cos[npr q[theta, v]] q[theta, v])/q[theta, v]^3;
fB[theta_, v_, npr_] := fBintr[theta, v, 0.75 npr] + fBextr[theta, v, npr];

As you can see, $fB[theta, v, npr]$ consists of two parts - internal and external, which describe a two-stage structure optical potential of the nano-object. Thus the nanoparticle consists of a solid core and a shell with different chemical composition and different optical interaction potential. The differential scattering cross-section is a square of the amplitude:

dsigma[theta_, v_, npr_] := Abs[fB[theta, v, npr]]^2;

Using the characteristic values for the task - $theta$ angle ranges [0, Pi], the velocity $v$ = 100 m/s, the radius of the nanoparticles $npr$ = 2.5 nm, we obtain the following graph:
And all looks good, but actually it is not. The differential scattering cross section can also be expressed using spherical Bessel functions of the first kind SphericalBesselJ that have instability at zero. It looks like this:

Further, if we normalize this distribution to use as PDF $dsigma[t, 100, 2.5*10^{-9}]/a$:

mydist[v_, npr_] := With[{a = NIntegrate[dsigma[tt, v, npr], {tt, 0, Pi}]},
ProbabilityDistribution[{"PDF", dsigma[t, v, npr]/a}, {t, 0, Pi},
Assumptions -> {v > 0, npr > 0}]
];
RandomVariate[mydist[100, 2.5*10^-9], {10}]

Then RandomVariate it just refuses to work and to generate random angles. I think that this is due to the instability near zero, and because $dsigma$ has a large number of $theta$ cycles. So Mathematica is unable to verify that I specified a valid PDF because it loses too much precision in the integration and the integral $dsigma[theta, 100, 2.5*10^{-9}]/a$ will not be equal to 1.

I have tried several options for generating random angles.

1. Using the RandomReal I generated uniformly distributed random value between the lower boundary $dsigma[theta, 100, 2.5 * 10 ^{-9}] / a$ (non-zero) to its upper boundary Limit[dsigma[theta, 100, 2.5*10^(-9)]/a, theta->0]. Then using FindRoot was looking for a matching $theta$. It works. But the problem is that such angles during the solution must generate millions, if not more. Accordingly, the method is very time-consuming in comparison with the generation of a large number of random variables using RandomVariate.
2. Had a look at the nice tutorial how to Create Your Own Distribution with Oleksandr Pavlyk (and actually more confused). Learned another way to generate.

RandomRandomRealRejection[ densityFunc, { hatFunc, hatRandomVector }, len, prec ] generates random variates from continuous univariate density. The density need not be normalized, as long as densityFunc <= hatFunc, and hatRandomVector draws from distribution with density ~ to hatFunc.

a = NIntegrate[dsigma[t, 100, 2.5*10^(-9)], {t, 0, Pi},
WorkingPrecision -> MachinePrecision];
b = 0.01 + Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a;
c = dsigma[Pi, 100, 2.5*10^(-9)]/a;

(* I am using piecewise *)
RandomRandomRealRejection[
Function[x,
Piecewise[{{Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a, x == 0},
{dsigma[x, 100, 2.5*10^(-9)]/a, x > 0}}], Listable],
{b &, Function[{len, prec}, RandomReal[Pi, len, WorkingPrecision -> prec]]},
10, MachinePrecision
]

Show[Histogram[%, Automatic, "PDF"],
Plot[Piecewise[{{Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a, x == 0},
{dsigma[x, 100, 2.5*10^(-9)]/a, x > 0}}], {x, 0, Pi}, PlotStyle -> Red,
PlotRange -> All], Plot[b, {x, 0, Pi}, PlotStyle -> Red, PlotRange -> All]
]

And in the end successfully generated random angles:
The problem is that, as I said earlier, these angles need to have tens of millions, and the generation of only $10^5$ was approximately 8 minutes (+/- 1 minute). Thus, this option is not very suitable for a large number of calculations.

The question I would like to ask is how to specify the probability distribution in my case? I suspect that it is not enough to use the correct method of random variables generating, but it is necessary to find a suitable algorithm for the distribution assignment.