Changing the Propagation Medium Density

[wjwoodley]

I (with @afedynitch) am trying to use PROPOSAL v7.4.2 to propagate muons through rock and I've noticed that when I change the density of the rock, the energy loss of the secondaries is identical, whereas at least some variation should be expected. This is essentially what I am running:

import proposal as pp

densities = [1, 5, 100] # [gcm^-3]

for i in densities:
    
    mu = pp.particle.MuMinusDef()
    rock = pp.medium.StandardRock()
    
    args = {
        "particle_def": mu,
        "target": rock,
        "interpolate": True,
        "cuts": pp.EnergyCutSettings(500, 0.05, True)
    }

    cross_sections = pp.crosssection.make_std_crosssection(**args)
    
    brems_param = pp.parametrization.bremsstrahlung.KelnerKokoulinPetrukhin(lpm=True, particle_def=mu, medium=rock)
    epair_param = pp.parametrization.pairproduction.KelnerKokoulinPetrukhin(lpm=True, particle_def=mu, medium=rock)
    ionis_param = pp.parametrization.ionization.BetheBlochRossi(energy_cuts=args["cuts"])
    shado_param = pp.parametrization.photonuclear.ShadowButkevichMikheyev()
    photo_param = pp.parametrization.photonuclear.AbramowiczLevinLevyMaor97(shadow_effect=shado_param)

    cross_sections[0] = pp.crosssection.make_crosssection(brems_param, **args)
    cross_sections[1] = pp.crosssection.make_crosssection(epair_param, **args)
    cross_sections[2] = pp.crosssection.make_crosssection(ionis_param, **args)
    cross_sections[3] = pp.crosssection.make_crosssection(photo_param, **args)

    collection = pp.PropagationUtilityCollection()
    collection.interaction = pp.make_interaction(cross_sections, True)
    collection.displacement = pp.make_displacement(cross_sections, True)
    collection.time = pp.make_time(cross_sections, mu, True)
    collection.decay = pp.make_decay(cross_sections, mu, True)
    pp.PropagationUtilityCollection.cont_rand = False

    utility = pp.PropagationUtility(collection = collection)
    detector = pp.geometry.Sphere(position = pp.Cartesian3D(0, 0, 0), radius = 10000000, inner_radius = 0)
    
    # Change the density
    
    density_distr = pp.density_distribution.density_homogeneous(mass_density=i)
    
    # Create the propagator

    propagator = pp.Propagator(mu, [(detector, utility, density_distr)])

    init_state = pp.particle.ParticleState()
    init_state.position = pp.Cartesian3D(0, 0, 0)
    init_state.direction = pp.Cartesian3D(0, 0, -1)
    init_state.energy = 1e5 + mu.mass

    # Propagate the muons

    pp.RandomGenerator.get().set_seed(500)

    # Propagate 3 km = 3e5 cm
    
    track = propagator.propagate(init_state, 3e5)

    # Print the secondaries' energies

    print(track.track_energies())

This outputs three identical lists:

[100105.6583745, 94190.74778033089, 89409.65106976482, 81521.10358506946, 80895.58736895806, 61667.89445004389, 58036.095812057545, 17981.00001545329, 17439.233219985068, 105.6583745]
[100105.6583745, 94190.74778033089, 89409.65106976482, 81521.10358506946, 80895.58736895806, 61667.89445004389, 58036.095812057545, 17981.00001545329, 17439.233219985068, 105.6583745]
[100105.6583745, 94190.74778033089, 89409.65106976482, 81521.10358506946, 80895.58736895806, 61667.89445004389, 58036.095812057545, 17981.00001545329, 17439.233219985068, 105.6583745]

Am I setting the density correctly? I seem to remember getting different results in the past (I unfortunately don't have a record of the exact PROPOSAL version I used, but it was earlier than v7.1.0); has something changed?

[Jean1995]

Hello! There are a few aspects that are important to note here.

Firstly, with the line pp.RandomGenerator.get().set_seed(500), you are resetting the random seed for every propagation process, so the same sequence of random number will be used for the different propagations (I believe this is intentional to compare the propagation process for the different mass densities, but I just wanted to mention it for the sake of completness).

Since PROPOSAL 7, we started to calculate everything internally in grammage. This means that changing the mass density does not have any effect on the calculations of energy losses, scattering, etc. (which is ok in first order approximation! more about this aspsect later...). The only difference for the three different propagations is therefore the (trivial) conversion from grammage to distance. This difference can be seen when you look, for example, at the propagated distances of the particles (i.e. print(track.track_propagated_distances())). Otherwise, all three propagation see the same cross sections, to the final energies are identical.

I've mentioned that for the calculation of energy losses, changing the mass density doesn't have an effect. This is only a valid approximation for linear effects, but doesn't work for non-linear effects, for example calculation of the LPM effect. To take these effects into account, this density correction needs to be explicitly passed to the parametrizations. You can do this by passing a density_correction argument, defining the correction factor to the standard density of the defined medium, to the parametrizations which consider the LPM effect.

See this adapted script:

import proposal as pp

densities = [1, 5, 100] # [gcm^-3]

for i in densities:
    
    mu = pp.particle.MuMinusDef()
    rock = pp.medium.StandardRock()

    
    args = {
        "particle_def": mu,
        "target": rock,
        "interpolate": True,
        "cuts": pp.EnergyCutSettings(500, 0.05, True)
    }

    cross_sections = pp.crosssection.make_std_crosssection(**args)
    
    density_correction = i / rock.mass_density # non-linear density corrections
    brems_param = pp.parametrization.bremsstrahlung.KelnerKokoulinPetrukhin(lpm=True, particle_def=mu, medium=rock, density_correction=density_correction)
    epair_param = pp.parametrization.pairproduction.KelnerKokoulinPetrukhin(lpm=True, particle_def=mu, medium=rock, density_correction=density_correction)
    ionis_param = pp.parametrization.ionization.BetheBlochRossi(energy_cuts=args["cuts"])
    shado_param = pp.parametrization.photonuclear.ShadowButkevichMikheyev()
    photo_param = pp.parametrization.photonuclear.AbramowiczLevinLevyMaor97(shadow_effect=shado_param)

    cross_sections[0] = pp.crosssection.make_crosssection(brems_param, **args)
    cross_sections[1] = pp.crosssection.make_crosssection(epair_param, **args)
    cross_sections[2] = pp.crosssection.make_crosssection(ionis_param, **args)
    cross_sections[3] = pp.crosssection.make_crosssection(photo_param, **args)

    collection = pp.PropagationUtilityCollection()
    collection.interaction = pp.make_interaction(cross_sections, True)
    collection.displacement = pp.make_displacement(cross_sections, True)
    collection.time = pp.make_time(cross_sections, mu, True)
    collection.decay = pp.make_decay(cross_sections, mu, True)
    pp.PropagationUtilityCollection.cont_rand = False

    utility = pp.PropagationUtility(collection = collection)
    detector = pp.geometry.Sphere(position = pp.Cartesian3D(0, 0, 0), radius = 10000000, inner_radius = 0)
    
    # Change the density
    
    density_distr = pp.density_distribution.density_homogeneous(mass_density=i)
    
    # Create the propagator

    propagator = pp.Propagator(mu, [(detector, utility, density_distr)])

    init_state = pp.particle.ParticleState()
    init_state.position = pp.Cartesian3D(0, 0, 0)
    init_state.direction = pp.Cartesian3D(0, 0, -1)
    init_state.energy = 1e5 + mu.mass

    # Propagate the muons

    pp.RandomGenerator.get().set_seed(500)

    # Propagate 3 km = 3e5 cm
    
    track = propagator.propagate(init_state, 3e5)

    # Print the secondaries' energies

    print(track.track_energies())

The corresponding output looks like this:

[100105.6583745, 94190.74779420358, 89409.65072571643, 81521.10325981908, 80895.5870354933, 61667.89418854576, 58036.09549228725, 17980.99964755952, 17439.232851920482, 105.6583745, 105.6583745]
[100105.6583745, 94190.74779109083, 89409.65045834142, 81521.10296014052, 80895.58674834007, 61667.89371871539, 58036.09530609559, 17980.999125967555, 17439.23233082569, 105.6583745, 105.6583745]
[100105.6583745, 94190.74785993651, 89409.65518622864, 81521.10791857465, 80895.59170423132, 61667.89928473408, 58036.10043307323, 17981.005846754393, 17439.239049203396, 105.6583745, 105.6583745]

Note that the differences are very small, since the non-linear corrections are negligible at the involved energies.