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

Last change on this file since 968 was 961, checked in by garnier, 17 years ago

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