From 9e84260824d2cc94161bfa7b8a6577d267e9fd40 Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Tue, 2 Jul 2024 09:03:01 +0000 Subject: [PATCH 01/22] Update digitizer_and_detector_modeling.rst Add spatial resolution 1D and 2D --- docs/digitizer_and_detector_modeling.rst | 50 ++++++++++++++++++++++-- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 3dd0326f0..5bd17f4d8 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -225,7 +225,7 @@ Five types of distribution are available in GATE, namely: * Gaussian distributions, defined by a mean value and a standard deviation. * Exponential distributions, defined by its power. * Manual distributions, defined by a discrete set of points specified in the GATE macro file. The data are linearly interpolated to define the function in a continuous range. -* File distribution, acting as the manual distribution, but where the points are defined in a separate ASCII file, whose name is given as a parameter. This method is appropriate for large numbers of points and allows to describe any distribution in a totally generic way. +* File distribution, acting as the manual distribution, but where the points are defined in a separate ASCII file, whose name is given as a parameter. This method is appropriate for large numbers of points and allows to describe any distribution in a totally generic way. Additionally, GATE supports reading 2D distributions from ASCII files where values are organized in matrices. A distribution is declared by specifying its name then by creating a new instance, with its type name:: @@ -297,7 +297,8 @@ The possible type name available corresponds to the five distributions described +----------------+--------------------------------------------------------------------------------+ | read | do read the file (should be called after specifying all the other parameters) | +----------------+--------------------------------------------------------------------------------+ - + | ReadMatrix2d | do read a data file that organizes its contents in a 2D matrix format | + +----------------+--------------------------------------------------------------------------------+ Singles Digitizers ------------------- @@ -518,11 +519,52 @@ In case if the position obtained after applying a Gaussian blurring exceeds the BEWARE: This relocation procedure is validated only for the first group level of crystals. **Example**:: - + /gate/digitizerMgr/crystal/SinglesDigitizer/Singles/insert spatialResolution /gate/digitizerMgr/crystal/SinglesDigitizer/Singles/spatialResolution/fwhm 1.0 mm /gate/digitizerMgr/crystal/SinglesDigitizer/Singles/spatialResolution/confineInsideOfSmallestElement 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 influence spatial resolution. +Here is an example of how to configure this in a macro file: + +**Example**:: + + + /gate/distributions/name my_distrib2D + /gate/distributions/insert File + /gate/distributions/my_distrib2D/setFileName Lut(X,Y).txt + /gate/distributions/my_distrib2D/readMatrix2d + + /gate/distributions/name my_distrib1D + /gate/distributions/insert File + /gate/distributions/my_distrib1D/setFileName macros/LutY.txt + /gate/distributions/my_distrib1D/read + + +Configure the spatial resolution for the digitizer using these distributions:: + +/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution +/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmXdistrib2D my_distrib2D +/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmYdistrib my_distrib1D +/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/confineInsideOfSmallestElement true + +These commands allow for more precise control over the spatial resolution by using predefined distributions for the X and Y axes. + +BEWARE : The file for 2D Distribution should be structured such that: + +-The first line contains the x values. + +-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 -28.50 -27.50 +-29.50 9.62 13.66 10.22 +-28.50 11.38 11.18 10.23 +-27.50 12.82 10.43 9.70 + Energy Framing ^^^^^^^^^^^^^^ *Previously Thresholder and Upholder* From 26287e6c741f193aeb019ea1c708ff631154d066 Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Tue, 2 Jul 2024 14:43:17 +0200 Subject: [PATCH 02/22] Update digitizer_and_detector_modeling.rst --- docs/digitizer_and_detector_modeling.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 5bd17f4d8..e3dfdcf22 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -225,7 +225,7 @@ Five types of distribution are available in GATE, namely: * Gaussian distributions, defined by a mean value and a standard deviation. * Exponential distributions, defined by its power. * Manual distributions, defined by a discrete set of points specified in the GATE macro file. The data are linearly interpolated to define the function in a continuous range. -* File distribution, acting as the manual distribution, but where the points are defined in a separate ASCII file, whose name is given as a parameter. This method is appropriate for large numbers of points and allows to describe any distribution in a totally generic way. Additionally, GATE supports reading 2D distributions from ASCII files where values are organized in matrices. +* File distribution, acting as the manual distribution, but where the points are defined in a separate ASCII file, whose name is given as a parameter. This method is appropriate for large numbers of points and allows to describe any distribution in a totally generic way. Now, GATE supports reading 2D distributions from ASCII files where values are organized in matrices. A distribution is declared by specifying its name then by creating a new instance, with its type name:: From d0992885a49c59d815acb0b02185a958211bcd9c Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Tue, 2 Jul 2024 16:05:51 +0200 Subject: [PATCH 03/22] Update digitizer_and_detector_modeling.rst --- docs/digitizer_and_detector_modeling.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index e3dfdcf22..04ddaabee 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -529,7 +529,7 @@ BEWARE: This relocation procedure is validated only for the first group level of This approach is particularly essential for monolithic crystal detectors, where factors like edge effects and interaction positions significantly influence spatial resolution. Here is an example of how to configure this in a macro file: -**Example**:: +**Example for 2D distribution**:: /gate/distributions/name my_distrib2D @@ -537,6 +537,8 @@ Here is an example of how to configure this in a macro file: /gate/distributions/my_distrib2D/setFileName Lut(X,Y).txt /gate/distributions/my_distrib2D/readMatrix2d +**Example for 2D distribution**:: + /gate/distributions/name my_distrib1D /gate/distributions/insert File /gate/distributions/my_distrib1D/setFileName macros/LutY.txt From 40d371ca26a3725f38b7764f1889744185c09fcb Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Wed, 3 Jul 2024 12:12:56 +0200 Subject: [PATCH 04/22] Update digitizer_and_detector_modeling.rst --- docs/digitizer_and_detector_modeling.rst | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 04ddaabee..73f611bfd 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -526,7 +526,7 @@ BEWARE: This relocation procedure is validated only for the first group level of **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 influence spatial resolution. +This approach is particularly essential for monolithic crystal detectors, where factors like edge effects and interaction positions significantly may influence spatial resolution. Here is an example of how to configure this in a macro file: **Example for 2D distribution**:: @@ -536,21 +536,20 @@ Here is an example of how to configure this in a macro file: /gate/distributions/insert File /gate/distributions/my_distrib2D/setFileName Lut(X,Y).txt /gate/distributions/my_distrib2D/readMatrix2d - -**Example for 2D distribution**:: + /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution + /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmXdistrib2D 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/read + /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution + /gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmYdistrib my_distrib1D + -Configure the spatial resolution for the digitizer using these distributions:: -/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/insert spatialResolution -/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmXdistrib2D my_distrib2D -/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/fwhmYdistrib my_distrib1D -/gate/digitizerMgr/crystalUnit/SinglesDigitizer/Singles/spatialResolution/confineInsideOfSmallestElement true These commands allow for more precise control over the spatial resolution by using predefined distributions for the X and Y axes. From 3dee074fa2873b5828e168b99c9eb2a502b1359e Mon Sep 17 00:00:00 2001 From: Thomas BAUDIER Date: Tue, 23 Jul 2024 09:08:55 +0200 Subject: [PATCH 05/22] Try debug Gate compilation on macos changing cache --- .github/workflows/main.yml | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index eede25bd1..47f126617 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -18,9 +18,12 @@ jobs: matrix: os: [ubuntu-latest, macos-latest] options: [none, rtk, torch, optic] + exclude: + - os: macos-latest + options: torch env: - ROOT_VERSION: 'v6-26-08' + ROOT_VERSION: 'v6-32-02' GEANT4_VERSION: 'v11.2.1' ITK_VERSION: 'v5.3.0' ROOT_DIR: $(HOME)/software/root @@ -38,8 +41,8 @@ jobs: uses: actions/cache@v3 with: path: ~/software - key: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-build4 - restore-keys: ${{ matrix.os }}-geant4-${{ env.GEANT4_VERSION }}-root-${{ env.ROOT_VERSION }}-build4 + 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 - name: Install dependencies run: | if [ "${{ matrix.os }}" == "ubuntu-latest" ]; then @@ -129,8 +132,8 @@ jobs: wget https://download.pytorch.org/libtorch/cpu/libtorch-cxx11-abi-shared-with-deps-1.10.1%2Bcpu.zip unzip libtorch-cxx11-abi-shared-with-deps-1.10.1+cpu.zip elif [ "${{ matrix.os }}" == 'macos-latest' ]; then - wget https://download.pytorch.org/libtorch/cpu/libtorch-macos-1.10.1.zip - unzip libtorch-macos-1.10.1.zip + wget https://download.pytorch.org/libtorch/cpu/libtorch-macos-arm64-2.2.0.zip + unzip libtorch-macos-arm64-2.2.0.zip fi fi if [ "${{ matrix.options }}" == 'rtk' ]; then From c61303a4a1a2c3856019164e77e224b674e819db Mon Sep 17 00:00:00 2001 From: kochebina Date: Fri, 31 May 2024 13:39:01 +0200 Subject: [PATCH 06/22] Common Output in GAteToTree --- source/digits_hits/include/GateToTree.hh | 5 ++ .../include/GateToTreeMessenger.hh | 4 ++ source/digits_hits/src/GateToTree.cc | 47 ++++++++++++++++--- source/digits_hits/src/GateToTreeMessenger.cc | 11 +++++ 4 files changed, 60 insertions(+), 7 deletions(-) diff --git a/source/digits_hits/include/GateToTree.hh b/source/digits_hits/include/GateToTree.hh index f4937d048..46de07d5c 100644 --- a/source/digits_hits/include/GateToTree.hh +++ b/source/digits_hits/include/GateToTree.hh @@ -65,6 +65,10 @@ public: static void SetOutputIDName(G4int id_system, const char * anOutputIDName, size_t depth); G4bool getHitsEnabled() const; void setHitsEnabled(G4bool mHitsEnabled); + + G4bool getHitsCommonOutputEnabled() const; + void setHitsCommonOutputEnabled(G4bool mHitsCommonOutputEnabled); + void addCollection(const std::string &str); //called by messenger //OK GND 2022 void setCCenabled(G4bool mCCenabled){m_cc_enabled=mCCenabled;}; @@ -135,6 +139,7 @@ private: std::vector m_listOfSinglesCollection; std::vector m_listOfCoincidencesCollection; G4bool m_hits_enabled; + G4bool m_hitsCommonOutput_enabled; G4String m_uselessFileName; //only for GiveNameOfFile which return a reference.. G4bool m_opticalData_enabled = false; diff --git a/source/digits_hits/include/GateToTreeMessenger.hh b/source/digits_hits/include/GateToTreeMessenger.hh index 27ece011e..2737d38d8 100644 --- a/source/digits_hits/include/GateToTreeMessenger.hh +++ b/source/digits_hits/include/GateToTreeMessenger.hh @@ -41,6 +41,10 @@ private: G4UIcmdWithoutParameter *m_enableHitsOutput; G4UIcmdWithoutParameter *m_disableHitsOutput; + G4UIcmdWithoutParameter *m_enableHitsCommonOutput; + G4UIcmdWithoutParameter *m_disableHitsCommonOutput; + + G4UIcmdWithoutParameter *m_enableOpticalDataOutput; G4UIcmdWithoutParameter *m_disableOpticalDataOutput; diff --git a/source/digits_hits/src/GateToTree.cc b/source/digits_hits/src/GateToTree.cc index 60e12700a..48e5f03e2 100644 --- a/source/digits_hits/src/GateToTree.cc +++ b/source/digits_hits/src/GateToTree.cc @@ -193,13 +193,21 @@ GateToTree::GateToTree(const G4String &name, GateOutputMgr *outputMgr, DigiMode m_messenger = new GateToTreeMessenger(this); m_hits_enabled = false; + m_hitsCommonOutput_enabled = false; + + + + } void GateToTree::RecordBeginOfAcquisition() { + if (m_hits_enabled && m_hitsCommonOutput_enabled) + GateError("Commands /hits/ and /hitsCommonOutput/ cannot be enabled at the same time"); + if (!this->IsEnabled()) return; - /* if (m_hits_enabled) { + if (m_hitsCommonOutput_enabled) { for (auto &&fileName: m_listOfFileName) { auto extension = getExtension(fileName); auto name = removeExtension(fileName); @@ -207,7 +215,8 @@ void GateToTree::RecordBeginOfAcquisition() { G4String hits_filename = name + ".hits." + extension; m_manager_hits.add_file(hits_filename, extension); } - */ + } + //OK GND 2022 for (auto &&m: m_mmanager_hits) { @@ -691,9 +700,15 @@ void GateToTree::RecordEndOfAcquisition() { //m_manager_optical.close(); //OK GND 2022 - + if(m_hitsCommonOutput_enabled) + { + m_manager_hits.close(); + } + else + { for (auto &&m: m_mmanager_hits) m.second.close(); + } for (auto &&m: m_mmanager_optical) m.second.close(); @@ -719,9 +734,10 @@ void GateToTree::RecordBeginOfEvent(const G4Event *event) { } void GateToTree::RecordEndOfEvent(const G4Event *event) { + auto fDM = G4DigiManager::GetDMpointer(); //OK GND 2022 - auto fDM = G4DigiManager::GetDMpointer(); + if (!m_hits_to_collectionID.size()) { for (auto &&m: m_mmanager_hits) { @@ -829,12 +845,19 @@ void GateToTree::RecordEndOfEvent(const G4Event *event) { m_energyIniT=hit->GetEnergyIniTrack(); } + if(m_hitsCommonOutput_enabled) + { + m_manager_hits.fill(); + } + else + { m.second.fill(); - //m_manager_hits.fill(); + } + } - } + } //auto fDM = G4DigiManager::GetDMpointer(); - + //} if (!m_singles_to_collectionID.size()) { for (auto &&m: m_mmanager_singles) { // auto collectionID = fDM->GetDigiCollectionID(m.first); @@ -1083,6 +1106,15 @@ void GateToTree::setHitsEnabled(G4bool mHitsEnabled) { } +G4bool GateToTree::getHitsCommonOutputEnabled() const { + return m_hitsCommonOutput_enabled; +} + +void GateToTree::setHitsCommonOutputEnabled(G4bool mHitsCommonOutputEnabled) { + m_hitsCommonOutput_enabled = mHitsCommonOutputEnabled; +} + + //OK GND 2022 void GateToTree::addHitsCollection(const std::string &str) { GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); @@ -1104,6 +1136,7 @@ void GateToTree::addHitsCollection(const std::string &str) { n_fileName = name + ".hits." + extension; else n_fileName = name + ".hits_" + str + "." + extension; + if(!m_hitsCommonOutput_enabled) m.add_file(n_fileName, extension); } diff --git a/source/digits_hits/src/GateToTreeMessenger.cc b/source/digits_hits/src/GateToTreeMessenger.cc index 6a4a2ea6c..2e709fb69 100644 --- a/source/digits_hits/src/GateToTreeMessenger.cc +++ b/source/digits_hits/src/GateToTreeMessenger.cc @@ -37,6 +37,10 @@ GateToTreeMessenger::GateToTreeMessenger(GateToTree *m) : m_enableHitsOutput = new G4UIcmdWithoutParameter("/gate/output/tree/hits/enable", this); m_disableHitsOutput = new G4UIcmdWithoutParameter("/gate/output/tree/hits/disable", this); + m_enableHitsCommonOutput = new G4UIcmdWithoutParameter("/gate/output/tree/hitsCommonOutput/enable", this); + m_disableHitsCommonOutput = new G4UIcmdWithoutParameter("/gate/output/tree/hitsCommonOutput/disable", this); + + m_enableOpticalDataOutput = new G4UIcmdWithoutParameter("/gate/output/tree/optical/enable", this); m_disableOpticalDataOutput = new G4UIcmdWithoutParameter("/gate/output/tree/optical/disable", this); @@ -103,6 +107,8 @@ GateToTreeMessenger::~GateToTreeMessenger() delete m_addOpticalCollectionCmd; delete m_enableHitsOutput; delete m_disableHitsOutput; + delete m_enableHitsCommonOutput; + delete m_disableHitsCommonOutput; } @@ -119,6 +125,11 @@ void GateToTreeMessenger::SetNewValue(G4UIcommand *icommand, G4String string) if(icommand == m_disableHitsOutput) m_gateToTree->setHitsEnabled(false); + if(icommand == m_enableHitsCommonOutput) + m_gateToTree->setHitsCommonOutputEnabled(true); + if(icommand == m_disableHitsCommonOutput) + m_gateToTree->setHitsCommonOutputEnabled(false); + if(icommand == m_enableOpticalDataOutput) m_gateToTree->setOpticalDataEnabled(true); if(icommand == m_disableOpticalDataOutput) From 89286f8634aa473fd2db3656026aa3496eaffaba Mon Sep 17 00:00:00 2001 From: kochebina Date: Fri, 7 Jun 2024 17:00:35 +0200 Subject: [PATCH 07/22] Branches Bug Correction --- source/digits_hits/src/GateToTree.cc | 154 ++++++++++++++++++++++++++- 1 file changed, 152 insertions(+), 2 deletions(-) diff --git a/source/digits_hits/src/GateToTree.cc b/source/digits_hits/src/GateToTree.cc index 48e5f03e2..f85f06eaa 100644 --- a/source/digits_hits/src/GateToTree.cc +++ b/source/digits_hits/src/GateToTree.cc @@ -215,6 +215,155 @@ void GateToTree::RecordBeginOfAcquisition() { G4String hits_filename = name + ".hits." + extension; m_manager_hits.add_file(hits_filename, extension); } + if (m_hitsParams_to_write.at("PDGEncoding").toSave()) + m_manager_hits.write_variable("PDGEncoding", &m_PDGEncoding); + + if (m_hitsParams_to_write.at("trackID").toSave()) + m_manager_hits.write_variable("trackID", &m_trackID); + + if (m_hitsParams_to_write.at("parentID").toSave()) + m_manager_hits.write_variable("parentID", &m_parentID); + + if (m_hitsParams_to_write.at("trackLocalTime").toSave()) + m_manager_hits.write_variable("trackLocalTime", &m_trackLocalTime); + + if (m_hitsParams_to_write.at("time").toSave()) + m_manager_hits.write_variable("time", &m_time[0]); + + if (m_hitsParams_to_write.at("runID").toSave()) + m_manager_hits.write_variable("runID", &m_runID); + + if (m_hitsParams_to_write.at("eventID").toSave()) + m_manager_hits.write_variable("eventID", &m_eventID[0]); + + if (m_hitsParams_to_write.at("sourceID").toSave()) + m_manager_hits.write_variable("sourceID", &m_sourceID[0]); + + if (m_hitsParams_to_write.at("primaryID").toSave()) + m_manager_hits.write_variable("primaryID", &m_primaryID); + + if (m_hitsParams_to_write.at("posX").toSave()) + m_manager_hits.write_variable("posX", &m_posX[0]); + + if (m_hitsParams_to_write.at("posY").toSave()) + m_manager_hits.write_variable("posY", &m_posY[0]); + + if (m_hitsParams_to_write.at("posZ").toSave()) + m_manager_hits.write_variable("posZ", &m_posZ[0]); + + if (m_hitsParams_to_write.at("localPosX").toSave()) + m_manager_hits.write_variable("localPosX", &m_localPosX); + + if (m_hitsParams_to_write.at("localPosY").toSave()) + m_manager_hits.write_variable("localPosY", &m_localPosY); + + if (m_hitsParams_to_write.at("localPosZ").toSave()) + m_manager_hits.write_variable("localPosZ", &m_localPosZ); + + if (m_hitsParams_to_write.at("momDirX").toSave()) + m_manager_hits.write_variable("momDirX", &m_momDirX); + + if (m_hitsParams_to_write.at("momDirY").toSave()) + m_manager_hits.write_variable("momDirY", &m_momDirY); + + if (m_hitsParams_to_write.at("momDirZ").toSave()) + m_manager_hits.write_variable("momDirZ", &m_momDirZ); + + if (m_hitsParams_to_write.at("edep").toSave()) + m_manager_hits.write_variable("edep", &m_edep[0]); + + if (m_hitsParams_to_write.at("stepLength").toSave()) + m_manager_hits.write_variable("stepLength", &m_stepLength); + + if (m_hitsParams_to_write.at("trackLength").toSave()) + m_manager_hits.write_variable("trackLength", &m_trackLength); + + if (m_hitsParams_to_write.at("rotationAngle").toSave()) + m_manager_hits.write_variable("rotationAngle", &m_rotationAngle); + + if (m_hitsParams_to_write.at("axialPos").toSave()) + m_manager_hits.write_variable("axialPos", &m_axialPos); + + if (m_hitsParams_to_write.at("processName").toSave()) + m_manager_hits.write_variable("processName", &m_processName, MAX_NB_CHARACTER); + + if (m_hitsParams_to_write.at("comptVolName").toSave()) + m_manager_hits.write_variable("comptVolName", &m_comptonVolumeName[0], MAX_NB_CHARACTER); + + if (m_hitsParams_to_write.at("RayleighVolName").toSave()) + m_manager_hits.write_variable("RayleighVolName", &m_RayleighVolumeName[0], MAX_NB_CHARACTER); + + if (m_hitsParams_to_write.at("volumeIDs").toSave()) { + for (auto i = 0; i < VOLUMEID_SIZE; ++i) { + std::stringstream ss; + ss << "volumeID[" << i << "]"; + m_manager_hits.write_variable(ss.str(), &m_volumeID[i]); + } + } + + if (m_hitsParams_to_write.at("sourcePosX").toSave()) + m_manager_hits.write_variable("sourcePosX", &m_sourcePosX[0]); + if (m_hitsParams_to_write.at("sourcePosY").toSave()) + m_manager_hits.write_variable("sourcePosY", &m_sourcePosY[0]); + if (m_hitsParams_to_write.at("sourcePosZ").toSave()) + m_manager_hits.write_variable("sourcePosZ", &m_sourcePosZ[0]); + + if (m_hitsParams_to_write.at("nPhantomCompton").toSave()) + m_manager_hits.write_variable("nPhantomCompton", &m_nPhantomCompton[0]); + + if (m_hitsParams_to_write.at("nCrystalCompton").toSave()) + m_manager_hits.write_variable("nCrystalCompton", &m_nCrystalCompton[0]); + + if (m_hitsParams_to_write.at("nPhantomRayleigh").toSave()) + m_manager_hits.write_variable("nPhantomRayleigh", &m_nPhantomRayleigh[0]); + + if (m_hitsParams_to_write.at("nCrystalRayleigh").toSave()) + m_manager_hits.write_variable("nCrystalRayleigh", &m_nCrystalRayleigh[0]); + + + if (m_hitsParams_to_write.at("componentsIDs").toSave()) { + + if (GateSystemListManager::GetInstance()->size() == 1) { + int k = 0; + for (auto depth = 0; depth < m_max_depth_system[k]; ++depth) { + std::stringstream ss; + ss << m_outputIDName[k][depth]; + m_manager_hits.write_variable(ss.str(), &m_outputID[0][k][depth]); + } + } else { + for (unsigned int k = 0; k < GateSystemListManager::GetInstance()->size(); ++k) { + auto system = GateSystemListManager::GetInstance()->GetSystem(k); + + for (auto depth = 0; depth < m_max_depth_system[k]; ++depth) { + if (!m_outputIDHasName[k][depth]) + continue; + std::stringstream ss; + ss << system->GetOwnName() << "/" << m_outputIDName[k][depth]; + m_manager_hits.write_variable(ss.str(), &m_outputID[0][k][depth]); + } + } + } + } + + if (m_hitsParams_to_write.at("photonID").toSave()) + m_manager_hits.write_variable("photonID", &m_photonID); + + if (m_hitsParams_to_write.at("systemID").toSave() && GateSystemListManager::GetInstance()->size() > 1) + m_manager_hits.write_variable("systemID", &m_systemID); + + if (m_hitsParams_to_write.at("sourceType").toSave()) + m_manager_hits.write_variable("sourceType", &m_sourceType); + + if (m_hitsParams_to_write.at("decayType").toSave()) + m_manager_hits.write_variable("decayType", &m_decayType); + + if (m_hitsParams_to_write.at("gammaType").toSave()) + m_manager_hits.write_variable("gammaType", &m_gammaType); + + + m_manager_hits.write_header(); + + } //OK GND 2022 @@ -1117,6 +1266,7 @@ void GateToTree::setHitsCommonOutputEnabled(G4bool mHitsCommonOutputEnabled) { //OK GND 2022 void GateToTree::addHitsCollection(const std::string &str) { + GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); G4String possibleValues = ""; @@ -1136,8 +1286,8 @@ void GateToTree::addHitsCollection(const std::string &str) { n_fileName = name + ".hits." + extension; else n_fileName = name + ".hits_" + str + "." + extension; - if(!m_hitsCommonOutput_enabled) - m.add_file(n_fileName, extension); + if(m_hits_enabled) + m.add_file(n_fileName, extension); } m_mmanager_hits.emplace(str, std::move(m)); From 0dec706eee8a26e23053dc493a7b5e45fba5a189 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Thu, 27 Jun 2024 09:30:30 +0200 Subject: [PATCH 08/22] Add option for spatial resolution 2D --- .../include/GateDistributionFile.hh | 8 +- .../include/GateDistributionFileMessenger.hh | 1 + .../include/GateDistributionMessenger.hh | 5 +- .../include/GateSpatialResolution.hh | 29 +++- .../include/GateSpatialResolutionMessenger.hh | 5 +- .../digits_hits/include/GateVDistribution.hh | 2 + .../include/GateVDistributionArray.hh | 24 ++- .../digits_hits/src/GateDistributionFile.cc | 65 ++++++- .../src/GateDistributionFileMessenger.cc | 9 +- .../src/GateDistributionMessenger.cc | 7 +- .../digits_hits/src/GateSpatialResolution.cc | 106 ++++++++++-- .../src/GateSpatialResolutionMessenger.cc | 48 +++++- source/digits_hits/src/GateVDistribution.cc | 8 + .../digits_hits/src/GateVDistributionArray.cc | 161 +++++++++++++++--- 14 files changed, 404 insertions(+), 74 deletions(-) diff --git a/source/digits_hits/include/GateDistributionFile.hh b/source/digits_hits/include/GateDistributionFile.hh index e8d4e73d2..cdf88cac0 100644 --- a/source/digits_hits/include/GateDistributionFile.hh +++ b/source/digits_hits/include/GateDistributionFile.hh @@ -32,7 +32,8 @@ class GateDistributionFile : public GateVDistributionArray inline G4int GetColumnX() const {return m_column_for_X;} inline G4int GetColumnY() const {return m_column_for_Y;} - void Read(); + void Read(); + void ReadMatrix2d(); virtual void DescribeMyself(size_t indent); private: @@ -41,6 +42,11 @@ class GateDistributionFile : public GateVDistributionArray G4int m_column_for_X; G4int m_column_for_Y; GateDistributionFileMessenger* m_messenger; + std::map, double> stddevMap; + + std::vector xValues; + std::vector yValues; + }; diff --git a/source/digits_hits/include/GateDistributionFileMessenger.hh b/source/digits_hits/include/GateDistributionFileMessenger.hh index 86ccc9426..281b88e3e 100644 --- a/source/digits_hits/include/GateDistributionFileMessenger.hh +++ b/source/digits_hits/include/GateDistributionFileMessenger.hh @@ -33,6 +33,7 @@ class GateDistributionFileMessenger: public GateDistributionArrayMessenger G4UIcmdWithAnInteger *setColXCmd; G4UIcmdWithAnInteger *setColYCmd; G4UIcmdWithoutParameter *readCmd; + G4UIcmdWithoutParameter *read2DCmd; G4UIcmdWithoutParameter *autoXCmd; }; diff --git a/source/digits_hits/include/GateDistributionMessenger.hh b/source/digits_hits/include/GateDistributionMessenger.hh index b33f3af72..747c6852a 100644 --- a/source/digits_hits/include/GateDistributionMessenger.hh +++ b/source/digits_hits/include/GateDistributionMessenger.hh @@ -30,9 +30,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger void SetNewValue(G4UIcommand* aCommand, G4String aString); void SetUnitX(const G4String& unitX); void SetUnitY(const G4String& unitY); + void SetUnitSigma(const G4String& unitSigma); inline G4String UnitCategoryX() const {return m_unitX.empty()?"":G4UIcommand::CategoryOf(m_unitX);} inline G4String UnitCategoryY() const {return m_unitY.empty()?"":G4UIcommand::CategoryOf(m_unitY);} - + inline G4String UnitCategorySigma() const {return m_unitSigma.empty()?"":G4UIcommand::CategoryOf(m_unitSigma);} private: G4String withUnity(G4double value,G4String category) const; G4UIcmdWithoutParameter *getMinX_Cmd ; @@ -41,8 +42,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger G4UIcmdWithoutParameter *getMaxY_Cmd ; G4UIcmdWithoutParameter *getRandom_Cmd ; G4UIcmdWithADoubleAndUnit *getValueCmd ; + G4UIcmdWithADoubleAndUnit *getSigmaCmd ; G4String m_unitX; G4String m_unitY; + G4String m_unitSigma; }; #endif diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh index 171bf855e..502bff961 100755 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/source/digits_hits/include/GateSpatialResolution.hh @@ -31,7 +31,10 @@ See LICENSE.md for further details #include "G4TouchableHistoryHandle.hh" #include "GateSinglesDigitizer.hh" +class GateVDistribution; + class GateSpatialResolution : public GateVDigitizerModule + { public: @@ -40,9 +43,15 @@ public: void Digitize() override; + + //! These functions return the resolution in use. G4double GetFWHM() { return m_fwhm; } - G4double GetFWHMx() { return m_fwhmX; } + GateVDistribution* GetFWHMxdistrib() { return m_fwhmXdistrib; } + GateVDistribution* GetFWHMydistrib() { return m_fwhmYdistrib; } + GateVDistribution* GetFWHMxdistrib2D() { return m_fwhmXdistrib2D; } + GateVDistribution* GetFWHMydistrib2D() { return m_fwhmYdistrib2D; } + G4double GetFWHMx() { return m_fwhmX; } G4double GetFWHMy() { return m_fwhmY; } G4double GetFWHMz() { return m_fwhmZ; } @@ -51,11 +60,14 @@ 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 SetFWHMxdistrib2D(GateVDistribution* dist) { m_fwhmXdistrib2D= dist; } + void SetFWHMydistrib2D(GateVDistribution* dist) { m_fwhmYdistrib2D = dist; } void SetFWHMx(G4double val) { m_fwhmX = val; } void SetFWHMy(G4double val) { m_fwhmY = val; } void SetFWHMz(G4double val) { m_fwhmZ = val; } - inline void ConfineInsideOfSmallestElement(const G4bool& value) { m_IsConfined = value; }; inline G4bool IsConfinedInsideOfSmallestElement() const { return m_IsConfined; } @@ -64,7 +76,6 @@ public: void UpdateVolumeID(); - //! Implementation of the pure virtual method declared by the base class GateClockDependent //! print-out the attributes specific of the blurring void DescribeMyself(size_t ); @@ -72,9 +83,16 @@ public: protected: G4double m_fwhm; + G4double m_fwhmX; + G4double m_fwhmY; G4double m_fwhmZ; + + GateVDistribution* m_fwhmXdistrib; + GateVDistribution* m_fwhmYdistrib; + GateVDistribution* m_fwhmXdistrib2D; + GateVDistribution* m_fwhmYdistrib2D; G4bool m_IsConfined; G4Navigator* m_Navigator; G4TouchableHistoryHandle m_Touchable; @@ -83,10 +101,9 @@ protected: private: - G4int m_systemDepth; - - GateDigi* m_outputDigi; + G4int m_systemDepth; + GateDigi* m_outputDigi;; GateSpatialResolutionMessenger *m_Messenger; GateDigiCollection* m_OutputDigiCollection; diff --git a/source/digits_hits/include/GateSpatialResolutionMessenger.hh b/source/digits_hits/include/GateSpatialResolutionMessenger.hh index 7b0d03c81..db588a6e3 100755 --- a/source/digits_hits/include/GateSpatialResolutionMessenger.hh +++ b/source/digits_hits/include/GateSpatialResolutionMessenger.hh @@ -27,7 +27,6 @@ See LICENSE.md for further details #include "GateClockDependentMessenger.hh" class GateSpatialResolution; class G4UIcmdWithAString; - class GateSpatialResolutionMessenger : public GateClockDependentMessenger { public: @@ -43,6 +42,10 @@ private: G4UIcmdWithADouble* spresolutionCmd; G4UIcmdWithADouble* spresolutionXCmd; + G4UIcmdWithAString *spresolutionXdistribCmd; + G4UIcmdWithAString *spresolutionYdistribCmd; + G4UIcmdWithAString *spresolutionXdistrib2DCmd; + G4UIcmdWithAString *spresolutionYdistrib2DCmd; G4UIcmdWithADouble* spresolutionYCmd; G4UIcmdWithADouble* spresolutionZCmd; G4UIcmdWithABool* confineCmd; diff --git a/source/digits_hits/include/GateVDistribution.hh b/source/digits_hits/include/GateVDistribution.hh index 7a428af02..ff524ccf2 100644 --- a/source/digits_hits/include/GateVDistribution.hh +++ b/source/digits_hits/include/GateVDistribution.hh @@ -41,6 +41,8 @@ class GateVDistribution : public GateNamedObject virtual G4double MaxX() const=0; virtual G4double MaxY() const=0; virtual G4double Value(G4double x) const=0; + virtual G4double Value2D(G4double x, G4double y) const; + // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const=0; diff --git a/source/digits_hits/include/GateVDistributionArray.hh b/source/digits_hits/include/GateVDistributionArray.hh index 226c23239..3353efb80 100644 --- a/source/digits_hits/include/GateVDistributionArray.hh +++ b/source/digits_hits/include/GateVDistributionArray.hh @@ -27,35 +27,43 @@ class GateVDistributionArray : public GateVDistribution inline void SetFactorX(G4double factor) {m_factorX=factor;} inline void SetFactorY(G4double factor) {m_factorY=factor;} - G4double Integral() const {return m_arrayRepartition.back();} + G4double Integral() const {return m_arrayRepartition.back();} virtual G4double MinX() const; virtual G4double MinY() const; virtual G4double MaxX() const; virtual G4double MaxY() const; virtual G4double Value(G4double x) const; + virtual G4double RepartitionValue(G4double x) const; + G4double BilinearInterpolation(G4double x, G4double y) const ; // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const; size_t GetSize() const {return m_arrayX.size();} - void Clear(); void SetAutoStart(G4int start) {m_autoStart=start;} + virtual G4double Value2D(G4double x, G4double y) const; protected: - void InsertPoint(G4double x,G4double y); - void InsertPoint(G4double y); + void InsertPoint(G4double x,G4double y, G4double sigma ); + void InsertPoint(G4double x, G4double y ); + void InsertPoint(G4double y ); //! private function G4int FindIdxBefore(G4double x ,const std::vector& array) const; - void FillRepartition(); + + + void FillRepartition() ; + std::vector& GetArrayX() {return m_arrayX;} std::vector& GetArrayY() {return m_arrayY;} + std::vector& GetArrayRepartition() {return m_arrayRepartition;} private: //! private members std::vector m_arrayX; std::vector m_arrayY; + G4double m_minX; G4double m_minY; G4double m_maxX; @@ -63,7 +71,13 @@ class GateVDistributionArray : public GateVDistribution std::vector m_arrayRepartition; //! repartition function calculus G4double m_factorX; G4double m_factorY; + G4int m_autoStart; + + std::map, double> stddevMap; + std::vector xValues; // Stocker les valeurs x + std::vector yValues; // Stocker les valeurs y + }; diff --git a/source/digits_hits/src/GateDistributionFile.cc b/source/digits_hits/src/GateDistributionFile.cc index 611c74092..973380e79 100644 --- a/source/digits_hits/src/GateDistributionFile.cc +++ b/source/digits_hits/src/GateDistributionFile.cc @@ -20,6 +20,7 @@ GateDistributionFile::GateDistributionFile(const G4String& itsName) , m_FileName() , m_column_for_X(0) , m_column_for_Y(1) + { m_messenger = new GateDistributionFileMessenger(this,itsName); } @@ -30,10 +31,10 @@ GateDistributionFile::~GateDistributionFile() //___________________________________________________________________ void GateDistributionFile::DescribeMyself(size_t indent) { - G4cout << GateTools::Indent(indent) - <<"File : " << m_FileName - <<'{' << m_column_for_X<<':'< xValues; + bool isFirstLine = true; + + while (std::getline(f, line)) { + std::stringstream iss(line); + + if (isFirstLine) { + // Process x values from the first line + G4double x; + + while (iss >> x) { + xValues.push_back(x); + if (iss.peek() == ',') iss.ignore(); + } + isFirstLine = false; + } else { + // Process y and stddev values for subsequent lines + G4double y; + iss >> y; + if (iss.peek() == ',') iss.ignore(); + // Read stddev values for each x value + for (size_t i = 0; i < xValues.size(); ++i) { + G4double stddev; + if (iss >> stddev) { + // Insert the (x, y) -> stddev pair into the map using InsertPoint + InsertPoint(xValues[i], y, stddev); + if (iss.peek() == ',') iss.ignore(); + } + } + } + } + + f.close(); + //G4cout << "Content of stddevMap:\n"; + //for (const auto& entry : stddevMap) { + // std::pair coordinates = entry.first; + // G4double stddev = entry.second; + //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl; + //} + + FillRepartition(); + } + diff --git a/source/digits_hits/src/GateDistributionFileMessenger.cc b/source/digits_hits/src/GateDistributionFileMessenger.cc index aae879075..107683d7b 100644 --- a/source/digits_hits/src/GateDistributionFileMessenger.cc +++ b/source/digits_hits/src/GateDistributionFileMessenger.cc @@ -45,7 +45,10 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil guidance = "Do read the file"; readCmd = new G4UIcmdWithoutParameter(cmdName,this); readCmd->SetGuidance(guidance); - + cmdName = GetDirectoryName()+"readMatrix2d"; + guidance = "Do read the file"; + read2DCmd = new G4UIcmdWithoutParameter(cmdName,this); + read2DCmd->SetGuidance(guidance); cmdName = GetDirectoryName()+"autoX"; guidance = "Set automatic mode for X"; autoXCmd = new G4UIcmdWithoutParameter(cmdName,this); @@ -58,6 +61,7 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil GateDistributionFileMessenger::~GateDistributionFileMessenger() { delete readCmd; + delete read2DCmd; delete autoXCmd; delete setColYCmd; delete setColXCmd; @@ -83,6 +87,9 @@ void GateDistributionFileMessenger::SetNewValue(G4UIcommand* command,G4String ne } else if( command==readCmd ) { GetDistributionFile()->Read(); } + else if( command==read2DCmd ){ + GetDistributionFile()->ReadMatrix2d(); + } else GateDistributionArrayMessenger::SetNewValue(command,newValue); } diff --git a/source/digits_hits/src/GateDistributionMessenger.cc b/source/digits_hits/src/GateDistributionMessenger.cc index 6272a88c3..0cd6bd444 100644 --- a/source/digits_hits/src/GateDistributionMessenger.cc +++ b/source/digits_hits/src/GateDistributionMessenger.cc @@ -54,6 +54,7 @@ GateDistributionMessenger::GateDistributionMessenger(GateVDistribution* itsDistr getValueCmd->SetGuidance(guidance); getValueCmd->SetParameterName("x",false); + } @@ -67,6 +68,7 @@ GateDistributionMessenger::~GateDistributionMessenger() delete getMaxY_Cmd; delete getRandom_Cmd; delete getValueCmd; + } G4String GateDistributionMessenger::withUnity(G4double value,G4String category) const { @@ -82,9 +84,8 @@ void GateDistributionMessenger::SetNewValue(G4UIcommand* command,G4String newVal { if ( command==getValueCmd ){ G4double x = getValueCmd->GetNewDoubleValue(newValue); - G4double y = GetDistribution()->Value(x); - G4cout<GetObjectName()<<'('<Value(x);} +else if( command==getMinX_Cmd ) { G4double x = GetDistribution()->MinX(); G4cout<GetObjectName()<<" MinX "<GetSD()->GetName()+"/SinglesDigitizer/"+digitizer->m_digitizerName+"/"+name,digitizer,digitizer->GetSD()), m_fwhm(0), m_fwhmX(0), + m_fwhmXdistrib(0), + m_fwhmYdistrib(0), + m_fwhmXdistrib2D(0), + m_fwhmYdistrib2D(0), m_fwhmY(0), m_fwhmZ(0), m_IsConfined(true), @@ -72,7 +81,7 @@ GateSpatialResolution::~GateSpatialResolution() void GateSpatialResolution::Digitize() { - if( m_fwhm!=0 && (m_fwhmX!=0 || m_fwhmY!=0 || m_fwhmZ!=0 ) ) + if( m_fwhm!=0 && (m_fwhmX!=0 || m_fwhmY!=0 || m_fwhmZ!=0 )) { GateError("***ERROR*** Spatial Resolution is ambiguous: you can set unique /fwhm for 3 axis OR set /fhwmX, /fhwmY, /fhwmZ!"); } @@ -81,17 +90,22 @@ void GateSpatialResolution::Digitize() G4double fwhmY; G4double fwhmZ; + if (m_fwhmX==0 && m_fwhmY==0 && m_fwhmZ==0 ) { fwhmX=m_fwhm; fwhmY=m_fwhm; fwhmZ=m_fwhm; + } else { fwhmX=m_fwhmX; fwhmY=m_fwhmY; fwhmZ=m_fwhmZ; + + + } @@ -139,11 +153,69 @@ void GateSpatialResolution::Digitize() G4double Px = P.x(); G4double Py = P.y(); G4double Pz = P.z(); - G4double PxNew = G4RandGauss::shoot(Px,fwhmX/GateConstants::fwhm_to_sigma); - G4double PyNew = G4RandGauss::shoot(Py,fwhmY/GateConstants::fwhm_to_sigma); - G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); //TC - + G4double stddevX, stddevY, stddevZ; + + if (m_fwhmXdistrib) { + // If the FWHM distribution for X is defined + if (m_fwhmYdistrib) { + // If the FWHM distribution for Y is also defined + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else if (m_fwhmY) { + // If the constant FWHM for Y is defined + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } // If the 2D distribution for Y is defined + else if (m_fwhmYdistrib2D) { + // Use the constant value for X and the 2D distribution for Y + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + stddevY = m_fwhmYdistrib2D->Value2D(P.x() * mm, P.y() * mm); + } + // If only the distribution for X is defined + else { + // Use the distribution for X and the constant value for Y + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } + + } else if (m_fwhmXdistrib2D) { + // If the 2D FWHM distribution for X is defined + if (m_fwhmYdistrib) { + // If the FWHM distribution for Y is defined + stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else if (m_fwhmY) { + // If the constant FWHM for Y is defined + stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } else { + // If neither the FWHM distribution nor constant for Y is defined + stddevX = stddevY = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); + } + } else if (m_fwhmYdistrib) { + // If the FWHM distribution for Y is defined + if (m_fwhmXdistrib) { + // If the FWHM distribution for X is also defined + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else if (m_fwhmX) { + // If the constant FWHM for X is defined + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else { + // If the 2D distribution for X is defined + stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } + } else { + // If neither the FWHM distributions for X nor Y are defined + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } + G4double PxNew = G4RandGauss::shoot(Px,stddevX); + G4double PyNew = G4RandGauss::shoot(Py,stddevY); + G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); if (m_IsConfined) { //set the position on the border of the crystal @@ -151,18 +223,22 @@ void GateSpatialResolution::Digitize() inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kXAxis, limits, at, Xmin, Xmax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kYAxis, limits, at, Ymin, Ymax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kZAxis, limits, at, Zmin, Zmax); + if(PxNewXmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - // + + + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC - //TC outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); + //TC + //outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); m_OutputDigiCollection->insert(m_outputDigi); - + //G4cout<<"PxNew"<Xmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - // - m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC //Getting world Volume // Do not use from TransportationManager as it is not recommended - G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); G4VPhysicalVolume *WorldVolume = navigator->GetWorldVolume(); m_Navigator = new G4Navigator(); m_Navigator->SetWorldVolume(WorldVolume); @@ -206,7 +281,7 @@ void GateSpatialResolution::Digitize() } } - } + } else { if (nVerboseLevel>1) @@ -237,5 +312,4 @@ void GateSpatialResolution::DescribeMyself(size_t indent ) if(m_fwhm) G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhm << Gateendl; else - G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhmX <<" "<< m_fwhmY<< " "<SetGuidance("Set the resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmX"; - spresolutionXCmd = new G4UIcmdWithADouble(cmdName,this); - spresolutionXCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); - + spresolutionXCmd= new G4UIcmdWithADouble(cmdName,this); + spresolutionXCmd->SetGuidance("Set the resolution "); cmdName = GetDirectoryName() + "fwhmY"; spresolutionYCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionYCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); @@ -41,7 +43,18 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol cmdName = GetDirectoryName() + "fwhmZ"; spresolutionZCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionZCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); - + cmdName = GetDirectoryName() + "fwhmXdistrib"; + spresolutionXdistribCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionXdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmYdistrib"; + spresolutionYdistribCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionYdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmXdistrib2D"; + spresolutionXdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionXdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmYdistrib2D"; + spresolutionYdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionYdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; confineCmd = new G4UIcmdWithABool(cmdName,this); confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); @@ -53,11 +66,14 @@ GateSpatialResolutionMessenger::~GateSpatialResolutionMessenger() { delete spresolutionCmd; delete spresolutionXCmd; + delete spresolutionXdistribCmd; + delete spresolutionYdistribCmd; + delete spresolutionXdistrib2DCmd; + delete spresolutionYdistrib2DCmd; delete spresolutionYCmd; delete spresolutionZCmd; delete confineCmd; - } @@ -65,8 +81,24 @@ void GateSpatialResolutionMessenger::SetNewValue(G4UIcommand * aCommand,G4String { if ( aCommand==spresolutionCmd ) { m_SpatialResolution->SetFWHM(spresolutionCmd->GetNewDoubleValue(newValue)); } - else if ( aCommand==spresolutionXCmd ) - { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } + else if ( aCommand==spresolutionXdistribCmd ) + { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib)m_SpatialResolution->SetFWHMxdistrib(distrib); +} + else if ( aCommand==spresolutionYdistribCmd ) + { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib)m_SpatialResolution->SetFWHMydistrib(distrib); + } + 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==spresolutionXCmd ) + { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionYCmd ) { m_SpatialResolution->SetFWHMy(spresolutionYCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionZCmd ) diff --git a/source/digits_hits/src/GateVDistribution.cc b/source/digits_hits/src/GateVDistribution.cc index ac9155771..400bdc904 100644 --- a/source/digits_hits/src/GateVDistribution.cc +++ b/source/digits_hits/src/GateVDistribution.cc @@ -19,3 +19,11 @@ GateVDistribution::GateVDistribution(const G4String& itsName) GateVDistribution::~GateVDistribution() { } + +G4double GateVDistribution::Value2D(G4double x, G4double y)const { + + + G4cerr << "Warning: GateVDistribution::Value2D always returns 0 for Gaussian, flat (uniform), and exponential distributions." < #include #include "GateTools.hh" - +#include GateVDistributionArray::GateVDistributionArray(const G4String& itsName) : GateVDistribution(itsName) @@ -60,9 +62,10 @@ void GateVDistributionArray::Clear() m_maxX=-DBL_MAX; m_maxX=-DBL_MAX; } -//___________________________________________________________________ +//_________________________________________________________________ G4double GateVDistributionArray::Value(G4double x) const { + if ( m_arrayX.empty() ) return 0; if ( m_arrayX.size() == 1 ) return m_arrayY[0]/m_factorY; x*=m_factorX; @@ -74,16 +77,26 @@ G4double GateVDistributionArray::Value(G4double x) const x2 = m_arrayX[1]; y2 = m_arrayY[0]; } else if (idx==(G4int)m_arrayX.size()-1){ x1 = m_arrayX[m_arrayX.size()-2]; y1 = m_arrayY[m_arrayX.size()-2]; - x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; + x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; } else { x1 = m_arrayX[idx] ; y1 = m_arrayY[idx]; x2 = m_arrayX[idx+1]; y2 = m_arrayY[idx+1]; } +//G4cout << "m_arrayX:" << Gateendl; + for (size_t i = 0; i < m_arrayX.size(); ++i) { + G4cout << m_arrayX[i] << Gateendl; + } - // return the linear interpolation + //G4cout << "m_arrayY:" << Gateendl; + for (size_t i = 0; i < m_arrayY.size(); ++i) { + G4cout << m_arrayY[i] << Gateendl; + } + + //return the linear interpolation return (y1 + (x-x1)*(y2-y1)/(x2-x1))/m_factorY; } -//___________________________________________________________________ + +//__________________________________________________________________ G4double GateVDistributionArray::RepartitionValue(G4double x) const { if ( m_arrayX.empty() ) return 0; @@ -119,24 +132,26 @@ G4double GateVDistributionArray::ShootRandom() const return (x1 + (y-y1)*(x2-x1)/(y2-y1))/m_factorX; } //___________________________________________________________________ -void GateVDistributionArray::FillRepartition() -{ - m_arrayRepartition.clear(); - if (m_arrayX.size()<2) return; - m_arrayRepartition.resize(m_arrayX.size()); - for (G4int i=0;i<(G4int)m_arrayX.size();++i){m_arrayRepartition[i]=0;} - for (G4int i=1;i<(G4int)m_arrayX.size();++i){ - G4double x1 = m_arrayX[i-1]; - G4double x2 = m_arrayX[i]; - G4double y1 = m_arrayY[i-1]; - G4double y2 = m_arrayY[i]; - m_arrayRepartition[i] = m_arrayRepartition[i-1] + (y1+y2)*(x2-x1)/2; - } - for (G4int i=0;i<(G4int)m_arrayX.size();++i){ - m_arrayRepartition[i]/=Integral(); -// G4cout<<"Repartition["<& array) const @@ -156,8 +171,8 @@ void GateVDistributionArray::InsertPoint(G4double x,G4double y) m_arrayY.push_back(y); } else if ( (i>=0) && m_arrayX[i] == x ) { //element replacement m_arrayY[i] = y; - G4cerr<<"[GateDistributionArray::InsertPoint] WARNING : replacing value for " - <m_maxY) m_maxY=y; } +void GateVDistributionArray::InsertPoint(G4double x, G4double y, G4double stddev) { + std::pair coordinates = std::make_pair(x, y); + stddevMap[coordinates] = stddev; + + + if (x < m_minX) m_minX = x; + if (x > m_maxX) m_maxX = x; + if (y < m_minY) m_minY = y; + if (y > m_maxY) m_maxY = y; + //G4cout << "Inserted: (x=" << x << ", y=" << y << ") -> stddev=" << stddev << std::endl; + + } + //___________________________________________________________________ void GateVDistributionArray::InsertPoint(G4double y) { if (GetSize()==0) InsertPoint(m_autoStart,y); else InsertPoint(m_maxX+1,y); } + + +G4double GateVDistributionArray::Value2D(double x, double y) const { + auto it = stddevMap.find(std::make_pair(x, y)); + if (it != stddevMap.end()) { + // If the point is found, return the standard deviation + G4double stddev= it->second; + // G4cout << "stddev for x=" << x << " and y=" << y << " is " << stddev<< G4endl; + return stddev; + } else { + //G4cout << "stddev not found for x=" << x << " and y=" << y << ". Performing bilinear interpolation." << G4endl; + G4double stddev = BilinearInterpolation(x, y); + return stddev; + } +} +G4double GateVDistributionArray::BilinearInterpolation(G4double x, G4double y) const { + auto lower_x_it = stddevMap.lower_bound(std::make_pair(x, -INFINITY)); + auto upper_x_it = stddevMap.upper_bound(std::make_pair(x, INFINITY)); + + if (lower_x_it == stddevMap.end() || upper_x_it == stddevMap.end()) { + G4cerr << "Error: Unable to find adjacent points in x direction." << G4endl; + return 0.0; + } + + // Adjust to make sure lower_x and upper_x are correct + G4double lower_x = (lower_x_it == stddevMap.begin()) ? lower_x_it->first.first : std::prev(lower_x_it)->first.first; + G4double upper_x = upper_x_it->first.first; + + if (lower_x == upper_x) { + G4cerr << "Error: Lower and upper x values are the same." << G4endl; + return 0.0; + } + + // Extract unique y values for these x coordinates + std::set y_values; + for (const auto& entry : stddevMap) { + if (entry.first.first == lower_x || entry.first.first == upper_x) { + y_values.insert(entry.first.second); + } + } + + if (y_values.empty()) { + G4cerr << "Error: No y values found for the given x values." << G4endl; + return 0.0; + } + + auto lower_y_it = y_values.lower_bound(y); + auto upper_y_it = y_values.upper_bound(y); + + if (lower_y_it == y_values.end() || upper_y_it == y_values.end()) { + G4cerr << "Error: Unable to find adjacent points in y direction." << G4endl; + return 0.0; + } + + G4double lower_y = (lower_y_it == y_values.begin()) ? *lower_y_it : *std::prev(lower_y_it); + G4double upper_y = *upper_y_it; + + if (lower_y == upper_y) { + G4cerr << "Error: Lower and upper y values are the same." << G4endl; + return 0.0; + } + + // Retrieve the sigma values for the four surrounding points + G4double stddev11 = stddevMap.at(std::make_pair(lower_x, lower_y)); + G4double stddev12 = stddevMap.at(std::make_pair(lower_x, upper_y)); + G4double stddev21 = stddevMap.at(std::make_pair(upper_x, lower_y)); + G4double stddev22 = stddevMap.at(std::make_pair(upper_x, upper_y)); + + //G4cout << "Interpolation points: (" << lower_x << ", " << lower_y << "), (" << upper_x << ", " << lower_y << "), (" + // << lower_x << ", " << upper_y << "), (" << upper_x << ", " << upper_y << ")" << G4endl; + //G4cout << "stddev values: " << stddev11 << ", " << stddev21 << ", " << stddev12 << ", " << stddev22 << G4endl; + + G4double interpolatedValue = (stddev11 * (upper_x - x) * (upper_y - y) + + stddev21 * (x - lower_x) * (upper_y - y) + + stddev12 * (upper_x - x) * (y - lower_y) + + stddev22 * (x - lower_x) * (y - lower_y)) / + ((upper_x - lower_x) * (upper_y - lower_y)); + + //G4cout << "Interpolated value: " << interpolatedValue << G4endl; + return interpolatedValue; + } + + From 9a2d4e248071925b4f046bcf3d2228b78c9cfa37 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Wed, 3 Jul 2024 14:54:30 +0200 Subject: [PATCH 09/22] Add modification --- .../include/GateDistributionFile.hh | 8 +- .../include/GateDistributionFileMessenger.hh | 1 - .../include/GateDistributionMessenger.hh | 5 +- .../include/GateSpatialResolution.hh | 29 +--- .../include/GateSpatialResolutionMessenger.hh | 5 +- .../digits_hits/include/GateVDistribution.hh | 2 - .../include/GateVDistributionArray.hh | 24 +-- .../digits_hits/src/GateDistributionFile.cc | 65 +------ .../src/GateDistributionFileMessenger.cc | 9 +- .../src/GateDistributionMessenger.cc | 7 +- .../digits_hits/src/GateSpatialResolution.cc | 106 ++---------- .../src/GateSpatialResolutionMessenger.cc | 48 +----- source/digits_hits/src/GateVDistribution.cc | 8 - .../digits_hits/src/GateVDistributionArray.cc | 161 +++--------------- 14 files changed, 74 insertions(+), 404 deletions(-) diff --git a/source/digits_hits/include/GateDistributionFile.hh b/source/digits_hits/include/GateDistributionFile.hh index cdf88cac0..e8d4e73d2 100644 --- a/source/digits_hits/include/GateDistributionFile.hh +++ b/source/digits_hits/include/GateDistributionFile.hh @@ -32,8 +32,7 @@ class GateDistributionFile : public GateVDistributionArray inline G4int GetColumnX() const {return m_column_for_X;} inline G4int GetColumnY() const {return m_column_for_Y;} - void Read(); - void ReadMatrix2d(); + void Read(); virtual void DescribeMyself(size_t indent); private: @@ -42,11 +41,6 @@ class GateDistributionFile : public GateVDistributionArray G4int m_column_for_X; G4int m_column_for_Y; GateDistributionFileMessenger* m_messenger; - std::map, double> stddevMap; - - std::vector xValues; - std::vector yValues; - }; diff --git a/source/digits_hits/include/GateDistributionFileMessenger.hh b/source/digits_hits/include/GateDistributionFileMessenger.hh index 281b88e3e..86ccc9426 100644 --- a/source/digits_hits/include/GateDistributionFileMessenger.hh +++ b/source/digits_hits/include/GateDistributionFileMessenger.hh @@ -33,7 +33,6 @@ class GateDistributionFileMessenger: public GateDistributionArrayMessenger G4UIcmdWithAnInteger *setColXCmd; G4UIcmdWithAnInteger *setColYCmd; G4UIcmdWithoutParameter *readCmd; - G4UIcmdWithoutParameter *read2DCmd; G4UIcmdWithoutParameter *autoXCmd; }; diff --git a/source/digits_hits/include/GateDistributionMessenger.hh b/source/digits_hits/include/GateDistributionMessenger.hh index 747c6852a..b33f3af72 100644 --- a/source/digits_hits/include/GateDistributionMessenger.hh +++ b/source/digits_hits/include/GateDistributionMessenger.hh @@ -30,10 +30,9 @@ class GateDistributionMessenger: public GateNamedObjectMessenger void SetNewValue(G4UIcommand* aCommand, G4String aString); void SetUnitX(const G4String& unitX); void SetUnitY(const G4String& unitY); - void SetUnitSigma(const G4String& unitSigma); inline G4String UnitCategoryX() const {return m_unitX.empty()?"":G4UIcommand::CategoryOf(m_unitX);} inline G4String UnitCategoryY() const {return m_unitY.empty()?"":G4UIcommand::CategoryOf(m_unitY);} - inline G4String UnitCategorySigma() const {return m_unitSigma.empty()?"":G4UIcommand::CategoryOf(m_unitSigma);} + private: G4String withUnity(G4double value,G4String category) const; G4UIcmdWithoutParameter *getMinX_Cmd ; @@ -42,10 +41,8 @@ class GateDistributionMessenger: public GateNamedObjectMessenger G4UIcmdWithoutParameter *getMaxY_Cmd ; G4UIcmdWithoutParameter *getRandom_Cmd ; G4UIcmdWithADoubleAndUnit *getValueCmd ; - G4UIcmdWithADoubleAndUnit *getSigmaCmd ; G4String m_unitX; G4String m_unitY; - G4String m_unitSigma; }; #endif diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh index 502bff961..171bf855e 100755 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/source/digits_hits/include/GateSpatialResolution.hh @@ -31,10 +31,7 @@ See LICENSE.md for further details #include "G4TouchableHistoryHandle.hh" #include "GateSinglesDigitizer.hh" -class GateVDistribution; - class GateSpatialResolution : public GateVDigitizerModule - { public: @@ -43,15 +40,9 @@ public: void Digitize() override; - - //! These functions return the resolution in use. G4double GetFWHM() { return m_fwhm; } - GateVDistribution* GetFWHMxdistrib() { return m_fwhmXdistrib; } - GateVDistribution* GetFWHMydistrib() { return m_fwhmYdistrib; } - GateVDistribution* GetFWHMxdistrib2D() { return m_fwhmXdistrib2D; } - GateVDistribution* GetFWHMydistrib2D() { return m_fwhmYdistrib2D; } - G4double GetFWHMx() { return m_fwhmX; } + G4double GetFWHMx() { return m_fwhmX; } G4double GetFWHMy() { return m_fwhmY; } G4double GetFWHMz() { return m_fwhmZ; } @@ -60,14 +51,11 @@ 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 SetFWHMxdistrib2D(GateVDistribution* dist) { m_fwhmXdistrib2D= dist; } - void SetFWHMydistrib2D(GateVDistribution* dist) { m_fwhmYdistrib2D = dist; } void SetFWHMx(G4double val) { m_fwhmX = val; } void SetFWHMy(G4double val) { m_fwhmY = val; } void SetFWHMz(G4double val) { m_fwhmZ = val; } + inline void ConfineInsideOfSmallestElement(const G4bool& value) { m_IsConfined = value; }; inline G4bool IsConfinedInsideOfSmallestElement() const { return m_IsConfined; } @@ -76,6 +64,7 @@ public: void UpdateVolumeID(); + //! Implementation of the pure virtual method declared by the base class GateClockDependent //! print-out the attributes specific of the blurring void DescribeMyself(size_t ); @@ -83,16 +72,9 @@ public: protected: G4double m_fwhm; - G4double m_fwhmX; - G4double m_fwhmY; G4double m_fwhmZ; - - GateVDistribution* m_fwhmXdistrib; - GateVDistribution* m_fwhmYdistrib; - GateVDistribution* m_fwhmXdistrib2D; - GateVDistribution* m_fwhmYdistrib2D; G4bool m_IsConfined; G4Navigator* m_Navigator; G4TouchableHistoryHandle m_Touchable; @@ -101,9 +83,10 @@ protected: private: - G4int m_systemDepth; + G4int m_systemDepth; + + GateDigi* m_outputDigi; - GateDigi* m_outputDigi;; GateSpatialResolutionMessenger *m_Messenger; GateDigiCollection* m_OutputDigiCollection; diff --git a/source/digits_hits/include/GateSpatialResolutionMessenger.hh b/source/digits_hits/include/GateSpatialResolutionMessenger.hh index db588a6e3..7b0d03c81 100755 --- a/source/digits_hits/include/GateSpatialResolutionMessenger.hh +++ b/source/digits_hits/include/GateSpatialResolutionMessenger.hh @@ -27,6 +27,7 @@ See LICENSE.md for further details #include "GateClockDependentMessenger.hh" class GateSpatialResolution; class G4UIcmdWithAString; + class GateSpatialResolutionMessenger : public GateClockDependentMessenger { public: @@ -42,10 +43,6 @@ private: G4UIcmdWithADouble* spresolutionCmd; G4UIcmdWithADouble* spresolutionXCmd; - G4UIcmdWithAString *spresolutionXdistribCmd; - G4UIcmdWithAString *spresolutionYdistribCmd; - G4UIcmdWithAString *spresolutionXdistrib2DCmd; - G4UIcmdWithAString *spresolutionYdistrib2DCmd; G4UIcmdWithADouble* spresolutionYCmd; G4UIcmdWithADouble* spresolutionZCmd; G4UIcmdWithABool* confineCmd; diff --git a/source/digits_hits/include/GateVDistribution.hh b/source/digits_hits/include/GateVDistribution.hh index ff524ccf2..7a428af02 100644 --- a/source/digits_hits/include/GateVDistribution.hh +++ b/source/digits_hits/include/GateVDistribution.hh @@ -41,8 +41,6 @@ class GateVDistribution : public GateNamedObject virtual G4double MaxX() const=0; virtual G4double MaxY() const=0; virtual G4double Value(G4double x) const=0; - virtual G4double Value2D(G4double x, G4double y) const; - // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const=0; diff --git a/source/digits_hits/include/GateVDistributionArray.hh b/source/digits_hits/include/GateVDistributionArray.hh index 3353efb80..226c23239 100644 --- a/source/digits_hits/include/GateVDistributionArray.hh +++ b/source/digits_hits/include/GateVDistributionArray.hh @@ -27,43 +27,35 @@ class GateVDistributionArray : public GateVDistribution inline void SetFactorX(G4double factor) {m_factorX=factor;} inline void SetFactorY(G4double factor) {m_factorY=factor;} - G4double Integral() const {return m_arrayRepartition.back();} + virtual G4double MinX() const; virtual G4double MinY() const; virtual G4double MaxX() const; virtual G4double MaxY() const; virtual G4double Value(G4double x) const; - virtual G4double RepartitionValue(G4double x) const; - G4double BilinearInterpolation(G4double x, G4double y) const ; // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const; size_t GetSize() const {return m_arrayX.size();} + void Clear(); void SetAutoStart(G4int start) {m_autoStart=start;} - virtual G4double Value2D(G4double x, G4double y) const; protected: - void InsertPoint(G4double x,G4double y, G4double sigma ); - void InsertPoint(G4double x, G4double y ); - void InsertPoint(G4double y ); + void InsertPoint(G4double x,G4double y); + void InsertPoint(G4double y); //! private function G4int FindIdxBefore(G4double x ,const std::vector& array) const; - - - void FillRepartition() ; - + void FillRepartition(); std::vector& GetArrayX() {return m_arrayX;} std::vector& GetArrayY() {return m_arrayY;} - std::vector& GetArrayRepartition() {return m_arrayRepartition;} private: //! private members std::vector m_arrayX; std::vector m_arrayY; - G4double m_minX; G4double m_minY; G4double m_maxX; @@ -71,13 +63,7 @@ class GateVDistributionArray : public GateVDistribution std::vector m_arrayRepartition; //! repartition function calculus G4double m_factorX; G4double m_factorY; - G4int m_autoStart; - - std::map, double> stddevMap; - std::vector xValues; // Stocker les valeurs x - std::vector yValues; // Stocker les valeurs y - }; diff --git a/source/digits_hits/src/GateDistributionFile.cc b/source/digits_hits/src/GateDistributionFile.cc index 973380e79..611c74092 100644 --- a/source/digits_hits/src/GateDistributionFile.cc +++ b/source/digits_hits/src/GateDistributionFile.cc @@ -20,7 +20,6 @@ GateDistributionFile::GateDistributionFile(const G4String& itsName) , m_FileName() , m_column_for_X(0) , m_column_for_Y(1) - { m_messenger = new GateDistributionFileMessenger(this,itsName); } @@ -31,10 +30,10 @@ GateDistributionFile::~GateDistributionFile() //___________________________________________________________________ void GateDistributionFile::DescribeMyself(size_t indent) { - //G4cout << GateTools::Indent(indent) - // <<"File : " << m_FileName - // <<'{' << m_column_for_X<<':'< xValues; - bool isFirstLine = true; - - while (std::getline(f, line)) { - std::stringstream iss(line); - - if (isFirstLine) { - // Process x values from the first line - G4double x; - - while (iss >> x) { - xValues.push_back(x); - if (iss.peek() == ',') iss.ignore(); - } - isFirstLine = false; - } else { - // Process y and stddev values for subsequent lines - G4double y; - iss >> y; - if (iss.peek() == ',') iss.ignore(); - // Read stddev values for each x value - for (size_t i = 0; i < xValues.size(); ++i) { - G4double stddev; - if (iss >> stddev) { - // Insert the (x, y) -> stddev pair into the map using InsertPoint - InsertPoint(xValues[i], y, stddev); - if (iss.peek() == ',') iss.ignore(); - } - } - } - } - - f.close(); - //G4cout << "Content of stddevMap:\n"; - //for (const auto& entry : stddevMap) { - // std::pair coordinates = entry.first; - // G4double stddev = entry.second; - //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl; - //} - - FillRepartition(); - } - diff --git a/source/digits_hits/src/GateDistributionFileMessenger.cc b/source/digits_hits/src/GateDistributionFileMessenger.cc index 107683d7b..aae879075 100644 --- a/source/digits_hits/src/GateDistributionFileMessenger.cc +++ b/source/digits_hits/src/GateDistributionFileMessenger.cc @@ -45,10 +45,7 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil guidance = "Do read the file"; readCmd = new G4UIcmdWithoutParameter(cmdName,this); readCmd->SetGuidance(guidance); - cmdName = GetDirectoryName()+"readMatrix2d"; - guidance = "Do read the file"; - read2DCmd = new G4UIcmdWithoutParameter(cmdName,this); - read2DCmd->SetGuidance(guidance); + cmdName = GetDirectoryName()+"autoX"; guidance = "Set automatic mode for X"; autoXCmd = new G4UIcmdWithoutParameter(cmdName,this); @@ -61,7 +58,6 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil GateDistributionFileMessenger::~GateDistributionFileMessenger() { delete readCmd; - delete read2DCmd; delete autoXCmd; delete setColYCmd; delete setColXCmd; @@ -87,9 +83,6 @@ void GateDistributionFileMessenger::SetNewValue(G4UIcommand* command,G4String ne } else if( command==readCmd ) { GetDistributionFile()->Read(); } - else if( command==read2DCmd ){ - GetDistributionFile()->ReadMatrix2d(); - } else GateDistributionArrayMessenger::SetNewValue(command,newValue); } diff --git a/source/digits_hits/src/GateDistributionMessenger.cc b/source/digits_hits/src/GateDistributionMessenger.cc index 0cd6bd444..6272a88c3 100644 --- a/source/digits_hits/src/GateDistributionMessenger.cc +++ b/source/digits_hits/src/GateDistributionMessenger.cc @@ -54,7 +54,6 @@ GateDistributionMessenger::GateDistributionMessenger(GateVDistribution* itsDistr getValueCmd->SetGuidance(guidance); getValueCmd->SetParameterName("x",false); - } @@ -68,7 +67,6 @@ GateDistributionMessenger::~GateDistributionMessenger() delete getMaxY_Cmd; delete getRandom_Cmd; delete getValueCmd; - } G4String GateDistributionMessenger::withUnity(G4double value,G4String category) const { @@ -84,8 +82,9 @@ void GateDistributionMessenger::SetNewValue(G4UIcommand* command,G4String newVal { if ( command==getValueCmd ){ G4double x = getValueCmd->GetNewDoubleValue(newValue); - G4double y = GetDistribution()->Value(x);} -else if( command==getMinX_Cmd ) { + G4double y = GetDistribution()->Value(x); + G4cout<GetObjectName()<<'('<MinX(); G4cout<GetObjectName()<<" MinX "<GetSD()->GetName()+"/SinglesDigitizer/"+digitizer->m_digitizerName+"/"+name,digitizer,digitizer->GetSD()), m_fwhm(0), m_fwhmX(0), - m_fwhmXdistrib(0), - m_fwhmYdistrib(0), - m_fwhmXdistrib2D(0), - m_fwhmYdistrib2D(0), m_fwhmY(0), m_fwhmZ(0), m_IsConfined(true), @@ -81,7 +72,7 @@ GateSpatialResolution::~GateSpatialResolution() void GateSpatialResolution::Digitize() { - if( m_fwhm!=0 && (m_fwhmX!=0 || m_fwhmY!=0 || m_fwhmZ!=0 )) + if( m_fwhm!=0 && (m_fwhmX!=0 || m_fwhmY!=0 || m_fwhmZ!=0 ) ) { GateError("***ERROR*** Spatial Resolution is ambiguous: you can set unique /fwhm for 3 axis OR set /fhwmX, /fhwmY, /fhwmZ!"); } @@ -90,22 +81,17 @@ void GateSpatialResolution::Digitize() G4double fwhmY; G4double fwhmZ; - if (m_fwhmX==0 && m_fwhmY==0 && m_fwhmZ==0 ) { fwhmX=m_fwhm; fwhmY=m_fwhm; fwhmZ=m_fwhm; - } else { fwhmX=m_fwhmX; fwhmY=m_fwhmY; fwhmZ=m_fwhmZ; - - - } @@ -153,69 +139,11 @@ void GateSpatialResolution::Digitize() G4double Px = P.x(); G4double Py = P.y(); G4double Pz = P.z(); - G4double stddevX, stddevY, stddevZ; - - if (m_fwhmXdistrib) { - // If the FWHM distribution for X is defined - if (m_fwhmYdistrib) { - // If the FWHM distribution for Y is also defined - stddevX = m_fwhmXdistrib->Value(P.x() * mm); - stddevY = m_fwhmYdistrib->Value(P.y() * mm); - } else if (m_fwhmY) { - // If the constant FWHM for Y is defined - stddevX = m_fwhmXdistrib->Value(P.x() * mm); - stddevY = fwhmY / GateConstants::fwhm_to_sigma; - } // If the 2D distribution for Y is defined - else if (m_fwhmYdistrib2D) { - // Use the constant value for X and the 2D distribution for Y - stddevX = fwhmX / GateConstants::fwhm_to_sigma; - stddevY = m_fwhmYdistrib2D->Value2D(P.x() * mm, P.y() * mm); - } - // If only the distribution for X is defined - else { - // Use the distribution for X and the constant value for Y - stddevX = m_fwhmXdistrib->Value(P.x() * mm); - stddevY = fwhmY / GateConstants::fwhm_to_sigma; - } - - } else if (m_fwhmXdistrib2D) { - // If the 2D FWHM distribution for X is defined - if (m_fwhmYdistrib) { - // If the FWHM distribution for Y is defined - stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); - stddevY = m_fwhmYdistrib->Value(P.y() * mm); - } else if (m_fwhmY) { - // If the constant FWHM for Y is defined - stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); - stddevY = fwhmY / GateConstants::fwhm_to_sigma; - } else { - // If neither the FWHM distribution nor constant for Y is defined - stddevX = stddevY = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); - } - } else if (m_fwhmYdistrib) { - // If the FWHM distribution for Y is defined - if (m_fwhmXdistrib) { - // If the FWHM distribution for X is also defined - stddevX = m_fwhmXdistrib->Value(P.x() * mm); - stddevY = m_fwhmYdistrib->Value(P.y() * mm); - } else if (m_fwhmX) { - // If the constant FWHM for X is defined - stddevX = fwhmX / GateConstants::fwhm_to_sigma; - stddevY = m_fwhmYdistrib->Value(P.y() * mm); - } else { - // If the 2D distribution for X is defined - stddevX = m_fwhmXdistrib2D->Value2D(P.x() * mm, P.y() * mm); - stddevY = m_fwhmYdistrib->Value(P.y() * mm); - } - } else { - // If neither the FWHM distributions for X nor Y are defined - stddevX = fwhmX / GateConstants::fwhm_to_sigma; - stddevY = fwhmY / GateConstants::fwhm_to_sigma; - } + G4double PxNew = G4RandGauss::shoot(Px,fwhmX/GateConstants::fwhm_to_sigma); + G4double PyNew = G4RandGauss::shoot(Py,fwhmY/GateConstants::fwhm_to_sigma); + G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); //TC + - G4double PxNew = G4RandGauss::shoot(Px,stddevX); - G4double PyNew = G4RandGauss::shoot(Py,stddevY); - G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); if (m_IsConfined) { //set the position on the border of the crystal @@ -223,22 +151,18 @@ void GateSpatialResolution::Digitize() inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kXAxis, limits, at, Xmin, Xmax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kYAxis, limits, at, Ymin, Ymax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kZAxis, limits, at, Zmin, Zmax); - if(PxNewXmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - - - + // m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC - //TC - //outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); + //TC outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); m_OutputDigiCollection->insert(m_outputDigi); - //G4cout<<"PxNew"<Xmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC + // + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC //Getting world Volume // Do not use from TransportationManager as it is not recommended - G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); G4VPhysicalVolume *WorldVolume = navigator->GetWorldVolume(); m_Navigator = new G4Navigator(); m_Navigator->SetWorldVolume(WorldVolume); @@ -281,7 +206,7 @@ void GateSpatialResolution::Digitize() } } - } + } else { if (nVerboseLevel>1) @@ -312,4 +237,5 @@ void GateSpatialResolution::DescribeMyself(size_t indent ) if(m_fwhm) G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhm << Gateendl; else - G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhmX <<" "<< m_fwhmY<< " "<SetGuidance("Set the resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmX"; - spresolutionXCmd= new G4UIcmdWithADouble(cmdName,this); - spresolutionXCmd->SetGuidance("Set the resolution "); + spresolutionXCmd = new G4UIcmdWithADouble(cmdName,this); + spresolutionXCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmY"; spresolutionYCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionYCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); @@ -43,18 +41,7 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol cmdName = GetDirectoryName() + "fwhmZ"; spresolutionZCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionZCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmXdistrib"; - spresolutionXdistribCmd = new G4UIcmdWithAString(cmdName,this); - spresolutionXdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmYdistrib"; - spresolutionYdistribCmd = new G4UIcmdWithAString(cmdName,this); - spresolutionYdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmXdistrib2D"; - spresolutionXdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); - spresolutionXdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmYdistrib2D"; - spresolutionYdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); - spresolutionYdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; confineCmd = new G4UIcmdWithABool(cmdName,this); confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); @@ -66,14 +53,11 @@ GateSpatialResolutionMessenger::~GateSpatialResolutionMessenger() { delete spresolutionCmd; delete spresolutionXCmd; - delete spresolutionXdistribCmd; - delete spresolutionYdistribCmd; - delete spresolutionXdistrib2DCmd; - delete spresolutionYdistrib2DCmd; delete spresolutionYCmd; delete spresolutionZCmd; delete confineCmd; + } @@ -81,24 +65,8 @@ void GateSpatialResolutionMessenger::SetNewValue(G4UIcommand * aCommand,G4String { if ( aCommand==spresolutionCmd ) { m_SpatialResolution->SetFWHM(spresolutionCmd->GetNewDoubleValue(newValue)); } - else if ( aCommand==spresolutionXdistribCmd ) - { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); - if (distrib)m_SpatialResolution->SetFWHMxdistrib(distrib); -} - else if ( aCommand==spresolutionYdistribCmd ) - { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); - if (distrib)m_SpatialResolution->SetFWHMydistrib(distrib); - } - 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==spresolutionXCmd ) - { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } + else if ( aCommand==spresolutionXCmd ) + { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionYCmd ) { m_SpatialResolution->SetFWHMy(spresolutionYCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionZCmd ) diff --git a/source/digits_hits/src/GateVDistribution.cc b/source/digits_hits/src/GateVDistribution.cc index 400bdc904..ac9155771 100644 --- a/source/digits_hits/src/GateVDistribution.cc +++ b/source/digits_hits/src/GateVDistribution.cc @@ -19,11 +19,3 @@ GateVDistribution::GateVDistribution(const G4String& itsName) GateVDistribution::~GateVDistribution() { } - -G4double GateVDistribution::Value2D(G4double x, G4double y)const { - - - G4cerr << "Warning: GateVDistribution::Value2D always returns 0 for Gaussian, flat (uniform), and exponential distributions." < #include #include "GateTools.hh" -#include + GateVDistributionArray::GateVDistributionArray(const G4String& itsName) : GateVDistribution(itsName) @@ -62,10 +60,9 @@ void GateVDistributionArray::Clear() m_maxX=-DBL_MAX; m_maxX=-DBL_MAX; } -//_________________________________________________________________ +//___________________________________________________________________ G4double GateVDistributionArray::Value(G4double x) const { - if ( m_arrayX.empty() ) return 0; if ( m_arrayX.size() == 1 ) return m_arrayY[0]/m_factorY; x*=m_factorX; @@ -77,26 +74,16 @@ G4double GateVDistributionArray::Value(G4double x) const x2 = m_arrayX[1]; y2 = m_arrayY[0]; } else if (idx==(G4int)m_arrayX.size()-1){ x1 = m_arrayX[m_arrayX.size()-2]; y1 = m_arrayY[m_arrayX.size()-2]; - x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; + x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; } else { x1 = m_arrayX[idx] ; y1 = m_arrayY[idx]; x2 = m_arrayX[idx+1]; y2 = m_arrayY[idx+1]; } -//G4cout << "m_arrayX:" << Gateendl; - for (size_t i = 0; i < m_arrayX.size(); ++i) { - G4cout << m_arrayX[i] << Gateendl; - } - //G4cout << "m_arrayY:" << Gateendl; - for (size_t i = 0; i < m_arrayY.size(); ++i) { - G4cout << m_arrayY[i] << Gateendl; - } - - //return the linear interpolation + // return the linear interpolation return (y1 + (x-x1)*(y2-y1)/(x2-x1))/m_factorY; } - -//__________________________________________________________________ +//___________________________________________________________________ G4double GateVDistributionArray::RepartitionValue(G4double x) const { if ( m_arrayX.empty() ) return 0; @@ -132,26 +119,24 @@ G4double GateVDistributionArray::ShootRandom() const return (x1 + (y-y1)*(x2-x1)/(y2-y1))/m_factorX; } //___________________________________________________________________ -void GateVDistributionArray:: FillRepartition() { - { - m_arrayRepartition.clear(); - if (m_arrayX.size()<2) return; - m_arrayRepartition.resize(m_arrayX.size()); - for (G4int i=0;i<(G4int)m_arrayX.size();++i){m_arrayRepartition[i]=0;} - for (G4int i=1;i<(G4int)m_arrayX.size();++i){ - G4double x1 = m_arrayX[i-1]; - G4double x2 = m_arrayX[i]; - G4double y1 = m_arrayY[i-1]; - G4double y2 = m_arrayY[i]; - m_arrayRepartition[i] = m_arrayRepartition[i-1] + (y1+y2)*(x2-x1)/2; - } - for (G4int i=0;i<(G4int)m_arrayX.size();++i){ - m_arrayRepartition[i]/=Integral(); - // G4cout<<"Repartition["<& array) const @@ -171,8 +156,8 @@ void GateVDistributionArray::InsertPoint(G4double x,G4double y) m_arrayY.push_back(y); } else if ( (i>=0) && m_arrayX[i] == x ) { //element replacement m_arrayY[i] = y; - //G4cerr<<"[GateDistributionArray::InsertPoint] WARNING : replacing value for " - // <m_maxY) m_maxY=y; } -void GateVDistributionArray::InsertPoint(G4double x, G4double y, G4double stddev) { - std::pair coordinates = std::make_pair(x, y); - stddevMap[coordinates] = stddev; - - - if (x < m_minX) m_minX = x; - if (x > m_maxX) m_maxX = x; - if (y < m_minY) m_minY = y; - if (y > m_maxY) m_maxY = y; - //G4cout << "Inserted: (x=" << x << ", y=" << y << ") -> stddev=" << stddev << std::endl; - - } - //___________________________________________________________________ void GateVDistributionArray::InsertPoint(G4double y) { if (GetSize()==0) InsertPoint(m_autoStart,y); else InsertPoint(m_maxX+1,y); } - - -G4double GateVDistributionArray::Value2D(double x, double y) const { - auto it = stddevMap.find(std::make_pair(x, y)); - if (it != stddevMap.end()) { - // If the point is found, return the standard deviation - G4double stddev= it->second; - // G4cout << "stddev for x=" << x << " and y=" << y << " is " << stddev<< G4endl; - return stddev; - } else { - //G4cout << "stddev not found for x=" << x << " and y=" << y << ". Performing bilinear interpolation." << G4endl; - G4double stddev = BilinearInterpolation(x, y); - return stddev; - } -} -G4double GateVDistributionArray::BilinearInterpolation(G4double x, G4double y) const { - auto lower_x_it = stddevMap.lower_bound(std::make_pair(x, -INFINITY)); - auto upper_x_it = stddevMap.upper_bound(std::make_pair(x, INFINITY)); - - if (lower_x_it == stddevMap.end() || upper_x_it == stddevMap.end()) { - G4cerr << "Error: Unable to find adjacent points in x direction." << G4endl; - return 0.0; - } - - // Adjust to make sure lower_x and upper_x are correct - G4double lower_x = (lower_x_it == stddevMap.begin()) ? lower_x_it->first.first : std::prev(lower_x_it)->first.first; - G4double upper_x = upper_x_it->first.first; - - if (lower_x == upper_x) { - G4cerr << "Error: Lower and upper x values are the same." << G4endl; - return 0.0; - } - - // Extract unique y values for these x coordinates - std::set y_values; - for (const auto& entry : stddevMap) { - if (entry.first.first == lower_x || entry.first.first == upper_x) { - y_values.insert(entry.first.second); - } - } - - if (y_values.empty()) { - G4cerr << "Error: No y values found for the given x values." << G4endl; - return 0.0; - } - - auto lower_y_it = y_values.lower_bound(y); - auto upper_y_it = y_values.upper_bound(y); - - if (lower_y_it == y_values.end() || upper_y_it == y_values.end()) { - G4cerr << "Error: Unable to find adjacent points in y direction." << G4endl; - return 0.0; - } - - G4double lower_y = (lower_y_it == y_values.begin()) ? *lower_y_it : *std::prev(lower_y_it); - G4double upper_y = *upper_y_it; - - if (lower_y == upper_y) { - G4cerr << "Error: Lower and upper y values are the same." << G4endl; - return 0.0; - } - - // Retrieve the sigma values for the four surrounding points - G4double stddev11 = stddevMap.at(std::make_pair(lower_x, lower_y)); - G4double stddev12 = stddevMap.at(std::make_pair(lower_x, upper_y)); - G4double stddev21 = stddevMap.at(std::make_pair(upper_x, lower_y)); - G4double stddev22 = stddevMap.at(std::make_pair(upper_x, upper_y)); - - //G4cout << "Interpolation points: (" << lower_x << ", " << lower_y << "), (" << upper_x << ", " << lower_y << "), (" - // << lower_x << ", " << upper_y << "), (" << upper_x << ", " << upper_y << ")" << G4endl; - //G4cout << "stddev values: " << stddev11 << ", " << stddev21 << ", " << stddev12 << ", " << stddev22 << G4endl; - - G4double interpolatedValue = (stddev11 * (upper_x - x) * (upper_y - y) + - stddev21 * (x - lower_x) * (upper_y - y) + - stddev12 * (upper_x - x) * (y - lower_y) + - stddev22 * (x - lower_x) * (y - lower_y)) / - ((upper_x - lower_x) * (upper_y - lower_y)); - - //G4cout << "Interpolated value: " << interpolatedValue << G4endl; - return interpolatedValue; - } - - From 6a52ee9b39bcabbc650e0b1c9ae2ca0398b915b8 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Wed, 3 Jul 2024 15:06:36 +0200 Subject: [PATCH 10/22] add comment --- .../digits_hits/src/GateSpatialResolution.cc | 122 +++++++++++++++--- 1 file changed, 101 insertions(+), 21 deletions(-) diff --git a/source/digits_hits/src/GateSpatialResolution.cc b/source/digits_hits/src/GateSpatialResolution.cc index 936cc2be5..7697f532e 100755 --- a/source/digits_hits/src/GateSpatialResolution.cc +++ b/source/digits_hits/src/GateSpatialResolution.cc @@ -1,4 +1,5 @@ + /*---------------------- Copyright (C): OpenGATE Collaboration @@ -21,9 +22,10 @@ - modified by Adrien Paillet 11/2022 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. + Implemented logic to determine standard deviations (stddevX, stddevY) based on defined 1D and 2D FWHM distributions for the X and Y axes. */ - #include "GateSpatialResolution.hh" #include "GateSpatialResolutionMessenger.hh" #include "GateDigi.hh" @@ -32,7 +34,6 @@ #include "GateObjectStore.hh" #include "GateConstants.hh" - #include "G4SystemOfUnits.hh" #include "G4EventManager.hh" #include "G4Event.hh" @@ -42,11 +43,18 @@ #include "G4UnitsTable.hh" #include "G4TransportationManager.hh" #include "G4Navigator.hh" +#include "GateVDistribution.hh" + + + GateSpatialResolution::GateSpatialResolution(GateSinglesDigitizer *digitizer, G4String name) :GateVDigitizerModule(name,"digitizerMgr/"+digitizer->GetSD()->GetName()+"/SinglesDigitizer/"+digitizer->m_digitizerName+"/"+name,digitizer,digitizer->GetSD()), m_fwhm(0), m_fwhmX(0), + m_fwhmXdistrib(0), + m_fwhmYdistrib(0), + m_fwhmXYdistrib2D(0), m_fwhmY(0), m_fwhmZ(0), m_IsConfined(true), @@ -70,28 +78,65 @@ GateSpatialResolution::~GateSpatialResolution() } -void GateSpatialResolution::Digitize() +void GateSpatialResolution::Digitize(){ + +if (m_fwhm != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmZ != 0 || m_fwhmXYdistrib2D != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 0)) { - if( m_fwhm!=0 && (m_fwhmX!=0 || m_fwhmY!=0 || m_fwhmZ!=0 ) ) - { - GateError("***ERROR*** Spatial Resolution is ambiguous: you can set unique /fwhm for 3 axis OR set /fhwmX, /fhwmY, /fhwmZ!"); - } + GateError("***ERROR*** Spatial Resolution is ambiguous: you can set a unique FWHM for all 3 axes OR set FWHM for X, Y, Z individually."); +} + + +if (m_fwhmXYdistrib2D != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 0)) +{ + GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM 2D distribution for (X,Y) OR set FWHM for X ,Y individually ."); +} + + +if (m_fwhmY != 0 && (m_fwhmYdistrib != 0)) +{ + GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM for Y distribution."); +} + +// Check if the spatial resolutions FWHM for the X axis are defined ambiguously +if (m_fwhmX != 0 && (m_fwhmXdistrib != 0)) +{ + GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM for X distribution."); +} + G4double fwhmX; G4double fwhmY; G4double fwhmZ; + if (m_fwhmX == 0 && m_fwhmY == 0 && m_fwhmZ == 0) + { + fwhmX = m_fwhm; + fwhmY = m_fwhm; + fwhmZ = m_fwhm; + } + else + { + fwhmX = m_fwhmX; + fwhmY = m_fwhmY; + fwhmZ = m_fwhmZ; + } + + if (m_fwhmX==0 && m_fwhmY==0 && m_fwhmZ==0 ) { fwhmX=m_fwhm; fwhmY=m_fwhm; fwhmZ=m_fwhm; + } else { fwhmX=m_fwhmX; fwhmY=m_fwhmY; fwhmZ=m_fwhmZ; + + + } @@ -139,11 +184,44 @@ void GateSpatialResolution::Digitize() G4double Px = P.x(); G4double Py = P.y(); G4double Pz = P.z(); - G4double PxNew = G4RandGauss::shoot(Px,fwhmX/GateConstants::fwhm_to_sigma); - G4double PyNew = G4RandGauss::shoot(Py,fwhmY/GateConstants::fwhm_to_sigma); - G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); //TC - + G4double stddevX, stddevY, stddevZ; + if (m_fwhmXYdistrib2D) { + // If the 2D FWHM distribution for X and Y is defined + + stddevX = m_fwhmXYdistrib2D->Value2D(P.x() * mm, P.y() * mm); + stddevY = stddevX; // Assuming the 2D distribution returns the same for both axes + } + else if (m_fwhmXdistrib) { + // If the FWHM distribution for X is defined + if (m_fwhmYdistrib) { + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else if (m_fwhmY) { + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } + else { + stddevX = m_fwhmXdistrib->Value(P.x() * mm); + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } + } else if (m_fwhmYdistrib) { + // If the FWHM distribution for Y is defined + if (m_fwhmX) { + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } else if (m_fwhmXdistrib) { + stddevX = m_fwhmXdistrib->Value(P.y() * mm); + stddevY = m_fwhmYdistrib->Value(P.y() * mm); + } + } else { + // If neither the FWHM distributions for X nor Y are defined + stddevX = fwhmX / GateConstants::fwhm_to_sigma; + stddevY = fwhmY / GateConstants::fwhm_to_sigma; + } + G4double PxNew = G4RandGauss::shoot(Px,stddevX); + G4double PyNew = G4RandGauss::shoot(Py,stddevY); + G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); if (m_IsConfined) { //set the position on the border of the crystal @@ -151,18 +229,22 @@ void GateSpatialResolution::Digitize() inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kXAxis, limits, at, Xmin, Xmax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kYAxis, limits, at, Ymin, Ymax); inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()->CalculateExtent(kZAxis, limits, at, Zmin, Zmax); + if(PxNewXmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - // + + + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC - //TC outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); + //TC + //outputPulse->SetGlobalPos(G4ThreeVector(PxNew,PyNew,PzNew)); m_OutputDigiCollection->insert(m_outputDigi); - + //G4cout<<"PxNew"<Xmax) PxNew=Xmax; if(PyNew>Ymax) PyNew=Ymax; if(PzNew>Zmax) PzNew=Zmax; - // - m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC + m_outputDigi->SetLocalPos(G4ThreeVector(PxNew,PyNew,PzNew)); //TC m_outputDigi->SetGlobalPos(m_outputDigi->GetVolumeID().MoveToAncestorVolumeFrame(m_outputDigi->GetLocalPos())); //TC //Getting world Volume // Do not use from TransportationManager as it is not recommended - G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); G4VPhysicalVolume *WorldVolume = navigator->GetWorldVolume(); m_Navigator = new G4Navigator(); m_Navigator->SetWorldVolume(WorldVolume); @@ -206,7 +287,7 @@ void GateSpatialResolution::Digitize() } } - } + } else { if (nVerboseLevel>1) @@ -237,5 +318,4 @@ void GateSpatialResolution::DescribeMyself(size_t indent ) if(m_fwhm) G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhm << Gateendl; else - G4cout << GateTools::Indent(indent) << "Spatial resolution : " << m_fwhmX <<" "<< m_fwhmY<< " "< Date: Wed, 3 Jul 2024 15:55:20 +0200 Subject: [PATCH 11/22] add modifications --- .../include/GateDistributionFile.hh | 8 +- .../include/GateDistributionFileMessenger.hh | 1 + .../include/GateDistributionMessenger.hh | 5 +- .../include/GateSpatialResolution.hh | 29 +++- .../include/GateSpatialResolutionMessenger.hh | 5 +- .../digits_hits/include/GateVDistribution.hh | 2 + .../include/GateVDistributionArray.hh | 24 ++- .../digits_hits/src/GateDistributionFile.cc | 69 +++++++- .../src/GateDistributionFileMessenger.cc | 9 +- .../src/GateDistributionMessenger.cc | 7 +- .../src/GateSpatialResolutionMessenger.cc | 46 ++++- source/digits_hits/src/GateVDistribution.cc | 8 + .../digits_hits/src/GateVDistributionArray.cc | 161 +++++++++++++++--- 13 files changed, 316 insertions(+), 58 deletions(-) diff --git a/source/digits_hits/include/GateDistributionFile.hh b/source/digits_hits/include/GateDistributionFile.hh index e8d4e73d2..cdf88cac0 100644 --- a/source/digits_hits/include/GateDistributionFile.hh +++ b/source/digits_hits/include/GateDistributionFile.hh @@ -32,7 +32,8 @@ class GateDistributionFile : public GateVDistributionArray inline G4int GetColumnX() const {return m_column_for_X;} inline G4int GetColumnY() const {return m_column_for_Y;} - void Read(); + void Read(); + void ReadMatrix2d(); virtual void DescribeMyself(size_t indent); private: @@ -41,6 +42,11 @@ class GateDistributionFile : public GateVDistributionArray G4int m_column_for_X; G4int m_column_for_Y; GateDistributionFileMessenger* m_messenger; + std::map, double> stddevMap; + + std::vector xValues; + std::vector yValues; + }; diff --git a/source/digits_hits/include/GateDistributionFileMessenger.hh b/source/digits_hits/include/GateDistributionFileMessenger.hh index 86ccc9426..281b88e3e 100644 --- a/source/digits_hits/include/GateDistributionFileMessenger.hh +++ b/source/digits_hits/include/GateDistributionFileMessenger.hh @@ -33,6 +33,7 @@ class GateDistributionFileMessenger: public GateDistributionArrayMessenger G4UIcmdWithAnInteger *setColXCmd; G4UIcmdWithAnInteger *setColYCmd; G4UIcmdWithoutParameter *readCmd; + G4UIcmdWithoutParameter *read2DCmd; G4UIcmdWithoutParameter *autoXCmd; }; diff --git a/source/digits_hits/include/GateDistributionMessenger.hh b/source/digits_hits/include/GateDistributionMessenger.hh index b33f3af72..747c6852a 100644 --- a/source/digits_hits/include/GateDistributionMessenger.hh +++ b/source/digits_hits/include/GateDistributionMessenger.hh @@ -30,9 +30,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger void SetNewValue(G4UIcommand* aCommand, G4String aString); void SetUnitX(const G4String& unitX); void SetUnitY(const G4String& unitY); + void SetUnitSigma(const G4String& unitSigma); inline G4String UnitCategoryX() const {return m_unitX.empty()?"":G4UIcommand::CategoryOf(m_unitX);} inline G4String UnitCategoryY() const {return m_unitY.empty()?"":G4UIcommand::CategoryOf(m_unitY);} - + inline G4String UnitCategorySigma() const {return m_unitSigma.empty()?"":G4UIcommand::CategoryOf(m_unitSigma);} private: G4String withUnity(G4double value,G4String category) const; G4UIcmdWithoutParameter *getMinX_Cmd ; @@ -41,8 +42,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger G4UIcmdWithoutParameter *getMaxY_Cmd ; G4UIcmdWithoutParameter *getRandom_Cmd ; G4UIcmdWithADoubleAndUnit *getValueCmd ; + G4UIcmdWithADoubleAndUnit *getSigmaCmd ; G4String m_unitX; G4String m_unitY; + G4String m_unitSigma; }; #endif diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh index 171bf855e..c28e7b37a 100755 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/source/digits_hits/include/GateSpatialResolution.hh @@ -31,7 +31,10 @@ See LICENSE.md for further details #include "G4TouchableHistoryHandle.hh" #include "GateSinglesDigitizer.hh" +class GateVDistribution; + class GateSpatialResolution : public GateVDigitizerModule + { public: @@ -40,9 +43,15 @@ public: void Digitize() override; + + //! These functions return the resolution in use. G4double GetFWHM() { return m_fwhm; } - G4double GetFWHMx() { return m_fwhmX; } + GateVDistribution* GetFWHMxdistrib() { return m_fwhmXdistrib; } + GateVDistribution* GetFWHMydistrib() { return m_fwhmYdistrib; } + GateVDistribution* GetFWHMxydistrib2D() { return m_fwhmXYdistrib2D; } + + G4double GetFWHMx() { return m_fwhmX; } G4double GetFWHMy() { return m_fwhmY; } G4double GetFWHMz() { return m_fwhmZ; } @@ -51,11 +60,14 @@ 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 SetFWHMxydistrib2D(GateVDistribution* dist) { m_fwhmXYdistrib2D= dist; } + void SetFWHMx(G4double val) { m_fwhmX = val; } void SetFWHMy(G4double val) { m_fwhmY = val; } void SetFWHMz(G4double val) { m_fwhmZ = val; } - inline void ConfineInsideOfSmallestElement(const G4bool& value) { m_IsConfined = value; }; inline G4bool IsConfinedInsideOfSmallestElement() const { return m_IsConfined; } @@ -64,7 +76,6 @@ public: void UpdateVolumeID(); - //! Implementation of the pure virtual method declared by the base class GateClockDependent //! print-out the attributes specific of the blurring void DescribeMyself(size_t ); @@ -72,9 +83,16 @@ public: protected: G4double m_fwhm; + G4double m_fwhmX; + G4double m_fwhmY; G4double m_fwhmZ; + + GateVDistribution* m_fwhmXdistrib; + GateVDistribution* m_fwhmYdistrib; + GateVDistribution* m_fwhmXYdistrib2D; + G4bool m_IsConfined; G4Navigator* m_Navigator; G4TouchableHistoryHandle m_Touchable; @@ -83,10 +101,9 @@ protected: private: - G4int m_systemDepth; - - GateDigi* m_outputDigi; + G4int m_systemDepth; + GateDigi* m_outputDigi;; GateSpatialResolutionMessenger *m_Messenger; GateDigiCollection* m_OutputDigiCollection; diff --git a/source/digits_hits/include/GateSpatialResolutionMessenger.hh b/source/digits_hits/include/GateSpatialResolutionMessenger.hh index 7b0d03c81..a140b97cd 100755 --- a/source/digits_hits/include/GateSpatialResolutionMessenger.hh +++ b/source/digits_hits/include/GateSpatialResolutionMessenger.hh @@ -27,7 +27,6 @@ See LICENSE.md for further details #include "GateClockDependentMessenger.hh" class GateSpatialResolution; class G4UIcmdWithAString; - class GateSpatialResolutionMessenger : public GateClockDependentMessenger { public: @@ -43,6 +42,10 @@ private: G4UIcmdWithADouble* spresolutionCmd; G4UIcmdWithADouble* spresolutionXCmd; + G4UIcmdWithAString *spresolutionXdistribCmd;// Command declaration for 1D X-resolution distribution + G4UIcmdWithAString *spresolutionYdistribCmd;// Command declaration for 1D Y-resolution distribution + G4UIcmdWithAString *spresolutionXYdistrib2DCmd; // Command declaration for 2D XY-resolution distribution + G4UIcmdWithADouble* spresolutionYCmd; G4UIcmdWithADouble* spresolutionZCmd; G4UIcmdWithABool* confineCmd; diff --git a/source/digits_hits/include/GateVDistribution.hh b/source/digits_hits/include/GateVDistribution.hh index 7a428af02..ff524ccf2 100644 --- a/source/digits_hits/include/GateVDistribution.hh +++ b/source/digits_hits/include/GateVDistribution.hh @@ -41,6 +41,8 @@ class GateVDistribution : public GateNamedObject virtual G4double MaxX() const=0; virtual G4double MaxY() const=0; virtual G4double Value(G4double x) const=0; + virtual G4double Value2D(G4double x, G4double y) const; + // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const=0; diff --git a/source/digits_hits/include/GateVDistributionArray.hh b/source/digits_hits/include/GateVDistributionArray.hh index 226c23239..3353efb80 100644 --- a/source/digits_hits/include/GateVDistributionArray.hh +++ b/source/digits_hits/include/GateVDistributionArray.hh @@ -27,35 +27,43 @@ class GateVDistributionArray : public GateVDistribution inline void SetFactorX(G4double factor) {m_factorX=factor;} inline void SetFactorY(G4double factor) {m_factorY=factor;} - G4double Integral() const {return m_arrayRepartition.back();} + G4double Integral() const {return m_arrayRepartition.back();} virtual G4double MinX() const; virtual G4double MinY() const; virtual G4double MaxX() const; virtual G4double MaxY() const; virtual G4double Value(G4double x) const; + virtual G4double RepartitionValue(G4double x) const; + G4double BilinearInterpolation(G4double x, G4double y) const ; // Returns a random number following the current distribution // should be optimised according to each distrbution type virtual G4double ShootRandom() const; size_t GetSize() const {return m_arrayX.size();} - void Clear(); void SetAutoStart(G4int start) {m_autoStart=start;} + virtual G4double Value2D(G4double x, G4double y) const; protected: - void InsertPoint(G4double x,G4double y); - void InsertPoint(G4double y); + void InsertPoint(G4double x,G4double y, G4double sigma ); + void InsertPoint(G4double x, G4double y ); + void InsertPoint(G4double y ); //! private function G4int FindIdxBefore(G4double x ,const std::vector& array) const; - void FillRepartition(); + + + void FillRepartition() ; + std::vector& GetArrayX() {return m_arrayX;} std::vector& GetArrayY() {return m_arrayY;} + std::vector& GetArrayRepartition() {return m_arrayRepartition;} private: //! private members std::vector m_arrayX; std::vector m_arrayY; + G4double m_minX; G4double m_minY; G4double m_maxX; @@ -63,7 +71,13 @@ class GateVDistributionArray : public GateVDistribution std::vector m_arrayRepartition; //! repartition function calculus G4double m_factorX; G4double m_factorY; + G4int m_autoStart; + + std::map, double> stddevMap; + std::vector xValues; // Stocker les valeurs x + std::vector yValues; // Stocker les valeurs y + }; diff --git a/source/digits_hits/src/GateDistributionFile.cc b/source/digits_hits/src/GateDistributionFile.cc index 611c74092..93602c3cd 100644 --- a/source/digits_hits/src/GateDistributionFile.cc +++ b/source/digits_hits/src/GateDistributionFile.cc @@ -20,6 +20,7 @@ GateDistributionFile::GateDistributionFile(const G4String& itsName) , m_FileName() , m_column_for_X(0) , m_column_for_Y(1) + { m_messenger = new GateDistributionFileMessenger(this,itsName); } @@ -30,16 +31,17 @@ GateDistributionFile::~GateDistributionFile() //___________________________________________________________________ void GateDistributionFile::DescribeMyself(size_t indent) { - G4cout << GateTools::Indent(indent) - <<"File : " << m_FileName - <<'{' << m_column_for_X<<':'< xValues; + bool isFirstLine = true; + + while (std::getline(f, line)) { + std::stringstream iss(line); + + if (isFirstLine) { + // Process x values from the first line + G4double x; + + while (iss >> x) { + xValues.push_back(x); + if (iss.peek() == ',') iss.ignore(); + } + isFirstLine = false; + } else { + // Process y and stddev values for subsequent lines + G4double y; + iss >> y; + if (iss.peek() == ',') iss.ignore(); + // Read stddev values for each x value + for (size_t i = 0; i < xValues.size(); ++i) { + G4double stddev; + if (iss >> stddev) { + // Insert the (x, y) -> stddev pair into the map using InsertPoint + InsertPoint(xValues[i], y, stddev); + if (iss.peek() == ',') iss.ignore(); + } + } + } + } + + f.close(); + //G4cout << "Content of stddevMap:\n"; + //for (const auto& entry : stddevMap) { + // std::pair coordinates = entry.first; + // G4double stddev = entry.second; + //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl; + //} + + FillRepartition(); + } + diff --git a/source/digits_hits/src/GateDistributionFileMessenger.cc b/source/digits_hits/src/GateDistributionFileMessenger.cc index aae879075..107683d7b 100644 --- a/source/digits_hits/src/GateDistributionFileMessenger.cc +++ b/source/digits_hits/src/GateDistributionFileMessenger.cc @@ -45,7 +45,10 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil guidance = "Do read the file"; readCmd = new G4UIcmdWithoutParameter(cmdName,this); readCmd->SetGuidance(guidance); - + cmdName = GetDirectoryName()+"readMatrix2d"; + guidance = "Do read the file"; + read2DCmd = new G4UIcmdWithoutParameter(cmdName,this); + read2DCmd->SetGuidance(guidance); cmdName = GetDirectoryName()+"autoX"; guidance = "Set automatic mode for X"; autoXCmd = new G4UIcmdWithoutParameter(cmdName,this); @@ -58,6 +61,7 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil GateDistributionFileMessenger::~GateDistributionFileMessenger() { delete readCmd; + delete read2DCmd; delete autoXCmd; delete setColYCmd; delete setColXCmd; @@ -83,6 +87,9 @@ void GateDistributionFileMessenger::SetNewValue(G4UIcommand* command,G4String ne } else if( command==readCmd ) { GetDistributionFile()->Read(); } + else if( command==read2DCmd ){ + GetDistributionFile()->ReadMatrix2d(); + } else GateDistributionArrayMessenger::SetNewValue(command,newValue); } diff --git a/source/digits_hits/src/GateDistributionMessenger.cc b/source/digits_hits/src/GateDistributionMessenger.cc index 6272a88c3..0cd6bd444 100644 --- a/source/digits_hits/src/GateDistributionMessenger.cc +++ b/source/digits_hits/src/GateDistributionMessenger.cc @@ -54,6 +54,7 @@ GateDistributionMessenger::GateDistributionMessenger(GateVDistribution* itsDistr getValueCmd->SetGuidance(guidance); getValueCmd->SetParameterName("x",false); + } @@ -67,6 +68,7 @@ GateDistributionMessenger::~GateDistributionMessenger() delete getMaxY_Cmd; delete getRandom_Cmd; delete getValueCmd; + } G4String GateDistributionMessenger::withUnity(G4double value,G4String category) const { @@ -82,9 +84,8 @@ void GateDistributionMessenger::SetNewValue(G4UIcommand* command,G4String newVal { if ( command==getValueCmd ){ G4double x = getValueCmd->GetNewDoubleValue(newValue); - G4double y = GetDistribution()->Value(x); - G4cout<GetObjectName()<<'('<Value(x);} +else if( command==getMinX_Cmd ) { G4double x = GetDistribution()->MinX(); G4cout<GetObjectName()<<" MinX "<SetGuidance("Set the resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmX"; - spresolutionXCmd = new G4UIcmdWithADouble(cmdName,this); - spresolutionXCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); - + spresolutionXCmd= new G4UIcmdWithADouble(cmdName,this); + spresolutionXCmd->SetGuidance("Set the resolution "); cmdName = GetDirectoryName() + "fwhmY"; spresolutionYCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionYCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); @@ -41,8 +43,17 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol cmdName = GetDirectoryName() + "fwhmZ"; spresolutionZCmd = new G4UIcmdWithADouble(cmdName,this); spresolutionZCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmXdistrib"; + spresolutionXdistribCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionXdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmYdistrib"; + spresolutionYdistribCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionYdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmXYdistrib2D"; + spresolutionXYdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); + spresolutionXYdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); + cmdName = GetDirectoryName() + "fwhmYdistrib2D"; - cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; confineCmd = new G4UIcmdWithABool(cmdName,this); confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); @@ -53,11 +64,14 @@ GateSpatialResolutionMessenger::~GateSpatialResolutionMessenger() { delete spresolutionCmd; delete spresolutionXCmd; + delete spresolutionXdistribCmd; + delete spresolutionYdistribCmd; + delete spresolutionXYdistrib2DCmd; + delete spresolutionYCmd; delete spresolutionZCmd; delete confineCmd; - } @@ -65,8 +79,24 @@ void GateSpatialResolutionMessenger::SetNewValue(G4UIcommand * aCommand,G4String { if ( aCommand==spresolutionCmd ) { m_SpatialResolution->SetFWHM(spresolutionCmd->GetNewDoubleValue(newValue)); } - else if ( aCommand==spresolutionXCmd ) - { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } + // Handle command for 1D X-distribution resolution + else if ( aCommand==spresolutionXdistribCmd ) + { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib)m_SpatialResolution->SetFWHMxdistrib(distrib); + } + // Handle command for 1D Y-distribution resolution + else if ( aCommand==spresolutionYdistribCmd ) + { GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib)m_SpatialResolution->SetFWHMydistrib(distrib); + } + // Handle command for 2D XY-distribution resolution + else if (aCommand == spresolutionXYdistrib2DCmd) + {GateVDistribution* distrib = (GateVDistribution*)GateDistributionListManager::GetInstance()->FindElementByBaseName(newValue); + if (distrib) m_SpatialResolution->SetFWHMxydistrib2D(distrib); + } + + else if ( aCommand==spresolutionXCmd ) + { m_SpatialResolution->SetFWHMx(spresolutionXCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionYCmd ) { m_SpatialResolution->SetFWHMy(spresolutionYCmd->GetNewDoubleValue(newValue)); } else if ( aCommand==spresolutionZCmd ) diff --git a/source/digits_hits/src/GateVDistribution.cc b/source/digits_hits/src/GateVDistribution.cc index ac9155771..400bdc904 100644 --- a/source/digits_hits/src/GateVDistribution.cc +++ b/source/digits_hits/src/GateVDistribution.cc @@ -19,3 +19,11 @@ GateVDistribution::GateVDistribution(const G4String& itsName) GateVDistribution::~GateVDistribution() { } + +G4double GateVDistribution::Value2D(G4double x, G4double y)const { + + + G4cerr << "Warning: GateVDistribution::Value2D always returns 0 for Gaussian, flat (uniform), and exponential distributions." < #include #include "GateTools.hh" - +#include GateVDistributionArray::GateVDistributionArray(const G4String& itsName) : GateVDistribution(itsName) @@ -60,9 +62,10 @@ void GateVDistributionArray::Clear() m_maxX=-DBL_MAX; m_maxX=-DBL_MAX; } -//___________________________________________________________________ +//_________________________________________________________________ G4double GateVDistributionArray::Value(G4double x) const { + if ( m_arrayX.empty() ) return 0; if ( m_arrayX.size() == 1 ) return m_arrayY[0]/m_factorY; x*=m_factorX; @@ -74,16 +77,26 @@ G4double GateVDistributionArray::Value(G4double x) const x2 = m_arrayX[1]; y2 = m_arrayY[0]; } else if (idx==(G4int)m_arrayX.size()-1){ x1 = m_arrayX[m_arrayX.size()-2]; y1 = m_arrayY[m_arrayX.size()-2]; - x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; + x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; } else { x1 = m_arrayX[idx] ; y1 = m_arrayY[idx]; x2 = m_arrayX[idx+1]; y2 = m_arrayY[idx+1]; } +//G4cout << "m_arrayX:" << Gateendl; + for (size_t i = 0; i < m_arrayX.size(); ++i) { + // G4cout << m_arrayX[i] << Gateendl; + } - // return the linear interpolation + //G4cout << "m_arrayY:" << Gateendl; + for (size_t i = 0; i < m_arrayY.size(); ++i) { + // G4cout << m_arrayY[i] << Gateendl; + } + + //return the linear interpolation return (y1 + (x-x1)*(y2-y1)/(x2-x1))/m_factorY; } -//___________________________________________________________________ + +//__________________________________________________________________ G4double GateVDistributionArray::RepartitionValue(G4double x) const { if ( m_arrayX.empty() ) return 0; @@ -119,24 +132,26 @@ G4double GateVDistributionArray::ShootRandom() const return (x1 + (y-y1)*(x2-x1)/(y2-y1))/m_factorX; } //___________________________________________________________________ -void GateVDistributionArray::FillRepartition() -{ - m_arrayRepartition.clear(); - if (m_arrayX.size()<2) return; - m_arrayRepartition.resize(m_arrayX.size()); - for (G4int i=0;i<(G4int)m_arrayX.size();++i){m_arrayRepartition[i]=0;} - for (G4int i=1;i<(G4int)m_arrayX.size();++i){ - G4double x1 = m_arrayX[i-1]; - G4double x2 = m_arrayX[i]; - G4double y1 = m_arrayY[i-1]; - G4double y2 = m_arrayY[i]; - m_arrayRepartition[i] = m_arrayRepartition[i-1] + (y1+y2)*(x2-x1)/2; - } - for (G4int i=0;i<(G4int)m_arrayX.size();++i){ - m_arrayRepartition[i]/=Integral(); -// G4cout<<"Repartition["<& array) const @@ -156,8 +171,8 @@ void GateVDistributionArray::InsertPoint(G4double x,G4double y) m_arrayY.push_back(y); } else if ( (i>=0) && m_arrayX[i] == x ) { //element replacement m_arrayY[i] = y; - G4cerr<<"[GateDistributionArray::InsertPoint] WARNING : replacing value for " - <m_maxY) m_maxY=y; } +void GateVDistributionArray::InsertPoint(G4double x, G4double y, G4double stddev) { + std::pair coordinates = std::make_pair(x, y); + stddevMap[coordinates] = stddev; + + + if (x < m_minX) m_minX = x; + if (x > m_maxX) m_maxX = x; + if (y < m_minY) m_minY = y; + if (y > m_maxY) m_maxY = y; + //G4cout << "Inserted: (x=" << x << ", y=" << y << ") -> stddev=" << stddev << std::endl; + + } + //___________________________________________________________________ void GateVDistributionArray::InsertPoint(G4double y) { if (GetSize()==0) InsertPoint(m_autoStart,y); else InsertPoint(m_maxX+1,y); } + + +G4double GateVDistributionArray::Value2D(double x, double y) const { + auto it = stddevMap.find(std::make_pair(x, y)); + if (it != stddevMap.end()) { + // If the point is found, return the standard deviation + G4double stddev= it->second; + // G4cout << "stddev for x=" << x << " and y=" << y << " is " << stddev<< G4endl; + return stddev; + } else { + //G4cout << "stddev not found for x=" << x << " and y=" << y << ". Performing bilinear interpolation." << G4endl; + G4double stddev = BilinearInterpolation(x, y); + return stddev; + } +} +G4double GateVDistributionArray::BilinearInterpolation(G4double x, G4double y) const { + auto lower_x_it = stddevMap.lower_bound(std::make_pair(x, -INFINITY)); + auto upper_x_it = stddevMap.upper_bound(std::make_pair(x, INFINITY)); + + if (lower_x_it == stddevMap.end() || upper_x_it == stddevMap.end()) { + G4cerr << "Error: Unable to find adjacent points in x direction." << G4endl; + return 0.0; + } + + // Adjust to make sure lower_x and upper_x are correct + G4double lower_x = (lower_x_it == stddevMap.begin()) ? lower_x_it->first.first : std::prev(lower_x_it)->first.first; + G4double upper_x = upper_x_it->first.first; + + if (lower_x == upper_x) { + G4cerr << "Error: Lower and upper x values are the same." << G4endl; + return 0.0; + } + + // Extract unique y values for these x coordinates + std::set y_values; + for (const auto& entry : stddevMap) { + if (entry.first.first == lower_x || entry.first.first == upper_x) { + y_values.insert(entry.first.second); + } + } + + if (y_values.empty()) { + G4cerr << "Error: No y values found for the given x values." << G4endl; + return 0.0; + } + + auto lower_y_it = y_values.lower_bound(y); + auto upper_y_it = y_values.upper_bound(y); + + if (lower_y_it == y_values.end() || upper_y_it == y_values.end()) { + G4cerr << "Error: Unable to find adjacent points in y direction." << G4endl; + return 0.0; + } + + G4double lower_y = (lower_y_it == y_values.begin()) ? *lower_y_it : *std::prev(lower_y_it); + G4double upper_y = *upper_y_it; + + if (lower_y == upper_y) { + G4cerr << "Error: Lower and upper y values are the same." << G4endl; + return 0.0; + } + + // Retrieve the sigma values for the four surrounding points + G4double stddev11 = stddevMap.at(std::make_pair(lower_x, lower_y)); + G4double stddev12 = stddevMap.at(std::make_pair(lower_x, upper_y)); + G4double stddev21 = stddevMap.at(std::make_pair(upper_x, lower_y)); + G4double stddev22 = stddevMap.at(std::make_pair(upper_x, upper_y)); + + //G4cout << "Interpolation points: (" << lower_x << ", " << lower_y << "), (" << upper_x << ", " << lower_y << "), (" + // << lower_x << ", " << upper_y << "), (" << upper_x << ", " << upper_y << ")" << G4endl; + //G4cout << "stddev values: " << stddev11 << ", " << stddev21 << ", " << stddev12 << ", " << stddev22 << G4endl; + + G4double interpolatedValue = (stddev11 * (upper_x - x) * (upper_y - y) + + stddev21 * (x - lower_x) * (upper_y - y) + + stddev12 * (upper_x - x) * (y - lower_y) + + stddev22 * (x - lower_x) * (y - lower_y)) / + ((upper_x - lower_x) * (upper_y - lower_y)); + + //G4cout << "Interpolated value: " << interpolatedValue << G4endl; + return interpolatedValue; + } + + From 1ed4af15e84549960e7a9808c511f9b8ce56e303 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Wed, 3 Jul 2024 16:14:13 +0200 Subject: [PATCH 12/22] add modifications --- .../include/GateDistributionMessenger.hh | 8 +- .../digits_hits/src/GateVDistributionArray.cc | 161 +++--------------- 2 files changed, 29 insertions(+), 140 deletions(-) diff --git a/source/digits_hits/include/GateDistributionMessenger.hh b/source/digits_hits/include/GateDistributionMessenger.hh index 747c6852a..74f2c7805 100644 --- a/source/digits_hits/include/GateDistributionMessenger.hh +++ b/source/digits_hits/include/GateDistributionMessenger.hh @@ -30,10 +30,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger void SetNewValue(G4UIcommand* aCommand, G4String aString); void SetUnitX(const G4String& unitX); void SetUnitY(const G4String& unitY); - void SetUnitSigma(const G4String& unitSigma); + inline G4String UnitCategoryX() const {return m_unitX.empty()?"":G4UIcommand::CategoryOf(m_unitX);} inline G4String UnitCategoryY() const {return m_unitY.empty()?"":G4UIcommand::CategoryOf(m_unitY);} - inline G4String UnitCategorySigma() const {return m_unitSigma.empty()?"":G4UIcommand::CategoryOf(m_unitSigma);} + private: G4String withUnity(G4double value,G4String category) const; G4UIcmdWithoutParameter *getMinX_Cmd ; @@ -42,10 +42,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger G4UIcmdWithoutParameter *getMaxY_Cmd ; G4UIcmdWithoutParameter *getRandom_Cmd ; G4UIcmdWithADoubleAndUnit *getValueCmd ; - G4UIcmdWithADoubleAndUnit *getSigmaCmd ; + G4String m_unitX; G4String m_unitY; - G4String m_unitSigma; + }; #endif diff --git a/source/digits_hits/src/GateVDistributionArray.cc b/source/digits_hits/src/GateVDistributionArray.cc index d9a59fe19..7134ba322 100644 --- a/source/digits_hits/src/GateVDistributionArray.cc +++ b/source/digits_hits/src/GateVDistributionArray.cc @@ -8,13 +8,11 @@ See LICENSE.md for further details #include "GateVDistributionArray.hh" -#include "GateDistributionFile.hh" - #include #include #include "GateTools.hh" -#include + GateVDistributionArray::GateVDistributionArray(const G4String& itsName) : GateVDistribution(itsName) @@ -62,10 +60,9 @@ void GateVDistributionArray::Clear() m_maxX=-DBL_MAX; m_maxX=-DBL_MAX; } -//_________________________________________________________________ +//___________________________________________________________________ G4double GateVDistributionArray::Value(G4double x) const { - if ( m_arrayX.empty() ) return 0; if ( m_arrayX.size() == 1 ) return m_arrayY[0]/m_factorY; x*=m_factorX; @@ -77,26 +74,16 @@ G4double GateVDistributionArray::Value(G4double x) const x2 = m_arrayX[1]; y2 = m_arrayY[0]; } else if (idx==(G4int)m_arrayX.size()-1){ x1 = m_arrayX[m_arrayX.size()-2]; y1 = m_arrayY[m_arrayX.size()-2]; - x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; + x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; } else { x1 = m_arrayX[idx] ; y1 = m_arrayY[idx]; x2 = m_arrayX[idx+1]; y2 = m_arrayY[idx+1]; } -//G4cout << "m_arrayX:" << Gateendl; - for (size_t i = 0; i < m_arrayX.size(); ++i) { - // G4cout << m_arrayX[i] << Gateendl; - } - //G4cout << "m_arrayY:" << Gateendl; - for (size_t i = 0; i < m_arrayY.size(); ++i) { - // G4cout << m_arrayY[i] << Gateendl; - } - - //return the linear interpolation + // return the linear interpolation return (y1 + (x-x1)*(y2-y1)/(x2-x1))/m_factorY; } - -//__________________________________________________________________ +//___________________________________________________________________ G4double GateVDistributionArray::RepartitionValue(G4double x) const { if ( m_arrayX.empty() ) return 0; @@ -132,26 +119,24 @@ G4double GateVDistributionArray::ShootRandom() const return (x1 + (y-y1)*(x2-x1)/(y2-y1))/m_factorX; } //___________________________________________________________________ -void GateVDistributionArray:: FillRepartition() { - { - m_arrayRepartition.clear(); - if (m_arrayX.size()<2) return; - m_arrayRepartition.resize(m_arrayX.size()); - for (G4int i=0;i<(G4int)m_arrayX.size();++i){m_arrayRepartition[i]=0;} - for (G4int i=1;i<(G4int)m_arrayX.size();++i){ - G4double x1 = m_arrayX[i-1]; - G4double x2 = m_arrayX[i]; - G4double y1 = m_arrayY[i-1]; - G4double y2 = m_arrayY[i]; - m_arrayRepartition[i] = m_arrayRepartition[i-1] + (y1+y2)*(x2-x1)/2; - } - for (G4int i=0;i<(G4int)m_arrayX.size();++i){ - m_arrayRepartition[i]/=Integral(); - // G4cout<<"Repartition["<& array) const @@ -171,8 +156,8 @@ void GateVDistributionArray::InsertPoint(G4double x,G4double y) m_arrayY.push_back(y); } else if ( (i>=0) && m_arrayX[i] == x ) { //element replacement m_arrayY[i] = y; - //G4cerr<<"[GateDistributionArray::InsertPoint] WARNING : replacing value for " - // <m_maxY) m_maxY=y; } -void GateVDistributionArray::InsertPoint(G4double x, G4double y, G4double stddev) { - std::pair coordinates = std::make_pair(x, y); - stddevMap[coordinates] = stddev; - - - if (x < m_minX) m_minX = x; - if (x > m_maxX) m_maxX = x; - if (y < m_minY) m_minY = y; - if (y > m_maxY) m_maxY = y; - //G4cout << "Inserted: (x=" << x << ", y=" << y << ") -> stddev=" << stddev << std::endl; - - } - //___________________________________________________________________ void GateVDistributionArray::InsertPoint(G4double y) { if (GetSize()==0) InsertPoint(m_autoStart,y); else InsertPoint(m_maxX+1,y); } - - -G4double GateVDistributionArray::Value2D(double x, double y) const { - auto it = stddevMap.find(std::make_pair(x, y)); - if (it != stddevMap.end()) { - // If the point is found, return the standard deviation - G4double stddev= it->second; - // G4cout << "stddev for x=" << x << " and y=" << y << " is " << stddev<< G4endl; - return stddev; - } else { - //G4cout << "stddev not found for x=" << x << " and y=" << y << ". Performing bilinear interpolation." << G4endl; - G4double stddev = BilinearInterpolation(x, y); - return stddev; - } -} -G4double GateVDistributionArray::BilinearInterpolation(G4double x, G4double y) const { - auto lower_x_it = stddevMap.lower_bound(std::make_pair(x, -INFINITY)); - auto upper_x_it = stddevMap.upper_bound(std::make_pair(x, INFINITY)); - - if (lower_x_it == stddevMap.end() || upper_x_it == stddevMap.end()) { - G4cerr << "Error: Unable to find adjacent points in x direction." << G4endl; - return 0.0; - } - - // Adjust to make sure lower_x and upper_x are correct - G4double lower_x = (lower_x_it == stddevMap.begin()) ? lower_x_it->first.first : std::prev(lower_x_it)->first.first; - G4double upper_x = upper_x_it->first.first; - - if (lower_x == upper_x) { - G4cerr << "Error: Lower and upper x values are the same." << G4endl; - return 0.0; - } - - // Extract unique y values for these x coordinates - std::set y_values; - for (const auto& entry : stddevMap) { - if (entry.first.first == lower_x || entry.first.first == upper_x) { - y_values.insert(entry.first.second); - } - } - - if (y_values.empty()) { - G4cerr << "Error: No y values found for the given x values." << G4endl; - return 0.0; - } - - auto lower_y_it = y_values.lower_bound(y); - auto upper_y_it = y_values.upper_bound(y); - - if (lower_y_it == y_values.end() || upper_y_it == y_values.end()) { - G4cerr << "Error: Unable to find adjacent points in y direction." << G4endl; - return 0.0; - } - - G4double lower_y = (lower_y_it == y_values.begin()) ? *lower_y_it : *std::prev(lower_y_it); - G4double upper_y = *upper_y_it; - - if (lower_y == upper_y) { - G4cerr << "Error: Lower and upper y values are the same." << G4endl; - return 0.0; - } - - // Retrieve the sigma values for the four surrounding points - G4double stddev11 = stddevMap.at(std::make_pair(lower_x, lower_y)); - G4double stddev12 = stddevMap.at(std::make_pair(lower_x, upper_y)); - G4double stddev21 = stddevMap.at(std::make_pair(upper_x, lower_y)); - G4double stddev22 = stddevMap.at(std::make_pair(upper_x, upper_y)); - - //G4cout << "Interpolation points: (" << lower_x << ", " << lower_y << "), (" << upper_x << ", " << lower_y << "), (" - // << lower_x << ", " << upper_y << "), (" << upper_x << ", " << upper_y << ")" << G4endl; - //G4cout << "stddev values: " << stddev11 << ", " << stddev21 << ", " << stddev12 << ", " << stddev22 << G4endl; - - G4double interpolatedValue = (stddev11 * (upper_x - x) * (upper_y - y) + - stddev21 * (x - lower_x) * (upper_y - y) + - stddev12 * (upper_x - x) * (y - lower_y) + - stddev22 * (x - lower_x) * (y - lower_y)) / - ((upper_x - lower_x) * (upper_y - lower_y)); - - //G4cout << "Interpolated value: " << interpolatedValue << G4endl; - return interpolatedValue; - } - - From 9ffa42da44cd71822c9b4e1fc0403d22246d8752 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Wed, 3 Jul 2024 16:23:03 +0200 Subject: [PATCH 13/22] add modification for GateVDistributionArray.cc --- .../digits_hits/src/GateVDistributionArray.cc | 163 +++++++++++++++--- 1 file changed, 137 insertions(+), 26 deletions(-) diff --git a/source/digits_hits/src/GateVDistributionArray.cc b/source/digits_hits/src/GateVDistributionArray.cc index 7134ba322..a934ae336 100644 --- a/source/digits_hits/src/GateVDistributionArray.cc +++ b/source/digits_hits/src/GateVDistributionArray.cc @@ -8,11 +8,13 @@ See LICENSE.md for further details #include "GateVDistributionArray.hh" +#include "GateDistributionFile.hh" + #include #include #include "GateTools.hh" - +#include GateVDistributionArray::GateVDistributionArray(const G4String& itsName) : GateVDistribution(itsName) @@ -60,9 +62,10 @@ void GateVDistributionArray::Clear() m_maxX=-DBL_MAX; m_maxX=-DBL_MAX; } -//___________________________________________________________________ +//_________________________________________________________________ G4double GateVDistributionArray::Value(G4double x) const { + if ( m_arrayX.empty() ) return 0; if ( m_arrayX.size() == 1 ) return m_arrayY[0]/m_factorY; x*=m_factorX; @@ -74,16 +77,26 @@ G4double GateVDistributionArray::Value(G4double x) const x2 = m_arrayX[1]; y2 = m_arrayY[0]; } else if (idx==(G4int)m_arrayX.size()-1){ x1 = m_arrayX[m_arrayX.size()-2]; y1 = m_arrayY[m_arrayX.size()-2]; - x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; + x2 = m_arrayX[m_arrayX.size()-1]; y2 = m_arrayY[m_arrayX.size()-1]; } else { x1 = m_arrayX[idx] ; y1 = m_arrayY[idx]; x2 = m_arrayX[idx+1]; y2 = m_arrayY[idx+1]; - } + } + //G4cout << "m_arrayX:" << Gateendl; + //for (size_t i = 0; i < m_arrayX.size(); ++i) { + // G4cout << m_arrayX[i] << Gateendl; + // } - // return the linear interpolation + //G4cout << "m_arrayY:" << Gateendl; + //for (size_t i = 0; i < m_arrayY.size(); ++i) { + // G4cout << m_arrayY[i] << Gateendl; + // } + + //return the linear interpolation return (y1 + (x-x1)*(y2-y1)/(x2-x1))/m_factorY; } -//___________________________________________________________________ + +//__________________________________________________________________ G4double GateVDistributionArray::RepartitionValue(G4double x) const { if ( m_arrayX.empty() ) return 0; @@ -119,24 +132,26 @@ G4double GateVDistributionArray::ShootRandom() const return (x1 + (y-y1)*(x2-x1)/(y2-y1))/m_factorX; } //___________________________________________________________________ -void GateVDistributionArray::FillRepartition() -{ - m_arrayRepartition.clear(); - if (m_arrayX.size()<2) return; - m_arrayRepartition.resize(m_arrayX.size()); - for (G4int i=0;i<(G4int)m_arrayX.size();++i){m_arrayRepartition[i]=0;} - for (G4int i=1;i<(G4int)m_arrayX.size();++i){ - G4double x1 = m_arrayX[i-1]; - G4double x2 = m_arrayX[i]; - G4double y1 = m_arrayY[i-1]; - G4double y2 = m_arrayY[i]; - m_arrayRepartition[i] = m_arrayRepartition[i-1] + (y1+y2)*(x2-x1)/2; - } - for (G4int i=0;i<(G4int)m_arrayX.size();++i){ - m_arrayRepartition[i]/=Integral(); -// G4cout<<"Repartition["<& array) const @@ -156,8 +171,8 @@ void GateVDistributionArray::InsertPoint(G4double x,G4double y) m_arrayY.push_back(y); } else if ( (i>=0) && m_arrayX[i] == x ) { //element replacement m_arrayY[i] = y; - G4cerr<<"[GateDistributionArray::InsertPoint] WARNING : replacing value for " - <m_maxY) m_maxY=y; } +void GateVDistributionArray::InsertPoint(G4double x, G4double y, G4double stddev) { + std::pair coordinates = std::make_pair(x, y); + stddevMap[coordinates] = stddev; + + + if (x < m_minX) m_minX = x; + if (x > m_maxX) m_maxX = x; + if (y < m_minY) m_minY = y; + if (y > m_maxY) m_maxY = y; + //G4cout << "Inserted: (x=" << x << ", y=" << y << ") -> stddev=" << stddev << std::endl; + + } + //___________________________________________________________________ void GateVDistributionArray::InsertPoint(G4double y) { if (GetSize()==0) InsertPoint(m_autoStart,y); else InsertPoint(m_maxX+1,y); } + + +G4double GateVDistributionArray::Value2D(double x, double y) const { + auto it = stddevMap.find(std::make_pair(x, y)); + if (it != stddevMap.end()) { + // If the point is found, return the standard deviation + G4double stddev= it->second; + // G4cout << "stddev for x=" << x << " and y=" << y << " is " << stddev<< G4endl; + return stddev; + } else { + //G4cout << "stddev not found for x=" << x << " and y=" << y << ". Performing bilinear interpolation." << G4endl; + G4double stddev = BilinearInterpolation(x, y); + return stddev; + } +} +G4double GateVDistributionArray::BilinearInterpolation(G4double x, G4double y) const { + auto lower_x_it = stddevMap.lower_bound(std::make_pair(x, -INFINITY)); + auto upper_x_it = stddevMap.upper_bound(std::make_pair(x, INFINITY)); + + if (lower_x_it == stddevMap.end() || upper_x_it == stddevMap.end()) { + G4cerr << "Error: Unable to find adjacent points in x direction." << G4endl; + return 0.0; + } + + // Adjust to make sure lower_x and upper_x are correct + G4double lower_x = (lower_x_it == stddevMap.begin()) ? lower_x_it->first.first : std::prev(lower_x_it)->first.first; + G4double upper_x = upper_x_it->first.first; + + if (lower_x == upper_x) { + G4cerr << "Error: Lower and upper x values are the same." << G4endl; + return 0.0; + } + + // Extract unique y values for these x coordinates + std::set y_values; + for (const auto& entry : stddevMap) { + if (entry.first.first == lower_x || entry.first.first == upper_x) { + y_values.insert(entry.first.second); + } + } + + if (y_values.empty()) { + G4cerr << "Error: No y values found for the given x values." << G4endl; + return 0.0; + } + + auto lower_y_it = y_values.lower_bound(y); + auto upper_y_it = y_values.upper_bound(y); + + if (lower_y_it == y_values.end() || upper_y_it == y_values.end()) { + G4cerr << "Error: Unable to find adjacent points in y direction." << G4endl; + return 0.0; + } + + G4double lower_y = (lower_y_it == y_values.begin()) ? *lower_y_it : *std::prev(lower_y_it); + G4double upper_y = *upper_y_it; + + if (lower_y == upper_y) { + G4cerr << "Error: Lower and upper y values are the same." << G4endl; + return 0.0; + } + + // Retrieve the sigma values for the four surrounding points + G4double stddev11 = stddevMap.at(std::make_pair(lower_x, lower_y)); + G4double stddev12 = stddevMap.at(std::make_pair(lower_x, upper_y)); + G4double stddev21 = stddevMap.at(std::make_pair(upper_x, lower_y)); + G4double stddev22 = stddevMap.at(std::make_pair(upper_x, upper_y)); + + //G4cout << "Interpolation points: (" << lower_x << ", " << lower_y << "), (" << upper_x << ", " << lower_y << "), (" + // << lower_x << ", " << upper_y << "), (" << upper_x << ", " << upper_y << ")" << G4endl; + //G4cout << "stddev values: " << stddev11 << ", " << stddev21 << ", " << stddev12 << ", " << stddev22 << G4endl; + + G4double interpolatedValue = (stddev11 * (upper_x - x) * (upper_y - y) + + stddev21 * (x - lower_x) * (upper_y - y) + + stddev12 * (upper_x - x) * (y - lower_y) + + stddev22 * (x - lower_x) * (y - lower_y)) / + ((upper_x - lower_x) * (upper_y - lower_y)); + + //G4cout << "Interpolated value: " << interpolatedValue << G4endl; + return interpolatedValue; + } + + From b244aa981c49186d160c2c18d3f7ec3c4169e168 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Thu, 4 Jul 2024 11:26:13 +0200 Subject: [PATCH 14/22] After verification, I discovered that I accidentally deleted confineInsideOfSmallestElement. --- source/digits_hits/src/GateSpatialResolutionMessenger.cc | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/source/digits_hits/src/GateSpatialResolutionMessenger.cc b/source/digits_hits/src/GateSpatialResolutionMessenger.cc index 6ef7b576d..0f3dd7f69 100755 --- a/source/digits_hits/src/GateSpatialResolutionMessenger.cc +++ b/source/digits_hits/src/GateSpatialResolutionMessenger.cc @@ -52,11 +52,10 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol cmdName = GetDirectoryName() + "fwhmXYdistrib2D"; spresolutionXYdistrib2DCmd = new G4UIcmdWithAString(cmdName,this); spresolutionXYdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); - cmdName = GetDirectoryName() + "fwhmYdistrib2D"; - - confineCmd = new G4UIcmdWithABool(cmdName,this); - confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); + cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; + confineCmd = new G4UIcmdWithABool(cmdName,this); + confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); } From 01ba9c51ba661e5389fac1ff6cab557e940a3501 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Fri, 5 Jul 2024 09:18:45 +0200 Subject: [PATCH 15/22] add 2D file verification in ReadMatrix2D --- .../digits_hits/src/GateDistributionFile.cc | 125 +++++++++++------- 1 file changed, 75 insertions(+), 50 deletions(-) diff --git a/source/digits_hits/src/GateDistributionFile.cc b/source/digits_hits/src/GateDistributionFile.cc index 93602c3cd..f0dd6312a 100644 --- a/source/digits_hits/src/GateDistributionFile.cc +++ b/source/digits_hits/src/GateDistributionFile.cc @@ -41,7 +41,7 @@ void GateDistributionFile::Read() { Clear(); G4cout<<"OPENING FILE "< xValues; - bool isFirstLine = true; - - while (std::getline(f, line)) { - std::stringstream iss(line); - - if (isFirstLine) { - // Process x values from the first line - G4double x; - - while (iss >> x) { - xValues.push_back(x); - if (iss.peek() == ',') iss.ignore(); - } - isFirstLine = false; - } else { - // Process y and stddev values for subsequent lines - G4double y; - iss >> y; - if (iss.peek() == ',') iss.ignore(); - // Read stddev values for each x value - for (size_t i = 0; i < xValues.size(); ++i) { - G4double stddev; - if (iss >> stddev) { - // Insert the (x, y) -> stddev pair into the map using InsertPoint - InsertPoint(xValues[i], y, stddev); - if (iss.peek() == ',') iss.ignore(); - } - } - } - } - - f.close(); - //G4cout << "Content of stddevMap:\n"; - //for (const auto& entry : stddevMap) { - // std::pair coordinates = entry.first; - // G4double stddev = entry.second; - //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl; - //} +void GateDistributionFile::ReadMatrix2d() { + Clear(); + G4cout << "OPENING FILE " << m_FileName << Gateendl; + G4cout << "This must be a 2D distribution" << Gateendl; + std::ifstream f(m_FileName); + if (!f.is_open()) { + G4cerr << "[GateDistributionFile::ReadMatrix2d] WARNING: File " << m_FileName << " can't be opened\n"; + return; + } + std::string line; + std::vector xValues; + bool isFirstLine = true; + bool is2D = true; + + while (std::getline(f, line)) { + std::stringstream iss(line); + + if (isFirstLine) { + // Process x values from the first line + G4double x; + size_t xCount = 0; + + while (iss >> x) { + xValues.push_back(x); + ++xCount; + if (iss.peek() == ',') iss.ignore(); + } + + // Check if we have at least two x values to be considered 2D + if (xCount < 3) { + is2D = false; + break; + } + + isFirstLine = false; + } else { + // Process y and stddev values for subsequent lines + G4double y; + iss >> y; + if (iss.peek() == ',') iss.ignore(); + size_t stddevCount = 0; + + // Read stddev values for each x value + for (size_t i = 0; i < xValues.size(); ++i) { + G4double stddev; + if (iss >> stddev) { + ++stddevCount; + // Insert the (x, y) -> stddev pair into the map using InsertPoint + InsertPoint(xValues[i], y, stddev); + if (iss.peek() == ',') iss.ignore(); + } + } + + // Check if the number of stddev values matches the number of x values + if (stddevCount != xValues.size()) { + is2D = false; + break; + } + } + } + + f.close(); + //G4cout << "Content of stddevMap:\n"; + //for (const auto& entry : stddevMap) { + // std::pair coordinates = entry.first; + // G4double stddev = entry.second; + //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl; + //} + + if (!is2D) { + G4cerr << "[GateDistributionFile::ReadMatrix2d] ERROR: File " << m_FileName << " is not a valid 2D distribution file\n"; + Clear(); + return; + } FillRepartition(); } From 21adda75cdcaee87f34b822236bd5eb8bfae46a5 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Tue, 9 Jul 2024 09:10:24 +0200 Subject: [PATCH 16/22] added SetSpatialResolutionParameters to check the parameter --- .../include/GateSpatialResolution.hh | 4 +- .../digits_hits/src/GateSpatialResolution.cc | 50 +++++++++---------- 2 files changed, 27 insertions(+), 27 deletions(-) mode change 100755 => 100644 source/digits_hits/include/GateSpatialResolution.hh mode change 100755 => 100644 source/digits_hits/src/GateSpatialResolution.cc diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh old mode 100755 new mode 100644 index c28e7b37a..1c27448e9 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/source/digits_hits/include/GateSpatialResolution.hh @@ -67,7 +67,7 @@ public: void SetFWHMx(G4double val) { m_fwhmX = val; } void SetFWHMy(G4double val) { m_fwhmY = val; } void SetFWHMz(G4double val) { m_fwhmZ = val; } - + void SetSpatialResolutionParameters(); inline void ConfineInsideOfSmallestElement(const G4bool& value) { m_IsConfined = value; }; inline G4bool IsConfinedInsideOfSmallestElement() const { return m_IsConfined; } @@ -109,7 +109,7 @@ private: GateDigiCollection* m_OutputDigiCollection; GateSinglesDigitizer *m_digitizer; - + G4bool m_IsFirstEntrance; G4VoxelLimits limits; G4double Xmin, Xmax, Ymin, Ymax, Zmin, Zmax; G4AffineTransform at; diff --git a/source/digits_hits/src/GateSpatialResolution.cc b/source/digits_hits/src/GateSpatialResolution.cc old mode 100755 new mode 100644 index 7697f532e..1621a0318 --- a/source/digits_hits/src/GateSpatialResolution.cc +++ b/source/digits_hits/src/GateSpatialResolution.cc @@ -62,6 +62,7 @@ GateSpatialResolution::GateSpatialResolution(GateSinglesDigitizer *digitizer, G4 m_Touchable(0), m_systemDepth(-1), m_outputDigi(0), + m_IsFirstEntrance(1), m_OutputDigiCollection(0), m_digitizer(digitizer) { @@ -76,32 +77,33 @@ GateSpatialResolution::~GateSpatialResolution() delete m_Messenger; } +void GateSpatialResolution::SetSpatialResolutionParameters() { + // Check FWHM parameters + if (m_fwhm != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmZ != 0 || m_fwhmXYdistrib2D != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 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_fwhmXYdistrib2D != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM 2D distribution for (X,Y) OR set FWHM for X, Y individually." << G4endl; + abort(); + } -void GateSpatialResolution::Digitize(){ - -if (m_fwhm != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmZ != 0 || m_fwhmXYdistrib2D != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 0)) -{ - GateError("***ERROR*** Spatial Resolution is ambiguous: you can set a unique FWHM for all 3 axes OR set FWHM for X, Y, Z individually."); -} - - -if (m_fwhmXYdistrib2D != 0 && (m_fwhmX != 0 || m_fwhmY != 0 || m_fwhmXdistrib != 0 || m_fwhmYdistrib != 0)) -{ - GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM 2D distribution for (X,Y) OR set FWHM for X ,Y individually ."); -} - + if (m_fwhmY != 0 && (m_fwhmYdistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM for Y distribution." << G4endl; + abort(); + } -if (m_fwhmY != 0 && (m_fwhmYdistrib != 0)) -{ - GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for Y OR set FWHM for Y distribution."); -} + if (m_fwhmX != 0 && (m_fwhmXdistrib != 0)) { + G4cout << "***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM for X distribution." << G4endl; + abort(); + }} -// Check if the spatial resolutions FWHM for the X axis are defined ambiguously -if (m_fwhmX != 0 && (m_fwhmXdistrib != 0)) -{ - GateError("***ERROR*** Spatial Resolution is ambiguous: you can set FWHM for X OR set FWHM for X distribution."); -} +void GateSpatialResolution::Digitize(){ + if (m_IsFirstEntrance) { + SetSpatialResolutionParameters(); + m_IsFirstEntrance = false; + } G4double fwhmX; @@ -222,7 +224,7 @@ if (m_fwhmX != 0 && (m_fwhmXdistrib != 0)) G4double PxNew = G4RandGauss::shoot(Px,stddevX); G4double PyNew = G4RandGauss::shoot(Py,stddevY); G4double PzNew = G4RandGauss::shoot(Pz,fwhmZ/GateConstants::fwhm_to_sigma); - if (m_IsConfined) + if (m_IsConfined) { //set the position on the border of the crystal //no need to update volume ID @@ -309,8 +311,6 @@ void GateSpatialResolution::UpdateVolumeID() G4int CopyNo = m_Touchable->GetReplicaNumber(m_systemDepth-1-i); m_outputDigi->ChangeVolumeIDAndOutputVolumeIDValue(i,CopyNo); } - - } void GateSpatialResolution::DescribeMyself(size_t indent ) From 2695629c11eaf8a3bfd3f33d12aad0107cb82582 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Thu, 11 Jul 2024 15:47:42 +0200 Subject: [PATCH 17/22] added units for spatial resolution --- .../include/GateSpatialResolutionMessenger.hh | 9 ++++--- .../src/GateSpatialResolutionMessenger.cc | 24 ++++++++++--------- 2 files changed, 17 insertions(+), 16 deletions(-) mode change 100755 => 100644 source/digits_hits/include/GateSpatialResolutionMessenger.hh mode change 100755 => 100644 source/digits_hits/src/GateSpatialResolutionMessenger.cc diff --git a/source/digits_hits/include/GateSpatialResolutionMessenger.hh b/source/digits_hits/include/GateSpatialResolutionMessenger.hh old mode 100755 new mode 100644 index a140b97cd..7d3a5ebd7 --- a/source/digits_hits/include/GateSpatialResolutionMessenger.hh +++ b/source/digits_hits/include/GateSpatialResolutionMessenger.hh @@ -40,14 +40,13 @@ public: private: GateSpatialResolution* m_SpatialResolution; - G4UIcmdWithADouble* spresolutionCmd; - G4UIcmdWithADouble* spresolutionXCmd; + G4UIcmdWithADoubleAndUnit* spresolutionCmd; + G4UIcmdWithADoubleAndUnit* spresolutionXCmd; + G4UIcmdWithADoubleAndUnit* spresolutionYCmd; + G4UIcmdWithADoubleAndUnit* spresolutionZCmd; G4UIcmdWithAString *spresolutionXdistribCmd;// Command declaration for 1D X-resolution distribution G4UIcmdWithAString *spresolutionYdistribCmd;// Command declaration for 1D Y-resolution distribution G4UIcmdWithAString *spresolutionXYdistrib2DCmd; // Command declaration for 2D XY-resolution distribution - - G4UIcmdWithADouble* spresolutionYCmd; - G4UIcmdWithADouble* spresolutionZCmd; G4UIcmdWithABool* confineCmd; diff --git a/source/digits_hits/src/GateSpatialResolutionMessenger.cc b/source/digits_hits/src/GateSpatialResolutionMessenger.cc old mode 100755 new mode 100644 index 0f3dd7f69..24ab69f5d --- a/source/digits_hits/src/GateSpatialResolutionMessenger.cc +++ b/source/digits_hits/src/GateSpatialResolutionMessenger.cc @@ -11,10 +11,9 @@ See LICENSE.md for further details #include "GateSpatialResolutionMessenger.hh" #include "GateSpatialResolution.hh" #include "GateDigitizerMgr.hh" - #include "G4SystemOfUnits.hh" #include "G4UIcmdWithADouble.hh" - +#include "G4UIcmdWithADoubleAndUnit.hh" #include "G4UIcmdWithABool.hh" #include "G4UIcmdWithAString.hh" #include "G4UIdirectory.hh" @@ -31,18 +30,21 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol G4String cmdName; cmdName = GetDirectoryName() + "fwhm"; - spresolutionCmd = new G4UIcmdWithADouble(cmdName,this); - spresolutionCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); + spresolutionCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this); + spresolutionCmd->SetGuidance("Set the resolution (in mm) in position for gaussian spblurring"); + spresolutionCmd->SetUnitCategory("Length"); cmdName = GetDirectoryName() + "fwhmX"; - spresolutionXCmd= new G4UIcmdWithADouble(cmdName,this); - spresolutionXCmd->SetGuidance("Set the resolution "); + spresolutionXCmd= new G4UIcmdWithADoubleAndUnit(cmdName,this); + spresolutionXCmd->SetGuidance("Set the resolution (in mm) in position for gaussian spblurring"); + spresolutionXCmd->SetUnitCategory("Length"); cmdName = GetDirectoryName() + "fwhmY"; - spresolutionYCmd = new G4UIcmdWithADouble(cmdName,this); + spresolutionYCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this); spresolutionYCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); - + spresolutionYCmd->SetUnitCategory("Length"); cmdName = GetDirectoryName() + "fwhmZ"; - spresolutionZCmd = new G4UIcmdWithADouble(cmdName,this); + spresolutionZCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this); spresolutionZCmd->SetGuidance("Set the resolution in position for gaussian spblurring"); + spresolutionZCmd->SetUnitCategory("Length"); cmdName = GetDirectoryName() + "fwhmXdistrib"; spresolutionXdistribCmd = new G4UIcmdWithAString(cmdName,this); spresolutionXdistribCmd->SetGuidance("Set the distribution resolution in position for gaussian spblurring"); @@ -54,8 +56,8 @@ GateSpatialResolutionMessenger::GateSpatialResolutionMessenger (GateSpatialResol spresolutionXYdistrib2DCmd->SetGuidance("Set the distribution 2D of spatial resolution in position for gaussian spblurring"); cmdName = GetDirectoryName() + "confineInsideOfSmallestElement"; - confineCmd = new G4UIcmdWithABool(cmdName,this); - confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); + confineCmd = new G4UIcmdWithABool(cmdName,this); + confineCmd->SetGuidance("To be set true, if you want to moves the outsiders of the crystal after spblurring inside the same crystal"); } From 24d61149e242f9117f50207a49fcda4bc92fe1af Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Mon, 29 Jul 2024 10:04:33 +0200 Subject: [PATCH 18/22] Update digitizer_and_detector_modeling.rst Documentation for Multiples Killer --- docs/digitizer_and_detector_modeling.rst | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 73f611bfd..19022836a 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -1422,9 +1422,9 @@ A presort buffer contains singles that have not yet been checked for coincidence Multiple coincidence removal ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -If the multiple coincidences are kept and not splitted into pairs (ie. if any of the **keepXXX** multiple coincidence policy is used), the multicoincidences could participate to dataflow occupancy, but could not be written to the disk. Unless otherwise specified, any multicoincidence is then cleared from data just before the disk writing. If needed, this clearing could be performed at any former coincidence processing step, by inserting the **multipleKiller** module at the required level. This module has no parameter and just kill the multicoincidence events. Multiple coincidences split into many pairs are not affected by this module and cannot be distinguished from the normal "simple" coincidences. To insert a multipleKiller, one has to use the syntax:: +If the multiple coincidences are kept and not split into pairs (i.e., if any of the **keepXXX** multiple coincidence policies are used), the multicoincidences could contribute to dataflow occupancy but cannot be written to the disk. Unless otherwise specified, any multicoincidence is then cleared from data just before the disk writing. If needed, this clearing could be performed at any earlier coincidence processing step by inserting the **multipleKiller** module at the required level. This module has no parameters and simply removes the multicoincidence events. Multiple coincidences split into many pairs are not affected by this module and cannot be distinguished from normal "simple" coincidences. To insert a multipleKiller, use the syntax:: +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert MultiplesKiller - /gate/digitizer/myCoincChain/insert multipleKiller Example of a digitizer setting ------------------------------ @@ -1529,6 +1529,7 @@ Example:: 87 /gate/digitizer/finalCoinc/buffer/setBufferSize 32 B 88 /gate/digitizer/finalCoinc/buffer/setReadFrequency 14.45 MHz 89 /gate/digitizer/finalCoinc/buffer/setMode 0 + Lines 1 to 15: The branch named "Singles" contains the result of applying the adder, readout, blurring, and threshold (50 keV) modules. From 6b054c4cc9db68f17abbafea769ac332f1fd921c Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Wed, 31 Jul 2024 09:41:16 +0200 Subject: [PATCH 19/22] Update digitizer_and_detector_modeling.rst --- docs/digitizer_and_detector_modeling.rst | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 19022836a..e80f9a04c 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -1154,7 +1154,7 @@ Then, the rejection can be set to the whole event or only to those pulses within Example:: - /gate/digitizerMgr/absorber/SinglesDigitizer/Singles/insert multipleRejection + /gate/digitizerMgr/absorber/SinglesDigitizer/Singles/insert multipleRejection /gate/digitizerMgr/absorber/SinglesDigitizer/Singles/multipleRejection/setMultipleDefinition volumeID /gate/digitizerMgr/absorber/SinglesDigitizer/Singles/multipleRejection/setEventRejection 1 @@ -1411,11 +1411,21 @@ The dead time for coincidences works in the same way as that acting on the *sing Coincidence buffers ~~~~~~~~~~~~~~~~~~~ +It simulates the operation of a detector by modeling coincidences, transfer speed limits, and data loss due to buffer capacity overflows. It manages a memory buffer for coincidence events, allowing for the modeling of data loss due to buffer overflow. It uses a read frequency and allows defining the buffer size and read mode, influencing how events are processed. There are two buffer operation modes. Mode 1 empties the entire buffer at each read clock tick, while Mode 0 reads events one by one. + +Example:: + +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert Buffer +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setBufferSize 64 B +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setReadFrequency 10 MHz +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setMode 1 + For a coincidence sorter user can chose a presort buffer with a following command: /gate/digitizer/Coincidences/setPresortBufferSize 256 + A presort buffer contains singles that have not yet been checked for coincidence with the already open coincidence windows. The default value is 256, the minimum value is 32. For more details, check https://iopscience.iop.org/article/10.1088/0031-9155/61/18/N522 @@ -1423,6 +1433,8 @@ Multiple coincidence removal ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If the multiple coincidences are kept and not split into pairs (i.e., if any of the **keepXXX** multiple coincidence policies are used), the multicoincidences could contribute to dataflow occupancy but cannot be written to the disk. Unless otherwise specified, any multicoincidence is then cleared from data just before the disk writing. If needed, this clearing could be performed at any earlier coincidence processing step by inserting the **multipleKiller** module at the required level. This module has no parameters and simply removes the multicoincidence events. Multiple coincidences split into many pairs are not affected by this module and cannot be distinguished from normal "simple" coincidences. To insert a multipleKiller, use the syntax:: + + /gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert MultiplesKiller From 70ddb9750038f7997382571aae1b331b8bb83994 Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Thu, 29 Aug 2024 15:57:29 +0200 Subject: [PATCH 20/22] modification --- source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc b/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc index b4b2210f6..2c8359942 100644 --- a/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc +++ b/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc @@ -124,7 +124,7 @@ void GateCoincidenceDigitizerMessenger::SetNewValue(G4UIcommand* command,G4Strin const G4String& GateCoincidenceDigitizerMessenger::DumpMap() { - static G4String theList = "deadtime MultiplesKiller";//readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger"; + static G4String theList = "deadtime multiplesKiller";//readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger"; return theList; } @@ -151,7 +151,7 @@ void GateCoincidenceDigitizerMessenger::DoInsertion(const G4String& childTypeNam m_CoinDigitizer->AddNewModule(newDM); } - else if (childTypeName=="MultiplesKiller") + else if (childTypeName=="multiplesKiller") { newDM = new GateCoincidenceMultiplesKiller(m_CoinDigitizer, DMname); m_CoinDigitizer->AddNewModule(newDM); From 88f543ea0720b22000efaa1146d98604ca2fc6fe Mon Sep 17 00:00:00 2001 From: Radia Oudihat Date: Fri, 30 Aug 2024 09:45:38 +0200 Subject: [PATCH 21/22] modification --- source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc b/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc index e8c151462..c6dd4c040 100644 --- a/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc +++ b/source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc @@ -123,7 +123,7 @@ void GateCoincidenceDigitizerMessenger::SetNewValue(G4UIcommand* command,G4Strin const G4String& GateCoincidenceDigitizerMessenger::DumpMap() { - static G4String theList = "deadtime MultiplesKiller Buffer";//readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger"; + static G4String theList = "deadtime multiplesKiller buffer";//readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger"; return theList; } @@ -155,7 +155,7 @@ void GateCoincidenceDigitizerMessenger::DoInsertion(const G4String& childTypeNam m_CoinDigitizer->AddNewModule(newDM); } - else if (childTypeName=="Buffer") + else if (childTypeName=="buffer") { newDM = new GateCoincidenceBuffer(m_CoinDigitizer, DMname); m_CoinDigitizer->AddNewModule(newDM); From 961ce782fffdffd562032c5e350a9e2a1c8dea6b Mon Sep 17 00:00:00 2001 From: Oudihat-Radia Date: Fri, 30 Aug 2024 10:51:07 +0200 Subject: [PATCH 22/22] Update digitizer_and_detector_modeling.rst --- docs/digitizer_and_detector_modeling.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index e80f9a04c..43c238524 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -1415,10 +1415,10 @@ It simulates the operation of a detector by modeling coincidences, transfer spee Example:: -/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert Buffer -/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setBufferSize 64 B -/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setReadFrequency 10 MHz -/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/Buffer/setMode 1 +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert buffer +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/buffer/setBufferSize 64 B +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/buffer/setReadFrequency 10 MHz +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/buffer/setMode 1 For a coincidence sorter user can chose a presort buffer with a following command: @@ -1435,7 +1435,7 @@ Multiple coincidence removal If the multiple coincidences are kept and not split into pairs (i.e., if any of the **keepXXX** multiple coincidence policies are used), the multicoincidences could contribute to dataflow occupancy but cannot be written to the disk. Unless otherwise specified, any multicoincidence is then cleared from data just before the disk writing. If needed, this clearing could be performed at any earlier coincidence processing step by inserting the **multipleKiller** module at the required level. This module has no parameters and simply removes the multicoincidence events. Multiple coincidences split into many pairs are not affected by this module and cannot be distinguished from normal "simple" coincidences. To insert a multipleKiller, use the syntax:: -/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert MultiplesKiller +/gate/digitizerMgr/CoincidenceDigitizer/finalCoinc/insert multiplesKiller Example of a digitizer setting