Hi Experts,
I am attempting to use Pol01 example to examine the angular distribution of initially polarized photons Compton scattering off G4_Si.
I have had to include several lines of code in the SteppingAction.cc which will print out information like, for example, the Compton scattering angle and the photon’s final kinetic energy.
The simulations seems to work well for horizontally and vertically polarized photons, linear polarized ± 45 degrees and right circularly polarized (RCP) photons. But when I try and run the simulation for left circuarly polarized (LCP) photons I get a segmentation fault part way through the printing of the information to my terminal window.
In macro code this is what I found
### RCP (works)
# /gun/polarization 0. 0. 1.
### LCP ( Segmentation fault )
/gun/polarization 0. 0. -1.
### vertical ( works )
# /gun/polarization -1. 0. 0.
### horizontal ( works )
# /gun/polarization 1. 0. 0.
### LP ( 45 degrees ) (works)
# /gun/polarization 0. 1. 0.
### LP ( -45 degrees ) ( works )
# /gun/polarization 0. -1. 0.
/gun/particle gamma
#
/gun/energy 511 keV
However, above this code there is the following instruction
/polarization/volume/set theBox 0. 0. 0.08
If I set 0.08 to 0. , then I do not get the segmentation fault. I would like help in figuring out why only LCP photons cause the segmentation fault?
FYI, the only extra code I included into Pol01 is as follows
// Inside this method you should access the track from the step object:
G4Track *aTrack = aStep->GetTrack();
// Then you should get the KineticEnergy of the track if the process is
// "Compton scattering" and if the volume remains the same in the next step.
// See the code lines below:
G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
const G4VProcess* CurrentProcess = preStepPoint->GetProcessDefinedStep();
if (CurrentProcess != 0)
{
// To implement the next line, one needs to include " #include "G4VProcess.hh" "
// (See SteppingAction.cc in TestEm6)
const G4String & StepProcessName = CurrentProcess->GetProcessName();
// Get track ID for current processes
G4double TrackID = aStep->GetTrack()->GetTrackID();
if( (StepProcessName == "pol-compt")&&(TrackID==1) )
{
// processing hit when entering the volume
G4double kineticEnergy = aStep->GetTrack()->GetKineticEnergy();
auto EkinUnits = G4BestUnit(kineticEnergy, "Energy");
// Get energy in units of keV
G4double En = kineticEnergy/keV;
// Determine angle in radians
G4double thetaRad = std::acos( 2- 511/En );
// Place condition to reject Compton scattering at 0 rad.
// (This is to avoid nan results when dividing through by 2pisin(theta))
if ( (thetaRad > 0) && (thetaRad < pi))
{
// fill histogram id = 13 (angular spread (in radians) as a function of scattering angle)
// analysisManager->FillH1(13, thetaRad);
// get particle name that's undergone the Compton scattering
G4String ParticleName = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
// Find position (c.f. /examples/extended/medical/GammaTherapy) in CheckVolumeSD.cc
// G4ThreeVector pos = aStep->GetPreStepPoint()->GetPosition();
// G4double z = pos.z();
// Check if volume is sensitive detector
G4String fromVolume = aTrack->GetNextVolume()->GetName();
// see https://geant4-forum.web.cern.ch/t/get-particle-info-in-steppingaction-if-particle-pass-through-specific-volume/2863/2
// for possible solution for finding which volume the gamma is in
// G4String SDname = aStep->GetPostStepPoint()->GetSensitiveDetector()->GetName();
// FOR CHECKING
std::cout << "fromVolume: " << fromVolume
<< " Particle: " << ParticleName
<< " Process: " << StepProcessName
<< " Energy: " << EkinUnits
<< " Angle: " << thetaRad/degree << " deg." << std::endl;
}
}
}
Thank you for your time.
Peter