G4EmCalculator.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4EmCalculator
00034 //
00035 // Author:        Vladimir Ivanchenko
00036 //
00037 // Creation date: 28.06.2004
00038 //
00039 // Modifications:
00040 // 12.09.2004 Add verbosity (V.Ivanchenko)
00041 // 17.11.2004 Change signature of methods, add new methods (V.Ivanchenko)
00042 // 08.04.2005 Major optimisation of internal interfaces (V.Ivantchenko)
00043 // 08.05.2005 Use updated interfaces (V.Ivantchenko)
00044 // 23.10.2005 Fix computations for ions (V.Ivantchenko)
00045 // 11.01.2006 Add GetCSDARange (V.Ivantchenko)
00046 // 26.01.2006 Rename GetRange -> GetRangeFromRestricteDEDX (V.Ivanchenko)
00047 // 14.03.2006 correction in GetCrossSectionPerVolume (mma)
00048 //            suppress GetCrossSectionPerAtom
00049 //            elm->GetA() in ComputeCrossSectionPerAtom
00050 // 22.03.2006 Add ComputeElectronicDEDX and ComputeTotalDEDX (V.Ivanchenko)
00051 // 13.05.2006 Add Corrections for ion stopping (V.Ivanchenko)
00052 // 29.09.2006 Uncomment computation of smoothing factor (V.Ivanchenko)
00053 // 27.10.2006 Change test energy to access lowEnergy model from 
00054 //            10 keV to 1 keV (V. Ivanchenko)
00055 // 15.03.2007 Add ComputeEnergyCutFromRangeCut methods (V.Ivanchenko)
00056 // 21.04.2008 Updated computations for ions (V.Ivanchenko)
00057 //
00058 // Class Description:
00059 //
00060 // -------------------------------------------------------------------
00061 //
00062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00064 
00065 #include "G4EmCalculator.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4LossTableManager.hh"
00068 #include "G4VEmProcess.hh"
00069 #include "G4VEnergyLossProcess.hh"
00070 #include "G4VMultipleScattering.hh"
00071 #include "G4Material.hh"
00072 #include "G4MaterialCutsCouple.hh"
00073 #include "G4ParticleDefinition.hh"
00074 #include "G4ParticleTable.hh"
00075 #include "G4PhysicsTable.hh"
00076 #include "G4ProductionCutsTable.hh"
00077 #include "G4ProcessManager.hh"
00078 #include "G4ionEffectiveCharge.hh"
00079 #include "G4RegionStore.hh"
00080 #include "G4Element.hh"
00081 #include "G4EmCorrections.hh"
00082 #include "G4GenericIon.hh"
00083 #include "G4ProcessVector.hh"
00084 #include "G4Gamma.hh"
00085 
00086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00087 
00088 G4EmCalculator::G4EmCalculator()
00089 {
00090   manager = G4LossTableManager::Instance();
00091   corr    = manager->EmCorrections();
00092   nLocalMaterials    = 0;
00093   verbose            = 0;
00094   currentCoupleIndex = 0;
00095   currentCouple      = 0;
00096   currentMaterial    = 0;
00097   currentParticle    = 0;
00098   lambdaParticle     = 0;
00099   baseParticle       = 0;
00100   currentLambda      = 0;
00101   currentModel       = 0;
00102   currentProcess     = 0;
00103   loweModel          = 0;
00104   chargeSquare       = 1.0;
00105   massRatio          = 1.0;
00106   mass               = 0.0;
00107   currentCut         = 0.0;
00108   currentParticleName= "";
00109   currentMaterialName= "";
00110   currentName        = "";
00111   lambdaName         = "";
00112   theGenericIon      = G4GenericIon::GenericIon();
00113   ionEffCharge       = new G4ionEffectiveCharge();
00114   isIon              = false;
00115   isApplicable       = false;
00116 }
00117 
00118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00119 
00120 G4EmCalculator::~G4EmCalculator()
00121 {
00122   delete ionEffCharge;
00123   for (G4int i=0; i<nLocalMaterials; ++i) {
00124     delete localCouples[i];
00125   }
00126 }
00127 
00128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00129 
00130 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p,
00131                                  const G4Material* mat, const G4Region* region)
00132 {
00133   G4double res = 0.0;
00134   const G4MaterialCutsCouple* couple = FindCouple(mat, region);
00135   if(couple && UpdateParticle(p, kinEnergy) ) {
00136     res = manager->GetDEDX(p, kinEnergy, couple);
00137     
00138     if(isIon) {
00139       if(FindEmModel(p, currentProcessName, kinEnergy)) {
00140         G4double length = CLHEP::nm;
00141         G4double eloss = res*length;
00142         //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res 
00143         //       << " de= " << eloss << G4endl;; 
00144         G4double niel  = 0.0;
00145         dynParticle.SetKineticEnergy(kinEnergy);
00146         currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
00147         currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
00148         res = eloss/length; 
00149         //G4cout << " de1= " << eloss << " res1= " << res 
00150         //       << " " << p->GetParticleName() <<G4endl;;
00151       }
00152     } 
00153     
00154     if(verbose>0) {
00155       G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
00156              << " DEDX(MeV/mm)= " << res*mm/MeV
00157              << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00158              << "  " <<  p->GetParticleName()
00159              << " in " <<  mat->GetName()
00160              << " isIon= " << isIon
00161              << G4endl;
00162     }
00163   }
00164   return res;
00165 }
00166 
00167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00168 
00169 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle,
00170                                  const G4String& material, const G4String& reg)
00171 {
00172   return GetDEDX(kinEnergy,FindParticle(particle),
00173                  FindMaterial(material),FindRegion(reg));
00174 }
00175 
00176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00177 
00178 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy, 
00179                                                    const G4ParticleDefinition* p,
00180                                                    const G4Material* mat,
00181                                                    const G4Region* region)
00182 {
00183   G4double res = 0.0;
00184   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00185   if(couple && UpdateParticle(p, kinEnergy)) {
00186     res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
00187     if(verbose>0) {
00188       G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00189              << " range(mm)= " << res/mm
00190              << "  " <<  p->GetParticleName()
00191              << " in " <<  mat->GetName()
00192              << G4endl;
00193     }
00194   }
00195   return res;
00196 }
00197 
00198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00199 
00200 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy, 
00201                                       const G4ParticleDefinition* p,
00202                                       const G4Material* mat, 
00203                                       const G4Region* region)
00204 {
00205   G4double res = 0.0;
00206   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00207   if(couple && UpdateParticle(p, kinEnergy)) {
00208     res = manager->GetCSDARange(p, kinEnergy, couple);
00209     if(verbose>0) {
00210       G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00211              << " range(mm)= " << res/mm
00212              << "  " <<  p->GetParticleName()
00213              << " in " <<  mat->GetName()
00214              << G4endl;
00215     }
00216   }
00217   return res;
00218 }
00219 
00220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00221 
00222 G4double G4EmCalculator::GetRange(G4double kinEnergy, 
00223                                   const G4ParticleDefinition* p,
00224                                   const G4Material* mat, 
00225                                   const G4Region* region)
00226 {
00227   G4double res = 0.0;
00228   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00229   if(couple && UpdateParticle(p, kinEnergy)) {
00230     res = manager->GetRange(p, kinEnergy, couple);
00231     if(verbose>0) {
00232       G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
00233              << " range(mm)= " << res/mm
00234              << "  " <<  p->GetParticleName()
00235              << " in " <<  mat->GetName()
00236              << G4endl;
00237     }
00238   }
00239   return res;
00240 }
00241 
00242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00243 
00244 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy, 
00245                                                    const G4String& particle,
00246                                                    const G4String& material, 
00247                                                    const G4String& reg)
00248 {
00249   return GetRangeFromRestricteDEDX(kinEnergy,FindParticle(particle),
00250                                    FindMaterial(material),FindRegion(reg));
00251 }
00252 
00253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00254 
00255 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy, 
00256                                       const G4String& particle,
00257                                       const G4String& material, 
00258                                       const G4String& reg)
00259 {
00260   return GetCSDARange(kinEnergy,FindParticle(particle),
00261                   FindMaterial(material),FindRegion(reg));
00262 }
00263 
00264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00265 
00266 G4double G4EmCalculator::GetRange(G4double kinEnergy, 
00267                                   const G4String& particle,
00268                                   const G4String& material, 
00269                                   const G4String& reg)
00270 {
00271   return GetRange(kinEnergy,FindParticle(particle),
00272                   FindMaterial(material),FindRegion(reg));
00273 }
00274 
00275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00276 
00277 G4double G4EmCalculator::GetKinEnergy(G4double range, 
00278                                       const G4ParticleDefinition* p,
00279                                       const G4Material* mat,
00280                                       const G4Region* region)
00281 {
00282   G4double res = 0.0;
00283   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00284   if(couple && UpdateParticle(p, 1.0*GeV)) {
00285     res = manager->GetEnergy(p, range, couple);
00286     if(verbose>0) {
00287       G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
00288              << " KinE(MeV)= " << res/MeV
00289              << "  " <<  p->GetParticleName()
00290              << " in " <<  mat->GetName()
00291              << G4endl;
00292     }
00293   }
00294   return res;
00295 }
00296 
00297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00298 
00299 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle,
00300                                       const G4String& material, const G4String& reg)
00301 {
00302   return GetKinEnergy(range,FindParticle(particle),
00303                       FindMaterial(material),FindRegion(reg));
00304 }
00305 
00306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00307 
00308 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
00309                                             const G4ParticleDefinition* p,
00310                                             const G4String& processName,
00311                                             const G4Material* mat,
00312                                             const G4Region* region)
00313 {
00314   G4double res = 0.0;
00315   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
00316 
00317   if(couple && UpdateParticle(p, kinEnergy)) {
00318     G4int idx = couple->GetIndex();
00319     FindLambdaTable(p, processName, kinEnergy);
00320 
00321     if(currentLambda) {
00322       G4double e = kinEnergy*massRatio;
00323       res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
00324     } else {
00325       res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, 
00326                                          kinEnergy);
00327     }
00328     if(verbose>0) {
00329       G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
00330              << " cross(cm-1)= " << res*cm
00331              << "  " <<  p->GetParticleName()
00332              << " in " <<  mat->GetName();
00333       if(verbose>1) 
00334         G4cout << "  idx= " << idx << "  Escaled((MeV)= " 
00335                << kinEnergy*massRatio 
00336                << "  q2= " << chargeSquare; 
00337       G4cout << G4endl;
00338     } 
00339   }
00340   return res;
00341 }
00342 
00343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00344 
00345 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
00346                                             const G4String& particle,
00347                                             const G4String& processName,
00348                                             const G4String& material,
00349                                             const G4String& reg)
00350 {
00351   return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
00352                                   FindMaterial(material),FindRegion(reg));
00353 }
00354 
00355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00356 
00357 G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom(
00358                                          const G4String& particle, 
00359                                          G4int Z, 
00360                                          G4AtomicShellEnumerator shell,
00361                                          G4double kinEnergy)
00362 {
00363   G4double res = 0.0;
00364   const G4ParticleDefinition* p = FindParticle(particle);
00365   G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
00366   if(p && ad) { 
00367     res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy); 
00368   }
00369   return res;
00370 }
00371 
00372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00373 
00374 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
00375                                          const G4ParticleDefinition* p,
00376                                          const G4String& processName,
00377                                          const G4Material* mat,
00378                                          const G4Region* region)
00379 {
00380   G4double res = DBL_MAX;
00381   G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
00382   if(x > 0.0) { res = 1.0/x; }
00383   if(verbose>1) {
00384     G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
00385            << " MFP(mm)= " << res/mm
00386            << "  " <<  p->GetParticleName()
00387            << " in " <<  mat->GetName()
00388            << G4endl;
00389   }
00390   return res;
00391 }
00392 
00393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00394 
00395 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
00396                                          const G4String& particle,
00397                                          const G4String& processName,
00398                                          const G4String& material,
00399                                          const G4String& reg)
00400 {
00401   return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
00402                          FindMaterial(material),FindRegion(reg));
00403 }
00404 
00405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00406 
00407 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
00408 {
00409   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00410   G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
00411   if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
00412 }
00413 
00414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00415 
00416 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
00417 {
00418   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00419   G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
00420   if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
00421 }
00422 
00423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00424 
00425 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
00426 {
00427   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
00428   G4cout << "### G4EmCalculator: Inverse Range Table for " 
00429          << p->GetParticleName() << G4endl;
00430   if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
00431 }
00432 
00433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00434 
00435 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
00436                                      const G4ParticleDefinition* p,
00437                                      const G4String& processName,
00438                                      const G4Material* mat,
00439                                            G4double cut)
00440 {
00441   currentMaterial = mat;
00442   currentMaterialName = mat->GetName();
00443   G4double res = 0.0;
00444   if(verbose > 1) {
00445     G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
00446            << " in " << currentMaterialName
00447            << " e(MeV)= " << kinEnergy/MeV << "  cut(MeV)= " << cut/MeV
00448            << G4endl;
00449   }
00450   if(UpdateParticle(p, kinEnergy)) {
00451     if(FindEmModel(p, processName, kinEnergy)) {
00452       G4double escaled = kinEnergy*massRatio;
00453       if(baseParticle) {
00454         res = currentModel->ComputeDEDXPerVolume(
00455               mat, baseParticle, escaled, cut) * chargeSquare;
00456         if(verbose > 1) {
00457           G4cout <<  baseParticle->GetParticleName()
00458                  << " Escaled(MeV)= " << escaled;
00459         }
00460       } else {
00461         res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
00462         if(verbose > 1) { G4cout <<  " no basePart E(MeV)= " << kinEnergy << " "; }
00463       }
00464       if(verbose > 1) {
00465         G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
00466                << " DEDX(MeV*cm^2/g)= "
00467                << res*gram/(MeV*cm2*mat->GetDensity())
00468                << G4endl;
00469       }
00470 
00471       // emulate smoothing procedure
00472       G4double eth = currentModel->LowEnergyLimit();
00473       // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
00474       if(loweModel) {
00475         G4double res0 = 0.0;
00476         G4double res1 = 0.0;
00477         if(baseParticle) {
00478           res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
00479                * chargeSquare;
00480           res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
00481                * chargeSquare;
00482         } else {
00483           res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
00484           res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
00485         }
00486         if(verbose > 1) {
00487           G4cout << "At boundary energy(MeV)= " << eth/MeV
00488                  << " DEDX(MeV/mm)= " << res1*mm/MeV
00489                  << G4endl;
00490         }
00491         
00492         //G4cout << "eth= " << eth << " escaled= " << escaled
00493         //       << " res0= " << res0 << " res1= "
00494         //       << res1 <<  "  q2= " << chargeSquare << G4endl;
00495         
00496         if(res1 > 0.0 && escaled > 0.0) {
00497           res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
00498         }
00499       } 
00500 
00501       // low energy correction for ions
00502       if(isIon) {
00503         G4double length = CLHEP::nm;
00504         const G4Region* r = 0;
00505         const G4MaterialCutsCouple* couple = FindCouple(mat, r);
00506         G4double eloss = res*length;
00507         G4double niel  = 0.0;
00508         dynParticle.SetKineticEnergy(kinEnergy);
00509         currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
00510         currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
00511         res = eloss/length; 
00512         
00513         if(verbose > 1) {
00514           G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
00515                  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00516                  << G4endl;
00517         }
00518       }
00519     }
00520       
00521     if(verbose > 0) {
00522       G4cout << "Sum: E(MeV)= " << kinEnergy/MeV
00523              << " DEDX(MeV/mm)= " << res*mm/MeV
00524              << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
00525              << " cut(MeV)= " << cut/MeV
00526              << "  " <<  p->GetParticleName()
00527              << " in " <<  currentMaterialName
00528              << " Zi^2= " << chargeSquare
00529              << " isIon=" << isIon
00530              << G4endl;
00531     }
00532   }
00533   return res;
00534 }
00535 
00536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00537 
00538 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
00539                                                const G4ParticleDefinition* part,
00540                                                const G4Material* mat,
00541                                                G4double cut)
00542 {
00543   currentMaterial = mat;
00544   currentMaterialName = mat->GetName();
00545   G4double dedx = 0.0;
00546   if(UpdateParticle(part, kinEnergy)) {
00547 
00548     G4LossTableManager* lManager = G4LossTableManager::Instance();
00549     const std::vector<G4VEnergyLossProcess*> vel =
00550       lManager->GetEnergyLossProcessVector();
00551     G4int n = vel.size();
00552 
00553     //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName() 
00554     //     << " n= " << n << G4endl;
00555  
00556     for(G4int i=0; i<n; ++i) {
00557       if(vel[i]) {
00558         G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
00559         if(ActiveForParticle(part, p)) {
00560           //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
00561           //     << "  " << (vel[i])->Particle()->GetParticleName() << G4endl; 
00562           dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
00563         }
00564       }
00565     }
00566   }
00567   return dedx;
00568 }
00569 
00570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00571 
00572 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy, const G4String& part,
00573                                                const G4String& mat, G4double cut)
00574 {
00575   return ComputeElectronicDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
00576 }
00577 
00578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00579 
00580 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 
00581                                           const G4ParticleDefinition* part,
00582                                           const G4Material* mat, 
00583                                           G4double cut)
00584 {
00585   G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
00586   if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
00587   return dedx;
00588 }
00589 
00590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00591 
00592 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 
00593                                           const G4String& part,
00594                                           const G4String& mat, 
00595                                           G4double cut)
00596 {
00597   return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
00598 }
00599 
00600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00601 
00602 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
00603                                      const G4String& particle,
00604                                      const G4String& processName,
00605                                      const G4String& material,
00606                                            G4double cut)
00607 {
00608   return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
00609                      FindMaterial(material),cut);
00610 }
00611 
00612 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00613 
00614 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
00615                                       const G4ParticleDefinition* p,
00616                                       const G4Material* mat)
00617 {
00618 
00619   G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false);
00620 
00621   if(verbose > 1) {
00622     G4cout <<  p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
00623            << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
00624            << " NuclearDEDX(MeV*cm^2/g)= "
00625            << res*gram/(MeV*cm2*mat->GetDensity())
00626            << G4endl;
00627   }
00628   return res;
00629 }
00630 
00631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00632 
00633 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
00634                                       const G4String& particle,
00635                                       const G4String& material)
00636 {
00637   return ComputeNuclearDEDX(kinEnergy,FindParticle(particle),
00638                             FindMaterial(material));
00639 }
00640 
00641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00642 
00643 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
00644                                                    G4double kinEnergy,
00645                                              const G4ParticleDefinition* p,
00646                                              const G4String& processName,
00647                                              const G4Material* mat,
00648                                                    G4double cut)
00649 {
00650   currentMaterial = mat;
00651   currentMaterialName = mat->GetName();
00652   G4double res = 0.0;
00653   if(UpdateParticle(p, kinEnergy)) {
00654     if(FindEmModel(p, processName, kinEnergy)) {
00655       G4double e = kinEnergy;
00656       if(baseParticle) {
00657         e *= kinEnergy*massRatio;
00658         res = currentModel->CrossSectionPerVolume(
00659               mat, baseParticle, e, cut, e) * chargeSquare;
00660       } else {
00661         res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e);
00662       }
00663       if(verbose>0) {
00664         G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
00665                << " cross(cm-1)= " << res*cm
00666                << " cut(keV)= " << cut/keV
00667                << "  " <<  p->GetParticleName()
00668                << " in " <<  mat->GetName()
00669                << G4endl;
00670       }
00671     }
00672   }
00673   return res;
00674 }
00675 
00676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00677 
00678 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
00679                                                    G4double kinEnergy,
00680                                              const G4String& particle,
00681                                              const G4String& processName,
00682                                              const G4String& material,
00683                                                    G4double cut)
00684 {
00685   return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),
00686                                       processName,
00687                                       FindMaterial(material),cut);
00688 }
00689 
00690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00691 
00692 G4double G4EmCalculator::ComputeCrossSectionPerAtom(
00693                                                    G4double kinEnergy,
00694                                              const G4ParticleDefinition* p,
00695                                              const G4String& processName,
00696                                                    G4double Z, G4double A,
00697                                                    G4double cut)
00698 {
00699   G4double res = 0.0;
00700   if(UpdateParticle(p, kinEnergy)) {
00701     if(FindEmModel(p, processName, kinEnergy)) {
00702       G4double e = kinEnergy;
00703       if(baseParticle) {
00704         e *= kinEnergy*massRatio;
00705         res = currentModel->ComputeCrossSectionPerAtom(
00706               baseParticle, e, Z, A, cut) * chargeSquare;
00707       } else {
00708         res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
00709       }
00710       if(verbose>0) {
00711         G4cout << "E(MeV)= " << kinEnergy/MeV
00712                << " cross(barn)= " << res/barn
00713                << "  " <<  p->GetParticleName()
00714                << " Z= " <<  Z << " A= " << A/(g/mole) << " g/mole"
00715                << G4endl;
00716       }
00717     }
00718   }
00719   return res;
00720 }
00721 
00722 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00723 
00724 G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
00725                                               const G4String& particle,
00726                                               const G4String& processName,
00727                                               const G4Element* elm,
00728                                                     G4double cut)
00729 {
00730   return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),
00731                                     processName,
00732                                     elm->GetZ(),elm->GetN(),cut);
00733 }
00734 
00735 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00736 
00737 G4double 
00738 G4EmCalculator::ComputeGammaAttenuationLength(G4double kinEnergy, 
00739                                               const G4Material* mat)
00740 {
00741   G4double res = 0.0;
00742   const G4ParticleDefinition* gamma = G4Gamma::Gamma();
00743   res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
00744   res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
00745   res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
00746   res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
00747   if(res > 0.0) { res = 1.0/res; }
00748   return res;
00749 }
00750 
00751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00752 
00753 G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom(
00754                                          const G4String& particle, 
00755                                          G4int Z, 
00756                                          G4AtomicShellEnumerator shell,
00757                                          G4double kinEnergy,
00758                                          const G4Material* mat)
00759 {
00760   G4double res = 0.0;
00761   const G4ParticleDefinition* p = FindParticle(particle);
00762   G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
00763   if(p && ad) { 
00764     res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell, 
00765                                                         kinEnergy, mat); 
00766   }
00767   return res;
00768 }
00769 
00770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00771 
00772 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
00773                                              const G4ParticleDefinition* p,
00774                                              const G4String& processName,
00775                                              const G4Material* mat,
00776                                                    G4double cut)
00777 {
00778   G4double mfp = DBL_MAX;
00779   G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
00780   if(x > 0.0) { mfp = 1.0/x; }
00781   if(verbose>1) {
00782     G4cout << "E(MeV)= " << kinEnergy/MeV
00783            << " MFP(mm)= " << mfp/mm
00784            << "  " <<  p->GetParticleName()
00785            << " in " <<  mat->GetName()
00786            << G4endl;
00787   }
00788   return mfp;
00789 }
00790 
00791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00792 
00793 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
00794                                              const G4String& particle,
00795                                              const G4String& processName,
00796                                              const G4String& material,
00797                                                    G4double cut)
00798 {
00799   return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
00800                              FindMaterial(material),cut);
00801 }
00802 
00803 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00804 
00805 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
00806                          G4double range, 
00807                          const G4ParticleDefinition* part,
00808                          const G4Material* mat)
00809 {
00810   return G4ProductionCutsTable::GetProductionCutsTable()->
00811     ConvertRangeToEnergy(part, mat, range);
00812 }
00813 
00814 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00815 
00816 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
00817                          G4double range, 
00818                          const G4String& particle,
00819                          const G4String& material)
00820 {
00821   return ComputeEnergyCutFromRangeCut(range,FindParticle(particle),
00822                                       FindMaterial(material));
00823 }
00824 
00825 
00826 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00827 
00828 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
00829                                       G4double kinEnergy)
00830 {
00831   if(p != currentParticle) {
00832 
00833     // new particle
00834     currentParticle = p;
00835     dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
00836     dynParticle.SetKineticEnergy(kinEnergy);
00837     baseParticle    = 0;
00838     currentParticleName = p->GetParticleName();
00839     massRatio       = 1.0;
00840     mass            = p->GetPDGMass();
00841     chargeSquare    = 1.0;
00842     currentProcess  = FindEnergyLossProcess(p);
00843     currentProcessName = "";
00844     isIon = false;
00845 
00846     // ionisation process exist
00847     if(currentProcess) {
00848       currentProcessName = currentProcess->GetProcessName();
00849       baseParticle = currentProcess->BaseParticle();
00850 
00851       // base particle is used
00852       if(baseParticle) {
00853         massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
00854         G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
00855         chargeSquare = q*q;
00856       } 
00857 
00858       if(p->GetParticleType()   == "nucleus" 
00859          && currentParticleName != "deuteron"  
00860          && currentParticleName != "triton"
00861          && currentParticleName != "alpha+"
00862          && currentParticleName != "helium"
00863          && currentParticleName != "hydrogen"
00864          ) {
00865         isIon = true;
00866         massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
00867         baseParticle = theGenericIon;
00868         //      G4cout << p->GetParticleName()
00869         // << " in " << currentMaterial->GetName()
00870         //       << "  e= " << kinEnergy << G4endl;
00871       }
00872     }
00873   }
00874 
00875   // Effective charge for ions
00876   if(isIon) {
00877     chargeSquare =
00878       corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
00879       * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
00880     if(currentProcess) {
00881       currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
00882       //G4cout << "NewP: massR= " << massRatio << "   q2= " << chargeSquare << G4endl;
00883     }
00884   }
00885   return true;
00886 }
00887 
00888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00889 
00890 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
00891 {
00892   const G4ParticleDefinition* p = 0;
00893   if(name != currentParticleName) {
00894     p = G4ParticleTable::GetParticleTable()->FindParticle(name);
00895     if(!p) {
00896       G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find " 
00897              << name << G4endl;
00898     }
00899   } else {
00900     p = currentParticle;
00901   }
00902   return p;
00903 }
00904 
00905 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00906 
00907 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
00908 {
00909   const G4ParticleDefinition* p = 
00910     G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
00911   return p;
00912 }
00913 
00914 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00915 
00916 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
00917 {
00918   if(name != currentMaterialName) {
00919     currentMaterial = G4Material::GetMaterial(name);
00920     currentMaterialName = name;
00921     if(!currentMaterial)
00922       G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find " 
00923              << name << G4endl;
00924   }
00925   return currentMaterial;
00926 }
00927 
00928 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00929 
00930 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
00931 {
00932   const G4Region* r = 0;
00933   if(reg != "" && reg != "world") {
00934     r = G4RegionStore::GetInstance()->GetRegion(reg);
00935   } else {
00936     r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
00937   }
00938   return r;
00939 }
00940 
00941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00942 
00943 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
00944                             const G4Material* material,
00945                             const G4Region* region)
00946 {
00947   if(!material) return 0;
00948   currentMaterial = material;
00949   currentMaterialName = material->GetName();
00950   // Access to materials
00951   const G4ProductionCutsTable* theCoupleTable=
00952         G4ProductionCutsTable::GetProductionCutsTable();
00953   const G4Region* r = region;
00954   if(!r) {
00955     r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
00956   }
00957   return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
00958 
00959 }
00960 
00961 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00962 
00963 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
00964 {
00965   if(!material) return false;
00966   currentMaterial = material;
00967   currentMaterialName = material->GetName();
00968   for (G4int i=0; i<nLocalMaterials; ++i) {
00969     if(material == localMaterials[i] && cut == localCuts[i]) {
00970       currentCouple = localCouples[i];
00971       currentCoupleIndex = currentCouple->GetIndex();
00972       currentCut = cut;
00973       return true;
00974     }
00975   }
00976   const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
00977   localMaterials.push_back(material);
00978   localCouples.push_back(cc);
00979   localCuts.push_back(cut);
00980   nLocalMaterials++;
00981   currentCouple = cc;
00982   currentCoupleIndex = currentCouple->GetIndex();
00983   currentCut = cut;
00984   return true;
00985 }
00986 
00987 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00988 
00989 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
00990                                      const G4String& processName,
00991                                      G4double kinEnergy)
00992 {
00993   // Search for the process
00994   if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
00995     lambdaName     = processName;
00996     currentLambda  = 0;
00997     lambdaParticle = p;
00998 
00999     const G4ParticleDefinition* part = p;
01000     if(isIon) { part = theGenericIon; }
01001 
01002     // Search for energy loss process
01003     currentName = processName;
01004     currentModel = 0;
01005     loweModel = 0;
01006 
01007     G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
01008     if(elproc) {
01009       currentLambda = elproc->LambdaTable();
01010       if(currentLambda) {
01011         isApplicable = true;
01012         if(verbose>1) { 
01013           G4cout << "G4VEnergyLossProcess is found out: " << currentName 
01014                  << G4endl;
01015         }
01016       }
01017       return;
01018     }
01019 
01020     // Search for discrete process 
01021     G4VEmProcess* proc = FindDiscreteProcess(part, processName);
01022     if(proc) {
01023       currentLambda = proc->LambdaTable();
01024       if(currentLambda) {
01025         isApplicable    = true;
01026         if(verbose>1) { 
01027           G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
01028         }
01029       }
01030       return;
01031     }
01032 
01033     // Search for msc process
01034     G4VMultipleScattering* msc = FindMscProcess(part, processName);
01035     if(msc) {
01036       currentModel = msc->SelectModel(kinEnergy,0);
01037       /*
01038       if(currentModel) {
01039         currentLambda = currentModel->GetCrossSectionTable();
01040         if(currentLambda) {
01041           isApplicable    = true;
01042           if(verbose>1) { 
01043             G4cout << "G4VMultipleScattering is found out: " << currentName 
01044                    << G4endl;
01045           }
01046         }
01047       }
01048       */
01049     }
01050   }
01051 }
01052 
01053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01054 
01055 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
01056                                    const G4String& processName,
01057                                          G4double kinEnergy)
01058 {
01059   isApplicable = false;
01060   if(!p) {
01061     G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined" 
01062            << G4endl;
01063     return isApplicable;
01064   }
01065   G4String partname =  p->GetParticleName();
01066   const G4ParticleDefinition* part = p;
01067   G4double scaledEnergy = kinEnergy*massRatio;
01068   if(isIon) { part = theGenericIon; } 
01069 
01070   if(verbose > 1) {
01071     G4cout << "## G4EmCalculator::FindEmModel for " << partname
01072            << " (type= " << p->GetParticleType()
01073            << ") and " << processName << " at E(MeV)= " << scaledEnergy 
01074            << G4endl;
01075     if(p != part) { G4cout << "  GenericIon is the base particle" << G4endl; }
01076   }
01077 
01078   // Search for energy loss process
01079   currentName = processName;
01080   currentModel = 0;
01081   loweModel = 0;
01082   size_t idx   = 0;
01083 
01084   G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
01085   if(elproc) {
01086     currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
01087     G4double eth = currentModel->LowEnergyLimit();
01088     if(eth > 0.0) {
01089       loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
01090       if(loweModel == currentModel) { loweModel = 0; }
01091     }
01092   }
01093 
01094   // Search for discrete process 
01095   if(!currentModel) {
01096     G4VEmProcess* proc = FindDiscreteProcess(part, processName);
01097     if(proc) {
01098       currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
01099       G4double eth = currentModel->LowEnergyLimit();
01100       if(eth > 0.0) {
01101         loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
01102         if(loweModel == currentModel) { loweModel = 0; }
01103       }
01104     }
01105   }
01106 
01107   // Search for msc process
01108   if(!currentModel) {
01109     G4VMultipleScattering* proc = FindMscProcess(part, processName);
01110     if(proc) {
01111       currentModel = proc->SelectModel(kinEnergy, idx);
01112       loweModel = 0;
01113     }
01114   }
01115   if(currentModel) {
01116     if(loweModel == currentModel) { loweModel = 0; }
01117     isApplicable = true;
01118     if(verbose > 1) {
01119       G4cout << "   Model <" << currentModel->GetName() 
01120              << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
01121              << " for " << part->GetParticleName();
01122       if(elproc) { 
01123         G4cout << " and " << elproc->GetProcessName() << "  " << elproc 
01124                << G4endl;
01125       }
01126       if(loweModel) { 
01127         G4cout << " LowEnergy model <" << loweModel->GetName() << ">"; 
01128       }
01129       G4cout << G4endl;
01130     } 
01131   }
01132   return isApplicable;
01133 }
01134 
01135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01136 
01137 G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
01138                       const G4ParticleDefinition* p)
01139 {
01140   G4VEnergyLossProcess* elp = 0;
01141   G4String partname =  p->GetParticleName();
01142   const G4ParticleDefinition* part = p;
01143   
01144   if(p->GetParticleType() == "nucleus" 
01145      && currentParticleName != "deuteron"  
01146      && currentParticleName != "triton"
01147      && currentParticleName != "alpha+"
01148      && currentParticleName != "helium"
01149      && currentParticleName != "hydrogen"
01150      ) { part = theGenericIon; } 
01151 
01152   elp = G4LossTableManager::Instance()->GetEnergyLossProcess(part);
01153   return elp;
01154 }
01155 
01156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01157 
01158 G4VEnergyLossProcess* 
01159 G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
01160                                   const G4String& processName)
01161 {
01162   G4VEnergyLossProcess* proc = 0;
01163   const std::vector<G4VEnergyLossProcess*> v = 
01164     G4LossTableManager::Instance()->GetEnergyLossProcessVector();
01165   G4int n = v.size();
01166   for(G4int i=0; i<n; ++i) {
01167     if((v[i])->GetProcessName() == processName) {
01168       G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01169       if(ActiveForParticle(part, p)) {
01170         proc = v[i];
01171         break;
01172       }
01173     }
01174   }
01175   return proc;
01176 }
01177 
01178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01179 
01180 G4VEmProcess* 
01181 G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
01182                                     const G4String& processName)
01183 {
01184   G4VEmProcess* proc = 0;
01185   const std::vector<G4VEmProcess*> v = 
01186     G4LossTableManager::Instance()->GetEmProcessVector();
01187   G4int n = v.size();
01188   for(G4int i=0; i<n; ++i) {
01189     if((v[i])->GetProcessName() == processName) {
01190       G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01191       if(ActiveForParticle(part, p)) {
01192         proc = v[i];
01193         break;
01194       }
01195     }
01196   }
01197   return proc;
01198 }
01199 
01200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01201 
01202 G4VMultipleScattering* 
01203 G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
01204                                const G4String& processName)
01205 {
01206   G4VMultipleScattering* proc = 0;
01207   const std::vector<G4VMultipleScattering*> v = 
01208     G4LossTableManager::Instance()->GetMultipleScatteringVector();
01209   G4int n = v.size();
01210   for(G4int i=0; i<n; ++i) {
01211     if((v[i])->GetProcessName() == processName) {
01212       G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
01213       if(ActiveForParticle(part, p)) {
01214         proc = v[i];
01215         break;
01216       }
01217     }
01218   }
01219   return proc;
01220 }
01221 
01222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01223 
01224 G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
01225                                          G4VProcess* proc)
01226 {
01227   G4ProcessManager* pm = part->GetProcessManager();
01228   G4ProcessVector* pv = pm->GetProcessList();
01229   G4int n = pv->size();
01230   G4bool res = false;
01231   for(G4int i=0; i<n; ++i) {
01232     if((*pv)[i] == proc) {
01233       if(pm->GetProcessActivation(i)) { res = true; }
01234       break;
01235     }
01236   }
01237   return res;
01238 }
01239 
01240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01241 
01242 void G4EmCalculator::SetVerbose(G4int verb)
01243 {
01244   verbose = verb;
01245 }
01246 
01247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01248 

Generated on Mon May 27 17:48:08 2013 for Geant4 by  doxygen 1.4.7