#include <G4UniversalFluctuation.hh>

Inheritance diagram for G4UniversalFluctuation:

Public Member Functions
	G4UniversalFluctuation (const G4String &nam="UniFluc")

	~G4UniversalFluctuation () override

G4double	SampleFluctuations (const G4MaterialCutsCouple , const G4DynamicParticle , const G4double, const G4double, const G4double, const G4double) override

G4double	Dispersion (const G4Material , const G4DynamicParticle , const G4double, const G4double, const G4double) override

void	InitialiseMe (const G4ParticleDefinition *) override

void	SetParticleAndCharge (const G4ParticleDefinition *, G4double q2) override

G4UniversalFluctuation &	operator= (const G4UniversalFluctuation &right)=delete

	G4UniversalFluctuation (const G4UniversalFluctuation &)=delete

Public Member Functions inherited from G4VEmFluctuationModel
	G4VEmFluctuationModel (const G4String &nam)

virtual	~G4VEmFluctuationModel ()

virtual G4double	SampleFluctuations (const G4MaterialCutsCouple , const G4DynamicParticle , const G4double tcut, const G4double tmax, const G4double length, const G4double meanLoss)=0

virtual G4double	Dispersion (const G4Material , const G4DynamicParticle , const G4double tcut, const G4double tmax, const G4double length)=0

virtual void	InitialiseMe (const G4ParticleDefinition *)

virtual void	SetParticleAndCharge (const G4ParticleDefinition *, G4double q2)

const G4String &	GetName () const

G4VEmFluctuationModel &	operator= (const G4VEmFluctuationModel &right)=delete

	G4VEmFluctuationModel (const G4VEmFluctuationModel &)=delete

Protected Member Functions
virtual G4double	SampleGlandz (CLHEP::HepRandomEngine rndm, const G4Material , const G4double tcut)

void	AddExcitation (CLHEP::HepRandomEngine *rndm, const G4double ax, const G4double ex, G4double &eav, G4double &eloss, G4double &esig2)

void	SampleGauss (CLHEP::HepRandomEngine *rndm, const G4double eav, const G4double esig2, G4double &eloss)

Protected Attributes
G4double	particleMass = 0.0

G4double	m_Inv_particleMass = DBL_MAX

G4double	m_massrate = DBL_MAX

G4double	chargeSquare = 1.0

G4double	ipotFluct = 0.0

G4double	ipotLogFluct = 0.0

G4double	e0 = 0.0

G4double	minNumberInteractionsBohr = 10.0

G4double	minLoss

G4double	nmaxCont = 8.0

G4double	rate = 0.56

G4double	fw = 4.0

G4double	a0 = 42.0

G4double	w2 = 0.0

G4double	meanLoss = 0.0

const G4ParticleDefinition *	particle = nullptr

G4double *	rndmarray = nullptr

G4int	sizearray = 30

Detailed Description

Definition at line 55 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

◆ G4UniversalFluctuation() [1/2]

G4UniversalFluctuation::G4UniversalFluctuation ( const G4String & nam = "UniFluc" )

explicit

Definition at line 58 of file G4UniversalFluctuation.cc.

 :G4VEmFluctuationModel(nam),
  minLoss(10.*CLHEP::eV)
{
  rndmarray = new G4double[sizearray];
}

◆ ~G4UniversalFluctuation()

G4UniversalFluctuation::~G4UniversalFluctuation ( )

override

Definition at line 67 of file G4UniversalFluctuation.cc.

{
  delete [] rndmarray;
}

◆ G4UniversalFluctuation() [2/2]

G4UniversalFluctuation::G4UniversalFluctuation ( const G4UniversalFluctuation & )

delete

Member Function Documentation

◆ AddExcitation()

void G4UniversalFluctuation::AddExcitation	(	CLHEP::HepRandomEngine *	rndm,
		const G4double	ax,
		const G4double	ex,
		G4double &	eav,
		G4double &	eloss,
		G4double &	esig2
	)

inlineprotected

Definition at line 129 of file G4UniversalFluctuation.hh.

{
  if(ax > nmaxCont) {
    eav  += ax*ex;
    esig2 += ax*ex*ex;
  } else {
    const G4int p = (G4int)G4Poisson(ax);
    if(p > 0) { eloss += ((p + 1) - 2.*rndm->flat())*ex; }
  }
}

Referenced by SampleGlandz(), and G4UrbanFluctuation::SampleGlandz().

◆ Dispersion()

G4double G4UniversalFluctuation::Dispersion	(	const G4Material *	material,
		const G4DynamicParticle *	dp,
		const G4double	tcut,
		const G4double	tmax,
		const G4double	length
	)

overridevirtual

Implements G4VEmFluctuationModel.

Definition at line 244 of file G4UniversalFluctuation.cc.

{
  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
  const G4double beta = dp->GetBeta();
  return (tmax/(beta*beta) - 0.5*tcut) * CLHEP::twopi_mc2_rcl2 * length
    * material->GetElectronDensity() * chargeSquare;
}

◆ InitialiseMe()

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition * part )

overridevirtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 74 of file G4UniversalFluctuation.cc.

{
  particle = part;
  particleMass = part->GetPDGMass();
  const G4double q = part->GetPDGCharge()/CLHEP::eplus;
 
  // Derived quantities
  m_Inv_particleMass = 1.0 / particleMass;
  m_massrate = CLHEP::electron_mass_c2 * m_Inv_particleMass;
  chargeSquare = q*q;
}

Referenced by Dispersion(), G4IonFluctuations::InitialiseMe(), and SampleFluctuations().

◆ operator=()

G4UniversalFluctuation & G4UniversalFluctuation::operator= ( const G4UniversalFluctuation & right )

delete

◆ SampleFluctuations()

G4double G4UniversalFluctuation::SampleFluctuations	(	const G4MaterialCutsCouple *	couple,
		const G4DynamicParticle *	dp,
		const G4double	tcut,
		const G4double	tmax,
		const G4double	length,
		const G4double	averageLoss
	)

overridevirtual

Implements G4VEmFluctuationModel.

Definition at line 89 of file G4UniversalFluctuation.cc.

{
  // 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;
  }
 
  auto ioni = material->GetIonisation();
  e0 = ioni->GetEnergy0fluct();
 
  // very small step or low-density material
  if(tcut <= e0) { return meanLoss; }
 
  ipotFluct = ioni->GetMeanExcitationEnergy();
  ipotLogFluct = ioni->GetLogMeanExcEnergy();
 
  // width correction for small cuts
  const G4double scaling = std::min(1.+0.5*CLHEP::keV/tcut, 1.50);
  meanLoss /= scaling;
 
  w2 = (tcut > ipotFluct) ? 
    G4Log(2.*CLHEP::electron_mass_c2*beta2*gam2)-beta2 : 0.0;
  return SampleGlandz(rndmEngineF, material, tcut)*scaling;
}

Referenced by G4IonFluctuations::SampleFluctuations().

◆ SampleGauss()

void G4UniversalFluctuation::SampleGauss	(	CLHEP::HepRandomEngine *	rndm,
		const G4double	eav,
		const G4double	esig2,
		G4double &	eloss
	)

inlineprotected

Definition at line 146 of file G4UniversalFluctuation.hh.

{
  G4double x = eav;
  const G4double sig = std::sqrt(esig2);
  if(eav < 0.25*sig) {
    x += (2.*rndm->flat() - 1.)*eav;
  } else {
    do { 
      x = G4RandGauss::shoot(rndm, eav, sig);
    } while (x < 0.0 || x > 2*eav);
    // Loop checking, 23-Feb-2016, Vladimir Ivanchenko
  }
  eloss += x;
} 

Referenced by SampleGlandz(), and G4UrbanFluctuation::SampleGlandz().

◆ SampleGlandz()

G4double G4UniversalFluctuation::SampleGlandz	(	CLHEP::HepRandomEngine *	rndm,
		const G4Material *	,
		const G4double	tcut
	)

protectedvirtual

Reimplemented in G4UrbanFluctuation.

Definition at line 173 of file G4UniversalFluctuation.cc.

{
  G4double a1(0.0), a3(0.0);
  G4double loss = 0.0;
  G4double e1 = ipotFluct;
 
  if(tcut > e1) {
    a1 = meanLoss*(1.-rate)/e1;
    if(a1 < a0) {
      const G4double fwnow = 0.1+(fw-0.1)*std::sqrt(a1/a0);
      a1 /= fwnow;
      e1 *= fwnow;
    } else {
      a1 /= fw;
      e1 *= fw;
    }   
  }
 
  const G4double w1 = tcut/e0;
  a3 = rate*meanLoss*(tcut - e0)/(e0*tcut*G4Log(w1));
  if(a1 <= 0.) { a3 /= rate; }
  
  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont  
  G4double emean = 0.;
  G4double sig2e = 0.;
 
  // excitation of type 1
  if(a1 > 0.0) { AddExcitation(rndmEngineF, a1, e1, emean, loss, sig2e); }
 
  if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
 
  // ionisation 
  if(a3 > 0.) {
    emean = 0.;
    sig2e = 0.;
    G4double p3 = a3;
    G4double alfa = 1.;
    if(a3 > nmaxCont) {
      alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
      const G4double alfa1  = alfa*G4Log(alfa)/(alfa-1.);
      const G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
      emean += namean*e0*alfa1;
      sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
      p3 = a3 - namean;
    }
 
    const G4double w3 = alfa*e0;
    if(tcut > w3) {
      const G4double w = (tcut-w3)/tcut;
      const G4int nnb = (G4int)G4Poisson(p3);
      if(nnb > 0) {
        if(nnb > sizearray) {
          sizearray = nnb;
          delete [] rndmarray;
          rndmarray = new G4double[nnb];
        }
        rndmEngineF->flatArray(nnb, rndmarray);
        for (G4int k=0; k<nnb; ++k) { loss += w3/(1.-w*rndmarray[k]); }
      }
    }
    if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
  }
  //G4cout << "### loss=" << loss << G4endl;
  return loss;
}

Referenced by SampleFluctuations().

◆ SetParticleAndCharge()

void G4UniversalFluctuation::SetParticleAndCharge	(	const G4ParticleDefinition *	part,
		G4double	q2
	)

overridevirtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 260 of file G4UniversalFluctuation.cc.

{
  if(part != particle) {
    particle = part;
    particleMass = part->GetPDGMass();
 
    // Derived quantities
    m_Inv_particleMass = 1.0 / particleMass;
    m_massrate = CLHEP::electron_mass_c2 * m_Inv_particleMass;
  }
  chargeSquare = q2;
}