Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmCalculator.cc
Go to the documentation of this file.
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$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4EmCalculator
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 28.06.2004
38 //
39 // Modifications:
40 // 12.09.2004 Add verbosity (V.Ivanchenko)
41 // 17.11.2004 Change signature of methods, add new methods (V.Ivanchenko)
42 // 08.04.2005 Major optimisation of internal interfaces (V.Ivantchenko)
43 // 08.05.2005 Use updated interfaces (V.Ivantchenko)
44 // 23.10.2005 Fix computations for ions (V.Ivantchenko)
45 // 11.01.2006 Add GetCSDARange (V.Ivantchenko)
46 // 26.01.2006 Rename GetRange -> GetRangeFromRestricteDEDX (V.Ivanchenko)
47 // 14.03.2006 correction in GetCrossSectionPerVolume (mma)
48 // suppress GetCrossSectionPerAtom
49 // elm->GetA() in ComputeCrossSectionPerAtom
50 // 22.03.2006 Add ComputeElectronicDEDX and ComputeTotalDEDX (V.Ivanchenko)
51 // 13.05.2006 Add Corrections for ion stopping (V.Ivanchenko)
52 // 29.09.2006 Uncomment computation of smoothing factor (V.Ivanchenko)
53 // 27.10.2006 Change test energy to access lowEnergy model from
54 // 10 keV to 1 keV (V. Ivanchenko)
55 // 15.03.2007 Add ComputeEnergyCutFromRangeCut methods (V.Ivanchenko)
56 // 21.04.2008 Updated computations for ions (V.Ivanchenko)
57 //
58 // Class Description:
59 //
60 // -------------------------------------------------------------------
61 //
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
65 #include "G4EmCalculator.hh"
66 #include "G4SystemOfUnits.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 #include "G4ProcessVector.hh"
84 #include "G4Gamma.hh"
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89 {
90  manager = G4LossTableManager::Instance();
91  corr = manager->EmCorrections();
92  nLocalMaterials = 0;
93  verbose = 0;
94  currentCoupleIndex = 0;
95  currentCouple = 0;
96  currentMaterial = 0;
97  currentParticle = 0;
98  lambdaParticle = 0;
99  baseParticle = 0;
100  currentLambda = 0;
101  currentModel = 0;
102  currentProcess = 0;
103  loweModel = 0;
104  chargeSquare = 1.0;
105  massRatio = 1.0;
106  mass = 0.0;
107  currentCut = 0.0;
108  currentParticleName= "";
109  currentMaterialName= "";
110  currentName = "";
111  lambdaName = "";
112  theGenericIon = G4GenericIon::GenericIon();
113  ionEffCharge = new G4ionEffectiveCharge();
114  isIon = false;
115  isApplicable = false;
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
121 {
122  delete ionEffCharge;
123  for (G4int i=0; i<nLocalMaterials; ++i) {
124  delete localCouples[i];
125  }
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
131  const G4Material* mat, const G4Region* region)
132 {
133  G4double res = 0.0;
134  const G4MaterialCutsCouple* couple = FindCouple(mat, region);
135  if(couple && UpdateParticle(p, kinEnergy) ) {
136  res = manager->GetDEDX(p, kinEnergy, couple);
137 
138  if(isIon) {
139  if(FindEmModel(p, currentProcessName, kinEnergy)) {
140  G4double length = CLHEP::nm;
141  G4double eloss = res*length;
142  //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res
143  // << " de= " << eloss << G4endl;;
144  G4double niel = 0.0;
145  dynParticle.SetKineticEnergy(kinEnergy);
146  currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
147  currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
148  res = eloss/length;
149  //G4cout << " de1= " << eloss << " res1= " << res
150  // << " " << p->GetParticleName() <<G4endl;;
151  }
152  }
153 
154  if(verbose>0) {
155  G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
156  << " DEDX(MeV/mm)= " << res*mm/MeV
157  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
158  << " " << p->GetParticleName()
159  << " in " << mat->GetName()
160  << " isIon= " << isIon
161  << G4endl;
162  }
163  }
164  return res;
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
170  const G4String& material, const G4String& reg)
171 {
172  return GetDEDX(kinEnergy,FindParticle(particle),
173  FindMaterial(material),FindRegion(reg));
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 
179  const G4ParticleDefinition* p,
180  const G4Material* mat,
181  const G4Region* region)
182 {
183  G4double res = 0.0;
184  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
185  if(couple && UpdateParticle(p, kinEnergy)) {
186  res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
187  if(verbose>0) {
188  G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
189  << " range(mm)= " << res/mm
190  << " " << p->GetParticleName()
191  << " in " << mat->GetName()
192  << G4endl;
193  }
194  }
195  return res;
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 
201  const G4ParticleDefinition* p,
202  const G4Material* mat,
203  const G4Region* region)
204 {
205  G4double res = 0.0;
206  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
207  if(couple && UpdateParticle(p, kinEnergy)) {
208  res = manager->GetCSDARange(p, kinEnergy, couple);
209  if(verbose>0) {
210  G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
211  << " range(mm)= " << res/mm
212  << " " << p->GetParticleName()
213  << " in " << mat->GetName()
214  << G4endl;
215  }
216  }
217  return res;
218 }
219 
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221 
223  const G4ParticleDefinition* p,
224  const G4Material* mat,
225  const G4Region* region)
226 {
227  G4double res = 0.0;
228  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
229  if(couple && UpdateParticle(p, kinEnergy)) {
230  res = manager->GetRange(p, kinEnergy, couple);
231  if(verbose>0) {
232  G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
233  << " range(mm)= " << res/mm
234  << " " << p->GetParticleName()
235  << " in " << mat->GetName()
236  << G4endl;
237  }
238  }
239  return res;
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243 
245  const G4String& particle,
246  const G4String& material,
247  const G4String& reg)
248 {
249  return GetRangeFromRestricteDEDX(kinEnergy,FindParticle(particle),
250  FindMaterial(material),FindRegion(reg));
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 
256  const G4String& particle,
257  const G4String& material,
258  const G4String& reg)
259 {
260  return GetCSDARange(kinEnergy,FindParticle(particle),
261  FindMaterial(material),FindRegion(reg));
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
267  const G4String& particle,
268  const G4String& material,
269  const G4String& reg)
270 {
271  return GetRange(kinEnergy,FindParticle(particle),
272  FindMaterial(material),FindRegion(reg));
273 }
274 
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
276 
278  const G4ParticleDefinition* p,
279  const G4Material* mat,
280  const G4Region* region)
281 {
282  G4double res = 0.0;
283  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
284  if(couple && UpdateParticle(p, 1.0*GeV)) {
285  res = manager->GetEnergy(p, range, couple);
286  if(verbose>0) {
287  G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
288  << " KinE(MeV)= " << res/MeV
289  << " " << p->GetParticleName()
290  << " in " << mat->GetName()
291  << G4endl;
292  }
293  }
294  return res;
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298 
300  const G4String& material, const G4String& reg)
301 {
302  return GetKinEnergy(range,FindParticle(particle),
303  FindMaterial(material),FindRegion(reg));
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
309  const G4ParticleDefinition* p,
310  const G4String& processName,
311  const G4Material* mat,
312  const G4Region* region)
313 {
314  G4double res = 0.0;
315  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
316 
317  if(couple && UpdateParticle(p, kinEnergy)) {
318  G4int idx = couple->GetIndex();
319  FindLambdaTable(p, processName, kinEnergy);
320 
321  if(currentLambda) {
322  G4double e = kinEnergy*massRatio;
323  res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
324  } else {
325  res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat,
326  kinEnergy);
327  }
328  if(verbose>0) {
329  G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
330  << " cross(cm-1)= " << res*cm
331  << " " << p->GetParticleName()
332  << " in " << mat->GetName();
333  if(verbose>1)
334  G4cout << " idx= " << idx << " Escaled((MeV)= "
335  << kinEnergy*massRatio
336  << " q2= " << chargeSquare;
337  G4cout << G4endl;
338  }
339  }
340  return res;
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344 
346  const G4String& particle,
347  const G4String& processName,
348  const G4String& material,
349  const G4String& reg)
350 {
351  return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
352  FindMaterial(material),FindRegion(reg));
353 }
354 
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
356 
358  const G4String& particle,
359  G4int Z,
361  G4double kinEnergy)
362 {
363  G4double res = 0.0;
364  const G4ParticleDefinition* p = FindParticle(particle);
365  G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
366  if(p && ad) {
367  res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy);
368  }
369  return res;
370 }
371 
372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373 
375  const G4ParticleDefinition* p,
376  const G4String& processName,
377  const G4Material* mat,
378  const G4Region* region)
379 {
380  G4double res = DBL_MAX;
381  G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
382  if(x > 0.0) { res = 1.0/x; }
383  if(verbose>1) {
384  G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
385  << " MFP(mm)= " << res/mm
386  << " " << p->GetParticleName()
387  << " in " << mat->GetName()
388  << G4endl;
389  }
390  return res;
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394 
396  const G4String& particle,
397  const G4String& processName,
398  const G4String& material,
399  const G4String& reg)
400 {
401  return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
402  FindMaterial(material),FindRegion(reg));
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
408 {
409  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
410  G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
411  if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
415 
417 {
418  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
419  G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
420  if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
421 }
422 
423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424 
426 {
427  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
428  G4cout << "### G4EmCalculator: Inverse Range Table for "
429  << p->GetParticleName() << G4endl;
430  if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
431 }
432 
433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434 
436  const G4ParticleDefinition* p,
437  const G4String& processName,
438  const G4Material* mat,
439  G4double cut)
440 {
441  currentMaterial = mat;
442  currentMaterialName = mat->GetName();
443  G4double res = 0.0;
444  if(verbose > 1) {
445  G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
446  << " in " << currentMaterialName
447  << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
448  << G4endl;
449  }
450  if(UpdateParticle(p, kinEnergy)) {
451  if(FindEmModel(p, processName, kinEnergy)) {
452  G4double escaled = kinEnergy*massRatio;
453  if(baseParticle) {
454  res = currentModel->ComputeDEDXPerVolume(
455  mat, baseParticle, escaled, cut) * chargeSquare;
456  if(verbose > 1) {
457  G4cout << baseParticle->GetParticleName()
458  << " Escaled(MeV)= " << escaled;
459  }
460  } else {
461  res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
462  if(verbose > 1) { G4cout << " no basePart E(MeV)= " << kinEnergy << " "; }
463  }
464  if(verbose > 1) {
465  G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
466  << " DEDX(MeV*cm^2/g)= "
467  << res*gram/(MeV*cm2*mat->GetDensity())
468  << G4endl;
469  }
470 
471  // emulate smoothing procedure
472  G4double eth = currentModel->LowEnergyLimit();
473  // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
474  if(loweModel) {
475  G4double res0 = 0.0;
476  G4double res1 = 0.0;
477  if(baseParticle) {
478  res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
479  * chargeSquare;
480  res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
481  * chargeSquare;
482  } else {
483  res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
484  res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
485  }
486  if(verbose > 1) {
487  G4cout << "At boundary energy(MeV)= " << eth/MeV
488  << " DEDX(MeV/mm)= " << res1*mm/MeV
489  << G4endl;
490  }
491 
492  //G4cout << "eth= " << eth << " escaled= " << escaled
493  // << " res0= " << res0 << " res1= "
494  // << res1 << " q2= " << chargeSquare << G4endl;
495 
496  if(res1 > 0.0 && escaled > 0.0) {
497  res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
498  }
499  }
500 
501  // low energy correction for ions
502  if(isIon) {
503  G4double length = CLHEP::nm;
504  const G4Region* r = 0;
505  const G4MaterialCutsCouple* couple = FindCouple(mat, r);
506  G4double eloss = res*length;
507  G4double niel = 0.0;
508  dynParticle.SetKineticEnergy(kinEnergy);
509  currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
510  currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
511  res = eloss/length;
512 
513  if(verbose > 1) {
514  G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
515  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
516  << G4endl;
517  }
518  }
519  }
520 
521  if(verbose > 0) {
522  G4cout << "Sum: E(MeV)= " << kinEnergy/MeV
523  << " DEDX(MeV/mm)= " << res*mm/MeV
524  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
525  << " cut(MeV)= " << cut/MeV
526  << " " << p->GetParticleName()
527  << " in " << currentMaterialName
528  << " Zi^2= " << chargeSquare
529  << " isIon=" << isIon
530  << G4endl;
531  }
532  }
533  return res;
534 }
535 
536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
537 
539  const G4ParticleDefinition* part,
540  const G4Material* mat,
541  G4double cut)
542 {
543  currentMaterial = mat;
544  currentMaterialName = mat->GetName();
545  G4double dedx = 0.0;
546  if(UpdateParticle(part, kinEnergy)) {
547 
549  const std::vector<G4VEnergyLossProcess*> vel =
550  lManager->GetEnergyLossProcessVector();
551  G4int n = vel.size();
552 
553  //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
554  // << " n= " << n << G4endl;
555 
556  for(G4int i=0; i<n; ++i) {
557  if(vel[i]) {
558  G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
559  if(ActiveForParticle(part, p)) {
560  //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
561  // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
562  dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
563  }
564  }
565  }
566  }
567  return dedx;
568 }
569 
570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571 
573  const G4String& mat, G4double cut)
574 {
575  return ComputeElectronicDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
576 }
577 
578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
579 
581  const G4ParticleDefinition* part,
582  const G4Material* mat,
583  G4double cut)
584 {
585  G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
586  if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
587  return dedx;
588 }
589 
590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
591 
593  const G4String& part,
594  const G4String& mat,
595  G4double cut)
596 {
597  return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
598 }
599 
600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
601 
603  const G4String& particle,
604  const G4String& processName,
605  const G4String& material,
606  G4double cut)
607 {
608  return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
609  FindMaterial(material),cut);
610 }
611 
612 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
613 
615  const G4ParticleDefinition* p,
616  const G4Material* mat)
617 {
618 
619  G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false);
620 
621  if(verbose > 1) {
622  G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
623  << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
624  << " NuclearDEDX(MeV*cm^2/g)= "
625  << res*gram/(MeV*cm2*mat->GetDensity())
626  << G4endl;
627  }
628  return res;
629 }
630 
631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
632 
634  const G4String& particle,
635  const G4String& material)
636 {
637  return ComputeNuclearDEDX(kinEnergy,FindParticle(particle),
638  FindMaterial(material));
639 }
640 
641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
642 
644  G4double kinEnergy,
645  const G4ParticleDefinition* p,
646  const G4String& processName,
647  const G4Material* mat,
648  G4double cut)
649 {
650  currentMaterial = mat;
651  currentMaterialName = mat->GetName();
652  G4double res = 0.0;
653  if(UpdateParticle(p, kinEnergy)) {
654  if(FindEmModel(p, processName, kinEnergy)) {
655  G4double e = kinEnergy;
656  if(baseParticle) {
657  e *= kinEnergy*massRatio;
658  res = currentModel->CrossSectionPerVolume(
659  mat, baseParticle, e, cut, e) * chargeSquare;
660  } else {
661  res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e);
662  }
663  if(verbose>0) {
664  G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
665  << " cross(cm-1)= " << res*cm
666  << " cut(keV)= " << cut/keV
667  << " " << p->GetParticleName()
668  << " in " << mat->GetName()
669  << G4endl;
670  }
671  }
672  }
673  return res;
674 }
675 
676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
677 
679  G4double kinEnergy,
680  const G4String& particle,
681  const G4String& processName,
682  const G4String& material,
683  G4double cut)
684 {
685  return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),
686  processName,
687  FindMaterial(material),cut);
688 }
689 
690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
691 
693  G4double kinEnergy,
694  const G4ParticleDefinition* p,
695  const G4String& processName,
696  G4double Z, G4double A,
697  G4double cut)
698 {
699  G4double res = 0.0;
700  if(UpdateParticle(p, kinEnergy)) {
701  if(FindEmModel(p, processName, kinEnergy)) {
702  G4double e = kinEnergy;
703  if(baseParticle) {
704  e *= kinEnergy*massRatio;
705  res = currentModel->ComputeCrossSectionPerAtom(
706  baseParticle, e, Z, A, cut) * chargeSquare;
707  } else {
708  res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
709  }
710  if(verbose>0) {
711  G4cout << "E(MeV)= " << kinEnergy/MeV
712  << " cross(barn)= " << res/barn
713  << " " << p->GetParticleName()
714  << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
715  << G4endl;
716  }
717  }
718  }
719  return res;
720 }
721 
722 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
723 
725  const G4String& particle,
726  const G4String& processName,
727  const G4Element* elm,
728  G4double cut)
729 {
730  return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),
731  processName,
732  elm->GetZ(),elm->GetN(),cut);
733 }
734 
735 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
736 
737 G4double
739  const G4Material* mat)
740 {
741  G4double res = 0.0;
742  const G4ParticleDefinition* gamma = G4Gamma::Gamma();
743  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
744  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
745  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
746  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
747  if(res > 0.0) { res = 1.0/res; }
748  return res;
749 }
750 
751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
752 
754  const G4String& particle,
755  G4int Z,
757  G4double kinEnergy,
758  const G4Material* mat)
759 {
760  G4double res = 0.0;
761  const G4ParticleDefinition* p = FindParticle(particle);
762  G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
763  if(p && ad) {
764  res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell,
765  kinEnergy, mat);
766  }
767  return res;
768 }
769 
770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
771 
773  const G4ParticleDefinition* p,
774  const G4String& processName,
775  const G4Material* mat,
776  G4double cut)
777 {
778  G4double mfp = DBL_MAX;
779  G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
780  if(x > 0.0) { mfp = 1.0/x; }
781  if(verbose>1) {
782  G4cout << "E(MeV)= " << kinEnergy/MeV
783  << " MFP(mm)= " << mfp/mm
784  << " " << p->GetParticleName()
785  << " in " << mat->GetName()
786  << G4endl;
787  }
788  return mfp;
789 }
790 
791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
792 
794  const G4String& particle,
795  const G4String& processName,
796  const G4String& material,
797  G4double cut)
798 {
799  return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
800  FindMaterial(material),cut);
801 }
802 
803 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
804 
806  G4double range,
807  const G4ParticleDefinition* part,
808  const G4Material* mat)
809 {
811  ConvertRangeToEnergy(part, mat, range);
812 }
813 
814 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
815 
817  G4double range,
818  const G4String& particle,
819  const G4String& material)
820 {
821  return ComputeEnergyCutFromRangeCut(range,FindParticle(particle),
822  FindMaterial(material));
823 }
824 
825 
826 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827 
828 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
829  G4double kinEnergy)
830 {
831  if(p != currentParticle) {
832 
833  // new particle
834  currentParticle = p;
835  dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
836  dynParticle.SetKineticEnergy(kinEnergy);
837  baseParticle = 0;
838  currentParticleName = p->GetParticleName();
839  massRatio = 1.0;
840  mass = p->GetPDGMass();
841  chargeSquare = 1.0;
842  currentProcess = FindEnergyLossProcess(p);
843  currentProcessName = "";
844  isIon = false;
845 
846  // ionisation process exist
847  if(currentProcess) {
848  currentProcessName = currentProcess->GetProcessName();
849  baseParticle = currentProcess->BaseParticle();
850 
851  // base particle is used
852  if(baseParticle) {
853  massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
854  G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
855  chargeSquare = q*q;
856  }
857 
858  if(p->GetParticleType() == "nucleus"
859  && currentParticleName != "deuteron"
860  && currentParticleName != "triton"
861  && currentParticleName != "alpha+"
862  && currentParticleName != "helium"
863  && currentParticleName != "hydrogen"
864  ) {
865  isIon = true;
866  massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
867  baseParticle = theGenericIon;
868  // G4cout << p->GetParticleName()
869  // << " in " << currentMaterial->GetName()
870  // << " e= " << kinEnergy << G4endl;
871  }
872  }
873  }
874 
875  // Effective charge for ions
876  if(isIon) {
877  chargeSquare =
878  corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
879  * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
880  if(currentProcess) {
881  currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
882  //G4cout << "NewP: massR= " << massRatio << " q2= " << chargeSquare << G4endl;
883  }
884  }
885  return true;
886 }
887 
888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
889 
891 {
892  const G4ParticleDefinition* p = 0;
893  if(name != currentParticleName) {
895  if(!p) {
896  G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
897  << name << G4endl;
898  }
899  } else {
900  p = currentParticle;
901  }
902  return p;
903 }
904 
905 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
906 
908 {
909  const G4ParticleDefinition* p =
911  return p;
912 }
913 
914 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
915 
917 {
918  if(name != currentMaterialName) {
919  currentMaterial = G4Material::GetMaterial(name);
920  currentMaterialName = name;
921  if(!currentMaterial)
922  G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
923  << name << G4endl;
924  }
925  return currentMaterial;
926 }
927 
928 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
929 
931 {
932  const G4Region* r = 0;
933  if(reg != "" && reg != "world") {
935  } else {
936  r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
937  }
938  return r;
939 }
940 
941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
942 
944  const G4Material* material,
945  const G4Region* region)
946 {
947  if(!material) return 0;
948  currentMaterial = material;
949  currentMaterialName = material->GetName();
950  // Access to materials
951  const G4ProductionCutsTable* theCoupleTable=
953  const G4Region* r = region;
954  if(!r) {
955  r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
956  }
957  return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
958 
959 }
960 
961 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
962 
963 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
964 {
965  if(!material) return false;
966  currentMaterial = material;
967  currentMaterialName = material->GetName();
968  for (G4int i=0; i<nLocalMaterials; ++i) {
969  if(material == localMaterials[i] && cut == localCuts[i]) {
970  currentCouple = localCouples[i];
971  currentCoupleIndex = currentCouple->GetIndex();
972  currentCut = cut;
973  return true;
974  }
975  }
976  const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
977  localMaterials.push_back(material);
978  localCouples.push_back(cc);
979  localCuts.push_back(cut);
980  nLocalMaterials++;
981  currentCouple = cc;
982  currentCoupleIndex = currentCouple->GetIndex();
983  currentCut = cut;
984  return true;
985 }
986 
987 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
988 
989 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
990  const G4String& processName,
991  G4double kinEnergy)
992 {
993  // Search for the process
994  if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
995  lambdaName = processName;
996  currentLambda = 0;
997  lambdaParticle = p;
998 
999  const G4ParticleDefinition* part = p;
1000  if(isIon) { part = theGenericIon; }
1001 
1002  // Search for energy loss process
1003  currentName = processName;
1004  currentModel = 0;
1005  loweModel = 0;
1006 
1007  G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
1008  if(elproc) {
1009  currentLambda = elproc->LambdaTable();
1010  if(currentLambda) {
1011  isApplicable = true;
1012  if(verbose>1) {
1013  G4cout << "G4VEnergyLossProcess is found out: " << currentName
1014  << G4endl;
1015  }
1016  }
1017  return;
1018  }
1019 
1020  // Search for discrete process
1021  G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1022  if(proc) {
1023  currentLambda = proc->LambdaTable();
1024  if(currentLambda) {
1025  isApplicable = true;
1026  if(verbose>1) {
1027  G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
1028  }
1029  }
1030  return;
1031  }
1032 
1033  // Search for msc process
1034  G4VMultipleScattering* msc = FindMscProcess(part, processName);
1035  if(msc) {
1036  currentModel = msc->SelectModel(kinEnergy,0);
1037  /*
1038  if(currentModel) {
1039  currentLambda = currentModel->GetCrossSectionTable();
1040  if(currentLambda) {
1041  isApplicable = true;
1042  if(verbose>1) {
1043  G4cout << "G4VMultipleScattering is found out: " << currentName
1044  << G4endl;
1045  }
1046  }
1047  }
1048  */
1049  }
1050  }
1051 }
1052 
1053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1054 
1055 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
1056  const G4String& processName,
1057  G4double kinEnergy)
1058 {
1059  isApplicable = false;
1060  if(!p) {
1061  G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined"
1062  << G4endl;
1063  return isApplicable;
1064  }
1065  G4String partname = p->GetParticleName();
1066  const G4ParticleDefinition* part = p;
1067  G4double scaledEnergy = kinEnergy*massRatio;
1068  if(isIon) { part = theGenericIon; }
1069 
1070  if(verbose > 1) {
1071  G4cout << "## G4EmCalculator::FindEmModel for " << partname
1072  << " (type= " << p->GetParticleType()
1073  << ") and " << processName << " at E(MeV)= " << scaledEnergy
1074  << G4endl;
1075  if(p != part) { G4cout << " GenericIon is the base particle" << G4endl; }
1076  }
1077 
1078  // Search for energy loss process
1079  currentName = processName;
1080  currentModel = 0;
1081  loweModel = 0;
1082  size_t idx = 0;
1083 
1084  G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
1085  if(elproc) {
1086  currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
1087  G4double eth = currentModel->LowEnergyLimit();
1088  if(eth > 0.0) {
1089  loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1090  if(loweModel == currentModel) { loweModel = 0; }
1091  }
1092  }
1093 
1094  // Search for discrete process
1095  if(!currentModel) {
1096  G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1097  if(proc) {
1098  currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
1099  G4double eth = currentModel->LowEnergyLimit();
1100  if(eth > 0.0) {
1101  loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1102  if(loweModel == currentModel) { loweModel = 0; }
1103  }
1104  }
1105  }
1106 
1107  // Search for msc process
1108  if(!currentModel) {
1109  G4VMultipleScattering* proc = FindMscProcess(part, processName);
1110  if(proc) {
1111  currentModel = proc->SelectModel(kinEnergy, idx);
1112  loweModel = 0;
1113  }
1114  }
1115  if(currentModel) {
1116  if(loweModel == currentModel) { loweModel = 0; }
1117  isApplicable = true;
1118  if(verbose > 1) {
1119  G4cout << " Model <" << currentModel->GetName()
1120  << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
1121  << " for " << part->GetParticleName();
1122  if(elproc) {
1123  G4cout << " and " << elproc->GetProcessName() << " " << elproc
1124  << G4endl;
1125  }
1126  if(loweModel) {
1127  G4cout << " LowEnergy model <" << loweModel->GetName() << ">";
1128  }
1129  G4cout << G4endl;
1130  }
1131  }
1132  return isApplicable;
1133 }
1134 
1135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1136 
1137 G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
1138  const G4ParticleDefinition* p)
1139 {
1140  G4VEnergyLossProcess* elp = 0;
1141  G4String partname = p->GetParticleName();
1142  const G4ParticleDefinition* part = p;
1143 
1144  if(p->GetParticleType() == "nucleus"
1145  && currentParticleName != "deuteron"
1146  && currentParticleName != "triton"
1147  && currentParticleName != "alpha+"
1148  && currentParticleName != "helium"
1149  && currentParticleName != "hydrogen"
1150  ) { part = theGenericIon; }
1151 
1153  return elp;
1154 }
1155 
1156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1157 
1159 G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
1160  const G4String& processName)
1161 {
1162  G4VEnergyLossProcess* proc = 0;
1163  const std::vector<G4VEnergyLossProcess*> v =
1165  G4int n = v.size();
1166  for(G4int i=0; i<n; ++i) {
1167  if((v[i])->GetProcessName() == processName) {
1168  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1169  if(ActiveForParticle(part, p)) {
1170  proc = v[i];
1171  break;
1172  }
1173  }
1174  }
1175  return proc;
1176 }
1177 
1178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1179 
1180 G4VEmProcess*
1181 G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
1182  const G4String& processName)
1183 {
1184  G4VEmProcess* proc = 0;
1185  const std::vector<G4VEmProcess*> v =
1187  G4int n = v.size();
1188  for(G4int i=0; i<n; ++i) {
1189  if((v[i])->GetProcessName() == processName) {
1190  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1191  if(ActiveForParticle(part, p)) {
1192  proc = v[i];
1193  break;
1194  }
1195  }
1196  }
1197  return proc;
1198 }
1199 
1200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1201 
1203 G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
1204  const G4String& processName)
1205 {
1206  G4VMultipleScattering* proc = 0;
1207  const std::vector<G4VMultipleScattering*> v =
1209  G4int n = v.size();
1210  for(G4int i=0; i<n; ++i) {
1211  if((v[i])->GetProcessName() == processName) {
1212  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1213  if(ActiveForParticle(part, p)) {
1214  proc = v[i];
1215  break;
1216  }
1217  }
1218  }
1219  return proc;
1220 }
1221 
1222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1223 
1224 G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
1225  G4VProcess* proc)
1226 {
1227  G4ProcessManager* pm = part->GetProcessManager();
1228  G4ProcessVector* pv = pm->GetProcessList();
1229  G4int n = pv->size();
1230  G4bool res = false;
1231  for(G4int i=0; i<n; ++i) {
1232  if((*pv)[i] == proc) {
1233  if(pm->GetProcessActivation(i)) { res = true; }
1234  break;
1235  }
1236  }
1237  return res;
1238 }
1239 
1240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1241 
1243 {
1244  verbose = verb;
1245 }
1246 
1247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1248