mPMT support

[JostMigenda]

Initial draft, still missing calibration. See discussion in #12.

[gpronost]

Thank you, I am testing it on sukap inside the AstroAnalysis framework

[gpronost]

It seems to work on a technical point of view, but there is a minor issue:
-> The number of PMT and mPMT in the default version are not the values from WCSim hybrid.
Fortunately the true number is smaller, so we can solve it by resetting the number of PMT and mPMT, but we should recompute the nearest neighbour after this. In order to do so, it would be nice to have a function allowing the overwrite directly, something like this:
inline void OverwriteNearestNeighbour() { GetNearestNeighbours(true); }

[JostMigenda]

Changing the number of (m)PMTs and then recomputing the neighbours might lead to problems right now. That’s because the name of the file that those neighbours are cached to only contains the detector geometry and neighbour distance. You’d need to remember to explicitly recalculate that file not just when you change the number of PMTs, but also at the beginning of every script you run; just in case the previous script you ran e.g. a week earlier used a different number of PMTs. So that seems very unsafe to me and I’d rather avoid doing that.

If you have a small number of different configurations (and always mask exactly the same PMT IDs), then it would probably be enough to add the number of PMTs to the file name. If that doesn’t work, could you explain your use case briefly? Hopefully we’ll be able to come up with another way.

[gpronost]

The issue is much simpler than what you are suggesting:

In WCSim Hybrid the number of Box and Line PMT and of 3inch PMT is not what you set as default:
[3]: INFO: BnL PMT: 19208
[3]: INFO: 3inch PMT: 182704
The difference is likely to come from the tank size difference between the Hybrid WCSim (which use the official HK tank dimensions) and the Official WCSim (which use outdated values).

Currently, FLOWER calculate the neighbour PMTs in the constructor function, which means that there is no opportunity to correct the number of Tubes beforehand. So we either we need to fix the values in FLOWER, either we need to overwrite the neighbour PMTs map (and this needs to be done only once, since the number of PMTs will not change frequently).

Edit: For masking PMT, this is a different issue, and I agree, recomputing neighbour PMTs map each time won't be good. However, I think that for the time being, we will use masking only on mPMTs, not on BnL PMTs, so the neighbour PMTs map should be affected. So, let's ignore this issue for now.

[JostMigenda]

Ah, thanks for clarifying! I’ll change the PMT numbers for the hybrid geometry in a moment. For masking, I agree—let’s ignore that for now and we’ll revisit it if/when necessary.

[gpronost]

For the MC Production, there are several configuration of PMT, I guess we need to implement them. I will make a Pull Request with the new numbers (for WCSim Hybrid).

[gpronost]

Here it is:
#16

[JostMigenda]

Note that in the recent commit, I use the same N_eff ⇔ Energy relation for the hybrid geometries with 3k, 5k and 10k mPMTs. In principle, the algorithm is independent of the number of (m)PMTs; though for major changes (e.g. 20% vs. 40% or B&L only vs. hybrid) it is better to recalibrate.

[JostMigenda]

One more thing I just remembered: @gpronost Is there a significant difference in photon detection efficiency (i.e. quantum efficiency * capture efficiency) between the B&L PMT and the 3" PMTs? In the original FLOWER (with just one PMT type), this was simply a constant factor that was absorbed into the N_eff ⇔ energy conversion; but if it is different between PMT types, that wouldn’t work as well (or would require recalibration for every geometry that changes the ratio of the PMT types.)

[JostMigenda]

WCSim gives the quantum efficiency as a function of wavelength of different PMT types in WCSimPMTObject.cc. Of course, we don’t know the wavelength of each photon in FLOWER, so we have to average over the Cherenkov spectrum.
The number of Cherenkov photons at a given wavelength is proportional to 1/λ² (derivation e.g. here), so integrating (with the Python code snippet below) gives a spectrum-averaged QE of:

B&L PMTs: 0.2194
mPMTs: 0.2268

So overall, there’s a 3.4% difference between both types. This is effectively absorbed in the Neff ⇔ Energy calibration, but it means that it wouldn’t be quite accurate to use the same calibration for detector configurations with different ratios of B&L PMTs to mPMTs. I’ll see if I can update the code in the next few days to separate the QE out.

from scipy import integrate, interpolate
wavelengths = range(280,680,20)
correctionFactor = 1/0.73
qe_bl = interpolate.pchip(wavelengths, [0.00*correctionFactor, .0008*correctionFactor, .1255*correctionFactor, .254962*correctionFactor, .2930*correctionFactor, .3127*correctionFactor, .3130*correctionFactor, .2994*correctionFactor, .2791*correctionFactor, .2491*correctionFactor, .2070*correctionFactor,  .1758*correctionFactor, .1384*correctionFactor, .0779*correctionFactor, .0473*correctionFactor, .0288*correctionFactor, .0149*correctionFactor, .0062*correctionFactor, .0002*correctionFactor, .0001*correctionFactor])
qe_mpmt = interpolate.pchip(wavelengths, [0.00*correctionFactor, .0787*correctionFactor, .1838*correctionFactor, .2401*correctionFactor, .2521*correctionFactor, .2695*correctionFactor, .2676*correctionFactor, .2593*correctionFactor, .2472*correctionFactor, .2276*correctionFactor, .1970*correctionFactor,  .1777*correctionFactor, .1547*correctionFactor, .1033*correctionFactor, .0727*correctionFactor, .0587*correctionFactor, .0470*correctionFactor, .0372*correctionFactor, .0285*correctionFactor, .0220*correctionFactor])
normalization = integrate.quad(lambda l: 1/l**2, 280, 680)[0]
print("B&L PMTs:", integrate.quad(lambda l: qe_bl(l)/l**2 / normalization, 280, 680)[0])
print("mPMTs:", integrate.quad(lambda l: qe_mpmt(l)/l**2 / normalization, 280, 680)[0])