From f8c38452c360bdceb71c0d258ba42c2f18ecf769 Mon Sep 17 00:00:00 2001 From: Jost Migenda Date: Sat, 12 Aug 2023 22:13:25 +0100 Subject: [PATCH 1/2] Revert "Remove o16nc support for this release tag" This reverts commit 7111bc7f72e15c4f22a84dc33c1df7135f460bea. --- doc/documentation.tex | 19 ++++ src/sntools/detectors.py | 6 +- src/sntools/genevts.py | 3 +- src/sntools/interaction_channels/o16nc_n.py | 112 ++++++++++++++++++++ src/sntools/interaction_channels/o16nc_p.py | 112 ++++++++++++++++++++ tests/test_o16nc_n.py | 51 +++++++++ tests/test_o16nc_p.py | 51 +++++++++ 7 files changed, 352 insertions(+), 2 deletions(-) create mode 100644 src/sntools/interaction_channels/o16nc_n.py create mode 100644 src/sntools/interaction_channels/o16nc_p.py create mode 100644 tests/test_o16nc_n.py create mode 100644 tests/test_o16nc_p.py diff --git a/doc/documentation.tex b/doc/documentation.tex index b56824a..317356d 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -220,6 +220,25 @@ \subsubsection{\texttt{o16e} and \texttt{o16eb}: Charged-Current Interactions on For a typical supernova neutrino flux, the difference in the resulting event spectra when using the four groups instead of all 42 nuclear states is very small. +\subsubsection{\texttt{o16nc\_n} and \texttt{o16nc\_p}: Neutral-Current Interactions on $^{16}$O} +In water, neutral-current interactions on $^{16}$O nuclei are a subdominant interaction channel that neutrinos and antineutrinos of all flavours contribute to equally. +The largest contribution to this channel comes from events where a single neutron or proton\footnote{ + Since the energy of the emitted proton is below the Cherenkov threshold, it cannot be detected in water Cherenkov detectors. + Therefore, sntools generates \texttt{o16nc\_p} events only for WbLS detectors. + However, \texttt{o16nc\_n} events \textit{are} generated also in water Cherenkov detectors, since the neutron capture signal is potentially detectable. +} is emitted, +\begin{align} +\nu + \/^{16}\text{O} \rightarrow \nu' + &^{16}\text{O}^*\\ + \Rightarrow &^{16}\text{O}^{*} \rightarrow X + n\\ +\nu + \/^{16}\text{O} \rightarrow \nu' + &^{16}\text{O}^*\\ + \Rightarrow &^{16}\text{O}^{*} \rightarrow X + p. +\end{align} + +In sntools, the cross sections for these interactions is based on the theoretical calculations tabulated in~\cite{Suzuki2018}. +The energy distribution of the outgoing particle is unknown. +In the current implementation, the energy is simply taken to be a random fraction of the available neutrino energy above threshold. + + \subsubsection{\texttt{c12e} and \texttt{c12eb}: Charged-Current Interactions on $^{12}$C} In liquid scintillator, charged-current interactions of \nue and \nuebar on $^{12}$C nuclei, \begin{align} diff --git a/src/sntools/detectors.py b/src/sntools/detectors.py index c191380..96f55f9 100644 --- a/src/sntools/detectors.py +++ b/src/sntools/detectors.py @@ -6,7 +6,7 @@ water = { "molecular_weight": 18.0153, # g/mol "density": 1.0, # g/cm^3 - "channel_weights": {"ibd": 2, "es": 10, "o16e": 1, "o16eb": 1}, # targets per molecule + "channel_weights": {"ibd": 2, "es": 10, "o16e": 1, "o16eb": 1, "o16nc_n": 1}, # targets per molecule } # liquid scintillator: approximated here as CH_2 @@ -39,6 +39,10 @@ def wbls(x): if not 0 <= x <= 1: raise ValueError("Fraction of Liquid Scintillator must be between 0 and 1!") + # The free proton from o16nc_p is undetectable in pure water, so not included above. + # Since it may be detectable in WbLS, we add it manually here. + water["channel_weights"]["o16nc_p"] = 1 + mw = x * ls["molecular_weight"] + (1 - x) * water["molecular_weight"] d = x * ls["density"] + (1 - x) * water["density"] cw = {} diff --git a/src/sntools/genevts.py b/src/sntools/genevts.py index 2e5cbce..6ba3c86 100644 --- a/src/sntools/genevts.py +++ b/src/sntools/genevts.py @@ -97,10 +97,11 @@ def parse_command_line_options(): help="Transformation between neutrino flux inside SN and flux in the detector on Earth. \ Choices: %(choices)s. Default: %(default)s.") - choices = ("ibd", "es", "ps", "o16e", "o16eb", "c12e", "c12eb", "c12nc") + choices = ("ibd", "es", "ps", "o16e", "o16eb", "o16nc_n", "o16nc_p", "c12e", "c12eb", "c12nc") parser.add_argument("-c", "--channel", metavar="INTCHANNEL", choices=choices, default="all", help="Interaction channels to consider. Currently supported: inverse beta decay (ibd), \ electron scattering (es), proton scattering (ps), nu_e + oxygen CC (o16e), nu_e-bar + oxygen CC (o16eb), \ + nu + oxygen NC with neutron emission (o16nc_n) or proton emission (o16nc_p), \ nu_e + carbon CC (c12e), nu_e-bar + carbon CC (c12eb) and nu + carbon NC (c12nc). \ Default: All channels available in selected detector.") diff --git a/src/sntools/interaction_channels/o16nc_n.py b/src/sntools/interaction_channels/o16nc_n.py new file mode 100644 index 0000000..6d34269 --- /dev/null +++ b/src/sntools/interaction_channels/o16nc_n.py @@ -0,0 +1,112 @@ +"""Implementation of: + nu + 16O -> nu' + 16O*, 16O* -> 15O + n. + +Based on data provided in Suzuki et al. 2018 (Phys. Rev. C 98, 034613). +A spline fit is used to obtain cross-section as a function of neutrino energy. +The threshold energy of the interaction is estimated from a by-eye fit as no +explicit threshold energy could be found. + +To determine the the differential cross section dSigma/dE (eNu, eE) from the +total cross section, we approximate a DiracDelta function with one that is +2*epsilon wide and 1/(2*epsilon) high, so that the integral is 1. +""" + +from sntools.event import Event +from sntools.interaction_channels import BaseChannel +from scipy.interpolate import interp1d +import random + +e_thr = 19.5 # approximate energy threshold of neutron emission +mN = 939.7 # neutron mass (MeV) +epsilon = 0.001 # for approximating DiracDelta distribution below + +# List of neutrino flavors ("e", "eb", "x", "xb") that interact in this channel. +possible_flavors = ("e", "eb", "x", "xb") + +# Energies (MeV) and cross-sections (10^-42 cm^2) for o16nc neutron emission taken from Suzuki et al 2018. +data = [[e_thr, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 80.0, 90.0, 100], + [0.0, 0.000193, 0.0163, 0.129, 0.48, 1.22, 2.49, 4.44, 7.17, 10.7, 15.2, 20.4, 32.8, 46.8, 61.1]] + +# Spline fit of the partial cross-section as a function of energy +fit = interp1d(data[0], data[1], kind='cubic', fill_value='extrapolate', bounds_error = False) + + +class Channel(BaseChannel): + def generate_event(self, eNu, dirx, diry, dirz): + """Return an event with the appropriate incoming/outgoing particles. + + Input: + eNu: neutrino energy + dirx, diry, dirz: direction of outgoing particle (normalized to 1) + """ + # Note: `self.flavor` is set during __init__ + nu_flv = {'e': 12, 'eb': -12, 'x': 14, 'xb': -14}[self.flavor] + + eE = self.get_eE(eNu, dirz) + + evt = Event(2008016 if nu_flv > 0 else -2008016) + evt.incoming_particles.append([nu_flv, eNu, 0, 0, 1]) # incoming nu + evt.incoming_particles.append((8016, 14900, 0, 0, 1)) # oxygen-16 nucleus at rest + evt.outgoing_particles.append([2112, eE, dirx, diry, dirz]) # emitted neutron + # evt.outgoing_particles.append([nu_flv, eNu-e_thr, 0, 0, 1]) # outgoing nu + return evt + + def bounds_eE(self, eNu, *args): + """Return kinematic bounds for integration over eE. + + Input: + eNu: neutrino energy (in MeV) + args: [ignore this] + Output: + list with minimum & maximum allowed energy of outgoing (detected) particle + """ + eE = self.get_eE(eNu) + return [eE - epsilon, eE + epsilon] + + def get_eE(self, eNu, cosT=0): + """Return energy (in MeV) of outgoing (detected) particle. + + Input: + eNu: neutrino energy (in MeV) + cosT: cosine of the angle between neutrino and outgoing (detected) particle + """ + eE = random.random()*(eNu-e_thr) + mN + return eE + + def dSigma_dE(self, eNu, eE): + """Return differential cross section in MeV^-2. + + Inputs: + eNu: neutrino energy + eE: energy of outgoing (detected) particle + """ + if eNu < e_thr: + # This should never happen, since we set bounds for eNu accordingly above + # ... but just in case: + return 0 + + sigma = fit(eNu)*10**(-42) # cross-section (in cm^2) at eNu from the fit of Suzuki et al. 2018 data + sigma *= (5.067731E10)**2 # convert cm^2 to MeV^-2: http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1) + return sigma / (2 * epsilon) # Ensure that integration over eE yields sigma + + def dSigma_dCosT(self, eNu, cosT): + """Return differential cross section in MeV^-2 as a function of the + emission angle of the outgoing (detected) particle. + + Input: + eNu: neutrino energy (MeV) + cosT: cosine of the angle between neutrino and outgoing (detected) particle + """ + # Energy dependence is unclear, so we use a constant value for now. + if abs(cosT) > 1: + return 0 + sigma = fit(eNu)*10**(-42) # cross-section (in cm^2) at eNu from the fit of Suzuki et al. 2018 data + sigma *= (5.067731E10)**2 # convert cm^2 to MeV^-2: http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1) + return 0.5*sigma + + # List with minimum & maximum energy of incoming neutrino. + bounds_eNu = (e_thr, 100) + + def _bounds_eNu(self, eE): + """Min/max neutrino energy that can produce a given neutron energy.""" + return (e_thr + eE - mN, 100) \ No newline at end of file diff --git a/src/sntools/interaction_channels/o16nc_p.py b/src/sntools/interaction_channels/o16nc_p.py new file mode 100644 index 0000000..f42996b --- /dev/null +++ b/src/sntools/interaction_channels/o16nc_p.py @@ -0,0 +1,112 @@ +"""Implementation of: + nu + 16O -> nu' + 16O*, 16O* -> 15N + p. + +Based on data provided in Suzuki et al. 2018 (Phys. Rev. C 98, 034613). +A spline fit is used to obtain cross-section as a function of neutrino energy. +The threshold energy of the interaction is estimated from a by-eye fit as no +explicit threshold energy could be found. + +To determine the the differential cross section dSigma/dE (eNu, eE) from the +total cross section, we approximate a DiracDelta function with one that is +2*epsilon wide and 1/(2*epsilon) high, so that the integral is 1. +""" + +from sntools.event import Event +from sntools.interaction_channels import BaseChannel +from scipy.interpolate import interp1d +import random + +e_thr = 14.0 # approximate energy threshold of proton emission +mP = 938.3 # proton mass (MeV) +epsilon = 0.001 # for approximating DiracDelta distribution below + +# List of neutrino flavors ("e", "eb", "x", "xb") that interact in this channel. +possible_flavors = ("e", "eb", "x", "xb") + +# Energies (MeV) and cross-sections (10^-42 cm^2) for o16nc proton emission taken from Suzuki et al 2018. +data = [[e_thr, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 80.0, 90.0, 100], + [0.0, 0.000605, 0.0308, 0.198, 0.776, 2.16, 4.76, 9.01, 15.2, 23.7, 34.5, 47.7, 62.9, 98.6, 138, 179]] + +# Spline fit of the partial cross-section as a function of energy +fit = interp1d(data[0], data[1], kind='cubic', fill_value='extrapolate', bounds_error = False) + + +class Channel(BaseChannel): + def generate_event(self, eNu, dirx, diry, dirz): + """Return an event with the appropriate incoming/outgoing particles. + + Input: + eNu: neutrino energy + dirx, diry, dirz: direction of outgoing particle (normalized to 1) + """ + # Note: `self.flavor` is set during __init__ + nu_flv = {'e': 12, 'eb': -12, 'x': 14, 'xb': -14}[self.flavor] + + eE = self.get_eE(eNu, dirz) + + evt = Event(2008016 if nu_flv > 0 else -2008016) + evt.incoming_particles.append([nu_flv, eNu, 0, 0, 1]) # incoming nu + evt.incoming_particles.append((8016, 14900, 0, 0, 1)) # oxygen-16 nucleus at rest + evt.outgoing_particles.append([2212, eE, dirx, diry, dirz]) # emitted proton + # evt.outgoing_particles.append([nu_flv, eNu-e_thr, 0, 0, 1]) # outgoing nu + return evt + + def bounds_eE(self, eNu, *args): + """Return kinematic bounds for integration over eE. + + Input: + eNu: neutrino energy (in MeV) + args: [ignore this] + Output: + list with minimum & maximum allowed energy of outgoing (detected) particle + """ + eE = self.get_eE(eNu) + return [eE - epsilon, eE + epsilon] + + def get_eE(self, eNu, cosT=0): + """Return energy (in MeV) of outgoing (detected) particle. + + Input: + eNu: neutrino energy (in MeV) + cosT: cosine of the angle between neutrino and outgoing (detected) particle + """ + eE = random.random()*(eNu-e_thr) + mP + return eE + + def dSigma_dE(self, eNu, eE): + """Return differential cross section in MeV^-2. + + Inputs: + eNu: neutrino energy + eE: energy of outgoing (detected) particle + """ + if eNu < e_thr: + # This should never happen, since we set bounds for eE and eNu accordingly above + # ... but just in case: + return 0 + + sigma = fit(eNu)*10**(-42) # cross-section (in cm^2) at eNu from the fit of Suzuki et al. 2018 data + sigma *= (5.067731E10)**2 # convert cm^2 to MeV^-2: http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1) + return sigma / (2 * epsilon) # Ensure that integration over eE yields sigma + + def dSigma_dCosT(self, eNu, cosT): + """Return differential cross section in MeV^-2 as a function of the + emission angle of the outgoing (detected) particle. + + Input: + eNu: neutrino energy (MeV) + cosT: cosine of the angle between neutrino and outgoing (detected) particle + """ + # Energy dependence is unclear, so we use a constant value for now. + if abs(cosT) > 1: + return 0 + sigma = fit(eNu)*10**(-42) # cross-section (in cm^2) at eNu from the fit of Suzuki et al. 2018 data + sigma *= (5.067731E10)**2 # convert cm^2 to MeV^-2: http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1) + return 0.5*sigma + + # List with minimum & maximum energy of incoming neutrino. + bounds_eNu = (e_thr, 100) + + def _bounds_eNu(self, eE): + """Min/max neutrino energy that can produce a given proton energy.""" + return (e_thr + eE - mP, 100) \ No newline at end of file diff --git a/tests/test_o16nc_n.py b/tests/test_o16nc_n.py new file mode 100644 index 0000000..ab59e9e --- /dev/null +++ b/tests/test_o16nc_n.py @@ -0,0 +1,51 @@ +import unittest + +from sntools.interaction_channels import o16nc_n +from ._crosssectiontest import CrossSectionTest + + +class O16ETest(CrossSectionTest): + c = o16nc_n.Channel('e') # ensure we can access interaction channel module as self.c + + # iterable with tuples (eNu, eE, dSigma_dE(eNu, eE)) + test_dSigma_dE_values = ( + (20, 0, 2.478303107626836e-22), + (20, 10, 2.478303107626836e-22), + (30, 10, 1.6564823879992844e-19), + (50, 10, 5.7013812424161415e-18), + (75, 10, 3.37836296864102e-17), + (100, 10, 7.845819682694285e-17), + ) + + # iterable with tuples (eNu, eE, dSigma_dE(eNu, eE)) + test_dSigma_dE_edgecases_values = ( + (19.49, 0, 0), # eNu too small + # dSigma_dE(eNu, eE) is independent of eE, so no edge cases for eE exist + ) + + # iterable with tuples (eNu, cosT, dSigma_dE(eNu, cosT)) + test_dSigma_dCosT_values = ( + (20, -0.99, 2.478303107626836e-25), + (20, 0, 2.478303107626836e-25), + (20, 0.99, 2.478303107626836e-25), + ) + + # iterable with tuples (eNu, cosT, get_eE(eNu, cosT)) + test_get_eE_values = ( + # cannot test this, since return value is uncertain + ) + + # iterable with tuples (eNu, bounds_eE(eNu)[0], bounds_eE(eNu)[1]) + test_bounds_eE_values = ( + # cannot test this, since return value is uncertain + ) + + # value of bounds_eNu[0] + test_bounds_eNu_minvalue = 19.5 + + +# ensure that unittest doesn't run tests in the base class, via https://stackoverflow.com/a/22836015 +del CrossSectionTest + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_o16nc_p.py b/tests/test_o16nc_p.py new file mode 100644 index 0000000..3d1f896 --- /dev/null +++ b/tests/test_o16nc_p.py @@ -0,0 +1,51 @@ +import unittest + +from sntools.interaction_channels import o16nc_p +from ._crosssectiontest import CrossSectionTest + + +class O16ETest(CrossSectionTest): + c = o16nc_p.Channel('e') # ensure we can access interaction channel module as self.c + + # iterable with tuples (eNu, eE, dSigma_dE(eNu, eE)) + test_dSigma_dE_values = ( + (20, 0, 3.955012213207594e-20), + (20, 10, 3.955012213207594e-20), + (30, 10, 9.96457622548407e-19), + (50, 10, 1.951824209115436e-17), + (75, 10, 1.0277407236510741e-16), + (100, 10, 2.2985298252083094e-16), + ) + + # iterable with tuples (eNu, eE, dSigma_dE(eNu, eE)) + test_dSigma_dE_edgecases_values = ( + (13.99, 0, 0), # eNu too small + # dSigma_dE(eNu, eE) is independent of eE, so no edge cases for eE exist + ) + + # iterable with tuples (eNu, cosT, dSigma_dE(eNu, cosT)) + test_dSigma_dCosT_values = ( + (20, -0.99, 3.955012213207594e-23), + (20, 0, 3.955012213207594e-23), + (20, 0.99, 3.955012213207594e-23), + ) + + # iterable with tuples (eNu, cosT, get_eE(eNu, cosT)) + test_get_eE_values = ( + # cannot test this, since return value is uncertain + ) + + # iterable with tuples (eNu, bounds_eE(eNu)[0], bounds_eE(eNu)[1]) + test_bounds_eE_values = ( + # cannot test this, since return value is uncertain + ) + + # value of bounds_eNu[0] + test_bounds_eNu_minvalue = 14.0 + + +# ensure that unittest doesn't run tests in the base class, via https://stackoverflow.com/a/22836015 +del CrossSectionTest + +if __name__ == "__main__": + unittest.main() From 9292e95f41647d04fd5e9f099b3d9c2c762fe4af Mon Sep 17 00:00:00 2001 From: Jost Migenda Date: Sat, 12 Aug 2023 22:24:03 +0100 Subject: [PATCH 2/2] =?UTF-8?q?Don=E2=80=99t=20list=20all=20channels=20exp?= =?UTF-8?q?licitly=20in=20introduction=20text=20Ensures=20this=20sentence?= =?UTF-8?q?=20can=E2=80=99t=20go=20out=20of=20sync=20when=20adding=20anoth?= =?UTF-8?q?er=20channel.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- doc/documentation.tex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/documentation.tex b/doc/documentation.tex index 317356d..c757788 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -168,7 +168,7 @@ \subsection{Detector Configurations}\label{sec:detector-configurations} \subsection{Interaction Channels} \label{sec:interaction-channels} sntools supports multiple different interaction channels described in this section. -By default, it will generate events across all channels that are available in the selected detector, but it can be restricted to a single channel by using the \texttt{--channel } command line argument, where \texttt{} can be one of \texttt{ibd}, \texttt{es}, \texttt{ps}, \texttt{o16e}, \texttt{o16eb}, \texttt{c12e}, \texttt{c12eb} and \texttt{c12nc}. +By default, it will generate events across all channels that are available in the selected detector, but it can be restricted to a single channel by using the \texttt{--channel } command line argument, with possible values (\texttt{ibd}, \texttt{es}, \texttt{ps}, \texttt{o16e}, etc.) described in the rest of this subsection. For all supported channels, sntools includes code to calculate the differential cross sections, outgoing particle energy as a function of neutrino energy and other quantities. These could be used as a library from other Python code, e.\,g. as follows:\\