Muon spectra

Author: Hans Dembinski (https://github.com/HDembinski)

• prompt: mother has lifetime < D0

#basic imports and ipython setup
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams["font.size"] = 15
rcParams["figure.figsize"] = (7, 6)
from scipy.interpolate import interp1d
from scipy.integrate import trapezoid
import numpy as np

#import solver related modules
from MCEq.core import MCEqRun
import mceq_config as config
#import primary model choices
import crflux.models as pm

interaction_models = (
    'DPMJet-III-19.1',
    'DPMJet-III-3.0.6',
    'EPOS-LHC',
    'QGSJet-II-04',
    'QGSJet-II-03',
    'QGSJet-01c',
    'SIBYLL-2.3c',
    'SIBYLL-2.1',
)

cr_energy = {"auger": 1e10, "icecube": 1e8}
h_obs = {"auger": 1425, "icecube": 2835}
theta_deg = {"auger": 67, "icecube": 13}

lepton_spectra = {}
height_grid = {}

for experiment in ("auger", "icecube"):
    if experiment not in lepton_spectra:
        lepton_spectra[experiment] = {}
    for interaction_model in interaction_models:
        if interaction_model in lepton_spectra[experiment]:
            print "skipping", interaction_model
            continue
        lepton_spectra[experiment][interaction_model] = {}
        for prim in ("proton", "iron"):
            print "=== processing", experiment, interaction_model, prim, "==="
            spectra = lepton_spectra[experiment][interaction_model][prim] = {}

            config.h_obs = h_obs[experiment]
            mceq_run = MCEqRun(
                interaction_model,
                primary_model=(pm.HillasGaisser2012, 'H3a'),
                theta_deg=theta_deg[experiment],
            )
            mceq_run.set_density_model(
                {
                    "icecube": ("MSIS00_IC", ("SouthPole", "January")),
                    "auger": ("CORSIKA", ("BK_USStd", None)),
                }[experiment]
            )

            a = {"proton": 1, "iron": 56}[prim]
            if prim == 'proton':
                mceq_run.set_single_primary_particle(
                    cr_energy[experiment], pdg_id = 2212)
            else:
                mceq_run.set_single_primary_particle(
                    cr_energy[experiment], corsika_id = 5626)

            Xvec = np.arange(1, mceq_run.density_model.max_X, 5)
            if experiment not in height_grid:
                height_grid[experiment] = (Xvec, mceq_run.density_model.s_lX2h(np.log(Xvec))/1e2)

            mceq_run.solve(int_grid=Xvec, grid_var="X")

            # muons
            spectra["egrid"] = mceq_run.e_grid
            for t in ("total", "conv", "pr"):
                spectra["mu_" + t] = (
                    mceq_run.get_solution(t + '_mu+', 0) +
                    mceq_run.get_solution(t + '_mu-', 0)
                )

            mu_long = np.empty_like(Xvec)
            for idx, X in enumerate(Xvec):
                y = (mceq_run.get_solution("total_mu+", 0, grid_idx=idx) +
                     mceq_run.get_solution("total_mu-", 0, grid_idx=idx))
                mu_long[idx] = trapezoid(y, mceq_run.e_grid)
            spectra["mu_long"] = mu_long

            if experiment == "icecube":
                for n in ("numu", "nue", "nutau"):
                    for t in ("total", "conv", "pr"):
                        # since there are no conventional tau neutrinos, prompt=total
                        if n == "nutau" and t != "total": continue
                        spectra[n + "_" + t] = (
                            mceq_run.get_solution(t + '_' + n, 0) +
                            mceq_run.get_solution(t + '_anti' + n, 0)
                        )

=== processing auger DPMJet-III-19.1 proton ===
MCEqRun::set_interaction_model(): DPMJETIII191
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)

/mnt/c/Users/afedy/OneDrive/devel/git/MCEq/MCEq/core.py:483: LinAlgWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number9.803333e-23
  emat, b_particle)

=== processing auger DPMJet-III-19.1 iron ===
MCEqRun::set_interaction_model(): DPMJETIII191
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)

/mnt/c/Users/afedy/OneDrive/devel/git/MCEq/MCEq/core.py:492: LinAlgWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number3.902777e-19
  emat, b_protons)
/mnt/c/Users/afedy/OneDrive/devel/git/MCEq/MCEq/core.py:498: LinAlgWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number3.902777e-19
  emat, b_neutrons)

=== processing auger DPMJet-III-3.0.6 proton ===
MCEqRun::set_interaction_model(): DPMJETIII306
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger DPMJet-III-3.0.6 iron ===
MCEqRun::set_interaction_model(): DPMJETIII306
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger EPOS-LHC proton ===
MCEqRun::set_interaction_model(): EPOSLHC
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger EPOS-LHC iron ===
MCEqRun::set_interaction_model(): EPOSLHC
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-II-04 proton ===
MCEqRun::set_interaction_model(): QGSJETII04
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-II-04 iron ===
MCEqRun::set_interaction_model(): QGSJETII04
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-II-03 proton ===
MCEqRun::set_interaction_model(): QGSJETII03
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-II-03 iron ===
MCEqRun::set_interaction_model(): QGSJETII03
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-01c proton ===
MCEqRun::set_interaction_model(): QGSJET01C
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger QGSJet-01c iron ===
MCEqRun::set_interaction_model(): QGSJET01C
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger SIBYLL-2.3c proton ===
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger SIBYLL-2.3c iron ===
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger SIBYLL-2.1 proton ===
MCEqRun::set_interaction_model(): SIBYLL21
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing auger SIBYLL-2.1 iron ===
MCEqRun::set_interaction_model(): SIBYLL21
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
MCEqRun::set_density_model(): Setting density profile to CORSIKA ('BK_USStd', None)
=== processing icecube DPMJet-III-19.1 proton ===
MCEqRun::set_interaction_model(): DPMJETIII191
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00

/mnt/c/Users/afedy/OneDrive/devel/git/MCEq/MCEq/core.py:483: LinAlgWarning: scipy.linalg.solve
Ill-conditioned matrix detected. Result is not guaranteed to be accurate.
Reciprocal condition number9.803333e-19
  emat, b_particle)

=== processing icecube DPMJet-III-19.1 iron ===
MCEqRun::set_interaction_model(): DPMJETIII191
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube DPMJet-III-3.0.6 proton ===
MCEqRun::set_interaction_model(): DPMJETIII306
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube DPMJet-III-3.0.6 iron ===
MCEqRun::set_interaction_model(): DPMJETIII306
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube EPOS-LHC proton ===
MCEqRun::set_interaction_model(): EPOSLHC
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube EPOS-LHC iron ===
MCEqRun::set_interaction_model(): EPOSLHC
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-II-04 proton ===
MCEqRun::set_interaction_model(): QGSJETII04
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-II-04 iron ===
MCEqRun::set_interaction_model(): QGSJETII04
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-II-03 proton ===
MCEqRun::set_interaction_model(): QGSJETII03
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-II-03 iron ===
MCEqRun::set_interaction_model(): QGSJETII03
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-01c proton ===
MCEqRun::set_interaction_model(): QGSJET01C
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube QGSJet-01c iron ===
MCEqRun::set_interaction_model(): QGSJET01C
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube SIBYLL-2.3c proton ===
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube SIBYLL-2.3c iron ===
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube SIBYLL-2.1 proton ===
MCEqRun::set_interaction_model(): SIBYLL21
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00
=== processing icecube SIBYLL-2.1 iron ===
MCEqRun::set_interaction_model(): SIBYLL21
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
MCEqRun::set_density_model(): Setting density profile to MSIS00_IC ('SouthPole', 'January')
MSIS00IceCubeCentered::set_theta(): latitude = -90.00 for zenith angle = 13.00

Plot with matplotlib

def model_color(name):
    model_colors = {"EPOS-LHC": "k",
                    "SIBYLL" : "b",
                    "DPMJet" : "g",
                    "QGSJet" : "r",
                    "QGSJet-01": "0.5"}
    for k,v in model_colors.items():
        if name.startswith(k):
            return v

def model_linestyle(name):
    if name in ("EPOS-LHC", "SIBYLL-2.3.5", "DPMJet-III-19.1", "QGSJet-II-04"):
        return "-"
    else:
        return ":"

# atmoshpere
for experiment, (X, h) in height_grid.items():
    plt.figure()
    plt.title({
        "auger": r"Auger: 1425 m, $\theta = 67^\circ$",
        "icecube": r"IceCube: 2835 m, $\theta = 13^\circ$",
    }[experiment])
    plt.plot(h/1e3, X, "k-")
    plt.axvspan(0, h_obs[experiment]/1e3, color="k", alpha=0.2)
    x_ground = interp1d(h, X)(h_obs[experiment] + 20.)
    print "%s %.0f gcm-2" % (experiment, x_ground)
    plt.xlim(0, 50)
    plt.xlabel("$h$/km")
    plt.ylabel("$X_\mathrm{atm}/\mathrm{g\,cm^{-2}}$")
    #plt.semilogy()

auger 2207 gcm-2
icecube 716 gcm-2

# longitudinal
h_mu_max = {"auger": 0.0, "icecube": 0.0}

for experiment in lepton_spectra:
    X, h = height_grid[experiment]
    for interaction_model in interaction_models:
        if interaction_model not in lepton_spectra[experiment]:
            continue
        prim_spectra = lepton_spectra[experiment][interaction_model]

        plt.figure()
        plt.title({
            "auger": r"Auger: 1425 m, $E_0 = 10^{10}\,$GeV $\theta = 67^\circ$",
            "icecube": r"IceCube: 2835 m, $E_0 = 10^8\,$GeV $\theta = 13^\circ$",
        }[experiment], x=0.55)
        for prim, spectra in prim_spectra.items():
            y = spectra["mu_long"]
            color = {"proton" : "r", "iron": "b"}[prim]
            label = {"proton" : "p", "iron": "Fe"}[prim]
            plt.plot(h/1e3, y, "k-", color=color, label=prim)
            hm = h[np.argmax(y)]
            h_mu_max[experiment] += hm
            print "%s %s h_mu_max %.1f km" % (experiment, prim, hm/1e3)
        print "ratio", prim_spectra["iron"]["mu_long"][-1]/prim_spectra["proton"]["mu_long"][-1]
        plt.xlim(0, 20)
        plt.legend(title=interaction_model)
        plt.axvspan(0, h_obs[experiment]/1e3, color="k", alpha=0.2, zorder=1)
        plt.xlabel("$h$/km")
        plt.ylabel("$N_\mu(E_\mu > 1\,\mathrm{GeV})$")

for k, v in h_mu_max.items():
    h_mu_max[k] = v/2
    print "%s avg h_mu_max %.1f km" % (experiment, v/2e3)

auger iron h_mu_max 9.3 km
auger proton h_mu_max 8.7 km
ratio 1.4867478709731174
auger iron h_mu_max 8.7 km
auger proton h_mu_max 54.4 km
ratio inf
auger iron h_mu_max 8.4 km
auger proton h_mu_max 7.8 km
ratio 1.362216523562472
auger iron h_mu_max 8.5 km
auger proton h_mu_max 8.0 km
ratio 1.3636377503286812
auger iron h_mu_max 8.6 km
auger proton h_mu_max 8.1 km
ratio 1.382826677404677
auger iron h_mu_max 9.1 km
auger proton h_mu_max 8.7 km
ratio 1.36522174742199
auger iron h_mu_max 8.1 km
auger proton h_mu_max 7.6 km
ratio 1.360661607932406

/home/afedyni/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:24: RuntimeWarning: divide by zero encountered in double_scalars

auger iron h_mu_max 8.7 km
auger proton h_mu_max 8.2 km
ratio 1.5078626483072022
icecube iron h_mu_max 3.2 km
icecube proton h_mu_max 2.9 km
ratio 1.567291165476035
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.6327638140302534
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.5468474737090006
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.5211633633107633
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.5922146014144292
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.4904957493431992
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.5442127338752043
icecube iron h_mu_max 2.9 km
icecube proton h_mu_max 2.9 km
ratio 1.61449446747557
icecube avg h_mu_max 90.5 km
icecube avg h_mu_max 23.0 km

# muon spectra 1
for experiment in lepton_spectra:
    for interaction_model in interaction_models:
        if interaction_model not in lepton_spectra[experiment]:
            continue
        prim_spectra = lepton_spectra[experiment][interaction_model]

        plt.figure()
        plt.title({
            "auger": r"Auger: 1425 m, $E_0 = 10^{10}\,$GeV $\theta = 67^\circ$",
            "icecube": r"IceCube: 2835 m, $E_0 = 10^8\,$GeV $\theta = 13^\circ$",
        }[experiment])
        for prim, spectra in prim_spectra.items():
            x = spectra["egrid"]
            y1 = spectra["mu_total"]
            y2 = spectra["mu_pr"]
            color = {"proton" : "r", "iron": "b"}[prim]
            label = {"proton" : "p", "iron": "Fe"}[prim]
            plt.plot(x, y1, "k-", color=color, lw=2, label=label + " total")
            plt.plot(x, y2, "k:", color=color, label=label + " prompt")
        plt.loglog()
        plt.xlim(.1, 1e8)
        plt.ylim(1e-14, 1e7)
        plt.xlabel(r"$E_\mu$/GeV")
        plt.legend(title=interaction_model)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-f6539dd31670> in <module>()
      1 # muon spectra 1
----> 2 for experiment in lepton_spectra:
      3     for interaction_model in interaction_models:
      4         if interaction_model not in lepton_spectra[experiment]:
      5             continue

NameError: name 'lepton_spectra' is not defined

# muon spectra 2
from collections import OrderedDict
for experiment in lepton_spectra:
    plt.figure(figsize=(10, 5))
    plt.title({
        "auger": r"Auger: 1425 m, p, $E_0 = 10^{10}\,$GeV, $\theta = 67^\circ$",
        "icecube": r"IceCube: 2835 m, p, $E_0 = 10^8\,$GeV, $\theta = 13^\circ$",
    }[experiment], x=0.5)

    y = OrderedDict()
    for interaction_model in interaction_models:
        if interaction_model not in lepton_spectra[experiment]:
            continue
        if model_linestyle(interaction_model) != "-":
            continue
        prim_spectra = lepton_spectra[experiment][interaction_model]
        spectra = prim_spectra["proton"]
        y[interaction_model] = spectra["mu_total"]
        x = spectra["egrid"]

    for interaction_model, yi in y.items():
        plt.plot(x, yi / y["EPOS-LHC"], "-",
                 color=model_color(interaction_model), label=interaction_model)
    plt.semilogx()
    plt.ylim(0., 2)
    plt.xlim(1, 1e3)
    plt.xlabel(r"$E_\mu$/GeV")
    plt.ylabel("muon spectrum relative to EPOS-LHC")
    plt.subplots_adjust(right=0.7)
    plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))

# muon spectra ratios
for experiment in lepton_spectra:
    ratios = []

    plt.figure(figsize=(10, 5))
    plt.title({
        "auger": r"Auger: 1425 m, $E_0 = 10^{10}\,$GeV $\theta = 67^\circ$",
        "icecube": r"IceCube: 2835 m, $E_0 = 10^8\,$GeV $\theta = 13^\circ$",
    }[experiment], x=0.55)

    for interaction_model in interaction_models:
        if interaction_model not in lepton_spectra[experiment]:
            continue
        prim_spectra = lepton_spectra[experiment][interaction_model]
        ratio = prim_spectra["iron"]["mu_total"] / prim_spectra["proton"]["mu_total"]
        xiron = prim_spectra["iron"]["egrid"]
        xproton = prim_spectra["proton"]["egrid"]
        assert(np.all(xiron == xproton))

        ratios.append(ratio)
        plt.plot(x, ratio,
                 color=model_color(interaction_model),
                 ls=model_linestyle(interaction_model),
                 label=interaction_model)
        plt.semilogx()
        plt.xlim(1, cr_energy[experiment]/1e2)
        plt.ylim(1, 3.2)
        plt.xlabel(r"$E_\mu$/GeV")
        plt.ylabel(r"iron / proton")
    plt.subplots_adjust(right=0.65)
    plt.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))

    plt.figure()
    plt.title({
        "auger": r"Auger: 1425 m, $E_0 = 10^{10}\,$GeV $\theta = 67^\circ$",
        "icecube": r"IceCube: 2835 m, $E_0 = 10^8\,$GeV $\theta = 13^\circ$",
    }[experiment])

    m = np.mean(ratios, axis=0)
    d = np.zeros_like(m)
    for r in ratios:
        a = np.abs(r-m)
        d[d<a] = a[d<a]
    s = np.var(ratios, axis=0) ** 0.5
    plt.plot(xproton, m/d)
    plt.semilogx()
    plt.xlim(1, cr_energy[experiment]/1e2)
    plt.xlabel(r"$E_\mu$/GeV")
    plt.ylabel(r"(average iron/proton ratio) / (model spread)")

C:\Users\afedy\Anaconda2\lib\site-packages\ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in divide
  from ipykernel import kernelapp as app
C:\Users\afedy\Anaconda2\lib\site-packages\ipykernel_launcher.py:43: RuntimeWarning: invalid value encountered in less
C:\Users\afedy\Anaconda2\lib\site-packages\ipykernel_launcher.py:45: RuntimeWarning: invalid value encountered in divide

# neutrino spectra
for interaction_model, prim_spectra in lepton_spectra["icecube"].items():
    plt.figure()
    plt.title(r"IceCube: 2835 m, $E_0 = 10^8\,$GeV $\theta = 13^\circ$")
    for prim, spectra in prim_spectra.items():
        x = spectra["egrid"]
        y1 = spectra["numu_total"]
        y2 = spectra["numu_pr"]
        color = {"proton" : "r", "iron": "b"}[prim]
        label = {"proton" : "p", "iron": "Fe"}[prim]
        plt.plot(x, y1, "k-", color=color, lw=2, label=label + " total")
        plt.plot(x, y2, "k:", color=color, label=label + " prompt")
    plt.loglog()
    plt.xlim(1, 1e8)
    plt.ylim(1e-14, 1e7)
    plt.xlabel(r"$E_{\nu_\mu}$/GeV")
    plt.legend(title=interaction_model)