scan data matching

common/framework/include/OMSimTools.hh

@@ -43,6 +43,10 @@ namespace Tools
     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);
+    double median(std::vector<double> p_vec);
+    double mean(const std::vector<double> &p_vec, const std::vector<double> &p_weights = {});
+    double std(const std::vector<double>& vec,  const std::vector<double> &p_weights = {});
+
     void throwError(const G4String& message);
     std::vector<G4String> splitStringByDelimiter(G4String const &p_string, char p_delim);
     std::vector<G4String> splitStringByDelimiter(char *p_char, char p_delim);

common/framework/src/OMSimTools.cc

@@ -482,5 +482,96 @@ namespace Tools
 		return splitStringByDelimiter(G4String(p_char), p_delim);
 	}
 
+	double median(std::vector<double> p_vec)
+	{
+		if (p_vec.empty())
+		{
+			throwError("Vector is empty");
+		}
+		auto n = p_vec.size() / 2;
+		std::nth_element(p_vec.begin(), p_vec.begin() + n, p_vec.end());
+		if (p_vec.size() % 2 == 0)
+		{
+			return (p_vec[n - 1] + p_vec[n]) / 2.0;
+		}
+		else
+		{
+			return p_vec[n];
+		}
+	}
+	double mean(const std::vector<double> &p_vec, const std::vector<double> &p_weights)
+	{
+		if (p_vec.empty())
+		{
+			throwError("Vector is empty");
+		}
+		if (!p_weights.empty() && p_vec.size() != p_weights.size())
+		{
+			throwError("Vectors and weights must have the same size");
+		}
+
+		if (p_weights.empty())
+		{
+			// Calculate the mean of p_vec without weights
+			return std::accumulate(p_vec.begin(), p_vec.end(), 0.0) / p_vec.size();
+		}
+		else
+		{
+			// Calculate the weighted mean
+			double weighted_sum = 0.0;
+			double weight_sum = 0.0;
+			for (size_t i = 0; i < p_vec.size(); ++i)
+			{
+				weighted_sum += p_vec[i] * p_weights[i];
+				weight_sum += p_weights[i];
+			}
+			return weighted_sum / weight_sum;
+		}
+	}
+
+	double std(const std::vector<double> &p_vec, const std::vector<double> &p_weights)
+	{
+		if (p_vec.size() < 2)
+		{
+			throwError("Vector must have at least two elements");
+		}
+
+		double meanV = mean(p_vec, p_weights);
+
+		if (p_weights.empty())
+		{
+			// Unweighted case
+			double sum_sq_diff = 0.0;
+			for (const auto &value : p_vec)
+			{
+				double diff = value - meanV;
+				sum_sq_diff += diff * diff;
+			}
+			return std::sqrt(sum_sq_diff / (p_vec.size() - 1));
+		}
+		else
+		{
+			// Weighted case
+			if (p_vec.size() != p_weights.size())
+			{
+				throwError("Vectors and weights must have the same size");
+			}
+
+			double sum_sq_diff = 0.0;
+			double sum_weights = 0.0;
+			double sum_weights_squared = 0.0;
+
+			for (size_t i = 0; i < p_vec.size(); ++i)
+			{
+				double diff = p_vec[i] - meanV;
+				sum_sq_diff += p_weights[i] * diff * diff;
+				sum_weights += p_weights[i];
+				sum_weights_squared += p_weights[i] * p_weights[i];
+			}
+
+			return std::sqrt(sum_sq_diff * sum_weights / (sum_weights * sum_weights - sum_weights_squared));
+		}
+	}
+
 	std::string visualisationURL = "https://icecube.github.io/OMSim/md_extra_doc_2_technicalities.html#autotoc_md30";
 }

documentation/extra_doc/2_technicalities.md

@@ -137,7 +137,7 @@ To calculate the quantum efficiency weight, we have to know how many photons are
  - In `OMSim_efficiency_calibration.cc` change the method `runQEbeamSimulation()` to use the newly defined beam instead of `runErlangenQEBeam`. Define there also the wavelength range to be simulated.
  - Run the simulation with the following command, changing the PMT number to the one you want to simulate (make sure that your PMT is in the list in `OMSimPMTConstruction::selectPMT`, or hard code it in `OMSimEffCaliDetector.cc`)
 ```bash
-./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 1 -n 500000 --threads 4 --output_file step1
+./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 1 -n 500000 --threads 4 --output_file step1 --detector_type 1
 ```
  - Check results and uncertainty. You want to determine the fraction pretty accurately, you may repeat the simulation at the UV region with higher statistics (change wavelength range in `runQEbeamSimulation()`).
 
@@ -178,7 +178,7 @@ Each PMT has its own derived class in `OMSimPMTResponse.cc`. If your PMT has not
 
 Now run the simulation again 
 ```bash
-./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 2 -n 100000 --threads 4 --output_file step2
+./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 2 -n 100000 --threads 4 --output_file step2 --detector_type 1
 ```
 
 and check that the weights are being calculated correctly
@@ -210,11 +210,11 @@ plt.legend()
 The last step is to create the collection efficiency weights to match the relative detection efficiency scans. For this the scan measurement is replicated in the simulation, scanning the PMT in a XY grid. The output file of the simulation of this step is a histogram with the position of absorbed photons for each beam position.
 
  - As before, we have to simulate the beam used during the scan measurement (see for example `Beam::configureXYZScan_PicoQuantSetup` and `Beam::runBeamPicoQuantSetup` for the beam used in Münster)
- - In in `runXYZfrontalScan()` of `OMSim_efficiency_calibration.cc` change the scan range (`lX` vector) and radius limit (`lRlim`) according to the diameter of your PMT
+ - In in `runXYZfrontalScan()` of `OMSim_efficiency_calibration.cc` change the scan range (`grid` vector) and radius limit (`rLim`) according to the diameter of your PMT
  - Also adjust the binning of the output histogram in `OMSimEffiCaliAnalyisis::writeHitPositionHistogram`
  - Run the XY grid simulation. 10000 photons per grid position should be enough, but you may increase / decrease statistics as you want
  ```bash
-./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 3 -n 10000 --threads 4 --output_file step3
+./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 3 -n 10000 --threads 4 --output_file step3 --detector_type 1
 ```
  - Fit the weights using simulation data. The analysis done for the mDOM can be found in the notebook located in `documentation/notebooks/detection_efficiency_matching/`.
 
@@ -224,10 +224,22 @@ The last step is to create the collection efficiency weights to match the relati
 
  ### Step 4: Matching gain / transit time scans
 
-The scan data must be corrected before use, as the coordinates of the beam do not necessarily correspond to the primary spot on the photocathode that is illuminated, since air-glass boundary refracts the beam.
+The scan data of transit time / gain must be corrected before use, as the coordinates of the beam do not necessarily correspond to the primary spot on the photocathode that is illuminated, since air-glass boundary refracts the beam. 
 
- - Set the photocathode as a 100% efficient by 
+ - As before, we have to simulate the beam used during the scan measurement (see for example `Beam::configureXYZScan_PicoQuantSetup` and `Beam::runBeamPicoQuantSetup` for the beam used in Münster)
+ - In in `runXYZfrontalScan()` of `OMSim_efficiency_calibration.cc` change the scan range (`grid` vector) and radius limit (`rLim`) according to the diameter of your PMT
+ - Adjust the binning of the output histogram in `OMSimEffiCaliAnalyisis::writeHitPositionHistogram` to a finner binning as before
+ - Run the XY grid simulation **using the simple PMT**. 10000 photons per grid position should be enough, but you may increase / decrease statistics as you want
 
+ ```bash
+./OMSim_efficiency_calibration --simple_PMT --pmt_model 0 --simulation_step 4 -n 10000 --threads 4 --output_file step4 --detector_type 1
+```
+ - Follow the analysis in the notebook `documentation/notebooks/scans_matching/` and save the created files in `common/data/PMTs/measurement_matching_data/scans/`.  Note that OMSim expects a naming convention for these files (see `OMSimPMTResponse::configureScansInterpolator`).
+ - Modify `getScannedWavelength` of your PMT response class for the wavelengths you scanned and add the `configureScansInterpolator` in the constructor.  
+ - Run the simulation in step 5 to check the output with the newly introduced files
+  ```bash
+./OMSim_efficiency_calibration --pmt_model 0 --simulation_step 5 -n 10000 --threads 4 --output_file step5 --detector_type 1
+```
 ---
 ---