| 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 // ClassName: G4EnergyLossForExtrapolator
30 //
31 // Description: This class provide calculation of energy loss, fluctuation,
32 // and msc angle
33 //
34 // Author: 09.12.04 V.Ivanchenko
35 //
36 // Modification:
37 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
38 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
39 // 21-03-06 Add verbosity defined in the constructor and Initialisation
40 // start only when first public method is called (V.Ivanchenko)
41 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
42 // 12-05-06 SEt linLossLimit=0.001 (VI)
43 //
44 //----------------------------------------------------------------------------
45 //
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
63
64 #ifdef G4MULTITHREADED
65 G4Mutex G4EnergyLossForExtrapolator::extrMutex = G4MUTEX_INITIALIZER;
66 #endif
67
68 G4TablesForExtrapolator* G4EnergyLossForExtrapolator::tables = nullptr;
69
70 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
71 : maxEnergyTransfer(DBL_MAX), verbose(verb)
72 {
73 emin = 1.*CLHEP::MeV;
74 emax = 100.*CLHEP::TeV;
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79 G4EnergyLossForExtrapolator::~G4EnergyLossForExtrapolator()
80 {
81 if(isMaster) {
82 delete tables;
83 tables = nullptr;
84 }
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89 G4double
90 G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
91 G4double stepLength,
92 const G4Material* mat,
93 const G4ParticleDefinition* part)
94 {
95 G4double kinEnergyFinal = kinEnergy;
96 if(SetupKinematics(part, mat, kinEnergy)) {
97 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
98 G4double r = ComputeRange(kinEnergy,part,mat);
99 if(r <= step) {
100 kinEnergyFinal = 0.0;
101 } else if(step < linLossLimit*r) {
102 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part,mat);
103 } else {
104 G4double r1 = r - step;
105 kinEnergyFinal = ComputeEnergy(r1,part,mat);
106 }
107 }
108 return kinEnergyFinal;
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112
113 G4double
114 G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
115 G4double stepLength,
116 const G4Material* mat,
117 const G4ParticleDefinition* part)
118 {
119 //G4cout << "G4EnergyLossForExtrapolator::EnergyBeforeStep" << G4endl;
120 G4double kinEnergyFinal = kinEnergy;
121
122 if(SetupKinematics(part, mat, kinEnergy)) {
123 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
124 G4double r = ComputeRange(kinEnergy,part,mat);
125
126 if(step < linLossLimit*r) {
127 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part,mat);
128 } else {
129 G4double r1 = r + step;
130 kinEnergyFinal = ComputeEnergy(r1,part,mat);
131 }
132 }
133 return kinEnergyFinal;
134 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137
138 G4double
139 G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
140 G4double stepLength,
141 const G4Material* mat,
142 const G4ParticleDefinition* part)
143 {
144 G4double res = stepLength;
145 //G4cout << "## G4EnergyLossForExtrapolator::TrueStepLength L= " << res
146 // << " " << part->GetParticleName() << G4endl;
147 if(SetupKinematics(part, mat, kinEnergy)) {
148 if(part == electron || part == positron) {
149 const G4double x = stepLength*
150 ComputeValue(kinEnergy, GetPhysicsTable(fMscElectron), mat->GetIndex());
151 //G4cout << " x= " << x << G4endl;
152 if(x < 0.2) { res *= (1.0 + 0.5*x + x*x/3.0); }
153 else if(x < 0.9999) { res = -G4Log(1.0 - x)*stepLength/x; }
154 else { res = ComputeRange(kinEnergy, part, mat); }
155 } else {
156 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
157 }
158 }
159 return res;
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163
164 G4bool
165 G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
166 const G4Material* mat,
167 G4double kinEnergy)
168 {
169 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
170 if(nullptr == part || nullptr == mat || kinEnergy < CLHEP::keV)
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 << "### G4EnergyLossForExtrapolator WARNING: material index i= "
181 << i << " above number of materials " << nmat << G4endl;
182 return false;
183 } else {
184 currentMaterial = mat;
185 electronDensity = mat->GetElectronDensity();
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/mass;
201 tmax = 2.0*bg2*CLHEP::electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
202 }
203 tmax = std::min(tmax, maxEnergyTransfer);
204 }
205 return true;
206 }
207
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209
210 const G4ParticleDefinition*
211 G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
212 {
213 currentParticle = G4ParticleTable::GetParticleTable()->FindParticle(name);
214 if(nullptr == currentParticle) {
215 G4cout << "### G4EnergyLossForExtrapolator WARNING: "
216 << "FindParticle() fails to find " << name << G4endl;
217 }
218 return currentParticle;
219 }
220
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222
223 G4double
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4double ekin,
225 const G4ParticleDefinition* part,
226 const G4Material* mat)
227 {
228 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
229 G4double x = 0.0;
230 if(part == electron) {
231 x = ComputeValue(ekin, GetPhysicsTable(fDedxElectron), mat->GetIndex());
232 } else if(part == positron) {
233 x = ComputeValue(ekin, GetPhysicsTable(fDedxPositron), mat->GetIndex());
234 } else if(part == muonPlus || part == muonMinus) {
235 x = ComputeValue(ekin, GetPhysicsTable(fDedxMuon), mat->GetIndex());
236 } else {
237 G4double e = ekin*CLHEP::proton_mass_c2/part->GetPDGMass();
238 G4double q = part->GetPDGCharge()/CLHEP::eplus;
239 x = ComputeValue(e, GetPhysicsTable(fDedxProton), mat->GetIndex())*q*q;
240 }
241 return x;
242 }
243
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245
246 G4double
247 G4EnergyLossForExtrapolator::ComputeRange(G4double ekin,
248 const G4ParticleDefinition* part,
249 const G4Material* mat)
250 {
251 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
252 G4double x = 0.0;
253 if(part == electron) {
254 x = ComputeValue(ekin, GetPhysicsTable(fRangeElectron), mat->GetIndex());
255 } else if(part == positron) {
256 x = ComputeValue(ekin, GetPhysicsTable(fRangePositron), mat->GetIndex());
257 } else if(part == muonPlus || part == muonMinus) {
258 x = ComputeValue(ekin, GetPhysicsTable(fRangeMuon), mat->GetIndex());
259 } else {
260 G4double massratio = CLHEP::proton_mass_c2/part->GetPDGMass();
261 G4double e = ekin*massratio;
262 G4double q = part->GetPDGCharge()/CLHEP::eplus;
263 x = ComputeValue(e, GetPhysicsTable(fRangeProton), mat->GetIndex())
264 /(q*q*massratio);
265 }
266 return x;
267 }
268
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270
271 G4double
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
273 const G4ParticleDefinition* part,
274 const G4Material* mat)
275 {
276 if(mat->GetNumberOfMaterials() != nmat) { Initialisation(); }
277 G4double x = 0.0;
278 if(part == electron) {
279 x = ComputeValue(range,GetPhysicsTable(fInvRangeElectron),mat->GetIndex());
280 } else if(part == positron) {
281 x = ComputeValue(range,GetPhysicsTable(fInvRangePositron),mat->GetIndex());
282 } else if(part == muonPlus || part == muonMinus) {
283 x = ComputeValue(range, GetPhysicsTable(fInvRangeMuon), mat->GetIndex());
284 } else {
285 G4double massratio = CLHEP::proton_mass_c2/part->GetPDGMass();
286 G4double q = part->GetPDGCharge()/CLHEP::eplus;
287 G4double r = range*massratio*q*q;
288 x = ComputeValue(r, GetPhysicsTable(fInvRangeProton), mat->GetIndex())/massratio;
289 }
290 return x;
291 }
292
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294
295 G4double
296 G4EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
297 G4double stepLength,
298 const G4Material* mat,
299 const G4ParticleDefinition* part)
300 {
301 G4double sig2 = 0.0;
302 if(SetupKinematics(part, mat, kinEnergy)) {
303 G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
304 sig2 = (1.0/beta2 - 0.5)
305 *CLHEP::twopi_mc2_rcl2*tmax*step*electronDensity*charge2;
306 }
307 return sig2;
308 }
309
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
311
312 G4double G4EnergyLossForExtrapolator::AverageScatteringAngle(
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*t)*(1.0 + 0.038*G4Log(y))
323 /(beta2*gam*part->GetPDGMass());
324 }
325 return theta;
326 }
327
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329
330 void G4EnergyLossForExtrapolator::Initialisation()
331 {
332 if(verbose>0) {
333 G4cout << "### G4EnergyLossForExtrapolator::Initialisation tables= "
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(verbose, nbins, emin, emax);
350 tables->Initialisation();
351 nmat = G4Material::GetNumberOfMaterials();
352 if(verbose > 0) {
353 G4cout << "### G4EnergyLossForExtrapolator::BuildTables for "
354 << nmat << " materials Nbins= "
355 << nbins << " Emin(MeV)= " << emin << " Emax(MeV)= " << emax
356 << G4endl;
357 }
358 #ifdef G4MULTITHREADED
359 }
360 G4MUTEXUNLOCK(&extrMutex);
361 #endif
362 }
363
364 // initialisation for the next run
365 if(isMaster && G4Material::GetNumberOfMaterials() != nmat) {
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........oooOO0OOooo........oooOO0OOooo....
378