update

common/data/Materials/Surf_Hamamatsu_R15458.dat

@@ -4,7 +4,54 @@
    "jType": "coated",
     "jCoatThicknessNM": 20,
    "jProperties": [
-      "RINDEX"
+      "RINDEX",
+      "IMAGINARYRINDEX"
+   ],
+   "jWavelength_IMAGINARYRINDEX": [
+      380,
+      395,
+      410,
+      425,
+      440,
+      455,
+      470,
+      485,
+      500,
+      515,
+      530,
+      545,
+      560,
+      575,
+      590,
+      605,
+      620,
+      635,
+      650,
+      665,
+      680
+   ],
+   "jValues_IMAGINARYRINDEX": [
+      1.69,
+      1.69,
+      1.71,
+      1.53,
+      1.50,
+      1.44,
+      1.34,
+      1.11,
+      1.06,
+      1.05,
+      0.86,
+      0.63,
+      0.53,
+      0.46,
+      0.42,
+      0.38,
+      0.37,
+      0.35,
+      0.34,
+      0.34,
+      0.33
    ],
    "jValues_RINDEX": [
       2.05939,

common/framework/include/OMSimOpBoundaryProcess.hh

@@ -311,6 +311,7 @@ private:
   // newly added
   size_t idx_abslength1 = 0;
   size_t idx_abslength2 = 0;
+  size_t idx_coatedimagindex = 0;
   size_t idx_coatedabslength = 0;
 
   G4bool fInvokeSD;

common/framework/include/OMSimTools.hh

@@ -31,7 +31,8 @@ namespace Tools
     std::vector<std::vector<double>> loadtxt(const std::string &pFilePath,
                                              bool pUnpack = true,
                                              size_t pSkipRows = 0,
-                                             char pDelimiter = ' ');
+                                             char pDelimiter = ' ',
+                                             char pComments = '#');
     std::vector<double> linspace(double start, double end, int num_points, bool endpoint = true);
     std::vector<double> logspace(double start, double end, int num_points, double base = 10.0, bool endpoint = true);
     void sortVectorByReference(std::vector<G4double> &referenceVector, std::vector<G4double> &sortVector);

common/framework/src/OMSim.cc

@@ -147,6 +147,7 @@ void OMSim::initialiseSimulation(OMSimDetectorConstruction* pDetectorConstructio
     mNavigator = std::make_unique<G4Navigator>();
 
     int nThreads = determineNumberOfThreads();
+
     mRunManager->SetNumberOfThreads(nThreads);
 
     mRunManager->SetUserInitialization(pDetectorConstruction);

common/framework/src/OMSimOpBoundaryProcess.cc

@@ -1602,7 +1602,7 @@ void G4OpBoundaryProcess::PhotocathodeCoatedComplex(const G4Track *pTrack)
 
   G4MaterialPropertiesTable *lMPT1 = fMaterial1->GetMaterialPropertiesTable();
   G4MaterialPropertyVector *lRIndexMPV1 = lMPT1 ? lMPT1->GetProperty(kRINDEX) : nullptr;
-  G4MaterialPropertyVector *lAbsLengthMPV1 = lMPT1? lMPT1->GetProperty(kABSLENGTH) : nullptr;
+  G4MaterialPropertyVector *lAbsLengthMPV1 = lMPT1 ? lMPT1->GetProperty(kABSLENGTH) : nullptr;
 
   fAbsorptionLength1 = lAbsLengthMPV1 ? lAbsLengthMPV1->Value(fPhotonMomentum, idx_abslength1) : 0;
 
@@ -1613,24 +1613,30 @@ void G4OpBoundaryProcess::PhotocathodeCoatedComplex(const G4Track *pTrack)
   fRindex2 = lRIndexMPV2 ? lRIndexMPV2->Value(fPhotonMomentum, idx_rindex2) : 0;
   fAbsorptionLength2 = lAbsLengthMPV2 ? lAbsLengthMPV2->Value(fPhotonMomentum, idx_abslength2) : 0;
 
-
   G4MaterialPropertiesTable *lMPTCoated = fOpticalSurface->GetMaterialPropertiesTable();
 
+  G4double lCoatedImagRIndex = 0;
   if (lMPTCoated)
   {
 
     G4MaterialPropertyVector *lPpCoated;
     if ((lPpCoated = lMPTCoated->GetProperty(kRINDEX)))
       lCoatedRindex = lPpCoated->Value(fPhotonMomentum, idx_coatedrindex);
-    if ((lPpCoated = lMPTCoated->GetProperty(kABSLENGTH)))
+    if ((lPpCoated = lMPTCoated->GetProperty(kIMAGINARYRINDEX)))
+    {
+      lCoatedImagRIndex = lPpCoated->Value(fPhotonMomentum, idx_coatedimagindex);
+    }
+    else if ((lPpCoated = lMPTCoated->GetProperty(kABSLENGTH)))
+    {
       lCoatedAbsLength = lPpCoated->Value(fPhotonMomentum, idx_coatedabslength);
+      lCoatedImagRIndex = lCoatedAbsLength != 0 ? lWavelength / (4 * pi * lCoatedAbsLength) : 0;
+    }
     lCoatedThickness = lMPTCoated->ConstPropertyExists(kCOATEDTHICKNESS) ? lMPTCoated->GetConstProperty(kCOATEDTHICKNESS) : 0;
   }
 
   // Convert absorption length to complex refractive index
   G4double limagRIndex1 = fAbsorptionLength1 != 0 ? lWavelength / (4 * pi * fAbsorptionLength1) : 0;
   G4double limagRIndex2 = fAbsorptionLength2 != 0 ? lWavelength / (4 * pi * fAbsorptionLength2) : 0;
-  G4double lCoatedImagRIndex = lCoatedAbsLength != 0 ? lWavelength / (4 * pi * lCoatedAbsLength) : 0;
 
   // Get material before the boundary
   G4VUserTrackInformation *lUserInformation = pTrack->GetUserInformation();
@@ -1730,7 +1736,6 @@ void G4OpBoundaryProcess::PhotocathodeCoatedComplex(const G4Track *pTrack)
       lCost1Complex = std::sqrt(G4complex(1.0, 0.0) - lSint1 * lSint1);
     }
 
-
     FresnelCoefficients lCoefficients1TL = CalculateFresnelCoefficientsComplex(lComplexRindex1, lComplexCoatedRindex, lCost1Complex, lCostTLComplex);
     FresnelCoefficients lCoefficientsTL2 = CalculateFresnelCoefficientsComplex(lComplexCoatedRindex, lComplexRindex2, lCostTLComplex, lCost2Complex);
     G4complex lBeta = lK0 * lComplexCoatedRindex * lCoatedThickness * lCostTLComplex;
@@ -1869,12 +1874,10 @@ OpticalLayerResult G4OpBoundaryProcess::Fresnel2ProbabilityComplex(FresnelCoeffi
   G4double lReflTE = std::norm(pCoefficients.rTE);
   G4double lReflTM = std::norm(pCoefficients.rTM);
 
-
   G4double prefactor = std::real((lComplexRindex2 * pCost2) / (lComplexRindex1 * pCost1));
   G4double lTransTE = prefactor * std::norm(pCoefficients.tTE);
   G4double lTransTM = prefactor * std::norm(pCoefficients.tTM);
 
-
   // real probabilities:
   lResult.Reflectivity = (lReflTE + lReflTM) * 0.5;
   lResult.Transmittance = (lTransTE + lTransTM) * 0.5;

common/framework/src/OMSimTools.cc

@@ -165,13 +165,15 @@ namespace Tools
 	 * the file. Default is 0.
 	 * @param pDelimiter The character used to separate values in each line of the
 	 * input file.
+	 * @param pComments Optional. The character used to indicate the start of a
+	 * comment. Default is '#'.
 	 * @return A 2D vector of doubles. The outer vector groups all columns (or
 	 * rows if 'unpack' is false), and each inner vector represents one of the
 	 * columns (or one of the rows if 'unpack' is false) of data file.
 	 * @throws std::runtime_error if the file cannot be opened.
 	 */
 	std::vector<std::vector<double>> loadtxt(const std::string &pFilePath, bool pUnpack,
-											 size_t pSkipRows, char pDelimiter)
+											 size_t pSkipRows, char pDelimiter, char pComments)
 	{
 		std::vector<std::vector<double>> lData;
 		std::ifstream lInFile(pFilePath);
@@ -189,6 +191,8 @@ namespace Tools
 		{
 			if (lRowCounter++ < pSkipRows)
 				continue;
+			if (!lLine.empty() && lLine[0] == pComments)
+            	continue;
 			std::vector<double> lRow;
 			std::stringstream lSs(lLine);
 			std::string lItem;

documentation/extra_doc/00_PMT.md

@@ -4,6 +4,7 @@
 
 ### QE calibration
 
+#### Step 0
 ```py
 import numpy as np
 from scipy.optimize import curve_fit, minimize
@@ -81,19 +82,14 @@ class CalibrationCurve:
         except RuntimeError:
             print("Error: Failed to fit single exponential function.")
 
-    def plot_single_exp(self):
+    def plot_single_exp(self, ax):
         """Plot the single exponential fit."""
         if self.popt_single is not None:
-            x = np.linspace(min(self.abs_length), max(self.abs_length), 100)
-            plt.figure(figsize=(10, 6))
-            plt.plot(x, single_exp(x * 1e9, *self.popt_single), 'b-', label='Single Exp Fit')
-            plt.plot(self.abs_length, self.fraction, 'ro', label='Data')
-            plt.xlabel('Absorption Length (m)')
-            plt.ylabel('Fraction')
-            plt.legend()
-            plt.title('Single Exponential Fit')
-            plt.show()
-
+            x = np.logspace(np.log10(np.amin(self.abs_length[self.mask])), -4, 100)
+            ax.plot(x, single_exp(x * 1e9, *self.popt_single), '--', label='Single Exp Fit')
+            ax.errorbar(self.abs_length[self.mask], self.fraction[self.mask], yerr = self.error[self.mask],
+                    fmt =  '.', label='Fitted Data')
+
     def get_needed_abs(self, QE: float) -> float:
         """Calculate the needed absorption length for a given quantum efficiency."""
         if self.popt_single is None:
@@ -115,8 +111,14 @@ plt.loglog()
 
 amplitude = [data[key].popt_single[0] for key in data]
 eff_thick = [data[key].popt_single[1] for key in data]
+max_QE = [curve.max_pos_y for curve in data.values()]
 wavelengths = [key for key in data]
 np.savetxt("mDOM_Hamamatsu_R15458_QE_matching_parameters.dat",
-           np.array([wavelengths_all, amp, eff_thick]).T,
-           delimiter="\t", header="Wavelength(nm) \t Amplitude \t Effective thickness (nm)")
-```
+           np.array([wavelengths_all, amp, eff_thick, max_QE]).T,
+           delimiter="\t", header="Wavelength(nm) \t Amplitude \t Effective thickness (nm) \t Max. QE")
+```
+
+
+#### Step 1
+
+#### Step 2