I am trying to calculate the energy/momentum distribution at the last step in the geometry boundary, just before particle exits the boundary.
I am using the following approach for that,
Firstly, I define a geometry in the Detector Construction, with a world volume and my detector geometry contained within (with Physical and Logical volume defined).
In the UserSteppingAction, I define the following, to get to the last step in volume.
const G4StepPoint* endPoint = step->GetPostStepPoint();
G4ParticleDefinition* prim_particle = fPrimary->GetParticleGun()->GetParticleDefinition();
G4bool condition = step->IsLastStepInVolume();
if ((prim_particle == G4MuonMinus::MuonMinus()) && condition == true)
{
G4cout << "Step is the last step in the volume" << G4endl;
energy_mu = endPoint->GetKineticEnergy();
pX_mu = endPoint->GetMomentum().x();
posX_mu = endPoint->GetPosition().x();
.....
}
The condition always comes out to be false (0 on printing it out), even though there are numerous primary particles being transported out of the volume. (I get the end of step energy and momentum values properly if I remove the condition).
Another approach I tried was using the ‘fGeomBoundary’ as follows:
In the context of TestEm3, I have checked your two tests, whit a geantino.
I attach the macro (can be run in batch or interactively with visualisation), the 2 outputs and the code added in SteppingAction. aditya.mac.txt (374 Bytes) aditya.lastStep.txt (2.3 KB) aditya.fGeomBoundary.txt (2.4 KB) aditya.stepAction.txt (448 Bytes)
I am still unable to debug the issue. My geometry is a simple rock geometry with particles (muons) being propagated via General Particle Source (GPS).
After printing out the tracking information, I get the following for one of the events…
In the stepping action, I use the command the same way as yours.
if (step->IsLastStepInVolume()) {
G4cout << "Step is the last step in the volume " << volume->GetName()<< "\n" << G4endl;
energy_mu = endPoint->GetKineticEnergy();
G4cout << "Muon energy on exit " << energy_mu << G4endl;
}
I try to get the output to the function call for IsLastStepInVolume, but it always comes out to be ‘0’
Similarily, fGeomBoundary doesn’t seem to work either. Maybe I’m missing something, I can’t quite figure out.
The particles are clearly transporting out of the volume as in the attached document. I don’t quite understand what I’m missing.
1- Did you try with a geantino ?
2- Are you sure that the volumes do not overlap ?
ps. May I suggest you to try to implement the SteppingVerbose of TestEm examples (but it will not solve the problem !).
Yes, I did try with a geantino, it passes straight through the geometry without taking any steps in it.
*********************************************************************************************************
* G4Track Information: Particle = geantino, Track ID = 1, Parent ID = 0
*********************************************************************************************************
Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName
0 -1.7e+05 -7.3e+04 5e+04 6.8e+03 0 0 0 RockBurden initStep
1 -1.59e+05 9.34e+03 -5e+04 6.8e+03 0 1.3e+05 1.3e+05 World Transportation
2 -5.32e+04 7.92e+05 -1e+06 6.8e+03 0 1.24e+06 1.37e+06 OutOfWorld Transportation
Rotated and Translated position (-44228.88199316448,-29116.20538633069,50000)
Final User Defined momentum vector (0,0,-1)
I am using the checkOverlaps boolean in the PVplacement which is fine.
Checking overlaps for volume RockBurden (G4Box) ... OK!
Checking overlaps for volume Scorer (G4Box) ... OK!
I did implement SteppingVerbose to see if it gives any further clues.
Do you think I could use this,
fUserDefinedLimit // Step defined by the user Step limit in the logical volume,
…to set a limit on the step coordinate, which when it reaches would give the energy/momenta at that step point?
I have attached my SteppingAction file.
Yes, none of the particles see the IsLastStepInVolume() condition, although without it, I can save the energy/momentum values to the Ntuple, but for all steps.
Also, for my purpose, I’m using General Particle Source (GPS) using position, energy and angular distribution to shoot particles through the geometry. Should that make any difference? (I’ll also attach the GPS macro)