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

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

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