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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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