ObservationPlanes show non-existing particles

This is a follow-up to the discussion on the slack development channel.

I've been simulating electromagnetic (1TeV initial electron) showers with the examples/corsika.cpp script. The only option I have changed is the emcut, which I changed to 20MeV.

If I look at the longitudinal profile (here for the charged particles), they look like this:

This means that I should expect zero particles on my observation plane. But if I look at the output from the ObservationPlane, I have 8931 entries in the parquet file! The energy distribution of these particles look like this:

Note that a lot of these particles have an energy very close to the cut energy, but not all of them.

For a 10 TeV shower, we would expect some charged particles according to the observation plane:

But the output of the observation plane shows 91676 particles!

This might, but doesn't have to be, related to issue #572.

added Bug Output labels

changed the description

Also from the slack discussion:

here it looks like the most of these particles are electrons and positrons. @JAlameddine Could you look at the distribution of ids?

Do they all have exactly the same energy?

no (see histogram above)

but a lot are at the cut energy

They don't have identical energies, no. Just a lot of them are very close to the cut, as Felix said

What about only those below the cut? The ones in Felix' copied table all have at least 5 common digits.

The issue is also present for photons if I remember correctly. But I will have a closer look later.

I just checked: 7202 out of 8931 charged particles have energies below the 20 MeV cut (0.0194890010539 GeV in kinetic energies).

They don't have all identical energies, but there are some accumulations. For example, there are 5924 particles which all have the (numerically) identical energy of 0.019486999139190 GeV. Note that they seem to have very different x and y coordinates on the observation plane:

Is this energy equal to 20 MeV - electron mass?

Which energy? The 0.019486999139190 GeV where we have a lot of particles? No, this is 2 keV below the cut.

Hmm so in my calculator 0.02 GeV - 0.019486999139190 GeV = 0.000513001 GeV which is rather close to the mass of the electron

doesn't explain anything though

Good point. All these 5924 particles don't have the exact cut value, but they are the particles with the energies closest to the cut value! Interesting.

@JAlameddine I assume this is current master? could you add the commit?

d356fdd1

This commit, so yes, it's the current master

I just re-checked: I don't see this unreasonable amount of photons! Only for charged particles!

I have now tried turning off the magnetic field and multiple scattering (which would be an obvious candidate for causing issues since these processes only affect charged particles), but this didn't solve the issue...

So maybe this is because the photons are massless??

just what I was thinking now.. maybe we stumbled over our confusing interface somewhere.

Nope energies are ok. It's the order of the secondary processes! Check this out:

[corsika:debug   (ParticleCut.inl:93)] ParticleCut: checking e+ (-11), E_kin= 0.019826012093476522 GeV, E=0.020337012\
093476523 GeV, m=0.000511 GeV                                                                                         
[corsika:info    (ParticleCut.inl:101)] removing low en. particle... pid=e+, eKin=0.019826012093476522GeV, eThreshold\
=0.02GeV                                                                                                              
[corsika:debug   (ObservationPlane.inl:58)] writing particle to output. pid=e+, eKin(post)=0.019486999923203692GeV, e\
Kin(pre)=0.019826012093476522GeV                                                                                      
[corsika:debug   (Cascade.inl:210)] Cascade: delete absorbed particle PID=e+ E=0.020337012093476523 GeV

looks like the particles are written to disk just before they are deleted!

That's very interesting... Does this change when one uses a different order in the process sequence? I thought the results would be independent of the order.

And so far I don't understand why this would happen with SO many particles, so why do SO many particles reach the Observation plane if they are not correctly cut? I could understanding this for particles close to the plane, but not for particles that are just anywhere in the shower

We mucked up the calculation. Here is the code from getMaxStepLength in proposal::ContinuousProcess. We are mixing kinetic and total energies which leads to the particles being below the cut by their mass times 0.99999 (or something).

    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be
    // a hyper parameter which must be adjusted.
    //
    auto const energy = vP.getEnergy();
    auto const energy_lim =
        std::max(energy * 0.9, // either 10% relative loss max., or
                 get_kinetic_energy_propagation_threshold(
                     code) // energy thresholds globally defined for individual particles
                     * 0.9999 // need to go slightly below global e-cut to assure removal
                              // in ParticleCut. This does not matter since at cut-time
                              // the entire energy is removed.
        );

fixing that the electrons end up at 0.019998 GeV as expected.

Is this maybe even related to #563 (closed) ?

I don't think so. I think the problem is that the cut and observation plane are both continuous processes but the actual deleting of particles happens at the higher level (in Cascade). So in the case that the particle absorption happens during the final step to the observation plane (it may be a big step!) they will also get written to file because the writing to file happens also during the last step. The deleting happens AFTER the step.

Nope that can't be it either. What I described above is only responsible for the peak at the cut value... If we fix it we'd still get thousands of particles arriving at the observation level, at least according to the energy spectrum shown above.

Ooops

I calculated the height of the particles when they are written to file inside ObservationPlane.

 LengthType const z =
          (step.getPositionPost() - plane_.getCenter()).dot(plane_.getNormal())

The result you see above. Btw this is a 100GeV shower (for performance reasons). There no photons arriving at the ground and also no electrons! BUT somehow many (all?) electrons (massive particles?) end up in the output.

I was expecting this to be the case, otherwise I wouldn't have an explanation why there are more particles on the observation plane than in the highest bin of the longitudinal profile

Do we know when this started to be wrong? I would guess with the implementation of the new Step object?

Here is my brake down of the problem:

• Observation level is a Continuous process
• every Continuous process is called in every step.
• this means ObservationLevel needs to figure out if it is the limiting process, i.e. if the step ends at the ObservationLevel. Only in this case it writes out the particles.
• which process is limiting the step is determined in Cascade by calling getMaxStepLength.
• the id of the limiting process is then passed to all Continuous processes in doContinuous.
• inside ProcessSequence the limiting id is compared with the current Continuous process' id doContinuous(step, limitId == ContinuousProcessIndex(IndexProcess2))

In the case of multiple nested process sequences eg. SwitchProcess in emContinuous, these id's seem to get confused, be ill defined, or something.

In any case this makes ObservationLevel think it is Proposal::ContinuousProcess

I imagine the id's look something like this (p1,p2,(p3,p4),p4).

A workaround is changing the order of LongitudinalProfile and ObservationLevel in the process sequence or rather not make ObservationLevel follow the continuous energy loss (Proposal and Bethe).

With a local fix showers at 1 TeV and emcut of 20, 10 and 5 MeV run through without any problem. There are tons of warnings from the proposal interface about negative step lengths though. I think we have to remove the warning. Not sure if setting the step length to zero itself causes any problems.

Testing 500 keV now.. there are lot's of warnings from proposal about a negative MeanFreePath. [proposal.interaction] [warning] Negative MeanFreePath detected at energy 1.2408481412992356 MeV. Returning INF instead.

The 500keV shower finished after an eternity.
I like the charge excess. I also do get 11 photons at ground level but that seems ok.

Before it is forgotten. Apart from the Observation level we learned that we need to fix

• label of column 'energy' in particle file should be 'kinetic_energy', or change output to 'total energy'
• fix the calculation of the energy threshold in proposal interface get_thres.. --> get_thres.. + mass
• writing out the particles that were removed by the energy cut could be interesting in the future ...
• a counter of the particles that reach the observation level and that is printed out at the conclusion of the shower is really useful

created branch 573-observationplanes-show-non-existing-particles to address this issue

mentioned in merge request !490 (merged)

mentioned in issue #579 (closed)

mentioned in commit 5e685c89

closed with merge request !490 (merged)

mentioned in issue #586 (closed)