| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //--------------------------------------------
28 //
29 // ClassName: G4EnergyLossForExtrapolator
30 //
31 // Description: This class provide calculatio
32 // and msc angle
33 //
34 // Author: 09.12.04 V.Ivanchenko
35 //
36 // Modification:
37 // 08-04-05 Rename Propogator -> Extrapolator
38 // 16-03-06 Add muon tables and fix bug in uni
39 // 21-03-06 Add verbosity defined in the const
40 // start only when first public metho
41 // 03-05-06 Remove unused pointer G4Material*
42 // 12-05-06 SEt linLossLimit=0.001 (VI)
43 //
44 //--------------------------------------------
45 //
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 #include "G4EnergyLossForExtrapolator.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh"
54 #include "G4MaterialCutsCouple.hh"
55 #include "G4Electron.hh"
56 #include "G4Positron.hh"
57 #include "G4Proton.hh"
58 #include "G4MuonPlus.hh"
59 #include "G4MuonMinus.hh"
60 #include "G4ParticleTable.hh"
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 #ifdef G4MULTITHREADED
65 G4Mutex G4EnergyLossForExtrapolator::extrMutex
66 #endif
67
68 G4TablesForExtrapolator* G4EnergyLossForExtrap
69
70 G4EnergyLossForExtrapolator::G4EnergyLossForEx
71 : maxEnergyTransfer(DBL_MAX), verbose(verb)
72 {
73 emin = 1.*CLHEP::MeV;
74 emax = 100.*CLHEP::TeV;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oo
78
79 G4EnergyLossForExtrapolator::~G4EnergyLossForE
80 {
81 if(isMaster) {
82 delete tables;
83 tables = nullptr;
84 }
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4double
90 G4EnergyLossForExtrapolator::EnergyAfterStep(G
91 G4double stepLength,
92 const G4Material* mat,
93 const G4ParticleDefinition* par
94 {
95 G4double kinEnergyFinal = kinEnergy;
96 if(SetupKinematics(part, mat, kinEnergy)) {
97 G4double step = TrueStepLength(kinEnergy,s
98 G4double r = ComputeRange(kinEnergy,part,
99 if(r <= step) {
100 kinEnergyFinal = 0.0;
101 } else if(step < linLossLimit*r) {
102 kinEnergyFinal -= step*ComputeDEDX(kinEn
103 } else {
104 G4double r1 = r - step;
105 kinEnergyFinal = ComputeEnergy(r1,part,m
106 }
107 }
108 return kinEnergyFinal;
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oo
112
113 G4double
114 G4EnergyLossForExtrapolator::EnergyBeforeStep(
115 G4double stepLength,
116 const G4Material* mat,
117 const G4ParticleDefinition* pa
118 {
119 //G4cout << "G4EnergyLossForExtrapolator::En
120 G4double kinEnergyFinal = kinEnergy;
121
122 if(SetupKinematics(part, mat, kinEnergy)) {
123 G4double step = TrueStepLength(kinEnergy,s
124 G4double r = ComputeRange(kinEnergy,part,m
125
126 if(step < linLossLimit*r) {
127 kinEnergyFinal += step*ComputeDEDX(kinEn
128 } else {
129 G4double r1 = r + step;
130 kinEnergyFinal = ComputeEnergy(r1,part,m
131 }
132 }
133 return kinEnergyFinal;
134 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oo
137
138 G4double
139 G4EnergyLossForExtrapolator::TrueStepLength(G4
140 G4double stepLength,
141 const G4Material* mat,
142 const G4ParticleDefinition* part
143 {
144 G4double res = stepLength;
145 //G4cout << "## G4EnergyLossForExtrapolator:
146 // << " " << part->GetParticleName() <<
147 if(SetupKinematics(part, mat, kinEnergy)) {
148 if(part == electron || part == positron) {
149 const G4double x = stepLength*
150 ComputeValue(kinEnergy, GetPhysicsTable(fMsc
151 //G4cout << " x= " << x << G4endl;
152 if(x < 0.2) { res *= (1.0 + 0.5*
153 else if(x < 0.9999) { res = -G4Log(1.0 -
154 else { res = ComputeRange(kinEnergy, par
155 } else {
156 res = ComputeTrueStep(mat,part,kinEnergy
157 }
158 }
159 return res;
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oo
163
164 G4bool
165 G4EnergyLossForExtrapolator::SetupKinematics(c
166 const G4Material* mat,
167 G4double kinEnergy)
168 {
169 if(mat->GetNumberOfMaterials() != nmat) { In
170 if(nullptr == part || nullptr == mat || kinE
171 { return false; }
172 if(part != currentParticle) {
173 currentParticle = part;
174 G4double q = part->GetPDGCharge()/eplus;
175 charge2 = q*q;
176 }
177 if(mat != currentMaterial) {
178 size_t i = mat->GetIndex();
179 if(i >= nmat) {
180 G4cout << "### G4EnergyLossForExtrapolat
181 << i << " above number of materials " <
182 return false;
183 } else {
184 currentMaterial = mat;
185 electronDensity = mat->GetElectronDensit
186 radLength = mat->GetRadlen();
187 }
188 }
189 if(kinEnergy != kineticEnergy) {
190 kineticEnergy = kinEnergy;
191 G4double mass = part->GetPDGMass();
192 G4double tau = kinEnergy/mass;
193
194 gam = tau + 1.0;
195 bg2 = tau * (tau + 2.0);
196 beta2 = bg2/(gam*gam);
197 tmax = kinEnergy;
198 if(part == electron) tmax *= 0.5;
199 else if(part != positron) {
200 G4double r = CLHEP::electron_mass_c2/mas
201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/(
202 }
203 tmax = std::min(tmax, maxEnergyTransfer);
204 }
205 return true;
206 }
207
208 //....oooOO0OOooo........oooOO0OOooo........oo
209
210 const G4ParticleDefinition*
211 G4EnergyLossForExtrapolator::FindParticle(cons
212 {
213 currentParticle = G4ParticleTable::GetPartic
214 if(nullptr == currentParticle) {
215 G4cout << "### G4EnergyLossForExtrapolator
216 << "FindParticle() fails to find " << nam
217 }
218 return currentParticle;
219 }
220
221 //....oooOO0OOooo........oooOO0OOooo........oo
222
223 G4double
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou
225 const G4ParticleDefinition* part,
226 const
227 {
228 if(mat->GetNumberOfMaterials() != nmat) { In
229 G4double x = 0.0;
230 if(part == electron) {
231 x = ComputeValue(ekin, GetPhysicsTable(fDe
232 } else if(part == positron) {
233 x = ComputeValue(ekin, GetPhysicsTable(fDe
234 } else if(part == muonPlus || part == muonMi
235 x = ComputeValue(ekin, GetPhysicsTable(fDe
236 } else {
237 G4double e = ekin*CLHEP::proton_mass_c2/pa
238 G4double q = part->GetPDGCharge()/CLHEP::e
239 x = ComputeValue(e, GetPhysicsTable(fDedxP
240 }
241 return x;
242 }
243
244 //....oooOO0OOooo........oooOO0OOooo........oo
245
246 G4double
247 G4EnergyLossForExtrapolator::ComputeRange(G4do
248 const G4ParticleDefinition* part,
249 const G4Material* mat)
250 {
251 if(mat->GetNumberOfMaterials() != nmat) { In
252 G4double x = 0.0;
253 if(part == electron) {
254 x = ComputeValue(ekin, GetPhysicsTable(fRa
255 } else if(part == positron) {
256 x = ComputeValue(ekin, GetPhysicsTable(fRa
257 } else if(part == muonPlus || part == muonMi
258 x = ComputeValue(ekin, GetPhysicsTable(fRa
259 } else {
260 G4double massratio = CLHEP::proton_mass_c2
261 G4double e = ekin*massratio;
262 G4double q = part->GetPDGCharge()/CLHEP::e
263 x = ComputeValue(e, GetPhysicsTable(fRange
264 /(q*q*massratio);
265 }
266 return x;
267 }
268
269 //....oooOO0OOooo........oooOO0OOooo........oo
270
271 G4double
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d
273 const G4ParticleDefinition* part,
274 const G4Material* mat)
275 {
276 if(mat->GetNumberOfMaterials() != nmat) { In
277 G4double x = 0.0;
278 if(part == electron) {
279 x = ComputeValue(range,GetPhysicsTable(fIn
280 } else if(part == positron) {
281 x = ComputeValue(range,GetPhysicsTable(fIn
282 } else if(part == muonPlus || part == muonMi
283 x = ComputeValue(range, GetPhysicsTable(fI
284 } else {
285 G4double massratio = CLHEP::proton_mass_c2
286 G4double q = part->GetPDGCharge()/CLHEP::e
287 G4double r = range*massratio*q*q;
288 x = ComputeValue(r, GetPhysicsTable(fInvRa
289 }
290 return x;
291 }
292
293 //....oooOO0OOooo........oooOO0OOooo........oo
294
295 G4double
296 G4EnergyLossForExtrapolator::EnergyDispersion(
297 G4double stepLength,
298 const G4Material* mat,
299 const G4ParticleDefinition* pa
300 {
301 G4double sig2 = 0.0;
302 if(SetupKinematics(part, mat, kinEnergy)) {
303 G4double step = ComputeTrueStep(mat,part,k
304 sig2 = (1.0/beta2 - 0.5)
305 *CLHEP::twopi_mc2_rcl2*tmax*step*electro
306 }
307 return sig2;
308 }
309
310 //....oooOO0OOooo........oooOO0OOooo........oo
311
312 G4double G4EnergyLossForExtrapolator::AverageS
313 G4double kinEnergy,
314 G4double stepLength,
315 const G4Material* mat,
316 const G4ParticleDefinition* part)
317 {
318 G4double theta = 0.0;
319 if(SetupKinematics(part, mat, kinEnergy)) {
320 G4double t = stepLength/radLength;
321 G4double y = std::max(0.001, t);
322 theta = 19.23*CLHEP::MeV*std::sqrt(charge2
323 /(beta2*gam*part->GetPDGMass());
324 }
325 return theta;
326 }
327
328 //....oooOO0OOooo........oooOO0OOooo........oo
329
330 void G4EnergyLossForExtrapolator::Initialisati
331 {
332 if(verbose>0) {
333 G4cout << "### G4EnergyLossForExtrapolator
334 << tables << G4endl;
335 }
336 electron = G4Electron::Electron();
337 positron = G4Positron::Positron();
338 proton = G4Proton::Proton();
339 muonPlus = G4MuonPlus::MuonPlus();
340 muonMinus= G4MuonMinus::MuonMinus();
341
342 // initialisation for the 1st run
343 if(nullptr == tables) {
344 #ifdef G4MULTITHREADED
345 G4MUTEXLOCK(&extrMutex);
346 if(nullptr == tables) {
347 #endif
348 isMaster = true;
349 tables = new G4TablesForExtrapolator(ver
350 tables->Initialisation();
351 nmat = G4Material::GetNumberOfMaterials(
352 if(verbose > 0) {
353 G4cout << "### G4EnergyLossForExtrapol
354 << nmat << " materials Nbins= "
355 << nbins << " Emin(MeV)= " << e
356 << G4endl;
357 }
358 #ifdef G4MULTITHREADED
359 }
360 G4MUTEXUNLOCK(&extrMutex);
361 #endif
362 }
363
364 // initialisation for the next run
365 if(isMaster && G4Material::GetNumberOfMateri
366 #ifdef G4MULTITHREADED
367 G4MUTEXLOCK(&extrMutex);
368 #endif
369 tables->Initialisation();
370 #ifdef G4MULTITHREADED
371 G4MUTEXUNLOCK(&extrMutex);
372 #endif
373 }
374 nmat = G4Material::GetNumberOfMaterials();
375 }
376
377 //....oooOO0OOooo........oooOO0OOooo........oo
378