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

Mott scattering #66

Open
wants to merge 17 commits into
base: main
Choose a base branch
from

Conversation

juniorpena
Copy link

This change is implementing Mott scattering of high energy electrons off of Helium nuclei. I uploaded the files from my Kassiopeia setup on my local machine into a repo I forked, and I made changes to the appropriate CMakeLists.txt files.

Copy link
Member

@zykure zykure left a comment

Choose a reason for hiding this comment

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

Thanks for the PR and for providing documentation with the code! I have some comments that need to be addressed before merging.

Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx Outdated Show resolved Hide resolved
Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx Outdated Show resolved Hide resolved
@zykure
Copy link
Member

zykure commented Jan 9, 2023

Is there any update on this PR @juniorpena ?

@juniorpena
Copy link
Author

juniorpena commented Feb 1, 2023 via email

@zykure zykure requested a review from 2xB February 2, 2023 08:19
@juniorpena
Copy link
Author

Hello, I've eliminated the use of ROOT in this scattering calculator. Now the GetTheta function uses inverse transform sampling and rejection sampling. I've also updated it to allow for both He and Ne Mott Scattering

@2xB
Copy link
Member

2xB commented May 7, 2024

Hey, thank you very much! Could you add documentation of the added objects to https://github.com/KATRIN-Experiment/Kassiopeia/blob/main/Kassiopeia/XML/Complete/Interactions.xml ? That is helpful both for users and for the review.

(Related, but independent issue: #106 )

@juniorpena
Copy link
Author

I've added the documentation

Copy link
Member

@2xB 2xB left a comment

Choose a reason for hiding this comment

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

After looking through the code, I had three remaining comments after which this looks good to merge!

Btw. if there will be a thesis describing a use case of this, one can also advertise it in the top comment of KSIntCalculatorMott.h, since it is always useful to see what the original use case for components of Kassiopeia was.

Comment on lines 69 to 74
std::vector<double> theta_arr;

// Populate theta_arr
for (int i = 0; i <= 1000; ++i) {
theta_arr.push_back(fThetaMin + i * (fThetaMax - fThetaMin) / 1000.0);
}
Copy link
Member

Choose a reason for hiding this comment

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

For a fixed-size array, std::array is better than std::vector. However, in this case, I think the array is not needed, as you just build it and directly afterwards use it - one could just calculate theta on the fly in the second loop. Also, it makes sense to document in a comment where the 1000 comes from and when it is valid?

Copy link
Author

Choose a reason for hiding this comment

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

It was arbitrarily chosen, but now there is a variable to control the size of the array and I added a comment about how it affects the simulation time

Copy link
Member

Choose a reason for hiding this comment

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

Thanks! I still believe one could avoid arrays/vectors completely along the lines of

#include <limits>

[...]

    // Find largest rescaling factor for rejection sampling over all angles
    double k = -std::numeric_limits<double>::infinity();

    for (int i = 0; i <= int(resolution); ++i) {
        const double current_theta = fThetaMin + i * (fThetaMax - fThetaMin) / resolution;
        const double current_k     = MDCS(theta, anEnergy) / Normalized_RDCS(theta);

        if (k < current_k) {
            k = current_k;
        }
    }

, but I'm not sure if that is relevant performance- or memory-wise at all.

Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx Outdated Show resolved Hide resolved
juniorpena and others added 2 commits June 10, 2024 11:23
…ilder.h

Co-authored-by: 2xB <31772910+2xB@users.noreply.github.com>
I've added code and comments for the choice of size of theta_arr responding to previous comment. I've also added comments wherever there are hard coded values, for where those values come from
@2xB
Copy link
Member

2xB commented Jul 4, 2024

Our Ubuntu 20.04 CI check failed with

/__w/Kassiopeia/Kassiopeia/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx:10:7: error: use of undeclared identifier 'KGeoBag'
using KGeoBag::KThreeVector;
      ^
/__w/Kassiopeia/Kassiopeia/Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx:43:5: error: unknown type name 'KThreeVector'; did you mean 'katrin::KThreeVector'?
    KThreeVector tPosition = anInitialParticle.GetPosition();
    ^~~~~~~~~~~~
    katrin::KThreeVector

, as can be seen by clicking on the check details...

@2xB
Copy link
Member

2xB commented Jul 4, 2024

Btw. additionally to fixing the above error, I believe that it would make sense to update your branch by merging our main branch into it, that would fix the fedora_37 test case.

Merged main branch, and fixed bug causing Ubuntu 20.04 CI check failure
Found a bug after running simulations, the momentum of the final state particle was being changed twice when executing the interaction
if (fNucleus == "He") {
amu = 4.0026; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002
} else if (fNucleus == "Ne") {
amu = 19.9924; // mass in amu from W.J. Huang et al 2021 Chinese Phys. C 45 030002
Copy link
Member

Choose a reason for hiding this comment

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

This is for 20Ne and the above for 4He, right? Is your entire code specific to those isotopes?

Kassiopeia/Interactions/Source/KSIntCalculatorMott.cxx Outdated Show resolved Hide resolved
@2xB
Copy link
Member

2xB commented Aug 26, 2024

@juniorpena What is the status here? I have added one question to resolve one of the two open threads here, the other I think is not that important but I leave it open just in case you want to comment on it as well. I think once at least the thread regarding the provided isotopes is solved, we should finally be able to merge this from my side. (CC @richeldichel )

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants