-
Notifications
You must be signed in to change notification settings - Fork 29
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fix hydrogen ionization scattering #124
base: main
Are you sure you want to change the base?
Conversation
Remove unrelated/unnecessary previous commit according to 2xB's suggestion. Co-authored-by: 2xB <[email protected]>
aReducedFinalEnergy = tReducedFinalEnergy; | ||
|
||
// convert to energy of primary by subtracting binding energy and secondary energy from input | ||
aReducedFinalEnergy = aReducedInitalEnergy-1.-tReducedFinalEnergy; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why the "-1"? And what is e_hlf, it also contains this "-1"? I see that the -1 also occurs at more places, but why?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Several Points:
-
the more efficient sampling is great
-
the angular distribution makes more sense of course
-
i vaguely remember having sampled only half the function intentionally. it is fundamentally symmetric because primary and secondary are indistinguishable, but basically you could just call the lower energy one the secondary. it's also a question of normalisation. if you need more input on that i have to take a deeper look.
-
the -1 should be the binding energy in the appropriate units used in the rudd paper.
cheers
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Several Points:
-
the more efficient sampling is great
-
the angular distribution makes more sense of course
-
i vaguely remember having sampled only half the function intentionally. it is fundamentally symmetric because primary and secondary are indistinguishable, but basically you could just call the lower energy one the secondary. it's also a question of normalisation. if you need more input on that i have to take a deeper look.
-
the -1 should be the binding energy in the appropriate units used in the rudd paper.
cheers
Modifications were made to correct mistakes in the hydrogen scattering implementation concerning ionization. Affected functions:
KSIntCalculatorHydrogenIonisation::ExecuteInteraction
linkKSIntCalculatorHydrogenIonisation::CalculateFinalEnergy
linkFixed issues:
1) Sampling from the singly differenctial cross section (SDCS), i.e. the ionization eloss function:
The sampling method originally implemented in

KSIntCalculatorHydrogenIonisation::CalculateFinalEnergy
to sample from the SDCS was erroneous. Firsty, when all outgoing electrons after ionization are taken into account (outgoing primary and secondary), the histogram of sampled energy losses should extend all the way down to zero. This was not the case in the previous implementation. Furthermore, it seems it was attempted to implement rejection-based sampling with a (inefficient) uniform proposal distribution. However, there are some mistakes and the distribution of sampled energies therfore does not follow the SDCS as intended.These issues are adressed by re-writing most parts of the function. For the new implementation, several new proposal functions for the rejection sampling were investigated, in order to make the new implementation much more efficient. An ideal proposal function must be as close to the target distribution, and one must ideally be able to create samples via inversion sampling. Since the SDCS is symmetrical around the middle, it suffices to define a proposal for one half. The second half is then autmatically taken into account, by generating the secondary with the energy E_secondary = E_input-E_binding-E_primary. The investigated proposal distributions are:

Out of these proposals, the custom powerlaw is closest and should be most efficient. A test in a standalone python implementation of the methods was used to assess the acceptance efficiency and the target trueness:
Target trueness:

Number of trials for successful sample generation:

Acceptance ratios:
Acceptance ratio (Kassiopeia sampling method) : 00.123 % 1e5 samples in 08:08 minutes
Acceptance ratio (Uniform rejection samping) : 00.121 % 1e5 samples in 08:01 minutes
Acceptance ratio (Exponential rejection samping) : 01.824 % 1e6 samples in 07:06 minutes
Acceptance ratio (Powerlaw rejection samping) : 53.575 % 1e6 samples in 00:19 minutes
To conclude, the custom powerlaw proposal distribution is true to the target and performs orders of magnitude better than the other methods, especially much better than the original implementation. Therefore, the powerlaw was chosen for the new Kassiopeia implementation.
2) Secondary electron angle
In the original implementation of
KSIntCalculatorHydrogenIonisation::ExecuteInteraction
, the phi and theta angles of the secondary electron from ionization were simply drawn from uniform distributions. However, due to the kinematics (momentum conservation) there should be a correlation between the scattering angles of the outgoing primary and secondary. For low electron energies this may be difficult to describe accurately. But for energies much larger than the binding energy (which is usally the case), one can assume free electron-electron scattering.Therefore, in the new implementation, the momentum direction (and thus the phi and theta angles) is simply computed by subtracting the momentum of the outgoing primary from the momentum of the incoming electron.
Both changes have been tested and are used by Julanan Songwadhana (SUT) for his KATRIN WGTS magnetic trap and source scattering simulations. No issues with the implementation were found so far.