Dependence on primary cosmic ray flux#
[1]:
#import primary model choices
import crflux.models as pm
import matplotlib.pyplot as plt
import numpy as np
import mceq_config as config
#import solver related modules
from MCEq.core import MCEqRun
Create an instance of an MCEqRun class. Most options are defined in the mceq_config module, and do not require change. Look into mceq_config.py or use the documentation.
If the initialization succeeds it will print out some information according to the debug level.
[2]:
mceq_run = MCEqRun(
#provide the string of the interaction model
interaction_model='SIBYLL2.3c',
#primary cosmic ray flux model
#support a tuple (primary model class (not instance!), arguments)
primary_model=(pm.HillasGaisser2012, "H3a"),
# Zenith angle in degrees. 0=vertical, 90=horizontal
theta_deg=0.0
)
MCEqRun::set_interaction_model(): SIBYLL23C
ParticleManager::_init_default_tracking(): Initializing default tracking categories (pi, K, mu)
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
MCEqRun::set_primary_model(): HillasGaisser2012 H3a
Solve and store results#
This code below computes fluxes of neutrinos, averaged over all directions, for different primary models.
[4]:
# Bump up the debug level to see what the calculation is doing
config.debug_level = 2
[5]:
#Define equidistant grid in cos(theta)
angles = np.arccos(np.linspace(1,0,11))*180./np.pi
#Power of energy to scale the flux
mag = 3
#obtain energy grid (nver changes) of the solution for the x-axis of the plots
e_grid = mceq_run.e_grid
p_spectrum_flux = []
#Initialize empty grid
for pmcount, pmodel in enumerate([(pm.HillasGaisser2012,'H3a'),
(pm.HillasGaisser2012,'H4a'),
(pm.GaisserStanevTilav,'3-gen'),
(pm.GaisserStanevTilav,'4-gen')]):
mceq_run.set_primary_model(*pmodel)
flux = {}
for frac in ['mu_conv','mu_pr','mu_total',
'numu_conv','numu_pr','numu_total',
'nue_conv','nue_pr','nue_total','nutau_pr']:
flux[frac] = np.zeros_like(e_grid)
#Sum fluxes, calculated for different angles
for theta in angles:
mceq_run.set_theta_deg(theta)
mceq_run.solve()
#_conv means conventional (mostly pions and kaons)
flux['mu_conv'] += (mceq_run.get_solution('conv_mu+', mag)
+ mceq_run.get_solution('conv_mu-', mag))
# _pr means prompt (the mother of the muon had a critical energy
# higher than a D meson. Includes all charm and direct resonance
# contribution)
flux['mu_pr'] += (mceq_run.get_solution('pr_mu+', mag)
+ mceq_run.get_solution('pr_mu-', mag))
# total means conventional + prompt
flux['mu_total'] += (mceq_run.get_solution('total_mu+', mag)
+ mceq_run.get_solution('total_mu-', mag))
# same meaning of prefixes for muon neutrinos as for muons
flux['numu_conv'] += (mceq_run.get_solution('conv_numu', mag)
+ mceq_run.get_solution('conv_antinumu', mag))
flux['numu_pr'] += (mceq_run.get_solution('pr_numu', mag)
+ mceq_run.get_solution('pr_antinumu', mag))
flux['numu_total'] += (mceq_run.get_solution('total_numu', mag)
+ mceq_run.get_solution('total_antinumu', mag))
# same meaning of prefixes for electron neutrinos as for muons
flux['nue_conv'] += (mceq_run.get_solution('conv_nue', mag)
+ mceq_run.get_solution('conv_antinue', mag))
flux['nue_pr'] += (mceq_run.get_solution('pr_nue', mag)
+ mceq_run.get_solution('pr_antinue', mag))
flux['nue_total'] += (mceq_run.get_solution('total_nue', mag)
+ mceq_run.get_solution('total_antinue', mag))
# since there are no conventional tau neutrinos, prompt=total
flux['nutau_pr'] += (mceq_run.get_solution('total_nutau', mag)
+ mceq_run.get_solution('total_antinutau', mag))
#average the results
for frac in ['mu_conv','mu_pr','mu_total',
'numu_conv','numu_pr','numu_total',
'nue_conv','nue_pr','nue_total','nutau_pr']:
flux[frac] = flux[frac]/float(len(angles))
p_spectrum_flux.append((flux,mceq_run.pmodel.sname,mceq_run.pmodel.name))
MCEqRun::set_primary_model(): HillasGaisser2012 H3a
MCEqRun::set_theta_deg(): Zenith angle 0.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1033.81g/cm2
MCEqRun::solve(): for 588 integration steps.
solv_MKL_sparse(): Performance: 2.77ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.64sec
MCEqRun::set_theta_deg(): Zenith angle 25.84
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1148.37g/cm2
MCEqRun::solve(): for 652 integration steps.
solv_MKL_sparse(): Performance: 2.74ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.79sec
MCEqRun::set_theta_deg(): Zenith angle 36.87
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1291.43g/cm2
MCEqRun::solve(): for 730 integration steps.
solv_MKL_sparse(): Performance: 2.64ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.94sec
MCEqRun::set_theta_deg(): Zenith angle 45.57
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1475.12g/cm2
MCEqRun::solve(): for 830 integration steps.
solv_MKL_sparse(): Performance: 2.78ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.31sec
MCEqRun::set_theta_deg(): Zenith angle 53.13
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1719.54g/cm2
MCEqRun::solve(): for 962 integration steps.
solv_MKL_sparse(): Performance: 3.02ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.91sec
MCEqRun::set_theta_deg(): Zenith angle 60.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2060.60g/cm2
MCEqRun::solve(): for 1142 integration steps.
solv_MKL_sparse(): Performance: 3.41ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.91sec
MCEqRun::set_theta_deg(): Zenith angle 66.42
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2569.28g/cm2
MCEqRun::solve(): for 1402 integration steps.
solv_MKL_sparse(): Performance: 2.83ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.97sec
MCEqRun::set_theta_deg(): Zenith angle 72.54
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 3407.46g/cm2
MCEqRun::solve(): for 1806 integration steps.
solv_MKL_sparse(): Performance: 2.73ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.94sec
MCEqRun::set_theta_deg(): Zenith angle 78.46
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 5037.05g/cm2
MCEqRun::solve(): for 2496 integration steps.
solv_MKL_sparse(): Performance: 2.77ms/iteration
MCEqRun::solve(): time elapsed during integration: 6.91sec
MCEqRun::set_theta_deg(): Zenith angle 84.26
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 9422.43g/cm2
MCEqRun::solve(): for 3837 integration steps.
solv_MKL_sparse(): Performance: 2.84ms/iteration
MCEqRun::solve(): time elapsed during integration: 10.91sec
MCEqRun::set_theta_deg(): Zenith angle 90.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 36583.25g/cm2
MCEqRun::solve(): for 7445 integration steps.
solv_MKL_sparse(): Performance: 5.36ms/iteration
MCEqRun::solve(): time elapsed during integration: 39.90sec
MCEqRun::set_primary_model(): HillasGaisser2012 H4a
MCEqRun::set_theta_deg(): Zenith angle 0.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1033.81g/cm2
MCEqRun::solve(): for 588 integration steps.
solv_MKL_sparse(): Performance: 3.34ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.97sec
MCEqRun::set_theta_deg(): Zenith angle 25.84
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1148.37g/cm2
MCEqRun::solve(): for 652 integration steps.
solv_MKL_sparse(): Performance: 3.06ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.01sec
MCEqRun::set_theta_deg(): Zenith angle 36.87
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1291.43g/cm2
MCEqRun::solve(): for 730 integration steps.
solv_MKL_sparse(): Performance: 3.05ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.24sec
MCEqRun::set_theta_deg(): Zenith angle 45.57
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1475.12g/cm2
MCEqRun::solve(): for 830 integration steps.
solv_MKL_sparse(): Performance: 3.02ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.51sec
MCEqRun::set_theta_deg(): Zenith angle 53.13
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1719.54g/cm2
MCEqRun::solve(): for 962 integration steps.
solv_MKL_sparse(): Performance: 4.42ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.26sec
MCEqRun::set_theta_deg(): Zenith angle 60.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2060.60g/cm2
MCEqRun::solve(): for 1142 integration steps.
solv_MKL_sparse(): Performance: 3.95ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.53sec
MCEqRun::set_theta_deg(): Zenith angle 66.42
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2569.28g/cm2
MCEqRun::solve(): for 1402 integration steps.
solv_MKL_sparse(): Performance: 3.61ms/iteration
MCEqRun::solve(): time elapsed during integration: 5.06sec
MCEqRun::set_theta_deg(): Zenith angle 72.54
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 3407.46g/cm2
MCEqRun::solve(): for 1806 integration steps.
solv_MKL_sparse(): Performance: 3.01ms/iteration
MCEqRun::solve(): time elapsed during integration: 5.45sec
MCEqRun::set_theta_deg(): Zenith angle 78.46
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 5037.05g/cm2
MCEqRun::solve(): for 2496 integration steps.
solv_MKL_sparse(): Performance: 3.19ms/iteration
MCEqRun::solve(): time elapsed during integration: 7.97sec
MCEqRun::set_theta_deg(): Zenith angle 84.26
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 9422.43g/cm2
MCEqRun::solve(): for 3837 integration steps.
solv_MKL_sparse(): Performance: 3.00ms/iteration
MCEqRun::solve(): time elapsed during integration: 11.52sec
MCEqRun::set_theta_deg(): Zenith angle 90.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 36583.25g/cm2
MCEqRun::solve(): for 7445 integration steps.
solv_MKL_sparse(): Performance: 5.83ms/iteration
MCEqRun::solve(): time elapsed during integration: 43.44sec
MCEqRun::set_primary_model(): GaisserStanevTilav 3-gen
MCEqRun::set_theta_deg(): Zenith angle 0.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1033.81g/cm2
MCEqRun::solve(): for 588 integration steps.
solv_MKL_sparse(): Performance: 2.93ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.73sec
MCEqRun::set_theta_deg(): Zenith angle 25.84
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1148.37g/cm2
MCEqRun::solve(): for 652 integration steps.
solv_MKL_sparse(): Performance: 3.61ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.37sec
MCEqRun::set_theta_deg(): Zenith angle 36.87
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1291.43g/cm2
MCEqRun::solve(): for 730 integration steps.
solv_MKL_sparse(): Performance: 3.29ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.41sec
MCEqRun::set_theta_deg(): Zenith angle 45.57
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1475.12g/cm2
MCEqRun::solve(): for 830 integration steps.
solv_MKL_sparse(): Performance: 4.35ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.63sec
MCEqRun::set_theta_deg(): Zenith angle 53.13
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1719.54g/cm2
MCEqRun::solve(): for 962 integration steps.
solv_MKL_sparse(): Performance: 3.48ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.36sec
MCEqRun::set_theta_deg(): Zenith angle 60.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2060.60g/cm2
MCEqRun::solve(): for 1142 integration steps.
solv_MKL_sparse(): Performance: 3.62ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.14sec
MCEqRun::set_theta_deg(): Zenith angle 66.42
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2569.28g/cm2
MCEqRun::solve(): for 1402 integration steps.
solv_MKL_sparse(): Performance: 3.40ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.77sec
MCEqRun::set_theta_deg(): Zenith angle 72.54
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 3407.46g/cm2
MCEqRun::solve(): for 1806 integration steps.
solv_MKL_sparse(): Performance: 3.18ms/iteration
MCEqRun::solve(): time elapsed during integration: 5.74sec
MCEqRun::set_theta_deg(): Zenith angle 78.46
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 5037.05g/cm2
MCEqRun::solve(): for 2496 integration steps.
solv_MKL_sparse(): Performance: 3.80ms/iteration
MCEqRun::solve(): time elapsed during integration: 9.50sec
MCEqRun::set_theta_deg(): Zenith angle 84.26
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 9422.43g/cm2
MCEqRun::solve(): for 3837 integration steps.
solv_MKL_sparse(): Performance: 4.20ms/iteration
MCEqRun::solve(): time elapsed during integration: 16.14sec
MCEqRun::set_theta_deg(): Zenith angle 90.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 36583.25g/cm2
MCEqRun::solve(): for 7445 integration steps.
solv_MKL_sparse(): Performance: 6.14ms/iteration
MCEqRun::solve(): time elapsed during integration: 45.69sec
MCEqRun::set_primary_model(): GaisserStanevTilav 4-gen
MCEqRun::set_theta_deg(): Zenith angle 0.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1033.81g/cm2
MCEqRun::solve(): for 588 integration steps.
solv_MKL_sparse(): Performance: 3.29ms/iteration
MCEqRun::solve(): time elapsed during integration: 1.94sec
MCEqRun::set_theta_deg(): Zenith angle 25.84
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1148.37g/cm2
MCEqRun::solve(): for 652 integration steps.
solv_MKL_sparse(): Performance: 3.37ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.21sec
MCEqRun::set_theta_deg(): Zenith angle 36.87
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1291.43g/cm2
MCEqRun::solve(): for 730 integration steps.
solv_MKL_sparse(): Performance: 3.76ms/iteration
MCEqRun::solve(): time elapsed during integration: 2.76sec
MCEqRun::set_theta_deg(): Zenith angle 45.57
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1475.12g/cm2
MCEqRun::solve(): for 830 integration steps.
solv_MKL_sparse(): Performance: 3.71ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.09sec
MCEqRun::set_theta_deg(): Zenith angle 53.13
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 1719.54g/cm2
MCEqRun::solve(): for 962 integration steps.
solv_MKL_sparse(): Performance: 3.16ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.05sec
MCEqRun::set_theta_deg(): Zenith angle 60.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2060.60g/cm2
MCEqRun::solve(): for 1142 integration steps.
solv_MKL_sparse(): Performance: 3.25ms/iteration
MCEqRun::solve(): time elapsed during integration: 3.72sec
MCEqRun::set_theta_deg(): Zenith angle 66.42
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 2569.28g/cm2
MCEqRun::solve(): for 1402 integration steps.
solv_MKL_sparse(): Performance: 3.38ms/iteration
MCEqRun::solve(): time elapsed during integration: 4.75sec
MCEqRun::set_theta_deg(): Zenith angle 72.54
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 3407.46g/cm2
MCEqRun::solve(): for 1806 integration steps.
solv_MKL_sparse(): Performance: 5.08ms/iteration
MCEqRun::solve(): time elapsed during integration: 9.19sec
MCEqRun::set_theta_deg(): Zenith angle 78.46
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 5037.05g/cm2
MCEqRun::solve(): for 2496 integration steps.
solv_MKL_sparse(): Performance: 4.47ms/iteration
MCEqRun::solve(): time elapsed during integration: 11.16sec
MCEqRun::set_theta_deg(): Zenith angle 84.26
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 9422.43g/cm2
MCEqRun::solve(): for 3837 integration steps.
solv_MKL_sparse(): Performance: 4.50ms/iteration
MCEqRun::solve(): time elapsed during integration: 17.30sec
MCEqRun::set_theta_deg(): Zenith angle 90.00
MCEqRun::solve(): Launching euler solver
MCEqRun::_calculate_integration_path(): X_surface = 36583.25g/cm2
MCEqRun::solve(): for 7445 integration steps.
solv_MKL_sparse(): Performance: 8.72ms/iteration
MCEqRun::solve(): time elapsed during integration: 64.94sec
Plot with matplotlib#
[6]:
#get path of the home directory + Desktop
desktop = os.path.join(os.path.expanduser("~"),'Desktop')
[8]:
for pref, lab in [('numu_',r'\nu_\mu'),
('mu_',r'\mu'),
('nue_',r'\nu_e')
]:
plt.figure(figsize=(4.5, 3.5))
for (flux, p_sname, p_name), col in zip(p_spectrum_flux,['k','r','g','b','c']):
plt.loglog(e_grid, flux[pref + 'total'], color=col, ls='-', lw=2.5,
label=p_sname, alpha=0.4)
plt.loglog(e_grid, flux[pref + 'conv'], color=col, ls='--', lw=1,
label='_nolabel_')
plt.loglog(e_grid, flux[pref + 'pr'], color=col,ls='-', lw=1,
label='_nolabel_')
plt.xlim(50,1e9)
plt.ylim(1e-5,1)
plt.xlabel(rf"$E_{{{lab}}}$ [GeV]")
plt.ylabel(r"$\Phi_{" + lab + "}$ (E/GeV)$^{" + str(mag) +" }$" +
"(cm$^{2}$ s sr GeV)$^{-1}$")
plt.legend(loc='upper right')
plt.tight_layout()
# Uncoment if you want to save the plot
# plt.savefig(os.path.join(desktop,pref + 'flux.pdf'))
Save as in ASCII file for other types of processing#
[12]:
for (flux, p_sname, p_name) in p_spectrum_flux:
np.savetxt(open(os.path.join(desktop, 'numu_flux_' + p_sname + '.txt'),'w'),
zip(e_grid,
flux['mu_conv'],flux['mu_pr'],flux['mu_total'],
flux['numu_conv'],flux['numu_pr'],flux['numu_total'],
flux['nue_conv'],flux['nue_pr'],flux['nue_total'],
flux['nutau_pr']),
fmt='%6.5E',
header=('lepton flux scaled with E**{0}. Order (E, mu_conv, mu_pr, mu_total, ' +
'numu_conv, numu_pr, numu_total, nue_conv, nue_pr, nue_total, ' +
'nutau_pr').format(mag)
)