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

Last change on this file since 1315 was 1315, checked in by garnier, 15 years ago

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