G4UniversalFluctuation Class Reference

#include <G4UniversalFluctuation.hh>

Inheritance diagram for G4UniversalFluctuation:

Public Member Functions
	G4UniversalFluctuation (const G4String &nam="UniFluc")
virtual	~G4UniversalFluctuation ()
virtual G4double	SampleFluctuations (const G4Material *, const G4DynamicParticle *, G4double &, G4double &, G4double &)
virtual G4double	Dispersion (const G4Material *, const G4DynamicParticle *, G4double &, G4double &)
virtual void	InitialiseMe (const G4ParticleDefinition *)
virtual void	SetParticleAndCharge (const G4ParticleDefinition *, G4double q2)
Public Member Functions inherited from G4VEmFluctuationModel
	G4VEmFluctuationModel (const G4String &nam)
virtual	~G4VEmFluctuationModel ()
virtual G4double	SampleFluctuations (const G4Material *, 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)
G4String	GetName () const

Detailed Description

Definition at line 62 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

G4UniversalFluctuation()

G4UniversalFluctuation::G4UniversalFluctuation

(

const G4String &

nam = "UniFluc"

)

Definition at line 83 of file G4UniversalFluctuation.cc.

 :G4VEmFluctuationModel(nam),
  particle(0),
  minNumberInteractionsBohr(10.0),
  theBohrBeta2(50.0*keV/proton_mass_c2),
  minLoss(10.*eV),
  nmaxCont(16.),
  rate(0.55),
  fw(4.)
{
  lastMaterial = 0;
 
  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct 
    = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall 
    = e1 = e2 = 0;
 
}

~G4UniversalFluctuation()

G4UniversalFluctuation::~G4UniversalFluctuation ( )

virtual

Definition at line 103 of file G4UniversalFluctuation.cc.

{}

Member Function Documentation

Dispersion()

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

virtual

Implements G4VEmFluctuationModel.

Definition at line 348 of file G4UniversalFluctuation.cc.

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

InitialiseMe()

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition *  part )

virtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 108 of file G4UniversalFluctuation.cc.

{
  particle       = part;
  particleMass   = part->GetPDGMass();
  G4double q     = part->GetPDGCharge()/eplus;
  chargeSquare   = q*q;
}

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

SampleFluctuations()

G4double G4UniversalFluctuation::SampleFluctuations ( const G4Material *  material, const G4DynamicParticle *  dp, G4double &  tmax, G4double &  length, G4double &  averageLoss  )

virtual

Implements G4VEmFluctuationModel.

Definition at line 118 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)
  //
  G4double meanLoss = averageLoss;
  G4double tkin  = dp->GetKineticEnergy();
  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
  if (meanLoss < minLoss) { return meanLoss; }
 
  if(!particle) { InitialiseMe(dp->GetDefinition()); }
  
  G4double tau   = tkin/particleMass;
  G4double gam   = tau + 1.0;
  G4double gam2  = gam*gam;
  G4double beta2 = tau*(tau + 2.0)/gam2;
 
  G4double loss(0.), siga(0.);
  
  // Gaussian regime
  // for heavy particles only and conditions
  // for Gauusian fluct. has been changed 
  //
  if ((particleMass > electron_mass_c2) &&
      (meanLoss >= minNumberInteractionsBohr*tmax))
  {
    G4double massrate = electron_mass_c2/particleMass ;
    G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
                        (1.+massrate*(2.*gam+massrate)) ;
    if (tmaxkine <= 2.*tmax)   
    {
      electronDensity = material->GetElectronDensity();
      siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
          * electronDensity * chargeSquare);
 
     
      G4double sn = meanLoss/siga;
  
      // thick target case 
      if (sn >= 2.0) {
 
    G4double twomeanLoss = meanLoss + meanLoss;
    do {
      loss = G4RandGauss::shoot(meanLoss,siga);
    } while (0.0 > loss || twomeanLoss < loss);
 
    // Gamma distribution
      } else {
 
    G4double neff = sn*sn;
    loss = meanLoss*CLHEP::RandGamma::shoot(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; }
 
  G4double losstot = 0.;
  G4int    nstep = 1;
  if(meanLoss < 25.*ipotFluct)
    {
      if(G4UniformRand() < 0.04*meanLoss/ipotFluct)
    { nstep = 1; }
      else
    { 
      nstep = 2;
      meanLoss *= 0.5; 
    }
    }
 
  for (G4int istep=0; istep < nstep; ++istep) {
    
    loss = 0.;
 
    G4double a1 = 0. , a2 = 0., a3 = 0. ;
 
    if(tmax > ipotFluct) {
      G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
 
      if(w2 > ipotLogFluct)  {
    G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
    a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
    if(w2 > e2LogFluct) {
      a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
    }
    if(a1 < nmaxCont) { 
      //small energy loss
      G4double sa1 = sqrt(a1);
      if(G4UniformRand() < exp(-sa1))
        {
          e1 = esmall;
          a1 = meanLoss*(1.-rate)/e1;
          a2 = 0.;
          e2 = e2Fluct;
        } 
      else
        {
          a1 = sa1 ;    
          e1 = sa1*e1Fluct;
          e2 = e2Fluct;
        }
 
    } else {
        //not small energy loss
        //correction to get better fwhm value
        a1 /= fw;
        e1 = fw*e1Fluct;
        e2 = e2Fluct;
    }
      }   
    }
 
    G4double w1 = tmax/e0;
    if(tmax > e0) {
      a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
    }
    //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont  
    G4double emean = 0.;
    G4double sig2e = 0., sige = 0.;
    G4double p1 = 0., p2 = 0., p3 = 0.;
 
    // excitation of type 1
    if(a1 > nmaxCont)
      {
    emean += a1*e1;
    sig2e += a1*e1*e1;
      }
    else if(a1 > 0.)
      {
    p1 = G4double(G4Poisson(a1));
    loss += p1*e1;
    if(p1 > 0.) {
      loss += (1.-2.*G4UniformRand())*e1;
    }
      }
 
 
    // excitation of type 2
    if(a2 > nmaxCont)
      {
    emean += a2*e2;
    sig2e += a2*e2*e2;
      }
    else if(a2 > 0.)
      {
    p2 = G4double(G4Poisson(a2));
    loss += p2*e2;
    if(p2 > 0.) 
      loss += (1.-2.*G4UniformRand())*e2;
      }
 
    if(emean > 0.)
      {
    sige   = sqrt(sig2e);
    loss += max(0.,G4RandGauss::shoot(emean,sige));
      }
 
    // ionisation 
    G4double lossc = 0.;
    if(a3 > 0.) {
      emean = 0.;
      sig2e = 0.;
      sige = 0.;
      p3 = a3;
      G4double alfa = 1.;
      if(a3 > nmaxCont)
    {
      alfa            = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
      G4double alfa1  = alfa*log(alfa)/(alfa-1.);
      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;
      G4double w  = (tmax-w2)/tmax;
      G4int nb = G4Poisson(p3);
      if(nb > 0) {
    for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
      }
    }
 
    if(emean > 0.)
      {
    sige   = sqrt(sig2e);
    lossc += max(0.,G4RandGauss::shoot(emean,sige));
      }
 
    loss += lossc;
 
    losstot += loss;
  }
  //G4cout << "Vavilov: " << losstot << "  Nstep= " << nstep << G4endl;
              
  return losstot;
 
}

Referenced by G4IonFluctuations::SampleFluctuations().

SetParticleAndCharge()

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

virtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 370 of file G4UniversalFluctuation.cc.

{
  if(part != particle) {
    particle       = part;
    particleMass   = part->GetPDGMass();
  }
  chargeSquare = q2;
}

Referenced by G4IonFluctuations::SetParticleAndCharge().

---

The documentation for this class was generated from the following files:

• G4UniversalFluctuation.hh
• G4UniversalFluctuation.cc