diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index be995ba51..975acea62 100644 --- a/.github/workflows/main.yml +++ b/.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,8 +56,10 @@ 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 \ @@ -65,24 +67,28 @@ jobs: 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 .. diff --git a/AUTHORS b/AUTHORS index 102b936ab..87495aa9a 100644 --- a/AUTHORS +++ b/AUTHORS @@ -25,6 +25,7 @@ Garcia Marie-Paule Garpebring Anders Germano Russo Gil Alex Vergara +Granado-Gonzalez Marc Grevillot Loic Groiselle Corinne Gueth Pierre diff --git a/CMakeLists.txt b/CMakeLists.txt index e009bc4b5..395a2798c 100644 --- a/CMakeLists.txt +++ b/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. diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 9b6f03e8a..d4d254131 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/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 ^^^^^^^^^^^^^^ diff --git a/source/digits_hits/include/GateDigi.hh b/source/digits_hits/include/GateDigi.hh index cf513cb44..95908f6df 100644 --- a/source/digits_hits/include/GateDigi.hh +++ b/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; diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh index 8fd30748a..549c19f44 100644 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/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; diff --git a/source/digits_hits/include/GateSpatialResolutionMessenger.hh b/source/digits_hits/include/GateSpatialResolutionMessenger.hh index 33d33213a..bfe534c3a 100644 --- a/source/digits_hits/include/GateSpatialResolutionMessenger.hh +++ b/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; diff --git a/source/digits_hits/src/GateDigi.cc b/source/digits_hits/src/GateDigi.cc index d6bea9e2f..78a09ad5d 100644 --- a/source/digits_hits/src/GateDigi.cc +++ b/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) { } diff --git a/source/digits_hits/src/GateDistributionTruncatedGaussian.cc b/source/digits_hits/src/GateDistributionTruncatedGaussian.cc index 03d22c17e..fae837e99 100644 --- a/source/digits_hits/src/GateDistributionTruncatedGaussian.cc +++ b/source/digits_hits/src/GateDistributionTruncatedGaussian.cc @@ -21,6 +21,7 @@ See LICENSE.md for further details //#include "GateDistributionTruncatedGaussianMessenger.hh" #include +#include //#include #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; } \ No newline at end of file diff --git a/source/digits_hits/src/GateSpatialResolution.cc b/source/digits_hits/src/GateSpatialResolution.cc index f5e795ff0..6c0b3573e 100644 --- a/source/digits_hits/src/GateSpatialResolution.cc +++ b/source/digits_hits/src/GateSpatialResolution.cc @@ -23,9 +23,16 @@ This blurring has been validated up to a given FWHM of 10mm. At higher FWHM, the number of "relocated" digis is no longer negligible. The blurring effect is then so compensated that resolution will improve compared to lower values of FWHM. -modified by Radia Oudihat 06/2024 - Added support for 1D and 2D FWHM distributions for X and Y, and applied Gaussian blurring. + Added support for 1D FWHM distribution for X and Y, and applied Gaussian blurring. Implemented logic to determine standard deviations (stddevX, stddevY) based on defined 1D and 2D FWHM distributions for the X and Y axes. -*/ +-modified by Marc Granado-Gonzalez 2025 + - Added 2D FWHM distribution for X, Y and Z, and applied Gaussian blurring. + Implemented logic to determine standard deviations (stddevX, stddevY, stddevZ) based on defined 2D FWHM distributions for the X, Y and Z axes. + - Added option to choose the axis pair for 2D distributions (nameAxis): "XZ" or "YZ" (default "YZ"). + - Added Truncated Gaussian option for confined and non-confined cases. + */ + + #include "GateSpatialResolution.hh" #include "GateSpatialResolutionMessenger.hh" #include "GateDigi.hh" @@ -55,11 +62,13 @@ GateSpatialResolution::GateSpatialResolution(GateSinglesDigitizer *digitizer, G4 m_fwhmX(0), m_fwhmY(0), m_fwhmZ(0), - m_fwhmXdistrib(0), - m_fwhmYdistrib(0), - m_fwhmZdistrib(0), - m_nameAxis("XY"), - m_fwhmDistrib2D(0), + m_fwhmXDistrib(0), + m_fwhmYDistrib(0), + m_fwhmZDistrib(0), + m_nameAxis("YZ"), + m_fwhmXDistrib2D(0), + m_fwhmYDistrib2D(0), + m_fwhmZDistrib2D(0), m_IsConfined(true), m_UseTruncatedGaussian(true), m_Navigator(0), @@ -83,38 +92,39 @@ GateSpatialResolution::~GateSpatialResolution() } void GateSpatialResolution::SetSpatialResolutionParameters() { // Check FWHM parameters - if (m_fwhm != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmZ != 0 || m_fwhmDistrib2D != 0 || m_fwhmXdistrib !=0 || m_fwhmYdistrib !=0 || m_fwhmZdistrib !=0 )) { + if (m_fwhm != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmZ != 0 || m_fwhmXDistrib2D != 0 || m_fwhmYDistrib2D != 0 || m_fwhmZDistrib2D != 0 || m_fwhmXDistrib !=0 || m_fwhmYDistrib !=0 || m_fwhmZDistrib !=0 )) { G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set a unique FWHM for all 3 axes OR set FWHM for X, Y, Z individually." << G4endl; abort(); } - if (m_fwhmDistrib2D) - { - if ((m_fwhmX != 0 || m_fwhmXdistrib !=0) && (m_nameAxis.find('X') != std::string::npos)) { - G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM for X distribution." << G4endl; - abort(); - } - - if ((m_fwhmY != 0|| m_fwhmYdistrib !=0) && (m_nameAxis.find('Y') != std::string::npos)) { - G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM for Y distribution." << G4endl; - abort(); - } - - if ((m_fwhmZ != 0 || m_fwhmZdistrib !=0) && (m_nameAxis.find('Z') != std::string::npos)) { - G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Z OR set FWHM for Z distribution." << G4endl; - abort(); - } - } + if (m_fwhmXDistrib2D || m_fwhmYDistrib2D || m_fwhmZDistrib2D) + { + // If a per-axis 2D distribution is provided for an axis, it is + // ambiguous to also provide a scalar FWHM or a 1D distribution for + // the same axis. Check per-axis instead of using nameAxis combinatorics. + if (m_fwhmXDistrib2D && (m_fwhmX != 0 || m_fwhmXDistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM distribution for X." << G4endl; + abort(); + } + if (m_fwhmYDistrib2D && (m_fwhmY != 0 || m_fwhmYDistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM distribution for Y." << G4endl; + abort(); + } + if (m_fwhmZDistrib2D && (m_fwhmZ != 0 || m_fwhmZDistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Z OR set FWHM distribution for Z." << G4endl; + abort(); + } + } - if (m_fwhmX != 0 && m_fwhmXdistrib !=0){ + if (m_fwhmX != 0 && m_fwhmXDistrib !=0){ G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM for Z distribution." << G4endl; abort(); } - if (m_fwhmY != 0 && m_fwhmYdistrib !=0){ + if (m_fwhmY != 0 && m_fwhmYDistrib !=0){ G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM for Z distribution." << G4endl; abort(); } - if (m_fwhmZ != 0 && m_fwhmZdistrib !=0){ + if (m_fwhmZ != 0 && m_fwhmZDistrib !=0){ G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Z OR set FWHM for Z distribution." << G4endl; abort(); } @@ -196,40 +206,30 @@ void GateSpatialResolution::Digitize(){ G4double Px = P.x(); G4double Py = P.y(); G4double Pz = P.z(); - G4double stddevX, stddevY, stddevZ; - - if (m_fwhmDistrib2D) { - if (m_nameAxis.size() != 2) { - GateError(" *** ERROR*** GateSpatialResolution::Digitize. " - "Attempt to use fwhmDistrib2D but the length of the named axis is not 2!\n"); - } else { - if (m_nameAxis == "XY") { - stddevX = m_fwhmDistrib2D->Value2D(P.x() * mm, P.y() * mm); - stddevY = m_fwhmDistrib2D->Value2D(P.x() * mm, P.y() * mm); - if (fwhmZ) stddevZ = fwhmZ / GateConstants::fwhm_to_sigma; - } else if (m_nameAxis == "XZ") { - stddevX = m_fwhmDistrib2D->Value2D(P.x() * mm, P.z() * mm); - stddevZ = m_fwhmDistrib2D->Value2D(P.x() * mm, P.z() * mm); - if (fwhmY) stddevY = fwhmY / GateConstants::fwhm_to_sigma; - } else if (m_nameAxis == "YZ") { - stddevY = m_fwhmDistrib2D->Value2D(P.y() * mm, P.z() * mm); - stddevZ = m_fwhmDistrib2D->Value2D(P.y() * mm, P.z() * mm); - if (fwhmX) stddevX = fwhmX / GateConstants::fwhm_to_sigma; - } else { - GateError(" *** ERROR*** GateSpatialResolution::Digitize. " - "Unrecognized axis configuration: " + m_nameAxis + "\n"); - } - } + G4double stddevX = 0., stddevY = 0., stddevZ = 0.; + + // Use the configured axis pair (m_nameAxis) to evaluate Value2D. + // Allowed pairs for PET: "XZ" or "YZ" (default "YZ"). + if (m_fwhmXDistrib2D || m_fwhmYDistrib2D || m_fwhmZDistrib2D) { + if (m_nameAxis == "XZ") { + if (m_fwhmXDistrib2D) stddevX = m_fwhmXDistrib2D->Value2D(P.x() * mm, P.z() * mm); + if (m_fwhmYDistrib2D) stddevY = m_fwhmYDistrib2D->Value2D(P.x() * mm, P.z() * mm); + if (m_fwhmZDistrib2D) stddevZ = m_fwhmZDistrib2D->Value2D(P.x() * mm, P.z() * mm); + } else { // YZ + if (m_fwhmXDistrib2D) stddevX = m_fwhmXDistrib2D->Value2D(P.y() * mm, P.z() * mm); + if (m_fwhmYDistrib2D) stddevY = m_fwhmYDistrib2D->Value2D(P.y() * mm, P.z() * mm); + if (m_fwhmZDistrib2D) stddevZ = m_fwhmZDistrib2D->Value2D(P.y() * mm, P.z() * mm); + } } else { - if (m_fwhmXdistrib) stddevX = m_fwhmXdistrib->Value(P.x() * mm); + if (m_fwhmXDistrib) stddevX = m_fwhmXDistrib->Value(P.x() * mm); else if (fwhmX) stddevX = fwhmX / GateConstants::fwhm_to_sigma; - if (m_fwhmYdistrib) stddevY = m_fwhmYdistrib->Value(P.y() * mm); + if (m_fwhmYDistrib) stddevY = m_fwhmYDistrib->Value(P.y() * mm); else if (fwhmY) stddevY = fwhmY / GateConstants::fwhm_to_sigma; - if (m_fwhmZdistrib) stddevZ = m_fwhmZdistrib->Value(P.z() * mm); + if (m_fwhmZDistrib) stddevZ = m_fwhmZDistrib->Value(P.z() * mm); else if (fwhmZ) stddevZ = fwhmZ / GateConstants::fwhm_to_sigma; } @@ -240,6 +240,11 @@ void GateSpatialResolution::Digitize(){ G4double PyNew ; G4double PzNew ; +// store the computed stddevs into the digi for later ROOT output + m_outputDigi->SetSpatialRes2DStdDevX(stddevX); + m_outputDigi->SetSpatialRes2DStdDevY(stddevY); + m_outputDigi->SetSpatialRes2DStdDevZ(stddevZ); + if (m_IsConfined) @@ -340,7 +345,13 @@ void GateSpatialResolution::Digitize(){ if (nVerboseLevel>1) G4cout << "[GateSpatialResolution::Digitize]: input digi collection is null -> nothing to do\n\n"; return; - } + // Ensure the chosen axis configuration is allowed for PET scanners + if (!(m_nameAxis == "XZ" || m_nameAxis == "YZ")) { + G4cout << "***ERROR*** GateSpatialResolution::SetSpatialResolutionParameters: " + "Only 'XZ' and 'YZ' are allowed as nameAxis values for 2D spatial resolution distributions.\n"; + abort(); + } + } StoreDigiCollection(m_OutputDigiCollection); } diff --git a/source/digits_hits/src/GateSpatialResolutionMessenger.cc b/source/digits_hits/src/GateSpatialResolutionMessenger.cc index 38e2711d7..413b271e2 100644 --- a/source/digits_hits/src/GateSpatialResolutionMessenger.cc +++ b/source/digits_hits/src/GateSpatialResolutionMessenger.cc @@ -57,15 +57,21 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol spresolutionZdistribCmd = new G4UIcmdWithAString(cmdName,this); spresolutionZdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmDistrib2D"; - spresolutionDistrib2DCmd = new G4UIcmdWithAString(cmdName,this); - spresolutionDistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring the name of the axis must be defined. default=XY"); + cmdName = GetDirectoryName() + "fwhmXDistrib2D"; + spresolutionXDistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionXDistrib2DCmd->SetGuidance("Set the 2D distribution for X axis (expects a 2D distribution object)"); + cmdName = GetDirectoryName() + "fwhmYDistrib2D"; + spresolutionYDistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionYDistrib2DCmd->SetGuidance("Set the 2D distribution for Y axis (expects a 2D distribution object)"); + cmdName = GetDirectoryName() + "fwhmZDistrib2D"; + spresolutionZDistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionZDistrib2DCmd->SetGuidance("Set the 2D distribution for Z axis (expects a 2D distribution object)"); - cmdName = GetDirectoryName()+"nameAxis"; - nameAxisCmd = new G4UIcmdWithAString(cmdName,this); - nameAxisCmd ->SetGuidance("Provide the name of the axis that will follow the spatial resolution distribution, X, Y, Z, XY, XZ, or YZ"); - nameAxisCmd ->SetCandidates("X Y Z XY XZ YZ"); + cmdName = GetDirectoryName()+"nameAxis"; + nameAxisCmd = new G4UIcmdWithAString(cmdName,this); + nameAxisCmd ->SetGuidance("Provide the coordinate pair used by 2D distributions: 'XZ' or 'YZ'. Default is 'YZ'."); + nameAxisCmd ->SetCandidates("XZ YZ"); cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; @@ -93,8 +99,12 @@ GateSpatialResolutionMessenger::~GateSpatialResolutionMessenger() delete spresolutionYdistribCmd; delete spresolutionZdistribCmd; + delete spresolutionXDistrib2DCmd; + delete spresolutionYDistrib2DCmd; + delete spresolutionZDistrib2DCmd; + delete nameAxisCmd; - delete spresolutionDistrib2DCmd; + delete confineCmd; delete useTruncatedGaussianCmd; @@ -123,13 +133,29 @@ void GateSpatialResolutionMessenger::SetNewValue(G4UIcommand * aCommand,G4String // Handle command for 2D-distribution resolution if (aCommand == nameAxisCmd) - { - m_SpatialResolution->SetNameAxis(newValue); - } - else if (aCommand == spresolutionDistrib2DCmd) - {GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); - if (distrib) m_SpatialResolution->SetFWHMDistrib2D(distrib); - } + { + // Only accept the PET-relevant options + if (newValue == "XZ" || newValue == "YZ") { + m_SpatialResolution->SetNameAxis(newValue); + } else { + G4cout << "***ERROR*** GateSpatialResolution::SetNewValue: nameAxis must be 'XZ' or 'YZ'\n"<< G4endl; + } + } + else if (aCommand == spresolutionXDistrib2DCmd) + { + GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib) m_SpatialResolution->SetFWHMXDistrib2D(distrib); + } + else if (aCommand == spresolutionYDistrib2DCmd) + { + GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib) m_SpatialResolution->SetFWHMYDistrib2D(distrib); + } + else if (aCommand == spresolutionZDistrib2DCmd) + { + GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib) m_SpatialResolution->SetFWHMZDistrib2D(distrib); + } else if ( aCommand==spresolutionXCmd ) { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } diff --git a/source/digits_hits/src/GateVirtualSegmentationSD.cc b/source/digits_hits/src/GateVirtualSegmentationSD.cc index 73107e648..125be63c5 100644 --- a/source/digits_hits/src/GateVirtualSegmentationSD.cc +++ b/source/digits_hits/src/GateVirtualSegmentationSD.cc @@ -426,7 +426,7 @@ void GateVirtualSegmentationSD::SetParameters() - if(digi_SpatialResolution->GetFWHMxdistrib()||digi_SpatialResolution->GetFWHMydistrib()||digi_SpatialResolution->GetFWHMDistrib2D()) + if(digi_SpatialResolution->GetFWHMxdistrib()||digi_SpatialResolution->GetFWHMydistrib()||digi_SpatialResolution->GetFWHMXDistrib2D()||digi_SpatialResolution->GetFWHMYDistrib2D()||digi_SpatialResolution->GetFWHMZDistrib2D()) { GateError("***ERROR*** No value of the target pitch has been provided and no value can be obtained from the spatial resolution distribution. /n Please provide a value for the pitch that is at least half of the minimum value of the distribution. "); } diff --git a/source/digits_hits/test/TestSpatialResStdDev.cc b/source/digits_hits/test/TestSpatialResStdDev.cc new file mode 100644 index 000000000..220cd958b --- /dev/null +++ b/source/digits_hits/test/TestSpatialResStdDev.cc @@ -0,0 +1,38 @@ +/* + Simple unit test to verify GateDigi stores spatial-resolution stddevs + This executable is optional and requires a Geant4-enabled build. +*/ +#include "GateDigi.hh" +#include +#include + +int main() +{ + GateDigi digi; + + const double x = 1.2345; + const double y = 2.3456; + const double z = 3.4567; + + digi.SetSpatialRes2DStdDevX(x); + digi.SetSpatialRes2DStdDevY(y); + digi.SetSpatialRes2DStdDevZ(z); + + const double eps = 1e-12; + + if (std::fabs(digi.GetSpatialRes2DStdDevX() - x) > eps) { + std::cerr << "Get/Set mismatch for StdDevX: got " << digi.GetSpatialRes2DStdDevX() << " expected " << x << "\n"; + return 1; + } + if (std::fabs(digi.GetSpatialRes2DStdDevY() - y) > eps) { + std::cerr << "Get/Set mismatch for StdDevY: got " << digi.GetSpatialRes2DStdDevY() << " expected " << y << "\n"; + return 2; + } + if (std::fabs(digi.GetSpatialRes2DStdDevZ() - z) > eps) { + std::cerr << "Get/Set mismatch for StdDevZ: got " << digi.GetSpatialRes2DStdDevZ() << " expected " << z << "\n"; + return 3; + } + + std::cout << "TestSpatialResStdDev: OK\n"; + return 0; +} diff --git a/source/general/include/GateRootDefs.hh b/source/general/include/GateRootDefs.hh index 747095c08..14727a544 100644 --- a/source/general/include/GateRootDefs.hh +++ b/source/general/include/GateRootDefs.hh @@ -338,6 +338,10 @@ class GateRootSingleBuffer Float_t energyIni; Int_t septalNb; //!< HDS : septal penetration + // Spatial resolution standard deviations (from 2D FWHM distribution) in mm + Float_t spatialRes2DStdDevX; + Float_t spatialRes2DStdDevY; + Float_t spatialRes2DStdDevZ; //@} }; @@ -434,6 +438,11 @@ class GateRootCoincBuffer Char_t comptonVolumeName2[40]; Char_t RayleighVolumeName2[40]; + // Spatial resolution standard deviations (averaged from both digis) in mm + Float_t spatialRes2DStdDevX; + Float_t spatialRes2DStdDevY; + Float_t spatialRes2DStdDevZ; + Float_t sinogramTheta; Float_t sinogramS; //@} diff --git a/source/general/src/GateRootDefs.cc b/source/general/src/GateRootDefs.cc index 140a8cfd3..bfc90339d 100644 --- a/source/general/src/GateRootDefs.cc +++ b/source/general/src/GateRootDefs.cc @@ -462,6 +462,11 @@ void GateRootSingleBuffer::Clear() // HDS : septal septalNb = 0; + + // initialize spatial resolution stddev fields + spatialRes2DStdDevX = 0.; + spatialRes2DStdDevY = 0.; + spatialRes2DStdDevZ = 0.; for ( d = 0 ; d < ROOT_VOLUMEIDSIZE ; ++d ) volumeID[d] = -1; @@ -490,6 +495,10 @@ void GateRootSingleBuffer::Clear() nCrystalCompt=0; nCrystalRayl=0; + spatialRes2DStdDevX = 0.; + spatialRes2DStdDevY = 0.; + spatialRes2DStdDevZ = 0.; + } @@ -551,6 +560,11 @@ void GateRootSingleBuffer::Fill(GateDigi* aDigi) } + // spatial resolution stddevs stored by digitizer (in internal units), convert to mm for output + spatialRes2DStdDevX = static_cast(aDigi->GetSpatialRes2DStdDevX()/mm); + spatialRes2DStdDevY = static_cast(aDigi->GetSpatialRes2DStdDevY()/mm); + spatialRes2DStdDevZ = static_cast(aDigi->GetSpatialRes2DStdDevZ()/mm); + aDigi->GetVolumeID().StoreDaughterIDs(volumeID,ROOT_VOLUMEIDSIZE); } @@ -630,8 +644,13 @@ void GateSingleTree::Init(GateRootSingleBuffer& buffer) } - //Initialized by default.TO DO: Mask option should be included or a flag - Branch("volumeID", (void *)buffer.volumeID,"volumeID[10]/I"); + // spatial resolution stddevs (2D distribution) - in mm + Branch("spatialRes2DStdDevX", &buffer.spatialRes2DStdDevX, "spatialRes2DStdDevX/F"); + Branch("spatialRes2DStdDevY", &buffer.spatialRes2DStdDevY, "spatialRes2DStdDevY/F"); + Branch("spatialRes2DStdDevZ", &buffer.spatialRes2DStdDevZ, "spatialRes2DStdDevZ/F"); + + //Initialized by default.TO DO: Mask option should be included or a flag + Branch("volumeID", (void *)buffer.volumeID,"volumeID[10]/I"); } @@ -682,6 +701,11 @@ void GateRootCoincBuffer::Clear() RayleighCrystal2 = -1; strcpy (comptonVolumeName2," "); strcpy (RayleighVolumeName2," "); + + // initialize spatial resolution stddev fields for coincidences + spatialRes2DStdDevX = 0.; + spatialRes2DStdDevY = 0.; + spatialRes2DStdDevZ = 0.; } @@ -736,6 +760,21 @@ void GateRootCoincBuffer::Fill(GateCoincidenceDigi* aDigi) strcpy (comptonVolumeName2,((aDigi->GetDigi(1))->GetComptonVolumeName()).c_str()); strcpy (RayleighVolumeName2,((aDigi->GetDigi(1))->GetRayleighVolumeName()).c_str()); + // spatial resolution stddevs: average the two constituent digis and convert to mm for output + { + G4double sX1 = (aDigi->GetDigi(0))->GetSpatialRes2DStdDevX(); + G4double sX2 = (aDigi->GetDigi(1))->GetSpatialRes2DStdDevX(); + spatialRes2DStdDevX = static_cast(((sX1 + sX2) / 2.0) / mm); + + G4double sY1 = (aDigi->GetDigi(0))->GetSpatialRes2DStdDevY(); + G4double sY2 = (aDigi->GetDigi(1))->GetSpatialRes2DStdDevY(); + spatialRes2DStdDevY = static_cast(((sY1 + sY2) / 2.0) / mm); + + G4double sZ1 = (aDigi->GetDigi(0))->GetSpatialRes2DStdDevZ(); + G4double sZ2 = (aDigi->GetDigi(1))->GetSpatialRes2DStdDevZ(); + spatialRes2DStdDevZ = static_cast(((sZ1 + sZ2) / 2.0) / mm); + } + sinogramTheta = ComputeSinogramTheta(); sinogramS = ComputeSinogramS(); } @@ -864,6 +903,11 @@ void GateCoincTree::Init(GateRootCoincBuffer& buffer) if ( GateCoincidenceDigi::GetCoincidenceASCIIMask(19) ) Branch("sinogramS", &buffer.sinogramS,"sinogramS/F"); + // spatial resolution stddevs (averaged per coincidence) - in mm + Branch("spatialRes2DStdDevX", &buffer.spatialRes2DStdDevX, "spatialRes2DStdDevX/F"); + Branch("spatialRes2DStdDevY", &buffer.spatialRes2DStdDevY, "spatialRes2DStdDevY/F"); + Branch("spatialRes2DStdDevZ", &buffer.spatialRes2DStdDevZ, "spatialRes2DStdDevZ/F"); + if ( GateCoincidenceDigi::GetCoincidenceASCIIMask(20) ) Branch("comptVolName1", (void *)buffer.comptonVolumeName1,"comptVolName1/C"); if ( GateCoincidenceDigi::GetCoincidenceASCIIMask(20) )