G4UniversalFluctuation Class Reference

#include <G4UniversalFluctuation.hh>

Inheritance diagram for G4UniversalFluctuation:

G4VEmFluctuationModel

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)

Detailed Description

Definition at line 62 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

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

Definition at line 83 of file G4UniversalFluctuation.cc.

00084  :G4VEmFluctuationModel(nam),
00085   particle(0),
00086   minNumberInteractionsBohr(10.0),
00087   theBohrBeta2(50.0*keV/proton_mass_c2),
00088   minLoss(10.*eV),
00089   nmaxCont(16.),
00090   rate(0.55),
00091   fw(4.)
00092 {
00093   lastMaterial = 0;
00094 
00095   particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct 
00096     = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall 
00097     = e1 = e2 = 0;
00098 
00099 }

G4UniversalFluctuation::~G4UniversalFluctuation (  )  [virtual]

Definition at line 103 of file G4UniversalFluctuation.cc.

00104 {}

Member Function Documentation

G4double G4UniversalFluctuation::Dispersion ( const G4Material *  ,
const G4DynamicParticle *  ,
G4double &  ,
G4double &   
) [virtual]

Implements G4VEmFluctuationModel.

Definition at line 348 of file G4UniversalFluctuation.cc.

References G4InuclParticleNames::gam, G4DynamicParticle::GetDefinition(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), and InitialiseMe().

00353 {
00354   if(!particle) { InitialiseMe(dp->GetDefinition()); }
00355 
00356   electronDensity = material->GetElectronDensity();
00357 
00358   G4double gam   = (dp->GetKineticEnergy())/particleMass + 1.0;
00359   G4double beta2 = 1.0 - 1.0/(gam*gam);
00360 
00361   G4double siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00362                  * electronDensity * chargeSquare;
00363 
00364   return siga;
00365 }

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition *   )  [virtual]

Reimplemented from G4VEmFluctuationModel.

Definition at line 108 of file G4UniversalFluctuation.cc.

References G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetPDGMass().

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

00109 {
00110   particle       = part;
00111   particleMass   = part->GetPDGMass();
00112   G4double q     = part->GetPDGCharge()/eplus;
00113   chargeSquare   = q*q;
00114 }

G4double G4UniversalFluctuation::SampleFluctuations ( const G4Material *  ,
const G4DynamicParticle *  ,
G4double &  ,
G4double &  ,
G4double &   
) [virtual]

Implements G4VEmFluctuationModel.

Definition at line 118 of file G4UniversalFluctuation.cc.

References G4Poisson(), G4UniformRand, G4InuclParticleNames::gam, G4DynamicParticle::GetDefinition(), G4Material::GetElectronDensity(), G4IonisParamMat::GetEnergy0fluct(), G4IonisParamMat::GetEnergy1fluct(), G4IonisParamMat::GetEnergy2fluct(), G4IonisParamMat::GetF1fluct(), G4IonisParamMat::GetF2fluct(), G4Material::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamMat::GetLogEnergy1fluct(), G4IonisParamMat::GetLogEnergy2fluct(), G4IonisParamMat::GetLogMeanExcEnergy(), G4IonisParamMat::GetMeanExcitationEnergy(), and InitialiseMe().

Referenced by G4IonFluctuations::SampleFluctuations().

00123 {
00124   // Calculate actual loss from the mean loss.
00125   // The model used to get the fluctuations is essentially the same
00126   // as in Glandz in Geant3 (Cern program library W5013, phys332).
00127   // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
00128 
00129   // shortcut for very small loss or from a step nearly equal to the range
00130   // (out of validity of the model)
00131   //
00132   G4double meanLoss = averageLoss;
00133   G4double tkin  = dp->GetKineticEnergy();
00134   //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
00135   if (meanLoss < minLoss) { return meanLoss; }
00136 
00137   if(!particle) { InitialiseMe(dp->GetDefinition()); }
00138   
00139   G4double tau   = tkin/particleMass;
00140   G4double gam   = tau + 1.0;
00141   G4double gam2  = gam*gam;
00142   G4double beta2 = tau*(tau + 2.0)/gam2;
00143 
00144   G4double loss(0.), siga(0.);
00145   
00146   // Gaussian regime
00147   // for heavy particles only and conditions
00148   // for Gauusian fluct. has been changed 
00149   //
00150   if ((particleMass > electron_mass_c2) &&
00151       (meanLoss >= minNumberInteractionsBohr*tmax))
00152   {
00153     G4double massrate = electron_mass_c2/particleMass ;
00154     G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
00155                         (1.+massrate*(2.*gam+massrate)) ;
00156     if (tmaxkine <= 2.*tmax)   
00157     {
00158       electronDensity = material->GetElectronDensity();
00159       siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00160                   * electronDensity * chargeSquare);
00161 
00162      
00163       G4double sn = meanLoss/siga;
00164   
00165       // thick target case 
00166       if (sn >= 2.0) {
00167 
00168         G4double twomeanLoss = meanLoss + meanLoss;
00169         do {
00170           loss = G4RandGauss::shoot(meanLoss,siga);
00171         } while (0.0 > loss || twomeanLoss < loss);
00172 
00173         // Gamma distribution
00174       } else {
00175 
00176         G4double neff = sn*sn;
00177         loss = meanLoss*CLHEP::RandGamma::shoot(neff,1.0)/neff;
00178       }
00179       //G4cout << "Gauss: " << loss << G4endl;
00180       return loss;
00181     }
00182   }
00183 
00184   // Glandz regime : initialisation
00185   //
00186   if (material != lastMaterial) {
00187     f1Fluct      = material->GetIonisation()->GetF1fluct();
00188     f2Fluct      = material->GetIonisation()->GetF2fluct();
00189     e1Fluct      = material->GetIonisation()->GetEnergy1fluct();
00190     e2Fluct      = material->GetIonisation()->GetEnergy2fluct();
00191     e1LogFluct   = material->GetIonisation()->GetLogEnergy1fluct();
00192     e2LogFluct   = material->GetIonisation()->GetLogEnergy2fluct();
00193     ipotFluct    = material->GetIonisation()->GetMeanExcitationEnergy();
00194     ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
00195     e0 = material->GetIonisation()->GetEnergy0fluct();
00196     esmall = 0.5*sqrt(e0*ipotFluct);  
00197     lastMaterial = material;   
00198   }
00199 
00200   // very small step or low-density material
00201   if(tmax <= e0) { return meanLoss; }
00202 
00203   G4double losstot = 0.;
00204   G4int    nstep = 1;
00205   if(meanLoss < 25.*ipotFluct)
00206     {
00207       if(G4UniformRand() < 0.04*meanLoss/ipotFluct)
00208         { nstep = 1; }
00209       else
00210         { 
00211           nstep = 2;
00212           meanLoss *= 0.5; 
00213         }
00214     }
00215 
00216   for (G4int istep=0; istep < nstep; ++istep) {
00217     
00218     loss = 0.;
00219 
00220     G4double a1 = 0. , a2 = 0., a3 = 0. ;
00221 
00222     if(tmax > ipotFluct) {
00223       G4double w2 = log(2.*electron_mass_c2*beta2*gam2)-beta2;
00224 
00225       if(w2 > ipotLogFluct)  {
00226         G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
00227         a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
00228         if(w2 > e2LogFluct) {
00229           a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
00230         }
00231         if(a1 < nmaxCont) { 
00232           //small energy loss
00233           G4double sa1 = sqrt(a1);
00234           if(G4UniformRand() < exp(-sa1))
00235             {
00236               e1 = esmall;
00237               a1 = meanLoss*(1.-rate)/e1;
00238               a2 = 0.;
00239               e2 = e2Fluct;
00240             } 
00241           else
00242             {
00243               a1 = sa1 ;    
00244               e1 = sa1*e1Fluct;
00245               e2 = e2Fluct;
00246             }
00247 
00248         } else {
00249             //not small energy loss
00250             //correction to get better fwhm value
00251             a1 /= fw;
00252             e1 = fw*e1Fluct;
00253             e2 = e2Fluct;
00254         }
00255       }   
00256     }
00257 
00258     G4double w1 = tmax/e0;
00259     if(tmax > e0) {
00260       a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*log(w1));
00261     }
00262     //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont  
00263     G4double emean = 0.;
00264     G4double sig2e = 0., sige = 0.;
00265     G4double p1 = 0., p2 = 0., p3 = 0.;
00266  
00267     // excitation of type 1
00268     if(a1 > nmaxCont)
00269       {
00270         emean += a1*e1;
00271         sig2e += a1*e1*e1;
00272       }
00273     else if(a1 > 0.)
00274       {
00275         p1 = G4double(G4Poisson(a1));
00276         loss += p1*e1;
00277         if(p1 > 0.) {
00278           loss += (1.-2.*G4UniformRand())*e1;
00279         }
00280       }
00281 
00282 
00283     // excitation of type 2
00284     if(a2 > nmaxCont)
00285       {
00286         emean += a2*e2;
00287         sig2e += a2*e2*e2;
00288       }
00289     else if(a2 > 0.)
00290       {
00291         p2 = G4double(G4Poisson(a2));
00292         loss += p2*e2;
00293         if(p2 > 0.) 
00294           loss += (1.-2.*G4UniformRand())*e2;
00295       }
00296 
00297     if(emean > 0.)
00298       {
00299         sige   = sqrt(sig2e);
00300         loss += max(0.,G4RandGauss::shoot(emean,sige));
00301       }
00302 
00303     // ionisation 
00304     G4double lossc = 0.;
00305     if(a3 > 0.) {
00306       emean = 0.;
00307       sig2e = 0.;
00308       sige = 0.;
00309       p3 = a3;
00310       G4double alfa = 1.;
00311       if(a3 > nmaxCont)
00312         {
00313           alfa            = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
00314           G4double alfa1  = alfa*log(alfa)/(alfa-1.);
00315           G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
00316           emean          += namean*e0*alfa1;
00317           sig2e          += e0*e0*namean*(alfa-alfa1*alfa1);
00318           p3              = a3-namean;
00319         }
00320 
00321       G4double w2 = alfa*e0;
00322       G4double w  = (tmax-w2)/tmax;
00323       G4int nb = G4Poisson(p3);
00324       if(nb > 0) {
00325         for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
00326       }
00327     }
00328 
00329     if(emean > 0.)
00330       {
00331         sige   = sqrt(sig2e);
00332         lossc += max(0.,G4RandGauss::shoot(emean,sige));
00333       }
00334 
00335     loss += lossc;
00336 
00337     losstot += loss;
00338   }
00339   //G4cout << "Vavilov: " << losstot << "  Nstep= " << nstep << G4endl;
00340               
00341   return losstot;
00342 
00343 }

void G4UniversalFluctuation::SetParticleAndCharge ( const G4ParticleDefinition *  ,
G4double  q2 
) [virtual]

Reimplemented from G4VEmFluctuationModel.

Definition at line 370 of file G4UniversalFluctuation.cc.

References G4ParticleDefinition::GetPDGMass().

Referenced by G4IonFluctuations::SetParticleAndCharge().

00372 {
00373   if(part != particle) {
00374     particle       = part;
00375     particleMass   = part->GetPDGMass();
00376   }
00377   chargeSquare = q2;
00378 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:37 2013 for Geant4 by  doxygen 1.4.7