Issues in implementing a custom EM process using G4VDiscreteProcess()

Physics Processes, Models and Cross Sections
1

Hello all,

I aim to implement a custom gamma process in a “Crystal” volume by replacing the Rayleigh scattering process (using /process/inactivate Rayl) while still keeping other low-energy processes (i.e., compt and phot). I am using a class derived from G4VDiscreteProcess for this. The gamma energies used in my code are limited from 4 keV to 200 keV (I am using G4EmPenelopePhysics). The issue: I see more singly scattered Compton events in the volume even if the mean free path proposed by my process is much smaller than that of the Compton process for the energy range used. I see my process dominating only for very low energies (e.g., less than 10 keV).

This plot shows an example of the attenuation coefficients taking part in my code’s physics:

example_mu
example_mu1920×940 45.7 KB

To debug this issue, I also tried to check for the mean free path of compt within the PostStepGetPhysicalInteractionLength() and always returned some much smaller value than it in the following way, but still see the same effect:

G4double myProcess::PostStepGetPhysicalInteractionLength(const G4Track& track,
                                                               G4double,
                                                               G4ForceCondition* condition)
{
    *condition = NotForced; 
    G4StepPoint* postStep = track.GetStep()->GetPostStepPoint();
    G4String postStepVolume = postStep->GetPhysicalVolume()->GetName();

    if(postStepVolume != "Crystal")
        return DBL_MAX;

    G4EmCalculator emCalculator;
    G4double energy = postStep->GetKineticEnergy();
    G4String material = track.GetMaterial()->GetName();
    G4double mfp_compt = emCalculator.GetMeanFreePath(energy, "gamma", "compt", material, "Crystal");

    G4double mfp_myProcess = mfp_compt/10;
    return mfp_myProcess;
}

in my header myProcess.hh, I also have this section:

class myProcess : public G4VDiscreteProcess
{
// ...
protected:
    virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition*)
    {return 0;}
// ...
}

What am I missing here?
Thank you for your time and expertise!

UPDATE: Looking at the Geant4 User’s Guide for Application Developers, there is a mention that:

G4double PostStepGetPhysicalInteractionLength( const G4Track& track,
G4double previousStepSize, G4ForceCondition* condition )

This method generates the step length allowed by its process. It also provides a flag to force the interaction to occur regardless of its step length.

So, instead of providing the mean free path, should I directly provide the step length, e.g., using: return CLHEP::RandExponential::shoot(mfp_myProcess);?