G4UniversalFluctuation Class Reference

#include <G4UniversalFluctuation.hh>

Inheritance diagram for G4UniversalFluctuation:

Collaboration diagram for G4UniversalFluctuation:

Public Member Functions
	G4UniversalFluctuation (const G4String &nam="UniFluc")
virtual	~G4UniversalFluctuation ()
virtual G4double	SampleFluctuations (const G4MaterialCutsCouple *, 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 ()
const G4String &	GetName () const

Private Member Functions
G4UniversalFluctuation &	operator= (const G4UniversalFluctuation &right)
	G4UniversalFluctuation (const G4UniversalFluctuation &)

Private Attributes
const G4ParticleDefinition *	particle
const G4Material *	lastMaterial
G4double	particleMass
G4double	m_Inv_particleMass
G4double	m_massrate
G4double	chargeSquare
G4double	ipotFluct
G4double	electronDensity
G4double	f1Fluct
G4double	f2Fluct
G4double	e1Fluct
G4double	e2Fluct
G4double	e1LogFluct
G4double	e2LogFluct
G4double	ipotLogFluct
G4double	e0
G4double	esmall
G4double	e1
G4double	e2
G4double	minNumberInteractionsBohr
G4double	theBohrBeta2
G4double	minLoss
G4double	nmaxCont
G4double	rate
G4double	fw
G4int	sizearray
G4double *	rndmarray

Detailed Description

Definition at line 62 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

G4UniversalFluctuation() [1/2]

G4UniversalFluctuation::G4UniversalFluctuation

(

const G4String &

nam = "UniFluc"

)

Definition at line 86 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;
  m_Inv_particleMass = m_massrate = DBL_MAX;
  sizearray = 30;
  rndmarray = new G4double[30];
}

~G4UniversalFluctuation()

G4UniversalFluctuation::~G4UniversalFluctuation ( )

virtual

Definition at line 108 of file G4UniversalFluctuation.cc.

{
  delete [] rndmarray;
}

G4UniversalFluctuation() [2/2]

G4UniversalFluctuation::G4UniversalFluctuation ( const G4UniversalFluctuation &  )

private

Member Function Documentation

Dispersion()

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

virtual

Implements G4VEmFluctuationModel.

Definition at line 370 of file G4UniversalFluctuation.cc.

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

  electronDensity = material->GetElectronDensity();

  G4double gam   = (dp->GetKineticEnergy())*m_Inv_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;
}

Here is the call graph for this function:

InitialiseMe()

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition *  part )

virtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 115 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

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

private

SampleFluctuations()

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

virtual

Implements G4VEmFluctuationModel.

Definition at line 130 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(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }

  CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
  
  G4double tau   = tkin * m_Inv_particleMass;            
  G4double gam   = tau + 1.0;
  G4double gam2  = gam*gam;
  G4double beta2 = tau*(tau + 2.0)/gam2;

  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);

      G4double sn = meanLoss/siga;
  
      // thick target case 
      if (sn >= 2.0) {

        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; }

  G4double losstot = 0.;
  G4int    nstep = 1;
  if(meanLoss < 25.*ipotFluct)
    {
      if(rndmEngineF->flat()*ipotFluct< 0.04*meanLoss)          
        { 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 = G4Log(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(rndmEngineF->flat() < G4Exp(-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*G4Log(w1));
      if(a1+a2 <= 0.) { 
        a3 /= rate; 
      }
    }
    //'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.*rndmEngineF->flat())*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.*rndmEngineF->flat())*e2;
      }
    if(emean > 0.)
      {
        sige   = sqrt(sig2e);
        loss += max(0.,G4RandGauss::shoot(rndmEngineF,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*G4Log(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;
      const G4int nb = G4Poisson(p3);
      if(nb > 0) {
        if(nb > sizearray) {
          sizearray = nb;
          delete [] rndmarray;
          rndmarray = new G4double[nb];
        }
        rndmEngineF->flatArray(nb, rndmarray);
        for (G4int k=0; k<nb; ++k) { lossc += w2/(1.-w*rndmarray[k]); }
      }

    if(emean > 0.)
      {
        sige   = sqrt(sig2e);
        lossc += max(0.,G4RandGauss::shoot(rndmEngineF,emean,sige));
      }
    }

    loss += lossc;

    losstot += loss;
  }
  //G4cout << "Vavilov: " << losstot << "  Nstep= " << nstep << G4endl;

  return losstot;

}

Here is the call graph for this function:

Here is the caller graph for this function:

SetParticleAndCharge()

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

virtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 392 of file G4UniversalFluctuation.cc.

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

    // Derived quantities
    if( particleMass != 0.0 ){
      m_Inv_particleMass = 1.0 / particleMass;
      m_massrate = electron_mass_c2 * m_Inv_particleMass ;
    }else{
      m_Inv_particleMass = DBL_MAX;
      m_massrate = DBL_MAX;
    }
  }
  chargeSquare = q2;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

chargeSquare

G4double G4UniversalFluctuation::chargeSquare

private

Definition at line 101 of file G4UniversalFluctuation.hh.

e0

G4double G4UniversalFluctuation::e0

private

Definition at line 114 of file G4UniversalFluctuation.hh.

e1

G4double G4UniversalFluctuation::e1

private

Definition at line 117 of file G4UniversalFluctuation.hh.

e1Fluct

G4double G4UniversalFluctuation::e1Fluct

private

Definition at line 109 of file G4UniversalFluctuation.hh.

e1LogFluct

G4double G4UniversalFluctuation::e1LogFluct

private

Definition at line 111 of file G4UniversalFluctuation.hh.

e2

G4double G4UniversalFluctuation::e2

private

Definition at line 117 of file G4UniversalFluctuation.hh.

e2Fluct

G4double G4UniversalFluctuation::e2Fluct

private

Definition at line 110 of file G4UniversalFluctuation.hh.

e2LogFluct

G4double G4UniversalFluctuation::e2LogFluct

private

Definition at line 112 of file G4UniversalFluctuation.hh.

electronDensity

G4double G4UniversalFluctuation::electronDensity

private

Definition at line 105 of file G4UniversalFluctuation.hh.

esmall

G4double G4UniversalFluctuation::esmall

private

Definition at line 115 of file G4UniversalFluctuation.hh.

f1Fluct

G4double G4UniversalFluctuation::f1Fluct

private

Definition at line 107 of file G4UniversalFluctuation.hh.

f2Fluct

G4double G4UniversalFluctuation::f2Fluct

private

Definition at line 108 of file G4UniversalFluctuation.hh.

fw

G4double G4UniversalFluctuation::fw

private

Definition at line 123 of file G4UniversalFluctuation.hh.

ipotFluct

G4double G4UniversalFluctuation::ipotFluct

private

Definition at line 104 of file G4UniversalFluctuation.hh.

ipotLogFluct

G4double G4UniversalFluctuation::ipotLogFluct

private

Definition at line 113 of file G4UniversalFluctuation.hh.

lastMaterial

const G4Material* G4UniversalFluctuation::lastMaterial

private

Definition at line 94 of file G4UniversalFluctuation.hh.

m_Inv_particleMass

G4double G4UniversalFluctuation::m_Inv_particleMass

private

Definition at line 99 of file G4UniversalFluctuation.hh.

m_massrate

G4double G4UniversalFluctuation::m_massrate

private

Definition at line 100 of file G4UniversalFluctuation.hh.

minLoss

G4double G4UniversalFluctuation::minLoss

private

Definition at line 121 of file G4UniversalFluctuation.hh.

minNumberInteractionsBohr

G4double G4UniversalFluctuation::minNumberInteractionsBohr

private

Definition at line 119 of file G4UniversalFluctuation.hh.

nmaxCont

G4double G4UniversalFluctuation::nmaxCont

private

Definition at line 122 of file G4UniversalFluctuation.hh.

particle

const G4ParticleDefinition* G4UniversalFluctuation::particle

private

Definition at line 93 of file G4UniversalFluctuation.hh.

particleMass

G4double G4UniversalFluctuation::particleMass

private

Definition at line 96 of file G4UniversalFluctuation.hh.

rate

G4double G4UniversalFluctuation::rate

private

Definition at line 123 of file G4UniversalFluctuation.hh.

rndmarray

G4double* G4UniversalFluctuation::rndmarray

private

Definition at line 126 of file G4UniversalFluctuation.hh.

sizearray

G4int G4UniversalFluctuation::sizearray

private

Definition at line 125 of file G4UniversalFluctuation.hh.

theBohrBeta2

G4double G4UniversalFluctuation::theBohrBeta2

private

Definition at line 120 of file G4UniversalFluctuation.hh.

---

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

• G4UniversalFluctuation.hh
• G4UniversalFluctuation.cc