| Geant4 Cross Reference |
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 "G4Material.hh"
51 #include "G4MaterialCutsCouple.hh"
52 #include "G4DynamicParticle.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Log.hh"
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58 G4UniversalFluctuation::G4UniversalFluctuation(const G4String& nam)
59 :G4VEmFluctuationModel(nam),
60 minLoss(10.*CLHEP::eV)
61 {
62 rndmarray = new G4double[sizearray];
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66
67 G4UniversalFluctuation::~G4UniversalFluctuation()
68 {
69 delete [] rndmarray;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 void G4UniversalFluctuation::InitialiseMe(const G4ParticleDefinition* part)
75 {
76 particle = part;
77 particleMass = part->GetPDGMass();
78 const G4double q = part->GetPDGCharge()/CLHEP::eplus;
79
80 // Derived quantities
81 m_Inv_particleMass = 1.0 / particleMass;
82 m_massrate = CLHEP::electron_mass_c2 * m_Inv_particleMass;
83 chargeSquare = q*q;
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 G4double
89 G4UniversalFluctuation::SampleFluctuations(const G4MaterialCutsCouple* couple,
90 const G4DynamicParticle* dp,
91 const G4double tcut,
92 const G4double tmax,
93 const G4double length,
94 const G4double averageLoss)
95 {
96 // Calculate actual loss from the mean loss.
97 // The model used to get the fluctuations is essentially the same
98 // as in Glandz in Geant3 (Cern program library W5013, phys332).
99 // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
100
101 // shortcut for very small loss or from a step nearly equal to the range
102 // (out of validity of the model)
103 //
104 if (averageLoss < minLoss) { return averageLoss; }
105 meanLoss = averageLoss;
106 const G4double tkin = dp->GetKineticEnergy();
107 //G4cout<< "Emean= "<< meanLoss<< " tmax= "<< tmax<< " L= "<<length<<G4endl;
108
109 if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
110
111 CLHEP::HepRandomEngine* rndmEngineF = G4Random::getTheEngine();
112
113 const G4double gam = tkin * m_Inv_particleMass + 1.0;
114 const G4double gam2 = gam*gam;
115 const G4double beta = dp->GetBeta();
116 const G4double beta2 = beta*beta;
117
118 G4double loss(0.), siga(0.);
119
120 const G4Material* material = couple->GetMaterial();
121
122 // Gaussian regime
123 // for heavy particles only and conditions
124 // for Gauusian fluct. has been changed
125 //
126 if (particleMass > CLHEP::electron_mass_c2 &&
127 meanLoss >= minNumberInteractionsBohr*tcut && tmax <= 2.*tcut) {
128
129 siga = std::sqrt((tmax/beta2 - 0.5*tcut)*CLHEP::twopi_mc2_rcl2*
130 length*chargeSquare*material->GetElectronDensity());
131 const G4double sn = meanLoss/siga;
132
133 // thick target case
134 if (sn >= 2.0) {
135
136 const G4double twomeanLoss = meanLoss + meanLoss;
137 do {
138 loss = G4RandGauss::shoot(rndmEngineF, meanLoss, siga);
139 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
140 } while (0.0 > loss || twomeanLoss < loss);
141
142 // Gamma distribution
143 } else {
144
145 const G4double neff = sn*sn;
146 loss = meanLoss*G4RandGamma::shoot(rndmEngineF, neff, 1.0)/neff;
147 }
148 //G4cout << "Gauss: " << loss << G4endl;
149 return loss;
150 }
151
152 auto ioni = material->GetIonisation();
153 e0 = ioni->GetEnergy0fluct();
154
155 // very small step or low-density material
156 if(tcut <= e0) { return meanLoss; }
157
158 ipotFluct = ioni->GetMeanExcitationEnergy();
159 ipotLogFluct = ioni->GetLogMeanExcEnergy();
160
161 // width correction for small cuts
162 const G4double scaling = std::min(1.+0.5*CLHEP::keV/tcut, 1.50);
163 meanLoss /= scaling;
164
165 w2 = (tcut > ipotFluct) ?
166 G4Log(2.*CLHEP::electron_mass_c2*beta2*gam2)-beta2 : 0.0;
167 return SampleGlandz(rndmEngineF, material, tcut)*scaling;
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171
172 G4double
173 G4UniversalFluctuation::SampleGlandz(CLHEP::HepRandomEngine* rndmEngineF,
174 const G4Material*,
175 const G4double tcut)
176 {
177 G4double a1(0.0), a3(0.0);
178 G4double loss = 0.0;
179 G4double e1 = ipotFluct;
180
181 if(tcut > e1) {
182 a1 = meanLoss*(1.-rate)/e1;
183 if(a1 < a0) {
184 const G4double fwnow = 0.1+(fw-0.1)*std::sqrt(a1/a0);
185 a1 /= fwnow;
186 e1 *= fwnow;
187 } else {
188 a1 /= fw;
189 e1 *= fw;
190 }
191 }
192
193 const G4double w1 = tcut/e0;
194 a3 = rate*meanLoss*(tcut - e0)/(e0*tcut*G4Log(w1));
195 if(a1 <= 0.) { a3 /= rate; }
196
197 //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont
198 G4double emean = 0.;
199 G4double sig2e = 0.;
200
201 // excitation of type 1
202 if(a1 > 0.0) { AddExcitation(rndmEngineF, a1, e1, emean, loss, sig2e); }
203
204 if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
205
206 // ionisation
207 if(a3 > 0.) {
208 emean = 0.;
209 sig2e = 0.;
210 G4double p3 = a3;
211 G4double alfa = 1.;
212 if(a3 > nmaxCont) {
213 alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3);
214 const G4double alfa1 = alfa*G4Log(alfa)/(alfa-1.);
215 const G4double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
216 emean += namean*e0*alfa1;
217 sig2e += e0*e0*namean*(alfa-alfa1*alfa1);
218 p3 = a3 - namean;
219 }
220
221 const G4double w3 = alfa*e0;
222 if(tcut > w3) {
223 const G4double w = (tcut-w3)/tcut;
224 const G4int nnb = (G4int)G4Poisson(p3);
225 if(nnb > 0) {
226 if(nnb > sizearray) {
227 sizearray = nnb;
228 delete [] rndmarray;
229 rndmarray = new G4double[nnb];
230 }
231 rndmEngineF->flatArray(nnb, rndmarray);
232 for (G4int k=0; k<nnb; ++k) { loss += w3/(1.-w*rndmarray[k]); }
233 }
234 }
235 if(sig2e > 0.0) { SampleGauss(rndmEngineF, emean, sig2e, loss); }
236 }
237 //G4cout << "### loss=" << loss << G4endl;
238 return loss;
239 }
240
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242
243
244 G4double G4UniversalFluctuation::Dispersion(
245 const G4Material* material,
246 const G4DynamicParticle* dp,
247 const G4double tcut,
248 const G4double tmax,
249 const G4double length)
250 {
251 if(dp->GetDefinition() != particle) { InitialiseMe(dp->GetDefinition()); }
252 const G4double beta = dp->GetBeta();
253 return (tmax/(beta*beta) - 0.5*tcut) * CLHEP::twopi_mc2_rcl2 * length
254 * material->GetElectronDensity() * chargeSquare;
255 }
256
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258
259 void
260 G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
261 G4double q2)
262 {
263 if(part != particle) {
264 particle = part;
265 particleMass = part->GetPDGMass();
266
267 // Derived quantities
268 m_Inv_particleMass = 1.0 / particleMass;
269 m_massrate = CLHEP::electron_mass_c2 * m_Inv_particleMass;
270 }
271 chargeSquare = q2;
272 }
273
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275