Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding Greisen Profile #15

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Open

Adding Greisen Profile #15

wants to merge 11 commits into from

Conversation

CDYMoore
Copy link

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
Copy link
Owner

vuillaut commented Sep 3, 2018

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

Copy link
Owner

@vuillaut vuillaut left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Requested changes (see comments in discussion)

@vuillaut
Copy link
Owner

vuillaut commented Sep 3, 2018

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants