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

TauRunner doesn't seem to be able to calculate the column depth for certain angles. #59

Open
williamluszczak opened this issue Jun 22, 2022 · 3 comments

Comments

@williamluszczak
Copy link

Certain choices of angle will result in TauRunner throwing an error when trying to calculat column depth. E.g:

    Earth    = construct_earth(layers=[(3., 1.0)])
    thetas_nadir = np.array([66.44148063158723])
    tracks = make_tracks(thetas_nadir)
    mytrack = tracks[thetas_nadir[0]]

    print("column depth", mytrack.total_column_depth(Earth))

Results in an error:

  File "generate_evts.py", line 116, in <module>
    print("column depth", mytrack.total_column_depth(Earth))
  File "/home/wluszczak/taurunner/taurunner/track/track.py", line 53, in total_column_depth
    set_spline(self, body)
  File "/home/wluszczak/taurunner/taurunner/track/utils/spline_column_depth.py", line 159, in set_spline
    X2x = construct_X2x(x2X)
  File "/home/wluszczak/taurunner/taurunner/track/utils/spline_column_depth.py", line 100, in construct_X2x
    raise ValueError('Constant X2x is only allow for tracks with no column depth')
ValueError: Constant X2x is only allow for tracks with no column depth```
@jlazar17
Copy link
Collaborator

TauRunner requires input in radians when running as an imported module. We have added a warning if the nadir angle is outside of [0, π].

@williamluszczak
Copy link
Author

I still seem to be getting some problems with values in that range:

thetas_nadir = np.array([2.323])
tracks = make_tracks(thetas_nadir)
mytrack = tracks[thetas_nadir[0]]

print("column depth", mytrack.total_column_depth(Earth))

Gives the same column depth error as before, but weirdly 2.32 radians works just fine

@tbertolez
Copy link

tbertolez commented Apr 3, 2024

Sorry for reopening this, but I seem to have the same problem.
Using a code very similar to the ones in the examples (but using an isotropic flux over the whole sky)

nevents  = 50 # number of events to simulate
eini     = 1e19 # initial energy in eV
th_min = 0
th_max = 180 # we want an isotropic flux over the whole 4pi solid angle
pid      = 16 # PDG MC Encoding particle ID (nutau)
xs_model = "CSMS" # neutrino cross section model 


Earth    = construct_earth(layers=[(4., 1.0)]) # Make Earth object with 4km water layer
xs       = CrossSections(xs_model)
rand  = np.random.RandomState(seed=7)

energies = make_initial_e(nevents, eini) # Return array of initial energies in eV
thetas   = make_initial_thetas(nevents, (th_min, th_max), rand = rand)

tau_prop = make_propagator(pid, Earth)

output = run_MC(energies, 
                thetas, 
                Earth, 
                xs,
                tau_prop, 
                rand,
         )

I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-4e57f02d4ef1> in <module>
     22                 xs,
     23                 tau_prop,
---> 24                 rand,
     25          )

[~/.local/lib/python3.7/site-packages/taurunner/main.py](https://vscode-remote+ssh-002dremote-002blogin-002erc-002efas-002eharvard-002eedu.vscode-resource.vscode-cdn.net/n/home11/tbertolez/IceCubeZoomEeV/~/.local/lib/python3.7/site-packages/taurunner/main.py) in run_MC(eini, thetas, body, xs, propagator, rand, no_secondaries, flavor, no_losses, condition, depth, track_type)
    210         )
    211 
--> 212         out      = Propagate(particle, my_track, body, condition=condition)
    213 
    214         if (out.survived==False):

[~/.local/lib/python3.7/site-packages/taurunner/Casino.py](https://vscode-remote+ssh-002dremote-002blogin-002erc-002efas-002eharvard-002eedu.vscode-resource.vscode-cdn.net/n/home11/tbertolez/IceCubeZoomEeV/~/.local/lib/python3.7/site-packages/taurunner/Casino.py) in Propagate(particle, track, body, condition)
     20     particle : taurunner Particle object after propagation
     21     '''
---> 22     total_column_depth = track.total_column_depth(body)
     23     total_distance     = track.x_to_d(1.-particle.position)*body.radius/units.km
     24     #keep iterating until final column depth is reached or a charged lepton is made

[~/.local/lib/python3.7/site-packages/taurunner/track/track.py](https://vscode-remote+ssh-002dremote-002blogin-002erc-002efas-002eharvard-002eedu.vscode-resource.vscode-cdn.net/n/home11/tbertolez/IceCubeZoomEeV/~/.local/lib/python3.7/site-packages/taurunner/track/track.py) in total_column_depth(self, body, safe_mode)
     51         hash_s = get_hash(self, body)
     52         if hash_s not in self._column_depth_lookup.keys():
---> 53             set_spline(self, body)
     54         return self._column_depth_lookup[hash_s][0]
     55 

[~/.local/lib/python3.7/site-packages/taurunner/track/utils/spline_column_depth.py](https://vscode-remote+ssh-002dremote-002blogin-002erc-002efas-002eharvard-002eedu.vscode-resource.vscode-cdn.net/n/home11/tbertolez/IceCubeZoomEeV/~/.local/lib/python3.7/site-packages/taurunner/track/utils/spline_column_depth.py) in set_spline(track, body, npad)
    150 
    151     # Construct the X2x on the fly
--> 152     X2x = construct_X2x(x2X)
    153 
    154     track._column_depth_lookup[hash_s] = (tot_X, x2X, X2x)

[~/.local/lib/python3.7/site-packages/taurunner/track/utils/spline_column_depth.py](https://vscode-remote+ssh-002dremote-002blogin-002erc-002efas-002eharvard-002eedu.vscode-resource.vscode-cdn.net/n/home11/tbertolez/IceCubeZoomEeV/~/.local/lib/python3.7/site-packages/taurunner/track/utils/spline_column_depth.py) in construct_X2x(x2X)
    105         cds        = [float(x2X(x)) for x in xx]
    106         cds[0]     = 0. # Sometimes spline support can be imperfect
--> 107         X2x_tck    = splrep(cds, xx)
    108         def X2x(X):
    109             return float(splev(X, X2x_tck))

/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_8_x86_64/lib/python3.7/site-packages/scipy/interpolate/fitpack.py in splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    287 
    288     """
--> 289     res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    290     return res
    291 

/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_8_x86_64/lib/python3.7/site-packages/scipy/interpolate/_fitpack_impl.py in splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    513         else:
    514             try:
--> 515                 raise _iermess[ier][1](_iermess[ier][0])
    516             except KeyError:
    517                 raise _iermess['unknown'][1](_iermess['unknown'][0])

ValueError: Error on input data

If I write th_max = 179.9, then I recover the same error message as William states. If I set th_max = 90, everything runs smoothly! Hope this helps :)

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

No branches or pull requests

3 participants