source: trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc@ 960

Last change on this file since 960 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 35.0 KB
RevLine 
[819]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// $Id: G4EmCalculator.cc,v 1.37 2007/08/16 15:55:42 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4EmCalculator
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 28.06.2004
39//
40// Modifications:
41// 12.09.2004 Add verbosity (V.Ivanchenko)
42// 17.11.2004 Change signature of methods, add new methods (V.Ivanchenko)
43// 08.04.2005 Major optimisation of internal interfaces (V.Ivantchenko)
44// 08.05.2005 Use updated interfaces (V.Ivantchenko)
45// 23.10.2005 Fix computations for ions (V.Ivantchenko)
46// 11.01.2006 Add GetCSDARange (V.Ivantchenko)
47// 26.01.2006 Rename GetRange -> GetRangeFromRestricteDEDX (V.Ivanchenko)
48// 14.03.2006 correction in GetCrossSectionPerVolume (mma)
49// suppress GetCrossSectionPerAtom
50// elm->GetA() in ComputeCrossSectionPerAtom
51// 22.03.2006 Add ComputeElectronicDEDX and ComputeTotalDEDX (V.Ivanchenko)
52// 13.05.2006 Add Corrections for ion stopping (V.Ivanchenko)
53// 29.09.2006 Uncomment computation of smoothing factor (V.Ivanchenko)
54// 27.10.2006 Change test energy to access lowEnergy model from
55// 10 keV to 1 keV (V. Ivanchenko)
56// 15.03.2007 Add ComputeEnergyCutFromRangeCut methods (V.Ivanchenko)
57//
58// Class Description:
59//
60// -------------------------------------------------------------------
61//
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65#include "G4EmCalculator.hh"
66#include "G4LossTableManager.hh"
67#include "G4VEmProcess.hh"
68#include "G4VEnergyLossProcess.hh"
69#include "G4VMultipleScattering.hh"
70#include "G4Material.hh"
71#include "G4MaterialCutsCouple.hh"
72#include "G4ParticleDefinition.hh"
73#include "G4ParticleTable.hh"
74#include "G4PhysicsTable.hh"
75#include "G4ProductionCutsTable.hh"
76#include "G4ProcessManager.hh"
77#include "G4ionEffectiveCharge.hh"
78#include "G4RegionStore.hh"
79#include "G4Element.hh"
80#include "G4EmCorrections.hh"
81#include "G4GenericIon.hh"
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85G4EmCalculator::G4EmCalculator()
86{
87 manager = G4LossTableManager::Instance();
88 corr = manager->EmCorrections();
89 nLocalMaterials = 0;
90 verbose = 0;
91 currentCoupleIndex = 0;
92 currentCouple = 0;
93 currentMaterial = 0;
94 currentParticle = 0;
95 baseParticle = 0;
96 currentLambda = 0;
97 currentModel = 0;
98 loweModel = 0;
99 chargeSquare = 1.0;
100 massRatio = 1.0;
101 currentParticleName= "";
102 currentMaterialName= "";
103 theGenericIon = G4GenericIon::GenericIon();
104 ionEffCharge = new G4ionEffectiveCharge();
105 isIon = false;
106 isApplicable = false;
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
111G4EmCalculator::~G4EmCalculator()
112{
113 delete ionEffCharge;
114 for (G4int i=0; i<nLocalMaterials; i++) {
115 delete localCouples[i];
116 }
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p,
122 const G4Material* mat, const G4Region* region)
123{
124 G4double res = 0.0;
125 const G4MaterialCutsCouple* couple = FindCouple(mat, region);
126 if(couple && UpdateParticle(p, kinEnergy) ) {
127 res = manager->GetDEDX(p, kinEnergy, couple);
128 if(verbose>0) {
129 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
130 << " DEDX(MeV/mm)= " << res*mm/MeV
131 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
132 << " " << p->GetParticleName()
133 << " in " << mat->GetName()
134 << G4endl;
135 }
136 }
137 return res;
138}
139
140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141
142G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle,
143 const G4String& material, const G4String& reg)
144{
145 return GetDEDX(kinEnergy,FindParticle(particle),
146 FindMaterial(material),FindRegion(reg));
147}
148
149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150
151G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
152 const G4ParticleDefinition* p,
153 const G4Material* mat,
154 const G4Region* region)
155{
156 G4double res = 0.0;
157 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
158 if(couple && UpdateParticle(p, kinEnergy)) {
159 res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
160 if(verbose>0) {
161 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
162 << " range(mm)= " << res/mm
163 << " " << p->GetParticleName()
164 << " in " << mat->GetName()
165 << G4endl;
166 }
167 }
168 return res;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
173G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
174 const G4ParticleDefinition* p,
175 const G4Material* mat,
176 const G4Region* region)
177{
178 G4double res = 0.0;
179 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
180 if(couple && UpdateParticle(p, kinEnergy)) {
181 res = manager->GetCSDARange(p, kinEnergy, couple);
182 if(verbose>0) {
183 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
184 << " range(mm)= " << res/mm
185 << " " << p->GetParticleName()
186 << " in " << mat->GetName()
187 << G4endl;
188 }
189 }
190 return res;
191}
192
193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194
195G4double G4EmCalculator::GetRange(G4double kinEnergy,
196 const G4ParticleDefinition* p,
197 const G4Material* mat,
198 const G4Region* region)
199{
200 G4double res = 0.0;
201 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
202 if(couple && UpdateParticle(p, kinEnergy)) {
203 res = manager->GetRange(p, kinEnergy, couple);
204 if(verbose>0) {
205 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
206 << " range(mm)= " << res/mm
207 << " " << p->GetParticleName()
208 << " in " << mat->GetName()
209 << G4endl;
210 }
211 }
212 return res;
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
217G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
218 const G4String& particle,
219 const G4String& material,
220 const G4String& reg)
221{
222 return GetRangeFromRestricteDEDX(kinEnergy,FindParticle(particle),
223 FindMaterial(material),FindRegion(reg));
224}
225
226//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227
228G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
229 const G4String& particle,
230 const G4String& material,
231 const G4String& reg)
232{
233 return GetCSDARange(kinEnergy,FindParticle(particle),
234 FindMaterial(material),FindRegion(reg));
235}
236
237//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
238
239G4double G4EmCalculator::GetRange(G4double kinEnergy,
240 const G4String& particle,
241 const G4String& material,
242 const G4String& reg)
243{
244 return GetRange(kinEnergy,FindParticle(particle),
245 FindMaterial(material),FindRegion(reg));
246}
247
248//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249
250G4double G4EmCalculator::GetKinEnergy(G4double range,
251 const G4ParticleDefinition* p,
252 const G4Material* mat,
253 const G4Region* region)
254{
255 G4double res = 0.0;
256 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
257 if(couple && UpdateParticle(p, 1.0*GeV)) {
258 res = manager->GetEnergy(p, range, couple);
259 if(verbose>0) {
260 G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
261 << " KinE(MeV)= " << res/MeV
262 << " " << p->GetParticleName()
263 << " in " << mat->GetName()
264 << G4endl;
265 }
266 }
267 return res;
268}
269
270//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
271
272G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle,
273 const G4String& material, const G4String& reg)
274{
275 return GetKinEnergy(range,FindParticle(particle),
276 FindMaterial(material),FindRegion(reg));
277}
278
279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280
281G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
282 const G4ParticleDefinition* p,
283 const G4String& processName,
284 const G4Material* mat,
285 const G4Region* region)
286{
287 G4double res = 0.0;
288 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
289
290 if(couple && UpdateParticle(p, kinEnergy)) {
291 G4int idx = couple->GetIndex();
292 FindLambdaTable(p, processName);
293 if(currentLambda) {
294 G4bool b;
295 G4double e = kinEnergy*massRatio;
296 res = (((*currentLambda)[idx])->GetValue(e,b))*chargeSquare;
297 if(verbose>0) {
298 G4cout << "E(MeV)= " << kinEnergy/MeV
299 << " cross(cm-1)= " << res*cm
300 << " " << p->GetParticleName()
301 << " in " << mat->GetName();
302 if(verbose>1)
303 G4cout << " idx= " << idx << " e(MeV)= " << e
304 << " q2= " << chargeSquare;
305 G4cout << G4endl;
306 }
307 }
308 }
309 return res;
310}
311
312//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
313
314G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
315 const G4String& particle,
316 const G4String& processName,
317 const G4String& material,
318 const G4String& reg)
319{
320 return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
321 FindMaterial(material),FindRegion(reg));
322}
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325
326G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
327 const G4ParticleDefinition* p,
328 const G4String& processName,
329 const G4Material* mat,
330 const G4Region* region)
331{
332 G4double res = DBL_MAX;
333 G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
334 if(x > 0.0) res = 1.0/x;
335 if(verbose>1) {
336 G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
337 << " MFP(mm)= " << res/mm
338 << " " << p->GetParticleName()
339 << " in " << mat->GetName()
340 << G4endl;
341 }
342 return res;
343}
344
345//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346
347G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
348 const G4String& particle,
349 const G4String& processName,
350 const G4String& material,
351 const G4String& reg)
352{
353 return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
354 FindMaterial(material),FindRegion(reg));
355}
356
357//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358
359void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
360{
361 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
362 G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
363 if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
364}
365
366//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367
368void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
369{
370 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
371 G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
372 if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
373}
374
375//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
376
377void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
378{
379 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
380 G4cout << "### G4EmCalculator: Inverse Range Table for "
381 << p->GetParticleName() << G4endl;
382 if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
388 const G4ParticleDefinition* p,
389 const G4String& processName,
390 const G4Material* mat,
391 G4double cut)
392{
393 currentMaterial = mat;
394 currentMaterialName = mat->GetName();
395 G4double res = 0.0;
396 if(verbose > 1) {
397 G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
398 << " in " << currentMaterialName
399 << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
400 << G4endl;
401 }
402 if(UpdateParticle(p, kinEnergy)) {
403 if(FindEmModel(p, processName, kinEnergy)) {
404 G4double escaled = kinEnergy*massRatio;
405 if(baseParticle) {
406 res = currentModel->ComputeDEDXPerVolume(
407 mat, baseParticle, escaled, cut) * chargeSquare;
408 if(verbose > 1)
409 G4cout << baseParticle->GetParticleName()
410 << " Escaled(MeV)= " << escaled;
411 } else {
412 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
413 if(verbose > 1) G4cout << " no basePart E(MeV)= " << kinEnergy;
414 }
415 if(verbose > 1)
416 G4cout << " DEDX(MeV/mm)= " << res*mm/MeV
417 << " DEDX(MeV*cm^2/g)= "
418 << res*gram/(MeV*cm2*mat->GetDensity())
419 << G4endl;
420
421 if(isIon) {
422 if(currentModel->HighEnergyLimit() > 100.*MeV)
423 res += corr->HighOrderCorrections(p,mat,kinEnergy);
424 else
425 res *= corr->EffectiveChargeCorrection(p,mat,kinEnergy);
426 if(verbose > 1)
427 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
428 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
429 << G4endl;
430 }
431
432
433 // emulate boundary region for different parameterisations
434 G4double eth = currentModel->LowEnergyLimit();
435 // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
436 if(eth > 0.05*MeV && eth < 10.*MeV && escaled > eth &&
437 loweModel != currentModel && loweModel) {
438 G4double res0 = 0.0;
439 G4double res1 = 0.0;
440 if(baseParticle) {
441 res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
442 * chargeSquare;
443 res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
444 * chargeSquare;
445 } else {
446 res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
447 res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
448 }
449 if(verbose > 1)
450 G4cout << "At boundary energy(MeV)= " << eth/MeV
451 << " DEDX(MeV/mm)= " << res1*mm/MeV
452 << G4endl;
453 if(isIon) res1 += corr->HighOrderCorrections(p,mat,eth/massRatio);
454 //G4cout << "eth= " << eth << " escaled= " << escaled
455 // << " res0= " << res0 << " res1= "
456 // << res1 << " q2= " << chargeSquare << G4endl;
457 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
458 }
459
460 if(verbose > 0) {
461 G4cout << "E(MeV)= " << kinEnergy/MeV
462 << " DEDX(MeV/mm)= " << res*mm/MeV
463 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
464 << " cut(MeV)= " << cut/MeV
465 << " " << p->GetParticleName()
466 << " in " << currentMaterialName
467 << " Zi^2= " << chargeSquare
468 << G4endl;
469 }
470 }
471 }
472 return res;
473}
474
475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476
477G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
478 const G4ParticleDefinition* part,
479 const G4Material* mat,
480 G4double cut)
481{
482 currentMaterial = mat;
483 currentMaterialName = mat->GetName();
484 G4double dedx = 0.0;
485 if(UpdateParticle(part, kinEnergy)) {
486 G4LossTableManager* lManager = G4LossTableManager::Instance();
487 const std::vector<G4VEnergyLossProcess*> vel =
488 lManager->GetEnergyLossProcessVector();
489 G4int n = vel.size();
490 for(G4int i=0; i<n; i++) {
491 const G4ParticleDefinition* p = (vel[i])->Particle();
492 if((!isIon && p == part) || (isIon && p == theGenericIon))
493 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
494 }
495 }
496 return dedx;
497}
498
499//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
500
501G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy, const G4String& part,
502 const G4String& mat, G4double cut)
503{
504 return ComputeElectronicDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
505}
506
507//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
508
509G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition* part,
510 const G4Material* mat, G4double cut)
511{
512 G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
513 if(mass > 700.*MeV) dedx += ComputeNuclearDEDX(kinEnergy,part,mat);
514 return dedx;
515}
516
517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
518
519G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, const G4String& part,
520 const G4String& mat, G4double cut)
521{
522 return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
523}
524
525//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526
527G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
528 const G4String& particle,
529 const G4String& processName,
530 const G4String& material,
531 G4double cut)
532{
533 return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
534 FindMaterial(material),cut);
535}
536
537//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
538
539G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
540 const G4ParticleDefinition* p,
541 const G4Material* mat)
542{
543
544 G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false);
545
546 if(verbose > 1) {
547 G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
548 << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
549 << " NuclearDEDX(MeV*cm^2/g)= "
550 << res*gram/(MeV*cm2*mat->GetDensity())
551 << G4endl;
552 }
553 return res;
554}
555
556//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
557
558G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
559 const G4String& particle,
560 const G4String& material)
561{
562 return ComputeNuclearDEDX(kinEnergy,FindParticle(particle),
563 FindMaterial(material));
564}
565
566//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
567
568G4double G4EmCalculator::ComputeCrossSectionPerVolume(
569 G4double kinEnergy,
570 const G4ParticleDefinition* p,
571 const G4String& processName,
572 const G4Material* mat,
573 G4double cut)
574{
575 currentMaterial = mat;
576 currentMaterialName = mat->GetName();
577 G4double res = 0.0;
578 if(UpdateParticle(p, kinEnergy)) {
579 if(FindEmModel(p, processName, kinEnergy)) {
580 G4double e = kinEnergy;
581 if(baseParticle) {
582 e *= kinEnergy*massRatio;
583 res = currentModel->CrossSectionPerVolume(
584 mat, baseParticle, e, cut, e) * chargeSquare;
585 } else {
586 res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e);
587 }
588 if(verbose>0) {
589 G4cout << "E(MeV)= " << kinEnergy/MeV
590 << " cross(cm-1)= " << res*cm
591 << " " << p->GetParticleName()
592 << " in " << mat->GetName()
593 << G4endl;
594 }
595 }
596 }
597 return res;
598}
599
600//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
601
602G4double G4EmCalculator::ComputeCrossSectionPerVolume(
603 G4double kinEnergy,
604 const G4String& particle,
605 const G4String& processName,
606 const G4String& material,
607 G4double cut)
608{
609 return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),
610 processName,
611 FindMaterial(material),cut);
612}
613
614//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
615
616G4double G4EmCalculator::ComputeCrossSectionPerAtom(
617 G4double kinEnergy,
618 const G4ParticleDefinition* p,
619 const G4String& processName,
620 G4double Z, G4double A,
621 G4double cut)
622{
623 G4double res = 0.0;
624 if(UpdateParticle(p, kinEnergy)) {
625 if(FindEmModel(p, processName, kinEnergy)) {
626 G4double e = kinEnergy;
627 if(baseParticle) {
628 e *= kinEnergy*massRatio;
629 res = currentModel->ComputeCrossSectionPerAtom(
630 baseParticle, e, Z, A, cut) * chargeSquare;
631 } else {
632 res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
633 }
634 if(verbose>0) {
635 G4cout << "E(MeV)= " << kinEnergy/MeV
636 << " cross(barn)= " << res/barn
637 << " " << p->GetParticleName()
638 << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
639 << G4endl;
640 }
641 }
642 }
643 return res;
644}
645
646//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
647
648G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
649 const G4String& particle,
650 const G4String& processName,
651 const G4Element* elm,
652 G4double cut)
653{
654 return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),
655 processName,
656 elm->GetZ(),elm->GetA(),cut);
657}
658
659//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
660
661G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
662 const G4ParticleDefinition* p,
663 const G4String& processName,
664 const G4Material* mat,
665 G4double cut)
666{
667 G4double mfp = DBL_MAX;
668 G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
669 if(x > 0.0) mfp = 1.0/x;
670 if(verbose>1) {
671 G4cout << "E(MeV)= " << kinEnergy/MeV
672 << " MFP(mm)= " << mfp/mm
673 << " " << p->GetParticleName()
674 << " in " << mat->GetName()
675 << G4endl;
676 }
677 return mfp;
678}
679
680//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
681
682G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
683 const G4String& particle,
684 const G4String& processName,
685 const G4String& material,
686 G4double cut)
687{
688 return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
689 FindMaterial(material),cut);
690}
691
692//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
693
694G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
695 G4double range,
696 const G4ParticleDefinition* part,
697 const G4Material* mat)
698{
699 return G4ProductionCutsTable::GetProductionCutsTable()->
700 ConvertRangeToEnergy(part, mat, range);
701}
702
703//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704
705G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
706 G4double range,
707 const G4String& particle,
708 const G4String& material)
709{
710 return ComputeEnergyCutFromRangeCut(range,FindParticle(particle),
711 FindMaterial(material));
712}
713
714
715//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716
717G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
718 G4double kinEnergy)
719{
720 if(p != currentParticle) {
721 currentParticle = p;
722 baseParticle = 0;
723 currentParticleName = p->GetParticleName();
724 massRatio = 1.0;
725 mass = p->GetPDGMass();
726 chargeSquare = 1.0;
727 currentProcess = FindEnergyLossProcess(p);
728 currentProcessName = "";
729 if(currentProcess) currentProcessName = currentProcess->GetProcessName();
730
731 if(p->GetParticleType() == "nucleus" &&
732 currentParticleName != "deuteron" && currentParticleName != "triton") {
733 baseParticle = theGenericIon;
734 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
735 isIon = true;
736 // G4cout << p->GetParticleName()
737 // << " in " << currentMaterial->GetName()
738 // << " e= " << kinEnergy << G4endl;
739 chargeSquare =
740 ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
741 //G4cout << "q2= " << chargeSquare << G4endl;
742 } else {
743 isIon = false;
744 if(currentProcess) {
745 baseParticle = currentProcess->BaseParticle();
746
747 if(baseParticle) {
748 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
749 G4double q = baseParticle->GetPDGCharge()/eplus;
750 chargeSquare /= (q*q);
751 }
752 }
753 }
754 }
755 if(isIon) {
756 chargeSquare =
757 ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
758 if(currentProcess)
759 currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
760 // G4cout << "massR= " << massRatio << " q2= " << chargeSquare << G4endl;
761 }
762 return true;
763}
764
765//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
766
767const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
768{
769 const G4ParticleDefinition* p = 0;
770 if(name != currentParticleName) {
771 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
772 if(!p) {
773 G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
774 << name << G4endl;
775 }
776 } else {
777 p = currentParticle;
778 }
779 return p;
780}
781
782//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
783
784const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
785{
786 if(name != currentMaterialName) {
787 currentMaterial = G4Material::GetMaterial(name);
788 currentMaterialName = name;
789 if(!currentMaterial)
790 G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
791 << name << G4endl;
792 }
793 return currentMaterial;
794}
795
796//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
797
798const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
799{
800 const G4Region* r = 0;
801 if(reg != "" || reg != "world")
802 r = G4RegionStore::GetInstance()->GetRegion(reg);
803 else
804 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
805 return r;
806}
807
808//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
809
810const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
811 const G4Material* material,
812 const G4Region* region)
813{
814 if(!material) return 0;
815 currentMaterial = material;
816 currentMaterialName = material->GetName();
817 // Access to materials
818 const G4ProductionCutsTable* theCoupleTable=
819 G4ProductionCutsTable::GetProductionCutsTable();
820 const G4Region* r = region;
821 if(!r)
822 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
823
824 return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
825
826}
827
828//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
829
830G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
831{
832 if(!material) return false;
833 currentMaterial = material;
834 currentMaterialName = material->GetName();
835 for (G4int i=0; i<nLocalMaterials; i++) {
836 if(material == localMaterials[i] && cut == localCuts[i]) {
837 currentCouple = localCouples[i];
838 currentCoupleIndex = currentCouple->GetIndex();
839 currentCut = cut;
840 return true;
841 }
842 }
843 const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
844 localMaterials.push_back(material);
845 localCouples.push_back(cc);
846 localCuts.push_back(cut);
847 nLocalMaterials++;
848 currentCouple = cc;
849 currentCoupleIndex = currentCouple->GetIndex();
850 currentCut = cut;
851 return true;
852}
853
854//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
855
856void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
857 const G4String& processName)
858{
859 // Search for the process
860 if (p != currentParticle || processName != currentName) {
861 currentName = processName;
862 currentLambda = 0;
863
864 G4String partname = p->GetParticleName();
865 const G4ParticleDefinition* part = p;
866 if(isIon) part = theGenericIon;
867
868 G4LossTableManager* lManager = G4LossTableManager::Instance();
869 const std::vector<G4VEnergyLossProcess*> vel =
870 lManager->GetEnergyLossProcessVector();
871 G4int n = vel.size();
872 for(G4int i=0; i<n; i++) {
873 if((vel[i])->GetProcessName() == currentName &&
874 (vel[i])->Particle() == part)
875 {
876 currentLambda = (vel[i])->LambdaTable();
877 isApplicable = true;
878 break;
879 }
880 }
881 if(!currentLambda) {
882 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
883 G4int n = vem.size();
884 for(G4int i=0; i<n; i++) {
885 if((vem[i])->GetProcessName() == currentName &&
886 (vem[i])->Particle() == part)
887 {
888 currentLambda = (vem[i])->LambdaTable();
889 isApplicable = true;
890 break;
891 }
892 }
893 }
894 if(!currentLambda) {
895 const std::vector<G4VMultipleScattering*> vmsc =
896 lManager->GetMultipleScatteringVector();
897 G4int n = vmsc.size();
898 for(G4int i=0; i<n; i++) {
899 if((vmsc[i])->GetProcessName() == currentName &&
900 (vmsc[i])->Particle() == part)
901 {
902 currentLambda = (vmsc[i])->LambdaTable();
903 isApplicable = true;
904 break;
905 }
906 }
907 }
908 }
909}
910
911//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
912
913G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
914 const G4String& processName,
915 G4double kinEnergy)
916{
917 G4bool res = false;
918 if(!p) {
919 G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined"
920 << G4endl;
921 return res;
922 }
923 G4String partname = p->GetParticleName();
924 const G4ParticleDefinition* part = p;
925 G4double scaledEnergy = kinEnergy*massRatio;
926 if(isIon) part = theGenericIon;
927
928 if(verbose > 1) {
929 G4cout << "G4EmCalculator::FindEmModel for " << partname
930 << " (type= " << p->GetParticleType()
931 << ") and " << processName << " at e(MeV)= " << scaledEnergy;
932 if(p != part) G4cout << " GenericIon is the base particle";
933 G4cout << G4endl;
934 }
935
936 // Search for the process
937 currentName = processName;
938 currentModel = 0;
939 loweModel = 0;
940 size_t idx = 0;
941 G4LossTableManager* lManager = G4LossTableManager::Instance();
942 const std::vector<G4VEnergyLossProcess*> vel =
943 lManager->GetEnergyLossProcessVector();
944 G4int n = vel.size();
945 for(G4int i=0; i<n; i++) {
946 // G4cout << "i= " << i << " part= "
947 // << (vel[i])->Particle()->GetParticleName()
948 // << " proc= " << (vel[i])->GetProcessName() << G4endl;
949 if((vel[i])->GetProcessName() == currentName &&
950 (vel[i])->Particle() == part)
951 {
952 const G4ParticleDefinition* bp = (vel[i])->BaseParticle();
953 // G4cout << "i= " << i << " bp= " << bp << G4endl;
954 if(!bp) {
955 currentModel = (vel[i])->SelectModelForMaterial(scaledEnergy, idx);
956 loweModel = (vel[i])->SelectModelForMaterial(keV, idx);
957 isApplicable = true;
958 break;
959 } else {
960 for(G4int j=0; j<n; j++) {
961 if((vel[j])->Particle() == bp) {
962 currentModel = (vel[j])->SelectModelForMaterial(scaledEnergy, idx);
963 loweModel = (vel[j])->SelectModelForMaterial(keV, idx);
964 isApplicable = true;
965 break;
966 }
967 }
968 }
969 }
970 }
971 if(!currentModel) {
972 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
973 G4int n = vem.size();
974 for(G4int i=0; i<n; i++) {
975 if((vem[i])->GetProcessName() == currentName &&
976 (vem[i])->Particle() == part)
977 {
978 currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx);
979 loweModel = (vem[i])->SelectModelForMaterial(keV, idx);
980 isApplicable = true;
981 break;
982 }
983 }
984 }
985 if(!currentModel) {
986 const std::vector<G4VMultipleScattering*> vmsc =
987 lManager->GetMultipleScatteringVector();
988 G4int n = vmsc.size();
989 for(G4int i=0; i<n; i++) {
990 if((vmsc[i])->GetProcessName() == currentName &&
991 (vmsc[i])->Particle() == part)
992 {
993 currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx);
994 loweModel = (vmsc[i])->SelectModelForMaterial(keV, idx);
995 isApplicable = true;
996 break;
997 }
998 }
999 }
1000 if(currentModel) res = true;
1001 return res;
1002}
1003
1004//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1005
1006G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
1007 const G4ParticleDefinition* p)
1008{
1009 G4VEnergyLossProcess* elp = 0;
1010 G4String partname = p->GetParticleName();
1011 const G4ParticleDefinition* part = p;
1012 if(p->GetParticleType() == "nucleus" &&
1013 partname != "deuteron" &&
1014 partname != "triton") part = G4GenericIon::GenericIon();
1015
1016 G4LossTableManager* lManager = G4LossTableManager::Instance();
1017 const std::vector<G4VEnergyLossProcess*> vel =
1018 lManager->GetEnergyLossProcessVector();
1019 G4int n = vel.size();
1020 for(G4int i=0; i<n; i++) {
1021 if((vel[i])->Particle() == part) {
1022 elp = vel[i];
1023 break;
1024 }
1025 }
1026 return elp;
1027}
1028
1029//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1030
1031void G4EmCalculator::SetVerbose(G4int verb)
1032{
1033 verbose = verb;
1034}
1035
1036//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1037
Note: See TracBrowser for help on using the repository browser.