Spatial resolution 2D distribution changes

.github/workflows/main.yml

@@ -23,9 +23,9 @@ jobs:
             options: torch
 
     env:
-      ROOT_VERSION: 'v6-32-02'
+      ROOT_VERSION: 'v6-38-00'
       ITK_VERSION: 'v5.4.2'
-      GEANT4_VERSION: 'v11.3.0'
+      GEANT4_VERSION: 'v11.4.0'
       ROOT_DIR: $(HOME)/software/root
       GEANT4_DIR: $(HOME)/software/geant4
 
@@ -41,8 +41,8 @@ jobs:
       uses: actions/cache@v3
       with:
         path: ~/software
-        key: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-build1
-        restore-keys: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-build1
+        key: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-qt6-root3800
+        restore-keys: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-qt6-root3800
     - name: Install dependencies
       run: |
         if [ "${{ matrix.os }}" == "ubuntu-latest" ]; then
@@ -56,33 +56,39 @@ jobs:
                                   libxpm-dev \
                                   libxft-dev \
                                   libxext-dev \
-                                  qtbase5-dev \
-                                  qt5-qmake \
+                                  qt6-base-dev \
+                                  qt6-base-dev-tools \
+                                  qt6-tools-dev \
+                                  qt6-tools-dev-tools \
                                   git \
                                   cmake \
                                   libssl-dev \
                                   python3.8-dev \
                                   ccache \
                                   libfftw3-dev
           gcc -v
-        elif [ "${{ matrix.os }}" == "macos-latest" ]; then
+        elif [ "${{ matrix.os }}" == "macos-latest" ]; then          
+          brew update || true
           brew install python@3.12 || true
-          brew link --overwrite python@3.12
-          #brew update
-          brew cleanup
-          brew config
-          #rm -rf /usr/local/bin/python3.11-config /usr/local/bin/2to3-3.11 /usr/local/bin/idle3.11 /usr/local/bin/pydoc3.11 /usr/local/bin/python3.11
-          #rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3
-          brew install --force --verbose --overwrite --debug qt5 \
+
+          brew link --overwrite --force python@3.12
+          export Python3_EXECUTABLE="$(brew --prefix python@3.12)/bin/python3"
+          brew install --force --verbose --overwrite --debug cmake \
+                                                             qt \
                                                              ccache \
                                                              tbb \
                                                              xrootd \
                                                              fftw
-          brew link qt5 --force && sudo ln -s /usr/local/opt/qt/mkspecs /usr/local/mkspecs && sudo ln -s /usr/local/opt/qt/plugins /usr/local/plugins
-          export PATH=/usr/local/opt/qt/bin:$PATH
-          export LDFLAGS="-L/usr/local/opt/qt/lib -L/usr/local/opt/llvm/lib"
-          export CPPFLAGS="-I/usr/local/opt/qt/include -I/usr/local/opt/llvm/include -fopenmp"
-        fi
+
+          brew cleanup
+          #LLVM_PREFIX="$(brew --prefix llvm || true)"
+          #export SDKROOT="$(xcrun --sdk macosx --show-sdk-path)"
+          export CC=/usr/bin/clang
+          export CXX=/usr/bin/clang++
+
+          cmake --version
+          which cmake
+        fi  
         cd $HOME/
         mkdir -p software option_dependencies
         cmake --version
@@ -95,13 +101,14 @@ jobs:
           mkdir src bin install
           git clone  --branch $ROOT_VERSION https://github.com/root-project/root.git --depth 1 src
           cd bin
+          export CC=/usr/bin/clang
+          export CXX=/usr/bin/clang++
           cmake ../src -DCMAKE_CXX_STANDARD=17 \
-                       -Dpython=OFF \
-                       -Dpyroot=OFF \
-                       -Dclad=OFF \
-                       -Dxrootd=OFF \
+                       -DCMAKE_C_COMPILER=$CC \
+                       -DCMAKE_CXX_COMPILER=$CXX \
+                       -DCMAKE_OSX_ARCHITECTURES=arm64 \
                        -DCMAKE_INSTALL_PREFIX=$HOME/software/root/install
-          make -j4 install
+          make install -j 4
           cd ..
           rm -rf bin src
         fi
@@ -116,8 +123,10 @@ jobs:
           cd bin
           cmake -DGEANT4_INSTALL_DATA=ON \
                 -DGEANT4_BUILD_MULTITHREADED=OFF \
+                -DGEANT4_USE_QT=ON \
                 -DGEANT4_INSTALL_DATADIR=$HOME/software/geant4/data \
                 -DCMAKE_INSTALL_PREFIX=$HOME/software/geant4/install \
+                -DQt6_DIR="$Qt6_DIR" \
                 ../src
           make -j4 install
           cd ..

AUTHORS

@@ -25,6 +25,7 @@ Garcia Marie-Paule
 Garpebring Anders
 Germano Russo
 Gil Alex Vergara
+Granado-Gonzalez Marc
 Grevillot Loic
 Groiselle Corinne
 Gueth Pierre

CMakeLists.txt

@@ -300,6 +300,7 @@ IF(GATE_COMPILE_GATEDIGIT)
     INSTALL(TARGETS Convert_CCMod2PETCoinc DESTINATION bin)
 ENDIF(GATE_COMPILE_GATEDIGIT)
 
+
 #=========================================================
 # We remove the warning option "shadow", because there are tons of
 # such warning related to clhep/g4 system of units.

docs/digitizer_and_detector_modeling.rst

@@ -533,32 +533,31 @@ BEWARE: The confined and the use of the truncated Gaussian are default options.
 /gate/digitizerMgr/crystal/SinglesDigitizer/Singles/spatialResolution/confineInsideOfSmallestElement true 
 /gate/digitizerMgr/pseudoCrystal/SinglesDigitizer/Singles/spatialResolution/useTruncatedGaussian 		true
 
-
 **Configuring Spatial Resolution with 1D and 2D Distributions**::
 
-This approach is particularly essential  for  monolithic crystal detectors, where factors like edge effects and interaction positions significantly  may influence spatial  resolution. The spatial distribution resolution allows to select what axis will be using such distribution. Currently only one distribution is allowed for all axis.
+This approach is particularly essential  for  monolithic crystal detectors, where factors such as edge effects and interaction positions may significantly influence spatial resolution.
+
+In order to map the 2D distribution with the crystal, the command "nameAxis" is used. The value of the axis provided will be the ones attributed to the 2 dimensions (columns and rows) of the distribution file. Assuming the crystal will be placed along the Z direction the options are "XZ" for a ring starting on top, and "YZ" for a ring starting on the sides. The default value for "nameAxis" is "YZ".
 
 Here is an example of how to configure this in a macro file:
 
 **Example for 2D distribution**::
 
 
-The axis are selected with a "nameAxis". IMPORTANT NOTE: Only the "XY" axis are available:
-
-/gate/distributions/name                                                                                        my_distrib2D
-/gate/distributions/insert                                                                                      File
-/gate/distributions/my_distrib2D/setFileName                                                                    data/my_stddev_distribuion_file.txt
-/gate/distributions/my_distrib2D/readMatrix2d
-/gate/digitizerMgr/pseudoCrystal/SinglesDigitizer/Singles/insert                                               spatialResolution
-/gate/digitizerMgr/pseudoCrystal/SinglesDigitizer/Singles/spatialResolution/nameAxis                           XY
-/gate/digitizerMgr/pseudoCrystal/SinglesDigitizer/Singles/spatialResolution/fwhmDistrib2D                      my_distrib2D
-
+  /gate/distributions/name    my_distrib2D
+  /gate/distributions/insert   File
+  /gate/distributions/my_distrib2D/setFileName    Lut_XY.txt
+  /gate/distributions/my_distrib2D/readMatrix2d
+  /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution
+  /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/nameAxis       YZ
+  /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmYdistrib2D my_distrib2D
+  /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmZdistrib2D my_distrib2D
 
 **Example for 1D distribution**::
 
   /gate/distributions/name   my_distrib1D
   /gate/distributions/insert  File
-  /gate/distributions/my_distrib1D/setFileName  macros/LutY.txt
+  /gate/distributions/my_distrib1D/setFileName  macros/Lut_Y.txt
   /gate/distributions/my_distrib1D/read
   /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution
   /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmYdistrib my_distrib1D
@@ -569,23 +568,18 @@ These commands allow for more precise control over the spatial resolution by usi
 
 BEWARE : The file for 2D Distribution  should be structured such that:
 
--The first line contains the values of the x position. See the example below where the first line has values from -29.5 to 29.5 in this example for a 59 mm crystal.
+-The first line contains the x values.
 
--Each subsequent line begins with the value of the y position followed by the standard deviation (stddev) values corresponding to each x value and y value pair.
+-Each subsequent line begins with a y value followed by the standard deviation (stddev) values corresponding to each x value and y value pair.
 
 **Example**::
 
--29.50 -22.94 -16.39 -9.83 -3.28 3.28 9.83 16.39 22.94 29.50
--29.50 8.05 5.1 4.24 4.23 4.23 4.19 4.56 4.6 5.02 7.68
--22.94 4.76 2.39 2.31 2.35 2.4 2.35 2.36 2.32 2.44 4.52
--16.39 4.53 2.45 2.28 2.39 2.35 2.36 2.46 2.46 2.55 4.52
--9.83 4.38 2.35 2.13 2.09 2.13 2.07 2.14 2.11 2.33 4.44
--3.28 4.18 2.18 1.97 1.96 1.97 1.97 2.03 2.02 2.19 4.2
-3.28 4.2 2.27 2.04 2.01 2.03 2.06 1.98 1.97 2.24 4.13
-9.83 4.24 2.4 2.23 2.24 2.21 2.28 2.22 2.21 2.35 4.36
-16.39 4.32 2.59 2.47 2.55 2.56 2.52 2.45 2.42 2.55 4.51
-22.94 4.07 2.32 2.33 2.34 2.34 2.27 2.36 2.29 2.4 4.3
-29.50 7.2 4.58 4.07 3.9 3.89 3.83 4.09 3.99 4.47 7.08
+
+  -30 -15 0 15 30
+  -15 9.31 7.25 6.22 7.31 9.73
+  0 9.42 6.25 3.25 6.22 9.72
+  15 9.42 6.53 3.15 6.32 9.71
+  30 9.42 7.45 6.25 7.32 9.74
 
 Energy Framing
 ^^^^^^^^^^^^^^

source/digits_hits/include/GateDigi.hh

@@ -173,6 +173,13 @@ public:
       inline void SetParentID(G4int parentID)          { m_parentID = parentID; }
       inline G4int GetParentID() const                { return m_parentID; }
 
+      inline void SetSpatialRes2DStdDevX(G4double v) { m_spatialRes2DStdDevX = v; }
+      inline G4double GetSpatialRes2DStdDevX() const { return m_spatialRes2DStdDevX; }
+      inline void SetSpatialRes2DStdDevY(G4double v) { m_spatialRes2DStdDevY = v; }
+      inline G4double GetSpatialRes2DStdDevY() const { return m_spatialRes2DStdDevY; }
+      inline void SetSpatialRes2DStdDevZ(G4double v) { m_spatialRes2DStdDevZ = v; }
+      inline G4double GetSpatialRes2DStdDevZ() const { return m_spatialRes2DStdDevZ; }
+
       //AE
       inline G4double GetEnergyError()   const                  	      { return m_energyError; }
       inline void SetEnergyError(G4double value)           	      { m_energyError = value; }
@@ -252,6 +259,11 @@ public:
   G4int m_nCrystalConv;
   //--------------------
 
+  // Spatial resolution info set by spatial resolution digitizer when using 2D distributions
+  G4double m_spatialRes2DStdDevX;
+  G4double m_spatialRes2DStdDevY;
+  G4double m_spatialRes2DStdDevZ;
+
 
   //! Pointer to the original crystal hit if known
   const void* m_mother;

source/digits_hits/include/GateSpatialResolution.hh

@@ -47,11 +47,15 @@ public:
 
   //! These functions return the resolution in use.
     G4double GetFWHM()				{ return m_fwhm; }
-    GateVDistribution* GetFWHMxdistrib()	{ return m_fwhmXdistrib; }
-    GateVDistribution* GetFWHMydistrib()    	{ return m_fwhmYdistrib; }
-    GateVDistribution* GetFWHMzdistrib()    	{ return m_fwhmZdistrib; }
+    GateVDistribution* GetFWHMxdistrib() 	{ return m_fwhmXDistrib; }
+    GateVDistribution* GetFWHMydistrib()    	{ return m_fwhmYDistrib; }
+    GateVDistribution* GetFWHMzdistrib()    	{ return m_fwhmZDistrib; }
     G4String GetNameAxis()				   { return m_nameAxis;     }
-    GateVDistribution* GetFWHMDistrib2D()	{ return m_fwhmDistrib2D; }
+    GateVDistribution* GetFWHMDistrib2D()	{ return m_fwhmXDistrib2D ? m_fwhmXDistrib2D : (m_fwhmYDistrib2D? m_fwhmYDistrib2D : m_fwhmZDistrib2D); }
+
+    GateVDistribution* GetFWHMXDistrib2D() { return m_fwhmXDistrib2D; }
+    GateVDistribution* GetFWHMYDistrib2D() { return m_fwhmYDistrib2D; }
+    GateVDistribution* GetFWHMZDistrib2D() { return m_fwhmZDistrib2D; }
 
     G4double GetFWHMx()         { return m_fwhmX; }
     G4double GetFWHMy()			{ return m_fwhmY; }
@@ -62,13 +66,18 @@ public:
       If you want a resolution of 10%, SetSpresolution(0.1)
     */
     void SetFWHM(G4double val)   { m_fwhm = val;  }
-    void SetFWHMxdistrib(GateVDistribution* dist)  { m_fwhmXdistrib= dist; }
-    void SetFWHMydistrib(GateVDistribution* dist)  { m_fwhmYdistrib = dist; }
-    void SetFWHMzdistrib(GateVDistribution* dist)  { m_fwhmZdistrib = dist; }
+    void SetFWHMxdistrib(GateVDistribution* dist)  { m_fwhmXDistrib= dist; }
+    void SetFWHMydistrib(GateVDistribution* dist)  { m_fwhmYDistrib = dist; }
+    void SetFWHMzdistrib(GateVDistribution* dist)  { m_fwhmZDistrib = dist; }
 
 
     void SetNameAxis(const G4String& name) {m_nameAxis=name;}
-    void SetFWHMDistrib2D(GateVDistribution* dist)  { m_fwhmDistrib2D= dist;}
+    // Backwards-compatible: sets all 3 axis 2D distributions to the same distribution
+    void SetFWHMDistrib2D(GateVDistribution* dist)  { m_fwhmXDistrib2D = m_fwhmYDistrib2D = m_fwhmZDistrib2D = dist; }
+
+    void SetFWHMXDistrib2D(GateVDistribution* dist) { m_fwhmXDistrib2D = dist; }
+    void SetFWHMYDistrib2D(GateVDistribution* dist) { m_fwhmYDistrib2D = dist; }
+    void SetFWHMZDistrib2D(GateVDistribution* dist) { m_fwhmZDistrib2D = dist; }
 
     void SetFWHMx(G4double val)   { m_fwhmX = val;  }
     void SetFWHMy(G4double val)   { m_fwhmY = val;  }
@@ -97,11 +106,13 @@ protected:
     G4double m_fwhmY;
     G4double m_fwhmZ;
 
-    GateVDistribution* m_fwhmXdistrib;
-    GateVDistribution* m_fwhmYdistrib;
-    GateVDistribution* m_fwhmZdistrib;
+    GateVDistribution* m_fwhmXDistrib;
+    GateVDistribution* m_fwhmYDistrib;
+    GateVDistribution* m_fwhmZDistrib;
 
-    GateVDistribution* m_fwhmDistrib2D;
+    GateVDistribution* m_fwhmXDistrib2D;
+    GateVDistribution* m_fwhmYDistrib2D;
+    GateVDistribution* m_fwhmZDistrib2D;
 
     G4String m_nameAxis;
 

source/digits_hits/include/GateSpatialResolutionMessenger.hh

@@ -48,7 +48,9 @@ private:
   G4UIcmdWithAString	*spresolutionYdistribCmd;// Command declaration for 1D Y-resolution distribution
   G4UIcmdWithAString	*spresolutionZdistribCmd;// Command declaration for 1D Y-resolution distribution
   G4UIcmdWithAString	*nameAxisCmd;
-  G4UIcmdWithAString	*spresolutionDistrib2DCmd; // Command declaration for 2D resolution distribution will be used with nameAxisCmd
+  G4UIcmdWithAString	*spresolutionXDistrib2DCmd; // per-axis 2D distribution commands will be used nameAxis to determine which axes to apply the distribution to
+  G4UIcmdWithAString	*spresolutionYDistrib2DCmd;
+  G4UIcmdWithAString	*spresolutionZDistrib2DCmd;
   G4UIcmdWithABool* 	confineCmd;
   G4UIcmdWithABool*		useTruncatedGaussianCmd;
 

source/digits_hits/src/GateDigi.cc

@@ -42,6 +42,9 @@ GateDigi::GateDigi(const void* itsMother):
       m_energyError(0.0),
       m_globalPosError(0.0),
       m_localPosError(0.0),
+      m_spatialRes2DStdDevX(0.0),
+      m_spatialRes2DStdDevY(0.0),
+      m_spatialRes2DStdDevZ(0.0),
       m_mother(itsMother)
 {
 }

source/digits_hits/src/GateDistributionTruncatedGaussian.cc

@@ -21,6 +21,7 @@ See LICENSE.md for further details
 
 //#include "GateDistributionTruncatedGaussianMessenger.hh"
 #include <math.h>
+#include <algorithm>
 //#include <CLHEP/Random/RandGauss.h>
 #include "Randomize.hh"
 #include "GateTools.hh"
@@ -144,7 +145,24 @@ G4double GateDistributionTruncatedGaussian::computeTruncatedSigmaStatic(G4double
     double mean_shift = (phi_lowLim - phi_hiLim) / Z;
     double variance_correction = 1 - (lowLim_std * phi_lowLim - hiLim_std * phi_hiLim) / Z - mean_shift * mean_shift;
 
-    return sigma * sqrt(variance_correction);
+    double truncatedSigma = sigma * sqrt(variance_correction);
+
+    // Edge correction: empirical fit to compensate for edge effects
+    constexpr double A = 3.61070188;
+    constexpr double B = -0.86538264;
+
+    double delta = std::min(mu - lowLimit, highLimit - mu);
+    delta = std::max(0., delta);
+
+    double delta_norm = delta / sigma;
+    double correction = 1;
+
+    // Cuts on delta_norm > 0.95 and delta_norm < 4 comply with the safe limits of the fit
+    delta_norm = std::max(.95, delta_norm);
+    if (delta_norm < 4)
+        correction = A * std::exp(B * delta_norm);
+
+    return truncatedSigma * correction;
 }
 
 //___________________________________________________________________
@@ -162,5 +180,22 @@ G4double GateDistributionTruncatedGaussian::computeTruncatedSigma() const{
     double mean_shift = (phi_lowLim - phi_hiLim) / Z;
     double variance_correction = 1 - (lowLim_std * phi_lowLim - hiLim_std * phi_hiLim) / Z - mean_shift * mean_shift;
 
-    return m_Sigma * sqrt(variance_correction);
+    double truncatedSigma = m_Sigma * sqrt(variance_correction);
+
+    // Edge correction: empirical fit to compensate for edge effects
+    constexpr double A = 3.61070188;
+    constexpr double B = -0.86538264;
+
+    double delta = std::min(m_Mu - m_lowLimit, m_highLimit - m_Mu);
+    delta = std::max(0., delta);
+
+    double delta_norm = delta / m_Sigma;
+    double correction = 1;
+
+    // Cuts on delta_norm > 0.95 and delta_norm < 4 comply with the safe limits of the fit
+    delta_norm = std::max(.95, delta_norm);
+    if (delta_norm < 4)
+        correction = A * std::exp(B * delta_norm);
+
+    return truncatedSigma * correction;
 }