Adding Greisen Profile

[CDYMoore]

Griesen profile added to sky_objects.py allowing for random distribution of particles along Z axis with a distribution dependent on the initial energy and interaction height of the shower.

X and Y profiles set by random Gaussian distributions scaled such that on average 90% of particles reside within 1 Moliere Radius of the shower axis (at sea level).

If shower = pschitt.sky_objects.shower:
-> Simulation using Greisen Function activated with "shower.Greisen_Profile()"
-> New shower input parameters needed in run script:
--> "shower.energy_primary = number[eV]" : sets energy of the primary gamma-ray photon.
--> "shower.height_of_first_interaction = height[m]" : sets the height at which the initial photon interacts with the atmosphere.
--> "shower.scale = number[-]" : sets the resolution of the distribution function produced from the Greisen Function. A lower number corresponds to a higher number of simulated particles. (Recommended range: 1-100).

Inputs no longer needed in run script:
-> "shower.number_of_particles" : determined by primary energy.
-> "shower_top/length/width : determined by preprogrammed functions (see top of comment).

[vuillaut]

Hi @CDYMoore and thank you very much for your contribution.
As this is a large PR, there will probably be a lot of questions and several checking steps before merging it.

First question, have you merged the last master version of pschitt in your version? If not, please do.

About the code:

• could you please use clear, self-explanatory names for all variables? (avoid the use of m, n, o e.g. that make reading the code difficult)
• why is a scale factor required? IMO, the variable shower.number_of_particles should give the scaling and number of particles. This is how other distributions have been designed and I think a user would expect it to work this way.
• the Greisen Formula should be defined as a proper function (not inside another function). If there are physical limits to variables, checks (assert) might be a good thing to add.
• as written, the Greisen Formula function works with x: numpy.ndarray as input. This is a very good thing that should be used as an advantage to avoid loops. As a general comment, one should avoid loops at all cost in Python and take advantage as much as possible of numpy vectorised operations. Please replace loops whenever possible (e.g. lines 676, 714, 717)
• files in pschitt/tests are unit test. Could you please add unit tests related to the functions you added in this PR? These tests should go in pschitt/tests/test_objects.py.
• please remove files pschitt/tests/CDYM_Tests/description.txt and pschitt/tests/CDYM_Tests/main_28818.py from PR

Example to avoid loops:

def Greisen_Formula(x, Tmax):
    T = x / Xo # Slant depth
    s = ( (3 * T) / (T + (2 * Tmax)) ) # Shower age
    Ne = ( 0.31 / np.sqrt(Tmax) ) * np.exp(T) * (s**((-3*T)/2))
    return Ne

Tmax = 10
X = np.array([1, 2, 3])
%timeit gf1 = np.array([Greisen_Formula(x, Tmax) for x in X])
gf1 = np.array([Greisen_Formula(x, Tmax) for x in X])
%timeit gf2 = Greisen_Formula(X, Tmax)
gf2 = Greisen_Formula(X, Tmax)
print((gf1 == gf2).all())

31.6 µs ± 2.22 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
18.5 µs ± 1.24 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)
True

[vuillaut]

Requested changes (see comments in discussion)

[vuillaut]

I now better understand the need for the scaling factor. However, I would suggest to remove it and have a number of particles = min(shower.number_of_particles, physical_number_of_particles) and explain this in the docstring of the shower.Greisen_Profile() function. Let me know your thoughts.