Why is sometimes Gaussian used in G4UniversalFluctuation energy straggling?

Physics Processes, Models and Cross Sections
1

Hi,

I am relatively new to Geant4. I was going through the source code to understand it better and I saw in G4UniversalFluctuation for calculating the energy straggling a combination of Gaussian and Gamma distribution is used depending on the mean to standard deviation ratio.

So my quesiton is, why isn’t the Gamma distribution used all the time? Is there big performance cost in using the Gamma distribution over the Gaussian distribution? Or is some other reason I’m not seeing?

Thanks,
Declan

The code snippet in question is:

G4double 
G4UniversalFluctuation::SampleFluctuations(const G4MaterialCutsCouple* couple,
                                           const G4DynamicParticle* dp,
                                           const G4double tcut,
                                           const G4double tmax,
                                           const G4double length,
                                           const G4double averageLoss)
{
  // Calculate actual loss from the mean loss.
  // The model used to get the fluctuations is essentially the same
  // as in Glandz in Geant3 (Cern program library W5013, phys332).
  // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual

  // shortcut for very small loss or from a step nearly equal to the range
  // (out of validity of the model)
  //
  if (averageLoss < minLoss) { return averageLoss; }
  meanLoss = averageLoss;
  const G4double tkin  = dp->GetKineticEnergy();
  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;

  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }

  CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
             
  const G4double gam   = tkin * m_Inv_particleMass + 1.0;
  const G4double gam2  = gam*gam;
  const G4double beta  = dp->GetBeta(); 
  const G4double beta2 = beta*beta;

  G4double loss(0.), siga(0.);

  const G4Material* material = couple->GetMaterial();
  
  // Gaussian regime
  // for heavy particles only and conditions
  // for Gauusian fluct. has been changed 
  //
  if (particleMass > CLHEP::electron_mass_c2 &&
      meanLoss >= minNumberInteractionsBohr*tcut && tmax <= 2.*tcut) {

    siga = std::sqrt((tmax/beta2 - 0.5*tcut)*CLHEP::twopi_mc2_rcl2* 
                      length*chargeSquare*material->GetElectronDensity());
    const G4double sn = meanLoss/siga;
  
    // thick target case 
    if (sn >= 2.0) {

      const G4double twomeanLoss = meanLoss + meanLoss;
      do {
	loss = G4RandGauss::shoot(rndmEngineF, meanLoss, siga);
	// Loop checking, 03-Aug-2015, Vladimir Ivanchenko
      } while  (0.0 > loss || twomeanLoss < loss);

      // Gamma distribution
    } else {

      const G4double neff = sn*sn;
      loss = meanLoss*G4RandGamma::shoot(rndmEngineF, neff, 1.0)/neff;
    }
    //G4cout << "Gauss: " << loss << G4endl;
    return loss;
  }
....

Geant4 Version: 11.2.0

2

My guess is that you are right about the computational speed (see this stack overflow page and comments) that is generally true. Gaussians have closed forms which make sampling from them significantly “easier”.

3

Ok that make sense. Thank you!