Theoretical Issue in Dosimetry

Greetings.

I used to perform dosimetry and energy deposition measurements using sensitive detectors. But since I found command-line macro easier, I decided to use that also to compare my results. I was expecting same results, but they are totally different.

I have an object as a mother volume. Inside of this package, there are lots of daughter volumes like water, plastic, glass, etc. I know the mass of the total package; it is about 2.41 kg. In my sensitive detector code, I used to calculate the total absorbed energy in MeV and then convert it to Joules. Then, I manually divided this value by the mass of the whole package to get dose in Gy. So,

Dose = Total Energy Deposition (J) / Total mass (kg)

With this method I got a dose of 6.34406e-8 Gy.

But when I used command-line and attached the scoring to the logical volumes as /score/create/realWorldLogVol <name of my logical volume> for each logical volume in the package, the calculated dose for the whole package is 1.522924e-5 Gy which is 1000 times greater. The total energy deposition from both code-based and command-line methods are kind of similar but the dose is different from my calculations.

So, I was wondering if I did not understand the idea of dosimetry properly. Is my manual division theoretically correct? The total energy deposition using code-based is 0.97 TeV, command-line is 0.511 TeV (Not same run! DIFFERENT RUNS). If I consider the calculated dose from command-line method, using the dose equation I provided, the mass of the whole package would be 5.38 g!!

I am very confused. Please assist me.

Do you have many daughter volumes? Scoring in that way will score each daughter volume with the same name separately. You also cannot use the mother volume for this. See slides 15 and 16 here.

Hi jrellin,

Please let me walk you through the process, forgive me if this is too much. Allow me to share my codes with you. The following code is for my SensitiveDetector.cc,

#include "SensitiveDetector.hh"

//Initialize static variable (this must remain outside of constructor)
G4double SensitiveDetector::fRunTotalEnergy = 0.;

SensitiveDetector::SensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{
    fTotalEnergy = 0.;
}

SensitiveDetector::~SensitiveDetector()
{
}

void SensitiveDetector::Initialize(G4HCofThisEvent*)
{
    // Reset for new event
    fTotalEnergy = 0.;

    // Resets run total energy at start of first event.
    G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    if(eventID == 0) {
        fRunTotalEnergy = 0.;
    }
}

G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*)
{
    // Get energy deposited in this step
    G4double edep = step->GetTotalEnergyDeposit();

    if(edep>0)
    {
    // Add to fTotalEnergy
    fTotalEnergy += edep;
    }
    
    return true;
}

void SensitiveDetector::EndOfEvent(G4HCofThisEvent*)
{

    //Accumulate to run total energy and check; if last event, print the result.
    fRunTotalEnergy += fTotalEnergy;
    
    G4int currentEvent = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    G4int totalEvents = G4RunManager::GetRunManager()->GetNumberOfEventsToBeProcessed();

    if (currentEvent == totalEvents - 1) {
    // Only print from the container detector
    if (GetName() == "ContainerLVDetector") {
        G4cout << "=== RUN COMPLETE ===" << G4endl;
        G4cout << "Total Events: " << totalEvents << G4endl;
        G4cout << "Total Energy Deposited in Run: " << fRunTotalEnergy << " MeV" << G4endl;
        G4cout << "===================" << G4endl;
    }
}

}

I attached this code to the ConstructSDandField in my DetectorConstruction as such,

void DetectorConstruction::ConstructSDandField()
{
    //SD manager must be created here only once
    G4SDManager* SDman = G4SDManager::GetSDMpointer();
    
    //Sensitive detector for containerLV
    SensitiveDetector *detector_containerLV = new SensitiveDetector("ContainerLVDetector");
    containerLV->SetSensitiveDetector(detector_containerLV);
    SDman->AddNewDetector(detector_containerLV);

    //Sensitive detector for Cardboard
    SensitiveDetector *detector_shellLogic = new SensitiveDetector("ShellLogicDetector");
    shellLogic->SetSensitiveDetector(detector_shellLogic);
    SDman->AddNewDetector(detector_shellLogic);

    //Sensitive detector for Foam
    SensitiveDetector *detector_foamLogic = new SensitiveDetector("FoamLogicDetector");
    foamLogic->SetSensitiveDetector(detector_foamLogic);
    SDman->AddNewDetector(detector_foamLogic);

    //Sensitive detector for Vial
    SensitiveDetector *detector_vialLV = new SensitiveDetector("VialLVDetector");
    vialLV->SetSensitiveDetector(detector_vialLV);
    SDman->AddNewDetector(detector_vialLV);

    //Sensitive detector for Water Phantom
    SensitiveDetector *detector_logicWater = new SensitiveDetector("LogicWaterDetector");
    logicWater->SetSensitiveDetector(detector_logicWater);
    SDman->AddNewDetector(detector_logicWater);
    
}

As you have already anticipated, I have many daughter volumes (2*256 volumes) that are replicated but the whole package consists of 5 logical volumes as above. This is the way I assigned my sensitive detector to the logical volumes and for every run, I always get 0.97 TeV energy deposition.

At the same time, I tried to examine the same thing using command-line and here is what I wrote,

# Score the air container (single volume)
/score/create/realWorldLogVol ContainerLV
/score/quantity/energyDeposit containerEnergy
/score/close

# Score the cardboard shell (single volume) 
/score/create/realWorldLogVol ShellLogical
/score/quantity/energyDeposit shellEnergy
/score/close

# Score the foam (2 pieces - will sum automatically)
/score/create/realWorldLogVol FoamLogical
/score/quantity/energyDeposit foamEnergy
/score/close

# Score all vial glass (all vials combined)
/score/create/realWorldLogVol VialLogical
/score/quantity/energyDeposit vialEnergy
/score/close

# Score all water (all water volumes combined)
/score/create/realWorldLogVol WaterLogical
/score/quantity/energyDeposit waterEnergy
/score/close

# Run simulation
/run/beamOn 1000000

# Dump all results to files
/score/dumpQuantityToFile ContainerLV containerEnergy container_energy.txt
/score/dumpQuantityToFile ShellLogical shellEnergy shell_energy.txt  
/score/dumpQuantityToFile FoamLogical foamEnergy foam_energy.txt
/score/dumpQuantityToFile VialLogical vialEnergy vial_energy.txt
/score/dumpQuantityToFile WaterLogical waterEnergy water_energy.txt

# Print summary to console
/score/list

And the reported results are OK; I know that I have 256 water samples inside the package and the text file for water_energy.txt has 256 rows, reporting the energy deposition of all the volumes. Same is true for other logical volumes I assigned.

But the funny thing is, when I sum up all the energy depositions, I get 0.511 TeV not 0.97 TeV. No matter how many times I run, I’ll always get the same results. I wonder what is the cause?

My Issue Regarding Dosimetry:

As I mentioned before, to calculate the dose, I used to convert the energy deposition I got from the simulation to Joules and divide it manually by the mass of my package (the mass calculation is done by Geant4 and I double-checked it manually, no problem).

Dose via my SensitiveDetector energy deposition report: 6.34e-8 Gy.
Dose via command-line energy deposition report: 3.34e-8 Gy

I decided to write another command-line for dose and this is the macro file,

/run/initialize

# Score the air container (single volume)
/score/create/realWorldLogVol ContainerLV
/score/quantity/doseDeposit containerDose
/score/close

# Score the cardboard shell (single volume) 
/score/create/realWorldLogVol ShellLogical
/score/quantity/doseDeposit shellDose
/score/close

# Score the foam (2 pieces - will sum automatically)
/score/create/realWorldLogVol FoamLogical
/score/quantity/doseDeposit foamDose
/score/close

# Score all vial glass (all vials combined)
/score/create/realWorldLogVol VialLogical
/score/quantity/doseDeposit vialDose
/score/close

# Score all water (all water volumes combined)
/score/create/realWorldLogVol WaterLogical
/score/quantity/doseDeposit waterDose
/score/close

# Run simulation
/run/beamOn 1000000

# Dump all results to files
/score/dumpQuantityToFile ContainerLV containerDose container_dose.txt
/score/dumpQuantityToFile ShellLogical shellDose shell_dose.txt  
/score/dumpQuantityToFile FoamLogical foamDose foam_dose.txt
/score/dumpQuantityToFile VialLogical vialDose vial_dose.txt
/score/dumpQuantityToFile WaterLogical waterDose water_dose.txt

# Print summary to console
/score/list

As you can see, same assignment, but this time for dose. This also gives peculiar results. The reported dose for the whole package is 1.52e-5!! this almost 1000 times larger. I am very confused and I don’t know which method is giving me the correct dose.

Please assist me.

Just as a point of clarification, are you running the SD at the same time as the scorer? I am pretty sure that behavior is unpredictable similar to trying to connect multiple SDs to the same volume.

No, not at the same time. Different runs. Actually, I tried to run both at the same time but somehow my SD reported 0 MeV energy absorption. I knew something is wrong, so I decided to take separate runs.

I would like to add an update. In my SD, I used to define the energy counter variable as static so that it adds up the energy deposition of all logical volumes I assigned. I thought, maybe this static variable is somehow over-counting the energy depositions. So I removed the static and allowed the SD to report the energy deposition on each logical volume separately.

The result is, ‘My SD was over-counting the energy depositions and now I have energy depositions similar to Scoring Manager results’. However, I found a systematic problem here:

Energy deposition results from SD:
ContainerLVDetector | 6.254381e+02 MeV
ShellLogicDetector | 1.023402e+04 MeV
FoamLogicDetector | 1.434552e+04 MeV
VialLVDetector | 1.937479e+05 MeV
LogicWaterDetector | 2.680398e+05 MeV

Energy deposition results from ScoringManager:
ContainerLV | 1.257005e+03 MeV
ShellLogical | 2.042153e+04 MeV
FoamLogical | 2.876773e+04 MeV
VialLogical | 1.939900e+05 MeV
WaterLogical | 2.673824e+05 MeV

As you can see, the results for Vial and Water are almost the same while the result for container, shell, and foam in ScoringManager is double the SD. There is a 2:1 ratio; this is a systematic problem. I still don’t know why this happens.

In addition, vial and water are cylindrical geometries while the rest are boxes.

This code is confusing:

  
    G4int currentEvent = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    G4int totalEvents = G4RunManager::GetRunManager()->GetNumberOfEventsToBeProcessed();

    if (currentEvent == totalEvents - 1) {
    // Only print from the container detector
    if (GetName() == "ContainerLVDetector") {
        G4cout << "=== RUN COMPLETE ===" << G4endl;
        G4cout << "Total Events: " << totalEvents << G4endl;
        G4cout << "Total Energy Deposited in Run: " << fRunTotalEnergy << " MeV" << G4endl;
        G4cout << "===================" << G4endl;
    }
}

This really should be printed inside the EndOfRun class. The processing order of events, tracks, and steps should not be assumed “in order” which would be compounded if you tried to use mulithreading, for example.

To eliminate that as a confounding cause, I’d refactor this into the End of Run similar to what is done with the “fEDepTarget” variable in rdecay02.

Dear @jrellin,

Thank you for the suggestion. I decided to take your approach. And here are the changes I made in my SD,

This is the SensitiveDetector.cc,

#include "SensitiveDetector.hh"
#include <iomanip>

SensitiveDetector::SensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{
    fTotalEnergy = 0.;
    fRunTotalEnergy = 0.;
}

SensitiveDetector::~SensitiveDetector()
{
}

void SensitiveDetector::Initialize(G4HCofThisEvent*)
{
    fTotalEnergy = 0.;  // Reset for new event
    
    G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    if(eventID == 0) {
        fRunTotalEnergy = 0.;  // Reset run total only once
    }
}

G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*)
{
    G4double edep = step->GetTotalEnergyDeposit();
    fTotalEnergy += edep;  // Accumulate for this event
    return true;
}

void SensitiveDetector::EndOfEvent(G4HCofThisEvent*)
{
    fRunTotalEnergy += fTotalEnergy;  //Just transfer to run total
    // NO RESET NEEDED HERE
}

Also, I added this line in my SensitiveDetector.hh in the public section,

G4double GetRunTotalEnergy() const { return fRunTotalEnergy; }

Then, I made my RunAction.cc and RunAction.hh as such,
This is RunAction.hh,

#ifndef RunAction_h
#define RunAction_h 1

#include "G4UserRunAction.hh"
#include "G4Run.hh"
#include "G4SDManager.hh"

class G4Run;

class RunAction : public G4UserRunAction
{
public:
    RunAction();
    virtual ~RunAction();

    virtual void BeginOfRunAction(const G4Run*);
    virtual void EndOfRunAction(const G4Run*);

private:
};

#endif

This is RunAction.cc,

#include "RunAction.hh"
#include "SensitiveDetector.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"
#include <iomanip>

RunAction::RunAction()
: G4UserRunAction()
{
    // Constructor
}

RunAction::~RunAction()
{
    // Destructor
}

void RunAction::BeginOfRunAction(const G4Run* run)
{
    G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
}

void RunAction::EndOfRunAction(const G4Run* run)
{
    G4int nofEvents = run->GetNumberOfEvent();
    if (nofEvents == 0) return;

    G4cout << "\n=== FINAL SD ENERGY RESULTS (Run " << run->GetRunID() << ") ===" << G4endl;
    G4cout << "Total Events: " << nofEvents << G4endl;
    G4cout << "=============================================" << G4endl;
    
    G4SDManager* SDman = G4SDManager::GetSDMpointer();
    
    // Container SD
    SensitiveDetector* containerSD = dynamic_cast<SensitiveDetector*>(
        SDman->FindSensitiveDetector("ContainerLVDetector"));
    if(containerSD) {
        G4cout << containerSD->GetName() << " Total Energy Deposition: " 
               << std::scientific << std::setprecision(6) 
               << containerSD->GetRunTotalEnergy()/CLHEP::MeV << " MeV" << G4endl;
    }
    
    // Shell SD
    SensitiveDetector* shellSD = dynamic_cast<SensitiveDetector*>(
        SDman->FindSensitiveDetector("ShellLogicDetector"));
    if(shellSD) {
        G4cout << shellSD->GetName() << " Total Energy Deposition: " 
               << std::scientific << std::setprecision(6) 
               << shellSD->GetRunTotalEnergy()/CLHEP::MeV << " MeV" << G4endl;
    }
    
    // Foam SD
    SensitiveDetector* foamSD = dynamic_cast<SensitiveDetector*>(
        SDman->FindSensitiveDetector("FoamLogicDetector"));
    if(foamSD) {
        G4cout << foamSD->GetName() << " Total Energy Deposition: " 
               << std::scientific << std::setprecision(6) 
               << foamSD->GetRunTotalEnergy()/CLHEP::MeV << " MeV" << G4endl;
    }
    
    // Vial SD
    SensitiveDetector* vialSD = dynamic_cast<SensitiveDetector*>(
        SDman->FindSensitiveDetector("VialLVDetector"));
    if(vialSD) {
        G4cout << vialSD->GetName() << " Total Energy Deposition: " 
               << std::scientific << std::setprecision(6) 
               << vialSD->GetRunTotalEnergy()/CLHEP::MeV << " MeV" << G4endl;
    }
    
    // Water SD
    SensitiveDetector* waterSD = dynamic_cast<SensitiveDetector*>(
        SDman->FindSensitiveDetector("LogicWaterDetector"));
    if(waterSD) {
        G4cout << waterSD->GetName() << " Total Energy Deposition: " 
               << std::scientific << std::setprecision(6) 
               << waterSD->GetRunTotalEnergy()/CLHEP::MeV << " MeV" << G4endl;
    }
    
    G4cout << "=============================================" << G4endl;
    
    G4cout << G4endl;
}

And because of Multi-threading, after I shot 1 million electrons, this is the result I got,

 SD ENERGY RESULTS (Run 0) ===
G4WT1 > Total Events: 499849
G4WT1 > =============================================
G4WT1 > ContainerLVDetector Total Energy Deposition: 6.318647e+02 MeV
G4WT1 > ShellLogicDetector Total Energy Deposition: 1.016732e+04 MeV
G4WT1 > FoamLogicDetector Total Energy Deposition: 1.437964e+04 MeV
G4WT1 > VialLVDetector Total Energy Deposition: 1.933064e+05 MeV
G4WT1 > LogicWaterDetector Total Energy Deposition: 2.665447e+05 MeV
G4WT1 > =============================================
G4WT1 > 
G4WT0 > 
=== FINAL SD ENERGY RESULTS (Run 0) ===
G4WT0 > Total Events: 500151
G4WT0 > =============================================
G4WT0 > ContainerLVDetector Total Energy Deposition: 6.322745e+02 MeV
G4WT0 > ShellLogicDetector Total Energy Deposition: 1.017765e+04 MeV
G4WT0 > FoamLogicDetector Total Energy Deposition: 1.440065e+04 MeV
G4WT0 > VialLVDetector Total Energy Deposition: 1.940219e+05 MeV
G4WT0 > LogicWaterDetector Total Energy Deposition: 2.669596e+05 MeV
G4WT0 > =============================================
G4WT0 > 

I manually added the energies but now we have another problem,

Combined Multithreading Results (My SD):

ContainerLV: 1.264139e+03 MeV

ShellLV: 2.034497e+04 MeV

FoamLV: 2.878029e+04 MeV

VialLV: 3.873283e+05 MeV

WaterLV: 5.335043e+05 MeV

Previous Results (Scoring Manager):

ContainerLV: 1.264139e+03 MeV

ShellLV: 2.034497e+04 MeV

FoamLV: 2.878029e+04 MeV

VialLV: 1.936261e+05 MeV

WaterLV: 2.666996e+05 MeV

Summary: Box geometries now match perfectly. But Cylindrical geometries show exactly 2× the expected values.

Something is wrong here that I cannot understand.

Do the bare minimal you need to swap the Vial and Water LV’s with roughly box equivalents to see if the problem persists. I haven’t used them extensively but IIRC some of these scorers do not behave well with objects that are not neatly cartesian like spheres and cylinders.

Thank you. Since the SD is giving me different results all the time, I decided to rely on command-line scoring.