Proton grazing angle simulation

Physics Processes, Models and Cross Sections
1

_Geant4 Version:_11.3.1
_Operating System:_Debian 6.1.133-1
_Compiler/Version:_CMake
_CMake Version:_3.25.1

I am trying to recreate a paper where a proton beam hits a gold mirror and then elasticly gets scattered to a detector. To simulate it, I shoot protons on a gold block from shallow angles (<5 deg).
The problem that I am running into is that the under a certain proton energy, it does not interact with the gold anymore. I have added the output for proton energies 1300-1400 keV with DumpInfo() about the used models and then the output given by verbose 2.
1400kev.txt (15.7 KB)
1300kev.txt (13.7 KB)

In the 1400 keV the proton does interact with the gold, but in the 1300 keV it does not, and I do not really understand why not. It feels like the models used have a cutoff at this energy range, but I cannot seem to find the documentation on it. I have recompiled G4 with the G4TENDL model and have exported the location into a new .bashrc.
I have also added MaxStepSizes to both the logic of the gold and that of the vacuum to make sure the proton does not pass through the gold block without interacting.

I am using the following in my physics list:
G4EmLivermorePhysics
G4OpticalPhysics
G4DecayPhysics
G4RadioactiveDecayPhysics
G4HadronPhysicsQGSP_BERT
G4HadronElasticPhysicsHP
G4IonElasticPhysics
G4IonPhysicsPHP
G4HadronPhysicsQGSP_BERT_HP

Please let me know if any code snippets would help in figuring this out!

2

Here are code snippets from the construction.cc of the gold block.

G4VPhysicalVolume* MyDetectorConstruction::Construct() {
G4NistManager* nist = G4NistManager::Instance();

// Define materials
G4Material* gold = nist->FindOrBuildMaterial("G4_Au");
G4Material* worldMat = nist->FindOrBuildMaterial("G4_Galactic"); 
G4Material* silicon = nist->FindOrBuildMaterial("G4_Si");

// Define world volume
G4double xWorld = 1 * m;
G4double yWorld = 1 * m;
G4double zWorld = 1 * m;

G4Box* solidWorld = new G4Box("solidWorld", xWorld, yWorld, zWorld);
G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, worldMat, "logicWorld");
G4VPhysicalVolume* physWorld = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), logicWorld, "physWorld", 0, false, 0, true);

// Add a step limiter to the vacuum (world volume)
G4double maxStepVacuum = 0.1 * mm; // Set a maximum step size for the vacuum
logicWorld->SetUserLimits(new G4UserLimits(maxStepVacuum));

// Place the gold block
G4Box* solidGoldBlock = new G4Box("solidGoldBlock", goldBlockWidth / 2, goldBlockThickness / 2, goldBlockHeight / 2);
G4LogicalVolume* logicGoldBlock = new G4LogicalVolume(solidGoldBlock, gold, "logicGoldBlock");
new G4PVPlacement(0, G4ThreeVector(0., -goldBlockThickness / 2, 0.), logicGoldBlock, "physGoldBlock", logicWorld, false, 0, true);

// Add a step limiter to the gold block
G4double maxStepGold = 0.00001 * mm; // Set a maximum step size for the gold block
logicGoldBlock->SetUserLimits(new G4UserLimits(maxStepGold));

// Place the detector
G4Tubs* siliconDetector = new G4Tubs("siliconDetector", 0, detectorRadius, detectorThickness / 2, 0, 2 * M_PI);
logicDetector = new G4LogicalVolume(siliconDetector, silicon, "logicDetector");
new G4PVPlacement(0, G4ThreeVector(0., 0., -detectorZPosition), logicDetector, "physDetector", logicWorld, false, 0, true);

return physWorld;
}

And here is part of the generator.cc

MyPrimaryGenerator::MyPrimaryGenerator()
{
fParticleGun = new G4ParticleGun(1); // Initialize particle gun with 1 particle per event

G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(0.1 * eV, 100 * GeV);

G4double lowEnergyCut = 0.1 * keV;
G4ProductionCuts* cuts = new G4ProductionCuts();
cuts->SetProductionCut(lowEnergyCut, "gamma");
cuts->SetProductionCut(lowEnergyCut, "e-");
cuts->SetProductionCut(lowEnergyCut, "e+");
cuts->SetProductionCut(lowEnergyCut, "proton");

G4Region* region = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
region->SetProductionCuts(cuts);
}

MyPrimaryGenerator::~MyPrimaryGenerator()
{
delete fParticleGun;
}

void MyPrimaryGenerator::GeneratePrimaries(G4Event* anEvent)
{
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

// Generate a proton
G4String particleName = "proton";
G4ParticleDefinition* particle = particleTable->FindParticle(particleName);

// Calculate the starting position based on the event ID
G4int eventID = anEvent->GetEventID();
G4double angle = eventID * (5.0 / 100.0) * deg; // Increment angle from 0 to 5 degrees over 100 events
if (angle > 5.0 * deg) {
angle = 5.0 * deg; 
}

// Set the particle's starting position
G4double yStart = 0.001 * m + std::tan(angle) * 1.0 * m; // Offset in y-axis based on angle
G4ThreeVector pos(0.0 * m, yStart, 0.1 * m); // start position
fParticleGun->SetParticlePosition(pos);

// Set the momentum direction to point toward the origin
G4ThreeVector mom(0.0, -yStart, -0.1 * m); // Momentum direction toward the origin
mom = mom.unit(); // Normalize the momentum vector
fParticleGun->SetParticleMomentumDirection(mom);

// Set the particle energy
fParticleGun->SetParticleMomentum(1300 * keV); 

// Set the particle definition
fParticleGun->SetParticleDefinition(particle);

// Generate the primary vertex
fParticleGun->GeneratePrimaryVertex(anEvent);
}

I am trying to recreate this paper.

3

Hi @Friso_Eigenhuis

Thank you for your post.

I believe the electromagnetic extended example TestEm5 (link) may help understand everyone the issue. Its geometry is quite similar to yours: a block of material placed in a world volume.

Could you please share a macro file that reproduces the issue using TestEm5?

Also, could you provide the bibliographic reference for the paper you mentioned?

Thank you for your time.

Best regards,
Alvaro

4

Thank you for your reply!

I have been trying to recreate the problem in the macro file of the TestEm5 example but without any success. If I adjust the hanson.mac to use low energy protons and the protons seem to interact with the gold as expected.

I have tried implementing the same physics list items in my own simulation it still does not work. These are G4EmStandardPhysics_option3, G4hPairProduction, G4hBremsstrahlung, G4hIonisation, G4hMultipleScattering. I have also tried to set /run/setCut 0.1 um the same as TestEm5 but without any difference.

My apologies about the wrong link, here is the paper I am trying to recreate as stepping stone to simulating a X-ray wolter type 1 nested mirror:
Diebold, S., Tenzer, C., Perinati, E. et al. Soft proton scattering efficiency measurements on x-ray mirror shells. Exp Astron 39, 343–365 (2015). Soft proton scattering efficiency measurements on x-ray mirror shells | Experimental Astronomy

So the issue still persists that under 1400 keV the protons just disappears at the intersection without any physics process. What are some steps I could take to solve this issue?

5

I have found the problem, at these low energies it becomes really important to use energy instead of momentum in the particle gun.
So instead of:

fParticleGun->SetParticleMomentum(1300 * keV);

it should be:

fParticleGun->SetParticleEnergy(1300* keV);

Bit of a oversight on my part.

6

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