Skip to content

Add updated HK detector geometry

The detector geometry currently implemented in FLOWER (PMT number, dark rate, hits-to-energy tuning, …) is from 2019. There have been significant updates since (ID dimensions have changed slightly, PMT number and arrangement has been fixed, dark rate of delivered PMTs has been measured, number and design of mPMTs has been fixed, …), so once GHOST is ready to simulate the latest geometry, we need to update FLOWER.

Since I’m about to move to a new position soon, I’ll document the necessary steps in some detail in the following.

Adding a new Detector Geometry

  • Add new entry to the kDetector_t enum in WCSimFLOWER.h and return that entry from DetectorEnumFromString in WCSimFLOWER.cpp
  • Define appropriate constants for this geometry in the WCSimFLOWER initialisation.
  • Check whether the corrections in GetNEff() require any modifications for the new detector geometry.
  • Add a fitting function for the new detector geometry to CorrectEnergy()—see below!
  • Optional: Delete outdated geometries that are no longer used.

Creating a Fitting Function That Fits a New Detector Geometry

At a very basic level, FLOWER works as follows:

  1. Get a trigger with a reconstructed event time and vertex from an external tool (e.g. from BONSAI)
  2. Iterate over detected hits in a short trigger time window and apply corrections for detector effects (scattering/absorption, dark noise, photocoverage, …) to estimate the total number of Cherenkov photons emitted by the particle, fNEff. A detailed description of this method (which is based on the one used in SK) and the individual corrections can be found in chapter 3.4.3 of my PhD thesis.
  3. Use a fitting function to match fNEff to a reconstructed particle energy fERec.

Determining a fitting fitting function in step 3 is the tricky bit.

  • Use WCSim/GHOST to generate sufficiently high statistics of monoenergetic events for sufficiently many electron energies
  • Reconstruct these datasets with a dummy fitting function (fERec = fNEff) and determine the average fNEff for each energy
  • Fit a polynomial to these data points, e.g. with the following sample code snippet:
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial

E = list(range(5, 10)) + list(range(10, 20, 2)) + list(range(20, 51, 5))  # list of energies in MeV
Neff = […]  # fill in average `fNEff` for these different energies

p = Polynomial.fit(Neff, E, 1, full=True)
print(f"Least squares linear fit: {p[0].convert()}")
print(f"Sum of squared residuals: {p[1][0]}")

fig, axs = plt.subplots(2)
axs[0].plot(E, Neff, 'x', label='MC')
axs[0].plot([p[0](n) for n in Neff], Neff, label=f"y = {p[0].convert()}")
axs[0].set_ylabel('Neff')
axs[0].legend()

axs[1].plot(E, [(p[0](n)-e) / e for e, n in zip(E, Neff)], 'x')
axs[1].set_ylabel('Relative error')

plt.xlabel('Energy [MeV]')
plt.show()

Note: This relation is expected to be approximately linear, but in practice, we observe slight deviations from linearity. (E.g. SK uses a fourth-order polynomial below 25 MeV and a linear fit above.) Depending on the required precision, consider using higher-order polynomials if needed and if you have sufficiently high statistics to be sure that you’re not overfitting on random noise. The fitting functions in the current version of FLOWER are based on fairly low statistics, because they were for an early version of the geometries where that level of fine-tuning was sufficient. Please do not rely on the current choices when deciding where to set the transition point between different polynomials or what order of polynomial to use!

mPMT-specific Changes

Currently, the main branch of FLOWER only includes detector geometries with a single photosensor type. Since the number of mPMTs is about 4% of the number of 50 cm PMTs (and the photosensitive area of an mPMT is smaller than that of a 50 cm PMT), ignoring mPMTs still gives a fairly good approximation of the low-energy reconstruction performance, where photocoverage is the limiting factor.

Initial support for hybrid geometries, which use a combination of both photosensor types, is pending in PR #14. Please see discussions in #12 to understand the limitations of this initial support.