From a90b7b37e4d1bdae33335ea87adca74c2c00c6b9 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Thu, 9 Feb 2023 15:53:47 +0100 Subject: [PATCH 01/12] Update .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 8a76c2668..54beb012a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +build/ _build/ _wiki/ *~ @@ -19,3 +20,4 @@ cmake-build-* /source/io_old/ /source/save_bug/ +.cache/ From c4102dc5de60c7be61fa153d5c03aadc978e3af9 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Thu, 9 Feb 2023 15:54:38 +0100 Subject: [PATCH 02/12] Use CMAKE_CXX_STANDARD 17 --- CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 777d05a1c..713924a46 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,6 +67,9 @@ ELSE() ENDIF() set(CMAKE_CXX_STANDARD ${Geant4_CXX_STANDARD}) +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + #========================================================= # Refuse to build with MT geant4, allow dev override IF(Geant4_multithreaded_FOUND) From 5850325c50954abb35e278ccc15e58007e72de8a Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Thu, 9 Feb 2023 15:55:02 +0100 Subject: [PATCH 03/12] Add GateBioDoseActor --- .../digits_hits/include/GateBioDoseActor.hh | 120 +++++++++ .../include/GateBioDoseActorMessenger.hh | 42 ++++ source/digits_hits/src/GateBioDoseActor.cc | 228 ++++++++++++++++++ .../src/GateBioDoseActorMessenger.cc | 60 +++++ 4 files changed, 450 insertions(+) create mode 100644 source/digits_hits/include/GateBioDoseActor.hh create mode 100644 source/digits_hits/include/GateBioDoseActorMessenger.hh create mode 100644 source/digits_hits/src/GateBioDoseActor.cc create mode 100644 source/digits_hits/src/GateBioDoseActorMessenger.cc diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh new file mode 100644 index 000000000..365eb8961 --- /dev/null +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -0,0 +1,120 @@ +#ifndef GATE_SOURCE_DIGITS_HITS_INCLUDE_GATEBIODOSEACTOR_HH +#define GATE_SOURCE_DIGITS_HITS_INCLUDE_GATEBIODOSEACTOR_HH + +/*---------------------- +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 +----------------------*/ + +/*! + \class GateBioDoseActor + \author Éloïse Salles, Alexis Pereda , Yasmine Ali +*/ + +#include + +#include "GateBioDoseActorMessenger.hh" +#include "GateImageWithStatistic.hh" +#include "GateVImageActor.hh" + +class GateBioDoseActor: public GateVImageActor +{ +public: + struct Deposited { + double alpha; + double beta; + double energy; + }; + + struct Coefficients { + double a, b; + }; + + struct AlphaBetaCoefficients { + Coefficients alpha; + Coefficients beta; + }; + + using VoxelIndex = int; + using DepositedMap = std::map; + + using Fragment = std::pair; + using AlphaBetaInterpolTable = std::map; + + using EnergyMaxForZ = std::map; + +public: + ~GateBioDoseActor() override = default; + + FCT_FOR_AUTO_CREATOR_ACTOR(GateBioDoseActor) + + //----------------------------------------------------------------------------- + // Constructs the sensor + void Construct() override; + + // G4 + void BeginOfRunAction(const G4Run* r) override; + void EndOfRunAction(const G4Run* r) override; + void BeginOfEventAction(const G4Event* event) override; + void UserSteppingActionInVoxel(const int index, const G4Step* step) override; + + // Do nothing but needed because pure virtual + void UserPreTrackActionInVoxel(const int, const G4Track*) override {} + void UserPostTrackActionInVoxel(const int, const G4Track*) override {} + + // Saves the data collected to the file + void SaveData() override; + void ResetData() override; + + // Scorer related + void Initialize(G4HCofThisEvent*) override {} + void EndOfEvent(G4HCofThisEvent*) override {} + + // Messenger + void SetAlphaRef(G4double alphaRef) { _alphaRef = alphaRef; } + void SetBetaRef(G4double betaRef) { _betaRef = betaRef; } + void SetBioDoseImageFilename(G4String s) { _bioDoseFilename = s; } + void SetCellLine(G4String s) { _cellLine = s; } + void SetBioPhysicalModel(G4String s) { _bioPhysicalModel = s; } + void SetSOBPWeight(G4double d) { _SOBPWeight = d; } + + // Input database + void BuildDatabase(); + Coefficients Interpol(double x1, double x2, double y1, double y2); + +protected: + GateBioDoseActor(G4String name, G4int depth = 0); + + void ApplyDeposit(int index, DepositedMap::iterator& it, double energyDep); + +private: + //Counters + int _currentEvent; + + GateBioDoseActorMessenger _messenger; + + EnergyMaxForZ _energyMaxForZ; + + // Initialisation + G4String _dataBase; + G4String _cellLine; + G4String _bioPhysicalModel; + double _alphaRef, _betaRef; //manual implanted + double _massOfVoxel; // G4_WATER + + // Maps + DepositedMap _depositedMap; + AlphaBetaInterpolTable _alphaBetaInterpolTable; + + // Images + G4double _SOBPWeight; + G4String _bioDoseFilename; + GateImageWithStatistic _bioDoseImage; +}; + +MAKE_AUTO_CREATOR_ACTOR(BioDoseActor, GateBioDoseActor) + +#endif diff --git a/source/digits_hits/include/GateBioDoseActorMessenger.hh b/source/digits_hits/include/GateBioDoseActorMessenger.hh new file mode 100644 index 000000000..a2e6f1531 --- /dev/null +++ b/source/digits_hits/include/GateBioDoseActorMessenger.hh @@ -0,0 +1,42 @@ +#ifndef GATE_SOURCE_DIGITS_HITS_INCLUDE_GATEBIODOSEACTORMESSENGER_HH +#define GATE_SOURCE_DIGITS_HITS_INCLUDE_GATEBIODOSEACTORMESSENGER_HH + +/*---------------------- +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 +----------------------*/ + +/*! + \class GateBioDoseActorMessenger + \author Éloïse Salles, Alexis Pereda , Yasmine Ali +*/ + +#include "G4UIcmdWithABool.hh" +#include "G4UIcmdWithADouble.hh" +#include "GateImageActorMessenger.hh" +#include + +class GateBioDoseActor; + +class GateBioDoseActorMessenger: public GateImageActorMessenger { +public: + GateBioDoseActorMessenger(GateBioDoseActor* sensor); + + void BuildCommands(G4String base) override; + void SetNewValue(G4UIcommand*, G4String) override; + +protected: + GateBioDoseActor* pBioDoseActor; + + std::unique_ptr pAlphaRefCmd; + std::unique_ptr pBetaRefCmd; + std::unique_ptr pImageFilenameCmd; + std::unique_ptr pCellLineCmd; + std::unique_ptr pBioPhysicalModelCmd; + std::unique_ptr pSOBPWeightCmd; +}; + +#endif diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc new file mode 100644 index 000000000..e5e4ba77e --- /dev/null +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -0,0 +1,228 @@ +/*---------------------- +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 +----------------------*/ + +#include "G4EmParameters.hh" +#include "GateBioDoseActor.hh" + +#define GATE_BUFFERSIZE + +//----------------------------------------------------------------------------- +GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): + GateVImageActor(name, depth), + _currentEvent(0), + _messenger(this), + _alphaRef(-1), + _betaRef(-1) +{ + GateDebugMessageInc("Actor", 4, "GateBioDoseActor() -- begin\n"); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::Construct() { + GateDebugMessageInc("Actor", 4, "GateBioDoseActor -- Construct - begin\n"); + GateVImageActor::Construct(); + + // Find G4_WATER + G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"); + // Find OtherMaterial + //G4NistManager::Instance()->FindOrBuildMaterial(mOtherMaterial); + G4EmParameters::Instance()->SetBuildCSDARange(true); + + // Enable callbacks BioDose + EnableBeginOfRunAction(true); + EnableEndOfRunAction(true); + EnableBeginOfEventAction(true); + EnablePreUserTrackingAction(false); + EnablePostUserTrackingAction(false); + EnableUserSteppingAction(true); + + // Output + _bioDoseFilename = G4String(removeExtension(mSaveFilename)) + "-BioDose." + G4String(getExtension(mSaveFilename)); + SetOriginTransformAndFlagToImage(_bioDoseImage); + _bioDoseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _bioDoseImage.Allocate(); + _bioDoseImage.SetFilename(_bioDoseFilename); + + //ResetData(); + + /////////////////////////////////////////////////////////////////////////////////////////// + //Just matrix information + G4cout << "Memory space to store physical dose into " << mResolution.x() * mResolution.y() * mResolution.z() << " voxels has been allocated " << G4endl; + + //Calculate the mass of voxel for dose estimation + _massOfVoxel = ((1 * (mVoxelSize.x() / cm * mVoxelSize.y() / cm * mVoxelSize.z() / cm)) * 1E-3); //G4_WATER (density : g.cm-3) + + // SOBP + if (_SOBPWeight == 0) { _SOBPWeight = 1; } + + //Building the cell line information + _dataBase = "data/" + _cellLine + "_" + _bioPhysicalModel + ".db"; + BuildDatabase(); + + if(_alphaRef < 0 || _betaRef < 0) + GateError("BioDoseActor " << GetName() << ": setAlphaRef and setBetaRef must be done"); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- // ok bio dose +void GateBioDoseActor::BuildDatabase() { + std::ifstream f(_dataBase); + if(!f) GateError("BioDoseActor " << GetName() << ": unable to open file '" << _dataBase << "'"); + + int nZ = 0; + double prevKineticEnergy = 1, prevAlpha = 1, prevBeta =1; + + for(std::string line; std::getline(f, line); ) { + std::istringstream iss(line); + std::string firstCol; + + iss >> firstCol; + + if(firstCol == "Fragment") { + if(nZ != 0) // prevKineticEnergy is the maximum kinetic energy for current nZ + _energyMaxForZ[nZ] = prevKineticEnergy; + + iss >> nZ; + prevKineticEnergy = 1; + prevAlpha = 1; + prevBeta = 1; + } else if(nZ != 0) { + double kineticEnergy, alpha, beta; + std::istringstream{firstCol} >> kineticEnergy; + iss >> alpha >> beta; + + auto alphaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevAlpha, alpha); + auto betaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevBeta, beta); + + // Saving the in the input databse + Fragment fragment{nZ, kineticEnergy}; + _alphaBetaInterpolTable[fragment] = {alphaCoeff, betaCoeff}; + + prevKineticEnergy = kineticEnergy; + prevAlpha = alpha; + prevBeta = beta; + } else { + GateError("BioDoseActor " << GetName() << ": bad database format in '" << _dataBase << "'"); + } + } + + if(nZ != 0) // last line read; prevKineticEnergy is the maximum kinetic energy for current nZ + _energyMaxForZ[nZ] = prevKineticEnergy; +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, double y1, double y2) { + //Function for a 1D linear interpolation. It returns a pair of a and b coefficients + double a = (y2 - y1) / (x2 - x1); + double b = y1 - x1 * a; + return {a, b}; +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +/// Save data +void GateBioDoseActor::SaveData() { + // Calculate Alpha Beta mix + // Fiding deposited alpha and beta stored in the maps for the right index i + for(auto const& [index, deposited]: _depositedMap) { + // Alpha Beta mix (final) + double alphamix = (deposited.alpha / deposited.energy); + double betamix = (deposited.beta / deposited.energy); + + // Calculate dose, dosebio + constexpr double MeVtoJoule = 1.60218e-13; + double dose = deposited.energy * MeVtoJoule / _massOfVoxel; + double dosebio = 0; + + if(dose != 0 && alphamix != 0 && betamix != 0) + dosebio = (-_alphaRef + std::sqrt((_alphaRef * _alphaRef) + 4 * _betaRef * (alphamix * dose + betamix * (dose * dose)))) / (2 * _betaRef); + + // Save data + _bioDoseImage.AddValue(index, dosebio); + } + //------------------------------------------------------------- + GateVActor::SaveData(); + _bioDoseImage.SaveData(_currentEvent); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::BeginOfRunAction(const G4Run* r) { + GateVActor::BeginOfRunAction(r); + GateDebugMessage("Actor", 3, "GateBioDoseActor -- Begin of Run\n"); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- // PAS DANS DOSE ACTOR +void GateBioDoseActor::EndOfRunAction(const G4Run* r) { + GateVActor::EndOfRunAction(r); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// Callback at each event +void GateBioDoseActor::BeginOfEventAction(const G4Event* e) { + GateVActor::BeginOfEventAction(e); + ++_currentEvent; +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* step) { + double const energyDep = step->GetTotalEnergyDeposit(); + + // Get information from step + // Particle + G4int nZ = step->GetTrack()->GetDefinition()->GetAtomicNumber(); + double kineticEnergyPerNucleon = (step->GetPreStepPoint()->GetKineticEnergy()) / (step->GetTrack()->GetDefinition()->GetAtomicMass()); //OK + + DepositedMap::iterator it = _depositedMap.find(index); + + if(energyDep > 0) { + if(it == std::end(_depositedMap)) { + _depositedMap[index] = {0, 0, 0}; + it = _depositedMap.find(index); + } + + auto& deposited = (*it).second; + + // Accumulate energy inconditionnaly + deposited.energy += energyDep; + + // Accumulation of alpha/beta if ion type if known + // -> check if the ion type is known + if(_energyMaxForZ.count(nZ)) { + //The max values in the database aren't being taking into account + //so for now it's coded like this to be sure the code takes them into account + double energyMax = _energyMaxForZ.at(nZ); + + AlphaBetaInterpolTable::const_iterator itr2; + if (kineticEnergyPerNucleon >= energyMax) { + // If the kinetic energy is the maximum value in the alpha beta tables, + // we have to use the a and b coefficient for this maximum value + Fragment fragmentKineticEnergyMax{nZ, energyMax}; + itr2 = _alphaBetaInterpolTable.find(fragmentKineticEnergyMax); + } else { + // We pair the ion type and the kinetic energy + Fragment fragmentKineticEnergy{nZ, kineticEnergyPerNucleon}; + itr2 = _alphaBetaInterpolTable.upper_bound(fragmentKineticEnergy); + } + + // Calculation of EZ, alphaDep and betaDep (K = a*EZ+b*E) + auto const& interpol = (*itr2).second; + + double ez = kineticEnergyPerNucleon * energyDep; + double alphaDep = interpol.alpha.a * ez + interpol.alpha.b * energyDep; + double betaDep = interpol.beta.a * ez + interpol.beta.b * energyDep; + + // Accumulate alpha/beta + deposited.alpha += alphaDep; + deposited.beta += betaDep; + } + } +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::ResetData() { + // _bioDoseImage.Reset(); +} +//----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc new file mode 100644 index 000000000..a7d52ee3c --- /dev/null +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -0,0 +1,60 @@ +/*---------------------- +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 +----------------------*/ + +#include "GateBioDoseActorMessenger.hh" +#include "GateBioDoseActor.hh" + +//----------------------------------------------------------------------------- +GateBioDoseActorMessenger::GateBioDoseActorMessenger(GateBioDoseActor* sensor): + GateImageActorMessenger(sensor), + pBioDoseActor(sensor) +{ + BuildCommands(baseName + sensor->GetObjectName()); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActorMessenger::BuildCommands(G4String base) { + G4String n; + + n = base + "/setAlphaRef"; + pAlphaRefCmd = std::make_unique(n, this); + pAlphaRefCmd->SetGuidance("See [...] for values from publications"); + + n = base + "/setBetaRef"; + pBetaRefCmd = std::make_unique(n, this); + pBetaRefCmd->SetGuidance("See [...] for values from publications"); + + n = base + "/setBioDoseImageFilename"; + pImageFilenameCmd = std::make_unique(n, this); + pImageFilenameCmd->SetGuidance("For an image for the biological dose as output"); + + n = base + "/setCellLine"; + pCellLineCmd = std::make_unique(n, this); + pCellLineCmd->SetGuidance("Pick a cell line to irradiate"); + + n = base + "/setBioPhysicalModel"; + pBioPhysicalModelCmd = std::make_unique(n, this); + pBioPhysicalModelCmd->SetGuidance("Pick a model to work with"); + + n = base + "/setWeightSOBP"; + pSOBPWeightCmd = std::make_unique(n, this); + pSOBPWeightCmd->SetGuidance("If passive SOBP, set bragg peak weight"); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String newValue) { + if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(newValue)); + else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(newValue)); + else if(cmd == pImageFilenameCmd.get()) pBioDoseActor->SetBioDoseImageFilename(newValue); + else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(newValue); + else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(newValue); + else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(newValue)); + + GateImageActorMessenger::SetNewValue(cmd, newValue); +} +//----------------------------------------------------------------------------- From 78db7905c97a9a653526c513ab68fc5a1eca26b9 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Mon, 13 Feb 2023 16:50:07 +0100 Subject: [PATCH 04/12] Add options to BioDoseActor to save dose, alpha-/beta-mix and RBE --- .../digits_hits/include/GateBioDoseActor.hh | 22 ++++- .../include/GateBioDoseActorMessenger.hh | 7 +- source/digits_hits/src/GateBioDoseActor.cc | 87 +++++++++++++++---- .../src/GateBioDoseActorMessenger.cc | 49 ++++++++--- 4 files changed, 133 insertions(+), 32 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index 365eb8961..f822f9c7f 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -76,11 +76,16 @@ public: // Messenger void SetAlphaRef(G4double alphaRef) { _alphaRef = alphaRef; } void SetBetaRef(G4double betaRef) { _betaRef = betaRef; } - void SetBioDoseImageFilename(G4String s) { _bioDoseFilename = s; } void SetCellLine(G4String s) { _cellLine = s; } void SetBioPhysicalModel(G4String s) { _bioPhysicalModel = s; } void SetSOBPWeight(G4double d) { _SOBPWeight = d; } + void SetEnableDose(bool e) { _enableDose = e; } + void SetEnableBioDose(bool e) { _enableBioDose = e; } + void SetEnableAlphaMix(bool e) { _enableAlphaMix = e; } + void SetEnableBetaMix(bool e) { _enableBetaMix = e; } + void SetEnableRBE(bool e) { _enableRBE = e; } + // Input database void BuildDatabase(); Coefficients Interpol(double x1, double x2, double y1, double y2); @@ -105,14 +110,25 @@ private: double _alphaRef, _betaRef; //manual implanted double _massOfVoxel; // G4_WATER + G4double _SOBPWeight; + // Maps DepositedMap _depositedMap; AlphaBetaInterpolTable _alphaBetaInterpolTable; // Images - G4double _SOBPWeight; - G4String _bioDoseFilename; GateImageWithStatistic _bioDoseImage; + GateImageWithStatistic _doseImage; + GateImageWithStatistic _alphaMixImage; + GateImageWithStatistic _betaMixImage; + GateImageWithStatistic _RBEImage; + + // Outputs + bool _enableDose; + bool _enableBioDose; + bool _enableAlphaMix; + bool _enableBetaMix; + bool _enableRBE; }; MAKE_AUTO_CREATOR_ACTOR(BioDoseActor, GateBioDoseActor) diff --git a/source/digits_hits/include/GateBioDoseActorMessenger.hh b/source/digits_hits/include/GateBioDoseActorMessenger.hh index a2e6f1531..0d2efba3c 100644 --- a/source/digits_hits/include/GateBioDoseActorMessenger.hh +++ b/source/digits_hits/include/GateBioDoseActorMessenger.hh @@ -33,10 +33,15 @@ protected: std::unique_ptr pAlphaRefCmd; std::unique_ptr pBetaRefCmd; - std::unique_ptr pImageFilenameCmd; std::unique_ptr pCellLineCmd; std::unique_ptr pBioPhysicalModelCmd; std::unique_ptr pSOBPWeightCmd; + + std::unique_ptr pEnableDoseCmd; + std::unique_ptr pEnableBioDoseCmd; + std::unique_ptr pEnableAlphaMixCmd; + std::unique_ptr pEnableBetaMixCmd; + std::unique_ptr pEnableRBECmd; }; #endif diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index e5e4ba77e..ae652b86b 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -17,7 +17,12 @@ GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): _currentEvent(0), _messenger(this), _alphaRef(-1), - _betaRef(-1) + _betaRef(-1), + _enableDose(false), + _enableBioDose(true), + _enableAlphaMix(false), + _enableBetaMix(false), + _enableRBE(false) { GateDebugMessageInc("Actor", 4, "GateBioDoseActor() -- begin\n"); } @@ -41,12 +46,56 @@ void GateBioDoseActor::Construct() { EnablePostUserTrackingAction(false); EnableUserSteppingAction(true); - // Output - _bioDoseFilename = G4String(removeExtension(mSaveFilename)) + "-BioDose." + G4String(getExtension(mSaveFilename)); - SetOriginTransformAndFlagToImage(_bioDoseImage); - _bioDoseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _bioDoseImage.Allocate(); - _bioDoseImage.SetFilename(_bioDoseFilename); + // Outputs + { + G4String basename = removeExtension(mSaveFilename); + G4String ext = getExtension(mSaveFilename); + + if(_enableDose) { + G4String filename = basename + "_dose." + ext; + + SetOriginTransformAndFlagToImage(_doseImage); + _doseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _doseImage.Allocate(); + _doseImage.SetFilename(filename); + } + + if(_enableBioDose) { + G4String filename = basename + "_biodose." + ext; + + SetOriginTransformAndFlagToImage(_bioDoseImage); + _bioDoseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _bioDoseImage.Allocate(); + _bioDoseImage.SetFilename(filename); + } + + if(_enableAlphaMix) { + G4String filename = basename + "_alphamix." + ext; + + SetOriginTransformAndFlagToImage(_alphaMixImage); + _alphaMixImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _alphaMixImage.Allocate(); + _alphaMixImage.SetFilename(filename); + } + + if(_enableBetaMix) { + G4String filename = basename + "_betamix." + ext; + + SetOriginTransformAndFlagToImage(_betaMixImage); + _betaMixImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _betaMixImage.Allocate(); + _betaMixImage.SetFilename(filename); + } + + if(_enableRBE) { + G4String filename = basename + "_rbe." + ext; + + SetOriginTransformAndFlagToImage(_RBEImage); + _RBEImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + _RBEImage.Allocate(); + _RBEImage.SetFilename(filename); + } + } //ResetData(); @@ -58,7 +107,7 @@ void GateBioDoseActor::Construct() { _massOfVoxel = ((1 * (mVoxelSize.x() / cm * mVoxelSize.y() / cm * mVoxelSize.z() / cm)) * 1E-3); //G4_WATER (density : g.cm-3) // SOBP - if (_SOBPWeight == 0) { _SOBPWeight = 1; } + if(_SOBPWeight == 0) { _SOBPWeight = 1; } //Building the cell line information _dataBase = "data/" + _cellLine + "_" + _bioPhysicalModel + ".db"; @@ -135,17 +184,26 @@ void GateBioDoseActor::SaveData() { // Calculate dose, dosebio constexpr double MeVtoJoule = 1.60218e-13; double dose = deposited.energy * MeVtoJoule / _massOfVoxel; - double dosebio = 0; + double biodose = 0; if(dose != 0 && alphamix != 0 && betamix != 0) - dosebio = (-_alphaRef + std::sqrt((_alphaRef * _alphaRef) + 4 * _betaRef * (alphamix * dose + betamix * (dose * dose)))) / (2 * _betaRef); - - // Save data - _bioDoseImage.AddValue(index, dosebio); + biodose = (-_alphaRef + std::sqrt((_alphaRef * _alphaRef) + 4 * _betaRef * (alphamix * dose + betamix * (dose * dose)))) / (2 * _betaRef); + + // Write data + if(_enableDose) _doseImage.AddValue(index, dose); + if(_enableBioDose) _bioDoseImage.AddValue(index, biodose); + if(_enableAlphaMix) _alphaMixImage.AddValue(index, alphamix); + if(_enableBetaMix) _betaMixImage.AddValue(index, betamix); + if(_enableRBE) _RBEImage.AddValue(index, dose > 0? biodose/dose : 0); } //------------------------------------------------------------- GateVActor::SaveData(); - _bioDoseImage.SaveData(_currentEvent); + + if(_enableDose) _doseImage.SaveData(_currentEvent); + if(_enableBioDose) _bioDoseImage.SaveData(_currentEvent); + if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); + if(_enableBetaMix) _betaMixImage.SaveData(_currentEvent); + if(_enableRBE) _RBEImage.SaveData(_currentEvent); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -223,6 +281,5 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActor::ResetData() { - // _bioDoseImage.Reset(); } //----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc index a7d52ee3c..c7577d164 100644 --- a/source/digits_hits/src/GateBioDoseActorMessenger.cc +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -8,6 +8,8 @@ See LICENSE.md for further details #include "GateBioDoseActorMessenger.hh" #include "GateBioDoseActor.hh" +#include +#include //----------------------------------------------------------------------------- GateBioDoseActorMessenger::GateBioDoseActorMessenger(GateBioDoseActor* sensor): @@ -29,10 +31,6 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { pBetaRefCmd = std::make_unique(n, this); pBetaRefCmd->SetGuidance("See [...] for values from publications"); - n = base + "/setBioDoseImageFilename"; - pImageFilenameCmd = std::make_unique(n, this); - pImageFilenameCmd->SetGuidance("For an image for the biological dose as output"); - n = base + "/setCellLine"; pCellLineCmd = std::make_unique(n, this); pCellLineCmd->SetGuidance("Pick a cell line to irradiate"); @@ -44,17 +42,42 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { n = base + "/setWeightSOBP"; pSOBPWeightCmd = std::make_unique(n, this); pSOBPWeightCmd->SetGuidance("If passive SOBP, set bragg peak weight"); + + // enable outputs + n = base + "/enableDose"; + pEnableDoseCmd = std::make_unique(n, this); + pEnableDoseCmd->SetGuidance("Enable dose output"); + + n = base + "/enableBioDose"; + pEnableBioDoseCmd = std::make_unique(n, this); + pEnableBioDoseCmd->SetGuidance("Enable biodose output (default: true)"); + + n = base + "/enableAlphaMix"; + pEnableAlphaMixCmd = std::make_unique(n, this); + pEnableAlphaMixCmd->SetGuidance("Enable alpha mix output"); + + n = base + "/enableBetaMix"; + pEnableBetaMixCmd = std::make_unique(n, this); + pEnableBetaMixCmd->SetGuidance("Enable beta mix output"); + + n = base + "/enableRBE"; + pEnableRBECmd = std::make_unique(n, this); + pEnableRBECmd->SetGuidance("Enable RBE output"); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String newValue) { - if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(newValue)); - else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(newValue)); - else if(cmd == pImageFilenameCmd.get()) pBioDoseActor->SetBioDoseImageFilename(newValue); - else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(newValue); - else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(newValue); - else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(newValue)); - - GateImageActorMessenger::SetNewValue(cmd, newValue); +void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { + if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(value)); + else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(value)); + else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(value); + else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(value); + else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(value)); + else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); + else if(cmd == pEnableBioDoseCmd.get()) pBioDoseActor->SetEnableBioDose(pEnableBioDoseCmd->GetNewBoolValue(value)); + else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); + else if(cmd == pEnableBetaMixCmd.get()) pBioDoseActor->SetEnableBetaMix(pEnableBetaMixCmd->GetNewBoolValue(value)); + else if(cmd == pEnableRBECmd.get()) pBioDoseActor->SetEnableRBE(pEnableRBECmd->GetNewBoolValue(value)); + + GateImageActorMessenger::SetNewValue(cmd, value); } //----------------------------------------------------------------------------- From f451cf602a29e9ef1140946a4b9f34c58fd8875d Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Tue, 7 Mar 2023 11:05:57 +0100 Subject: [PATCH 05/12] Fix mass calculation per voxel in biodose actor --- source/digits_hits/include/GateBioDoseActor.hh | 3 ++- source/digits_hits/src/GateBioDoseActor.cc | 13 ++++++++----- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index f822f9c7f..30a27662f 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -27,6 +27,8 @@ public: double alpha; double beta; double energy; + + double mass; }; struct Coefficients { @@ -108,7 +110,6 @@ private: G4String _cellLine; G4String _bioPhysicalModel; double _alphaRef, _betaRef; //manual implanted - double _massOfVoxel; // G4_WATER G4double _SOBPWeight; diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index ae652b86b..d18bfd1ab 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -103,9 +103,6 @@ void GateBioDoseActor::Construct() { //Just matrix information G4cout << "Memory space to store physical dose into " << mResolution.x() * mResolution.y() * mResolution.z() << " voxels has been allocated " << G4endl; - //Calculate the mass of voxel for dose estimation - _massOfVoxel = ((1 * (mVoxelSize.x() / cm * mVoxelSize.y() / cm * mVoxelSize.z() / cm)) * 1E-3); //G4_WATER (density : g.cm-3) - // SOBP if(_SOBPWeight == 0) { _SOBPWeight = 1; } @@ -175,7 +172,6 @@ GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, /// Save data void GateBioDoseActor::SaveData() { // Calculate Alpha Beta mix - // Fiding deposited alpha and beta stored in the maps for the right index i for(auto const& [index, deposited]: _depositedMap) { // Alpha Beta mix (final) double alphamix = (deposited.alpha / deposited.energy); @@ -183,9 +179,12 @@ void GateBioDoseActor::SaveData() { // Calculate dose, dosebio constexpr double MeVtoJoule = 1.60218e-13; - double dose = deposited.energy * MeVtoJoule / _massOfVoxel; + double dose = 0; double biodose = 0; + if(deposited.mass != 0) + dose = deposited.energy * MeVtoJoule / deposited.mass; + if(dose != 0 && alphamix != 0 && betamix != 0) biodose = (-_alphaRef + std::sqrt((_alphaRef * _alphaRef) + 4 * _betaRef * (alphamix * dose + betamix * (dose * dose)))) / (2 * _betaRef); @@ -246,6 +245,10 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* // Accumulate energy inconditionnaly deposited.energy += energyDep; + auto current_material = step->GetPreStepPoint()->GetMaterial(); + double density = current_material->GetDensity(); + deposited.mass = ((1 * (mVoxelSize.x() / cm * mVoxelSize.y() / cm * mVoxelSize.z() / cm)) * density); + // Accumulation of alpha/beta if ion type if known // -> check if the ion type is known if(_energyMaxForZ.count(nZ)) { From b1e49500c53d78e47a134b1c30d1b19d42e34dd9 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Mon, 24 Apr 2023 10:53:43 +0200 Subject: [PATCH 06/12] Improve code, add edep output and fix computation --- .../digits_hits/include/GateBioDoseActor.hh | 13 +- .../include/GateBioDoseActorMessenger.hh | 1 + source/digits_hits/src/GateBioDoseActor.cc | 225 +++++++++--------- .../src/GateBioDoseActorMessenger.cc | 5 + 4 files changed, 132 insertions(+), 112 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index 30a27662f..a5193b4d3 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -25,10 +25,9 @@ class GateBioDoseActor: public GateVImageActor public: struct Deposited { double alpha; - double beta; + double sqrtBeta; double energy; - - double mass; + double dose; }; struct Coefficients { @@ -37,7 +36,7 @@ public: struct AlphaBetaCoefficients { Coefficients alpha; - Coefficients beta; + Coefficients sqrtBeta; }; using VoxelIndex = int; @@ -82,6 +81,7 @@ public: void SetBioPhysicalModel(G4String s) { _bioPhysicalModel = s; } void SetSOBPWeight(G4double d) { _SOBPWeight = d; } + void SetEnableEdep(bool e) { _enableEdep = e; } void SetEnableDose(bool e) { _enableDose = e; } void SetEnableBioDose(bool e) { _enableBioDose = e; } void SetEnableAlphaMix(bool e) { _enableAlphaMix = e; } @@ -119,17 +119,22 @@ private: // Images GateImageWithStatistic _bioDoseImage; + GateImageWithStatistic _edepImage; GateImageWithStatistic _doseImage; GateImageWithStatistic _alphaMixImage; GateImageWithStatistic _betaMixImage; GateImageWithStatistic _RBEImage; // Outputs + bool _enableEdep; bool _enableDose; bool _enableBioDose; bool _enableAlphaMix; bool _enableBetaMix; bool _enableRBE; + + int _eventCount = 0; + int _eventWithKnownIonCount = 0; }; MAKE_AUTO_CREATOR_ACTOR(BioDoseActor, GateBioDoseActor) diff --git a/source/digits_hits/include/GateBioDoseActorMessenger.hh b/source/digits_hits/include/GateBioDoseActorMessenger.hh index 0d2efba3c..235cc1836 100644 --- a/source/digits_hits/include/GateBioDoseActorMessenger.hh +++ b/source/digits_hits/include/GateBioDoseActorMessenger.hh @@ -37,6 +37,7 @@ protected: std::unique_ptr pBioPhysicalModelCmd; std::unique_ptr pSOBPWeightCmd; + std::unique_ptr pEnableEdepCmd; std::unique_ptr pEnableDoseCmd; std::unique_ptr pEnableBioDoseCmd; std::unique_ptr pEnableAlphaMixCmd; diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index d18bfd1ab..6eb49118e 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -8,6 +8,9 @@ See LICENSE.md for further details #include "G4EmParameters.hh" #include "GateBioDoseActor.hh" +#include "GateImage.hh" +#include "GateImageWithStatistic.hh" +#include #define GATE_BUFFERSIZE @@ -18,6 +21,7 @@ GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): _messenger(this), _alphaRef(-1), _betaRef(-1), + _enableEdep(false), _enableDose(false), _enableBioDose(true), _enableAlphaMix(false), @@ -51,53 +55,24 @@ void GateBioDoseActor::Construct() { G4String basename = removeExtension(mSaveFilename); G4String ext = getExtension(mSaveFilename); - if(_enableDose) { - G4String filename = basename + "_dose." + ext; - - SetOriginTransformAndFlagToImage(_doseImage); - _doseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _doseImage.Allocate(); - _doseImage.SetFilename(filename); - } - - if(_enableBioDose) { - G4String filename = basename + "_biodose." + ext; - - SetOriginTransformAndFlagToImage(_bioDoseImage); - _bioDoseImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _bioDoseImage.Allocate(); - _bioDoseImage.SetFilename(filename); - } - - if(_enableAlphaMix) { - G4String filename = basename + "_alphamix." + ext; - - SetOriginTransformAndFlagToImage(_alphaMixImage); - _alphaMixImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _alphaMixImage.Allocate(); - _alphaMixImage.SetFilename(filename); - } - - if(_enableBetaMix) { - G4String filename = basename + "_betamix." + ext; - - SetOriginTransformAndFlagToImage(_betaMixImage); - _betaMixImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _betaMixImage.Allocate(); - _betaMixImage.SetFilename(filename); - } - - if(_enableRBE) { - G4String filename = basename + "_rbe." + ext; - - SetOriginTransformAndFlagToImage(_RBEImage); - _RBEImage.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - _RBEImage.Allocate(); - _RBEImage.SetFilename(filename); - } + auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix) { + G4String filename = basename + "_" + suffix + "." + ext; + + SetOriginTransformAndFlagToImage(image); + image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + image.Allocate(); + image.SetFilename(filename); + }; + + if(_enableEdep) setupImage(_edepImage, "edep"); + if(_enableDose) setupImage(_doseImage, "dose"); + if(_enableBioDose) setupImage(_bioDoseImage, "biodose"); + if(_enableAlphaMix) setupImage(_alphaMixImage, "alphamix"); + if(_enableBetaMix) setupImage(_betaMixImage, "betamix"); + if(_enableRBE) setupImage(_RBEImage, "rbe"); } - //ResetData(); + ResetData(); /////////////////////////////////////////////////////////////////////////////////////////// //Just matrix information @@ -142,11 +117,11 @@ void GateBioDoseActor::BuildDatabase() { iss >> alpha >> beta; auto alphaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevAlpha, alpha); - auto betaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevBeta, beta); + auto sqrtBetaCoeff = Interpol(prevKineticEnergy, kineticEnergy, std::sqrt(prevBeta), std::sqrt(beta)); // Saving the in the input databse Fragment fragment{nZ, kineticEnergy}; - _alphaBetaInterpolTable[fragment] = {alphaCoeff, betaCoeff}; + _alphaBetaInterpolTable[fragment] = {alphaCoeff, sqrtBetaCoeff}; prevKineticEnergy = kineticEnergy; prevAlpha = alpha; @@ -169,35 +144,44 @@ GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -/// Save data void GateBioDoseActor::SaveData() { + GateDebugMessageInc("Actor", 4, "GateBioDoseActor::SaveData() known ion events / total events: " << _eventWithKnownIonCount << " / " << _eventCount << "\n"); + // Calculate Alpha Beta mix for(auto const& [index, deposited]: _depositedMap) { // Alpha Beta mix (final) - double alphamix = (deposited.alpha / deposited.energy); - double betamix = (deposited.beta / deposited.energy); + double alphamix = 0; + double betamix = 0; + + if(deposited.energy != 0) { + alphamix = (deposited.alpha / deposited.energy); + betamix = (deposited.sqrtBeta / deposited.energy) * (deposited.sqrtBeta / deposited.energy); + } - // Calculate dose, dosebio - constexpr double MeVtoJoule = 1.60218e-13; - double dose = 0; - double biodose = 0; + // Calculate biological dose and RBE + double biodose = 0; + double rbe = 0; - if(deposited.mass != 0) - dose = deposited.energy * MeVtoJoule / deposited.mass; + if(deposited.dose > 0 && alphamix != 0 && betamix != 0) { + auto const sqAlphaRef = _alphaRef * _alphaRef; + auto const sqDose = deposited.dose * deposited.dose; + biodose = (-_alphaRef + std::sqrt(sqAlphaRef + 4 * _betaRef * (alphamix * deposited.dose + betamix * sqDose))) / (2 * _betaRef); - if(dose != 0 && alphamix != 0 && betamix != 0) - biodose = (-_alphaRef + std::sqrt((_alphaRef * _alphaRef) + 4 * _betaRef * (alphamix * dose + betamix * (dose * dose)))) / (2 * _betaRef); + rbe = biodose / deposited.dose; + } // Write data - if(_enableDose) _doseImage.AddValue(index, dose); - if(_enableBioDose) _bioDoseImage.AddValue(index, biodose); - if(_enableAlphaMix) _alphaMixImage.AddValue(index, alphamix); - if(_enableBetaMix) _betaMixImage.AddValue(index, betamix); - if(_enableRBE) _RBEImage.AddValue(index, dose > 0? biodose/dose : 0); + if(_enableEdep) _edepImage.SetValue(index, deposited.energy); + if(_enableDose) _doseImage.SetValue(index, deposited.dose); + if(_enableBioDose) _bioDoseImage.SetValue(index, biodose); + if(_enableAlphaMix) _alphaMixImage.SetValue(index, alphamix); + if(_enableBetaMix) _betaMixImage.SetValue(index, betamix); + if(_enableRBE) _RBEImage.SetValue(index, rbe); } //------------------------------------------------------------- GateVActor::SaveData(); + if(_enableEdep) _edepImage.SaveData(_currentEvent); if(_enableDose) _doseImage.SaveData(_currentEvent); if(_enableBioDose) _bioDoseImage.SaveData(_currentEvent); if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); @@ -225,64 +209,89 @@ void GateBioDoseActor::BeginOfEventAction(const G4Event* e) { //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* step) { - double const energyDep = step->GetTotalEnergyDeposit(); + double const weight = step->GetTrack()->GetWeight(); + double const energyDep = step->GetTotalEnergyDeposit() * weight; + + if(energyDep == 0) return; + if(index < 0) return; + + DepositedMap::iterator it = _depositedMap.find(index); + if(it == std::end(_depositedMap)) { + _depositedMap[index] = {0, 0, 0, 0}; + it = _depositedMap.find(index); + } + + auto& deposited = (*it).second; + + // Accumulate energy inconditionnaly + deposited.energy += energyDep; + + if(_enableDose || _enableBioDose || _enableRBE) { + decltype(_doseImage)* image = nullptr; + if(_enableDose) image = &_doseImage; + else if(_enableBioDose) image = &_bioDoseImage; + else if(_enableRBE) image = &_RBEImage; + + auto currentMaterial = step->GetPreStepPoint()->GetMaterial(); + double density = currentMaterial->GetDensity(); + double mass = image->GetVoxelVolume() * density; + + deposited.dose += energyDep / mass / CLHEP::gray; + } // Get information from step // Particle G4int nZ = step->GetTrack()->GetDefinition()->GetAtomicNumber(); double kineticEnergyPerNucleon = (step->GetPreStepPoint()->GetKineticEnergy()) / (step->GetTrack()->GetDefinition()->GetAtomicMass()); //OK - DepositedMap::iterator it = _depositedMap.find(index); + ++_eventCount; - if(energyDep > 0) { - if(it == std::end(_depositedMap)) { - _depositedMap[index] = {0, 0, 0}; - it = _depositedMap.find(index); - } + // Accumulation of alpha/beta if ion type if known + // -> check if the ion type is known + if(_energyMaxForZ.count(nZ)) { + ++_eventWithKnownIonCount; - auto& deposited = (*it).second; - - // Accumulate energy inconditionnaly - deposited.energy += energyDep; - - auto current_material = step->GetPreStepPoint()->GetMaterial(); - double density = current_material->GetDensity(); - deposited.mass = ((1 * (mVoxelSize.x() / cm * mVoxelSize.y() / cm * mVoxelSize.z() / cm)) * density); - - // Accumulation of alpha/beta if ion type if known - // -> check if the ion type is known - if(_energyMaxForZ.count(nZ)) { - //The max values in the database aren't being taking into account - //so for now it's coded like this to be sure the code takes them into account - double energyMax = _energyMaxForZ.at(nZ); - - AlphaBetaInterpolTable::const_iterator itr2; - if (kineticEnergyPerNucleon >= energyMax) { - // If the kinetic energy is the maximum value in the alpha beta tables, - // we have to use the a and b coefficient for this maximum value - Fragment fragmentKineticEnergyMax{nZ, energyMax}; - itr2 = _alphaBetaInterpolTable.find(fragmentKineticEnergyMax); - } else { - // We pair the ion type and the kinetic energy - Fragment fragmentKineticEnergy{nZ, kineticEnergyPerNucleon}; - itr2 = _alphaBetaInterpolTable.upper_bound(fragmentKineticEnergy); - } - - // Calculation of EZ, alphaDep and betaDep (K = a*EZ+b*E) - auto const& interpol = (*itr2).second; - - double ez = kineticEnergyPerNucleon * energyDep; - double alphaDep = interpol.alpha.a * ez + interpol.alpha.b * energyDep; - double betaDep = interpol.beta.a * ez + interpol.beta.b * energyDep; - - // Accumulate alpha/beta - deposited.alpha += alphaDep; - deposited.beta += betaDep; + //The max values in the database aren't being taking into account + //so for now it's coded like this to be sure the code takes them into account + double energyMax = _energyMaxForZ.at(nZ); + + AlphaBetaInterpolTable::const_iterator itr2; + if (kineticEnergyPerNucleon >= energyMax) { + // If the kinetic energy is the maximum value in the alpha beta tables, + // we have to use the a and b coefficient for this maximum value + Fragment fragmentKineticEnergyMax{nZ, energyMax}; + itr2 = _alphaBetaInterpolTable.find(fragmentKineticEnergyMax); + } else { + // We pair the ion type and the kinetic energy + Fragment fragmentKineticEnergy{nZ, kineticEnergyPerNucleon}; + itr2 = _alphaBetaInterpolTable.upper_bound(fragmentKineticEnergy); } + + // Calculation of EZ, alphaDep and betaDep (K = a*EZ+b*E) + auto const& interpol = (*itr2).second; + + double alpha = (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b); + double sqrtBeta = (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b); + + if(alpha < 0) alpha = 0; + if(sqrtBeta < 0) sqrtBeta = 0; + + double alphaDep = alpha * energyDep; + double sqrtBetaDep = sqrtBeta * energyDep; + + // Accumulate alpha/beta + deposited.alpha += alphaDep; + deposited.sqrtBeta += sqrtBetaDep; } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActor::ResetData() { + if(_enableEdep) _edepImage.Reset(); + if(_enableDose) _doseImage.Reset(); + if(_enableBioDose) _bioDoseImage.Reset(); + if(_enableAlphaMix) _alphaMixImage.Reset(); + if(_enableBetaMix) _betaMixImage.Reset(); + if(_enableRBE) _RBEImage.Reset(); } //----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc index c7577d164..dc3822514 100644 --- a/source/digits_hits/src/GateBioDoseActorMessenger.cc +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -44,6 +44,10 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { pSOBPWeightCmd->SetGuidance("If passive SOBP, set bragg peak weight"); // enable outputs + n = base + "/enableEdep"; + pEnableEdepCmd = std::make_unique(n, this); + pEnableEdepCmd->SetGuidance("Enable deposited energy output"); + n = base + "/enableDose"; pEnableDoseCmd = std::make_unique(n, this); pEnableDoseCmd->SetGuidance("Enable dose output"); @@ -72,6 +76,7 @@ void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(value); else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(value); else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(value)); + else if(cmd == pEnableEdepCmd.get()) pBioDoseActor->SetEnableEdep(pEnableEdepCmd->GetNewBoolValue(value)); else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); else if(cmd == pEnableBioDoseCmd.get()) pBioDoseActor->SetEnableBioDose(pEnableBioDoseCmd->GetNewBoolValue(value)); else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); From 47838566024d0c1608c4fce1d4686f9e77a49ede Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Fri, 1 Dec 2023 16:10:07 +0100 Subject: [PATCH 07/12] Statistical uncertainty (biodose), fix formulas --- .../digits_hits/include/GateBioDoseActor.hh | 52 ++- .../include/GateBioDoseActorMessenger.hh | 6 +- source/digits_hits/src/GateBioDoseActor.cc | 423 ++++++++++++++---- .../src/GateBioDoseActorMessenger.cc | 20 +- 4 files changed, 392 insertions(+), 109 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index a5193b4d3..c3ed226c1 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -41,6 +41,7 @@ public: using VoxelIndex = int; using DepositedMap = std::map; + using VoxelIndices = std::set; using Fragment = std::pair; using AlphaBetaInterpolTable = std::map; @@ -48,8 +49,6 @@ public: using EnergyMaxForZ = std::map; public: - ~GateBioDoseActor() override = default; - FCT_FOR_AUTO_CREATOR_ACTOR(GateBioDoseActor) //----------------------------------------------------------------------------- @@ -60,6 +59,7 @@ public: void BeginOfRunAction(const G4Run* r) override; void EndOfRunAction(const G4Run* r) override; void BeginOfEventAction(const G4Event* event) override; + void EndOfEventAction(const G4Event* event) override; void UserSteppingActionInVoxel(const int index, const G4Step* step) override; // Do nothing but needed because pure virtual @@ -75,28 +75,28 @@ public: void EndOfEvent(G4HCofThisEvent*) override {} // Messenger + void SetDoseScaleFactor(G4double doseScaleFactor) { _doseScaleFactor = doseScaleFactor; } void SetAlphaRef(G4double alphaRef) { _alphaRef = alphaRef; } void SetBetaRef(G4double betaRef) { _betaRef = betaRef; } - void SetCellLine(G4String s) { _cellLine = s; } - void SetBioPhysicalModel(G4String s) { _bioPhysicalModel = s; } - void SetSOBPWeight(G4double d) { _SOBPWeight = d; } + void SetCellLine(G4String s) { _cellLine = std::move(s); } + void SetBioPhysicalModel(G4String s) { _bioPhysicalModel = std::move(s); } + void SetSOBPWeight(G4double d) { _sobpWeight = d; } void SetEnableEdep(bool e) { _enableEdep = e; } void SetEnableDose(bool e) { _enableDose = e; } void SetEnableBioDose(bool e) { _enableBioDose = e; } void SetEnableAlphaMix(bool e) { _enableAlphaMix = e; } - void SetEnableBetaMix(bool e) { _enableBetaMix = e; } + void SetEnableSqrtBetaMix(bool e) { _enableSqrtBetaMix = e; } void SetEnableRBE(bool e) { _enableRBE = e; } + void SetEnableUncertainties(bool e) { _enableUncertainties = e; } // Input database void BuildDatabase(); - Coefficients Interpol(double x1, double x2, double y1, double y2); + static Coefficients Interpol(double x1, double x2, double y1, double y2); protected: GateBioDoseActor(G4String name, G4int depth = 0); - void ApplyDeposit(int index, DepositedMap::iterator& it, double energyDep); - private: //Counters int _currentEvent; @@ -109,32 +109,52 @@ private: G4String _dataBase; G4String _cellLine; G4String _bioPhysicalModel; - double _alphaRef, _betaRef; //manual implanted + double _alphaRef, _betaRef; + double _doseScaleFactor = 1.; - G4double _SOBPWeight; + G4double _sobpWeight; // Maps DepositedMap _depositedMap; AlphaBetaInterpolTable _alphaBetaInterpolTable; + VoxelIndices _eventVoxelIndices; + // Images + GateImageWithStatistic _eventCountImage; + GateImageWithStatistic _bioDoseImage; GateImageWithStatistic _edepImage; GateImageWithStatistic _doseImage; GateImageWithStatistic _alphaMixImage; - GateImageWithStatistic _betaMixImage; - GateImageWithStatistic _RBEImage; + GateImageWithStatistic _sqrtBetaMixImage; + GateImageWithStatistic _rbeImage; + + GateImageWithStatistic _biodoseUncertaintyImage; + GateImageWithStatistic _eventEdepImage; + GateImageWithStatistic _eventDoseImage; + GateImageWithStatistic _squaredDoseImage; + GateImageWithStatistic _eventAlphaImage; + GateImageWithStatistic _squaredAlphaMixImage; + GateImageWithStatistic _eventSqrtBetaImage; + GateImageWithStatistic _squaredSqrtBetaMixImage; + + GateImageWithStatistic _alphaMixSqrtBetaMixImage; + GateImageWithStatistic _alphaMixDoseImage; + GateImageWithStatistic _sqrtBetaMixDoseImage; // Outputs bool _enableEdep; bool _enableDose; bool _enableBioDose; bool _enableAlphaMix; - bool _enableBetaMix; + bool _enableSqrtBetaMix; bool _enableRBE; + bool _enableUncertainties; - int _eventCount = 0; - int _eventWithKnownIonCount = 0; + // Extra information + int _stepCount = 0; + int _stepWithKnownIonCount = 0; }; MAKE_AUTO_CREATOR_ACTOR(BioDoseActor, GateBioDoseActor) diff --git a/source/digits_hits/include/GateBioDoseActorMessenger.hh b/source/digits_hits/include/GateBioDoseActorMessenger.hh index 235cc1836..ec388fcfc 100644 --- a/source/digits_hits/include/GateBioDoseActorMessenger.hh +++ b/source/digits_hits/include/GateBioDoseActorMessenger.hh @@ -28,9 +28,10 @@ public: void BuildCommands(G4String base) override; void SetNewValue(G4UIcommand*, G4String) override; -protected: +private: GateBioDoseActor* pBioDoseActor; + std::unique_ptr pDoseScaleFactorCmd; std::unique_ptr pAlphaRefCmd; std::unique_ptr pBetaRefCmd; std::unique_ptr pCellLineCmd; @@ -41,8 +42,9 @@ protected: std::unique_ptr pEnableDoseCmd; std::unique_ptr pEnableBioDoseCmd; std::unique_ptr pEnableAlphaMixCmd; - std::unique_ptr pEnableBetaMixCmd; + std::unique_ptr pEnableSqrtBetaMixCmd; std::unique_ptr pEnableRBECmd; + std::unique_ptr pEnableUncertaintyCmd; }; #endif diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index 6eb49118e..8889dbef1 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -1,4 +1,4 @@ -/*---------------------- +/*-- Copyright (C): OpenGATE Collaboration This software is distributed under the terms @@ -16,17 +16,19 @@ See LICENSE.md for further details //----------------------------------------------------------------------------- GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): - GateVImageActor(name, depth), + GateVImageActor(std::move(name), depth), _currentEvent(0), _messenger(this), _alphaRef(-1), _betaRef(-1), + _sobpWeight(0), _enableEdep(false), - _enableDose(false), + _enableDose(true), _enableBioDose(true), _enableAlphaMix(false), - _enableBetaMix(false), - _enableRBE(false) + _enableSqrtBetaMix(false), + _enableRBE(false), + _enableUncertainties(true) { GateDebugMessageInc("Actor", 4, "GateBioDoseActor() -- begin\n"); } @@ -46,6 +48,7 @@ void GateBioDoseActor::Construct() { EnableBeginOfRunAction(true); EnableEndOfRunAction(true); EnableBeginOfEventAction(true); + EnableEndOfEventAction(true); EnablePreUserTrackingAction(false); EnablePostUserTrackingAction(false); EnableUserSteppingAction(true); @@ -55,21 +58,42 @@ void GateBioDoseActor::Construct() { G4String basename = removeExtension(mSaveFilename); G4String ext = getExtension(mSaveFilename); - auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix) { - G4String filename = basename + "_" + suffix + "." + ext; - + auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { SetOriginTransformAndFlagToImage(image); image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); image.Allocate(); - image.SetFilename(filename); + + if(!suffix.empty()) { + G4String filename = basename + "_" + suffix + "." + ext; + image.SetFilename(filename); + } }; - if(_enableEdep) setupImage(_edepImage, "edep"); - if(_enableDose) setupImage(_doseImage, "dose"); - if(_enableBioDose) setupImage(_bioDoseImage, "biodose"); - if(_enableAlphaMix) setupImage(_alphaMixImage, "alphamix"); - if(_enableBetaMix) setupImage(_betaMixImage, "betamix"); - if(_enableRBE) setupImage(_RBEImage, "rbe"); + setupImage(_eventCountImage); + + if(_enableEdep) setupImage(_edepImage, "edep"); + if(_enableDose) setupImage(_doseImage, "dose"); + setupImage(_bioDoseImage, "biodose"); + setupImage(_alphaMixImage, "alphamix"); + setupImage(_sqrtBetaMixImage, "sqrtbetamix"); + if(_enableRBE) setupImage(_rbeImage, "rbe"); + + setupImage(_eventEdepImage); + setupImage(_eventDoseImage); + setupImage(_eventAlphaImage); + setupImage(_eventSqrtBetaImage); + + if(_enableUncertainties) { + setupImage(_biodoseUncertaintyImage, "biodose_uncertainty"); + + setupImage(_squaredDoseImage); + setupImage(_squaredAlphaMixImage); + setupImage(_squaredSqrtBetaMixImage); + + setupImage(_alphaMixSqrtBetaMixImage); + setupImage(_alphaMixDoseImage); + setupImage(_sqrtBetaMixDoseImage); + } } ResetData(); @@ -79,7 +103,7 @@ void GateBioDoseActor::Construct() { G4cout << "Memory space to store physical dose into " << mResolution.x() * mResolution.y() * mResolution.z() << " voxels has been allocated " << G4endl; // SOBP - if(_SOBPWeight == 0) { _SOBPWeight = 1; } + if(_sobpWeight == 0) { _sobpWeight = 1; } //Building the cell line information _dataBase = "data/" + _cellLine + "_" + _bioPhysicalModel + ".db"; @@ -89,13 +113,15 @@ void GateBioDoseActor::Construct() { GateError("BioDoseActor " << GetName() << ": setAlphaRef and setBetaRef must be done"); } //----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- // ok bio dose +//----------------------------------------------------------------------------- void GateBioDoseActor::BuildDatabase() { std::ifstream f(_dataBase); if(!f) GateError("BioDoseActor " << GetName() << ": unable to open file '" << _dataBase << "'"); int nZ = 0; - double prevKineticEnergy = 1, prevAlpha = 1, prevBeta =1; + double prevKineticEnergy = 1; + double prevAlpha = 1; + double prevBeta =1; for(std::string line; std::getline(f, line); ) { std::istringstream iss(line); @@ -112,9 +138,12 @@ void GateBioDoseActor::BuildDatabase() { prevAlpha = 1; prevBeta = 1; } else if(nZ != 0) { - double kineticEnergy, alpha, beta; + double kineticEnergy = 0; + double alpha = 0; + double beta = 0; std::istringstream{firstCol} >> kineticEnergy; - iss >> alpha >> beta; + iss >> alpha; + iss >> beta; auto alphaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevAlpha, alpha); auto sqrtBetaCoeff = Interpol(prevKineticEnergy, kineticEnergy, std::sqrt(prevBeta), std::sqrt(beta)); @@ -147,64 +176,264 @@ GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, void GateBioDoseActor::SaveData() { GateDebugMessageInc("Actor", 4, "GateBioDoseActor::SaveData() known ion events / total events: " << _eventWithKnownIonCount << " / " << _eventCount << "\n"); - // Calculate Alpha Beta mix - for(auto const& [index, deposited]: _depositedMap) { - // Alpha Beta mix (final) - double alphamix = 0; - double betamix = 0; + /* ******************* */ + // TODO remove + G4String basename = removeExtension(mSaveFilename); + G4String ext = getExtension(mSaveFilename); + auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { + SetOriginTransformAndFlagToImage(image); + image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + image.Allocate(); - if(deposited.energy != 0) { - alphamix = (deposited.alpha / deposited.energy); - betamix = (deposited.sqrtBeta / deposited.energy) * (deposited.sqrtBeta / deposited.energy); + if(!suffix.empty()) { + G4String filename = basename + "_" + suffix + "." + ext; + image.SetFilename(filename); } + }; + + GateImageWithStatistic partAlphaMixImg; + GateImageWithStatistic partSqrtBetaMixImg; + GateImageWithStatistic partDoseImg; + GateImageWithStatistic partAlphaMixSqrtBetaMixImg; + GateImageWithStatistic partAlphaMixDoseImg; + GateImageWithStatistic partSqrtBetaMixDoseImg; + GateImageWithStatistic partSquaredBioDoseImg; + + GateImageWithStatistic squaredPdBioDoseAlphaMixImg; + GateImageWithStatistic squaredPdBioDoseSqrtBetaMixImg; + GateImageWithStatistic squaredPdBioDoseDoseImg; + GateImageWithStatistic varAlphaMixImg; + GateImageWithStatistic varSqrtBetaMixImg; + GateImageWithStatistic varDoseImg; + + setupImage(partAlphaMixImg, "part_alphamix"); + setupImage(partSqrtBetaMixImg, "part_sqrtbetamix"); + setupImage(partDoseImg, "part_dose"); + setupImage(partAlphaMixSqrtBetaMixImg, "part_alphamixsqrtbetamix"); + setupImage(partAlphaMixDoseImg, "part_alphamixdose"); + setupImage(partSqrtBetaMixDoseImg, "part_sqrtbetamixdose"); + setupImage(partSquaredBioDoseImg, "part_squaredbiodose"); + + setupImage(squaredPdBioDoseAlphaMixImg, "squared_pdbiodosealphamix"); + setupImage(squaredPdBioDoseSqrtBetaMixImg, "squared_pdbiodosesqrtbetamix"); + setupImage(squaredPdBioDoseDoseImg, "squared_pdbiodosedose"); + setupImage(varAlphaMixImg, "var_alphamix"); + setupImage(varSqrtBetaMixImg, "var_sqrtbetamix"); + setupImage(varDoseImg, "var_dose"); + /* ******************* */ + + auto const sqAlphaRef = _alphaRef * _alphaRef; + double const n = _currentEvent; + + for(auto const& [index, deposited]: _depositedMap) { + auto const eventCount = _eventCountImage.GetValue(index); + + auto const alphaMixMean = _alphaMixImage.GetValue(index) / eventCount; + auto const sqrtBetaMixMean = _sqrtBetaMixImage.GetValue(index) / eventCount; + auto const dose = deposited.dose; + auto const scaledDose = _doseScaleFactor * dose; + auto const sqScaledDose = scaledDose * scaledDose; + auto const delta = sqAlphaRef + 4 * _betaRef * + (alphaMixMean * scaledDose + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDose); + + double sqrtDelta = 0; + if(delta >= 0) + sqrtDelta = std::sqrt(delta); // Calculate biological dose and RBE double biodose = 0; double rbe = 0; - if(deposited.dose > 0 && alphamix != 0 && betamix != 0) { - auto const sqAlphaRef = _alphaRef * _alphaRef; - auto const sqDose = deposited.dose * deposited.dose; - biodose = (-_alphaRef + std::sqrt(sqAlphaRef + 4 * _betaRef * (alphamix * deposited.dose + betamix * sqDose))) / (2 * _betaRef); + if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0) + biodose = (-_alphaRef + sqrtDelta) / (2 * _betaRef); + if(biodose < 0) biodose = 0; // TODO improve + + if(scaledDose > 0) + rbe = biodose / scaledDose; + + if(index == 80) { // TODO remove + std::ofstream of{"/tmp/eoea", std::ios_base::app}; + of << "alphaRef: " << _alphaRef << '\n'; + of << "sqrtBetaRef: " << _betaRef << '\n'; + of << "s: " << _doseScaleFactor << '\n'; + of << "eventCount: " << eventCount << '\n'; + of << "alphaMixMean: " << _alphaMixImage.GetValue(index) << " / " << eventCount << " = " << alphaMixMean << '\n'; + of << "sqrtBetaMixMean: " << sqrtBetaMixMean << '\n'; + of << "delta: " << delta << '\n'; + of << "final sum dose: " << dose << '\n'; + of << "final sum squared dose: " << _squaredDoseImage.GetValue(index) << '\n'; + of << "biodose: " << biodose << '\n'; + of << "rbe: " << rbe << '\n'; + } - rbe = biodose / deposited.dose; + if(_enableUncertainties) { + if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0 && sqrtDelta > 0 && _currentEvent > 0) { + double const doseMean = dose / n; + + double sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); + double sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); + double sumSquaredDose = _squaredDoseImage.GetValue(index); + + double pdBiodoseAlphaMix = scaledDose / sqrtDelta; + double pdBiodoseSqrtBetaMix = 2 * sqScaledDose * sqrtBetaMixMean / sqrtDelta; + double pdBiodoseDose = ( + (alphaMixMean * _doseScaleFactor) + + 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose + ) / sqrtDelta; + + double varAlphaMix = sumSquaredAlphaMix / eventCount - alphaMixMean * alphaMixMean; + double varSqrtBetaMix = sumSquaredSqrtBetaMix / eventCount - sqrtBetaMixMean * sqrtBetaMixMean; + double varDose = sumSquaredDose / n - doseMean * doseMean; + + double sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); + double sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); + double sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); + double covAlphaMixSqrtBetaMix = sumAlphaMixSqrtBetaMix / eventCount - alphaMixMean * sqrtBetaMixMean; + double covAlphaMixDose = sumAlphaMixDose / n - alphaMixMean * doseMean; + double covSqrtBetaMixDose = sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean; + + double partAlphaMix = pdBiodoseAlphaMix * pdBiodoseAlphaMix * varAlphaMix; + double partSqrtBetaMix = pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix * varSqrtBetaMix; + double partDose = pdBiodoseDose * pdBiodoseDose * varDose; + double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMix * pdBiodoseSqrtBetaMix * covAlphaMixSqrtBetaMix; + double partAlphaMixDose = 2 * pdBiodoseAlphaMix * pdBiodoseDose * covAlphaMixDose; + double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMix * pdBiodoseDose * covSqrtBetaMixDose; + + { // TODO remove + partAlphaMixImg.SetValue(index, partAlphaMix); + partSqrtBetaMixImg.SetValue(index, partSqrtBetaMix); + partDoseImg.SetValue(index, partDose); + partAlphaMixDoseImg.SetValue(index, partAlphaMixDose); + partAlphaMixSqrtBetaMixImg.SetValue(index, partAlphaMixSqrtBetaMix); + partSqrtBetaMixDoseImg.SetValue(index, partSqrtBetaMixDose); + partSquaredBioDoseImg.SetValue(index, biodose * biodose); + + squaredPdBioDoseAlphaMixImg.SetValue(index, pdBiodoseAlphaMix * pdBiodoseAlphaMix); + squaredPdBioDoseSqrtBetaMixImg.SetValue(index, pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix); + squaredPdBioDoseDoseImg.SetValue(index, pdBiodoseDose * pdBiodoseDose); + varAlphaMixImg.SetValue(index, varAlphaMix); + varSqrtBetaMixImg.SetValue(index, varSqrtBetaMix); + varDoseImg.SetValue(index, varDose); + } + + double uncertaintyBiodose = std::sqrt( + partAlphaMix + partSqrtBetaMix + partDose + + partAlphaMixSqrtBetaMix + partAlphaMixDose + partSqrtBetaMixDose + ) / biodose; + + _biodoseUncertaintyImage.SetValue(index, uncertaintyBiodose); + } else { + _biodoseUncertaintyImage.SetValue(index, 1); + } } // Write data - if(_enableEdep) _edepImage.SetValue(index, deposited.energy); - if(_enableDose) _doseImage.SetValue(index, deposited.dose); - if(_enableBioDose) _bioDoseImage.SetValue(index, biodose); - if(_enableAlphaMix) _alphaMixImage.SetValue(index, alphamix); - if(_enableBetaMix) _betaMixImage.SetValue(index, betamix); - if(_enableRBE) _RBEImage.SetValue(index, rbe); + if(_enableEdep) _edepImage.SetValue(index, deposited.energy); + if(_enableDose) _doseImage.SetValue(index, scaledDose); + _bioDoseImage.SetValue(index, biodose); + if(_enableRBE) _rbeImage.SetValue(index, rbe); } - //------------------------------------------------------------- + GateVActor::SaveData(); - if(_enableEdep) _edepImage.SaveData(_currentEvent); - if(_enableDose) _doseImage.SaveData(_currentEvent); - if(_enableBioDose) _bioDoseImage.SaveData(_currentEvent); - if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); - if(_enableBetaMix) _betaMixImage.SaveData(_currentEvent); - if(_enableRBE) _RBEImage.SaveData(_currentEvent); + if(_enableEdep) _edepImage.SaveData(_currentEvent); + if(_enableDose) _doseImage.SaveData(_currentEvent); + if(_enableBioDose) _bioDoseImage.SaveData(_currentEvent); + if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); + if(_enableSqrtBetaMix) _sqrtBetaMixImage.SaveData(_currentEvent); + if(_enableRBE) _rbeImage.SaveData(_currentEvent); + if(_enableUncertainties) _biodoseUncertaintyImage.SaveData(_currentEvent); + + { // TODO remove + partAlphaMixImg.SaveData(_currentEvent); + partSqrtBetaMixImg.SaveData(_currentEvent); + partDoseImg.SaveData(_currentEvent); + partAlphaMixDoseImg.SaveData(_currentEvent); + partAlphaMixSqrtBetaMixImg.SaveData(_currentEvent); + partSqrtBetaMixDoseImg.SaveData(_currentEvent); + partSquaredBioDoseImg.SaveData(_currentEvent); + + squaredPdBioDoseAlphaMixImg.SaveData(_currentEvent); + squaredPdBioDoseSqrtBetaMixImg.SaveData(_currentEvent); + squaredPdBioDoseDoseImg.SaveData(_currentEvent); + varAlphaMixImg.SaveData(_currentEvent); + varSqrtBetaMixImg.SaveData(_currentEvent); + varDoseImg.SaveData(_currentEvent); + } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActor::BeginOfRunAction(const G4Run* r) { GateVActor::BeginOfRunAction(r); GateDebugMessage("Actor", 3, "GateBioDoseActor -- Begin of Run\n"); + + ResetData(); } //----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- // PAS DANS DOSE ACTOR +//----------------------------------------------------------------------------- void GateBioDoseActor::EndOfRunAction(const G4Run* r) { GateVActor::EndOfRunAction(r); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -// Callback at each event void GateBioDoseActor::BeginOfEventAction(const G4Event* e) { GateVActor::BeginOfEventAction(e); ++_currentEvent; + + _eventVoxelIndices.clear(); + + _eventEdepImage.Reset(); + _eventDoseImage.Reset(); + _eventAlphaImage.Reset(); + _eventSqrtBetaImage.Reset(); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::EndOfEventAction(const G4Event* e) { + GateVActor::EndOfEventAction(e); + + for(auto const& index: _eventVoxelIndices) { + auto const eventEdep = _eventEdepImage.GetValue(index); + auto const eventDose = _eventDoseImage.GetValue(index); + auto const eventAlphaMix = _eventAlphaImage.GetValue(index) / eventEdep; + auto const eventSqrtBetaMix = _eventSqrtBetaImage.GetValue(index) / eventEdep; + + _eventCountImage.AddValue(index, 1); + + _alphaMixImage.AddValue(index, eventAlphaMix); + _sqrtBetaMixImage.AddValue(index, eventSqrtBetaMix); + + if(_enableUncertainties) { + _squaredDoseImage.AddValue(index, eventDose * eventDose); + _squaredAlphaMixImage.AddValue(index, eventAlphaMix * eventAlphaMix); + _squaredSqrtBetaMixImage.AddValue(index, eventSqrtBetaMix * eventSqrtBetaMix); + + _alphaMixSqrtBetaMixImage.AddValue(index, eventAlphaMix * eventSqrtBetaMix); + _alphaMixDoseImage.AddValue(index, eventAlphaMix * eventDose); + _sqrtBetaMixDoseImage.AddValue(index, eventSqrtBetaMix * eventDose); + } + + if(index == 80) { // TODO remove + std::ofstream of{"/tmp/alpha", std::ios_base::app}; + of << "EVENT alpha: " << _eventAlphaImage.GetValue(index) << '\n'; + of << "EVENT sqrtBeta: " << _eventSqrtBetaImage.GetValue(index) << '\n'; + of << "EVENT energyDep: " << eventEdep << '\n'; + of << "--------------------\n"; + } + + std::ofstream of{"/tmp/check", std::ios_base::app}; + of << "index: " << index << '\n'; + if(index == 80 || index == 320) { + std::ofstream of{"/tmp/eoea", std::ios_base::app}; + of << "index: " << index << '\n'; + of << "eventEdep: " << eventEdep << '\n'; + of << "eventDose: " << eventDose << '\n'; + of << "eventDose²: " << eventDose * eventDose << '\n'; + of << "current sum dose: " << _depositedMap.at(index).dose << '\n'; + of << "current sum squared dose: " << _squaredDoseImage.GetValue(index) << '\n'; + of << "------------------------\n"; + } + } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -215,7 +444,7 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* if(energyDep == 0) return; if(index < 0) return; - DepositedMap::iterator it = _depositedMap.find(index); + auto it = _depositedMap.find(index); if(it == std::end(_depositedMap)) { _depositedMap[index] = {0, 0, 0, 0}; it = _depositedMap.find(index); @@ -227,71 +456,93 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* deposited.energy += energyDep; if(_enableDose || _enableBioDose || _enableRBE) { - decltype(_doseImage)* image = nullptr; - if(_enableDose) image = &_doseImage; - else if(_enableBioDose) image = &_bioDoseImage; - else if(_enableRBE) image = &_RBEImage; - - auto currentMaterial = step->GetPreStepPoint()->GetMaterial(); + auto* currentMaterial = step->GetPreStepPoint()->GetMaterial(); double density = currentMaterial->GetDensity(); - double mass = image->GetVoxelVolume() * density; + double mass = _bioDoseImage.GetVoxelVolume() * density; + double dose = energyDep / mass / CLHEP::gray; - deposited.dose += energyDep / mass / CLHEP::gray; + deposited.dose += dose; + _eventDoseImage.AddValue(index, dose); } // Get information from step // Particle G4int nZ = step->GetTrack()->GetDefinition()->GetAtomicNumber(); - double kineticEnergyPerNucleon = (step->GetPreStepPoint()->GetKineticEnergy()) / (step->GetTrack()->GetDefinition()->GetAtomicMass()); //OK + double kineticEnergyPerNucleon = (step->GetPreStepPoint()->GetKineticEnergy()) / (step->GetTrack()->GetDefinition()->GetAtomicMass()); - ++_eventCount; + ++_stepCount; // Accumulation of alpha/beta if ion type if known // -> check if the ion type is known - if(_energyMaxForZ.count(nZ)) { - ++_eventWithKnownIonCount; + if(_energyMaxForZ.count(nZ) != 0) { + ++_stepWithKnownIonCount; - //The max values in the database aren't being taking into account - //so for now it's coded like this to be sure the code takes them into account double energyMax = _energyMaxForZ.at(nZ); - AlphaBetaInterpolTable::const_iterator itr2; - if (kineticEnergyPerNucleon >= energyMax) { - // If the kinetic energy is the maximum value in the alpha beta tables, - // we have to use the a and b coefficient for this maximum value + AlphaBetaInterpolTable::const_iterator itInterpol; + if(kineticEnergyPerNucleon >= energyMax) { Fragment fragmentKineticEnergyMax{nZ, energyMax}; - itr2 = _alphaBetaInterpolTable.find(fragmentKineticEnergyMax); + itInterpol = _alphaBetaInterpolTable.find(fragmentKineticEnergyMax); } else { - // We pair the ion type and the kinetic energy Fragment fragmentKineticEnergy{nZ, kineticEnergyPerNucleon}; - itr2 = _alphaBetaInterpolTable.upper_bound(fragmentKineticEnergy); + itInterpol = _alphaBetaInterpolTable.upper_bound(fragmentKineticEnergy); } - // Calculation of EZ, alphaDep and betaDep (K = a*EZ+b*E) - auto const& interpol = (*itr2).second; + // Calculation of alphaDep and betaDep (K = (a*Z+b)*E) + auto const& interpol = (*itInterpol).second; + + double alpha = (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b) * energyDep; + double sqrtBeta = (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b) * energyDep; + + if(index == 80) { // TODO remove + std::ofstream of{"/tmp/alpha", std::ios_base::app}; + of << "alpha: " << (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b) << '\n'; + of << "sqrtBeta: " << (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b) << '\n'; + of << "energyDep: " << energyDep << '\n'; + of << "--------------------\n"; + } - double alpha = (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b); - double sqrtBeta = (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b); - if(alpha < 0) alpha = 0; if(sqrtBeta < 0) sqrtBeta = 0; - double alphaDep = alpha * energyDep; - double sqrtBetaDep = sqrtBeta * energyDep; - // Accumulate alpha/beta - deposited.alpha += alphaDep; - deposited.sqrtBeta += sqrtBetaDep; + deposited.alpha += alpha; + deposited.sqrtBeta += sqrtBeta; + + _eventEdepImage.AddValue(index, energyDep); + _eventAlphaImage.AddValue(index, alpha); + _eventSqrtBetaImage.AddValue(index, sqrtBeta); + + _eventVoxelIndices.insert(index); } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActor::ResetData() { - if(_enableEdep) _edepImage.Reset(); - if(_enableDose) _doseImage.Reset(); - if(_enableBioDose) _bioDoseImage.Reset(); - if(_enableAlphaMix) _alphaMixImage.Reset(); - if(_enableBetaMix) _betaMixImage.Reset(); - if(_enableRBE) _RBEImage.Reset(); + _eventCountImage.Reset(); + + if(_enableEdep) _edepImage.Reset(); + if(_enableDose) _doseImage.Reset(); + _bioDoseImage.Reset(); + _alphaMixImage.Reset(); + _sqrtBetaMixImage.Reset(); + if(_enableRBE) _rbeImage.Reset(); + + _eventEdepImage.Reset(); + _eventDoseImage.Reset(); + _eventAlphaImage.Reset(); + _eventSqrtBetaImage.Reset(); + + if(_enableUncertainties) { + _biodoseUncertaintyImage.Reset(); + + _squaredDoseImage.Reset(); + _squaredAlphaMixImage.Reset(); + _squaredSqrtBetaMixImage.Reset(); + + _alphaMixSqrtBetaMixImage.Reset(); + _alphaMixDoseImage.Reset(); + _sqrtBetaMixDoseImage.Reset(); + } } //----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc index dc3822514..eb1b335ee 100644 --- a/source/digits_hits/src/GateBioDoseActorMessenger.cc +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -23,6 +23,10 @@ GateBioDoseActorMessenger::GateBioDoseActorMessenger(GateBioDoseActor* sensor): void GateBioDoseActorMessenger::BuildCommands(G4String base) { G4String n; + n = base + "/setDoseScaleFactor"; + pDoseScaleFactorCmd = std::make_unique(n, this); + pDoseScaleFactorCmd->SetGuidance("Set (physical) dose scale factor (default: 1.0)"); + n = base + "/setAlphaRef"; pAlphaRefCmd = std::make_unique(n, this); pAlphaRefCmd->SetGuidance("See [...] for values from publications"); @@ -60,18 +64,23 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { pEnableAlphaMixCmd = std::make_unique(n, this); pEnableAlphaMixCmd->SetGuidance("Enable alpha mix output"); - n = base + "/enableBetaMix"; - pEnableBetaMixCmd = std::make_unique(n, this); - pEnableBetaMixCmd->SetGuidance("Enable beta mix output"); + n = base + "/enableSqrtBetaMix"; + pEnableSqrtBetaMixCmd = std::make_unique(n, this); + pEnableSqrtBetaMixCmd->SetGuidance("Enable sqrt(beta) mix output"); n = base + "/enableRBE"; pEnableRBECmd = std::make_unique(n, this); pEnableRBECmd->SetGuidance("Enable RBE output"); + + n = base + "/enableUncertainty"; + pEnableUncertaintyCmd = std::make_unique(n, this); + pEnableUncertaintyCmd->SetGuidance("Enable uncertainty outputs"); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { - if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(value)); + if(cmd == pDoseScaleFactorCmd.get()) pBioDoseActor->SetDoseScaleFactor(pDoseScaleFactorCmd->GetNewDoubleValue(value)); + else if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(value)); else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(value)); else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(value); else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(value); @@ -80,8 +89,9 @@ void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); else if(cmd == pEnableBioDoseCmd.get()) pBioDoseActor->SetEnableBioDose(pEnableBioDoseCmd->GetNewBoolValue(value)); else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); - else if(cmd == pEnableBetaMixCmd.get()) pBioDoseActor->SetEnableBetaMix(pEnableBetaMixCmd->GetNewBoolValue(value)); + else if(cmd == pEnableSqrtBetaMixCmd.get()) pBioDoseActor->SetEnableSqrtBetaMix(pEnableSqrtBetaMixCmd->GetNewBoolValue(value)); else if(cmd == pEnableRBECmd.get()) pBioDoseActor->SetEnableRBE(pEnableRBECmd->GetNewBoolValue(value)); + else if(cmd == pEnableUncertaintyCmd.get()) pBioDoseActor->SetEnableUncertainties(pEnableUncertaintyCmd->GetNewBoolValue(value)); GateImageActorMessenger::SetNewValue(cmd, value); } From 6b26603117c245b8a03f55ba83b7f613ef2a0b54 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Wed, 14 Feb 2024 18:08:55 +0100 Subject: [PATCH 08/12] Prepare for release, cleanup code --- .../digits_hits/include/GateBioDoseActor.hh | 36 +- source/digits_hits/src/GateBioDoseActor.cc | 402 ++++++------------ .../src/GateBioDoseActorMessenger.cc | 5 - 3 files changed, 146 insertions(+), 297 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index c3ed226c1..3d3365021 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -23,13 +23,6 @@ See LICENSE.md for further details class GateBioDoseActor: public GateVImageActor { public: - struct Deposited { - double alpha; - double sqrtBeta; - double energy; - double dose; - }; - struct Coefficients { double a, b; }; @@ -40,7 +33,6 @@ public: }; using VoxelIndex = int; - using DepositedMap = std::map; using VoxelIndices = std::set; using Fragment = std::pair; @@ -84,19 +76,18 @@ public: void SetEnableEdep(bool e) { _enableEdep = e; } void SetEnableDose(bool e) { _enableDose = e; } - void SetEnableBioDose(bool e) { _enableBioDose = e; } void SetEnableAlphaMix(bool e) { _enableAlphaMix = e; } void SetEnableSqrtBetaMix(bool e) { _enableSqrtBetaMix = e; } void SetEnableRBE(bool e) { _enableRBE = e; } void SetEnableUncertainties(bool e) { _enableUncertainties = e; } - // Input database - void BuildDatabase(); - static Coefficients Interpol(double x1, double x2, double y1, double y2); - protected: GateBioDoseActor(G4String name, G4int depth = 0); + void updateData(); + void buildDatabase(); + static Coefficients interpol(double x1, double x2, double y1, double y2); + private: //Counters int _currentEvent; @@ -114,31 +105,31 @@ private: G4double _sobpWeight; - // Maps - DepositedMap _depositedMap; AlphaBetaInterpolTable _alphaBetaInterpolTable; VoxelIndices _eventVoxelIndices; + VoxelIndices _voxelIndices; // Images - GateImageWithStatistic _eventCountImage; + GateImageWithStatistic _hitEventCountImage; + + GateImageWithStatistic _eventEdepImage; + GateImageWithStatistic _eventDoseImage; + GateImageWithStatistic _eventAlphaImage; + GateImageWithStatistic _eventSqrtBetaImage; - GateImageWithStatistic _bioDoseImage; GateImageWithStatistic _edepImage; GateImageWithStatistic _doseImage; + GateImageWithStatistic _scaledDoseImage; GateImageWithStatistic _alphaMixImage; GateImageWithStatistic _sqrtBetaMixImage; + GateImageWithStatistic _bioDoseImage; GateImageWithStatistic _rbeImage; GateImageWithStatistic _biodoseUncertaintyImage; - GateImageWithStatistic _eventEdepImage; - GateImageWithStatistic _eventDoseImage; GateImageWithStatistic _squaredDoseImage; - GateImageWithStatistic _eventAlphaImage; GateImageWithStatistic _squaredAlphaMixImage; - GateImageWithStatistic _eventSqrtBetaImage; GateImageWithStatistic _squaredSqrtBetaMixImage; - GateImageWithStatistic _alphaMixSqrtBetaMixImage; GateImageWithStatistic _alphaMixDoseImage; GateImageWithStatistic _sqrtBetaMixDoseImage; @@ -146,7 +137,6 @@ private: // Outputs bool _enableEdep; bool _enableDose; - bool _enableBioDose; bool _enableAlphaMix; bool _enableSqrtBetaMix; bool _enableRBE; diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index 8889dbef1..f496fbff0 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -8,7 +8,6 @@ See LICENSE.md for further details #include "G4EmParameters.hh" #include "GateBioDoseActor.hh" -#include "GateImage.hh" #include "GateImageWithStatistic.hh" #include @@ -24,7 +23,6 @@ GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): _sobpWeight(0), _enableEdep(false), _enableDose(true), - _enableBioDose(true), _enableAlphaMix(false), _enableSqrtBetaMix(false), _enableRBE(false), @@ -57,7 +55,6 @@ void GateBioDoseActor::Construct() { { G4String basename = removeExtension(mSaveFilename); G4String ext = getExtension(mSaveFilename); - auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { SetOriginTransformAndFlagToImage(image); image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); @@ -69,27 +66,26 @@ void GateBioDoseActor::Construct() { } }; - setupImage(_eventCountImage); - - if(_enableEdep) setupImage(_edepImage, "edep"); - if(_enableDose) setupImage(_doseImage, "dose"); - setupImage(_bioDoseImage, "biodose"); - setupImage(_alphaMixImage, "alphamix"); - setupImage(_sqrtBetaMixImage, "sqrtbetamix"); - if(_enableRBE) setupImage(_rbeImage, "rbe"); + setupImage(_hitEventCountImage); setupImage(_eventEdepImage); setupImage(_eventDoseImage); setupImage(_eventAlphaImage); setupImage(_eventSqrtBetaImage); + if(_enableEdep) setupImage(_edepImage, "edep"); + setupImage(_doseImage); // dose output can be scaled, see SaveData() + if(_enableDose) setupImage(_scaledDoseImage, "dose"); + setupImage(_alphaMixImage, "alphamix"); + setupImage(_sqrtBetaMixImage, "sqrtbetamix"); + setupImage(_bioDoseImage, "biodose"); + if(_enableRBE) setupImage(_rbeImage, "rbe"); + if(_enableUncertainties) { setupImage(_biodoseUncertaintyImage, "biodose_uncertainty"); - setupImage(_squaredDoseImage); setupImage(_squaredAlphaMixImage); setupImage(_squaredSqrtBetaMixImage); - setupImage(_alphaMixSqrtBetaMixImage); setupImage(_alphaMixDoseImage); setupImage(_sqrtBetaMixDoseImage); @@ -107,14 +103,14 @@ void GateBioDoseActor::Construct() { //Building the cell line information _dataBase = "data/" + _cellLine + "_" + _bioPhysicalModel + ".db"; - BuildDatabase(); + buildDatabase(); if(_alphaRef < 0 || _betaRef < 0) GateError("BioDoseActor " << GetName() << ": setAlphaRef and setBetaRef must be done"); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -void GateBioDoseActor::BuildDatabase() { +void GateBioDoseActor::buildDatabase() { std::ifstream f(_dataBase); if(!f) GateError("BioDoseActor " << GetName() << ": unable to open file '" << _dataBase << "'"); @@ -145,8 +141,8 @@ void GateBioDoseActor::BuildDatabase() { iss >> alpha; iss >> beta; - auto alphaCoeff = Interpol(prevKineticEnergy, kineticEnergy, prevAlpha, alpha); - auto sqrtBetaCoeff = Interpol(prevKineticEnergy, kineticEnergy, std::sqrt(prevBeta), std::sqrt(beta)); + auto alphaCoeff = interpol(prevKineticEnergy, kineticEnergy, prevAlpha, alpha); + auto sqrtBetaCoeff = interpol(prevKineticEnergy, kineticEnergy, std::sqrt(prevBeta), std::sqrt(beta)); // Saving the in the input databse Fragment fragment{nZ, kineticEnergy}; @@ -165,7 +161,7 @@ void GateBioDoseActor::BuildDatabase() { } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, double y1, double y2) { +GateBioDoseActor::Coefficients GateBioDoseActor::interpol(double x1, double x2, double y1, double y2) { //Function for a 1D linear interpolation. It returns a pair of a and b coefficients double a = (y2 - y1) / (x2 - x1); double b = y1 - x1 * a; @@ -173,196 +169,6 @@ GateBioDoseActor::Coefficients GateBioDoseActor::Interpol(double x1, double x2, } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -void GateBioDoseActor::SaveData() { - GateDebugMessageInc("Actor", 4, "GateBioDoseActor::SaveData() known ion events / total events: " << _eventWithKnownIonCount << " / " << _eventCount << "\n"); - - /* ******************* */ - // TODO remove - G4String basename = removeExtension(mSaveFilename); - G4String ext = getExtension(mSaveFilename); - auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { - SetOriginTransformAndFlagToImage(image); - image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - image.Allocate(); - - if(!suffix.empty()) { - G4String filename = basename + "_" + suffix + "." + ext; - image.SetFilename(filename); - } - }; - - GateImageWithStatistic partAlphaMixImg; - GateImageWithStatistic partSqrtBetaMixImg; - GateImageWithStatistic partDoseImg; - GateImageWithStatistic partAlphaMixSqrtBetaMixImg; - GateImageWithStatistic partAlphaMixDoseImg; - GateImageWithStatistic partSqrtBetaMixDoseImg; - GateImageWithStatistic partSquaredBioDoseImg; - - GateImageWithStatistic squaredPdBioDoseAlphaMixImg; - GateImageWithStatistic squaredPdBioDoseSqrtBetaMixImg; - GateImageWithStatistic squaredPdBioDoseDoseImg; - GateImageWithStatistic varAlphaMixImg; - GateImageWithStatistic varSqrtBetaMixImg; - GateImageWithStatistic varDoseImg; - - setupImage(partAlphaMixImg, "part_alphamix"); - setupImage(partSqrtBetaMixImg, "part_sqrtbetamix"); - setupImage(partDoseImg, "part_dose"); - setupImage(partAlphaMixSqrtBetaMixImg, "part_alphamixsqrtbetamix"); - setupImage(partAlphaMixDoseImg, "part_alphamixdose"); - setupImage(partSqrtBetaMixDoseImg, "part_sqrtbetamixdose"); - setupImage(partSquaredBioDoseImg, "part_squaredbiodose"); - - setupImage(squaredPdBioDoseAlphaMixImg, "squared_pdbiodosealphamix"); - setupImage(squaredPdBioDoseSqrtBetaMixImg, "squared_pdbiodosesqrtbetamix"); - setupImage(squaredPdBioDoseDoseImg, "squared_pdbiodosedose"); - setupImage(varAlphaMixImg, "var_alphamix"); - setupImage(varSqrtBetaMixImg, "var_sqrtbetamix"); - setupImage(varDoseImg, "var_dose"); - /* ******************* */ - - auto const sqAlphaRef = _alphaRef * _alphaRef; - double const n = _currentEvent; - - for(auto const& [index, deposited]: _depositedMap) { - auto const eventCount = _eventCountImage.GetValue(index); - - auto const alphaMixMean = _alphaMixImage.GetValue(index) / eventCount; - auto const sqrtBetaMixMean = _sqrtBetaMixImage.GetValue(index) / eventCount; - auto const dose = deposited.dose; - auto const scaledDose = _doseScaleFactor * dose; - auto const sqScaledDose = scaledDose * scaledDose; - auto const delta = sqAlphaRef + 4 * _betaRef * - (alphaMixMean * scaledDose + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDose); - - double sqrtDelta = 0; - if(delta >= 0) - sqrtDelta = std::sqrt(delta); - - // Calculate biological dose and RBE - double biodose = 0; - double rbe = 0; - - if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0) - biodose = (-_alphaRef + sqrtDelta) / (2 * _betaRef); - if(biodose < 0) biodose = 0; // TODO improve - - if(scaledDose > 0) - rbe = biodose / scaledDose; - - if(index == 80) { // TODO remove - std::ofstream of{"/tmp/eoea", std::ios_base::app}; - of << "alphaRef: " << _alphaRef << '\n'; - of << "sqrtBetaRef: " << _betaRef << '\n'; - of << "s: " << _doseScaleFactor << '\n'; - of << "eventCount: " << eventCount << '\n'; - of << "alphaMixMean: " << _alphaMixImage.GetValue(index) << " / " << eventCount << " = " << alphaMixMean << '\n'; - of << "sqrtBetaMixMean: " << sqrtBetaMixMean << '\n'; - of << "delta: " << delta << '\n'; - of << "final sum dose: " << dose << '\n'; - of << "final sum squared dose: " << _squaredDoseImage.GetValue(index) << '\n'; - of << "biodose: " << biodose << '\n'; - of << "rbe: " << rbe << '\n'; - } - - if(_enableUncertainties) { - if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0 && sqrtDelta > 0 && _currentEvent > 0) { - double const doseMean = dose / n; - - double sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); - double sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); - double sumSquaredDose = _squaredDoseImage.GetValue(index); - - double pdBiodoseAlphaMix = scaledDose / sqrtDelta; - double pdBiodoseSqrtBetaMix = 2 * sqScaledDose * sqrtBetaMixMean / sqrtDelta; - double pdBiodoseDose = ( - (alphaMixMean * _doseScaleFactor) + - 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose - ) / sqrtDelta; - - double varAlphaMix = sumSquaredAlphaMix / eventCount - alphaMixMean * alphaMixMean; - double varSqrtBetaMix = sumSquaredSqrtBetaMix / eventCount - sqrtBetaMixMean * sqrtBetaMixMean; - double varDose = sumSquaredDose / n - doseMean * doseMean; - - double sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); - double sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); - double sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); - double covAlphaMixSqrtBetaMix = sumAlphaMixSqrtBetaMix / eventCount - alphaMixMean * sqrtBetaMixMean; - double covAlphaMixDose = sumAlphaMixDose / n - alphaMixMean * doseMean; - double covSqrtBetaMixDose = sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean; - - double partAlphaMix = pdBiodoseAlphaMix * pdBiodoseAlphaMix * varAlphaMix; - double partSqrtBetaMix = pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix * varSqrtBetaMix; - double partDose = pdBiodoseDose * pdBiodoseDose * varDose; - double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMix * pdBiodoseSqrtBetaMix * covAlphaMixSqrtBetaMix; - double partAlphaMixDose = 2 * pdBiodoseAlphaMix * pdBiodoseDose * covAlphaMixDose; - double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMix * pdBiodoseDose * covSqrtBetaMixDose; - - { // TODO remove - partAlphaMixImg.SetValue(index, partAlphaMix); - partSqrtBetaMixImg.SetValue(index, partSqrtBetaMix); - partDoseImg.SetValue(index, partDose); - partAlphaMixDoseImg.SetValue(index, partAlphaMixDose); - partAlphaMixSqrtBetaMixImg.SetValue(index, partAlphaMixSqrtBetaMix); - partSqrtBetaMixDoseImg.SetValue(index, partSqrtBetaMixDose); - partSquaredBioDoseImg.SetValue(index, biodose * biodose); - - squaredPdBioDoseAlphaMixImg.SetValue(index, pdBiodoseAlphaMix * pdBiodoseAlphaMix); - squaredPdBioDoseSqrtBetaMixImg.SetValue(index, pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix); - squaredPdBioDoseDoseImg.SetValue(index, pdBiodoseDose * pdBiodoseDose); - varAlphaMixImg.SetValue(index, varAlphaMix); - varSqrtBetaMixImg.SetValue(index, varSqrtBetaMix); - varDoseImg.SetValue(index, varDose); - } - - double uncertaintyBiodose = std::sqrt( - partAlphaMix + partSqrtBetaMix + partDose + - partAlphaMixSqrtBetaMix + partAlphaMixDose + partSqrtBetaMixDose - ) / biodose; - - _biodoseUncertaintyImage.SetValue(index, uncertaintyBiodose); - } else { - _biodoseUncertaintyImage.SetValue(index, 1); - } - } - - // Write data - if(_enableEdep) _edepImage.SetValue(index, deposited.energy); - if(_enableDose) _doseImage.SetValue(index, scaledDose); - _bioDoseImage.SetValue(index, biodose); - if(_enableRBE) _rbeImage.SetValue(index, rbe); - } - - GateVActor::SaveData(); - - if(_enableEdep) _edepImage.SaveData(_currentEvent); - if(_enableDose) _doseImage.SaveData(_currentEvent); - if(_enableBioDose) _bioDoseImage.SaveData(_currentEvent); - if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); - if(_enableSqrtBetaMix) _sqrtBetaMixImage.SaveData(_currentEvent); - if(_enableRBE) _rbeImage.SaveData(_currentEvent); - if(_enableUncertainties) _biodoseUncertaintyImage.SaveData(_currentEvent); - - { // TODO remove - partAlphaMixImg.SaveData(_currentEvent); - partSqrtBetaMixImg.SaveData(_currentEvent); - partDoseImg.SaveData(_currentEvent); - partAlphaMixDoseImg.SaveData(_currentEvent); - partAlphaMixSqrtBetaMixImg.SaveData(_currentEvent); - partSqrtBetaMixDoseImg.SaveData(_currentEvent); - partSquaredBioDoseImg.SaveData(_currentEvent); - - squaredPdBioDoseAlphaMixImg.SaveData(_currentEvent); - squaredPdBioDoseSqrtBetaMixImg.SaveData(_currentEvent); - squaredPdBioDoseDoseImg.SaveData(_currentEvent); - varAlphaMixImg.SaveData(_currentEvent); - varSqrtBetaMixImg.SaveData(_currentEvent); - varDoseImg.SaveData(_currentEvent); - } -} -//----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- void GateBioDoseActor::BeginOfRunAction(const G4Run* r) { GateVActor::BeginOfRunAction(r); GateDebugMessage("Actor", 3, "GateBioDoseActor -- Begin of Run\n"); @@ -398,8 +204,12 @@ void GateBioDoseActor::EndOfEventAction(const G4Event* e) { auto const eventAlphaMix = _eventAlphaImage.GetValue(index) / eventEdep; auto const eventSqrtBetaMix = _eventSqrtBetaImage.GetValue(index) / eventEdep; - _eventCountImage.AddValue(index, 1); + _voxelIndices.insert(index); + _hitEventCountImage.AddValue(index, 1); + + if(_enableEdep) _edepImage.AddValue(index, eventEdep); + _doseImage.AddValue(index, eventDose); _alphaMixImage.AddValue(index, eventAlphaMix); _sqrtBetaMixImage.AddValue(index, eventSqrtBetaMix); @@ -412,27 +222,6 @@ void GateBioDoseActor::EndOfEventAction(const G4Event* e) { _alphaMixDoseImage.AddValue(index, eventAlphaMix * eventDose); _sqrtBetaMixDoseImage.AddValue(index, eventSqrtBetaMix * eventDose); } - - if(index == 80) { // TODO remove - std::ofstream of{"/tmp/alpha", std::ios_base::app}; - of << "EVENT alpha: " << _eventAlphaImage.GetValue(index) << '\n'; - of << "EVENT sqrtBeta: " << _eventSqrtBetaImage.GetValue(index) << '\n'; - of << "EVENT energyDep: " << eventEdep << '\n'; - of << "--------------------\n"; - } - - std::ofstream of{"/tmp/check", std::ios_base::app}; - of << "index: " << index << '\n'; - if(index == 80 || index == 320) { - std::ofstream of{"/tmp/eoea", std::ios_base::app}; - of << "index: " << index << '\n'; - of << "eventEdep: " << eventEdep << '\n'; - of << "eventDose: " << eventDose << '\n'; - of << "eventDose²: " << eventDose * eventDose << '\n'; - of << "current sum dose: " << _depositedMap.at(index).dose << '\n'; - of << "current sum squared dose: " << _squaredDoseImage.GetValue(index) << '\n'; - of << "------------------------\n"; - } } } //----------------------------------------------------------------------------- @@ -444,26 +233,15 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* if(energyDep == 0) return; if(index < 0) return; - auto it = _depositedMap.find(index); - if(it == std::end(_depositedMap)) { - _depositedMap[index] = {0, 0, 0, 0}; - it = _depositedMap.find(index); - } - - auto& deposited = (*it).second; - // Accumulate energy inconditionnaly - deposited.energy += energyDep; + _eventEdepImage.AddValue(index, energyDep); - if(_enableDose || _enableBioDose || _enableRBE) { - auto* currentMaterial = step->GetPreStepPoint()->GetMaterial(); - double density = currentMaterial->GetDensity(); - double mass = _bioDoseImage.GetVoxelVolume() * density; - double dose = energyDep / mass / CLHEP::gray; + auto* currentMaterial = step->GetPreStepPoint()->GetMaterial(); + double density = currentMaterial->GetDensity(); + double mass = _bioDoseImage.GetVoxelVolume() * density; + double dose = energyDep / mass / CLHEP::gray; - deposited.dose += dose; - _eventDoseImage.AddValue(index, dose); - } + _eventDoseImage.AddValue(index, dose); // Get information from step // Particle @@ -494,22 +272,10 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* double alpha = (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b) * energyDep; double sqrtBeta = (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b) * energyDep; - if(index == 80) { // TODO remove - std::ofstream of{"/tmp/alpha", std::ios_base::app}; - of << "alpha: " << (interpol.alpha.a * kineticEnergyPerNucleon + interpol.alpha.b) << '\n'; - of << "sqrtBeta: " << (interpol.sqrtBeta.a * kineticEnergyPerNucleon + interpol.sqrtBeta.b) << '\n'; - of << "energyDep: " << energyDep << '\n'; - of << "--------------------\n"; - } - if(alpha < 0) alpha = 0; if(sqrtBeta < 0) sqrtBeta = 0; // Accumulate alpha/beta - deposited.alpha += alpha; - deposited.sqrtBeta += sqrtBeta; - - _eventEdepImage.AddValue(index, energyDep); _eventAlphaImage.AddValue(index, alpha); _eventSqrtBetaImage.AddValue(index, sqrtBeta); @@ -518,28 +284,126 @@ void GateBioDoseActor::UserSteppingActionInVoxel(const int index, const G4Step* } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -void GateBioDoseActor::ResetData() { - _eventCountImage.Reset(); +void GateBioDoseActor::updateData() { + auto const sqAlphaRef = _alphaRef * _alphaRef; + double const n = _currentEvent; - if(_enableEdep) _edepImage.Reset(); - if(_enableDose) _doseImage.Reset(); - _bioDoseImage.Reset(); - _alphaMixImage.Reset(); - _sqrtBetaMixImage.Reset(); - if(_enableRBE) _rbeImage.Reset(); + for(auto const& index: _voxelIndices) { + auto const hitEventCount = _hitEventCountImage.GetValue(index); + + auto const alphaMixMean = _alphaMixImage.GetValue(index) / hitEventCount; + auto const sqrtBetaMixMean = _sqrtBetaMixImage.GetValue(index) / hitEventCount; + auto const dose = _doseImage.GetValue(index); + auto const scaledDose = _doseScaleFactor * dose; + auto const sqScaledDose = scaledDose * scaledDose; + auto const delta = sqAlphaRef + 4 * _betaRef * + (alphaMixMean * scaledDose + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDose); + + double sqrtDelta = 0; + if(delta >= 0) + sqrtDelta = std::sqrt(delta); + + // Calculate biological dose and RBE + double biodose = 0; + double rbe = 0; + + if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0) + biodose = (-_alphaRef + sqrtDelta) / (2 * _betaRef); + if(biodose < 0) biodose = 0; // TODO improve + + if(scaledDose > 0) + rbe = biodose / scaledDose; + + if(_enableUncertainties) { + if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0 && sqrtDelta > 0 && _currentEvent > 0) { + double const doseMean = dose / n; + + double sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); + double sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); + double sumSquaredDose = _squaredDoseImage.GetValue(index); + + double pdBiodoseAlphaMix = scaledDose / sqrtDelta; + double pdBiodoseSqrtBetaMix = 2 * sqScaledDose * sqrtBetaMixMean / sqrtDelta; + double pdBiodoseDose = ( + (alphaMixMean * _doseScaleFactor) + + 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose + ) / sqrtDelta; + + double varAlphaMix = sumSquaredAlphaMix / hitEventCount - alphaMixMean * alphaMixMean; + double varSqrtBetaMix = sumSquaredSqrtBetaMix / hitEventCount - sqrtBetaMixMean * sqrtBetaMixMean; + double varDose = sumSquaredDose / n - doseMean * doseMean; + + double sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); + double sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); + double sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); + double covAlphaMixSqrtBetaMix = sumAlphaMixSqrtBetaMix / hitEventCount - alphaMixMean * sqrtBetaMixMean; + double covAlphaMixDose = sumAlphaMixDose / n - alphaMixMean * doseMean; + double covSqrtBetaMixDose = sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean; + + double partAlphaMix = pdBiodoseAlphaMix * pdBiodoseAlphaMix * varAlphaMix; + double partSqrtBetaMix = pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix * varSqrtBetaMix; + double partDose = pdBiodoseDose * pdBiodoseDose * varDose; + double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMix * pdBiodoseSqrtBetaMix * covAlphaMixSqrtBetaMix; + double partAlphaMixDose = 2 * pdBiodoseAlphaMix * pdBiodoseDose * covAlphaMixDose; + double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMix * pdBiodoseDose * covSqrtBetaMixDose; + + double uncertaintyBiodose = std::sqrt( + partAlphaMix + partSqrtBetaMix + partDose + + partAlphaMixSqrtBetaMix + partAlphaMixDose + partSqrtBetaMixDose + ) / biodose; + + _biodoseUncertaintyImage.SetValue(index, uncertaintyBiodose); + } else { + _biodoseUncertaintyImage.SetValue(index, 1); + } + } + + // Write data + if(_enableDose) _scaledDoseImage.SetValue(index, scaledDose); + _bioDoseImage.SetValue(index, biodose); + if(_enableRBE) _rbeImage.SetValue(index, rbe); + } +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::SaveData() { + GateDebugMessageInc("Actor", 4, "GateBioDoseActor::SaveData() known ion events / total events: " << _eventWithKnownIonCount << " / " << _eventCount << "\n"); + + updateData(); + + GateVActor::SaveData(); + + if(_enableEdep) _edepImage.SaveData(_currentEvent); + if(_enableDose) _scaledDoseImage.SaveData(_currentEvent); + if(_enableAlphaMix) _alphaMixImage.SaveData(_currentEvent); + if(_enableSqrtBetaMix) _sqrtBetaMixImage.SaveData(_currentEvent); + _bioDoseImage.SaveData(_currentEvent); + if(_enableRBE) _rbeImage.SaveData(_currentEvent); + if(_enableUncertainties) _biodoseUncertaintyImage.SaveData(_currentEvent); +} +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +void GateBioDoseActor::ResetData() { + _hitEventCountImage.Reset(); _eventEdepImage.Reset(); _eventDoseImage.Reset(); _eventAlphaImage.Reset(); _eventSqrtBetaImage.Reset(); + if(_enableEdep) _edepImage.Reset(); + _doseImage.Reset(); + if(_enableDose) _scaledDoseImage.Reset(); + _alphaMixImage.Reset(); + _sqrtBetaMixImage.Reset(); + _bioDoseImage.Reset(); + if(_enableRBE) _rbeImage.Reset(); + if(_enableUncertainties) { _biodoseUncertaintyImage.Reset(); - _squaredDoseImage.Reset(); _squaredAlphaMixImage.Reset(); _squaredSqrtBetaMixImage.Reset(); - _alphaMixSqrtBetaMixImage.Reset(); _alphaMixDoseImage.Reset(); _sqrtBetaMixDoseImage.Reset(); diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc index eb1b335ee..b6b979e35 100644 --- a/source/digits_hits/src/GateBioDoseActorMessenger.cc +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -56,10 +56,6 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { pEnableDoseCmd = std::make_unique(n, this); pEnableDoseCmd->SetGuidance("Enable dose output"); - n = base + "/enableBioDose"; - pEnableBioDoseCmd = std::make_unique(n, this); - pEnableBioDoseCmd->SetGuidance("Enable biodose output (default: true)"); - n = base + "/enableAlphaMix"; pEnableAlphaMixCmd = std::make_unique(n, this); pEnableAlphaMixCmd->SetGuidance("Enable alpha mix output"); @@ -87,7 +83,6 @@ void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(value)); else if(cmd == pEnableEdepCmd.get()) pBioDoseActor->SetEnableEdep(pEnableEdepCmd->GetNewBoolValue(value)); else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); - else if(cmd == pEnableBioDoseCmd.get()) pBioDoseActor->SetEnableBioDose(pEnableBioDoseCmd->GetNewBoolValue(value)); else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); else if(cmd == pEnableSqrtBetaMixCmd.get()) pBioDoseActor->SetEnableSqrtBetaMix(pEnableSqrtBetaMixCmd->GetNewBoolValue(value)); else if(cmd == pEnableRBECmd.get()) pBioDoseActor->SetEnableRBE(pEnableRBECmd->GetNewBoolValue(value)); From 80496b4d0e7d14c6ad6edd65079f15aa594f7b70 Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Mon, 19 Feb 2024 16:04:18 +0100 Subject: [PATCH 09/12] Use dose mean for uncertainty, adapt calculations to it --- .../digits_hits/include/GateBioDoseActor.hh | 1 + source/digits_hits/src/GateBioDoseActor.cc | 127 +++++++++++++++--- 2 files changed, 111 insertions(+), 17 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index 3d3365021..4cba1d89c 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -126,6 +126,7 @@ private: GateImageWithStatistic _bioDoseImage; GateImageWithStatistic _rbeImage; + GateImageWithStatistic _doseUncertaintyImage; GateImageWithStatistic _biodoseUncertaintyImage; GateImageWithStatistic _squaredDoseImage; GateImageWithStatistic _squaredAlphaMixImage; diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index f496fbff0..269c210b8 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -10,6 +10,7 @@ See LICENSE.md for further details #include "GateBioDoseActor.hh" #include "GateImageWithStatistic.hh" #include +#include #define GATE_BUFFERSIZE @@ -82,6 +83,7 @@ void GateBioDoseActor::Construct() { if(_enableRBE) setupImage(_rbeImage, "rbe"); if(_enableUncertainties) { + setupImage(_doseUncertaintyImage, "dose_uncertainty"); setupImage(_biodoseUncertaintyImage, "biodose_uncertainty"); setupImage(_squaredDoseImage); setupImage(_squaredAlphaMixImage); @@ -288,6 +290,52 @@ void GateBioDoseActor::updateData() { auto const sqAlphaRef = _alphaRef * _alphaRef; double const n = _currentEvent; + /* ******************* */ + // TODO remove + G4String basename = removeExtension(mSaveFilename); + G4String ext = getExtension(mSaveFilename); + auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { + SetOriginTransformAndFlagToImage(image); + image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); + image.Allocate(); + + if(!suffix.empty()) { + G4String filename = basename + "_" + suffix + "." + ext; + image.SetFilename(filename); + } + }; + + GateImageWithStatistic partAlphaMixImg; + GateImageWithStatistic partSqrtBetaMixImg; + GateImageWithStatistic partDoseImg; + GateImageWithStatistic partAlphaMixSqrtBetaMixImg; + GateImageWithStatistic partAlphaMixDoseImg; + GateImageWithStatistic partSqrtBetaMixDoseImg; + GateImageWithStatistic partSquaredBioDoseImg; + + GateImageWithStatistic squaredPdBioDoseAlphaMixImg; + GateImageWithStatistic squaredPdBioDoseSqrtBetaMixImg; + GateImageWithStatistic squaredPdBioDoseDoseImg; + GateImageWithStatistic varAlphaMixImg; + GateImageWithStatistic varSqrtBetaMixImg; + GateImageWithStatistic varDoseImg; + + setupImage(partAlphaMixImg, "part_alphamix"); + setupImage(partSqrtBetaMixImg, "part_sqrtbetamix"); + setupImage(partDoseImg, "part_dose"); + setupImage(partAlphaMixSqrtBetaMixImg, "part_alphamixsqrtbetamix"); + setupImage(partAlphaMixDoseImg, "part_alphamixdose"); + setupImage(partSqrtBetaMixDoseImg, "part_sqrtbetamixdose"); + setupImage(partSquaredBioDoseImg, "part_squaredbiodose"); + + setupImage(squaredPdBioDoseAlphaMixImg, "squared_pdbiodosealphamix"); + setupImage(squaredPdBioDoseSqrtBetaMixImg, "squared_pdbiodosesqrtbetamix"); + setupImage(squaredPdBioDoseDoseImg, "squared_pdbiodosedose"); + setupImage(varAlphaMixImg, "var_alphamix"); + setupImage(varSqrtBetaMixImg, "var_sqrtbetamix"); + setupImage(varDoseImg, "var_dose"); + /* ******************* */ + for(auto const& index: _voxelIndices) { auto const hitEventCount = _hitEventCountImage.GetValue(index); @@ -317,43 +365,69 @@ void GateBioDoseActor::updateData() { if(_enableUncertainties) { if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0 && sqrtDelta > 0 && _currentEvent > 0) { double const doseMean = dose / n; + double const biodoseMean = biodose / n; + auto const deltaMean = sqAlphaRef + 4 * _betaRef * + (alphaMixMean * scaledDose/n + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDose/n/n); + double sqrtDeltaMean = 0; + if(deltaMean >= 0) + sqrtDeltaMean = std::sqrt(deltaMean); double sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); double sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); double sumSquaredDose = _squaredDoseImage.GetValue(index); - double pdBiodoseAlphaMix = scaledDose / sqrtDelta; - double pdBiodoseSqrtBetaMix = 2 * sqScaledDose * sqrtBetaMixMean / sqrtDelta; + double pdBiodoseAlphaMixMean = scaledDose/n / sqrtDeltaMean; + double pdBiodoseSqrtBetaMixMean = 2 * sqScaledDose/n/n * sqrtBetaMixMean / sqrtDeltaMean; double pdBiodoseDose = ( (alphaMixMean * _doseScaleFactor) + - 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose - ) / sqrtDelta; + 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose/n + ) / sqrtDeltaMean; - double varAlphaMix = sumSquaredAlphaMix / hitEventCount - alphaMixMean * alphaMixMean; - double varSqrtBetaMix = sumSquaredSqrtBetaMix / hitEventCount - sqrtBetaMixMean * sqrtBetaMixMean; - double varDose = sumSquaredDose / n - doseMean * doseMean; + double varAlphaMixMean = (sumSquaredAlphaMix / hitEventCount - alphaMixMean * alphaMixMean) / hitEventCount; + double varSqrtBetaMixMean = (sumSquaredSqrtBetaMix / hitEventCount - sqrtBetaMixMean * sqrtBetaMixMean) / hitEventCount; + double varDose = (sumSquaredDose / n - doseMean * doseMean) / n; double sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); double sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); double sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); - double covAlphaMixSqrtBetaMix = sumAlphaMixSqrtBetaMix / hitEventCount - alphaMixMean * sqrtBetaMixMean; - double covAlphaMixDose = sumAlphaMixDose / n - alphaMixMean * doseMean; - double covSqrtBetaMixDose = sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean; + double covAlphaMixMeanSqrtBetaMixMean = (sumAlphaMixSqrtBetaMix / hitEventCount - alphaMixMean * sqrtBetaMixMean) / hitEventCount; + double covAlphaMixMeanDose = (sumAlphaMixDose / n - alphaMixMean * doseMean) / n; + double covSqrtBetaMixMeanDose = (sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean) / n; - double partAlphaMix = pdBiodoseAlphaMix * pdBiodoseAlphaMix * varAlphaMix; - double partSqrtBetaMix = pdBiodoseSqrtBetaMix * pdBiodoseSqrtBetaMix * varSqrtBetaMix; + double partAlphaMix = pdBiodoseAlphaMixMean * pdBiodoseAlphaMixMean * varAlphaMixMean; + double partSqrtBetaMix = pdBiodoseSqrtBetaMixMean * pdBiodoseSqrtBetaMixMean * varSqrtBetaMixMean; double partDose = pdBiodoseDose * pdBiodoseDose * varDose; - double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMix * pdBiodoseSqrtBetaMix * covAlphaMixSqrtBetaMix; - double partAlphaMixDose = 2 * pdBiodoseAlphaMix * pdBiodoseDose * covAlphaMixDose; - double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMix * pdBiodoseDose * covSqrtBetaMixDose; + double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMixMean * pdBiodoseSqrtBetaMixMean * covAlphaMixMeanSqrtBetaMixMean; + double partAlphaMixDose = 2 * pdBiodoseAlphaMixMean * pdBiodoseDose * covAlphaMixMeanDose; + double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMixMean * pdBiodoseDose * covSqrtBetaMixMeanDose; + double uncertaintyDose = std::sqrt(varDose) / doseMean; double uncertaintyBiodose = std::sqrt( partAlphaMix + partSqrtBetaMix + partDose + partAlphaMixSqrtBetaMix + partAlphaMixDose + partSqrtBetaMixDose - ) / biodose; - + ) / (biodose/n); + + { // TODO remove + partAlphaMixImg.SetValue(index, partAlphaMix); + partSqrtBetaMixImg.SetValue(index, partSqrtBetaMix); + partDoseImg.SetValue(index, partDose); + partAlphaMixDoseImg.SetValue(index, partAlphaMixDose); + partAlphaMixSqrtBetaMixImg.SetValue(index, partAlphaMixSqrtBetaMix); + partSqrtBetaMixDoseImg.SetValue(index, partSqrtBetaMixDose); + partSquaredBioDoseImg.SetValue(index, biodose * biodose); + + squaredPdBioDoseAlphaMixImg.SetValue(index, pdBiodoseAlphaMixMean * pdBiodoseAlphaMixMean); + squaredPdBioDoseSqrtBetaMixImg.SetValue(index, pdBiodoseSqrtBetaMixMean * pdBiodoseSqrtBetaMixMean); + squaredPdBioDoseDoseImg.SetValue(index, pdBiodoseDose * pdBiodoseDose); + varAlphaMixImg.SetValue(index, varAlphaMixMean); + varSqrtBetaMixImg.SetValue(index, varSqrtBetaMixMean); + varDoseImg.SetValue(index, varDose); + } + + _doseUncertaintyImage.SetValue(index, uncertaintyDose); _biodoseUncertaintyImage.SetValue(index, uncertaintyBiodose); } else { + _doseUncertaintyImage.SetValue(index, 1); _biodoseUncertaintyImage.SetValue(index, 1); } } @@ -363,6 +437,23 @@ void GateBioDoseActor::updateData() { _bioDoseImage.SetValue(index, biodose); if(_enableRBE) _rbeImage.SetValue(index, rbe); } + + { // TODO remove + partAlphaMixImg.SaveData(_currentEvent); + partSqrtBetaMixImg.SaveData(_currentEvent); + partDoseImg.SaveData(_currentEvent); + partAlphaMixDoseImg.SaveData(_currentEvent); + partAlphaMixSqrtBetaMixImg.SaveData(_currentEvent); + partSqrtBetaMixDoseImg.SaveData(_currentEvent); + partSquaredBioDoseImg.SaveData(_currentEvent); + + squaredPdBioDoseAlphaMixImg.SaveData(_currentEvent); + squaredPdBioDoseSqrtBetaMixImg.SaveData(_currentEvent); + squaredPdBioDoseDoseImg.SaveData(_currentEvent); + varAlphaMixImg.SaveData(_currentEvent); + varSqrtBetaMixImg.SaveData(_currentEvent); + varDoseImg.SaveData(_currentEvent); + } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -379,6 +470,7 @@ void GateBioDoseActor::SaveData() { if(_enableSqrtBetaMix) _sqrtBetaMixImage.SaveData(_currentEvent); _bioDoseImage.SaveData(_currentEvent); if(_enableRBE) _rbeImage.SaveData(_currentEvent); + if(_enableUncertainties) _doseUncertaintyImage.SaveData(_currentEvent); if(_enableUncertainties) _biodoseUncertaintyImage.SaveData(_currentEvent); } //----------------------------------------------------------------------------- @@ -400,6 +492,7 @@ void GateBioDoseActor::ResetData() { if(_enableRBE) _rbeImage.Reset(); if(_enableUncertainties) { + _doseUncertaintyImage.Reset(); _biodoseUncertaintyImage.Reset(); _squaredDoseImage.Reset(); _squaredAlphaMixImage.Reset(); From df47b5d21cc39b48fabe81b16eeae28a00fa3ece Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Wed, 21 Feb 2024 18:14:04 +0100 Subject: [PATCH 10/12] Add outputs for uncertainty details and hit-events --- .../digits_hits/include/GateBioDoseActor.hh | 35 ++- .../include/GateBioDoseActorMessenger.hh | 2 + source/digits_hits/src/GateBioDoseActor.cc | 224 +++++++----------- .../src/GateBioDoseActorMessenger.cc | 35 ++- 4 files changed, 142 insertions(+), 154 deletions(-) diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index 4cba1d89c..9bc150ddb 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -79,7 +79,9 @@ public: void SetEnableAlphaMix(bool e) { _enableAlphaMix = e; } void SetEnableSqrtBetaMix(bool e) { _enableSqrtBetaMix = e; } void SetEnableRBE(bool e) { _enableRBE = e; } - void SetEnableUncertainties(bool e) { _enableUncertainties = e; } + void SetEnableUncertainty(bool e) { _enableUncertainty = e; } + void SetEnableUncertaintyDetails(bool e) { _enableUncertaintyDetails = e; } + void SetEnableHitEventCount(bool e) { _enableHitEventCount = e; } protected: GateBioDoseActor(G4String name, G4int depth = 0); @@ -90,7 +92,7 @@ protected: private: //Counters - int _currentEvent; + int _currentEvent = 0; GateBioDoseActorMessenger _messenger; @@ -100,10 +102,11 @@ private: G4String _dataBase; G4String _cellLine; G4String _bioPhysicalModel; - double _alphaRef, _betaRef; + double _alphaRef = -1; + double _betaRef = -1; double _doseScaleFactor = 1.; - G4double _sobpWeight; + G4double _sobpWeight = 0; AlphaBetaInterpolTable _alphaBetaInterpolTable; @@ -135,13 +138,25 @@ private: GateImageWithStatistic _alphaMixDoseImage; GateImageWithStatistic _sqrtBetaMixDoseImage; + GateImageWithStatistic _pdBiodoseAlphaMixMean; + GateImageWithStatistic _pdBiodoseSqrtBetaMixMean; + GateImageWithStatistic _pdBiodoseDoseMean; + GateImageWithStatistic _varAlphaMixMeanImage; + GateImageWithStatistic _varSqrtBetaMixMeanImage; + GateImageWithStatistic _varDoseMeanImage; + GateImageWithStatistic _covAlphaMixMeanSqrtBetaMixMeanImage; + GateImageWithStatistic _covAlphaMixMeanDoseMean; + GateImageWithStatistic _covSqrtBetaMixMeanDoseMean; + // Outputs - bool _enableEdep; - bool _enableDose; - bool _enableAlphaMix; - bool _enableSqrtBetaMix; - bool _enableRBE; - bool _enableUncertainties; + bool _enableEdep = false; + bool _enableDose = true; + bool _enableAlphaMix = false; + bool _enableSqrtBetaMix = false; + bool _enableRBE = false; + bool _enableUncertainty = true; + bool _enableUncertaintyDetails = false; + bool _enableHitEventCount = false; // Extra information int _stepCount = 0; diff --git a/source/digits_hits/include/GateBioDoseActorMessenger.hh b/source/digits_hits/include/GateBioDoseActorMessenger.hh index ec388fcfc..07f959a0b 100644 --- a/source/digits_hits/include/GateBioDoseActorMessenger.hh +++ b/source/digits_hits/include/GateBioDoseActorMessenger.hh @@ -45,6 +45,8 @@ private: std::unique_ptr pEnableSqrtBetaMixCmd; std::unique_ptr pEnableRBECmd; std::unique_ptr pEnableUncertaintyCmd; + std::unique_ptr pEnableUncertaintyDetailsCmd; + std::unique_ptr pEnableHitEventCountCmd; }; #endif diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index 269c210b8..b6191ff06 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -17,17 +17,7 @@ See LICENSE.md for further details //----------------------------------------------------------------------------- GateBioDoseActor::GateBioDoseActor(G4String name, G4int depth): GateVImageActor(std::move(name), depth), - _currentEvent(0), - _messenger(this), - _alphaRef(-1), - _betaRef(-1), - _sobpWeight(0), - _enableEdep(false), - _enableDose(true), - _enableAlphaMix(false), - _enableSqrtBetaMix(false), - _enableRBE(false), - _enableUncertainties(true) + _messenger(this) { GateDebugMessageInc("Actor", 4, "GateBioDoseActor() -- begin\n"); } @@ -37,10 +27,7 @@ void GateBioDoseActor::Construct() { GateDebugMessageInc("Actor", 4, "GateBioDoseActor -- Construct - begin\n"); GateVImageActor::Construct(); - // Find G4_WATER G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"); - // Find OtherMaterial - //G4NistManager::Instance()->FindOrBuildMaterial(mOtherMaterial); G4EmParameters::Instance()->SetBuildCSDARange(true); // Enable callbacks BioDose @@ -67,7 +54,7 @@ void GateBioDoseActor::Construct() { } }; - setupImage(_hitEventCountImage); + setupImage(_hitEventCountImage, "hitevent_count"); setupImage(_eventEdepImage); setupImage(_eventDoseImage); @@ -82,7 +69,7 @@ void GateBioDoseActor::Construct() { setupImage(_bioDoseImage, "biodose"); if(_enableRBE) setupImage(_rbeImage, "rbe"); - if(_enableUncertainties) { + if(_enableUncertainty) { setupImage(_doseUncertaintyImage, "dose_uncertainty"); setupImage(_biodoseUncertaintyImage, "biodose_uncertainty"); setupImage(_squaredDoseImage); @@ -91,6 +78,18 @@ void GateBioDoseActor::Construct() { setupImage(_alphaMixSqrtBetaMixImage); setupImage(_alphaMixDoseImage); setupImage(_sqrtBetaMixDoseImage); + + if(_enableUncertaintyDetails) { + setupImage(_pdBiodoseAlphaMixMean, "pd_biodose_alphamixmean"); + setupImage(_pdBiodoseSqrtBetaMixMean, "pd_biodose_sqrtbetamixmean"); + setupImage(_pdBiodoseDoseMean, "pd_biodose_dosemean"); + setupImage(_varAlphaMixMeanImage, "var_alphamixmean"); + setupImage(_varSqrtBetaMixMeanImage, "var_sqrtbetamixmean"); + setupImage(_varDoseMeanImage, "var_dosemean"); + setupImage(_covAlphaMixMeanSqrtBetaMixMeanImage, "cov_alphamixmean_sqrtbetamixmean"); + setupImage(_covAlphaMixMeanDoseMean, "cov_alphamixmean_dosemean"); + setupImage(_covSqrtBetaMixMeanDoseMean, "cov_sqrtbetamixmean_dosemean"); + } } } @@ -215,7 +214,7 @@ void GateBioDoseActor::EndOfEventAction(const G4Event* e) { _alphaMixImage.AddValue(index, eventAlphaMix); _sqrtBetaMixImage.AddValue(index, eventSqrtBetaMix); - if(_enableUncertainties) { + if(_enableUncertainty) { _squaredDoseImage.AddValue(index, eventDose * eventDose); _squaredAlphaMixImage.AddValue(index, eventAlphaMix * eventAlphaMix); _squaredSqrtBetaMixImage.AddValue(index, eventSqrtBetaMix * eventSqrtBetaMix); @@ -290,52 +289,6 @@ void GateBioDoseActor::updateData() { auto const sqAlphaRef = _alphaRef * _alphaRef; double const n = _currentEvent; - /* ******************* */ - // TODO remove - G4String basename = removeExtension(mSaveFilename); - G4String ext = getExtension(mSaveFilename); - auto setupImage = [&](GateImageWithStatistic& image, std::string const& suffix = "") { - SetOriginTransformAndFlagToImage(image); - image.SetResolutionAndHalfSize(mResolution, mHalfSize, mPosition); - image.Allocate(); - - if(!suffix.empty()) { - G4String filename = basename + "_" + suffix + "." + ext; - image.SetFilename(filename); - } - }; - - GateImageWithStatistic partAlphaMixImg; - GateImageWithStatistic partSqrtBetaMixImg; - GateImageWithStatistic partDoseImg; - GateImageWithStatistic partAlphaMixSqrtBetaMixImg; - GateImageWithStatistic partAlphaMixDoseImg; - GateImageWithStatistic partSqrtBetaMixDoseImg; - GateImageWithStatistic partSquaredBioDoseImg; - - GateImageWithStatistic squaredPdBioDoseAlphaMixImg; - GateImageWithStatistic squaredPdBioDoseSqrtBetaMixImg; - GateImageWithStatistic squaredPdBioDoseDoseImg; - GateImageWithStatistic varAlphaMixImg; - GateImageWithStatistic varSqrtBetaMixImg; - GateImageWithStatistic varDoseImg; - - setupImage(partAlphaMixImg, "part_alphamix"); - setupImage(partSqrtBetaMixImg, "part_sqrtbetamix"); - setupImage(partDoseImg, "part_dose"); - setupImage(partAlphaMixSqrtBetaMixImg, "part_alphamixsqrtbetamix"); - setupImage(partAlphaMixDoseImg, "part_alphamixdose"); - setupImage(partSqrtBetaMixDoseImg, "part_sqrtbetamixdose"); - setupImage(partSquaredBioDoseImg, "part_squaredbiodose"); - - setupImage(squaredPdBioDoseAlphaMixImg, "squared_pdbiodosealphamix"); - setupImage(squaredPdBioDoseSqrtBetaMixImg, "squared_pdbiodosesqrtbetamix"); - setupImage(squaredPdBioDoseDoseImg, "squared_pdbiodosedose"); - setupImage(varAlphaMixImg, "var_alphamix"); - setupImage(varSqrtBetaMixImg, "var_sqrtbetamix"); - setupImage(varDoseImg, "var_dose"); - /* ******************* */ - for(auto const& index: _voxelIndices) { auto const hitEventCount = _hitEventCountImage.GetValue(index); @@ -362,66 +315,63 @@ void GateBioDoseActor::updateData() { if(scaledDose > 0) rbe = biodose / scaledDose; - if(_enableUncertainties) { + if(_enableUncertainty) { if(scaledDose > 0 && alphaMixMean != 0 && sqrtBetaMixMean != 0 && sqrtDelta > 0 && _currentEvent > 0) { - double const doseMean = dose / n; - double const biodoseMean = biodose / n; + auto const doseMean = dose / n; + auto const scaledDoseMean = scaledDose / n; + auto const sqScaledDoseMean = sqScaledDose / n / n; + auto const biodoseMean = biodose / n; auto const deltaMean = sqAlphaRef + 4 * _betaRef * - (alphaMixMean * scaledDose/n + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDose/n/n); - double sqrtDeltaMean = 0; + (alphaMixMean * scaledDoseMean + sqrtBetaMixMean * sqrtBetaMixMean * sqScaledDoseMean); + double sqrtDeltaMean = 0.; if(deltaMean >= 0) sqrtDeltaMean = std::sqrt(deltaMean); - double sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); - double sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); - double sumSquaredDose = _squaredDoseImage.GetValue(index); + auto sumSquaredAlphaMix = _squaredAlphaMixImage.GetValue(index); + auto sumSquaredSqrtBetaMix = _squaredSqrtBetaMixImage.GetValue(index); + auto sumSquaredDose = _squaredDoseImage.GetValue(index); - double pdBiodoseAlphaMixMean = scaledDose/n / sqrtDeltaMean; - double pdBiodoseSqrtBetaMixMean = 2 * sqScaledDose/n/n * sqrtBetaMixMean / sqrtDeltaMean; - double pdBiodoseDose = ( + auto pdBiodoseAlphaMixMean = scaledDoseMean / sqrtDeltaMean; + auto pdBiodoseSqrtBetaMixMean = 2 * sqScaledDoseMean * sqrtBetaMixMean / sqrtDeltaMean; + auto pdBiodoseDoseMean = ( (alphaMixMean * _doseScaleFactor) + - 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDose/n + 2 * sqrtBetaMixMean * sqrtBetaMixMean * _doseScaleFactor * scaledDoseMean ) / sqrtDeltaMean; - double varAlphaMixMean = (sumSquaredAlphaMix / hitEventCount - alphaMixMean * alphaMixMean) / hitEventCount; - double varSqrtBetaMixMean = (sumSquaredSqrtBetaMix / hitEventCount - sqrtBetaMixMean * sqrtBetaMixMean) / hitEventCount; - double varDose = (sumSquaredDose / n - doseMean * doseMean) / n; - - double sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); - double sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); - double sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); - double covAlphaMixMeanSqrtBetaMixMean = (sumAlphaMixSqrtBetaMix / hitEventCount - alphaMixMean * sqrtBetaMixMean) / hitEventCount; - double covAlphaMixMeanDose = (sumAlphaMixDose / n - alphaMixMean * doseMean) / n; - double covSqrtBetaMixMeanDose = (sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean) / n; - - double partAlphaMix = pdBiodoseAlphaMixMean * pdBiodoseAlphaMixMean * varAlphaMixMean; - double partSqrtBetaMix = pdBiodoseSqrtBetaMixMean * pdBiodoseSqrtBetaMixMean * varSqrtBetaMixMean; - double partDose = pdBiodoseDose * pdBiodoseDose * varDose; - double partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMixMean * pdBiodoseSqrtBetaMixMean * covAlphaMixMeanSqrtBetaMixMean; - double partAlphaMixDose = 2 * pdBiodoseAlphaMixMean * pdBiodoseDose * covAlphaMixMeanDose; - double partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMixMean * pdBiodoseDose * covSqrtBetaMixMeanDose; - - double uncertaintyDose = std::sqrt(varDose) / doseMean; - double uncertaintyBiodose = std::sqrt( + auto varAlphaMixMean = (sumSquaredAlphaMix / hitEventCount - alphaMixMean * alphaMixMean) / hitEventCount; + auto varSqrtBetaMixMean = (sumSquaredSqrtBetaMix / hitEventCount - sqrtBetaMixMean * sqrtBetaMixMean) / hitEventCount; + auto varDoseMean = (sumSquaredDose / n - doseMean * doseMean) / n; + + auto sumAlphaMixSqrtBetaMix = _alphaMixSqrtBetaMixImage.GetValue(index); + auto sumAlphaMixDose = _alphaMixDoseImage.GetValue(index); + auto sumSqrtBetaMixDose = _sqrtBetaMixDoseImage.GetValue(index); + auto covAlphaMixMeanSqrtBetaMixMean = (sumAlphaMixSqrtBetaMix / hitEventCount - alphaMixMean * sqrtBetaMixMean) / hitEventCount; + auto covAlphaMixMeanDoseMean = (sumAlphaMixDose / n - alphaMixMean * doseMean) / n; + auto covSqrtBetaMixMeanDoseMean = (sumSqrtBetaMixDose / n - sqrtBetaMixMean * doseMean) / n; + + auto partAlphaMix = pdBiodoseAlphaMixMean * pdBiodoseAlphaMixMean * varAlphaMixMean; + auto partSqrtBetaMix = pdBiodoseSqrtBetaMixMean * pdBiodoseSqrtBetaMixMean * varSqrtBetaMixMean; + auto partDose = pdBiodoseDoseMean * pdBiodoseDoseMean * varDoseMean; + auto partAlphaMixSqrtBetaMix = 2 * pdBiodoseAlphaMixMean * pdBiodoseSqrtBetaMixMean * covAlphaMixMeanSqrtBetaMixMean; + auto partAlphaMixDose = 2 * pdBiodoseAlphaMixMean * pdBiodoseDoseMean * covAlphaMixMeanDoseMean; + auto partSqrtBetaMixDose = 2 * pdBiodoseSqrtBetaMixMean * pdBiodoseDoseMean * covSqrtBetaMixMeanDoseMean; + + auto uncertaintyDose = std::sqrt(varDoseMean) / doseMean; + auto uncertaintyBiodose = std::sqrt( partAlphaMix + partSqrtBetaMix + partDose + partAlphaMixSqrtBetaMix + partAlphaMixDose + partSqrtBetaMixDose - ) / (biodose/n); - - { // TODO remove - partAlphaMixImg.SetValue(index, partAlphaMix); - partSqrtBetaMixImg.SetValue(index, partSqrtBetaMix); - partDoseImg.SetValue(index, partDose); - partAlphaMixDoseImg.SetValue(index, partAlphaMixDose); - partAlphaMixSqrtBetaMixImg.SetValue(index, partAlphaMixSqrtBetaMix); - partSqrtBetaMixDoseImg.SetValue(index, partSqrtBetaMixDose); - partSquaredBioDoseImg.SetValue(index, biodose * biodose); - - squaredPdBioDoseAlphaMixImg.SetValue(index, pdBiodoseAlphaMixMean * pdBiodoseAlphaMixMean); - squaredPdBioDoseSqrtBetaMixImg.SetValue(index, pdBiodoseSqrtBetaMixMean * pdBiodoseSqrtBetaMixMean); - squaredPdBioDoseDoseImg.SetValue(index, pdBiodoseDose * pdBiodoseDose); - varAlphaMixImg.SetValue(index, varAlphaMixMean); - varSqrtBetaMixImg.SetValue(index, varSqrtBetaMixMean); - varDoseImg.SetValue(index, varDose); + ) / biodoseMean; + + if(_enableUncertaintyDetails) { + _pdBiodoseAlphaMixMean.SetValue(index, pdBiodoseAlphaMixMean); + _pdBiodoseSqrtBetaMixMean.SetValue(index, pdBiodoseSqrtBetaMixMean); + _pdBiodoseDoseMean.SetValue(index, pdBiodoseDoseMean); + _varAlphaMixMeanImage.SetValue(index, varAlphaMixMean); + _varSqrtBetaMixMeanImage.SetValue(index, varSqrtBetaMixMean); + _varDoseMeanImage.SetValue(index, varDoseMean); + _covAlphaMixMeanSqrtBetaMixMeanImage.SetValue(index, covAlphaMixMeanSqrtBetaMixMean); + _covAlphaMixMeanDoseMean.SetValue(index, covAlphaMixMeanDoseMean); + _covSqrtBetaMixMeanDoseMean.SetValue(index, covSqrtBetaMixMeanDoseMean); } _doseUncertaintyImage.SetValue(index, uncertaintyDose); @@ -437,23 +387,6 @@ void GateBioDoseActor::updateData() { _bioDoseImage.SetValue(index, biodose); if(_enableRBE) _rbeImage.SetValue(index, rbe); } - - { // TODO remove - partAlphaMixImg.SaveData(_currentEvent); - partSqrtBetaMixImg.SaveData(_currentEvent); - partDoseImg.SaveData(_currentEvent); - partAlphaMixDoseImg.SaveData(_currentEvent); - partAlphaMixSqrtBetaMixImg.SaveData(_currentEvent); - partSqrtBetaMixDoseImg.SaveData(_currentEvent); - partSquaredBioDoseImg.SaveData(_currentEvent); - - squaredPdBioDoseAlphaMixImg.SaveData(_currentEvent); - squaredPdBioDoseSqrtBetaMixImg.SaveData(_currentEvent); - squaredPdBioDoseDoseImg.SaveData(_currentEvent); - varAlphaMixImg.SaveData(_currentEvent); - varSqrtBetaMixImg.SaveData(_currentEvent); - varDoseImg.SaveData(_currentEvent); - } } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -470,8 +403,23 @@ void GateBioDoseActor::SaveData() { if(_enableSqrtBetaMix) _sqrtBetaMixImage.SaveData(_currentEvent); _bioDoseImage.SaveData(_currentEvent); if(_enableRBE) _rbeImage.SaveData(_currentEvent); - if(_enableUncertainties) _doseUncertaintyImage.SaveData(_currentEvent); - if(_enableUncertainties) _biodoseUncertaintyImage.SaveData(_currentEvent); + if(_enableUncertainty) { + _doseUncertaintyImage.SaveData(_currentEvent); + _biodoseUncertaintyImage.SaveData(_currentEvent); + + if(_enableUncertaintyDetails) { + _pdBiodoseAlphaMixMean.SaveData(_currentEvent); + _pdBiodoseSqrtBetaMixMean.SaveData(_currentEvent); + _pdBiodoseDoseMean.SaveData(_currentEvent); + _varAlphaMixMeanImage.SaveData(_currentEvent); + _varSqrtBetaMixMeanImage.SaveData(_currentEvent); + _varDoseMeanImage.SaveData(_currentEvent); + _covAlphaMixMeanSqrtBetaMixMeanImage.SaveData(_currentEvent); + _covAlphaMixMeanDoseMean.SaveData(_currentEvent); + _covSqrtBetaMixMeanDoseMean.SaveData(_currentEvent); + } + } + if(_enableHitEventCount) _hitEventCountImage.SaveData(_currentEvent); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -491,7 +439,7 @@ void GateBioDoseActor::ResetData() { _bioDoseImage.Reset(); if(_enableRBE) _rbeImage.Reset(); - if(_enableUncertainties) { + if(_enableUncertainty) { _doseUncertaintyImage.Reset(); _biodoseUncertaintyImage.Reset(); _squaredDoseImage.Reset(); @@ -500,6 +448,18 @@ void GateBioDoseActor::ResetData() { _alphaMixSqrtBetaMixImage.Reset(); _alphaMixDoseImage.Reset(); _sqrtBetaMixDoseImage.Reset(); + + if(_enableUncertaintyDetails) { + _pdBiodoseAlphaMixMean.Reset(); + _pdBiodoseSqrtBetaMixMean.Reset(); + _pdBiodoseDoseMean.Reset(); + _varAlphaMixMeanImage.Reset(); + _varSqrtBetaMixMeanImage.Reset(); + _varDoseMeanImage.Reset(); + _covAlphaMixMeanSqrtBetaMixMeanImage.Reset(); + _covAlphaMixMeanDoseMean.Reset(); + _covSqrtBetaMixMeanDoseMean.Reset(); + } } } //----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateBioDoseActorMessenger.cc b/source/digits_hits/src/GateBioDoseActorMessenger.cc index b6b979e35..517cb783b 100644 --- a/source/digits_hits/src/GateBioDoseActorMessenger.cc +++ b/source/digits_hits/src/GateBioDoseActorMessenger.cc @@ -71,22 +71,33 @@ void GateBioDoseActorMessenger::BuildCommands(G4String base) { n = base + "/enableUncertainty"; pEnableUncertaintyCmd = std::make_unique(n, this); pEnableUncertaintyCmd->SetGuidance("Enable uncertainty outputs"); + + n = base + "/enableUncertaintyDetails"; + pEnableUncertaintyDetailsCmd = std::make_unique(n, this); + pEnableUncertaintyDetailsCmd->SetGuidance("Enable uncertainty intermediary outputs"); + + n = base + "/enableHitEventCount"; + pEnableHitEventCountCmd = std::make_unique(n, this); + pEnableHitEventCountCmd->SetGuidance("Enable hit-event count output"); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- void GateBioDoseActorMessenger::SetNewValue(G4UIcommand* cmd, G4String value) { - if(cmd == pDoseScaleFactorCmd.get()) pBioDoseActor->SetDoseScaleFactor(pDoseScaleFactorCmd->GetNewDoubleValue(value)); - else if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(value)); - else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(value)); - else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(value); - else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(value); - else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(value)); - else if(cmd == pEnableEdepCmd.get()) pBioDoseActor->SetEnableEdep(pEnableEdepCmd->GetNewBoolValue(value)); - else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); - else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); - else if(cmd == pEnableSqrtBetaMixCmd.get()) pBioDoseActor->SetEnableSqrtBetaMix(pEnableSqrtBetaMixCmd->GetNewBoolValue(value)); - else if(cmd == pEnableRBECmd.get()) pBioDoseActor->SetEnableRBE(pEnableRBECmd->GetNewBoolValue(value)); - else if(cmd == pEnableUncertaintyCmd.get()) pBioDoseActor->SetEnableUncertainties(pEnableUncertaintyCmd->GetNewBoolValue(value)); + if(cmd == pDoseScaleFactorCmd.get()) pBioDoseActor->SetDoseScaleFactor(pDoseScaleFactorCmd->GetNewDoubleValue(value)); + else if(cmd == pAlphaRefCmd.get()) pBioDoseActor->SetAlphaRef(pAlphaRefCmd->GetNewDoubleValue(value)); + else if(cmd == pBetaRefCmd.get()) pBioDoseActor->SetBetaRef(pBetaRefCmd->GetNewDoubleValue(value)); + else if(cmd == pCellLineCmd.get()) pBioDoseActor->SetCellLine(value); + else if(cmd == pBioPhysicalModelCmd.get()) pBioDoseActor->SetBioPhysicalModel(value); + else if(cmd == pSOBPWeightCmd.get()) pBioDoseActor->SetSOBPWeight(pSOBPWeightCmd->GetNewDoubleValue(value)); + else if(cmd == pEnableEdepCmd.get()) pBioDoseActor->SetEnableEdep(pEnableEdepCmd->GetNewBoolValue(value)); + else if(cmd == pEnableDoseCmd.get()) pBioDoseActor->SetEnableDose(pEnableDoseCmd->GetNewBoolValue(value)); + else if(cmd == pEnableAlphaMixCmd.get()) pBioDoseActor->SetEnableAlphaMix(pEnableAlphaMixCmd->GetNewBoolValue(value)); + else if(cmd == pEnableSqrtBetaMixCmd.get()) pBioDoseActor->SetEnableSqrtBetaMix(pEnableSqrtBetaMixCmd->GetNewBoolValue(value)); + else if(cmd == pEnableRBECmd.get()) pBioDoseActor->SetEnableRBE(pEnableRBECmd->GetNewBoolValue(value)); + else if(cmd == pEnableUncertaintyCmd.get()) pBioDoseActor->SetEnableUncertainty(pEnableUncertaintyCmd->GetNewBoolValue(value)); + else if(cmd == pEnableUncertaintyDetailsCmd.get()) + pBioDoseActor->SetEnableUncertaintyDetails(pEnableUncertaintyDetailsCmd->GetNewBoolValue(value)); + else if(cmd == pEnableHitEventCountCmd.get()) pBioDoseActor->SetEnableHitEventCount(pEnableHitEventCountCmd->GetNewBoolValue(value)); GateImageActorMessenger::SetNewValue(cmd, value); } From 90ee6139185399d67746d6f92a84ee56e363ee1c Mon Sep 17 00:00:00 2001 From: Alexis Pereda Date: Mon, 18 Mar 2024 16:06:23 +0100 Subject: [PATCH 11/12] Add documentation and biological data --- AUTHORS | 1 + data/HSG_MMKM.db | 431 ++++++++++++++++++ data/HSG_NanOx.db | 416 +++++++++++++++++ ...to_interact_with_the_simulation_actors.rst | 84 ++++ .../digits_hits/include/GateBioDoseActor.hh | 14 +- source/digits_hits/src/GateBioDoseActor.cc | 40 +- 6 files changed, 959 insertions(+), 27 deletions(-) create mode 100644 data/HSG_MMKM.db create mode 100644 data/HSG_NanOx.db diff --git a/AUTHORS b/AUTHORS index e1dcc9c7e..102b936ab 100644 --- a/AUTHORS +++ b/AUTHORS @@ -48,6 +48,7 @@ Miklavec Mojca Mouton Claire Padilla Fatima Patay Gergely +Pereda Alexis Perez Hector Pommranz Christian Pike Lucy diff --git a/data/HSG_MMKM.db b/data/HSG_MMKM.db new file mode 100644 index 000000000..ce14d86c9 --- /dev/null +++ b/data/HSG_MMKM.db @@ -0,0 +1,431 @@ +Fragment 1; +0.1 1.57949019869911 0.047 +0.125 1.54247543384665 0.047 +0.15 1.4972852739187 0.047 +0.175 1.44966916924569 0.047 +0.2 1.38612514094072 0.047 +0.225 1.35761278567687 0.047 +0.25 1.31545384310958 0.047 +0.275 1.25179827306997 0.047 +0.3 1.20974390659241 0.047 +0.325 1.17270042701559 0.047 +0.35 1.13316814861987 0.047 +0.375 1.10341388783488 0.047 +0.4 1.083952291255 0.047 +0.425 1.03774318175064 0.047 +0.45 1.0281585942452 0.047 +0.475 0.988043948524032 0.047 +0.5 0.97039156529744 0.047 +0.525 0.941908178578341 0.047 +0.55 0.927571467254184 0.047 +0.6 0.893462336251923 0.047 +0.625 0.867341543482448 0.047 +0.65 0.846482254096548 0.047 +0.675 0.834568478438686 0.047 +0.7 0.821647475911647 0.047 +0.725 0.819826822713914 0.047 +0.75 0.805696154497498 0.047 +0.775 0.773649475617298 0.047 +0.8 0.762286113558253 0.047 +0.825 0.76343052337703 0.047 +0.85 0.749618054526763 0.047 +0.875 0.740088873111887 0.047 +0.9 0.722906452814284 0.047 +0.925 0.702859335660836 0.047 +0.95 0.694479602708279 0.047 +0.975 0.680325415638843 0.047 +0.988 0.70064262322935 0.047 +1 0.684390916768456 0.047 +1.375 0.666682549091409 0.047 +1.438 0.668564320475712 0.047 +1.5 0.622320144027226 0.047 +1.562 0.624851409142543 0.047 +1.625 0.603719538295307 0.047 +1.688 0.57716356055289 0.047 +1.75 0.574089085279603 0.047 +1.875 0.549991735425752 0.047 +2 0.526984423813353 0.047 +2.125 0.516236634533239 0.047 +2.25 0.487862571348046 0.047 +2.375 0.49214023309534 0.047 +2.5 0.459952553606391 0.047 +2.75 0.447349344901264 0.047 +3 0.444231765576353 0.047 +3.25 0.404942320848148 0.047 +3.5 0.392024007503326 0.047 +3.875 0.368160854149416 0.047 +4.25 0.360048787799576 0.047 +4.625 0.34873737425516 0.047 +5 0.331398700870807 0.047 +6 0.309435708584424 0.047 +7 0.293318209651091 0.047 +7.5 0.286659900758886 0.047 +8 0.282276169047493 0.047 +9 0.276902199013123 0.047 +10 0.271999732662595 0.047 +13 0.255211087252365 0.047 +14 0.241710674898587 0.047 +14.5 0.247787843853242 0.047 +15 0.243816674391015 0.047 +16 0.242464451793647 0.047 +17 0.243104936903852 0.047 +18.5 0.235094920292335 0.047 +20 0.225703564512416 0.047 +22.5 0.224795303783935 0.047 +25 0.231892670927331 0.047 +30 0.22226532371744 0.047 +35 0.219972180339683 0.047 +40 0.216583659599233 0.047 +42.5 0.20906029509513 0.047 +45 0.207829415228691 0.047 +50 0.202964450152741 0.047 +60 0.20627190303056 0.047 +70 0.209188910942346 0.047 +72.5 0.206332704243841 0.047 +75 0.208513246599234 0.047 +80 0.20166199491309 0.047 +85 0.204394018440942 0.047 +87.5 0.205084904351221 0.047 +90 0.196876372738078 0.047 +100 0.202595018458597 0.047 +110 0.203738028455536 0.047 +115 0.200152488987945 0.047 +120 0.201925583783688 0.047 +125 0.200833355702528 0.047 +130 0.200799746817376 0.047 +132.5 0.203244637136372 0.047 +135 0.201657368322226 0.047 +140 0.204013567966513 0.047 +160 0.199029925040146 0.047 +165 0.196711394763513 0.047 +170 0.198645019275942 0.047 +175 0.20001442557956 0.047 +180 0.196907511474943 0.047 +190 0.193554689166491 0.047 +200 0.197126711379263 0.047 +212.5 0.195925922147248 0.047 +225 0.195891925842979 0.047 +237.5 0.196353101850613 0.047 +250 0.199361141429766 0.047 +275 0.19614345428252 0.047 +300 0.196771041884298 0.047 +Fragment 2; +0.1 1.24815984258466 0.047 +0.115 1.22746833860529 0.047 +0.125 1.22568511984478 0.047 +0.135 1.22120019513435 0.047 +0.15 1.24004341046125 0.047 +0.175 1.27041403877105 0.047 +0.2 1.31284762557775 0.047 +0.225 1.34915114542436 0.047 +0.25 1.29916973452767 0.047 +0.275 1.3465188663924 0.047 +0.3 1.38853140901389 0.047 +0.325 1.42088580637717 0.047 +0.35 1.4651409822187 0.047 +0.375 1.49847344090863 0.047 +0.4 1.53231440665997 0.047 +0.425 1.55732901298041 0.047 +0.45 1.58887323257286 0.047 +0.475 1.61656993380224 0.047 +0.5 1.63477051521759 0.047 +0.525 1.64792086696669 0.047 +0.55 1.66367944667746 0.047 +0.56 1.67113349322127 0.047 +0.575 1.6782992316412 0.047 +0.58 1.68293757263274 0.047 +0.6 1.63902947990545 0.047 +0.625 1.65875805730164 0.047 +0.65 1.66760265894815 0.047 +0.675 1.67681552296571 0.047 +0.7 1.69362129846441 0.047 +0.725 1.69910794646881 0.047 +0.75 1.70121063944025 0.047 +0.775 1.70745567180774 0.047 +0.8 1.71339902777767 0.047 +0.825 1.71126047429641 0.047 +0.85 1.71475229369128 0.047 +0.875 1.71500818224546 0.047 +0.9 1.71253335416299 0.047 +0.925 1.71132512825384 0.047 +0.95 1.71042246319104 0.047 +0.975 1.70675238369539 0.047 +1 1.70283361040198 0.047 +1.125 1.67711028793932 0.047 +1.25 1.64115820323465 0.047 +1.312 1.61601704313548 0.047 +1.375 1.58515430782202 0.047 +1.438 1.57607563637513 0.047 +1.5 1.5591009949567 0.047 +1.562 1.5277574022372 0.047 +1.625 1.48072600340654 0.047 +1.688 1.47294230638862 0.047 +1.75 1.44849869172935 0.047 +1.8 1.42273654594659 0.047 +1.875 1.41554337214412 0.047 +1.9 1.36718645244497 0.047 +2 1.3484944522396 0.047 +2.125 1.29311494810033 0.047 +2.25 1.26544839044617 0.047 +2.375 1.22545364518765 0.047 +2.5 1.19956144020413 0.047 +2.75 1.10928181163848 0.047 +3 1.04033857934838 0.047 +3.25 0.995540726050855 0.047 +3.5 0.939271653992763 0.047 +3.75 0.898083473945884 0.047 +3.875 0.871295000837716 0.047 +4.25 0.816071445456982 0.047 +4.625 0.783037807932898 0.047 +4.8 0.767416691831474 0.047 +5 0.726291169041016 0.047 +5.5 0.68394019561239 0.047 +6 0.673790001122744 0.047 +7 0.601740849822646 0.047 +7.5 0.564825587003732 0.047 +8 0.542893492290651 0.047 +9 0.507043541377444 0.047 +9.5 0.499657624288272 0.047 +10 0.473240404403522 0.047 +12 0.439991881379019 0.047 +13 0.420059338583068 0.047 +14 0.413250267461632 0.047 +14.5 0.391318882034672 0.047 +15 0.39654305362603 0.047 +16 0.378167329233946 0.047 +17 0.371518217044929 0.047 +18.5 0.356832473293254 0.047 +20 0.337732188538681 0.047 +22.5 0.343088845534802 0.047 +25 0.309216389064121 0.047 +30 0.304885443714819 0.047 +35 0.282440098717072 0.047 +37.5 0.282404904154488 0.047 +40 0.27468125989886 0.047 +42.5 0.272361322449968 0.047 +45 0.26413866400718 0.047 +50 0.261870895854828 0.047 +60 0.252706707697897 0.047 +65 0.254545600880529 0.047 +70 0.242680258778023 0.047 +72.5 0.236189063172561 0.047 +75 0.24217093029468 0.047 +80 0.241967783234024 0.047 +85 0.237417901026164 0.047 +87.5 0.224328291081684 0.047 +90 0.234060411847778 0.047 +100 0.232027163414251 0.047 +110 0.227104795366579 0.047 +120 0.225454628016225 0.047 +130 0.223244894695277 0.047 +132.5 0.222853212506675 0.047 +135 0.217177779416561 0.047 +140 0.220034117529259 0.047 +160 0.21958561786659 0.047 +165 0.22058133341671 0.047 +170 0.219412021777746 0.047 +175 0.214836787141238 0.047 +180 0.216984920563624 0.047 +185 0.215153173813492 0.047 +190 0.210218241826287 0.047 +195 0.21383991264429 0.047 +200 0.216711620340268 0.047 +212.5 0.211443948581367 0.047 +225 0.215000201019383 0.047 +237.5 0.214196466693849 0.047 +250 0.206178971320345 0.047 +275 0.214644281664695 0.047 +300 0.206571356500144 0.047 +400 0.205153803569437 0.047 +500 0.205572062099992 0.047 +600 0.203623117877382 0.047 +700 0.200803428881801 0.047 +800 0.203894585196959 0.047 +900 0.199942243041688 0.047 +1000 0.201760692595498 0.047 +Fragment 6; +0.15 0.46771870915923 0.047 +0.175 0.45111001174474 0.047 +0.2 0.442579578522224 0.047 +0.225 0.43993529644597 0.047 +0.25 0.441768729083835 0.047 +0.275 0.429387809673272 0.047 +0.3 0.431473559281486 0.047 +0.35 0.439964885784134 0.047 +0.375 0.4453294251114 0.047 +0.4 0.451464155151481 0.047 +0.45 0.463408551092484 0.047 +0.475 0.46873158590075 0.047 +0.5 0.47613389481539 0.047 +0.525 0.456852622000955 0.047 +0.53 0.459498587392891 0.047 +0.55 0.460728170810213 0.047 +0.575 0.468338539236766 0.047 +0.6 0.475326384659887 0.047 +0.625 0.479916921469282 0.047 +0.65 0.486957860339879 0.047 +0.675 0.492065332768845 0.047 +0.7 0.498205945650057 0.047 +0.725 0.504736486846607 0.047 +0.75 0.502410693215626 0.047 +0.775 0.505867264051584 0.047 +0.8 0.514459720257531 0.047 +0.825 0.518473905794478 0.047 +0.85 0.52388804064196 0.047 +0.875 0.531222964603523 0.047 +0.9 0.53590714837395 0.047 +0.925 0.541569578249899 0.047 +0.95 0.548719272413843 0.047 +0.975 0.553286544219464 0.047 +0.988 0.556281712211522 0.047 +1 0.562795749287014 0.047 +1.125 0.591653025355066 0.047 +1.25 0.591137207023418 0.047 +1.312 0.60733679768358 0.047 +1.375 0.622003460017602 0.047 +1.438 0.635098466581411 0.047 +1.5 0.65299583066629 0.047 +1.562 0.66946775657227 0.047 +1.625 0.67977620005098 0.047 +1.688 0.693572811017881 0.047 +1.75 0.712159887126403 0.047 +1.875 0.727856610759531 0.047 +2 0.888962880223317 0.047 +2.125 0.759156873161599 0.047 +2.25 0.793361306942189 0.047 +2.5 0.849708323816713 0.047 +2.75 0.910260069524955 0.047 +3 0.976484322022195 0.047 +3.25 1.03816132450917 0.047 +3.5 1.09473276781447 0.047 +3.875 1.15288725796274 0.047 +4.25 1.22456704960088 0.047 +4.625 1.32422535930487 0.047 +5 1.37995585636382 0.047 +6 1.52956520927953 0.047 +7 1.63993938234722 0.047 +7.5 1.66828351827068 0.047 +8 1.69427749920672 0.047 +9 1.71119675143086 0.047 +10 1.70798646477594 0.047 +13 1.6106151202602 0.047 +14 1.58210303587246 0.047 +14.5 1.54805692116174 0.047 +15 1.53932394539696 0.047 +16 1.50354560433042 0.047 +17 1.43526086238606 0.047 +18.5 1.399256827658 0.047 +20 1.34314618467131 0.047 +22.5 1.2419292102451 0.047 +25 1.17541954239796 0.047 +27 1.12778867078749 0.047 +30 1.06823883613599 0.047 +32 1.02448408520427 0.047 +35 0.95657931321957 0.047 +37.5 0.927925494033008 0.047 +40 0.88692286564128 0.047 +42.5 0.842987024527601 0.047 +45 0.82114992786009 0.047 +50 0.778912417769797 0.047 +60 0.695630575871119 0.047 +65 0.652122921046971 0.047 +70 0.63331278944267 0.047 +72.5 0.611366721110614 0.047 +75 0.610805552932955 0.047 +80 0.592297511546729 0.047 +85 0.561119937378219 0.047 +87.5 0.560842376874587 0.047 +90 0.551197241208067 0.047 +100 0.523065078596558 0.047 +110 0.435238717853221 0.047 +120 0.474587966653229 0.047 +130 0.459646867165241 0.047 +132.5 0.453614530001866 0.047 +135 0.450636736824129 0.047 +140 0.445401074156315 0.047 +150 0.433357233053043 0.047 +160 0.424317118814647 0.047 +165 0.418435518267906 0.047 +175 0.408097739693348 0.047 +180 0.406822210554081 0.047 +185 0.399617574107289 0.047 +190 0.3960727715728 0.047 +195 0.391721887925199 0.047 +200 0.390341788246792 0.047 +212.5 0.3832824910486 0.047 +225 0.37363717529352 0.047 +237.5 0.367024071361367 0.047 +250 0.356661599411494 0.047 +275 0.355088605178929 0.047 +300 0.343722067778925 0.047 +Fragment 8; +0.1 0.532461736696031 0.047 +0.2 0.371901005763506 0.047 +0.25 0.363828363120222 0.047 +0.3 0.362712366880964 0.047 +0.35 0.366673649486981 0.047 +0.4 0.371538103944604 0.047 +0.45 0.377892108746988 0.047 +0.5 0.385642537689374 0.047 +0.55 0.371077107579352 0.047 +0.6 0.378556192415654 0.047 +0.7 0.393462400899198 0.047 +0.75 0.400098795861757 0.047 +0.8 0.407720935506561 0.047 +0.85 0.415379268610058 0.047 +0.9 0.422353470426836 0.047 +0.95 0.429477298199205 0.047 +1 0.437998976622219 0.047 +1.25 0.464336584114938 0.047 +1.5 0.498419168587147 0.047 +1.75 0.529578500134331 0.047 +2.25 0.578520279357553 0.047 +2.5 0.60638189324565 0.047 +2.75 0.638646625796802 0.047 +3 0.677427211519294 0.047 +3.5 0.751691858959154 0.047 +5 0.925608093458013 0.047 +6 1.05940180232884 0.047 +7 1.15948732484172 0.047 +8 1.26845445214636 0.047 +10 1.46833895356548 0.047 +13 1.64849230555153 0.047 +15 1.69646944214935 0.047 +20 1.69024842312232 0.047 +25 1.59903475226544 0.047 +30 1.50137089080467 0.047 +40 1.30292309231338 0.047 +50 1.13696207641347 0.047 +60 1.0290955710318 0.047 +70 0.926617553332013 0.047 +80 0.872197538842629 0.047 +90 0.805902827211741 0.047 +100 0.762288996264791 0.047 +110 0.735713914820247 0.047 +120 0.687503191709931 0.047 +130 0.659619799078658 0.047 +140 0.628200237768366 0.047 +150 0.610227714446274 0.047 +160 0.595197820960409 0.047 +180 0.562616568527284 0.047 +200 0.533018309573731 0.047 +250 0.491422841286475 0.047 +300 0.452001385193884 0.047 +350 0.433833355818465 0.047 +400 0.416712460291429 0.047 +Fragment 10; +0.75 0.346146147055066 0.047 +1.5 0.404248605191343 0.047 +1.7 0.426190820611571 0.047 +2.8 0.512684103166857 0.047 +5 0.671397806503662 0.047 +7 0.830628497205261 0.047 +8 0.917148437188257 0.047 +10 1.0681266249944 0.047 +13 1.28613529356722 0.047 +15 1.41025072440687 0.047 +80 1.17458249228727 0.047 +85 1.13390665507559 0.047 +90 1.09350611577311 0.047 +95 1.0655081455022 0.047 +100 1.04158596712227 0.047 \ No newline at end of file diff --git a/data/HSG_NanOx.db b/data/HSG_NanOx.db new file mode 100644 index 000000000..c70127047 --- /dev/null +++ b/data/HSG_NanOx.db @@ -0,0 +1,416 @@ +Fragment 1; +0.1 3.52785 0.0586794 +0.125 3.58379 0.0219491 +0.15 3.64192 0.0977552 +0.175 3.64134 0.045627 +0.2 3.59205 0.0522845 +0.225 3.48742 0.0763334 +0.25 3.38711 0.0486352 +0.275 3.23556 0.0140717 +0.3 3.10038 0.0564686 +0.325 2.92819 0.0736448 +0.35 2.74536 0.0440922 +0.375 2.64766 0.068873 +0.4 2.50822 0.0519665 +0.425 2.35826 0.0598118 +0.45 2.24049 0.0523686 +0.475 2.11282 0.0327819 +0.5 2.00902 0.0619912 +0.525 1.90411 0.048899 +0.55 1.81201 0.03943 +0.6 1.64405 0.0482093 +0.625 1.5553 0.0315918 +0.65 1.50294 0.0514883 +0.675 1.42643 0.0530718 +0.7 1.37338 0.0520192 +0.725 1.30831 0.0516787 +0.75 1.27779 0.0520332 +0.775 1.22202 0.0610424 +0.8 1.17571 0.0540186 +0.825 1.14052 0.0496843 +0.85 1.10919 0.0522072 +0.875 1.07663 0.0538031 +0.9 1.04062 0.0541082 +0.925 1.0123 0.0436477 +0.95 0.987964 0.0546012 +0.975 0.954701 0.052234 +0.9875 0.956628 0.0588447 +1 0.931771 0.058776 +1.25 0.922846 0.0596302 +1.375 0.829386 0.057546 +1.4375 0.811162 0.0635889 +1.5 0.77844 0.0576283 +1.5625 0.743542 0.0617835 +1.625 0.724542 0.0590762 +1.6875 0.705097 0.0629384 +1.75 0.690471 0.0610498 +1.875 0.657967 0.0689717 +2 0.630247 0.0609551 +2.125 0.59985 0.0558018 +2.25 0.576881 0.0605689 +2.375 0.560703 0.0620355 +2.5 0.548153 0.0558827 +2.75 0.521374 0.0553158 +3 0.49762 0.0574269 +3.25 0.487405 0.0658839 +3.5 0.467097 0.0576851 +3.875 0.459079 0.0697264 +4.25 0.436973 0.0635334 +4.625 0.425273 0.0644498 +5 0.420139 0.0625937 +6 0.406463 0.0625907 +7 0.389055 0.0639294 +7.5 0.390886 0.0518132 +8 0.377757 0.0648603 +9 0.380409 0.0416042 +10 0.375543 0.0629273 +13 0.36375 0.0637803 +14 0.355133 0.0663663 +14.5 0.363962 0.039987 +15 0.36261 0.0653386 +16 0.36199 0.0194867 +17 0.358723 0.0634626 +18.5 0.354093 0.0679557 +20 0.361473 0.0655274 +22.5 0.349285 0.0663717 +25 0.348734 0.0648933 +30 0.340981 0.0686984 +35 0.338952 0.0647853 +40 0.33787 0.0675714 +42.5 0.336445 0.0712075 +45 0.344386 0.0727018 +50 0.344973 0.0997508 +60 0.341337 0.0779679 +70 0.330335 0.0788438 +72.5 0.339633 0.083321 +75 0.341254 0.0831872 +80 0.343223 0.088238 +85 0.347429 0.0845721 +87.5 0.340158 0.0884021 +90 0.343359 0.089741 +100 0.3399 0.0930667 +110 0.344508 0.0972525 +115 0.341763 0.010713 +120 0.33716 0.093344 +125 0.348853 0.096485 +130 0.348952 0.0958178 +132.5 0.342216 0.108084 +135 0.339428 0.109774 +140 0.337172 0.0978914 +160 0.345544 0.101325 +165 0.341937 0.0939792 +170 0.351141 0.0839463 +175 0.346595 0.0956681 +180 0.331431 0.100256 +190 0.334884 0.104186 +200 0.338724 0.105441 +212.5 0.333464 0.094867 +225 0.33002 0.129717 +237.5 0.342825 0.0705913 +250 0.346535 0.109787 +275 0.339406 0.0998982 +300 0.339219 0.108748 +Fragment 2; +0.1 1.35471 0.0207273 +0.115 1.34943 0.068132 +0.125 1.36143 0.0561437 +0.135 1.36702 0.0479669 +0.15 1.40449 0.0523761 +0.175 1.46176 0.0186978 +0.2 1.52903 0.0422706 +0.225 1.54966 0.0685532 +0.25 1.55482 0.0423817 +0.275 1.62378 0.0624645 +0.3 1.69296 0.0560868 +0.325 1.75868 0.0827488 +0.35 1.82571 0.09662 +0.375 1.89151 0.0417015 +0.4 1.95804 0.0406013 +0.425 2.01788 0.0406869 +0.45 2.07856 0.0427398 +0.475 2.13803 0.0509513 +0.5 2.19514 0.0360138 +0.525 2.24906 0.0560736 +0.55 2.30165 0.0719733 +0.56 2.31162 0.0589844 +0.575 2.35233 0.0283324 +0.58 2.3359 0.0610067 +0.6 2.24179 0.0365905 +0.625 2.27969 0.048764 +0.65 2.32832 0.0571975 +0.675 2.36932 0.0582657 +0.7 2.40823 0.0107764 +0.725 2.44624 0.0803496 +0.75 2.47837 0.0209884 +0.775 2.51364 0.0847122 +0.8 2.54292 0.0519947 +0.825 2.57208 0.0609218 +0.85 2.59087 0.0639041 +0.875 2.61384 0.0825088 +0.9 2.64201 0.0602317 +0.925 2.65902 0.0399342 +0.95 2.67573 0.0100122 +0.975 2.68649 0.0900614 +1 2.69945 0.0460438 +1.125 2.72759 0.0649546 +1.25 2.68623 0.0592119 +1.3125 2.67842 0.0437095 +1.375 2.72147 0.0511202 +1.4375 2.6023 0.0481748 +1.5 2.55898 0.0135152 +1.5625 2.51015 0.0251847 +1.625 2.51876 0.0820981 +1.6875 2.40387 0.0796317 +1.75 2.349 0.0442119 +1.8 2.30672 0.0467401 +1.875 2.28272 0.0412346 +1.9 2.23023 0.0968624 +2 2.12334 0.0421089 +2.125 2.05126 0.0369635 +2.25 1.91357 0.0728045 +2.375 1.82345 0.0653356 +2.5 1.73184 0.0416488 +2.75 1.562 0.0705918 +3 1.43353 0.0454924 +3.25 1.31698 0.0449153 +3.5 1.21521 0.0409911 +3.75 1.13259 0.0515363 +3.875 1.09196 0.0562185 +4.25 0.99491 0.0515848 +4.625 0.913462 0.074134 +4.8 0.894276 0.0399684 +5 0.845962 0.0545597 +5.5 0.788932 0.0548308 +6 0.729322 0.0589305 +7 0.65191 0.0584992 +7.5 0.618191 0.0546201 +8 0.594936 0.0590312 +9 0.552776 0.0556349 +9.5 0.537816 0.0270089 +10 0.521669 0.0646037 +12 0.492294 0.0631951 +13 0.471319 0.0600745 +14 0.457 0.0540569 +14.5 0.454425 0.0647896 +15 0.447441 0.0639881 +16 0.442108 0.0585427 +17 0.427663 0.054866 +18.5 0.416323 0.0629227 +20 0.408267 0.0635374 +22.5 0.405099 0.0675533 +25 0.390136 0.0661453 +30 0.377674 0.0652145 +35 0.375329 0.0698207 +37.5 0.376945 0.068874 +40 0.370289 0.0662104 +42.5 0.372768 0.0683548 +45 0.364365 0.0703231 +50 0.364275 0.0681571 +60 0.355372 0.0747071 +65 0.355306 0.0760348 +70 0.357653 0.0778917 +72.5 0.351378 0.0813151 +75 0.354749 0.0843505 +80 0.356145 0.0820511 +85 0.351307 0.087815 +87.5 0.359332 0.0900896 +90 0.353549 0.0857066 +100 0.352662 0.0869735 +110 0.350016 0.0897816 +120 0.352662 0.0924707 +130 0.349242 0.0993062 +132.5 0.350188 0.0977779 +135 0.351315 0.101791 +140 0.351914 0.0996269 +160 0.352059 0.103381 +165 0.348055 0.108882 +170 0.349539 0.104125 +175 0.349366 0.097993 +180 0.350419 0.104505 +185 0.346769 0.107044 +190 0.342395 0.0732676 +195 0.349724 0.100912 +200 0.346889 0.116881 +212.5 0.348807 0.116918 +225 0.347231 0.0997025 +237.5 0.353986 0.110258 +250 0.34613 0.112309 +275 0.34762 0.107218 +300 0.347978 0.101295 +400 0.344634 0.102051 +500 0.34295 0.109641 +600 0.345637 0.10687 +700 0.342868 0.115623 +800 0.341284 0.0974119 +900 0.3453 0.103232 +1000 0.340597 0.103002 +Fragment 6; +0.1 0.554922 0.0174572 +0.15 0.507327 0.0501012 +0.175 0.496762 0.0185161 +0.2 0.498712 0.0128507 +0.225 0.507064 0.0398217 +0.25 0.517163 0.0321235 +0.275 0.524741 0.028797 +0.3 0.527277 0.024745 +0.35 0.562364 0.0325134 +0.375 0.581553 0.04897 +0.4 0.600128 0.022651 +0.45 0.63878 0.0493113 +0.475 0.658849 0.0709208 +0.5 0.678314 0.0311766 +0.525 0.651331 0.0378454 +0.53 0.654257 0.0387549 +0.55 0.668569 0.0323383 +0.575 0.685794 0.0958369 +0.6 0.698747 0.0925473 +0.625 0.71445 0.0123321 +0.65 0.729095 0.018834 +0.675 0.743807 0.033274 +0.7 0.757473 0.032298 +0.725 0.771899 0.029902 +0.75 0.7715 0.0258204 +0.775 0.780193 0.0222819 +0.8 0.792595 0.0115817 +0.825 0.804823 0.0254497 +0.85 0.818199 0.0192019 +0.875 0.828681 0.019282 +0.9 0.838734 0.0572171 +0.925 0.847032 0.0501417 +0.95 0.857266 0.0462406 +0.975 0.868805 0.0614917 +0.9875 0.875099 0.0422347 +1 0.875006 0.08374 +1.125 0.918784 0.0602353 +1.25 0.920182 0.0925326 +1.3125 0.937178 0.0795297 +1.375 0.949605 0.0347631 +1.4375 0.964492 0.0309748 +1.5 0.981289 0.0190316 +1.5625 0.993225 0.0224706 +1.625 1.00636 0.0147146 +1.6875 1.0206 0.039587 +1.75 1.0308 0.0626443 +1.875 1.03536 0.0595288 +2 1.03757 0.140597 +2.125 1.05661 0.0580031 +2.25 1.07951 0.047584 +2.5 1.11612 0.0134392 +2.75 1.15096 0.0417291 +3 1.18637 0.0549517 +3.25 1.2192 0.0958209 +3.5 1.25362 0.034796 +3.875 1.27999 0.0602744 +4.25 1.3268 0.0283756 +4.625 1.37202 0.0533761 +5 1.41833 0.0934593 +6 1.53289 0.0576861 +7 1.6441 0.0352824 +7.5 1.6905 0.0484951 +8 1.73733 0.0206558 +9 1.81111 0.0270153 +10 1.86519 0.0145637 +13 1.9013 0.0735242 +14 1.88 0.0524471 +14.5 1.86598 0.0753743 +15 1.84803 0.0247954 +16 1.80727 0.0386726 +17 1.75589 0.127889 +18.5 1.6866 0.0289437 +20 1.59615 0.0269112 +22.5 1.46482 0.02547 +25 1.3439 0.0377534 +27 1.25757 0.0699588 +30 1.14384 0.063285 +32 1.08349 0.0503975 +35 0.995382 0.0563748 +37.5 0.937827 0.0468086 +40 0.883756 0.0536156 +42.5 0.843352 0.06574 +45 0.806053 0.060595 +50 0.742661 0.0559711 +60 0.656602 0.0738968 +65 0.622154 0.0275371 +70 0.587718 0.0599057 +72.5 0.581034 0.0643803 +75 0.579367 0.0263508 +80 0.570936 0.0643563 +85 0.557118 0.0420902 +87.5 0.55305 0.064244 +90 0.549996 0.060972 +100 0.529454 0.044553 +110 0.514361 0.0641633 +120 0.500774 0.0656887 +130 0.490693 0.0572786 +132.5 0.486477 0.0687478 +135 0.484113 0.0336333 +140 0.479973 0.0652569 +150 0.471674 0.0623355 +160 0.464491 0.0678329 +165 0.461521 0.0694929 +175 0.455044 0.0702457 +180 0.454354 0.0690091 +185 0.450706 0.0807745 +190 0.446442 0.0477704 +195 0.445182 0.0765487 +200 0.443174 0.075495 +212.5 0.436772 0.0750927 +225 0.432475 0.0477424 +237.5 0.429575 0.0770921 +250 0.428672 0.0863672 +275 0.419821 0.0642036 +300 0.414074 0.0833777 +Fragment 8; +0.1 0.413736 0.0142512 +0.2 0.389411 0.0289218 +0.25 0.390787 0.0132718 +0.3 0.408917 0.0151436 +0.35 0.432501 0.0102363 +0.4 0.457564 0.0140088 +0.45 0.487057 0.0416276 +0.5 0.512859 0.029806 +0.55 0.497537 0.0513708 +0.6 0.519719 0.023812 +0.7 0.561975 0.0147102 +0.75 0.585803 0.0239479 +0.8 0.608447 0.00467075 +0.85 0.618215 0.0265444 +0.9 0.636862 0.0320603 +0.95 0.65226 0.0291653 +1 0.670432 0.0378584 +1.25 0.710918 0.0462263 +1.5 0.7642 0.0888616 +1.75 0.794974 0.0780676 +2.25 0.836185 0.0696362 +2.5 0.864835 0.0403956 +2.75 0.880924 0.0148689 +3 0.894613 0.00887251 +3.5 0.941321 0.0402146 +5 1.04036 0.0320689 +6 1.09991 0.0179879 +7 1.17086 0.0553748 +8 1.23791 0.03446 +10 1.36301 0.0780824 +13 1.53735 0.0380272 +15 1.63487 0.0208637 +20 1.77359 0.0114724 +25 1.78011 0.0369413 +30 1.69153 0.0243579 +40 1.43087 0.0687465 +50 1.20946 0.0680245 +60 1.044 0.0593522 +70 0.93624 0.0608797 +80 0.854476 0.0671706 +90 0.803408 0.0549333 +100 0.758682 0.0548792 +110 0.720772 0.0656161 +120 0.689713 0.0541748 +130 0.659153 0.0666365 +140 0.641074 0.0580988 +150 0.621408 0.0605031 +160 0.602909 0.0606604 +180 0.576667 0.0703917 +200 0.556955 0.0680523 +250 0.517284 0.0828523 +300 0.493476 0.0861731 +350 0.47441 0.0914356 diff --git a/docs/tools_to_interact_with_the_simulation_actors.rst b/docs/tools_to_interact_with_the_simulation_actors.rst index c052eea73..c1314cd3a 100644 --- a/docs/tools_to_interact_with_the_simulation_actors.rst +++ b/docs/tools_to_interact_with_the_simulation_actors.rst @@ -242,6 +242,90 @@ The output of the TetMeshDoseActor is a csv-file tabulating the results, e.g.:: Each row corresponds to one tetrahedron. The region marker column identifies to which macroscopic structure a tetrahedron belongs to -- it is equal to the region attribute defined for this tetrahedron in the '.ele' file the TetMeshBox is constructed from. +.. _biodose_measurement_biodoseactor-label: + +Biological dose measurement (BioDoseActor) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The BioDoseActor builds 3D images of:: + + the energy deposited (edep) + physical and biological doses deposited + alphamix and sqrtbetamix deposited + the number of hits + +in a given box volume. +The biological dose is computed as:: + + :math:`\frac{-\alpha_{ref} + \sqrt{\alpha_{ref}^2 + 4\,\beta_{ref}\,({\alpha_{mix}}\,Z + ({\sqrt{\beta_{mix}}}\,Z)^2 )}}{2\,\beta_{ref}}` + +It takes into account the weight of particles. +It can store multiple information into a 3D grid, each information can be enabled by using:: + + /gate/actor/[Actor Name]/enableEdep true + /gate/actor/[Actor Name]/enableDose true + /gate/actor/[Actor Name]/enableAlphaMix true + /gate/actor/[Actor Name]/enableSqrtBetaMix true + /gate/actor/[Actor Name]/enableRBE true + /gate/actor/[Actor Name]/enableUncertainty true + /gate/actor/[Actor Name]/enableUncertaintyDetails true + /gate/actor/[Actor Name]/enableHitEventCount true + +Informations can be disabled by using "false" (default) instead of "true". +The unit of edep is MeV, the unit of dose is Gy, the unit of biological dose is Gy RBE. +The uncertainty outputs are relative statistical uncertainties. + +For the output, the suffixes:: + + _edep + _dose + _alphamix + _sqrtbetamix + _biodose + _rbe + _hitevent_count + _dose_uncertainty + _biodose_uncertainty + +are added to the output file name given by the user. + +The user must provide a cell line and a biophysical model:: + + /gate/actor/[Actor Name]/setCellLine HSG + /gate/actor/[Actor Name]/setBioPhysicalModel NanOx + +which will be combined to make a filename: data/{CellLine}_{BioPhysicalModel}.db, for example data/HSG_NanOx.db. +Only HSG (Human Salivary Gland) cell line is provided. +Two biophysical models are available: NanOx (carbon, proton) and MMKM (carbon). +User can provide its own data files (cell line associated to a biophysical model database). + +The user must also provide the alpha and beta reference values corresponding to the cell line:: + + /gate/actor/[Actor Name]/setAlphaRef 0.313 + /gate/actor/[Actor Name]/setBetaRef 0.0615 + +If the user wants to reach a large dose value output to +The user can use a dose scale factor in order to reach a large dose value without needing to run higher number of primary particles. +The dose scale factor is set with this command (default 1):: + + /gate/actor/[Actor Name]/setDoseScaleFactor 1e3 + +It will multiply the dose value (and affect the biological dose value). +Note that it will not affect the relative uncertainty, so it can be used once some uncertainty threshold has been reached. + +Usage example:: + + /gate/actor/addActor BioDoseActor MyBio + /gate/actor/MyBio/attachTo MyVolume + /gate/actor/MyBio/setVoxelSize 200 200 1 mm + /gate/actor/MyBio/setPosition 0 0 0 + /gate/actor/MyBio/setCellLine HSG + /gate/actor/MyBio/setBioPhysicalModel NanOx + /gate/actor/MyBio/setAlphaRef 0.313 + /gate/actor/MyBio/setBetaRef 0.0615 + /gate/actor/MyBio/enableDose true + /gate/actor/MyBio/save output/{particleName}.mhd + .. _kill_track-label: Kill track diff --git a/source/digits_hits/include/GateBioDoseActor.hh b/source/digits_hits/include/GateBioDoseActor.hh index 9bc150ddb..5bf161b82 100644 --- a/source/digits_hits/include/GateBioDoseActor.hh +++ b/source/digits_hits/include/GateBioDoseActor.hh @@ -138,23 +138,23 @@ private: GateImageWithStatistic _alphaMixDoseImage; GateImageWithStatistic _sqrtBetaMixDoseImage; - GateImageWithStatistic _pdBiodoseAlphaMixMean; - GateImageWithStatistic _pdBiodoseSqrtBetaMixMean; - GateImageWithStatistic _pdBiodoseDoseMean; + GateImageWithStatistic _pdBiodoseAlphaMixMeanImage; + GateImageWithStatistic _pdBiodoseSqrtBetaMixMeanImage; + GateImageWithStatistic _pdBiodoseDoseMeanImage; GateImageWithStatistic _varAlphaMixMeanImage; GateImageWithStatistic _varSqrtBetaMixMeanImage; GateImageWithStatistic _varDoseMeanImage; GateImageWithStatistic _covAlphaMixMeanSqrtBetaMixMeanImage; - GateImageWithStatistic _covAlphaMixMeanDoseMean; - GateImageWithStatistic _covSqrtBetaMixMeanDoseMean; + GateImageWithStatistic _covAlphaMixMeanDoseMeanImage; + GateImageWithStatistic _covSqrtBetaMixMeanDoseMeanImage; // Outputs bool _enableEdep = false; - bool _enableDose = true; + bool _enableDose = false; bool _enableAlphaMix = false; bool _enableSqrtBetaMix = false; bool _enableRBE = false; - bool _enableUncertainty = true; + bool _enableUncertainty = false; bool _enableUncertaintyDetails = false; bool _enableHitEventCount = false; diff --git a/source/digits_hits/src/GateBioDoseActor.cc b/source/digits_hits/src/GateBioDoseActor.cc index b6191ff06..e8ba9bc72 100644 --- a/source/digits_hits/src/GateBioDoseActor.cc +++ b/source/digits_hits/src/GateBioDoseActor.cc @@ -80,15 +80,15 @@ void GateBioDoseActor::Construct() { setupImage(_sqrtBetaMixDoseImage); if(_enableUncertaintyDetails) { - setupImage(_pdBiodoseAlphaMixMean, "pd_biodose_alphamixmean"); - setupImage(_pdBiodoseSqrtBetaMixMean, "pd_biodose_sqrtbetamixmean"); - setupImage(_pdBiodoseDoseMean, "pd_biodose_dosemean"); + setupImage(_pdBiodoseAlphaMixMeanImage, "pd_biodose_alphamixmean"); + setupImage(_pdBiodoseSqrtBetaMixMeanImage, "pd_biodose_sqrtbetamixmean"); + setupImage(_pdBiodoseDoseMeanImage, "pd_biodose_dosemean"); setupImage(_varAlphaMixMeanImage, "var_alphamixmean"); setupImage(_varSqrtBetaMixMeanImage, "var_sqrtbetamixmean"); setupImage(_varDoseMeanImage, "var_dosemean"); setupImage(_covAlphaMixMeanSqrtBetaMixMeanImage, "cov_alphamixmean_sqrtbetamixmean"); - setupImage(_covAlphaMixMeanDoseMean, "cov_alphamixmean_dosemean"); - setupImage(_covSqrtBetaMixMeanDoseMean, "cov_sqrtbetamixmean_dosemean"); + setupImage(_covAlphaMixMeanDoseMeanImage, "cov_alphamixmean_dosemean"); + setupImage(_covSqrtBetaMixMeanDoseMeanImage, "cov_sqrtbetamixmean_dosemean"); } } } @@ -363,15 +363,15 @@ void GateBioDoseActor::updateData() { ) / biodoseMean; if(_enableUncertaintyDetails) { - _pdBiodoseAlphaMixMean.SetValue(index, pdBiodoseAlphaMixMean); - _pdBiodoseSqrtBetaMixMean.SetValue(index, pdBiodoseSqrtBetaMixMean); - _pdBiodoseDoseMean.SetValue(index, pdBiodoseDoseMean); + _pdBiodoseAlphaMixMeanImage.SetValue(index, pdBiodoseAlphaMixMean); + _pdBiodoseSqrtBetaMixMeanImage.SetValue(index, pdBiodoseSqrtBetaMixMean); + _pdBiodoseDoseMeanImage.SetValue(index, pdBiodoseDoseMean); _varAlphaMixMeanImage.SetValue(index, varAlphaMixMean); _varSqrtBetaMixMeanImage.SetValue(index, varSqrtBetaMixMean); _varDoseMeanImage.SetValue(index, varDoseMean); _covAlphaMixMeanSqrtBetaMixMeanImage.SetValue(index, covAlphaMixMeanSqrtBetaMixMean); - _covAlphaMixMeanDoseMean.SetValue(index, covAlphaMixMeanDoseMean); - _covSqrtBetaMixMeanDoseMean.SetValue(index, covSqrtBetaMixMeanDoseMean); + _covAlphaMixMeanDoseMeanImage.SetValue(index, covAlphaMixMeanDoseMean); + _covSqrtBetaMixMeanDoseMeanImage.SetValue(index, covSqrtBetaMixMeanDoseMean); } _doseUncertaintyImage.SetValue(index, uncertaintyDose); @@ -408,15 +408,15 @@ void GateBioDoseActor::SaveData() { _biodoseUncertaintyImage.SaveData(_currentEvent); if(_enableUncertaintyDetails) { - _pdBiodoseAlphaMixMean.SaveData(_currentEvent); - _pdBiodoseSqrtBetaMixMean.SaveData(_currentEvent); - _pdBiodoseDoseMean.SaveData(_currentEvent); + _pdBiodoseAlphaMixMeanImage.SaveData(_currentEvent); + _pdBiodoseSqrtBetaMixMeanImage.SaveData(_currentEvent); + _pdBiodoseDoseMeanImage.SaveData(_currentEvent); _varAlphaMixMeanImage.SaveData(_currentEvent); _varSqrtBetaMixMeanImage.SaveData(_currentEvent); _varDoseMeanImage.SaveData(_currentEvent); _covAlphaMixMeanSqrtBetaMixMeanImage.SaveData(_currentEvent); - _covAlphaMixMeanDoseMean.SaveData(_currentEvent); - _covSqrtBetaMixMeanDoseMean.SaveData(_currentEvent); + _covAlphaMixMeanDoseMeanImage.SaveData(_currentEvent); + _covSqrtBetaMixMeanDoseMeanImage.SaveData(_currentEvent); } } if(_enableHitEventCount) _hitEventCountImage.SaveData(_currentEvent); @@ -450,15 +450,15 @@ void GateBioDoseActor::ResetData() { _sqrtBetaMixDoseImage.Reset(); if(_enableUncertaintyDetails) { - _pdBiodoseAlphaMixMean.Reset(); - _pdBiodoseSqrtBetaMixMean.Reset(); - _pdBiodoseDoseMean.Reset(); + _pdBiodoseAlphaMixMeanImage.Reset(); + _pdBiodoseSqrtBetaMixMeanImage.Reset(); + _pdBiodoseDoseMeanImage.Reset(); _varAlphaMixMeanImage.Reset(); _varSqrtBetaMixMeanImage.Reset(); _varDoseMeanImage.Reset(); _covAlphaMixMeanSqrtBetaMixMeanImage.Reset(); - _covAlphaMixMeanDoseMean.Reset(); - _covSqrtBetaMixMeanDoseMean.Reset(); + _covAlphaMixMeanDoseMeanImage.Reset(); + _covSqrtBetaMixMeanDoseMeanImage.Reset(); } } } From 431dddafa5eed5881fa8122778f13790288f462c Mon Sep 17 00:00:00 2001 From: tbaudier Date: Wed, 27 Mar 2024 18:43:35 +0100 Subject: [PATCH 12/12] Add t34_biodose_actor test --- .github/workflows/main.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index ce0a1ba64..bac6a8dd1 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -283,7 +283,8 @@ jobs: t30_dna, t31_vpgTLE-tt, t32_isotopes, - t33_invert_filter] + t33_invert_filter, + t34_biodose_actor] steps: - name: Checkout github repo