Cascade: Problems with Multiple Scattering (in combination with tracking)

In this issue, I describe a current problem with Cascade on the master branch in detail. Afterwards, I share some thought about both a possible short-term solution for this problem, as well as the implications on a long-term solution that we should find.

What is currently wrong in Cascade

In a very simplified way, the order of the steps that are currently done in Cascade::step look like this:

1. Sample the grammage until the next interaction occurs
2. Determine the track that the particle will take to the next interaction (Taking into account magnetic fields and possible intersections with other media)
3. Apply continuous processes (doContinuous) on the calculated track
4. Sample and calculate interaction (interaction) that happens at the end of the track

Now, let's look at the way that step 3 is currently done in Cascade::step (L207-L223):

    // apply all continuous processes on particle + track
    if (sequence_.doContinuous(particle, step, limitingId) ==
        ProcessReturn::ParticleAbsorbed) {
      CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
                        particle.getPID(), particle.getEnergy() / 1_GeV);
      if (particle.isErased()) {
        CORSIKA_LOG_WARN(
            "Particle marked as Absorbed in doContinuous, but prematurely erased. This "
            "may be bug. Check.");
      } else {
        particle.erase();
      }
      return; // particle is gone -> return
    }
    particle.setTime(particle.getTime() + step.getDuration());
    particle.setPosition(step.getPosition(1));
    particle.setDirection(step.getDirection(1));

We see that after doContinuous has been called, the time, position and direction of our particle are set using the information provided by the track.

Let's look at what proposal::ContinuousProcess::doContinuous is doing:

  inline ProcessReturn ContinuousProcess<TOutput>::doContinuous(TParticle& vP,
                                                                TTrajectory const& track,
                                                                bool const) {
    // ...
    auto final_energy = (c->second).disp->UpperLimitTrackIntegral(
                            vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
                        1_MeV;
    auto dE = vP.getEnergy() - final_energy;

    // if the particle has a charge take multiple scattering into account
    if (vP.getChargeNumber() != 0) scatter(vP, dE, dX);
    vP.setEnergy(final_energy); // on the stack, this is just kinetic energy, E-m
    // ..
  }

One step is that PROPOSAL calculates continuous losses and sets the particle energy accordingly. But more importantly, PROPOSAL calls proposal::scatter. This is how proposal::scatter looks like:

  inline void ContinuousProcess<TOutput>::scatter(TParticle& particle,
                                                  HEPEnergyType const& loss,
                                                  GrammageType const& grammage) {

    // ...

    // Cast corsika vector to proposal vector
    auto particle_dir = particle.getDirection();
    auto d = particle_dir.getComponents();
    auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
                                           d.getZ().magnitude());

    // ...

    [[maybe_unused]] auto [unused1, final_direction] =
        PROPOSAL::multiple_scattering::ScatterInitialDirection(direction, deflection);

    // update particle direction after continuous loss caused by multiple
    // scattering
    particle.setDirection(
        {particle_dir.getCoordinateSystem(),
         {final_direction.GetX(), final_direction.GetY(), final_direction.GetZ()}});
  }

In the last three lines, PROPOSAL changes the direction of our particle to account for the effects of multiple scattering. But if we look back at the code of Cascade::step, this change of the particle direction is actually overwritten by this line:

    particle.setDirection(step.getDirection(1));

Conclusion 1: The direction changes due to multiple scattering, calculated by PROPOSAL, are not applied on the current master branch!

Why did this happen?

To understand what has caused this error, let's look at the doContinuous call inside Cascade::step again:

    // apply all continuous processes on particle + track
    if (sequence_.doContinuous(particle, step, limitingId)

Both particle and step carry information about the direction (and position) of the particle. Firstly, it is not clear whether doContinuous should change the direction of the particle using the track object or using the particle object (or whether doContinuous is even allowed to change these objects?). Until now, PROPOSAL::doContinuous only touched particle to change the directional information.

Secondly, is is not clear which directional information particle is carrying here: The direction of the particle at the beginning of the track, or the direction of the particle at the end of the track? (Exactly this question has already been raised in issue #391). The problem is that this behavior has changed in this commit - Before this commit, the direction of particle passed to doContinuous actually corresponded to the direction at the end of the track. However, the PROPOSAL interface has not been adapted to this change, which finally broke the treatment of multiple scattering in cascade.

Conclusion 2: Not adapting the PROPOSAL interface after a major change in Cascade, as well as ambiguities in the directional information, caused this bug.

Why is it not trivial to fix this?

Let's assume that we don't want to change the code in Cascade::step. This means that the continuous processes can only change the direction (and position) of the particle by altering the track object. However, the direction of the particle at the end of the track can not be changed this way: track is basically just a connection of two points (either by a straight or a curved line), and just changing the direction at the end of the track while not altering the rest of the trajectory is geometrically impossible.

How can this problem still be fixed?

This is an issue on the current master branch, and it makes predictions of particle distributions unphysical (see, for example, slide 2 of this talk to see the effect of neglecting multiple scattering on the lateral particle distributions). This means that we should try to find a quick fix as soon as possible.

But furthermore, there are several conceptional issues coming along with this issue, that will need a long term solution to be solved.

Quick fix

As a quick fix, it would be desirable to somehow reproduce the former behavior in Cascade. This means that doContinuous can influence the final particle direction after a continuous step.

For this, we probably need to "hack" a solution inside Cascade (while also making some changes to proposal::doContinous). This could, for example, look something like this:

if (particle.getDirection() != step.getDirection(1)) {
   // we know that something inside `doContinuous changed,
   // so we don't want to overwrite these changes
} else {
   particle.setDirection(step.getDirection(1));
}

@Nikos and I will try to find such a quick-fix solution and open a MR for it, so we can discuss it there.

Long-term solution

I would like to highlight two problems with the current treatment of multiple scattering:

1. We can't consider lateral displacement in multiple scattering right now

Even if we apply a quick-fix, multiple scattering is currently applied by changing the particle direction at the end of a tracking step. However, for a complete treatment of multiple scattering, we would also need to take into account a lateral displacement of the particle. This is visualized in this sketch from the EGS5 manual:

We have directional changes, denoted by the angles Θ and Φ (which we take into account), but also a lateral displacement, denoted by Δx and Δy (which we do not take into account). Furthermore, it should be noted that there is also going to be a difference between s and t (i.e. the length traveled by the particle, and therefore also the time that passed. This is also not taken into account).

However, if we want to take these effects into account, we also have to think of an efficent mechanism that avoids that our track leaves the current medium (without us noticing) due to multiple scattering effects.

If we neglect the lateral displacement, this means that we need to limit the maximum possible step sizes to still achieve physical results (This has been done in EGS4. In their manual, they also state under which conditions this approximation is valid). This might not be a big problem for air shower simulations, but might be problematic if we want to go to other applications, for example simulations in more dense media. (For more information about this topic, see for example the article about the Presta algorithm.)

2. We need to take into account the influence of magnetic fields on the track

We also need to consider that multiple scattering needs to be treated in parallel to the magnetic deflections. This is non-trivial, even for the sole directional changes, since we are not working with straight-lines anymore.

Conclusion concerning a long-term solution

Both points lead to the conclusion that it is inevitable to have a combined treatment of Multiple Scattering and Tracking. Separating both, as it is done right now, is very error-prone and can easily lead to inconsistencies (as it has right now). When thinking about such a concept, this should be combined with the ideas on a redesign of the continuous processes, as described in issue #391.