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

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

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