I have solved the problem, with the help of this excellent explanation.
https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4
The histogram below shows the agreement I get when I compare the angular distribution compared to Klein-Nishina theory. The title should read 0.5 mm cube of Si.
(The adjustment to my ProccessHits code + some explanation follows. I hope this can help others.)
G4bool SensitiveDetectorCPET04::ProcessHits(G4Step* aStep, G4TouchableHistory*)
{
// Instantiate G4Analysismanager
auto analysisManager = G4AnalysisManager::Instance();
// c.f. https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4
// 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 == "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 = 0 (angular spread (in radians) as a function of scattering angle)
analysisManager->FillH1(0, thetaRad); // (This produces the data in the above graph) ---> I implement the theoretical curve in my root script.
// fill of histogram id = 1 (angular spread (in radians) as a function of soild angle)
// c.f. TestEm5 (TrackingAction.cc) and also (corrections at)
// https://geant4-forum.web.cern.ch/t/physics-validation/2893
G4double dthetaRad1 = ( analysisManager->GetH1Width(1) );
G4double unit1 = analysisManager->GetH1Unit(1);
G4double weight1 = (unit1)/( twopi*std::sin(thetaRad)*dthetaRad1 );
analysisManager->FillH1(1, thetaRad, weight1); **GRAPH NOT SHOWM**
// 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 x = pos.x();
// G4double y = pos.y();
G4double z = pos.z();
// Fill NTuples
analysisManager->FillNtupleDColumn(0, En/keV);
analysisManager->FillNtupleDColumn(1, thetaRad/deg);
// analysisManager->FillNtupleDColumn(2, x/um);
// analysisManager->FillNtupleDColumn(3, y/um);
analysisManager->FillNtupleDColumn(2, z/um);
analysisManager->AddNtupleRow(0);
// FOR CHECKING
std::cout << "fromVolume: " << fromVolume
// << " unit: " << unit
<< " Particle: " << ParticleName
<< " Process: " << StepProcessName
<< " Energy: " << EkinUnits
<< " Angle: " << thetaRad/degree << " deg." << std::endl;
// << " Scattering site: " << pos << std::endl;
}
}
}
return true;
}
Best
Peter