lots of cleanup

* minor documentation updates
* add __repr__ to BaseChannel class
* add documentation.pdf to .gitignore file
* follow flake8 conventions for code formatting
* delete undocumented HK-specific energy threshold output (if --verbose)

.gitignore

@@ -11,6 +11,7 @@ dist/
 *.egg-info/
 
 # LaTeX documentation
+documentation.pdf
 *.aux
 *.bbl
 *.blg

doc/documentation.tex

@@ -169,7 +169,8 @@ \subsection{Interaction Channels} \label{sec:interaction-channels}
 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:\\
 \texttt{>> from sntools.interaction\_channels import ibd}\\
-\texttt{>> ibd.bounds\_eE(eNu=20)  \# min/max energy of outgoing positron in MeV}\\
+\texttt{>> c = ibd.Channel('eb')}\\
+\texttt{>> c.bounds\_eE(eNu=20)  \# min/max energy of outgoing positron in MeV}\\
 \texttt{[17.941220926704954, 18.70577898514169]}
 
 Type \texttt{help(ibd)} on the Python command line for documentation or see \href{https://github.com/SNOwGLoBES/snowglobes/pull/12/commits/f1cf96d8c326b99ac35474bb6a36ef4d1fde809d}{this SNOwGLoBES pull request} for a full usage example.
@@ -445,23 +446,23 @@ \subsection{Add a New Detector Material}\label{sec:new-detector-material}
 \subsection{Add a New Interaction Channel}\label{sec:new-interaction-channel}
 
 \begin{itemize}
-\item In the \texttt{fluxes/} folder, create a new file for your interaction channel. It must define two things:
+\item In the \texttt{interaction\_channels/} folder, create a new file for your interaction channel. It must define two things:
 	\begin{itemize}
 		\item A list of neutrino flavors that can interact in this channel (e.g. if all flavors interact in your channel: \texttt{possible\_flavors = ["e", "eb", "x", "xb"]})
-		\item A class named \texttt{Channel}, which inherits from the \texttt{BaseChannel} class defined in \texttt{interaction\_channels/\_\_init\_\_.py}. Overwrite all functions that are marked \texttt{@abstractmethod} in the \texttt{BaseChannel} class to implement the relevant physics.
-		(See docstrings for explanations of all variables and functions.)
+		\item A class named \texttt{Channel}, which inherits from the \texttt{BaseChannel} class defined in \texttt{interaction\_channels/\_\_init\_\_.py}. Overwrite all functions and properties that are marked \texttt{@abstractmethod} in the \texttt{BaseChannel} class to implement the relevant physics.
+		(See docstrings for explanations.)
 	\end{itemize}
 \item In \texttt{genevts.py}, add the channel in the \texttt{parse\_command\_line\_options()} function.
 \item In \texttt{detectors.py}, add the channel to all relevant detector materials.
 \end{itemize}
 
 
 \subsection{Add a New Input Format}\label{sec:new-input-format}
-Ask whoever provided you with that input file, whether they could make the fluxes available in the \texttt{gamma} format described in section~\ref{sec:format-gamma}.
-If that’s not possible:
+First, check if it is straightforward to transform the input files into one of the supported formats described in section~\ref{sec:input-formats}.
+If not:
 \begin{itemize}
-\item In the \texttt{formats/} folder, copy \texttt{gamma.py} to create a new module for your interaction channel and implement the three functions in there.
-(See the doc strings in the source code for details.)
+\item In the \texttt{formats/} folder, create a new file for your interaction channel. It must contain a class named \texttt{Flux}, which inherits from the \texttt{BaseFlux} class defined in \texttt{formats/\_\_init\_\_.py}. Overwrite all functions that are marked \texttt{@abstractmethod} in the \texttt{BaseFlux} class.
+(See docstrings for explanations.)
 \item In \texttt{genevts.py}, add the format in the \texttt{parse\_command\_line\_options()} function.
 \end{itemize}
 

fluxes/sample-gamma.txt

@@ -7,9 +7,10 @@
 #   Column 4: luminosity of nu_e (erg/s)
 #   Column 5-7: like 2-4, but for anti-nu_e
 #   Column 8-10: like 2-4, but for nu_x (i.e. one of nu_mu, nu_tau or their antiparticles)
+# Note: nu_x luminosity in column 10 is for one neutrino type, not the sum of all four.)
 
 # Example usage:
-# $ python genevts.py fluxes/sample-gamma.txt --format gamma --starttime 110 --endtime 120 -v
+# $ python sntools/genevts.py fluxes/sample-gamma.txt --format gamma --starttime 110 --endtime 120 -v
 
 # The nu_e sample data was generated using
 #   time = 0.1 + 0.0005 * i + 0.00001 * (random() - 0.5))

sntools/channel.py

@@ -32,8 +32,6 @@ def gen_evts(_channel, _flux, scale, seed, verbose):
     # To save time, we cache results in a dictionary.
     cached_flux = {}
 
-    thr_e = 3.511  # detection threshold in HK: 3 MeV kinetic energy + rest mass
-
     # integrate over eE and then eNu to obtain the event rate at time t
     raw_nevts = [scale * integrate.nquad(ddEventRate, [channel.bounds_eE, channel.bounds_eNu], args=[t], opts=[channel._opts, {}])[0]
                  for t in flux.raw_times]
@@ -55,14 +53,6 @@ def gen_evts(_channel, _flux, scale, seed, verbose):
     binned_nevt = np.random.poisson(binned_nevt_th)  # Get random number of events in each bin from Poisson distribution
     flux.prepare_evt_gen(binned_t)  # give flux script a chance to pre-compute values
 
-    if verbose:  # compute events above threshold energy `thr_e`
-        thr_bounds_eE = lambda _eNu, *args: [max(thr_e, channel.bounds_eE(_eNu)[0]), max(thr_e, channel.bounds_eE(_eNu)[1])]
-        thr_raw_nevts = [scale * integrate.nquad(ddEventRate, [thr_bounds_eE, channel.bounds_eNu], args=[t], opts=[channel._opts, {}])[0]
-                         for t in flux.raw_times]
-        thr_event_rate = interpolate.pchip(flux.raw_times, thr_raw_nevts)
-        thr_binned_nevt_th = thr_event_rate(binned_t)
-        thr_nevt = sum(binned_nevt)
-
     events = []
     for i in range(n_bins):
         t0 = flux.starttime + i * bin_width
@@ -78,13 +68,7 @@ def gen_evts(_channel, _flux, scale, seed, verbose):
             evt.time = t0 + random.random() * bin_width
             events.append(evt)
 
-            if verbose and evt.outgoing_particles[0][1] < thr_e:
-                thr_nevt -= 1
-
     print(f"[{tag}] Generated {sum(binned_nevt)} particles (expected: {sum(binned_nevt_th):.2f} particles)")
-    if verbose:
-        print(f"[{tag}] -> above threshold of {thr_e} MeV: {thr_nevt} particles (expected: {sum(thr_binned_nevt_th):.2f})")
-        print("**************************************")
 
     return events
 
@@ -129,9 +113,8 @@ def rejection_sample(dist, min_val, max_val, n_bins=100):
 
 def get_eNu(time):
     """Get energy of interacting neutrino using rejection sampling."""
-    dist = lambda _eNu: integrate.quad(
-        ddEventRate, *channel.bounds_eE(_eNu), args=(_eNu, time), points=channel._opts(_eNu)["points"]
-    )[0]
+    def dist(eNu):
+        return integrate.quad(ddEventRate, *channel.bounds_eE(eNu), args=(eNu, time), points=channel._opts(eNu)["points"])[0]
     eNu = rejection_sample(dist, *channel.bounds_eNu, n_bins=200)
     return eNu
 
@@ -140,7 +123,8 @@ def get_direction(eNu):
     """Get direction of outgoing particle using rejection sampling.
     (Assumes that incoming neutrino with energy eNu moves in z direction.)
     """
-    dist = lambda _cosT: channel.dSigma_dCosT(eNu, _cosT)
+    def dist(cosT):
+        return channel.dSigma_dCosT(eNu, cosT)
     cosT = rejection_sample(dist, -1, 1, 200)
     sinT = sin(acos(cosT))
     phi = 2 * pi * random.random()  # randomly distributed in [0, 2 pi)

sntools/event.py

@@ -47,11 +47,11 @@ def ratpac_string(self, i, events):
         GeV = 0.001   # convert from MeV
         mm = 10       # convert from cm
         ns = 1000000  # convert from ms
-        
+
         dt = self.time
         if i > 0:
-            dt -= events[i-1].time
-        
+            dt -= events[i - 1].time
+
         s = f"{len(self.outgoing_particles)}\n"
         for idx, (pid, e, dirx, diry, dirz) in enumerate(self.outgoing_particles):
             mass = 0.0
@@ -63,9 +63,9 @@ def ratpac_string(self, i, events):
                 mass = 939.56563
             p2 = (e**2) - (mass**2)
             p = p2**0.5
-            px = dirx*p
-            py = diry*p
-            pz = dirz*p
+            px = dirx * p
+            py = diry * p
+            pz = dirz * p
             if idx > 0:
                 dt = 0.0
             s += f"1 {pid} 0 0 {px * GeV:.8e} {py * GeV:.8e} {pz * GeV:.8e} {mass * GeV:.8e} {dt * ns:.5e} {self.vertex[0] * mm:.5e} {self.vertex[1] * mm:.5e} {self.vertex[2] * mm:.5e}\n"

sntools/formats/__init__.py

@@ -54,7 +54,7 @@ def nu_emission(self, eNu, time):
 class WeightedFlux(BaseFlux):
     def __init__(self, unweighted_flux, scale, distance=10.0) -> None:
         """Initialize a WeightedFlux.
-        
+
         Arguments:
         unweighted_flux -- BaseFlux instance corresponding to original flux (at supernova)
         scale -- scaling factor from flux transformation, e.g. 1.0 for NoTransformation

sntools/formats/gamma.py

@@ -46,11 +46,10 @@ def parse_input(self, input, inflv, starttime, endtime):
 
         # save mean energy, mean squared energy, luminosity to dictionary to look up in nu_emission() below
         self.flux = {}
+        # input files contain time in column 1 and data for nu_e (cols 2-4), anti-nu_e (cols 5-7) and nu_x (cols 8-10)
+        offset = {"e": 1, "eb": 4, "x": 7, "xb": 7}[inflv]
         for timebin in indata:
-            # input files contain information for nu_e in columns 1-3, for
-            # anti-nu_e in cols 4-6 and for nu_x in columns 7-9
-            offset = {"e": 1, "eb": 4, "x": 7, "xb": 7}[inflv]
-            (mean_e, mean_e_sq, lum) = timebin[offset : offset + 3]
+            (mean_e, mean_e_sq, lum) = timebin[offset: offset + 3]
             t = timebin[0]
             self.flux[t] = (mean_e, mean_e_sq, lum * 624.151)  # convert lum from erg/s to MeV/ms
 

sntools/formats/nakazato.py

@@ -26,7 +26,7 @@ def parse_input(self, input, inflv, starttime, endtime):
             indata = [list(map(float, line.split())) for line in infile if not (line.startswith("#") or line.isspace())]
 
         # 21 lines per time bin
-        chunks = [indata[21 * i : 21 * (i + 1)] for i in range(len(indata) // 21)]
+        chunks = [indata[21 * i: 21 * (i + 1)] for i in range(len(indata) // 21)]
 
         # input files contain information for e, eb & x in neighbouring columns,
         # so depending on the flavor, we might need an offset

sntools/formats/princeton.py

@@ -43,7 +43,7 @@ def parse_input(self, input, inflv, starttime, endtime):
             self.times.append(time)
 
             diff_number_flux = [0]  # Set flux at 0 MeV to 0
-            for emean, diff_lum in zip(ebins[1:-1], line[offset : offset + 20]):
+            for emean, diff_lum in zip(ebins[1:-1], line[offset: offset + 20]):
                 diff_lum *= 1e50  # file gives spectral luminosity in 10^50 erg/s/MeV
                 diff_lum *= 624.151  # convert erg/s/MeV to MeV/ms/MeV
                 if offset == 41:

sntools/formats/totani.py

@@ -119,13 +119,13 @@ def _parse(self, input, format, flv):
         if format == "early":
             # 42 lines per time bin, 26 bins in wilson-early.txt
             for i in range(26):
-                chunks.append(raw_indata[42 * i : 42 * (i + 1)])
+                chunks.append(raw_indata[42 * i: 42 * (i + 1)])
             line_N = 6
             range_egroup = range(19, 39)
         elif format == "late":
             # 46 lines per time bin, 36 bins in wilson-late.txt
             for i in range(36):
-                chunks.append(raw_indata[46 * i : 46 * (i + 1)])
+                chunks.append(raw_indata[46 * i: 46 * (i + 1)])
             line_N = 8
             range_egroup = range(21, 41)
 
@@ -171,7 +171,7 @@ def _parse_nb(self, input):
 
         # 26 lines per time bin, 99 bins in wilson-nb.txt. Bin 6 is equivalent to
         # 40ms post-bounce & bin 56 is 50ms, so we only select that range:
-        chunks = [raw_indata[26 * i : 26 * (i + 1)] for i in range(6, 57)]
+        chunks = [raw_indata[26 * i: 26 * (i + 1)] for i in range(6, 57)]
 
         # for each time bin, save data to dictionaries to look up later
         for chunk in chunks:

sntools/interaction_channels/__init__.py

@@ -11,6 +11,9 @@ class BaseChannel(ABC):
     def __init__(self, flv) -> None:
         self.flavor = flv
 
+    def __repr__(self) -> str:
+        return f"{self.__module__}.Channel('{self.flavor}')"
+
     @abstractmethod
     def generate_event(self, eNu, dirx, diry, dirz):
         """Return an event with the appropriate incoming/outgoing particles.

sntools/interaction_channels/c12e.py

@@ -69,9 +69,9 @@ def dSigma_dE(self, eNu, eE):
             # ... but just in case:
             return 0
 
-        sigma0 = 3.439E-44 * (5.067731E10)**2  # convert cm^2 to MeV^-2, see http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1)
         a1, a2, a3 = 10.164, -0.4666, 0.0546
-        sigma = sigma0 * (a1 * (eNu - e_thr) + a2 * (eNu - e_thr)**2 + a3 * (eNu - e_thr)**3)
+        sigma = 3.439e-44 * (a1 * (eNu - e_thr) + a2 * (eNu - e_thr)**2 + a3 * (eNu - e_thr)**3)
+        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):

sntools/interaction_channels/c12eb.py

@@ -69,9 +69,9 @@ def dSigma_dE(self, eNu, eE):
             # ... but just in case:
             return 0
 
-        sigma0 = 9.55E-46 * (5.067731E10)**2  # convert cm^2 to MeV^-2, see http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1)
         a1, a2, a3 = -106.3, 25.15, 0.3697
-        sigma = sigma0 * (a1 * (eNu - e_thr) + a2 * (eNu - e_thr)**2 + a3 * (eNu - e_thr)**3)
+        sigma = 9.55e-46 * (a1 * (eNu - e_thr) + a2 * (eNu - e_thr)**2 + a3 * (eNu - e_thr)**3)
+        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):

sntools/interaction_channels/c12nc.py

@@ -73,9 +73,8 @@ def dSigma_dE(self, eNu, eE):
             # ... but just in case:
             return 0
 
-        sigma0 = 1.08E-38 * (5.067731E10)**2  # convert cm^2 to MeV^-2, see http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1)
-        sigma = sigma0 * ((eNu - e_thr) / 939)**2
-        sigma *= 1.11**2  # coupling constant \beta (see discussion above)
+        sigma = 1.08E-38 * ((eNu - e_thr) / 939)**2 * 1.11**2  # Assume beta * kappa = 1.11 (see discussion above)
+        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):

sntools/interaction_channels/o16e.py

@@ -102,7 +102,7 @@ def dSigma_dE(self, eNu, eE):
         for g in range(1, 5):
             sigma += self.partial_dSigma_dE(eNu, eE, g)
 
-        sigma *= (5.067731E10)**2  # convert cm^2 to MeV^-2, see http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1)
+        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 partial_dSigma_dE(self, eNu, eE, g):  # eq. (4) of arXiv:1809.08398

sntools/interaction_channels/o16eb.py

@@ -102,7 +102,7 @@ def dSigma_dE(self, eNu, eE):
         for g in range(1, 5):
             sigma += self.partial_dSigma_dE(eNu, eE, g)
 
-        sigma *= (5.067731E10)**2  # convert cm^2 to MeV^-2, see http://www.wolframalpha.com/input/?i=cm%2F(hbar+*+c)+in+MeV%5E(-1)
+        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 partial_dSigma_dE(self, eNu, eE, g):  # eq. (4) of arXiv:1809.08398

sntools/transformation.py

@@ -87,7 +87,7 @@ def __init__(self, name):
 
         if xf.prob_ee(0, 10) != xf.prob_ee(0, 20) or xf.prob_ee(0, 10) != xf.prob_ee(0.2, 10):
             # Probability appears to depend on time and/or energy
-            raise ValueError(f"The transformation '{name}' from SNEWPy appears to be time- or energy-dependent. This is not currently supported.")
+            raise ValueError(f"The transformation 'SNEWPY-{name}' is time- or energy-dependent. This is not yet supported.")
 
         transformation = (
             ("e", float(xf.prob_ee(0, 0)), "e"),  # nu_e that originated as nu_e