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

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

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