Skip to content
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

Add Wentzel macro xs calculator and fix a_sq_factor #1274

Merged
merged 10 commits into from
Jun 22, 2024

Conversation

amandalund
Copy link
Contributor

@amandalund amandalund commented Jun 14, 2024

  • Fix an error in the calculation of a_sq_factor for Wentzel OK&VI
  • Fix calculation of inv_mass_cbrt_sq
  • Fix a_sq_factor units
  • Add a Wentzel total xs calculator (to be used in the Wentzel VI msc model)
  • Move some of the common code for the Urban and Wentzel VI tests into a msc test base.

@amandalund amandalund added bug Something isn't working enhancement New feature or request physics Particles, processes, and stepping algorithms labels Jun 14, 2024
@amandalund amandalund requested review from sethrj and whokion June 14, 2024 00:28
@amandalund amandalund force-pushed the wentzel-xs-calculator branch from 0fbb7dc to 8a51eca Compare June 14, 2024 00:29
Copy link
Member

@sethrj sethrj left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM, but I'll let @whokion take a look before merging. I'm going to remove the bug label since this is mostly an enhancement and the bug hasn't been seen in non-testing use cases.

@@ -36,7 +36,7 @@ struct CoulombParameters
//! Factor used to calculate the maximum scattering angle off of a nucleus
real_type a_sq_factor{real_type(0.5)
* ipow<2>(constants::hbar_planck * constants::c_light
* units::femtometer)};
/ units::femtometer)};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe a Quantity (FmSq?) would help elucidate this factor? The divide vs multiply has bitten me several times...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would also be grateful if someone (maybe @whokion?) could take a second look at this factor and where it's used (in geant4 and celeritas) and make sure everything looks ok. I'm still not really sure where it comes from or if I have all the units right...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mo problem, I will take a second look for this today.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, it looks all good to me. This factor is related tothe screening parameter, $A \equiv \frac{1}{4}(\frac{\hbar}{pR})^2$ (where $p$ is the particle momentum) which is originated from the Wentzel SS differential cross section based on the first Born approximation using a simple potential $V(r) \sim \frac{exp(-r/R)}{r}$ where $R \sim 0.885 Z^{1/3} a_o$ is the screening radius in terms of the Bohr radius $a_o$ and the atomic (charge) number $Z$ (i.e, based on Thomas-Fermi model) and then is also a parameter of the angular distribution of the single scattering - please see eq 23, 26 of the paper by J.M. Fernrindez-Varea et al.. When the SS is combined with the WentzelVI multiple scattering, the maximum angle of the scattering is therefore limited by $1-A$ (i.e., the line that you pointed in Geant4).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the explanation and reference! That makes more sense now. I'll take another look and see if I can understand the differences between the screening parameter and the value that's used in the code to limit the angle.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok I fixed another bug (in the calculation of the effective atomic mass^(-2/3)), but the quantity that we're subtracting here is still about 4 orders of magnitude smaller than the screening parameter. Should they be roughly the same? The screening parameter itself is also very small (like O(1e-10)) which means the maximum angle is going to be limited to near zero instead of the polar angle limit... so I think I'm still missing something.

Copy link
Contributor

@whokion whokion Jun 14, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I guess that the magnitude, 1-factorA^2/ [(mass number)^{2/3} * p^2] should be close to 1, but not too close to 1 - based on my rough calculation, it is about 0.95/0.98 for electrons of 100/200 MeV (kinetic) energy on Pb (A~200).

@@ -77,7 +77,7 @@ WentzelOKVIParams::WentzelOKVIParams(SPConstMaterials materials,
host_data.params.a_sq_factor
= real_type(0.5)
* ipow<2>(options.angle_limit_factor * constants::hbar_planck
* constants::c_light * units::femtometer);
* constants::c_light / units::femtometer);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the factor's always being set here (it looks like??) maybe just default the CoulombParameters to zero and reduce the duplication here.

@sethrj sethrj removed the bug Something isn't working label Jun 14, 2024
Copy link
Contributor

@whokion whokion left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM too.

@amandalund
Copy link
Contributor Author

One more fix... I think a_sq_factor should have units MeV^2. The cos(theta) limit looks more reasonable now.

Copy link
Member

@sethrj sethrj left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good work!

@sethrj sethrj merged commit c87c167 into celeritas-project:develop Jun 22, 2024
29 checks passed
@amandalund amandalund deleted the wentzel-xs-calculator branch June 22, 2024 21:35
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request physics Particles, processes, and stepping algorithms
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants