Muon bundles with MCEq?

[pkalaczynski]

Yet another question:

in arXiv:1806.04140v1 in Fig. 4.6. you show a plot done for muon bundles. If I understand correctly all plots in the paper are done with MCEq. How do you do muon bundle plots instead of single muons with MCEq? Do you somehow divide the flux by the total number of particles?

[afedynitch]

You can use MCEqRun::set_single_primary_particle(energy, pdg_id=2212) or so to compute the muon spectrum dNmu/dE from a single CR primary. The spectrum would correspond to the average spectrum of a muon bundle. The mean multiplicity is the zero-th moment = int dE dNmu/dE. Of course you can not get the fluctuations of the muon number or any type of conditional probabilities in numerics.

Edit: Previously the formula shown was for the mean energy = (int dE E dNmu/dE) / (int dE dNmu/dE)

[pkalaczynski]

Ok, will it work correctly if I just put mag=1 in the first integration like this?

muon_bundle_flux = 
(mceq.get_solution('mu+', 0) +mceq.get_solution('mu-', 0))
*(mceq.get_solution('mu+', 1, integrate=True) +mceq.get_solution('mu-', 1, integrate=True))
/(mceq.get_solution('mu+', 0, integrate=True) +mceq.get_solution('mu-', 0, integrate=True))

The plot looks rather okay-ish, but I want to make sure I'm doing it right.

I do not really need the spectrum for a single primary, but rather the total one, so I skip adding pdg_id. I could evaluate systematic uncertainties for each of the 5 primaries I use for CORSIKA simulations and then add them up but that would be superfluous I think. Total spectrum should do just fine imo.

EDIT: the y label is wrong, here I used mag=0, not 3.

[afedynitch]

The "muon bundle" is a term from underground physics and is not commonly defined. In particular it is not defined for a spectrum of CR primaries. The term inclusive flux implies that you integrate over "infite" time = over all initial CR energies. This flux would be observed by an experiment counting each single muon including single counting of muons in a bundle. Dunno if this helps, but you can look at the intro sections of https://arxiv.org/abs/1806.04140.

One bundle comes always from a single interacting cosmic ray. Of course you can compute the spectrum of the bundle muons (from set_singleprimary...). The flux of all bundles (irrespective of their muon multiplicity or spectrum) is equal to the cosmic ray "rate", because all cosmic rays produce some sort of bundle.

So you have to really think carefully about what you need. You might be after counting the number of bundles of multiplicity N that your experiment will observe at depth x and angle y. For that you need something like a survival probability, depending on the minimal detection threshold of muons Emu_min. This will be the lower integration boundary for the multiplicity computation. For instance, a bundle will be seen by ORCA if N_mu(E_mu > threshold) >= 1 (or defined by the det. efficiency from MC).

The first step would be to compute for each CR energy and incident angle dN_mu/dE at the detector. By using a weighted sum for R_bckgr(theta) = int dphi dE_CR dN_cr/dE P_detection(E_CR, theta) A_proj(theta) you can get something like the background rate. This is an approximation. The proper way would be to use an effective area A_eff(E_mu, theta), that would translate the flux at the surface to the underground flux. Fragments of what I wrote above give you some alternative to avoid a full MC eff. area. As you will see this approach will be very sensitive to the det. threshold and linearly dependent on A. But you can adjust these to data that you may already have. These values will evolve with detector construction. Just some thoughts here...

Edit: What I'm suggesting above is actually a computation of the effective area with MCEq. There are some subtle caveats, but as a first estimate it is better then nothing. The "systematics" that you refer to only affect fluxes at the surface. So it's an independent approach because Phi_mu X A_eff is independent. The approach above gives you some sort of approximate handle on A_eff systematics, for instance E_mu_min. But this stuff is usually done in MC for a reason.

[pkalaczynski]

Thanks for an elaborate answer, will need a while to digest all the information ;) Btw, I think I have found a typo in units on y-axis in FIG. 4.3. in https://arxiv.org/pdf/1806.04140.pdf. It should be (cm^2 s sr GeV), not (cm^2 s sr GeV^2), right?

[afedynitch]

Cheers. Nope, the flux is multiplied with E^3 that gives you another GeV^3 factor.

[pkalaczynski]

Oh, so it's the unit for the whole expression E^3*phi? But is then the unit of phi (GeV^2)/(cm^2 s sr)? I would expect 1/(GeV (cm^2) s sr), especially comparing to Fig. 4.4

[pkalaczynski]

Okay, after being occupied with some other stuff, I finally managed to produce this plot:

It shows a comparison with CORSIKA at sea level (same hadronic model etc.). I would say this looks quite ok especially given that MCEq muon bundles are averaged (and it's obviously not really true that their average multiplicity is constant with respect to energy). Bands I just got from varying the zenith angle.

This time I just did (that's for some1 else looking for how to do it in the future I guess):

from scipy.integrate import simps
muon_flux = (mceq.get_solution('mu+', mag) + mceq.get_solution('mu-', mag))
muon_flux_integrated = simps(muon_flux, mceq.e_grid)
muon_flux_E_integrated = simps(muon_flux*mceq.e_grid, mceq.e_grid)
muon_bundle_flux = muon_flux*muon_flux_integrated/muon_flux_E_integrated

[pkalaczynski]

Hmm, I'm having some afterthoughts about the units. From the formula for the mean multiplicity you gave me: mean multiplicity = (int dE E dNmu/dE) / (int dE dNmu/dE) it seems this mean multiplicity will be in units of energy, so I guess I should in fact still divide by the bin width? After testing it by just doing mean multiplicity = ((int dE E dNmu/dE) / (int dE dNmu/dE) )/E , I get this: , which is more similar to what we have with CORSIKA (bundles dominate the high- and single muons the low-energy flux).

[afedynitch]

Hi!

hmm, actually I was not thinking properly. The formula is the mean energy of course. The mean "bundle" multiplicity is just the denominator = \int_0^\inf dE dN/dE for a single primary. Multiplicity is a scalar and not a spectrum.

The agreement between CORSIKA and MCEq for muon spectra and multiplicity should be within ~few %, or depending on the model <1%.

I don't understand what you mean by the "bundle spectrum" on the plot. Is it as a function of primary cosmic ray energy, or for a certain fixed energy the spectrum of secondary muons?

[pkalaczynski]

Hey,

this is however then just mean multiplicity over all the muon energies, not exactly what I need I would say.

Well, for CORSIKA I may be weighting wrong, I have not done flux plots before, only rates and unweighted, so there may be some detail I'm missing.

By bundle spectrum I mean bundle flux as a function of bundle energy:

Right now I'm having some doubts if this is really calculable from muon spectra. For every bundle we have:

(by <M> I mean the average multiplicity)

and

It's not clear to me how to go from to

[pkalaczynski]

Btw, using set_single_primary_particle would require me to scan over a range of primary energies for all 5 primaries I use (p,He,C,O,Fe) and I don't see the gain, since I just use GST3 model as my default, which anyway only uses those 5 primaries, so I'd expect I'll get more or less the same just by setting primary_model = (crf.GaisserStanevTilav,'3-gen') (unless I miss something)?

[afedynitch]

Again, a bundle is only a meaningful definition for a single primary. The rate of bundles is identical to the rate of cosmic rays because each cosmic ray produces one bundle when interacting.

You can try computing the mean multiplicity and then multiply it by the mean energy as previously defined. You can compute all the values from a spectrum. Mean multiplicity, mean energy and total energy. And compare <"M"> with \sum E_\mu, if this is a good approximation or not.

Edit: The M was not shown before...

[afedynitch]

Hi! Is there anything left to discuss on this? If not, I will close this issue

[pkalaczynski]

I think not for now. After some considerations, I decided to just use single muon flux, assuming the muon vs muon bundle differences will cancel out (I use it to obtain the relative systematic uncertainties). I will reopen if anything changes.

[afedynitch]

Makes sense!

[pkalaczynski]

Hi Anatoli,

I have tried to actually reproduce the few % agreement (at least for single muons), which you mentioned, but I fail to do that:

What I do for MCEq:

from MCEq.core import config, MCEqRun
import crflux.models as crf
import matplotlib.pyplot as plt
from scipy.integrate import simps

# Initalize MCEq by creating the user interface object MCEqRun
mceq = MCEqRun(
    e_min=1e2,e_max=1e12,
    interaction_model='SIBYLL23C', # High-energy hadronic interaction model
    primary_model = (crf.GaisserStanevTilav,'3-gen'), # cosmic ray flux at the top of the atmosphere
    density_model = ('CORSIKA',("KM3NeT", None)), # atmosphere model
    theta_deg = 0. # zenith angle
)

# Solve the equation system
# mceq.set_single_primary_particle(E=1e8,pdg_id=1000260560) # proton: 2212, iron: 1000260560
mceq.solve()

# Obtain the result
# Multiply fluxes be E**mag to resolve the features of the steep spectrum
mag = 0
muon_flux = mceq.get_solution('mu+', mag) + mceq.get_solution('mu-', mag)
muon_flux_integrated = simps(muon_flux, mceq.e_grid)
muon_flux_E_integrated = simps(muon_flux*mceq.e_grid, mceq.e_grid)
muon_bundle_flux = muon_flux*(muon_flux_integrated/muon_flux_E_integrated)*mceq.e_grid

# plotting:
plt.loglog(mceq.e_grid, muon_flux*1e4, label='muons (MCEq)') # *1e4 because MCEq output normally is in 1/cm^2
plt.loglog(mceq.e_grid, muon_bundle_flux*1e4, label='bundles (MCEq)')

For CORSIKA, I extract the info with km3net-specific software, but what I do is:

• I save the energy at sea (be it of individual muons or for the bundles)
• I compute the flux weight by computing the CR flux (GST3) from the primary energy & primary id, then divide the flux by the number of generated showers per primary (from all the runs I use)
• Do the weighted histogram

I use CORSIKA v7.7410 with SIBYLL2.3d and UrQMD-1.3c.

[afedynitch]

Hi.

You’re still computing the multiplicity in MCEq incorrectly and it is the mean energy or something like this. You need to do the same thing as in CORSIKA. Use set_single_primary_particle and compute the muon yields at the surface or underwater for each primary energy and species. The integral of get_solution will give you the number of muons per primary. Then you can weight it with your CR flux model.

[pkalaczynski]

Yeah, I have to update that.

But still, what about the (individual) muon flux? This should match CORSIKA pretty much out of the box as I understand.

[pkalaczynski]

By the way, I am a bit confused by what I get when doing set_single_primary_particle:

I would expect the predicted flux of muons created in showers caused by 1EeV protons to be much lower than the overall muon flux prediction ...

What I do to obtain the orange curve is I just add

mceq.set_single_primary_particle(E=1e9,pdg_id=2212)

right before mceq.solve() and proceed the same way as for the regular muon flux.

[pkalaczynski]

I've tested the multiplicity evaluation and this seems to indeed work rather nicely: (I tested for proton showers only)

[afedynitch]

Indeed. Seems to look fine. On the log scale it is not easy to see that it hits indeed the mean but the distributions look quite wide and the mean may not have a very rigorous meaning.

Is this underground, or surface?

Regarding your question above,

1EeV protons to be much lower than the overall muon flux prediction

The distributions show 2 different things, the units are different (1/cm2 GeV sr cm2) for the flux and (1/GeV) for the result from set_single_primary, which gives the average yield of muons produced by 1 primary CR. To get the bundle multiplicity, you need to multiply each yield by the CR flux evaluated at the set_single_pr.. energy times the bin width of this the energy grid. Then sum over the all weighted yields to obtain the flux.

This is shown, for example, in Eq. 1 of this paper. The w is the CR flux integrated between the bin edges, the Y would be the yield obtained from MCEq.

[pkalaczynski]

Is this underground, or surface?

That's at the sea surface. I did not want to include additional bias due to propagation etc. in the check.

The distributions show 2 different things, the units are different (1/cm2 GeV sr cm2) for the flux and (1/GeV) for the result from set_single_primary, which gives the average yield of muons produced by 1 primary CR. To get the bundle multiplicity, you need to multiply each yield by the CR flux evaluated at the set_single_pr.. energy times the bin width of this the energy grid. Then sum over the all weighted yields to obtain the flux.

Oh, ok, if it's in different units, then it makes sense. I will try doing that.

[pkalaczynski]

It seems that I pretty much reproduce the muon flux by multiplying by the CR flux (I use GST3) and Eprim binwidth and summing over all Eprim bins:

[afedynitch]

Seems ok. Now, the rate of bundles at the surface is equal to the CR all-particle flux because at these high energies each cosmic ray produces one or more muons. But at underground/-water, only a few muons survive and the rate of bundles becomes smaller.

If you have tools to propagate the muon spectrum per primary, then you would take the muon yield (let’s call it spectrum) from the surface, calculate the spectrum at a certain depth and integrate over it to get the expected mean number of muons, this will give you the average bundle multiplicity.

Then, it depends of what you wanna know and what the goal of the calculation is.

Convolving the underground multiplicities with the CR flux (similar to what you did above) will give you a multiplicity spectrum at a certain depth. Or, calculating the rate of bundles would be a similar integration, but there you would introduce weights like (w=1 for Nmu>=2 and otherwise 0) to compute bundle rates. Note that the multiplicity spectrum drops very rapidly and is subject to all kinds of fluctuations. Therefore, mean numbers may not represent well the physics and the observations.

[pkalaczynski]

Alright, I think this is solved, closing (hopefully for the last time 😅 )