Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) 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 ()
 
const G4StringGetName () const
 

Detailed Description

Definition at line 64 of file G4UniversalFluctuation.hh.

Constructor & Destructor Documentation

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

Definition at line 86 of file G4UniversalFluctuation.cc.

88  particle(nullptr),
89  minNumberInteractionsBohr(10.0),
90  minLoss(10.*eV),
91  nmaxCont(16.),
92  rate(0.55),
93  fw(4.)
94 {
95  lastMaterial = 0;
96 
97  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
98  = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
99  = e1 = e2 = 0;
100  m_Inv_particleMass = m_massrate = DBL_MAX;
101  sizearray = 30;
102  rndmarray = new G4double[30];
103 }
static constexpr double eV
Definition: G4SIunits.hh:215
double G4double
Definition: G4Types.hh:76
G4VEmFluctuationModel(const G4String &nam)
#define DBL_MAX
Definition: templates.hh:83
G4UniversalFluctuation::~G4UniversalFluctuation ( )
virtual

Definition at line 107 of file G4UniversalFluctuation.cc.

108 {
109  delete [] rndmarray;
110 }

Member Function Documentation

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

Implements G4VEmFluctuationModel.

Definition at line 335 of file G4UniversalFluctuation.cc.

340 {
341  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
342 
343  electronDensity = material->GetElectronDensity();
344 
345  G4double gam = (dp->GetKineticEnergy())*m_Inv_particleMass + 1.0;
346  G4double beta2 = 1.0 - 1.0/(gam*gam);
347 
348  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
349  * electronDensity * chargeSquare;
350 
351  return siga;
352 }
G4double GetKineticEnergy() const
G4ParticleDefinition * GetDefinition() const
G4double GetElectronDensity() const
Definition: G4Material.hh:217
virtual void InitialiseMe(const G4ParticleDefinition *) final
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4UniversalFluctuation::InitialiseMe ( const G4ParticleDefinition part)
finalvirtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 114 of file G4UniversalFluctuation.cc.

115 {
116  particle = part;
117  particleMass = part->GetPDGMass();
118  G4double q = part->GetPDGCharge()/eplus;
119 
120  // Derived quantities
121  m_Inv_particleMass = 1.0 / particleMass;
122  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
123  chargeSquare = q*q;
124 }
static constexpr double eplus
Definition: G4SIunits.hh:199
float electron_mass_c2
Definition: hepunit.py:274
G4double GetPDGMass() const
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const

Here is the call graph for this function:

Here is the caller graph for this function:

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

Implements G4VEmFluctuationModel.

Definition at line 129 of file G4UniversalFluctuation.cc.

134 {
135  // Calculate actual loss from the mean loss.
136  // The model used to get the fluctuations is essentially the same
137  // as in Glandz in Geant3 (Cern program library W5013, phys332).
138  // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
139 
140  // shortcut for very small loss or from a step nearly equal to the range
141  // (out of validity of the model)
142  //
143  G4double meanLoss = averageLoss;
144  G4double tkin = dp->GetKineticEnergy();
145  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
146  if (meanLoss < minLoss) { return meanLoss; }
147 
148  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
149 
150  CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
151 
152  G4double tau = tkin * m_Inv_particleMass;
153  G4double gam = tau + 1.0;
154  G4double gam2 = gam*gam;
155  G4double beta2 = tau*(tau + 2.0)/gam2;
156 
157  G4double loss(0.), siga(0.);
158 
159  const G4Material* material = couple->GetMaterial();
160 
161  // Gaussian regime
162  // for heavy particles only and conditions
163  // for Gauusian fluct. has been changed
164  //
165  if ((particleMass > electron_mass_c2) &&
166  (meanLoss >= minNumberInteractionsBohr*tmax))
167  {
168  G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
169  (1.+m_massrate*(2.*gam+m_massrate)) ;
170  if (tmaxkine <= 2.*tmax)
171  {
172  electronDensity = material->GetElectronDensity();
173  siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
174  * electronDensity * chargeSquare);
175 
176  G4double sn = meanLoss/siga;
177 
178  // thick target case
179  if (sn >= 2.0) {
180 
181  G4double twomeanLoss = meanLoss + meanLoss;
182  do {
183  loss = G4RandGauss::shoot(rndmEngineF,meanLoss,siga);
184  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
185  } while (0.0 > loss || twomeanLoss < loss);
186 
187  // Gamma distribution
188  } else {
189 
190  G4double neff = sn*sn;
191  loss = meanLoss*G4RandGamma::shoot(rndmEngineF,neff,1.0)/neff;
192  }
193  //G4cout << "Gauss: " << loss << G4endl;
194  return loss;
195  }
196  }
197 
198  // Glandz regime : initialisation
199  //
200  if (material != lastMaterial) {
201  f1Fluct = material->GetIonisation()->GetF1fluct();
202  f2Fluct = material->GetIonisation()->GetF2fluct();
203  e1Fluct = material->GetIonisation()->GetEnergy1fluct();
204  e2Fluct = material->GetIonisation()->GetEnergy2fluct();
205  e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
206  e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
207  ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
208  ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
209  e0 = material->GetIonisation()->GetEnergy0fluct();
210  esmall = 0.5*sqrt(e0*ipotFluct);
211  lastMaterial = material;
212  }
213 
214  // very small step or low-density material
215  if(tmax <= e0) { return meanLoss; }
216 
217  G4double losstot = 0.;
218  G4int nstep = 1;
219  if(meanLoss < 25.*ipotFluct)
220  {
221  if(rndmEngineF->flat()*ipotFluct< 0.04*meanLoss)
222  { nstep = 1; }
223  else
224  {
225  nstep = 2;
226  meanLoss *= 0.5;
227  }
228  }
229 
230  G4double a1, a2, a3;
231  for (G4int istep=0; istep < nstep; ++istep) {
232 
233  loss = a1 = a2 = a3 = 0.;
234  e1 = e1Fluct;
235  e2 = e2Fluct;
236 
237  if(tmax > ipotFluct) {
238  G4double w2 = G4Log(2.*electron_mass_c2*beta2*gam2)-beta2;
239 
240  if(w2 > ipotLogFluct) {
241  if(w2 > e2LogFluct) {
242  G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
243  a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
244  a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
245  } else {
246  a1 = meanLoss*(1.-rate)/e1;
247  }
248 
249  if(a1 < nmaxCont) {
250  //small energy loss
251  G4double sa1 = sqrt(a1);
252  if(rndmEngineF->flat() < G4Exp(-sa1))
253  {
254  e1 = esmall;
255  a1 = meanLoss*(1.-rate)/e1;
256  a2 = 0.;
257  }
258  else
259  {
260  a1 = sa1 ;
261  e1 = sa1*e1Fluct;
262  }
263 
264  } else {
265  //not small energy loss
266  //correction to get better fwhm value
267  a1 /= fw;
268  e1 = fw*e1Fluct;
269  }
270  }
271  }
272 
273  G4double w1 = tmax/e0;
274  if(tmax > e0) {
275  a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*G4Log(w1));
276  if(a1+a2 <= 0.) {
277  a3 /= rate;
278  }
279  }
280  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
281  G4double emean = 0.;
282  G4double sig2e = 0.;
283 
284  // excitation of type 1
285  AddExcitation(rndmEngineF, a1, e1, emean, loss, sig2e);
286 
287  // excitation of type 2
288  AddExcitation(rndmEngineF, a2, e2, emean, loss, sig2e);
289 
290  if(emean > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
291 
292  // ionisation
293  if(a3 > 0.) {
294  emean = 0.;
295  sig2e = 0.;
296  G4double p3 = a3;
297  G4double alfa = 1.;
298  if(a3 > nmaxCont)
299  {
300  alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
301  G4double alfa1 = alfa*G4Log(alfa)/(alfa-1.);
302  G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
303  emean += namean*e0*alfa1;
304  sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
305  p3 = a3-namean;
306  }
307 
308  G4double w2 = alfa*e0;
309  if(tmax > w2) {
310  G4double w = (tmax-w2)/tmax;
311  G4int nb = G4Poisson(p3);
312  if(nb > 0) {
313  if(nb > sizearray) {
314  sizearray = nb;
315  delete [] rndmarray;
316  rndmarray = new G4double[nb];
317  }
318  rndmEngineF->flatArray(nb, rndmarray);
319  for (G4int k=0; k<nb; ++k) { loss += w2/(1.-w*rndmarray[k]); }
320  }
321  }
322  if(emean > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
323  }
324  losstot += loss;
325  }
326  //G4cout << "Vavilov: " << losstot << " Nstep= " << nstep << G4endl;
327 
328  return losstot;
329 
330 }
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetKineticEnergy() const
G4double GetEnergy2fluct() const
G4double GetLogEnergy2fluct() const
G4ParticleDefinition * GetDefinition() const
double C(double temp)
G4double GetLogMeanExcEnergy() const
int G4int
Definition: G4Types.hh:78
string material
Definition: eplot.py:19
G4double GetEnergy0fluct() const
G4double GetElectronDensity() const
Definition: G4Material.hh:217
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual void InitialiseMe(const G4ParticleDefinition *) final
G4double GetLogEnergy1fluct() const
G4double GetMeanExcitationEnergy() const
G4double GetF2fluct() const
double G4double
Definition: G4Types.hh:76
G4double GetF1fluct() const
const G4Material * GetMaterial() const
G4double GetEnergy1fluct() const

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4VEmFluctuationModel.

Definition at line 357 of file G4UniversalFluctuation.cc.

359 {
360  if(part != particle) {
361  particle = part;
362  particleMass = part->GetPDGMass();
363 
364  // Derived quantities
365  if( particleMass != 0.0 ){
366  m_Inv_particleMass = 1.0 / particleMass;
367  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
368  }else{
369  m_Inv_particleMass = DBL_MAX;
370  m_massrate = DBL_MAX;
371  }
372  }
373  chargeSquare = q2;
374 }
float electron_mass_c2
Definition: hepunit.py:274
G4double GetPDGMass() const
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

Here is the caller graph for this function:


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