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

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

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