add calibration factor to DetectionEfficiencies (#165)

* add calibration factorand text field to DetectionEfficiencies

* add calibration_factor to helpers

---------

Co-authored-by: Casper da Costa-Luis <[email protected]>

Author: KrisThielemans

cpp/helpers/include/petsird_helpers.h

@@ -101,9 +101,10 @@ make_detection_bin(const ScannerInformation& scanner, const TypeOfModule& type_o
 
 inline float
 get_detection_efficiency(const ScannerInformation& scanner, const TypeOfModulePair& type_of_module_pair,
-                         const DetectionBin& detection_bin_1, const DetectionBin& detection_bin_2)
+                         const DetectionBin& detection_bin_1, const DetectionBin& detection_bin_2,
+                         bool with_calibration_factor = true)
 {
-  float eff = 1.0F;
+  float eff = with_calibration_factor ? scanner.detection_efficiencies.calibration_factor : 1.0F;
   const auto& detection_bin_efficiencies = scanner.detection_efficiencies.detection_bin_efficiencies;
   if (detection_bin_efficiencies)
     {
@@ -141,9 +142,10 @@ get_detection_efficiency(const ScannerInformation& scanner, const TypeOfModulePa
 
 inline float
 get_detection_efficiency(const ScannerInformation& scanner, const TypeOfModulePair& type_of_module_pair,
-                         const CoincidenceEvent& event)
+                         const CoincidenceEvent& event, bool with_calibration_factor = true)
 {
-  return get_detection_efficiency(scanner, type_of_module_pair, event.detection_bins[0], event.detection_bins[1]);
+  return get_detection_efficiency(scanner, type_of_module_pair, event.detection_bins[0], event.detection_bins[1],
+                                  with_calibration_factor);
 }
 
 } // namespace petsird_helpers

cpp/helpers/include/petsird_helpers/create.h

@@ -49,6 +49,9 @@ construct_2D_nested_vector(std::size_t size0, std::size_t size1)
 
   Elements will be constructed via the default constructors, so you will still have
   to fill in the actual values.
+
+  To prevent dramatic errors, the calibration_factor is set to 1, but this factor has to be set
+  correctly afterwards.
 */
 inline void
 initialize_scanner_information_dimensions(petsird::ScannerInformation& scanner, const std::size_t num_module_types,
@@ -60,6 +63,8 @@ initialize_scanner_information_dimensions(petsird::ScannerInformation& scanner,
   scanner.energy_resolution_at_511 = construct_vector<float>(num_module_types);
 
   scanner.detection_efficiencies = petsird::DetectionEfficiencies();
+  // set this to 1 to avoid the case where the user forgets to set it, in which case it might be zero.
+  scanner.detection_efficiencies.calibration_factor = 1.F;
 
   if (allocate_detection_bin_efficiencies)
     {

cpp/helpers/petsird_analysis.cpp

@@ -135,6 +135,7 @@ main(int argc, char const* argv[])
       std::cout << "Event energy mid points: " << energy_mid_points << std::endl;
       all_energy_mid_points.push_back(energy_mid_points);
 
+      std::cout << "Calibration factor: " << scanner.detection_efficiencies.calibration_factor << std::endl;
       std::cout << "Singles Histogram Level: ";
       switch (scanner.singles_histogram_level)
         {

cpp/helpers/petsird_generator.cpp

@@ -113,6 +113,9 @@ get_scanner_geometry()
 void
 set_detection_efficiencies(petsird::ScannerInformation& scanner)
 {
+  // use some non-physical value for the calibration factor
+  scanner.detection_efficiencies.calibration_factor = 42.F;
+
   const auto num_module_types = scanner.scanner_geometry.NumberOfModuleTypes();
   // only 1 type of module in the current scanner
   assert(num_module_types == 1);

model/DetectionEfficiencies.yml

@@ -48,14 +48,31 @@ ModulePairEfficienciesVector: ModulePairEfficiencies*
 # 3. SGID = modulePairSGIDLUT[type_of_module1][type_of_module2][module_index1, module_index2]
 # 4. if (SGID < 0) return 0
 # 5. module_pair_efficiencies = modulePairEfficienciesVector[type_of_module1][type_of_module2][SGID]
-# 6. return detectionBinEfficiencies[type_of_module1](detection_bin1) * detectionBinEfficiencies[type_of_module2](detection_bin2)
+# 6. return calibrationFactor
+#       * detectionBinEfficiencies[type_of_module1](detection_bin1) * detectionBinEfficiencies[type_of_module2](detection_bin2)
 #       * module_pair_efficiencies[detection_bin_in_module1, detection_bin_in_module2]
 #
 # If either of the components is not present, its value is considered to be 1.
 #
+# The calibrationFactor is defined in terms of measured/expected "true unscattered" coincidences.
+# Moreover, it is (in principle) defined in terms of an idealised measurement of large uniform source.
+# Therefore, "expected" counts should take into account the
+# - activity (as estimated via a dose calibrator)
+# - positron emission branching fraction
+# - detectionBinEfficiencies and modulePairEfficiencies
+# - the line integral (in millimeter) through the source
+# This definition has the effect that for a given pair of detection_bins, we should have
+#   mean_precorrected_true_coincidences =
+#     detection_effiency * line_integral(number_of_positron_annihilations_in_source_as_estimated_from_dose_calibrator)
+#   (with the precorrection taking into scatter and attenuation)
+#
 # Note that computing a detection efficiency for a triple coincidence is left to the user.
 DetectionEfficiencies: !record
   fields:
+    # Text field describing how these efficiencies were determined
+    methodDescription: string
+    # Calibration factor for "absolute" efficiency
+    calibrationFactor: float
     # List of detection efficiencies for every detection bin (one for each module-type).
     # Constraint: size(detectionBinEfficiencies) == ScannerGeometry.numberOfModuleTypes()
     detectionBinEfficiencies: DetectionBinEfficiencies*?

python/petsird/helpers/__init__.py

@@ -87,28 +87,31 @@ def make_detection_bin(
 def get_detection_efficiency(scanner: petsird.ScannerInformation,
                              type_of_module_pair: petsird.TypeOfModulePair,
                              detection_bin_1: petsird.DetectionBin,
-                             detection_bin_2: petsird.DetectionBin) -> float:
+                             detection_bin_2: petsird.DetectionBin,
+                             with_calibration_factor: bool = ...) -> float:
     """Compute the detection efficiency for a pair of detectors"""
     ...
 
 
 @typing.overload
 def get_detection_efficiency(scanner: petsird.ScannerInformation,
                              type_of_module_pair: petsird.TypeOfModulePair,
-                             event: petsird.CoincidenceEvent) -> float:
+                             event: petsird.CoincidenceEvent,
+                             with_calibration_factor: bool = ...) -> float:
     """Compute the detection efficiency for a coincidence event"""
     ...
 
 
 _DetectionBinUnannotated = typing.get_args(petsird.DetectionBin)[0]
 
 
-def get_detection_efficiency(
-        scanner: petsird.ScannerInformation,
-        type_of_module_pair: petsird.TypeOfModulePair,
-        event_or_detection_bin_1: typing.Union[petsird.CoincidenceEvent,
-                                               petsird.DetectionBin],
-        detection_bin_2: petsird.DetectionBin = None) -> float:
+def get_detection_efficiency(scanner: petsird.ScannerInformation,
+                             type_of_module_pair: petsird.TypeOfModulePair,
+                             event_or_detection_bin_1: typing.Union[
+                                 petsird.CoincidenceEvent,
+                                 petsird.DetectionBin],
+                             detection_bin_2: petsird.DetectionBin = None,
+                             with_calibration_factor: bool = True) -> float:
     """Compute the detection efficiency"""
     if isinstance(event_or_detection_bin_1, _DetectionBinUnannotated):
         detection_bin_1 = event_or_detection_bin_1

python/petsird/helpers/analysis.py

@@ -75,6 +75,9 @@ def parserCreator():
                                  event_energy_bin_edges[1:]) / 2
             all_energy_mid_points.append(energy_mid_points)
             print("Event energy mid points: ", energy_mid_points)
+
+        print("Calibration factor: ",
+              scanner.detection_efficiencies.calibration_factor)
         if scanner.detection_efficiencies.module_pair_sgidlut is not None:
             for type_of_module0 in range(num_module_types):
                 for type_of_module1 in range(num_module_types):

python/petsird/helpers/generator.py

@@ -343,6 +343,8 @@ def get_detection_efficiencies(
 ) -> petsird.DetectionEfficiencies:
     """return some (non-physical) detection efficiencies"""
 
+    # use some non-physical value for the calibration factor
+    calibration_factor = 42.
     num_types_of_modules = scanner.scanner_geometry.number_of_module_types()
     # single-level list for detection_bin_efficiencies
     all_detection_bin_efficiencies = []
@@ -377,6 +379,7 @@ def get_detection_efficiencies(
             module_pair_efficiencies_vectors0)
 
     return petsird.DetectionEfficiencies(
+        calibration_factor=calibration_factor,
         detection_bin_efficiencies=all_detection_bin_efficiencies,
         module_pair_sgidlut=all_LUTs,
         module_pair_efficiencies_vectors=all_module_pair_efficiencies_vectors,