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

Adding o16nc interaction channel #44

Merged
merged 38 commits into from
Jan 20, 2023
Merged

Adding o16nc interaction channel #44

merged 38 commits into from
Jan 20, 2023

Conversation

SWakely
Copy link
Contributor

@SWakely SWakely commented Oct 19, 2022

o16nc_p.py and o16nc_n.py have been written to model neutral current interactions on O-16 (proton and neutron emission channels).
The cross-section as a function of energy has been fitted from data presented in (https://link.aps.org/doi/10.1103/PhysRevC.98.034613) and the threshold energy estimated from these fits.
The energy spectrum of the outgoing particle is assumed to be uniform as no information could be found for the correct shape of the spectrum. The outgoing energy is therefore generated at random from the available energy range.

SWakely and others added 24 commits August 24, 2022 15:07
This is the in progress code for neutral current interactions on O-16. The dominant neutral current interactions in the energy range of interest are those which result in a nucleon knockout (Suzuki et al. 2018). As such there are two code files, for neutron and proton knock-out (o16nc_n_code and o16nc_p_code respectively).
A spline function is used to produce a fit 'function' for cross-section as a function of neutrino energy. The data this spline uses is taken from Suzuki et al. 2018 with additional estimated data points (estimated from by-eye fits) to ensure a smooth spline fit. The threshold energies for the interactions have also been estimated from a by-eye fit as no theoretical or experimental research could be found into the exact values.
If the neutrino energy is large enough, the remnant nucleus (after nucleon knock-out) has ~70% chance of being left in an excited state (T2K 2019 paper). In this case two potentially detectable particles are produced, the knocked-out nucleon and a de-excitation photon.
A function, random(), generates a random number between 0 and 1 which is used to determine whether a de-excitation photon is produced. This has been written in the o16nc_n_code and will be added to the o16nc_p_code. I will probably adapt this function to a gamma(eNu) function which will return true if a de-excitation photon is produced.
A function to generate the kinematics of the knocked-out nucleon has been added to each code, however it has not been written in full yet.
…urns True if a gamma emission occurs. That is to say when the neutrino energy is larger than the threshold energy for leaving the remnant nucleus is the excited state and if a randomly generated number between 0-1 is larger than the probability of it being left in the ground state. The function gamma(eNu) will return True ~70% of the time when eNu > e_thr_g. The generate_event function has been adjusted to use the gamma(eNu) function instead of random().
…a de-excitation emission rather than a knock-out process. The neutrino excites the O-16 nucleus which de-excites via the emission of a nucleon. This emission may leave the daughter nucleus in an excited state. If this is the case the daughter nucleus then further de-excites via the emission of a gamma ray. Not considering the nucleon as being knocked out by the neutrino reconciles the apparent contradiction between Suzuki et al. 2018 and the T2K paper. Suzuki et al. 2018 shows nucleon final states (after NC neutrino O-16 interactions) as dominant in the 15-90 MeV range whereas the T2K paper states that the nucleon knock-out process only dominates above 200 MeV. If the Suzuki et al. 2018 data is for nucleon emission rather than knock-out, then the nucleon emission can dominate at low energies while knock-out only dominates at >200 MeV. Suzuki et al. don't explicitly say whether the nucleons are de-excitation emissions, however this would agree with the neutrino 'seeing' the nucleus as a whole in this energy range as the de Broglie wavelength is larger than a nucleon.

> This has some implications for the kinematics of the nucleon. If knocked-out, the direction of the nucleon's travel would depend on that of the neutrino, the kinematics would therefore be dealing with a three body final state (neutrino, nucleon and nucleus). However, if the neutrino excites the nucleus which subsequently emitts the nucleon, then the nucleon emission is isotropic (assuming spin effects of the nucleus average out over the whole of the detector).
> If a photon is produced as well, it's emission is also isotropic. As such both the nucleon and the photon are emitted in a random direction. A random (normalised) direction is given as an argument for the generate_event() function. I have assigned this to the nucleon independent of photon emission. I have therefore created a get_photon_direction() function which will generate a random (normalised) direction for the photon emission.
> The energy of the nucleon (eE) is chosen at random from the energy range 0 to [neutrino energy - threshold energy] where the threshold energy is changed to the threshold energy of gamma production if a gamma is produced.
…ls were correct in the relevant code and all other code was identical.

Improved commentary and docstrings for code.
…between ''' ''' relates to nucleon/photon emission and can be added again at a later point if desired. Re-reading Suzuki et al. makes it clear that the cross-sections given relate to nucleon ONLY emission.
… for given emission and commented out estimated data points for other emissions.
…o start of code. Threshold energy in data set to equal the value given prior. Should the threshold energy need changing it now only need be changed once in each document.
…nucleon from its momentum and mass. The numeric code for the interaction channel was changed to 2008016 and (-2008016). I could not find a reference for what this number should be however looking at the neutral current interaction channel for c12 and the charged current interaction channels for c12 and o16, it appeared that 200 signaled neutral current, while 8016 is the particle ID for the O16 nucleus.
… neutrino plus the rest mass of the nucleon.
…w returns the bounds of the delta function around eE
… energy required to produce a particle of given energy is e_thr + eE - (nucleon mass)
Copy link
Member

@JostMigenda JostMigenda left a comment

Choose a reason for hiding this comment

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

Many thanks for your contributions @SWakely! At first glance, this looks really good already. I have a few minor comments to clean the code up a bit; we can go over those at our meeting next week. I will take a closer look at the paper next week and double check the physics, too.

Before I merge this, it also needs unit tests and documentation. Let me know if you are able to work on those or if you’d like me to take over and write them!

sntools/detectors.py Outdated Show resolved Hide resolved
sntools/genevts.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
JostMigenda and others added 3 commits October 27, 2022 15:57
…ion channels with data of relevant emission channel only. The fit()(eNu) function currently returns an array of the cross-section value rather than just a number.
@JostMigenda JostMigenda linked an issue Nov 14, 2022 that may be closed by this pull request
@JostMigenda
Copy link
Member

Hi @SWakely! We talked about the random.random() call in get_eE(eNu, cosT) at our meeting recently and noticed that the current implementation would be a problem for bounds_eE since a different random value would be used for the lower and upper bound. (To fix that, we can call get_eE just once and then use that single value to calculate the lower and upper bounds.)
I’ve also said I would check if this could lead to problems elsewhere in the code. In the implementation in channel.py, there are a few occasions where get_eE() gets called indirectly; here are my thoughts about those:

  • In line 40 & line 69, via bounds_eE: This is then used as integration range for dSigma_dE(eNu, eE)—but since that doesn’t depend on eE, having a random value isn’t a problem here as long as the width of the integration range is two epsilon.
  • In line 71, via generate_event(): This is where the value that ends up in the output file is generated. Here, it’s totally fine to have a random value.

In summary: I think once bounds_eE is fixed, everything else should be alright.

Copy link
Member

@JostMigenda JostMigenda left a comment

Choose a reason for hiding this comment

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

I just looked through the implementation in more detailed; here are some changes I would suggest. They’re fairly straightforward and shouldn’t take much time. Hope it’s all clear, but let me know if I’ve misunderstood something!

(Note: All comments apply analogously to o16nc_p.py.)

sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
sntools/interaction_channels/o16nc_n.py Outdated Show resolved Hide resolved
SWakely and others added 9 commits December 3, 2022 13:15
…iable. eE = self.get_eE(eNu) is used in bounds_eE() and the return uses eE to generate the boundaries, this ensures the same value of eE is used for the two bounds. Unnecessary if statement removed from dSigma_dE() function. dSigma_dCosT() function adjusted to return 0.5*sigma. All these change have been made to both o16nc_p and o16nc_n code.
…iable. eE = self.get_eE(eNu) is used in bounds_eE() and the return uses eE to generate the boundaries, this ensures the same value of eE is used for the two bounds. Unnecessary if statement removed from dSigma_dE() function. dSigma_dCosT() function adjusted to return 0.5*sigma. All these change have been made to both o16nc_p and o16nc_n code.
These are not visible in the arXiv version due to a typesetting issue
@JostMigenda
Copy link
Member

Added a brief section in the documentation; the only thing that’s missing before this PR can be merged are unit tests for both channels. I’ll put that on my todo list for next week.

@JostMigenda
Copy link
Member

Added some basic tests for both subchannels, so this is now ready to merge. Once again, @SWakely, thanks for the contributions and congratulations on your first successful PR!

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.

NC interactions on oxygen
2 participants