notebook.community

Edit and run

Determination of partial contributions of hadrons (needs to be reworked for MCEq V1.0+)

This notebook reproduces Figure 9 from the proceedings arXiv:1503.00544.

This updated version uses SIBYLL 2.3c. Results might be therefore slightly different.





In [1]:



    


#basic imports and ipython setup
%load_ext autoreload
%matplotlib inline
%autoreload 2
import os
import matplotlib.pyplot as plt
import numpy as np
import sys
sys.path.append('../MCEq')
os.chdir('..')

#import solver related modules
from MCEq.core import MCEqRun
from mceq_config import config, mceq_config_without
#import primary model choices
import CRFluxModels as pm



    











In [2]:



    


config['compact_mode'] = False
config['low_energy_extension']['enable'] = True



    











In [3]:



    


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.Thunman, None),
# Zenith angle in degrees. 0=vertical, 90=horizontal
theta_deg=60.0,
#expand the rest of the options from mceq_config.py
**config
)



    


















    






InteractionYields::_load(): Looking for /afs/ifh.de/group/that/work-af/git/MCEq/data/SIBYLL23C_yields_ledpm.bz2
extend_to_low_energies(): Low energy extension requested: /afs/ifh.de/group/that/work-af/git/MCEq/data/SIBYLL23C_yields_ledpm.bz2
extend_to_low_energies(): Saving /afs/ifh.de/group/that/work-af/git/MCEq/data/SIBYLL23C_yields_ledpm.bz2
DecayYields:_load():: Loading file /afs/ifh.de/group/that/work-af/git/MCEq/data/decays_v1.ppd
DecayYields::_load(): Warning: no decay distributions for 2112 found.
DecayYields::_load(): Warning: no decay distributions for -2112 found.

Hadrons and stable particles:

"gamma", "p", "p-bar", "n", "n-bar"

Mixed:

"pi+", "pi-", "K0L", "K-", "K+", "Lambda0", "Lambda0-bar", "Xi0-bar", 
"Xi0", "K0S", "Xi-", "Xi--bar", "Sigma-", "Sigma--bar", "Sigma+", "Sigma+-bar", 
"Omega-", "Omega--bar", "D-", "D+", "Ds-", "Ds+", "D0-bar", "D0", 
"XiC+", "XiC+-bar", "tau-", "tau+", "LambdaC+", "LambdaC+-bar", "XiC0", "XiC0-bar", 
"OmegaC0", "OmegaC0-bar", "pi0"

Resonances:

"eta", "Sigma0", "Sigma0-bar", "D*-", "D*+", "eta*", "jpsi", "phi", 
"SigmaC0", "SigmaC0-bar", "SigmaC+", "SigmaC+-bar", "SigmaC++", "SigmaC++-bar", "omega", "Xi*0", 
"Xi*0-bar", "Xi*-", "Xi*--bar", "SigmaC*0", "SigmaC*0-bar", "SigmaC*++", "SigmaC*++-bar", "SigmaC*+", 
"SigmaC*+-bar", "K*0-bar", "K*0", "K*-", "K*+", "Sigma*+", "Sigma*+-bar", "Sigma*0", 
"Sigma*0-bar", "Sigma*-", "Sigma*--bar", "etaC", "rho0", "rho+", "rho-", "Delta0", 
"Delta-", "Delta+", "Delta++", "Delta++-bar", "Delta+-bar", "Delta--bar", "Delta0-bar", "D*0", 
"D*0-bar", "Ds*-", "XiC*0", "XiC*+", "XiC*+-bar", "XiC*0-bar", "Ds*+"

Leptons:

"e-", "nue", "numu", "nutau", "antinutau", "antinumu", "antinue", "e+", 
"mu-", "mu+"

Aliases:
"k_nue", "k_numu", "k_nutau", "pi_antinutau", "pi_antinumu", "pi_antinue", "obs_nue", "obs_numu", 
"obs_nutau", "pr_antinumu", "pr_antinue", "pr_nutau", "pr_nue", "pr_numu", "obs_antinutau", "obs_antinumu", 
"obs_antinue", "pr_antinutau", "pi_nue", "pi_numu", "pi_nutau", "k_antinutau", "k_antinumu", "k_antinue", 
"k_mu-", "pi_mu+", "obs_mu-", "pr_mu+", "pr_mu-", "obs_mu+", "pi_mu-", "k_mu+"

Total number of species: 82
MCEqRun::set_interaction_model():  SIBYLL23C
InteractionYields:set_interaction_model():: Model SIBYLL23C already loaded.
InteractionYields:set_interaction_model():: Model SIBYLL23C already loaded.
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.07%
    shape      : 7216 x 7216
    nnz        : 1079188
D Matrix info:
    density    : 1.04%
    shape      : 7216 x 7216
    nnz        : 542903
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::set_density_model():  CORSIKA ('BK_USStd', None)
MCEqRun::set_theta_deg():  60.0
CorsikaAtmosphere::calculate_density_spline(): Calculating spline of rho(X) for zenith 60.0 degrees.
.. took 0.03s
MCEqRun::set_primary_model():  Thunman None

















In [4]:



    


res_groups = [
    (["D+", "D-", "D0", "D0-bar", "D*+", "D*-", "D*0", "D*0-bar"], r"$D^\pm + D^0$"),
    (["Ds+", "Ds-"], r"$D_s$"),
    (["K0S"], r"$K^0_s$"),
    (["K0L"], r"$K^0_L$"),
    (["LambdaC+", "LambdaC+-bar"], r'$\Lambda_C$'),
    (["OmegaC0", "OmegaC0-bar", "XiC+", 
      "XiC+-bar", "XiC0", "XiC0-bar", 
      "SigmaC*+", "SigmaC*++", "tau+", "tau-", 
      "SigmaC*++-bar", "SigmaC*+-bar", "SigmaC*0", 
      "SigmaC*0-bar", "SigmaC+", "SigmaC++", 
      "SigmaC++-bar", "SigmaC+-bar", "SigmaC0", 
      "SigmaC0-bar"], r"other prompt"),
    (["K0S", "K0L", "n", "n-bar", "p", "p-bar", 
      "Lambda0", "Lambda0-bar", "Sigma-", 
      "Sigma--bar", "Xi-", "Xi--bar", "Xi0", 
      "Xi0-bar", "Delta-", "Delta--bar", "Delta0", 
      "Delta0-bar", "K*+", "K*-", "K*0", "K*0-bar", 
      "Omega-", "Omega--bar", "Sigma*+", "Sigma*+-bar", 
      "Sigma*-", "Sigma*--bar", "Sigma*0", "Sigma*0-bar", 
      "Sigma+", "Sigma+-bar", "Sigma0", "Sigma0-bar", 
      "Delta+", "Delta++", "Delta++-bar", "Delta+-bar", 
      "Xi*-", "Xi*--bar", "Xi*0", "Xi*0-bar"],r"other conv. mu"),
    (["K0S", "K0L", "n", "n-bar", "p", "p-bar", 
      "Lambda0", "Lambda0-bar", "Sigma-", "Sigma--bar", 
      "Xi-", "Xi--bar", "Xi0", "Xi0-bar", "Delta-", 
      "Delta--bar", "Delta0", "Delta0-bar", "K*+", 
      "K*-", "K*0", "K*0-bar", "Omega-", "Omega--bar", 
      "Sigma*+", "Sigma*+-bar", "Sigma*-", "Sigma*--bar", 
      "Sigma*0", "Sigma*0-bar", "Sigma+", "Sigma+-bar", 
      "Sigma0", "Sigma0-bar", "Delta+", "Delta++", 
      "Delta++-bar", "Delta+-bar", "Xi*-", "Xi*--bar", 
      "Xi*0", "Xi*0-bar", "eta", "eta*", "omega", 
      "phi", "pi0", "rho+", "rho-", "rho0"],r"other conv. numu"),
    (["n", "n-bar", "p", "p-bar", "Lambda0", 
      "Lambda0-bar", "Sigma-", "Sigma--bar", 
      "Xi-", "Xi--bar", "Xi0", "Xi0-bar", "Delta-", 
      "Delta--bar", "Delta0", "Delta0-bar", 
      "K*+", "K*-", "K*0", "K*0-bar", "Omega-", 
      "Omega--bar", "Sigma*+", "Sigma*+-bar", "Sigma*-", 
      "Sigma*--bar", "Sigma*0", "Sigma*0-bar", "Sigma+", 
      "Sigma+-bar", "Sigma0", "Sigma0-bar", "Delta+", 
      "Delta++", "Delta++-bar", "Delta+-bar", "Xi*-", 
      "Xi*--bar", "Xi*0", "Xi*0-bar", "eta", "eta*", 
      "omega", "phi", "pi0", "rho+", "rho-", "rho0"],r"other conv. nue"),
    (["eta", "eta*", "omega", "phi", 
      "pi0", "rho+", "rho-", "rho0",
      "etaC", "jpsi"], r"unflavored"),
    (["tau+","tau-"], r"tau decay"),
    (["mu+", "mu-", "pi_mu+", "pi_mu-", 
      "k_mu+", "k_mu-", "pr_mu+", "pr_mu-"], r"$\mu$ decay")]

color_spectrum = ['b', 'r', 'g', 'orange', 'cyan', 'violet',
                  'brown', 'pink', 'yellow', 'lightblue']



    











In [5]:



    


mag = 3
mu_obs = {}
numu_obs = {}
nue_obs = {}
nutau_obs = {}
for res_group, res_title in res_groups:
    mceq_run.set_obs_particles(res_group)
    mceq_run.solve()
    mu_obs[res_title] = mceq_run.get_solution('obs_mu+', mag) + \
                        mceq_run.get_solution('obs_mu-', mag)
    numu_obs[res_title] = mceq_run.get_solution('obs_numu', mag) + \
                          mceq_run.get_solution('obs_antinumu', mag)
    nue_obs[res_title] = mceq_run.get_solution('obs_nue', mag) + \
                         mceq_run.get_solution('obs_antinue', mag)
    nutau_obs[res_title] = mceq_run.get_solution('obs_nutau', mag) + \
                           mceq_run.get_solution('obs_antinutau', mag)



    


















    






MCEqRun::set_obs_particles(): Converted names:D+, D-, D0, D0-bar, D*+, D*-, D*0, D*0-bar
to: 411, -411, 421, -421, 413, -413, 10421, -10421
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.48%
    shape      : 7216 x 7216
    nnz        : 1293665
D Matrix info:
    density    : 1.09%
    shape      : 7216 x 7216
    nnz        : 567537
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_calculate_integration_path(): X_surface = 2060.60278351
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   1.12ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 1.1231918335 sec
MCEqRun::set_obs_particles(): Converted names:Ds+, Ds-
to: 431, -431
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.49%
    shape      : 7216 x 7216
    nnz        : 1295104
D Matrix info:
    density    : 1.07%
    shape      : 7216 x 7216
    nnz        : 556719
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   1.11ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 1.11393809319 sec
MCEqRun::set_obs_particles(): Converted names:K0S
to: 310
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.15%
    shape      : 7216 x 7216
    nnz        : 1119114
D Matrix info:
    density    : 1.05%
    shape      : 7216 x 7216
    nnz        : 547024
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.64ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.641925096512 sec
MCEqRun::set_obs_particles(): Converted names:K0L
to: 130
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.15%
    shape      : 7216 x 7216
    nnz        : 1120604
D Matrix info:
    density    : 1.07%
    shape      : 7216 x 7216
    nnz        : 557779
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.63ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.637985944748 sec
MCEqRun::set_obs_particles(): Converted names:LambdaC+, LambdaC+-bar
to: 4122, -4122
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.40%
    shape      : 7216 x 7216
    nnz        : 1248887
D Matrix info:
    density    : 1.06%
    shape      : 7216 x 7216
    nnz        : 551811
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.68ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.688147068024 sec
MCEqRun::set_obs_particles(): Converted names:OmegaC0, OmegaC0-bar, XiC+, XiC+-bar, XiC0, XiC0-bar, SigmaC*+, SigmaC*++, tau+, tau-, SigmaC*++-bar, SigmaC*+-bar, SigmaC*0, SigmaC*0-bar, SigmaC+, SigmaC++, SigmaC++-bar, SigmaC+-bar, SigmaC0, SigmaC0-bar
to: 4332, -4332, 4232, -4232, 4132, -4132, 4214, 4224, -15, 15, -4224, -4214, 4114, -4114, 4212, 4222, -4222, -4212, 4112, -4112
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.52%
    shape      : 7216 x 7216
    nnz        : 1313834
D Matrix info:
    density    : 1.13%
    shape      : 7216 x 7216
    nnz        : 588218
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.76ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.764829874039 sec
MCEqRun::set_obs_particles(): Converted names:K0S, K0L, n, n-bar, p, p-bar, Lambda0, Lambda0-bar, Sigma-, Sigma--bar, Xi-, Xi--bar, Xi0, Xi0-bar, Delta-, Delta--bar, Delta0, Delta0-bar, K*+, K*-, K*0, K*0-bar, Omega-, Omega--bar, Sigma*+, Sigma*+-bar, Sigma*-, Sigma*--bar, Sigma*0, Sigma*0-bar, Sigma+, Sigma+-bar, Sigma0, Sigma0-bar, Delta+, Delta++, Delta++-bar, Delta+-bar, Xi*-, Xi*--bar, Xi*0, Xi*0-bar
to: 310, 130, 2112, -2112, 2212, -2212, 3122, -3122, 3112, -3112, 3312, -3312, 3322, -3322, 1114, -1114, 2114, -2114, 323, -323, 313, -313, 3334, -3334, 3224, -3224, 3114, -3114, 3214, -3214, 3222, -3222, 3212, -3212, 2214, 2224, -2224, -2214, 3314, -3314, 3324, -3324
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.32%
    shape      : 7216 x 7216
    nnz        : 1207520
D Matrix info:
    density    : 1.16%
    shape      : 7216 x 7216
    nnz        : 603044
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.68ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.681896209717 sec
MCEqRun::set_obs_particles(): Converted names:K0S, K0L, n, n-bar, p, p-bar, Lambda0, Lambda0-bar, Sigma-, Sigma--bar, Xi-, Xi--bar, Xi0, Xi0-bar, Delta-, Delta--bar, Delta0, Delta0-bar, K*+, K*-, K*0, K*0-bar, Omega-, Omega--bar, Sigma*+, Sigma*+-bar, Sigma*-, Sigma*--bar, Sigma*0, Sigma*0-bar, Sigma+, Sigma+-bar, Sigma0, Sigma0-bar, Delta+, Delta++, Delta++-bar, Delta+-bar, Xi*-, Xi*--bar, Xi*0, Xi*0-bar, eta, eta*, omega, phi, pi0, rho+, rho-, rho0
to: 310, 130, 2112, -2112, 2212, -2212, 3122, -3122, 3112, -3112, 3312, -3312, 3322, -3322, 1114, -1114, 2114, -2114, 323, -323, 313, -313, 3334, -3334, 3224, -3224, 3114, -3114, 3214, -3214, 3222, -3222, 3212, -3212, 2214, 2224, -2224, -2214, 3314, -3314, 3324, -3324, 221, 331, 223, 333, 111, 213, -213, 113
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.39%
    shape      : 7216 x 7216
    nnz        : 1242887
D Matrix info:
    density    : 1.16%
    shape      : 7216 x 7216
    nnz        : 603044
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.70ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.703466176987 sec
MCEqRun::set_obs_particles(): Converted names:n, n-bar, p, p-bar, Lambda0, Lambda0-bar, Sigma-, Sigma--bar, Xi-, Xi--bar, Xi0, Xi0-bar, Delta-, Delta--bar, Delta0, Delta0-bar, K*+, K*-, K*0, K*0-bar, Omega-, Omega--bar, Sigma*+, Sigma*+-bar, Sigma*-, Sigma*--bar, Sigma*0, Sigma*0-bar, Sigma+, Sigma+-bar, Sigma0, Sigma0-bar, Delta+, Delta++, Delta++-bar, Delta+-bar, Xi*-, Xi*--bar, Xi*0, Xi*0-bar, eta, eta*, omega, phi, pi0, rho+, rho-, rho0
to: 2112, -2112, 2212, -2212, 3122, -3122, 3112, -3112, 3312, -3312, 3322, -3322, 1114, -1114, 2114, -2114, 323, -323, 313, -313, 3334, -3334, 3224, -3224, 3114, -3114, 3214, -3214, 3222, -3222, 3212, -3212, 2214, 2224, -2224, -2214, 3314, -3314, 3324, -3324, 221, 331, 223, 333, 111, 213, -213, 113
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.38%
    shape      : 7216 x 7216
    nnz        : 1240823
D Matrix info:
    density    : 1.12%
    shape      : 7216 x 7216
    nnz        : 584047
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.69ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.695334911346 sec
MCEqRun::set_obs_particles(): Converted names:eta, eta*, omega, phi, pi0, rho+, rho-, rho0, etaC, jpsi
to: 221, 331, 223, 333, 111, 213, -213, 113, 441, 443
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.30%
    shape      : 7216 x 7216
    nnz        : 1196495
D Matrix info:
    density    : 1.04%
    shape      : 7216 x 7216
    nnz        : 542903
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.65ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.650470018387 sec
MCEqRun::set_obs_particles(): Converted names:tau+, tau-
to: -15, 15
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.51%
    shape      : 7216 x 7216
    nnz        : 1306133
D Matrix info:
    density    : 1.07%
    shape      : 7216 x 7216
    nnz        : 559040
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.73ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.734627008438 sec
MCEqRun::set_obs_particles(): Converted names:mu+, mu-, pi_mu+, pi_mu-, k_mu+, k_mu-, pr_mu+, pr_mu-
to: -13, 13, -7113, 7113, -7213, 7213, -7013, 7013
MCEqRun::_init_default_matrices():Start filling matrices. Skip_D_matrix = False
MCEqRun::_convert_to_sparse():Converting to sparse (CSR) matrix format.
C Matrix info:
    density    : 2.07%
    shape      : 7216 x 7216
    nnz        : 1079188
D Matrix info:
    density    : 1.15%
    shape      : 7216 x 7216
    nnz        : 598923
MCEqRun::_init_default_matrices():Done filling matrices.
MCEqRun::_forward_euler(): Solver will perform 1004 integration steps.
Performance:   0.61ms/iteration

MCEqRun::_forward_euler(): time elapsed during integration: 0.609716176987 sec

















In [6]:



    


mu_conv = mceq_run.get_solution('conv_mu+', mag) + \
          mceq_run.get_solution('conv_mu-', mag)

numu_conv = mceq_run.get_solution('conv_numu', mag) + \
            mceq_run.get_solution('conv_antinumu', mag)

nue_conv = mceq_run.get_solution('conv_nue', mag) + \
           mceq_run.get_solution('conv_antinue', mag)

mu_prompt = mceq_run.get_solution('pr_mu+', mag) + \
          mceq_run.get_solution('pr_mu-', mag)

numu_prompt = mceq_run.get_solution('pr_numu', mag) + \
            mceq_run.get_solution('pr_antinumu', mag)

nue_prompt = mceq_run.get_solution('pr_nue', mag) + \
           mceq_run.get_solution('pr_antinue', mag)


mu_total = mceq_run.get_solution('total_mu+', mag) + \
           mceq_run.get_solution('total_mu-', mag)

numu_total = mceq_run.get_solution('total_numu', mag) + \
            mceq_run.get_solution('total_antinumu', mag)

nue_total = mceq_run.get_solution('total_nue', mag) + \
            mceq_run.get_solution('total_antinue', mag)

nutau_total = mceq_run.get_solution('total_nutau', mag) + \
              mceq_run.get_solution('total_antinutau', mag)

mu_pi = mceq_run.get_solution('pi_mu+', mag) + \
           mceq_run.get_solution('pi_mu-', mag)

numu_pi = mceq_run.get_solution('pi_numu', mag) + \
            mceq_run.get_solution('pi_antinumu', mag)

nue_pi = mceq_run.get_solution('pi_nue', mag) + \
            mceq_run.get_solution('pi_antinue', mag)

mu_k = mceq_run.get_solution('k_mu+', mag) + \
           mceq_run.get_solution('k_mu-', mag)

numu_k = mceq_run.get_solution('k_numu', mag) + \
            mceq_run.get_solution('k_antinumu', mag)

nue_k = mceq_run.get_solution('k_nue', mag) + \
            mceq_run.get_solution('k_antinue', mag)



    











In [13]:



    


color_spectrum = ['b', 'r', 'g', 'orange', 'cyan', 'violet',
                  'brown', 'pink', 'yellow', 'lightblue']
fig = plt.figure(figsize=(11, 7.5))
fig.set_tight_layout(dict(rect=[0.01, 0.01, 0.97, 0.95]))
e_grid = mceq_run.e_grid
# plt.suptitle('{0} atmospheric lepton fluxes, {1} primary flux model.'.format(mceq_run.iamodel_name,
#                                                                              mceq_run.pmodel.sname))
plt.subplot(221)        
plt.plot(e_grid, mu_total, ls='-', lw=2, color='black', label='total')
plt.plot(e_grid, mu_conv, ls='dotted', lw=2, color='black', label=r'total conv.')
plt.plot(e_grid, mu_prompt, ls='--', lw=2, color='black', label='total prompt')
plt.plot(e_grid, mu_obs['other prompt'], ls='--', lw=1.5, color='lightblue', label=r'other prompt')
plt.plot(e_grid, mu_obs['other conv. mu'], ls='--', lw=1.5, color='grey', label=r'other conv.')
plt.plot(e_grid, mu_pi, ls='-', lw=1.5, color='red', label=r'$\pi$')
plt.plot(e_grid, mu_k, ls='-', lw=1.5, color='green', label=r'$K$')
plt.plot(e_grid, mu_obs['$D^\\pm + D^0$'], ls='-', lw=1.5, color='orange', label=r'$D^\pm, D^0$')
plt.plot(e_grid, mu_obs['$D_s$'], ls='-', lw=1.5, color='cyan', label=r'$D_s$')
plt.plot(e_grid, mu_obs['$\\Lambda_C$'], ls='-', lw=1.5, color='violet', label=r'$\Lambda_C$')
plt.plot(e_grid, mu_obs['unflavored'], ls='-', lw=1.5, color='pink', label=r'unflavored')

plt.loglog()
plt.xlabel(r"$E_{\mu}$ [GeV]")
plt.ylabel(r"$\Phi_{\mu}$ (E/GeV)$^{" + str(mag) + "}$ (cm$^{2}$ s sr GeV)$^{-1}$")
plt.xlim([10, 1e9])
plt.ylim([1e-6, 1e3])
plt.legend(loc='upper right', frameon=False, numpoints=1, ncol=2, fontsize='small')

plt.subplot(222)
plt.plot(e_grid, numu_total, ls='-', lw=2, color='black', label='total')
plt.plot(e_grid, numu_conv, ls='dotted', lw=2, color='black', label=r'total conv.')
plt.plot(e_grid, numu_prompt, ls='--', lw=2, color='black', label='total prompt')
plt.plot(e_grid, numu_obs['other conv. numu'], ls='--', lw=1.5, color='grey', label=r'other conv.')
plt.plot(e_grid, numu_obs['other prompt'], ls='--', lw=1.5, color='lightblue', label=r'other prompt')
plt.plot(e_grid, numu_pi, ls='-', lw=1.5, color='red', label=r'$\pi$')
plt.plot(e_grid, numu_k, ls='-', lw=1.5, color='green', label=r'$K$')
plt.plot(e_grid, numu_obs['$D^\\pm + D^0$'], ls='-', lw=1.5, color='orange', label=r'$D^\pm, D^0$')
plt.plot(e_grid, numu_obs['$D_s$'], ls='-', lw=1.5, color='cyan', label=r'$D_s$')
plt.plot(e_grid, numu_obs['$\\Lambda_C$'], ls='-', lw=1.5, color='violet', label=r'$\Lambda_C$')
plt.loglog()
plt.xlabel(r"$E_{\nu_{\mu}}$ [GeV]")
plt.ylabel(r"$\Phi_{\nu_{\mu}}$ (E/GeV)$^{" + str(mag) + "}$ (cm$^{2}$ s sr GeV)$^{-1}$")
plt.xlim([10, 1e9])
plt.ylim([1e-7, 1e1])
plt.legend(loc='upper right', frameon=False, numpoints=1, ncol=2, fontsize='small')

plt.subplot(223)
plt.plot(e_grid, nue_total, ls='-', lw=2, color='black', label='total')
plt.plot(e_grid, nue_conv, ls='dotted', lw=2, color='black', label=r'total conv.')
plt.plot(e_grid, nue_prompt, ls='--', lw=2, color='black', label='total prompt')
plt.plot(e_grid, nue_obs['other conv. nue'], ls='--', lw=1.5, color='grey', label=r'other conv.')
plt.plot(e_grid, nue_obs['other prompt'], ls='--', lw=1.5, color='lightblue', label=r'other prompt')
plt.plot(e_grid, nue_pi, ls='-', lw=1.5, color='red', label=r'$\pi$')
plt.plot(e_grid, nue_k, ls='-', lw=1.5, color='green', label=r'$K$')
plt.plot(e_grid, nue_obs['$K^0_s$'], ls='-', lw=1.5, color='brown', label=r'$K^0_S$')
plt.plot(e_grid, nue_obs['$K^0_L$'], ls='-', lw=1.5, color='blue', label=r'$K^0_L$')
plt.plot(e_grid, nue_obs['$\\mu$ decay'], ls='-', lw=1.5, color='pink', label=r'$\mu$ decay')
plt.plot(e_grid, nue_obs['$D^\\pm + D^0$'], ls='-', lw=1.5, color='orange', label=r'$D^\pm, D^0$')
plt.plot(e_grid, nue_obs['$D_s$'], ls='-', lw=1.5, color='cyan', label=r'$D_s$')
plt.plot(e_grid, nue_obs['$\\Lambda_C$'], ls='-', lw=1.5, color='violet', label=r'$\Lambda_C$')
plt.plot(e_grid, nue_obs["tau decay"], ls='-', lw=1.5, color='darkblue', label=r'$\tau$')

plt.loglog()
plt.xlabel(r"$E_{\nu_{e}}$ [GeV]")
plt.ylabel(r"$\Phi_{\nu_{e}}$ (E/GeV)$^{" + str(mag) + "}$ (cm$^{2}$ s sr GeV)$^{-1}$")
plt.xlim([10, 1e9])
plt.ylim([1e-7, 1e1])
plt.legend(loc=1, frameon=False, numpoints=1, ncol=3, fontsize='small')

plt.subplot(224)
plt.plot(e_grid, nutau_total, ls='-', lw=2, color='black', alpha=0.6, label='total')
plt.plot(e_grid, nutau_obs["tau decay"], ls='-', lw=1.5, color='darkblue', label=r'$\tau$')
plt.plot(e_grid, nutau_obs['$D^\\pm + D^0$'], ls='-', lw=1.5, color='orange', label=r'$D^\pm, D^0$')
plt.plot(e_grid, nutau_obs['$D_s$'], ls='-', lw=1.5, color='cyan', label=r'$D_s$')


plt.loglog()
plt.xlabel(r"$E_{\nu_{\tau}}$ [GeV]")
plt.ylabel(r"$\Phi_{\nu_{\tau}}$ (E/GeV)$^{" + str(mag) + "}$ (cm$^{2}$ s sr GeV)$^{-1}$")
plt.xlim([10, 1e9])
plt.ylim([1e-9, 1e-2])
plt.legend(loc=1, frameon=False, numpoints=1, ncol=2, fontsize='small')
plt.tight_layout(pad=0.1)
# plt.savefig('detailed_flux.pdf')



    


















    
























In [8]:

Content source: afedynitch/MCEq

Similar notebooks:

notebook.community | gallery | about