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

Last change on this file since 1303 was 1228, checked in by garnier, 16 years ago

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