Error when setting low particle cuts for em showers
An error that terminates the shower occurs when setting low particle cuts (10 MeV) and electromagnetic interaction. The error is not present when setting high cuts (10 GeV). The error message received was:
terminate called after throwing an instance of 'boost::wrapexcept<boost::math::evaluation_error>'
what(): Error in function boost::math::tools::bisect<double>: No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = -9.8713280286356291).
Aborted (core dumped)Previously several messages of:
No node was found where the function value was greater than zero. Either there are no continuous losses for the cut, or something went wrong.
coming from crosssection/CrossSectionDEDX/AxisBuilderDEDX.h are shown.
The proposal_new_installation_and_bugfixes branch was used (latest commit 60f00801) and PROPOSAL 7.0.2 is running.
I leave here the examples used and the ouput to reproduce it. Vertical EAS (modified to start with an electron and have electromagnetic interaction) was run with ./bin/vertical_EAS 0 1 1000 1 and em shower with ./bin/em_shower 1000. em_shower.cpp vertical_EAS.cpp output_em_shower output_vertical_EAS
Activity
added EM label
I think that issue was caused pulling the current main branch here. I didn't notice it because the CI is currently not working and i haven't checked each line
The energy above which the particles are considered stochastic is set to the particle mass. This does not work for leptons because the crosssection divergence for small energies. Instead, a value in the order of the ParticleCut is usually a reasonable choice.
I have overlooked the pull request !303 (merged). (Sorry @riehn) I think @JAlameddine and I have to think about this. The quick fix would be to set the threshold for each electromagnetic particle with
set_energy_thresholdto an energy near the particle cut or simply checkout the previous commit.I will have another look at it later this afternoon.
Edited by Maximilian Sackeladded Bug Critical In Progress labels
It took longer than expected, but I could back trace the error. If I run the std em_shower example with an initial energy from 1000 GeV and a particle cut from 10 MeV I could reproduce the failure. If I back trace where it comes from and log the geomMaxLength I get a small negative value (have a look here). I think it comes from the leap frog algorithm. A negative distance is not well-defined in PROPOSAL and the error is thrown.
I'm not sure if a small negative distance is expected or maybe a bug? @rulrich I could introduce a geometric imprecision, but then particles that travel a distance less than the imprecision would no longer lose energy via the electromagnetic process. I would like to avoid this if possible, because otherwise they lose systematically too little energy. Maybe it's related to the little stochastic cut, but I'm not sure jet. The particle only make very small step sizes because of the
get_energy_thresholdinitialization.I am considering introducing a particle related variable that can set an energy per type, above which losses are treated stochastically. This prevents that with very high energetic particle cuts, the step size is extremely reduced and a large part of the particles is immediately deleted from the stack by the cut. Any opinions about this?
Edited by Maximilian Sackel
About the stepsizes: negative steps should not appear. This is hopefully fixed in !271 (merged) . Could you make a test if the problem is fixed in !271 (merged)? This would be great.
This MR fixes the issue mentioned here also. Should I pull it into this branch, or should I wait till it's on the main?
Another question: Since I have not found it in a hurry I wanted to ask how to set the logging in corsika actually. A reference to the appropriate place is enough for me.
Edited by Maximilian Sackel-
For the MR: don't pull this in. It is a fix for something specific and it is sufficient to have both independently merged on master during the next days (asap).
For the logging:
https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/wikis/Logging-Guidelines
and good example for how-to log is:
https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/blob/master/corsika/modules/OnShellCheck.hpp
Doing the merge with 280-assumption-of-no-continuous-between-two-stochastic-losses solved the issue for 1 TeV showers, of 50 not one failed. When doing 100 TeV showers only 11 finished. Most of them with the previous error
Error in function boost::math::tools::bisect<double>: No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the intervalbut the error
Error in function boost::math::tools::newton_raphson_iterate<double>: There appears to be no root to be found in boost::math::tools::newton_raphson_iterate, perhaps we have a local minima near current best guess of 0.970703125also showed up. So it seems to solve it only partially.
I will have a closer look as soon as !271 (merged) has been merged.
I am still getting these boost errors with these settings (100 TeV + 10 MeV cut), although these are rather rare cases (compared to what Juan has been experiencing). I will investigate this..
Edited by Jean-Marco AlameddineOk, now I am actually unable to reproduce the runtime errors from boost. I've been able to propagate ~30 100 TeV showers (10 MeV PaticleCut) without any problems, and even some 1000 TeV showers.
The only problem I encountered were some other Arithmetic Exceptions thrown by PROPOSAL, but I introduced a quick fix in PROPOSAL which will be included in the next PROPOSAL release (version 7.0.5). Until proposal/7.0.5 is included in corsika, one can simply disable the floating point exceptions if these errors occur.
assigned to @rulrich
assigned to @JAlameddine and unassigned @rulrich
@JAlameddine your last comment sounds like the issue is solved. If you can confirm this please close the issue.
