Difference between hand calculation and Simulation result

HyeonSeok-Do · February 25, 2025, 3:33am

Dear geant4 and physics experts,

Hello,
I am calculating the depth at which 95% of 100 keV photons are absorbed when they are incident on a silicon target.
In other words, I’m determining the Si target depth required for 95% containment using 100keV mono energy photon.

First, I calculated it by my hands using “Exponential Attenuation law” and “nist dataset”.

.
Compare to this simple caculation result, I made a Geant4 code.

My environment is like this.
[50,50,50]m huge Si target.
I fired 100keV mono energy photon at [0,0,-30]m to z position.

G4PhysListFactory physicsFactory;
G4VModularPhysicsList* physicsList = physicsFactory.GetReferencePhysList("FTFP_BERT_EMZ");
runManager->SetUserInitialization(physicsList);

‘’’
void SteppingAction::UserSteppingAction(const G4Step* step) {
auto preStepPoint = step->GetPreStepPoint();
auto postStepPoint = step->GetPostStepPoint();

// Get position information
G4ThreeVector prePosition = preStepPoint->GetPosition();
G4ThreeVector postPosition = postStepPoint->GetPosition();

// Get volume information
const G4VPhysicalVolume* preVolume = preStepPoint->GetTouchableHandle()->GetVolume();
const G4VPhysicalVolume* postVolume = postStepPoint->GetTouchableHandle()->GetVolume();

// Only consider primary particle (TrackID == 1)
if (step->GetTrack()->GetTrackID() == 1){

    // Record the end step when particle is killed
    if (step->GetTrack()->GetTrackStatus() == fStopAndKill) {
        G4ThreeVector finalpostPosition = (postVolume && postVolume->GetName() == "Silicon") ? postPosition : G4ThreeVector(-999,-999,-9999);

        // position, detectorID 저장 및 출력
        fEventAction->RecordFinalStep(finalpostPosition);
    }

}

I used this SteppingAction Code, and take final position of Z which absorbs in Si target.

I made 1000 runs , and in each run, I made 10000events.
So total 10,000,000 events.
About 20% goes outOfWorld so I rejected.

For each 1000 run, I calculated 95%containment depth using 10000events, and I calculated the average.

The result is like below.

I really wonder why the results are different.
Actually I tested another physicslists for low EM interaction, but the result was same.
I wish my method of making the code is wrong.

Thank you in advance for your help.

Best regards,
Hyeonseok

jrellin · February 25, 2025, 6:30am

It is far easier to do this with PreUserTrackingAction and PostUserTrackingAction since both of these are called exactly once. The first at “creation” and the second when “killed”. It also helps you avoid odd exceptions like you do here since particles that escape will be at the edge of the world volume anyway:

 G4ThreeVector finalpostPosition = (postVolume && postVolume->GetName() == "Silicon") ? postPosition : G4ThreeVector(-999,-999,-9999);

Lastly, two questions. How are you defining your Silicon block? Second, how are you defining depth. The distances you care about are in cm but everything is in meters…

About 20% goes outOfWorld so I rejected.

How is that possible if the range is 10 cm in a 50 m block?

HyeonSeok-Do · February 25, 2025, 7:33am

Thank you for your reply.

I also used TrackingAction Class before, but the result was the same.

First, my detectorConstruction is like this.
‘’’
G4VPhysicalVolume* DetectorConstruction::Construct()
{
ConstructMaterials();
G4Material* vacuum = G4Material::GetMaterial(“G4_Galactic”);
G4Material* silicon = G4Material::GetMaterial(“G4_Si”);
G4bool checkOverlaps = true;

//==============================
// World volume
//==============================
G4double worldSizeX = 100 * m;
G4double worldSizeY = 100 * m;
G4double worldSizeZ = 100 * m;
G4VSolid* worldSolid = new G4Box("world", worldSizeX/2, worldSizeY/2, worldSizeZ/2);
G4LogicalVolume* worldLogical = new G4LogicalVolume(worldSolid, vacuum, "worldLogical");
G4VPhysicalVolume* worldPhysical = new G4PVPlacement(
            nullptr,
            G4ThreeVector(),
            worldLogical,
            "World",
            nullptr,
            false,
            0,
            checkOverlaps
        );

//==============================
// Xray Scanner (silicon detector)
//==============================
G4double SizeX = 50 * m;
G4double SizeY = 50 * m;
G4double SizeZ = 50 * m;
// Pixel 크기 (각 Pixel의 크기)
G4Box* siliconSolid = new G4Box("silicon", SizeX / 2, SizeY / 2, SizeZ / 2);
G4LogicalVolume* siliconLV = new G4LogicalVolume(siliconSolid, silicon, "siliconLogical");

new G4PVPlacement(
        nullptr,
        G4ThreeVector(0, 0, 0),
        siliconLV,
        "Silicon",
        worldLogical,
        false,
        0,
        checkOverlaps
    );

‘’’
When I visualize this(simply for 200 events),

The second image is zooming in on the interaction point.

The reason for using such a huge target is that I wanted to grab whole events that live long.
About a total of 10,000,000 events, the longest event(by Z direction) goes about 50cm.

I think the events that go outOfWorld are backscattered by Compton scattering.
So I rejected these events and only used the events which finally absorbed in silicon by the photoelectric effect.

I just sorted the z positions of events absorbed by the silicon in ascending order and then truncated the list at the 95% mark.
Is my method of getting 95% containment depth wrong?
Do I need to consider not only the Z direction but also the X, and Y directions?

I appreciate your help.

pico · February 25, 2025, 10:36am

Hello,

The exponential attenuation is only an approximation for a collimated monoenergetic beam. From the picture I am not convinced the beam is collimated (or at least collinear).

My two cents.

/Pico

jrellin · February 25, 2025, 2:56pm

Shrink the dimensions of everything to speed up computation. Also that 20% that you throw out points to a flaw in how you are assessing range. You should be able to capture these values since if you took a weighted average of near 0 with the 80% you are keeping you get much closer to the “hand calculation”.

One trick would be to keep the stepping action that you have and just record the Z everywhere. Initialize some “zmax” variable at -9999. With every step check if the Z is greater than zmax, than make it the new zmax if so. Handles all the edges cases and is independent of volume/boundary issues.

maire1 · February 25, 2025, 5:38pm

1- Geant4 data set are printed in example TestEm0
Run the attached macro (em0.mac.txt) and compare results with NIST/xcom values
em0.mac.txt (236 Bytes)

2- Absorption - tranmission ratios are simulated in example TestEm13
Take the attached macro (em13.mac.txt) and play with size to find the 5% transmission
PS. in PhysListEmStandard::ConstructProcess(),do not forget to uncomment G4RayleighScattering process
em13.mac.txt (343 Bytes)

HyeonSeok-Do · February 27, 2025, 7:48am

Dear Michel,

I appreciate your help.
It was the thing that I wanted to know.
Maybe there was some misunderstanding about the exponential attenuation formula.
I think my simulation method was wrong.

Plus, thank you to all the other experts who helped me.

Best regards,
Hyeonseok

Best

maire1 · February 27, 2025, 10:37am

To test the exponential law, you must measure the track length of the primary particle from its creation to its first interaction only.
Compton or Rayleigh processes may have several successive interactions; they must not be taken into account.
In TestEm13, we kill an event at end of the first step of the primary particle …

HyeonSeok-Do · February 27, 2025, 11:10am

Dear Michel,

Thank you for your comment.
I studied the experimental setup for the exponential formula thanks to you.

My understanding is that this formula is focused on “Transmittance”, not absorption.
I mean, is the expression ‘5% transmitted depth’ more appropriate than ‘95% absorbed depth’?

Lastly, does TestEm13 consider an event where a photon interacts even once with the material as being absorbed, even if the photon could eventually pass through the material?

Best regards,
Hyeonseok.

maire1 · February 27, 2025, 4:18pm

The answer to your last question is yes.
In TesTem13 printout, we prefer the expression unaltered . A particle which has interacted - whatever happen after this interaction - is no more unaltered.
It was also the meaning of Pico answer : collimated monoenergertic beam

maire1 · March 6, 2025, 8:50am

1- TestEm13 is an exercise to count the number of incident monoenergetic particles which have not yet interacted after a given distance. This number is governed by an exponential law.

2- it is a different exercise to count the number of incident particles which remain confined within a given volume. This number is not directly governed by an exponential law.

Your computation by hand correspond to the first situation, while your simulation is the second situation.

maire1 · March 7, 2025, 2:11pm

In my previous post, I did not answered to one of your question.

The purpose of TestEm13 is to check if a primary particle interact or not. Once we get this information, the rest of the event has no interest. We drop it only to save time ! there is no other physics reason behind …

If you prefer, you can keep the event, but with some additional protection in SteppingAction. Something like :
if ((itrack==1)&&(istep==1)) run->countProcesses(procName)

(Indeed, I did not check this …)

HyeonSeok-Do · March 8, 2025, 5:11am

I appreciate your kind reply.
Thanks to you, I can understand the real meaning of the formula and check for misunderstandings.
Thank you again.

system · March 15, 2025, 5:11am

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