ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UniversalFluctuation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UniversalFluctuation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4UniversalFluctuation
33 //
34 // Author: V. Ivanchenko for Laszlo Urban
35 //
36 // Creation date: 03.01.2002
37 //
38 // Modifications:
39 //
40 //
41 
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 
45 #include "G4UniversalFluctuation.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh"
49 #include "G4Poisson.hh"
50 #include "G4Step.hh"
51 #include "G4Material.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4DynamicParticle.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4Log.hh"
56 #include "G4Exp.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 using namespace std;
61 
62 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
63  :G4VEmFluctuationModel(nam),
64  particle(nullptr),
65  minNumberInteractionsBohr(10.0),
66  minLoss(10.*eV),
67  nmaxCont(16.),
68  rate(0.56),
69  a0(50.),
70  fw(4.00)
71 {
72  lastMaterial = nullptr;
73  particleMass = chargeSquare = ipotFluct = electronDensity = f1Fluct = f2Fluct
74  = e1Fluct = e2Fluct = e1LogFluct = e2LogFluct = ipotLogFluct = e0 = esmall
75  = e1 = e2 = 0.0;
76  m_Inv_particleMass = m_massrate = DBL_MAX;
77  sizearray = 30;
78  rndmarray = new G4double[30];
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
83 G4UniversalFluctuation::~G4UniversalFluctuation()
84 {
85  delete [] rndmarray;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
90 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
91 {
92  particle = part;
93  particleMass = part->GetPDGMass();
94  G4double q = part->GetPDGCharge()/eplus;
95 
96  // Derived quantities
97  m_Inv_particleMass = 1.0 / particleMass;
98  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
99  chargeSquare = q*q;
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 G4double
105 G4UniversalFluctuation::SampleFluctuations(const G4MaterialCutsCouple* couple,
106  const G4DynamicParticle* dp,
107  G4double tmax,
108  G4double length,
109  G4double averageLoss)
110 {
111  // Calculate actual loss from the mean loss.
112  // The model used to get the fluctuations is essentially the same
113  // as in Glandz in Geant3 (Cern program library W5013, phys332).
114  // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
115 
116  // shortcut for very small loss or from a step nearly equal to the range
117  // (out of validity of the model)
118  //
119  G4double meanLoss = averageLoss;
120  G4double tkin = dp->GetKineticEnergy();
121  //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
122  if (meanLoss < minLoss) { return meanLoss; }
123 
124  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
125 
126  CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
127 
128  G4double tau = tkin * m_Inv_particleMass;
129  G4double gam = tau + 1.0;
130  G4double gam2 = gam*gam;
131  G4double beta2 = tau*(tau + 2.0)/gam2;
132 
133  G4double loss(0.), siga(0.);
134 
135  const G4Material* material = couple->GetMaterial();
136 
137  // Gaussian regime
138  // for heavy particles only and conditions
139  // for Gauusian fluct. has been changed
140  //
141  if ((particleMass > electron_mass_c2) &&
142  (meanLoss >= minNumberInteractionsBohr*tmax))
143  {
144  G4double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
145  (1.+m_massrate*(2.*gam+m_massrate)) ;
146  if (tmaxkine <= 2.*tmax)
147  {
148  electronDensity = material->GetElectronDensity();
149  siga = sqrt((1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
150  * electronDensity * chargeSquare);
151 
152  G4double sn = meanLoss/siga;
153 
154  // thick target case
155  if (sn >= 2.0) {
156 
157  G4double twomeanLoss = meanLoss + meanLoss;
158  do {
159  loss = G4RandGauss::shoot(rndmEngineF,meanLoss,siga);
160  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
161  } while (0.0 > loss || twomeanLoss < loss);
162 
163  // Gamma distribution
164  } else {
165 
166  G4double neff = sn*sn;
167  loss = meanLoss*G4RandGamma::shoot(rndmEngineF,neff,1.0)/neff;
168  }
169  //G4cout << "Gauss: " << loss << G4endl;
170  return loss;
171  }
172  }
173 
174  // Glandz regime : initialisation
175  //
176  if (material != lastMaterial) {
177  f1Fluct = material->GetIonisation()->GetF1fluct();
178  f2Fluct = material->GetIonisation()->GetF2fluct();
179  e1Fluct = material->GetIonisation()->GetEnergy1fluct();
180  e2Fluct = material->GetIonisation()->GetEnergy2fluct();
181  e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
182  e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
183  ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy();
184  ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
185  e0 = material->GetIonisation()->GetEnergy0fluct();
186  esmall = 0.5*sqrt(e0*ipotFluct);
187  lastMaterial = material;
188  }
189 
190  // very small step or low-density material
191  if(tmax <= e0) { return meanLoss; }
192 
193  // width correction for small cuts
194  G4double scaling = std::min(1.+0.5*CLHEP::keV/tmax,1.50);
195  meanLoss /= scaling;
196 
197  G4double a1(0.0), a2(0.0), a3(0.0);
198 
199  loss = 0.0;
200 
201  e1 = e1Fluct;
202  e2 = e2Fluct;
203 
204  if(tmax > ipotFluct) {
205  G4double w2 = G4Log(2.*electron_mass_c2*beta2*gam2)-beta2;
206 
207  if(w2 > ipotLogFluct) {
208  if(w2 > e2LogFluct) {
209  G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
210  a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
211  a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
212  } else {
213  a1 = meanLoss*(1.-rate)/e1;
214  }
215  if(a1 < a0) {
216  G4double fwnow = 0.5+(fw-0.5)*sqrt(a1/a0);
217  a1 /= fwnow;
218  e1 *= fwnow;
219  } else {
220  a1 /= fw;
221  e1 = fw*e1Fluct;
222  }
223  }
224  }
225 
226  G4double w1 = tmax/e0;
227  if(tmax > e0) {
228  a3 = rate*meanLoss*(tmax-e0)/(e0*tmax*G4Log(w1));
229  if(a1+a2 <= 0.) {
230  a3 /= rate;
231  }
232  }
233  //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
234  G4double emean = 0.;
235  G4double sig2e = 0.;
236 
237  // excitation of type 1
238  if(a1 > 0.0) { AddExcitation(rndmEngineF, a1, e1, emean, loss, sig2e); }
239 
240  // excitation of type 2
241  if(a2 > 0.0) { AddExcitation(rndmEngineF, a2, e2, emean, loss, sig2e); }
242 
243  if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
244 
245  // ionisation
246  if(a3 > 0.) {
247  emean = 0.;
248  sig2e = 0.;
249  G4double p3 = a3;
250  G4double alfa = 1.;
251  if(a3 > nmaxCont)
252  {
253  alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
254  G4double alfa1 = alfa*G4Log(alfa)/(alfa-1.);
255  G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
256  emean += namean*e0*alfa1;
257  sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
258  p3 = a3-namean;
259  }
260 
261  G4double w2 = alfa*e0;
262  if(tmax > w2) {
263  G4double w = (tmax-w2)/tmax;
264  G4int nnb = G4Poisson(p3);
265  if(nnb > 0) {
266  if(nnb > sizearray) {
267  sizearray = nnb;
268  delete [] rndmarray;
269  rndmarray = new G4double[nnb];
270  }
271  rndmEngineF->flatArray(nnb, rndmarray);
272  for (G4int k=0; k<nnb; ++k) { loss += w2/(1.-w*rndmarray[k]); }
273  }
274  }
275  if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
276  }
277 
278  loss *= scaling;
279 
280  return loss;
281 
282 }
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285 
286 
287 G4double G4UniversalFluctuation::Dispersion(
288  const G4Material* material,
289  const G4DynamicParticle* dp,
290  G4double tmax,
291  G4double length)
292 {
293  if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
294 
295  electronDensity = material->GetElectronDensity();
296 
297  G4double gam = (dp->GetKineticEnergy())*m_Inv_particleMass + 1.0;
298  G4double beta2 = 1.0 - 1.0/(gam*gam);
299 
300  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
301  * electronDensity * chargeSquare;
302 
303  return siga;
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
308 void
309 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
310  G4double q2)
311 {
312  if(part != particle) {
313  particle = part;
314  particleMass = part->GetPDGMass();
315 
316  // Derived quantities
317  if( particleMass != 0.0 ){
318  m_Inv_particleMass = 1.0 / particleMass;
319  m_massrate = electron_mass_c2 * m_Inv_particleMass ;
320  }else{
321  m_Inv_particleMass = DBL_MAX;
322  m_massrate = DBL_MAX;
323  }
324  }
325  chargeSquare = q2;
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......