Getting lower energy depostion than expected using isotope production


I am trying to score energy depostion in organs of human phantoms for uniform distribution of I-131 isotope via macro using commands,

Example macro file to run in the batch mode

Set the verbose

/run/verbose 2

Set the number of threads for multi-threading mode

/run/numberOfThreads 30


/gun/particle ion
/gun/ion 53 131 0 0

Set the nps

/run/beamOn 1000000

I am getting correct energy deposition in cross-irradiation cases (when source and target organs are different), However or slef-irrdiation cases energy depostion is somewhat lower (about 50% only) of what is expected.
I am using
// EM physics
RegisterPhysics(new G4EmLivermorePhysics());

// Decay
RegisterPhysics(new G4DecayPhysics());

// Radioactive decay
RegisterPhysics(new G4RadioactiveDecayPhysics());

physics and // Set the secondary production cut for photons to 1 µm
G4double photonCut = 1.0 * CLHEP::um; // 1 µm
SetCutValue(photonCut, “gamma”);

// Set the secondary production cut for electrons to 1 µm
G4double electronCut = 1.0 * CLHEP::um;  // 1 µm
SetCutValue(electronCut, "e-");

Production cut.
Any Idea what is wrong in this??


Hi I would suggest to use the advanced example ICRP145Phantom and check if you get the same results.

Hello @guatelli
I checked ICRP145Phantoms advance example and compared with ICRP provided GEANT4 files for same. In advance example internal source sampling scheme is not there also instead of G4EmLivermorePhysics QGSP_BIC_HP physics is used. Is there any other difference between these two files?

I agree you have to add the internal source sampling scheme, I would suggest to use the General Particle Source.There is difference between G4EMLivermorePhysics and QGSP_BIC_HP. Please, check the documentation of Geant4. The first models EM physics processes only. QGSP_BIC_HP is a complete physics list of Geant4, modelling both EM and hadronic physics interactions (please check the G4 doc).
If you model the radioactive decay, make sure to have the Radioactive Decay activated in your physics list. We recently published a paper on a project of internal dosimetry. Please, look here, it may have some useful information for you:


I used gps to set internal source sampling scheme in ICRP145phantoms example using same scheme as provided in ICRP example (that is via rejection method). Also physics I set is

// physics list
G4VModularPhysicsList* physicsList = new QGSP_BIC_HP(); // or another physics list
physicsList->RegisterPhysics(new G4DecayPhysics());
physicsList->RegisterPhysics(new G4RadioactiveDecayPhysics());

After this I am using macro to set particle



/gps/particle ion
/gps/ion 53 131
However the results are same. The energy deposited in self target case is still lower than expected.

How do you record the energy deposition?

This is SD and edep scorer in detector construction class,

void TETDetectorConstruction::ConstructSDandField()
 // Define detector (Phantom SD) and scorer (eDep)
 G4SDManager* pSDman = G4SDManager::GetSDMpointer();
 G4String phantomSDname = "PhantomSD";

 // MultiFunctional detector
 auto* MFDet = new G4MultiFunctionalDetector(phantomSDname);
 pSDman->AddNewDetector( MFDet );

 // scorer for energy depositon in each organ
 MFDet->RegisterPrimitive(new TETPSEnergyDeposit("eDep", fTetData));

 // attach the detector to logical volume for parameterised geometry (phantom geometry)
 SetSensitiveDetector(fTetLogic, MFDet);

The energy Deposition class is ,

#include "TETPSEnergyDeposit.hh"

TETPSEnergyDeposit::TETPSEnergyDeposit(G4String name, TETModelImport* _tetData)
  :G4PSEnergyDeposit(name), fTetData(_tetData)

G4int TETPSEnergyDeposit::GetIndex(G4Step* aStep)
  // return the organ ID (= material index)
  G4int copyNo = aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber();
  return fTetData->GetMaterialIndex(copyNo);

and finally;

event by event scoring of energy in organs is

void TETRun::RecordEvent(const G4Event* event)
  fCollID = G4SDManager::GetSDMpointer()->GetCollectionID("PhantomSD/eDep");

 // Hits collections
 G4HCofThisEvent* HCE = event->GetHCofThisEvent();
 if(!HCE) return;
 auto* evtMap = static_cast<G4THitsMap<G4double>*>(HCE->GetHC(fCollID));
 // sum up the energy deposition and the square of it
 for (auto itr : *evtMap->GetMap()) {
		fEdepMap[itr.first].first  += *itr.second;                   //sum
		fEdepMap[itr.first].second += (*itr.second) * (*itr.second); //sum square

The overall setup seems to work fine when simulating energy deposition for electrons and photons individually, However, for isotope decay setup problem arises…


I all predefined Physics Lists DecayPhysics is already included. In all HP PhysicsList the RadioactiveDecay is also included. You need to check if these processes are not instantiated twice.


1 Like

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.