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

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

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