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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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