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

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

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