diff --git a/docs/digitizer_and_detector_modeling.rst b/docs/digitizer_and_detector_modeling.rst index 2bcae1a4f..e847203f1 100644 --- a/docs/digitizer_and_detector_modeling.rst +++ b/docs/digitizer_and_detector_modeling.rst @@ -1159,6 +1159,107 @@ Example:: /gate/digitizerMgr/absorber/SinglesDigitizer/Singles/multipleRejection/setEventRejection 1 + + + + + + +Virtual segmentation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + +In traditional PET image reconstruction, software like CASToR utilizes crystal IDs instead of the position of the interaction. This approach has sufficed because the small crystals' size inherently defines the spatial resolution. However, new PET scanner systems are exploring the use of monolithic crystals, which can reconstruct interaction positions within the crystal with a specific resolution. + +The Virtual Segmentation Digitizer module provides a mechanism to generate an ID based on a virtual segmentation of the monolithic crystal, aligned with its spatial resolution. Ideally, the pitch size should be at least half of the position resolution. This virtual segmentation occurs post-simulation, ensuring that the simulation speed remains uncompromised even when dealing with large systems and numerous crystals such as total-body PET scans. + +A GateTool associated with this digitizer allows users to create a new geometry macro with the segmented geometry, suitable for use in image reconstruction software. + +**Geometry Requirements:** + +To utilize this digitizer, the cylindricalPET geometry must be configured as follows: +The size of the crystal should be defined at the Submodule level using Air as the material. +The crystal and layer0 levels should both reflect the crystal size and use the crystal's material. +This setup allows new virtual IDs for the XYZ axes to be assigned at the Layer, Crystal, and Submodule levels, respectively. + +Example:: + + # CRYSTAL + /gate/rsector/daughters/name crystal + /gate/rsector/daughters/insert box + /gate/crystal/geometry/setXLength 10. mm + /gate/crystal/geometry/setYLength 59. mm + /gate/crystal/geometry/setZLength 59. mm + /gate/crystal/setMaterial Air + + # Level saved for the virtual COLUMN + /gate/crystal/daughters/name column + /gate/crystal/daughters/insert box + /gate/column/geometry/setXLength 10. mm + /gate/column/geometry/setYLength 59. mm + /gate/column/geometry/setZLength 59. mm + /gate/column/setMaterial Air + + + #Level saved for the virtual ROW + /gate/column/daughters/name row + /gate/column/daughters/insert box + /gate/row/geometry/setXLength 10. mm + /gate/row/geometry/setYLength 59. mm + /gate/row/geometry/setZLength 59. mm + /gate/row/setMaterial Air + + #Level saved for the virtual Pseudo-Crystal (now the real material needs to be used) + /gate/row/daughters/name pseudo-crystal + /gate/row/daughters/insert box + /gate/pseudo-crystal/geometry/setXLength 10. mm + /gate/pseudo-crystal/geometry/setYLength 59. mm + /gate/pseudo-crystal/geometry/setZLength 59. mm + /gate/pseudo-crystal/setMaterial PWO + + +**Commands** + +*"nameAxis"* Enable users to specify which axes require discretization. The axis selected can be any combination of "XYZ", "XY", "XZ" or "YZ", "X", "Y" or "Z". + +*"pitch, pitchX, pitchY, pitchZ"* allow users to specify the desired pitch size for all axis, or specific values for X,Y and Z axis. If no pitch size is provided, the digitizer will use the spatial resolution value to compute the optimal pitch size ensuring that the number of bins, defined as (crystal size)/(pitch), is integer. Note that the spatial resolution must be a single value for each axis; if a the spatial resolution is defined with a distribution and no pitch value is provided, the digitizer will not function correctly. + +*"useMacroGenerator"* boolean flag that determines if the user wants to generate the new geometry macro with the segmentation. + + + + +Example providing spatial resolution:: + +/gate/digitizerMgr//SinglesDigitizer/Singles/insert spatialResolution +/gate/digitizerMgr//SinglesDigitizer/Singles/spatialResolution/fwhm 2. mm + +/gate/digitizerMgr//SinglesDigitizer/Singles/insert virtualSegmentation +/gate/digitizerMgr//SinglesDigitizer/Singles/virtualSegmentation/nameAxis XYZ +/gate/digitizerMgr//SinglesDigitizer/Singles/virtualSegmentation/useMacroGenerator true + +In this case, a value for the FWHM for the spatial resolution was provided but no pitch size was given to the virtual segmentation. The module will read the value of spatial resolution and will generate a pitch size that is, at least, half of the value of the FWHM in spatial resolution. + + +Example providing the pitch:: + +/gate/digitizerMgr//SinglesDigitizer/Singles/insert virtualSegmentation +/gate/digitizerMgr//SinglesDigitizer/Singles/virtualSegmentation/nameAxis XYZ +/gate/digitizerMgr//SinglesDigitizer/Singles/virtualSegmentation/pitch 1.0 mm +/gate/digitizerMgr//SinglesDigitizer/Singles/virtualSegmentation/useMacroGenerator true + +In this case, the pitch size is provided by the user and it will be used regardless of any values of FWHM provided in the spatial resolution module. + + + + + + + + + + + .. _digitizer_multiple_processor_chains-label: Examples of multiple Single Digitizers @@ -1766,3 +1867,4 @@ Finally, we call the 'triCoincProcessor' module and we plug on it the second sys /gate/digitizer/TriCoinc/triCoincProcessor/setWindow 15 ns /gate/digitizer/TriCoinc/triCoincProcessor/setSinglesBufferSize 40 + diff --git a/source/digits_hits/include/GateDigi.hh b/source/digits_hits/include/GateDigi.hh index 7be01d0ae..cf513cb44 100644 --- a/source/digits_hits/include/GateDigi.hh +++ b/source/digits_hits/include/GateDigi.hh @@ -125,6 +125,7 @@ public: inline G4double GetScannerRotAngle() const { return m_scannerRotAngle; } inline void SetOutputVolumeID(const GateOutputVolumeID& outputVolumeID) { m_outputVolumeID = outputVolumeID; } + inline void SetOutputVolumeID(const G4int copyNum,G4int depth) { m_outputVolumeID[depth] = copyNum; } inline const GateOutputVolumeID& GetOutputVolumeID() const { return m_outputVolumeID; } inline G4int GetComponentID(size_t depth) const { return (m_outputVolumeID.size()>depth) ? m_outputVolumeID[depth] : -1; } diff --git a/source/digits_hits/include/GateDiscretizer.hh b/source/digits_hits/include/GateDiscretizer.hh new file mode 100644 index 000000000..d6f807b8a --- /dev/null +++ b/source/digits_hits/include/GateDiscretizer.hh @@ -0,0 +1,53 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +// OK GND 2022 +/*! + \class GateDiscretizer (by marc.granado@universite-paris-saclay.fr) + \brief for discretizing the position within a monolithic crystal + + - For each volume the local position of the interactions within the crystal are discretized. + the X,Y,Z vector is translated into the IDs of virtual divisions within the crystal + ocuppying the levels of SubmoduleID, CrystalID and LayerID. + +*/ + +#ifndef GateDiscretizer_h +#define GateDiscretizer_h 1 + +#include "GateVDigitizerModule.hh" +#include "GateDigi.hh" +#include "GateClockDependent.hh" +#include "GateCrystalSD.hh" + +#include "globals.hh" + +#include "GateSinglesDigitizer.hh" + + +class GateDiscretizer : public GateVDigitizerModule +{ +public: + + GateDiscretizer(GateSinglesDigitizer *digitizer, G4String name); + ~GateDiscretizer(); + adder_policy_t m_positionPolicy; + +private: + GateAdderMessenger *m_Messenger; + + GateDigi* m_outputDigi; + GateDigiCollection* m_OutputDigiCollection; + GateSinglesDigitizer *m_digitizer; + + + + +}; + +#endif diff --git a/source/digits_hits/include/GateSpatialResolution.hh b/source/digits_hits/include/GateSpatialResolution.hh index 1c27448e9..778ca3afa 100644 --- a/source/digits_hits/include/GateSpatialResolution.hh +++ b/source/digits_hits/include/GateSpatialResolution.hh @@ -46,14 +46,14 @@ public: //! These functions return the resolution in use. - G4double GetFWHM() { return m_fwhm; } - GateVDistribution* GetFWHMxdistrib() { return m_fwhmXdistrib; } - GateVDistribution* GetFWHMydistrib() { return m_fwhmYdistrib; } - GateVDistribution* GetFWHMxydistrib2D() { return m_fwhmXYdistrib2D; } + G4double GetFWHM() { return m_fwhm; } + 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; } + G4double GetFWHMx() { return m_fwhmX; } + G4double GetFWHMy() { return m_fwhmY; } + G4double GetFWHMz() { return m_fwhmZ; } //! These functions set the spresolution of a gaussian spblurring. /*! @@ -76,6 +76,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 ); diff --git a/source/digits_hits/include/GateVirtualSegmentationSD.hh b/source/digits_hits/include/GateVirtualSegmentationSD.hh new file mode 100644 index 000000000..005154d50 --- /dev/null +++ b/source/digits_hits/include/GateVirtualSegmentationSD.hh @@ -0,0 +1,135 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +// OK GND 2022 +/*This class is not used by GATE ! + The purpose of this class is to help to create new users digitizer module(DM). + Please, check GateVirtualSegmentationSD.cc for more detals + */ + + +/*! \class GateVirtualSegmentationSD + \brief GateVirtualSegmentationSD - For each volume the local position of the interactions within the crystal are virtually segmented. + the X,Y,Z vector is translated into the IDs of virtual divisions within the crystal + + + - GateVirtualSegmentationSD - by marc.granado@universite-paris-saclay.fr + + \sa GateVirtualSegmentationSD, GateVirtualSegmentationSDMessenger +*/ + +#ifndef GateVirtualSegmentationSD_h +#define GateVirtualSegmentationSD_h 1 + +#include "GateVDigitizerModule.hh" +#include "GateDigi.hh" +#include "GateClockDependent.hh" +#include "GateCrystalSD.hh" + +#include "globals.hh" + +#include "GateVirtualSegmentationSDMessenger.hh" +#include "GateSinglesDigitizer.hh" + + +class GateVirtualSegmentationSD : public GateVDigitizerModule +{ +public: + + GateVirtualSegmentationSD(GateSinglesDigitizer *digitizer, G4String name); + + ~GateVirtualSegmentationSD(); + + void Digitize() override; + + //! These functions return the pitch in use. + G4String GetNameAxis() { return m_nameAxis; } + G4double GetPitch() { return m_pitch; } + G4double GetPitchX() { return m_pitchX; } + G4double GetPitchY() { return m_pitchY; } + G4double GetPitchZ() { return m_pitchZ; } + + G4double calculatePitch(G4double crystal_size, G4double target_pitch); + + + //Set Parameter + void SetNameAxis(const G4String&); + + void SetUseMacroGenerator(G4bool val) {useMacroGenerator=val;} + + //Set Variables + void SetPitch(G4double val) { m_pitch = val; } + void SetPitchX(G4double val) { m_pitchX = val; } + void SetPitchY(G4double val) { m_pitchY = val; } + void SetPitchZ(G4double val) { m_pitchZ = val; } + + void SetParameters(); + + void SetVirtualID(int nBins,double pitch, G4double pos, int depth); + void DescribeMyself(size_t ); + +protected: + // *******implement your parameters here + G4double m_pitch; + + G4String m_nameAxis; + G4double m_pitchX; + G4double m_pitchY; + G4double m_pitchZ; + + G4int nBinsX; + G4int nBinsY; + G4int nBinsZ; + + G4double pitchX; + G4double pitchY; + G4double pitchZ; + + G4int depthX; + G4int depthY; + G4int depthZ; + + G4double xLength; + G4double yLength; + G4double zLength; + + //These are in Spatial pitch but aren't there rom touchable? + //We'll need to check + + //G4Navigator* m_Navigator; + //G4TouchableHistoryHandle m_Touchable; + +private: + + G4int m_systemDepth; + + G4int m_IsFirstEntry; + + G4bool useMacroGenerator; + + GateDigi* m_outputDigi; + + GateVirtualSegmentationSDMessenger *m_Messenger; + + GateDigiCollection* m_OutputDigiCollection; + + GateSinglesDigitizer *m_digitizer; + + +}; + + +#endif + + + + + + + + diff --git a/source/digits_hits/include/GateVirtualSegmentationSDMessenger.hh b/source/digits_hits/include/GateVirtualSegmentationSDMessenger.hh new file mode 100644 index 000000000..9a9fde1dc --- /dev/null +++ b/source/digits_hits/include/GateVirtualSegmentationSDMessenger.hh @@ -0,0 +1,65 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +// OK GND 2022 +/*This class is not used by GATE ! + The purpose of this class is to help to create new users digitizer module(DM). + Please, check GateVirtualSegmentationSD.cc for more detals + */ + + +/*! \class GateVirtualSegmentationSDMessenger + \brief Messenger for the GateVirtualSegmentationSD + + - GateVirtualSegmentationSD - by marc.granado@universite-paris-saclay.fr + + \sa GateVirtualSegmentationSD, GateVirtualSegmentationSDMessenger +*/ + + +#ifndef GateVirtualSegmentationSDMessenger_h +#define GateVirtualSegmentationSDMessenger_h 1 + +#include "G4UImessenger.hh" +#include "globals.hh" + +#include "GateClockDependentMessenger.hh" +class GateVirtualSegmentationSD; +class G4UIcmdWithAString; + +class GateVirtualSegmentationSDMessenger : public GateClockDependentMessenger +{ +public: + + GateVirtualSegmentationSDMessenger(GateVirtualSegmentationSD*); + ~GateVirtualSegmentationSDMessenger(); + + void SetNewValue(G4UIcommand*, G4String); + + +private: + GateVirtualSegmentationSD* m_VirtualSegmentationSD; + + G4UIcmdWithAString* nameAxisCmd; + + G4UIcmdWithABool* useMacroGeneratorCmd; + G4UIcmdWithADoubleAndUnit* pitchCmd; + G4UIcmdWithADoubleAndUnit* pitchXCmd; + G4UIcmdWithADoubleAndUnit* pitchYCmd; + G4UIcmdWithADoubleAndUnit* pitchZCmd; +}; + +#endif + + + + + + + + diff --git a/source/digits_hits/src/GateDigi.cc b/source/digits_hits/src/GateDigi.cc index e7462243a..d6bea9e2f 100644 --- a/source/digits_hits/src/GateDigi.cc +++ b/source/digits_hits/src/GateDigi.cc @@ -190,6 +190,9 @@ void GateDigi::ChangeVolumeIDAndOutputVolumeIDValue(size_t depth, G4int copyNo) delete volSelector; // Finally change the outputVolumeID accordingly m_outputVolumeID[depth] = copyNo; + + //With the new method the line above could be changed to: + //SetOutputVolumeID(copyNo,depth) } diff --git a/source/digits_hits/src/GateDummyDigitizerModule.cc b/source/digits_hits/src/GateDummyDigitizerModule.cc index a6a6ac6d5..8df34bf49 100755 --- a/source/digits_hits/src/GateDummyDigitizerModule.cc +++ b/source/digits_hits/src/GateDummyDigitizerModule.cc @@ -15,7 +15,7 @@ - Create your DM by coping this class and GateDummyDigitizerMessenger class for your DM messenger - Places to change are marked with // ****** comment and called with "dummy" names - - Include your module to GateSinglesDigitizerMessenger in the method DoInsertion(..) + - Include your module to GateSinglesDigitizerMessenger.cc in the method DoInsertion(..) If you adapting some already exiting class from Old Gate Digitizer here is some of the tips - Digitize () is a fusion of GateVPulseProcessor::ProcessPulseList and GateXXX::ProcessOnePulse @@ -28,7 +28,7 @@ To create new Digitizer Module (DM), please, follow the steps: 1) Copy .cc and .hh of GateDummyDigitizerModule, GateDummyDigitizerModuleMessenger to GateYourNewDigitizerModule and GateYourNewDigitizerModuleMessenger 2) Replace in these new files : DummyDigitizerModule -> YourNewDigitizerModule - 3) Compile 1st time (so not forget to redo ccmake to) + 3) Compile 1st time (so not forget to redo ccmake too) 4) Adapt GateYourNewDigitizerModuleMessenger.cc 5) Adapt GateYourNewDigitizerModuleMessenger.hh (!!!! DO NOT FORGET TO WRITE A SHORT EXPLANATION ON WHAT DOES YOUR DM !!!!) 6) Adapt GateYourNewDigitizerModule.hh @@ -46,7 +46,7 @@ inputPulse -> inputDigi outputPulse -> m_outputDigi + correct the first declaration (as in this Dummy module) outputPulseList.push_back(outputPulse) -> m_OutputDigiCollection->insert(m_outputDigi); - 10) Add YourDigitizerModule to GateSinglesDigitizer.cc + 10) Add YourDigitizerModule to GateSinglesDigitizerMessenger - #include "YourDigitizerModule.hh" - in DumpMap() method in static G4String theList = " ...." diff --git a/source/digits_hits/src/GateSinglesDigitizerMessenger.cc b/source/digits_hits/src/GateSinglesDigitizerMessenger.cc index d1fd0bd38..b603cd077 100755 --- a/source/digits_hits/src/GateSinglesDigitizerMessenger.cc +++ b/source/digits_hits/src/GateSinglesDigitizerMessenger.cc @@ -9,7 +9,7 @@ See LICENSE.md for further details /*! \class GateSinglesDigitizerMessenger - Last modification (Adaptation to GND): May 2023 by Mohamed-Jordan Soumano mjsoumano@yahoo.com + Last modification (Adaptation to GND): July 2024 by Marc Granado-Gonzalez marc.granado@universite-paris-saclay.fr */ #include "GateSinglesDigitizer.hh" @@ -55,6 +55,7 @@ See LICENSE.md for further details #include "GateAdderComptPhotIdeal.hh" #include "GateClustering.hh" #include "GateTimeDelay.hh" +#include "GateVirtualSegmentationSD.hh" /* #include "GateLocalTimeDelay.hh" @@ -126,7 +127,7 @@ const G4String& GateSinglesDigitizerMessenger::DumpMap() { - static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection"; + static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection virtualSegmentation"; return theList; } @@ -267,6 +268,13 @@ else if (childTypeName=="multipleRejection") m_digitizer->AddNewModule(newDM); } + else if (childTypeName=="virtualSegmentation") + { + newDM = new GateVirtualSegmentationSD(m_digitizer, DMname); + m_digitizer->AddNewModule(newDM); + } + + /* diff --git a/source/digits_hits/src/GateSpatialResolution.cc b/source/digits_hits/src/GateSpatialResolution.cc index 1621a0318..f6a19ed06 100644 --- a/source/digits_hits/src/GateSpatialResolution.cc +++ b/source/digits_hits/src/GateSpatialResolution.cc @@ -313,6 +313,9 @@ void GateSpatialResolution::UpdateVolumeID() } } + + + void GateSpatialResolution::DescribeMyself(size_t indent ) { if(m_fwhm) diff --git a/source/digits_hits/src/GateVirtualSegmentationSD.cc b/source/digits_hits/src/GateVirtualSegmentationSD.cc new file mode 100644 index 000000000..afe840117 --- /dev/null +++ b/source/digits_hits/src/GateVirtualSegmentationSD.cc @@ -0,0 +1,559 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +// OK GND 2022 +/*! + \class GateVirtualSegmentationSD (by marc.granado@universite-paris-saclay.fr) + \brief for discretizing the position through a virtual segmentation of a monolithic crystal + + - For each volume the local position of the interactions within the crystal are virtually segmented. + the X,Y,Z vector is translated into the IDs of virtual divisions within the crystal + occupying the levels of SubmoduleID, CrystalID and LayerID. + +*/ + + + +#include "GateVirtualSegmentationSD.hh" +#include "GateVirtualSegmentationSDMessenger.hh" +#include "GateDigi.hh" + +#include "GateSpatialResolution.hh" +#include "GateSpatialResolutionMessenger.hh" + +#include "GateBoxComponent.hh" + +#include "GateDigitizerMgr.hh" + +#include "G4Box.hh" +#include "G4SystemOfUnits.hh" +#include "G4EventManager.hh" +#include "G4Event.hh" +#include "G4SDManager.hh" +#include "G4DigiManager.hh" +#include "G4ios.hh" +#include "G4UnitsTable.hh" + +#include +#include +#include +#include + + + + + +GateVirtualSegmentationSD::GateVirtualSegmentationSD(GateSinglesDigitizer *digitizer, G4String name) + :GateVDigitizerModule(name,"digitizerMgr/"+digitizer->GetSD()->GetName()+"/SinglesDigitizer/"+digitizer->m_digitizerName+"/"+name,digitizer,digitizer->GetSD()), + + m_IsFirstEntry(1), + m_nameAxis("XYZ"), + useMacroGenerator(0), + + m_pitch(0), + m_pitchX(0), + m_pitchY(0), + m_pitchZ(0), + + m_outputDigi(0), + m_OutputDigiCollection(0), + m_digitizer(digitizer), + + pitchX(0), + pitchY(0), + pitchZ(0), + + nBinsX(0), + nBinsY(0), + nBinsZ(0), + + depthX(5), + depthY(4), + depthZ(3), + + xLength(0), + yLength(0), + zLength(0) + + +{ + G4String colName = digitizer->GetOutputName() ; + collectionName.push_back(colName); + m_Messenger = new GateVirtualSegmentationSDMessenger(this); +} + + + + +GateVirtualSegmentationSD::~GateVirtualSegmentationSD() +{ + delete m_Messenger; + +} + + + + + + +void GateVirtualSegmentationSD::Digitize() +{ + + GateVSystem* m_system = ((GateSinglesDigitizer*)this->GetDigitizer())->GetSystem(); + + if (m_system==NULL) G4Exception( "GateVirtualSegmentationSD::Digitize", "Digitize", FatalException, + "Failed to get the system corresponding to that digitizer. Abort.\n"); + + if (!m_system->CheckIfEnoughLevelsAreDefined()) + { + GateError( " *** ERROR*** GateVirtualSegmentationSD::Digitize. Not all defined geometry levels has their mother levels defined." + "(Ex.: for cylindricalPET, the levels are: rsector, module, submodule, crystal). If you have defined submodule, you have to have resector and module defined as well." + "Please, add them to your geometry macro in /gate/systems/cylindricalPET/XXX/attach YYY. Abort.\n"); + } + + + + m_systemDepth = m_system->GetTreeDepth(); + + G4String digitizerName = m_digitizer->m_digitizerName; + G4String outputCollName = m_digitizer-> GetOutputName(); + + m_OutputDigiCollection = new GateDigiCollection(GetName(),outputCollName); // to create the Digi Collection + + G4DigiManager* DigiMan = G4DigiManager::GetDMpointer(); + + + + GateDigiCollection* IDC = 0; + IDC = (GateDigiCollection*) (DigiMan->GetDigiCollection(m_DCID)); + + GateDigi* inputDigi; + std::vector< GateDigi* >* OutputDigiCollectionVector = m_OutputDigiCollection->GetVector (); + std::vector::iterator iter; + + + + + if (IDC) + { + G4int n_digi = IDC->entries(); + + //loop over input digits + for (G4int i=0;iGetVolumeID().size()){ + G4Box* box = dynamic_cast(inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid()); + xLength = 2*box->GetXHalfLength(); + yLength = 2*box->GetYHalfLength(); + zLength = 2*box->GetZHalfLength(); + } + else GateError( " *** ERROR*** No volumeID in inputDigi!"); + + + + + if (m_IsFirstEntry){ + + if(pitchX > xLength){ pitchX=xLength; + GateWarning("The size provided for pitchX is bigger than the crystal syze! The length of the crystal has been selected as the new pitchX");} + if(pitchY > yLength){ pitchY=yLength; + GateWarning("The size provided for pitchY is bigger than the crystal syze! The length of the crystal has been selected as the new pitchY");} + if(pitchZ > zLength){ pitchZ=zLength; + GateWarning("The size provided for pitchZ is bigger than the crystal syze! The length of the crystal has been selected as the new pitchX");} + + SetParameters(); + + + + if (useMacroGenerator){ + /* + std::string path_to_script = "/home/granado/GATE_projects/MonoCrystals/tests/mac/"; + std::string path_to_macros = "/home/granado/GATE_projects/MonoCrystals/tests/mac/"; + //Once this is implemented as a GateTool it will be: + //std::string command = "macro_converter"; + std::string command = "python3 "+path_to_script+"macro_converter.py"; + command +=" -d "+path_to_macros+"digitizer.mac"; + command +=" -g "+path_to_macros+"geometry_pseudo-crystal.mac"; + + if (pitchX) command +=" -x "+std::to_string(pitchX); + if (pitchY) command +=" -y "+std::to_string(pitchY); + if (pitchZ) command +=" -z "+std::to_string(pitchZ); + + if (nBinsX) command +=" --binsX "+std::to_string(nBinsX); + if (nBinsY) command +=" --binsY "+std::to_string(nBinsY); + if (nBinsZ) command +=" --binsZ "+std::to_string(nBinsZ); + + int result = system(command.c_str()); + if (result ==0) std::cout<<"The macro converter script has been executed corectly"<1) + { + G4cout << "Merged previous pulse for volume " << inputDigi->GetVolumeID() + << " with new pulse of energy " << G4BestUnit(inputDigi->GetEnergy(),"Energy") <<".\n" + << "Resulting pulse is: \n" + << **iter << Gateendl << Gateendl ; + } + + + m_outputDigi = new GateDigi(*inputDigi); + m_outputDigi->SetEnergyIniTrack(-1); + m_outputDigi->SetEnergyFin(-1); + + + //Accessing the local position of the hit to define and set the new VirtualID + + G4ThreeVector localPos = inputDigi->GetLocalPos(); + + if(m_nameAxis.find('X') != std::string::npos){ + SetVirtualID(nBinsX,pitchX,localPos.getX(),depthX); + } + if(m_nameAxis.find('Y') != std::string::npos){ + SetVirtualID(nBinsY,pitchY,localPos.getY(),depthY); + } + if(m_nameAxis.find('Z') != std::string::npos){ + SetVirtualID(nBinsZ,pitchZ,localPos.getZ(),depthZ); + } + + + + + + + + if (nVerboseLevel>1) + G4cout << "Created new pulse for volume " << inputDigi->GetVolumeID() << ".\n" + << "Resulting pulse is: \n" + << *m_outputDigi << Gateendl << Gateendl ; + /// !!!!!! The following line should be kept !!!! -> inserts the outputdigi to collection + m_OutputDigiCollection->insert(m_outputDigi); + + + ////// ** End of the part from ProcessOnePulse + + + if (nVerboseLevel==1) { + G4cout << "[GateVirtualSegmentationSD::Digitize]: returning output pulse-list with " << OutputDigiCollectionVector->size() << " entries\n"; + for (iter=OutputDigiCollectionVector->begin(); iter!= OutputDigiCollectionVector->end() ; ++iter) + G4cout << **iter << Gateendl; + G4cout << Gateendl; + } + /// *** End of the part from ProcessPulseList + } //loop over input digits + } //IDC + else + { + if (nVerboseLevel>1) + G4cout << "[GateVirtualSegmentationSD::Digitize]: input digi collection is null -> nothing to do\n\n"; + return; + } + StoreDigiCollection(m_OutputDigiCollection); + + +} + + + + + + + + + +void GateVirtualSegmentationSD::SetNameAxis(const G4String ¶m) +{ + m_nameAxis=param; +} + + +/////////////////////////////////////////// +////////////// Methods of DM ////////////// +/////////////////////////////////////////// + + +//TODO: This function could be further optimised (but only used once) + +G4double GateVirtualSegmentationSD::calculatePitch(G4double crystal_size, G4double target_pitch) { + // Target pitch size + + double best_pitch = 0; + double min_diff = 999; + int num_pitches = 1; + + while (true) { + double pitch = crystal_size / num_pitches; + + // Check if pitch meets both conditions + if (pitch <= target_pitch && std::fmod(crystal_size, pitch) == 0) { + double diff = std::fabs(pitch - target_pitch); + if (diff <= min_diff) { + min_diff = diff; + best_pitch = pitch; + } + } + + // Break if pitch is smaller than a threshold, to prevent infinite loop + if ((best_pitch <= target_pitch)&&(best_pitch!=0)) { + break; + } + else if (num_pitches>1000){ + std::cout<<"ERROR! No pitch found in 1000 iterations! Size = "<(m_digitizer->FindDigitizerModule("digitizerMgr/" + +m_digitizer->GetSD()->GetName() + +"/SinglesDigitizer/" + + m_digitizer->GetName() + + + "/spatialResolution")); + + + //Setting the pitch size provided by the user + + if(m_pitch){ + if(m_nameAxis.find('X') != std::string::npos) + {pitchX = m_pitch;} + + if(m_nameAxis.find('Y') != std::string::npos) + {pitchY = m_pitch;} + + if(m_nameAxis.find('Z') != std::string::npos) + {pitchZ = m_pitch;} + + + if (m_pitchX!=0 || m_pitchY!=0 ||m_pitchZ!=0) + {GateWarning("The values provided for target pitch are ambiguous. Only the value of 'pitch' was taken, any value provided for 'pitchX', 'pitchY' and 'pitchZ' was ignored");} + } + + else if (m_pitchX!=0 || m_pitchY!=0 ||m_pitchZ!=0){ + + if(m_pitchX && m_nameAxis.find('X') != std::string::npos) + {pitchX = m_pitchX;} + + if(m_pitchY && m_nameAxis.find('Y') != std::string::npos) + {pitchY = m_pitchY;} + + if(m_pitchZ && m_nameAxis.find('Z') != std::string::npos) + {pitchZ = m_pitchZ;} + + } + + else if (digi_SpatialResolution){ + if (digi_SpatialResolution->GetFWHM()) + { + + if(m_nameAxis.find('X') != std::string::npos) + {pitchX = 0.5*digi_SpatialResolution->GetFWHM();} + + if(m_nameAxis.find('Y') != std::string::npos) + {pitchY = 0.5*digi_SpatialResolution->GetFWHM();} + + if(m_nameAxis.find('Z') != std::string::npos) + {pitchZ = 0.5*digi_SpatialResolution->GetFWHM();} + } + + + else if(digi_SpatialResolution->GetFWHMx() || digi_SpatialResolution->GetFWHMy()|| digi_SpatialResolution->GetFWHMz()){ + + + + if(digi_SpatialResolution->GetFWHMxdistrib()||digi_SpatialResolution->GetFWHMydistrib()||digi_SpatialResolution->GetFWHMxydistrib2D()) + { + GateError("***ERROR*** No value of the target pitch has been provided and no value can be obtained from the spatial resolution distribution. /n Please provide a value for the pitch that is at least half of the minimum value of the distribution. "); + } + + + + + if(digi_SpatialResolution->GetFWHMx() && m_nameAxis.find('X') != std::string::npos) + {pitchX = 0.5*digi_SpatialResolution->GetFWHMx();} + + if(digi_SpatialResolution->GetFWHMy() && m_nameAxis.find('Y') != std::string::npos) + {pitchY = 0.5*digi_SpatialResolution->GetFWHMy();} + + if(digi_SpatialResolution->GetFWHMz() && m_nameAxis.find('Z') != std::string::npos) + {pitchZ = 0.5*digi_SpatialResolution->GetFWHMz();} + } + } + + if(pitchX == 0 && m_nameAxis.find('X') != std::string::npos) + {GateError("***ERROR*** Virtual setmentation in X axis has been selected but no value for FWHM was provided."); + } + else if (pitchX!=0){ + + std::cout<<"Calculating pitchX with = "<GetDigitizer())->GetSystem(); + + + GateSystemComponent* m_submoduleComponent = m_system->FindComponent("submodule"); + GateSystemComponent* m_crystalComponent = m_system->FindComponent("crystal"); + GateSystemComponent* m_layer0Component= m_system->FindComponent("layer0"); + std::cout<<"System depth is: "< GetVolumeNumber(); + if (repeaterNumber > 1) + GateError( " *** ERROR*** GateVirtualSegmentationSD::Digitize. Crystal level needs to be empty (no repeaters) to use VirtualSegmentationSD with XYZ segmentation"); + + + + if (m_nameAxis.size()>=2) + { + repeaterNumber = m_crystalComponent-> GetVolumeNumber(); + if (repeaterNumber > 1) + GateError( " *** ERROR*** GateVirtualSegmentationSD::Digitize. Crystal level needs to be empty (no repeaters) to use VirtualSegmentationSD with XYZ segmentation"); + + } + + + + + if (m_nameAxis.size()==3) + { + repeaterNumber = m_submoduleComponent-> GetVolumeNumber(); + if (repeaterNumber >1) + GateError( " *** ERROR*** GateVirtualSegmentationSD::Digitize. Submodule level needs to be empty (no repeaters) to use VirtualSegmentationSD with XYZ segmentation"); + + } + +} + + + + + + diff --git a/source/digits_hits/src/GateVirtualSegmentationSDMessenger.cc b/source/digits_hits/src/GateVirtualSegmentationSDMessenger.cc new file mode 100644 index 000000000..4e18140d3 --- /dev/null +++ b/source/digits_hits/src/GateVirtualSegmentationSDMessenger.cc @@ -0,0 +1,135 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +// OK GND 2022 +/*This class is not used by GATE ! + The purpose of this class is to help to create new users digitizer module(DM). + Please, check GateVirtualSegmentationSD.cc for more detals + */ + + + +/*! + \class GateVirtualSegmentationSD (by marc.granado@universite-paris-saclay.fr) + \brief for discretizing the position through a virtual segmentation of a monolithic crystal + + - For each volume the local position of the interactions within the crystal are virtually segmented. + the X,Y,Z vector is translated into the IDs of virtual divisions within the crystal + occupying the levels of SubmoduleID, CrystalID and LayerID. +*/ + + +#include "GateVirtualSegmentationSDMessenger.hh" +#include "GateVirtualSegmentationSD.hh" +#include "GateDigitizerMgr.hh" + +#include "G4SystemOfUnits.hh" +#include "G4UIcmdWithADouble.hh" +#include "G4UIcmdWithADoubleAndUnit.hh" +#include "G4UIcmdWithABool.hh" +#include "G4UIcmdWithAString.hh" +#include "G4UIdirectory.hh" + + + +GateVirtualSegmentationSDMessenger::GateVirtualSegmentationSDMessenger (GateVirtualSegmentationSD* VirtualSegmentationSD) +:GateClockDependentMessenger(VirtualSegmentationSD), + m_VirtualSegmentationSD(VirtualSegmentationSD) +{ + G4String guidance; + G4String cmdName; + + + cmdName = GetDirectoryName()+"nameAxis"; + nameAxisCmd = new G4UIcmdWithAString(cmdName,this); + nameAxisCmd ->SetGuidance("Provide the number of axis that need to be virtually segmented, XYZ, XY, XZ, or YZ"); + nameAxisCmd ->SetCandidates("XYZ XY XZ YZ"); + + cmdName = GetDirectoryName() + "useMacroGenerator"; + useMacroGeneratorCmd = new G4UIcmdWithABool(cmdName,this); + useMacroGeneratorCmd->SetGuidance("To be set true, if there is a need to create a new geometry macro for CASToR"); + + + cmdName = GetDirectoryName() + "pitch"; + pitchCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this);//new G4UIcmdWithADouble(cmdName,this); + pitchCmd->SetGuidance("Set the size of the pitch for all the selected dimensions"); + + + + cmdName = GetDirectoryName()+"pitchX"; + pitchXCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this);//new G4UIcmdWithADouble(cmdName,this); + pitchXCmd->SetGuidance("Set the size of the pitch for the X dimension"); + + cmdName = GetDirectoryName()+"pitchY"; + pitchYCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this);//new G4UIcmdWithADouble(cmdName,this); + pitchYCmd->SetGuidance("Set the size of the pitch for the Y dimension"); + + cmdName = GetDirectoryName()+"pitchZ"; + pitchZCmd = new G4UIcmdWithADoubleAndUnit(cmdName,this);//new G4UIcmdWithADouble(cmdName,this); + pitchZCmd->SetGuidance("Set the size of the pitch for the Z dimension"); + +} + + +GateVirtualSegmentationSDMessenger::~GateVirtualSegmentationSDMessenger() +{ + delete nameAxisCmd; + delete useMacroGeneratorCmd; + delete pitchCmd; + delete pitchXCmd; + delete pitchYCmd; + delete pitchZCmd; +} + + +void GateVirtualSegmentationSDMessenger::SetNewValue(G4UIcommand * aCommand,G4String newValue) +{ + + if (aCommand == nameAxisCmd) + { + m_VirtualSegmentationSD->SetNameAxis(newValue); + } + + else if (aCommand == pitchCmd) + { + m_VirtualSegmentationSD->SetPitch(pitchCmd->GetNewDoubleValue(newValue)); + } + else if (aCommand == pitchXCmd) + { + m_VirtualSegmentationSD->SetPitchX(pitchCmd->GetNewDoubleValue(newValue)); + + } + + else if (aCommand == pitchYCmd) + { + m_VirtualSegmentationSD->SetPitchY(pitchCmd->GetNewDoubleValue(newValue)); + } + else if (aCommand == pitchZCmd) + { + m_VirtualSegmentationSD->SetPitchZ(pitchCmd->GetNewDoubleValue(newValue)); + } + else if ( aCommand==useMacroGeneratorCmd ) + { m_VirtualSegmentationSD->SetUseMacroGenerator(useMacroGeneratorCmd->GetNewBoolValue(newValue)); } + else + { + GateClockDependentMessenger::SetNewValue(aCommand,newValue); + } +} + + + + + + + + + + + + +