#include <G4UniversalFluctuation.hh>

Inheritance diagram for G4UniversalFluctuation:

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

virtual	~G4UniversalFluctuation ()

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

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

virtual void	InitialiseMe (const G4ParticleDefinition *) final

virtual void	SetParticleAndCharge (const G4ParticleDefinition *, G4double q2) final

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

virtual	~G4VEmFluctuationModel ()

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

virtual G4double	Dispersion (const G4Material , const G4DynamicParticle , G4double tmax, 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

Detailed Description

Definition at line 57 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

◆ G4UniversalFluctuation()

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

explicit

Definition at line 62 of file G4UniversalFluctuation.cc.

 :G4VEmFluctuationModel(nam),
  particle(nullptr),
  minNumberInteractionsBohr(10.0),
  minLoss(10.*eV),
  nmaxCont(8.),
  rate(0.56),
  a0(42),
  fw(4.00)
{
  lastMaterial = nullptr;
  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct 
    = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall 
    = e1 = e2 = 0.0;
  m_Inv_particleMass = m_massrate = DBL_MAX;
  sizearray = 30;
  rndmarray = new G4double[30];
}

◆ ~G4UniversalFluctuation()

G4UniversalFluctuation::~G4UniversalFluctuation ( )

virtual

Definition at line 83 of file G4UniversalFluctuation.cc.

{
  delete [] rndmarray;
}

Member Function Documentation

◆ Dispersion()

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

overridevirtual

Implements G4VEmFluctuationModel.

Definition at line 287 of file G4UniversalFluctuation.cc.

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

◆ InitialiseMe()

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition * part )

finalvirtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 90 of file G4UniversalFluctuation.cc.

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

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

◆ SampleFluctuations()

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

overridevirtual

Implements G4VEmFluctuationModel.

Definition at line 105 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; }
  G4double 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 > electron_mass_c2) &&
      (meanLoss >= minNumberInteractionsBohr*tmax))
  {
    G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
                        (1.+m_massrate*(2.*gam+m_massrate)) ;
    if (tmaxkine <= 2.*tmax)   
    {
      electronDensity = material->GetElectronDensity();
      siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
                  * electronDensity * chargeSquare);
 
      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 {
 
        G4double neff = sn*sn;
        loss = meanLoss*G4RandGamma::shoot(rndmEngineF,neff,1.0)/neff;
      }
      //G4cout << "Gauss: " << loss << G4endl;
      return loss;
    }
  }
 
  // Glandz regime : initialisation
  //
  if (material != lastMaterial) {
    f1Fluct      = material->GetIonisation()->GetF1fluct();
    f2Fluct      = material->GetIonisation()->GetF2fluct();
    e1Fluct      = material->GetIonisation()->GetEnergy1fluct();
    e2Fluct      = material->GetIonisation()->GetEnergy2fluct();
    e1LogFluct   = material->GetIonisation()->GetLogEnergy1fluct();
    e2LogFluct   = material->GetIonisation()->GetLogEnergy2fluct();
    ipotFluct    = material->GetIonisation()->GetMeanExcitationEnergy();
    ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
    e0 = material->GetIonisation()->GetEnergy0fluct();
    esmall = 0.5*sqrt(e0*ipotFluct);  
    lastMaterial = material;   
  }
 
  // very small step or low-density material
  if(tmax <= e0) { return meanLoss; }
 
  // width correction for small cuts
  const G4double scaling = std::min(1.+0.5*CLHEP::keV/tmax,1.50);
  meanLoss /= scaling;
 
  G4double a1(0.0), a2(0.0), a3(0.0);
    
  loss = 0.0;
 
  e1 = e1Fluct;
  e2 = e2Fluct;
 
  if(tmax > ipotFluct) {
    G4double w2 = G4Log(2.*electron_mass_c2*beta2*gam2)-beta2;
 
    if(w2 > ipotLogFluct)  {
      if(w2 > e2LogFluct) {
    const G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
    a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
    a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
      } else {
    a1 = meanLoss*(1.-rate)/e1;
      }
      if(a1 < a0) { 
        const G4double fwnow = 0.5+(fw-0.5)*sqrt(a1/a0);
        a1 /= fwnow;
        e1 *= fwnow;
      } else {
        a1 /= fw;
        e1 = fw*e1Fluct;
      }
    }   
  }
 
  G4double w1 = tmax/e0;
  if(tmax > e0) {
    a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*G4Log(w1));
    if(a1+a2 <= 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); }
 
  // excitation of type 2
  if(a2 > 0.0) { AddExcitation(rndmEngineF, a2, e2, 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;
      }
 
    G4double w2 = alfa*e0;
    if(tmax > w2) {
      G4double w  = (tmax-w2)/tmax;
      G4int nnb = 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 += w2/(1.-w*rndmarray[k]); }
      }
    }
    if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
  }
 
  loss *= scaling;
 
  return loss;
 
}

Referenced by G4IonFluctuations::SampleFluctuations().

◆ SetParticleAndCharge()

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