Geant4 Cross Reference |
1 // << 1 // 2 // ******************************************* << 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * License and Disclaimer * 4 // * << 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided << 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License << 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ << 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. << 9 // * include a list of copyright holders. * 10 // * << 10 // * * 11 // * Neither the authors of this software syst << 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin << 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran << 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum << 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio << 16 // * for the full disclaimer and the limitation of liability. * 17 // * << 17 // * * 18 // * This code implementation is the result << 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio << 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag << 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati << 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof << 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* << 24 // ******************************************************************** 25 // << 25 // 26 // << 26 // $Id: G4EmCalculator.cc,v 1.59 2011-01-03 19:34:03 vnivanch Exp $ 27 // ------------------------------------------- << 27 // GEANT4 tag $Name: geant4-09-04-patch-02 $ 28 // << 28 // 29 // GEANT4 Class file << 29 // ------------------------------------------------------------------- 30 // << 30 // 31 // << 31 // GEANT4 Class file 32 // File name: G4EmCalculator << 32 // 33 // << 33 // 34 // Author: Vladimir Ivanchenko << 34 // File name: G4EmCalculator 35 // << 35 // 36 // Creation date: 28.06.2004 << 36 // Author: Vladimir Ivanchenko 37 // << 37 // 38 // << 38 // Creation date: 28.06.2004 39 // Class Description: V.Ivanchenko & M.Novak << 39 // 40 // << 40 // Modifications: 41 // ------------------------------------------- << 41 // 12.09.2004 Add verbosity (V.Ivanchenko) 42 // << 42 // 17.11.2004 Change signature of methods, add new methods (V.Ivanchenko) 43 //....oooOO0OOooo........oooOO0OOooo........oo << 43 // 08.04.2005 Major optimisation of internal interfaces (V.Ivantchenko) 44 //....oooOO0OOooo........oooOO0OOooo........oo << 44 // 08.05.2005 Use updated interfaces (V.Ivantchenko) 45 << 45 // 23.10.2005 Fix computations for ions (V.Ivantchenko) 46 #include "G4EmCalculator.hh" << 46 // 11.01.2006 Add GetCSDARange (V.Ivantchenko) 47 #include "G4SystemOfUnits.hh" << 47 // 26.01.2006 Rename GetRange -> GetRangeFromRestricteDEDX (V.Ivanchenko) 48 #include "G4LossTableManager.hh" << 48 // 14.03.2006 correction in GetCrossSectionPerVolume (mma) 49 #include "G4EmParameters.hh" << 49 // suppress GetCrossSectionPerAtom 50 #include "G4NistManager.hh" << 50 // elm->GetA() in ComputeCrossSectionPerAtom 51 #include "G4DynamicParticle.hh" << 51 // 22.03.2006 Add ComputeElectronicDEDX and ComputeTotalDEDX (V.Ivanchenko) 52 #include "G4VEmProcess.hh" << 52 // 13.05.2006 Add Corrections for ion stopping (V.Ivanchenko) 53 #include "G4VEnergyLossProcess.hh" << 53 // 29.09.2006 Uncomment computation of smoothing factor (V.Ivanchenko) 54 #include "G4VMultipleScattering.hh" << 54 // 27.10.2006 Change test energy to access lowEnergy model from 55 #include "G4Material.hh" << 55 // 10 keV to 1 keV (V. Ivanchenko) 56 #include "G4MaterialCutsCouple.hh" << 56 // 15.03.2007 Add ComputeEnergyCutFromRangeCut methods (V.Ivanchenko) 57 #include "G4ParticleDefinition.hh" << 57 // 21.04.2008 Updated computations for ions (V.Ivanchenko) 58 #include "G4ParticleTable.hh" << 58 // 59 #include "G4IonTable.hh" << 59 // Class Description: 60 #include "G4PhysicsTable.hh" << 60 // 61 #include "G4ProductionCutsTable.hh" << 61 // ------------------------------------------------------------------- 62 #include "G4ProcessManager.hh" << 62 // 63 #include "G4ionEffectiveCharge.hh" << 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 64 #include "G4RegionStore.hh" << 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 65 #include "G4Element.hh" << 65 66 #include "G4EmCorrections.hh" << 66 #include "G4EmCalculator.hh" 67 #include "G4GenericIon.hh" << 67 #include "G4LossTableManager.hh" 68 #include "G4ProcessVector.hh" << 68 #include "G4VEmProcess.hh" 69 #include "G4Gamma.hh" << 69 #include "G4VEnergyLossProcess.hh" 70 #include "G4Electron.hh" << 70 #include "G4VMultipleScattering.hh" 71 #include "G4Positron.hh" << 71 #include "G4Material.hh" 72 #include "G4EmUtility.hh" << 72 #include "G4MaterialCutsCouple.hh" 73 << 73 #include "G4ParticleDefinition.hh" 74 //....oooOO0OOooo........oooOO0OOooo........oo << 74 #include "G4ParticleTable.hh" 75 << 75 #include "G4PhysicsTable.hh" 76 G4EmCalculator::G4EmCalculator() << 76 #include "G4ProductionCutsTable.hh" 77 { << 77 #include "G4ProcessManager.hh" 78 manager = G4LossTableManager::Instance(); << 78 #include "G4ionEffectiveCharge.hh" 79 nist = G4NistManager::Instance(); << 79 #include "G4RegionStore.hh" 80 theParameters = G4EmParameters::Instance(); << 80 #include "G4Element.hh" 81 corr = manager->EmCorrections(); << 81 #include "G4EmCorrections.hh" 82 cutenergy[0] = cutenergy[1] = cutenergy[2] = << 82 #include "G4GenericIon.hh" 83 theGenericIon = G4GenericIon::GenericIon(); << 83 84 ionEffCharge = new G4ionEffectiveCharge(); << 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 dynParticle = new G4DynamicParticle(); << 85 86 ionTable = G4ParticleTable::GetParticle << 86 G4EmCalculator::G4EmCalculator() 87 } << 87 { 88 << 88 manager = G4LossTableManager::Instance(); 89 //....oooOO0OOooo........oooOO0OOooo........oo << 89 corr = manager->EmCorrections(); 90 << 90 nLocalMaterials = 0; 91 G4EmCalculator::~G4EmCalculator() << 91 verbose = 0; 92 { << 92 currentCoupleIndex = 0; 93 delete ionEffCharge; << 93 currentCouple = 0; 94 delete dynParticle; << 94 currentMaterial = 0; 95 for (G4int i=0; i<nLocalMaterials; ++i) { << 95 currentParticle = 0; 96 delete localCouples[i]; << 96 lambdaParticle = 0; 97 } << 97 baseParticle = 0; 98 } << 98 currentLambda = 0; 99 << 99 currentModel = 0; 100 //....oooOO0OOooo........oooOO0OOooo........oo << 100 currentProcess = 0; 101 << 101 loweModel = 0; 102 G4double G4EmCalculator::GetDEDX(G4double kinE << 102 chargeSquare = 1.0; 103 const G4Parti << 103 massRatio = 1.0; 104 const G4Mater << 104 mass = 0.0; 105 const G4Regio << 105 currentCut = 0.0; 106 { << 106 currentParticleName= ""; 107 G4double res = 0.0; << 107 currentMaterialName= ""; 108 const G4MaterialCutsCouple* couple = FindCou << 108 currentName = ""; 109 if(nullptr != couple && UpdateParticle(p, ki << 109 lambdaName = ""; 110 res = manager->GetDEDX(p, kinEnergy, coupl << 110 theGenericIon = G4GenericIon::GenericIon(); 111 << 111 ionEffCharge = new G4ionEffectiveCharge(); 112 if(isIon) { << 112 isIon = false; 113 if(FindEmModel(p, currentProcessName, ki << 113 isApplicable = false; 114 G4double length = CLHEP::nm; << 114 } 115 G4double eloss = res*length; << 115 116 //G4cout << "### GetDEDX: E= " << kinE << 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 117 // << " de= " << eloss << G4endl << 117 118 dynParticle->SetKineticEnergy(kinEnerg << 118 G4EmCalculator::~G4EmCalculator() 119 currentModel->GetChargeSquareRatio(p, << 119 { 120 currentModel->CorrectionsAlongStep(cou << 120 delete ionEffCharge; 121 res = eloss/length; << 121 for (G4int i=0; i<nLocalMaterials; ++i) { 122 //G4cout << " de1= " << eloss << << 122 delete localCouples[i]; 123 // << " " << p->GetParticleName( << 123 } 124 } << 124 } 125 } << 125 126 << 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 127 if(verbose>0) { << 127 128 G4cout << "G4EmCalculator::GetDEDX: E(Me << 128 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p, 129 << " DEDX(MeV/mm)= " << res*mm/Me << 129 const G4Material* mat, const G4Region* region) 130 << " DEDX(MeV*cm^2/g)= " << res*g << 130 { 131 << " " << p->GetParticleName() << 131 G4double res = 0.0; 132 << " in " << mat->GetName() << 132 const G4MaterialCutsCouple* couple = FindCouple(mat, region); 133 << " isIon= " << isIon << 133 if(couple && UpdateParticle(p, kinEnergy) ) { 134 << G4endl; << 134 res = manager->GetDEDX(p, kinEnergy, couple); 135 } << 135 136 } << 136 if(isIon) { 137 return res; << 137 if(FindEmModel(p, currentProcessName, kinEnergy)) { 138 } << 138 G4double length = CLHEP::nm; 139 << 139 G4double eloss = res*length; 140 //....oooOO0OOooo........oooOO0OOooo........oo << 140 //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res 141 << 141 // << " de= " << eloss << G4endl;; 142 G4double G4EmCalculator::GetRangeFromRestricte << 142 G4double niel = 0.0; 143 << 143 dynParticle.SetKineticEnergy(kinEnergy); 144 << 144 currentModel->GetChargeSquareRatio(p, mat, kinEnergy); 145 << 145 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length); 146 { << 146 res = eloss/length; 147 G4double res = 0.0; << 147 //G4cout << " de1= " << eloss << " res1= " << res 148 const G4MaterialCutsCouple* couple = FindCou << 148 // << " " << p->GetParticleName() <<G4endl;; 149 if(couple && UpdateParticle(p, kinEnergy)) { << 149 } 150 res = manager->GetRangeFromRestricteDEDX(p << 150 } 151 if(verbose>1) { << 151 152 G4cout << " G4EmCalculator::GetRangeFrom << 152 if(verbose>0) { 153 << kinEnergy/MeV << 153 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV 154 << " range(mm)= " << res/mm << 154 << " DEDX(MeV/mm)= " << res*mm/MeV 155 << " " << p->GetParticleName() << 155 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 156 << " in " << mat->GetName() << 156 << " " << p->GetParticleName() 157 << G4endl; << 157 << " in " << mat->GetName() 158 } << 158 << " isIon= " << isIon 159 } << 159 << G4endl; 160 return res; << 160 } 161 } << 161 } 162 << 162 return res; 163 //....oooOO0OOooo........oooOO0OOooo........oo << 163 } 164 << 164 165 G4double G4EmCalculator::GetCSDARange(G4double << 165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 166 const G4 << 166 167 const G4 << 167 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle, 168 const G4 << 168 const G4String& material, const G4String& reg) 169 { << 169 { 170 G4double res = 0.0; << 170 return GetDEDX(kinEnergy,FindParticle(particle), 171 if(!theParameters->BuildCSDARange()) { << 171 FindMaterial(material),FindRegion(reg)); 172 G4ExceptionDescription ed; << 172 } 173 ed << "G4EmCalculator::GetCSDARange: CSDA << 173 174 << " use UI command: /process/eLoss/CSD << 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 175 G4Exception("G4EmCalculator::GetCSDARange" << 175 176 JustWarning, ed); << 176 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy, 177 return res; << 177 const G4ParticleDefinition* p, 178 } << 178 const G4Material* mat, 179 << 179 const G4Region* region) 180 const G4MaterialCutsCouple* couple = FindCou << 180 { 181 if(nullptr != couple && UpdateParticle(p, ki << 181 G4double res = 0.0; 182 res = manager->GetCSDARange(p, kinEnergy, << 182 const G4MaterialCutsCouple* couple = FindCouple(mat,region); 183 if(verbose>1) { << 183 if(couple && UpdateParticle(p, kinEnergy)) { 184 G4cout << " G4EmCalculator::GetCSDARange << 184 res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple); 185 << " range(mm)= " << res/mm << 185 if(verbose>0) { 186 << " " << p->GetParticleName() << 186 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV 187 << " in " << mat->GetName() << 187 << " range(mm)= " << res/mm 188 << G4endl; << 188 << " " << p->GetParticleName() 189 } << 189 << " in " << mat->GetName() 190 } << 190 << G4endl; 191 return res; << 191 } 192 } << 192 } 193 << 193 return res; 194 //....oooOO0OOooo........oooOO0OOooo........oo << 194 } 195 << 195 196 G4double G4EmCalculator::GetRange(G4double kin << 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 197 const G4Part << 197 198 const G4Mate << 198 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy, 199 const G4Regi << 199 const G4ParticleDefinition* p, 200 { << 200 const G4Material* mat, 201 G4double res = 0.0; << 201 const G4Region* region) 202 if(theParameters->BuildCSDARange()) { << 202 { 203 res = GetCSDARange(kinEnergy, p, mat, regi << 203 G4double res = 0.0; 204 } else { << 204 const G4MaterialCutsCouple* couple = FindCouple(mat,region); 205 res = GetRangeFromRestricteDEDX(kinEnergy, << 205 if(couple && UpdateParticle(p, kinEnergy)) { 206 } << 206 res = manager->GetCSDARange(p, kinEnergy, couple); 207 return res; << 207 if(verbose>0) { 208 } << 208 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV 209 << 209 << " range(mm)= " << res/mm 210 //....oooOO0OOooo........oooOO0OOooo........oo << 210 << " " << p->GetParticleName() 211 << 211 << " in " << mat->GetName() 212 G4double G4EmCalculator::GetKinEnergy(G4double << 212 << G4endl; 213 const G4 << 213 } 214 const G4 << 214 } 215 const G4 << 215 return res; 216 { << 216 } 217 G4double res = 0.0; << 217 218 const G4MaterialCutsCouple* couple = FindCou << 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 219 if(nullptr != couple && UpdateParticle(p, 1. << 219 220 res = manager->GetEnergy(p, range, couple) << 220 G4double G4EmCalculator::GetRange(G4double kinEnergy, 221 if(verbose>0) { << 221 const G4ParticleDefinition* p, 222 G4cout << "G4EmCalculator::GetKinEnergy: << 222 const G4Material* mat, 223 << " KinE(MeV)= " << res/MeV << 223 const G4Region* region) 224 << " " << p->GetParticleName() << 224 { 225 << " in " << mat->GetName() << 225 G4double res = 0.0; 226 << G4endl; << 226 const G4MaterialCutsCouple* couple = FindCouple(mat,region); 227 } << 227 if(couple && UpdateParticle(p, kinEnergy)) { 228 } << 228 res = manager->GetRange(p, kinEnergy, couple); 229 return res; << 229 if(verbose>0) { 230 } << 230 G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV 231 << 231 << " range(mm)= " << res/mm 232 //....oooOO0OOooo........oooOO0OOooo........oo << 232 << " " << p->GetParticleName() 233 << 233 << " in " << mat->GetName() 234 G4double G4EmCalculator::GetCrossSectionPerVol << 234 << G4endl; 235 co << 235 } 236 co << 236 } 237 co << 237 return res; 238 co << 238 } 239 { << 239 240 G4double res = 0.0; << 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 241 const G4MaterialCutsCouple* couple = FindCou << 241 242 << 242 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy, 243 if(nullptr != couple && UpdateParticle(p, ki << 243 const G4String& particle, 244 if(FindEmModel(p, processName, kinEnergy)) << 244 const G4String& material, 245 G4int idx = couple->GetIndex(); << 245 const G4String& reg) 246 G4int procType = -1; << 246 { 247 FindLambdaTable(p, processName, kinEnerg << 247 return GetRangeFromRestricteDEDX(kinEnergy,FindParticle(particle), 248 << 248 FindMaterial(material),FindRegion(reg)); 249 G4VEmProcess* emproc = FindDiscreteProce << 249 } 250 if(nullptr != emproc) { << 250 251 res = emproc->GetCrossSection(kinEnergy, cou << 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 252 } else if(currentLambda) { << 252 253 // special tables are built for Msc mo << 253 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy, 254 // procType is set in FindLambdaTable << 254 const G4String& particle, 255 if(procType==2) { << 255 const G4String& material, 256 auto mscM = static_cast<G4VMscModel* << 256 const G4String& reg) 257 mscM->SetCurrentCouple(couple); << 257 { 258 G4double tr1Mfp = mscM->GetTransport << 258 return GetCSDARange(kinEnergy,FindParticle(particle), 259 if (tr1Mfp<DBL_MAX) { << 259 FindMaterial(material),FindRegion(reg)); 260 res = 1./tr1Mfp; << 260 } 261 } << 261 262 } else { << 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 263 G4double e = kinEnergy*massRatio; << 263 264 res = (((*currentLambda)[idx])->Valu << 264 G4double G4EmCalculator::GetRange(G4double kinEnergy, 265 } << 265 const G4String& particle, 266 } else { << 266 const G4String& material, 267 res = ComputeCrossSectionPerVolume(kin << 267 const G4String& reg) 268 } << 268 { 269 if(verbose>0) { << 269 return GetRange(kinEnergy,FindParticle(particle), 270 G4cout << "G4EmCalculator::GetXSPerVol << 270 FindMaterial(material),FindRegion(reg)); 271 << " cross(cm-1)= " << res*cm << 271 } 272 << " " << p->GetParticleName( << 272 273 << " in " << mat->GetName(); << 273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 274 if(verbose>1) << 274 275 G4cout << " idx= " << idx << " Esc << 275 G4double G4EmCalculator::GetKinEnergy(G4double range, 276 << kinEnergy*massRatio << 276 const G4ParticleDefinition* p, 277 << " q2= " << chargeSquare; << 277 const G4Material* mat, 278 G4cout << G4endl; << 278 const G4Region* region) 279 } << 279 { 280 } << 280 G4double res = 0.0; 281 } << 281 const G4MaterialCutsCouple* couple = FindCouple(mat,region); 282 return res; << 282 if(couple && UpdateParticle(p, 1.0*GeV)) { 283 } << 283 res = manager->GetEnergy(p, range, couple); 284 << 284 if(verbose>0) { 285 //....oooOO0OOooo........oooOO0OOooo........oo << 285 G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm 286 << 286 << " KinE(MeV)= " << res/MeV 287 G4double G4EmCalculator::GetShellIonisationCro << 287 << " " << p->GetParticleName() 288 const << 288 << " in " << mat->GetName() 289 G4int << 289 << G4endl; 290 G4Ato << 290 } 291 G4dou << 291 } 292 { << 292 return res; 293 G4double res = 0.0; << 293 } 294 const G4ParticleDefinition* p = FindParticle << 294 295 G4VAtomDeexcitation* ad = manager->AtomDeexc << 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 296 if(nullptr != p && nullptr != ad) { << 296 297 res = ad->GetShellIonisationCrossSectionPe << 297 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle, 298 } << 298 const G4String& material, const G4String& reg) 299 return res; << 299 { 300 } << 300 return GetKinEnergy(range,FindParticle(particle), 301 << 301 FindMaterial(material),FindRegion(reg)); 302 //....oooOO0OOooo........oooOO0OOooo........oo << 302 } 303 << 303 304 G4double G4EmCalculator::GetMeanFreePath(G4dou << 304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 305 const << 305 306 const << 306 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy, 307 const << 307 const G4ParticleDefinition* p, 308 const << 308 const G4String& processName, 309 { << 309 const G4Material* mat, 310 G4double res = DBL_MAX; << 310 const G4Region* region) 311 G4double x = GetCrossSectionPerVolume(kinEne << 311 { 312 if(x > 0.0) { res = 1.0/x; } << 312 G4double res = 0.0; 313 if(verbose>1) { << 313 const G4MaterialCutsCouple* couple = FindCouple(mat,region); 314 G4cout << "G4EmCalculator::GetMeanFreePath << 314 315 << " MFP(mm)= " << res/mm << 315 if(couple && UpdateParticle(p, kinEnergy)) { 316 << " " << p->GetParticleName() << 316 G4int idx = couple->GetIndex(); 317 << " in " << mat->GetName() << 317 FindLambdaTable(p, processName); 318 << G4endl; << 318 319 } << 319 if(currentLambda) { 320 return res; << 320 G4double e = kinEnergy*massRatio; 321 } << 321 res = (((*currentLambda)[idx])->Value(e))*chargeSquare; 322 << 322 if(verbose>0) { 323 //....oooOO0OOooo........oooOO0OOooo........oo << 323 G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV 324 << 324 << " cross(cm-1)= " << res*cm 325 void G4EmCalculator::PrintDEDXTable(const G4Pa << 325 << " " << p->GetParticleName() 326 { << 326 << " in " << mat->GetName(); 327 const G4VEnergyLossProcess* elp = manager->G << 327 if(verbose>1) 328 G4cout << "##### DEDX Table for " << p->GetP << 328 G4cout << " idx= " << idx << " Escaled((MeV)= " << e 329 if(nullptr != elp) G4cout << *(elp->DEDXTabl << 329 << " q2= " << chargeSquare; 330 } << 330 G4cout << G4endl; 331 << 331 } 332 //....oooOO0OOooo........oooOO0OOooo........oo << 332 } 333 << 333 } 334 void G4EmCalculator::PrintRangeTable(const G4P << 334 return res; 335 { << 335 } 336 const G4VEnergyLossProcess* elp = manager->G << 336 337 G4cout << "##### Range Table for " << p->Get << 337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 338 if(nullptr != elp) G4cout << *(elp->RangeTab << 338 339 } << 339 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy, 340 << 340 const G4String& particle, 341 //....oooOO0OOooo........oooOO0OOooo........oo << 341 const G4String& processName, 342 << 342 const G4String& material, 343 void G4EmCalculator::PrintInverseRangeTable(co << 343 const G4String& reg) 344 { << 344 { 345 const G4VEnergyLossProcess* elp = manager->G << 345 return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName, 346 G4cout << "### G4EmCalculator: Inverse Range << 346 FindMaterial(material),FindRegion(reg)); 347 << p->GetParticleName() << G4endl; << 347 } 348 if(nullptr != elp) G4cout << *(elp->InverseR << 348 349 } << 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 350 << 350 351 //....oooOO0OOooo........oooOO0OOooo........oo << 351 G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom( 352 << 352 const G4String& particle, 353 G4double G4EmCalculator::ComputeDEDX(G4double << 353 G4int Z, 354 const G4P << 354 G4AtomicShellEnumerator shell, 355 const G4S << 355 G4double kinEnergy) 356 const G4M << 356 { 357 G4d << 357 G4double res = 0.0; 358 { << 358 const G4ParticleDefinition* p = FindParticle(particle); 359 SetupMaterial(mat); << 359 G4VAtomDeexcitation* ad = manager->AtomDeexcitation(); 360 G4double res = 0.0; << 360 if(p && ad) { 361 if(verbose > 1) { << 361 res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy); 362 G4cout << "### G4EmCalculator::ComputeDEDX << 362 } 363 << " in " << currentMaterialName << 363 return res; 364 << " e(MeV)= " << kinEnergy/MeV << << 364 } 365 << G4endl; << 365 366 } << 366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 367 if(UpdateParticle(p, kinEnergy)) { << 367 368 if(FindEmModel(p, processName, kinEnergy)) << 368 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy, 369 G4double escaled = kinEnergy*massRatio; << 369 const G4ParticleDefinition* p, 370 if(nullptr != baseParticle) { << 370 const G4String& processName, 371 res = currentModel->ComputeDEDXPerVolume(mat << 371 const G4Material* mat, 372 << 372 const G4Region* region) 373 if(verbose > 1) { << 373 { 374 G4cout << "Particle: " << p->GetParticleNa << 374 G4double res = DBL_MAX; 375 << " E(MeV)=" << kinEnergy << 375 G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region); 376 << " Base particle: " << baseParticle->Ge << 376 if(x > 0.0) { res = 1.0/x; } 377 << " Escaled(MeV)= " << escaled << 377 if(verbose>1) { 378 << " q2=" << chargeSquare << G4endl; << 378 G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV 379 } << 379 << " MFP(mm)= " << res/mm 380 } else { << 380 << " " << p->GetParticleName() 381 res = currentModel->ComputeDEDXPerVolume(mat << 381 << " in " << mat->GetName() 382 if(verbose > 1) { << 382 << G4endl; 383 G4cout << "Particle: " << p->GetParticleNa << 383 } 384 << " E(MeV)=" << kinEnergy << G4endl; << 384 return res; 385 } << 385 } 386 } << 386 387 if(verbose > 1) { << 387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 388 G4cout << currentModel->GetName() << ": DEDX << 388 389 << " DEDX(MeV*cm^2/g)= " << 389 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy, 390 << res*gram/(MeV*cm2*mat->GetDensity( << 390 const G4String& particle, 391 << G4endl; << 391 const G4String& processName, 392 } << 392 const G4String& material, 393 // emulate smoothing procedure << 393 const G4String& reg) 394 if(applySmoothing && nullptr != loweMode << 394 { 395 G4double eth = currentModel->LowEnergyLimit( << 395 return GetMeanFreePath(kinEnergy,FindParticle(particle),processName, 396 G4double res0 = 0.0; << 396 FindMaterial(material),FindRegion(reg)); 397 G4double res1 = 0.0; << 397 } 398 if(nullptr != baseParticle) { << 398 399 res1 = chargeSquare* << 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 400 currentModel->ComputeDEDXPerVolume(mat, << 400 401 res0 = chargeSquare* << 401 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p) 402 loweModel->ComputeDEDXPerVolume(mat, bas << 402 { 403 } else { << 403 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p); 404 res1 = currentModel->ComputeDEDXPerVolume( << 404 G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl; 405 res0 = loweModel->ComputeDEDXPerVolume(mat << 405 if(elp) G4cout << *(elp->DEDXTable()) << G4endl; 406 } << 406 } 407 if(res1 > 0.0 && escaled > 0.0) { << 407 408 res *= (1.0 + (res0/res1 - 1.0)*eth/escale << 408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 409 } << 409 410 if(verbose > 1) { << 410 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p) 411 G4cout << "At boundary energy(MeV)= " << e << 411 { 412 << " DEDX(MeV/mm)= " << res0*mm/MeV << " << 412 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p); 413 << " after correction DEDX(MeV/mm)=" << r << 413 G4cout << "##### Range Table for " << p->GetParticleName() << G4endl; 414 } << 414 if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl; 415 } << 415 } 416 // correction for ions << 416 417 if(isIon) { << 417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 418 const G4double length = CLHEP::nm; << 418 419 if(UpdateCouple(mat, cut)) { << 419 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p) 420 G4double eloss = res*length; << 420 { 421 dynParticle->SetKineticEnergy(kinEnergy); << 421 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p); 422 currentModel->CorrectionsAlongStep(current << 422 G4cout << "### G4EmCalculator: Inverse Range Table for " 423 l << 423 << p->GetParticleName() << G4endl; 424 res = eloss/length; << 424 if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl; 425 << 425 } 426 if(verbose > 1) { << 426 427 G4cout << "After Corrections: DEDX(MeV/m << 427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 428 << " DEDX(MeV*cm^2/g)= " << 428 429 << res*gram/(MeV*cm2*mat->GetDensity()) << 429 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy, 430 } << 430 const G4ParticleDefinition* p, 431 } << 431 const G4String& processName, 432 } << 432 const G4Material* mat, 433 if(verbose > 0) { << 433 G4double cut) 434 G4cout << "## E(MeV)= " << kinEnergy/MeV << 434 { 435 << " DEDX(MeV/mm)= " << res*mm/MeV << 435 currentMaterial = mat; 436 << " DEDX(MeV*cm^2/g)= " << res*gram/ << 436 currentMaterialName = mat->GetName(); 437 << " cut(MeV)= " << cut/MeV << 437 G4double res = 0.0; 438 << " " << p->GetParticleName() << 438 if(verbose > 1) { 439 << " in " << currentMaterialName << 439 G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName() 440 << " Zi^2= " << chargeSquare << 440 << " in " << currentMaterialName 441 << " isIon=" << isIon << 441 << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV 442 << G4endl; << 442 << G4endl; 443 } << 443 } 444 } << 444 if(UpdateParticle(p, kinEnergy)) { 445 } << 445 if(FindEmModel(p, processName, kinEnergy)) { 446 return res; << 446 G4double escaled = kinEnergy*massRatio; 447 } << 447 if(baseParticle) { 448 << 448 res = currentModel->ComputeDEDXPerVolume( 449 //....oooOO0OOooo........oooOO0OOooo........oo << 449 mat, baseParticle, escaled, cut) * chargeSquare; 450 << 450 if(verbose > 1) { 451 G4double G4EmCalculator::ComputeElectronicDEDX << 451 G4cout << baseParticle->GetParticleName() 452 << 452 << " Escaled(MeV)= " << escaled; 453 << 453 } 454 << 454 } else { 455 { << 455 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut); 456 SetupMaterial(mat); << 456 if(verbose > 1) G4cout << " no basePart E(MeV)= " << kinEnergy; 457 G4double dedx = 0.0; << 457 } 458 if(UpdateParticle(part, kinEnergy)) { << 458 if(verbose > 1) { 459 << 459 G4cout << " DEDX(MeV/mm)= " << res*mm/MeV 460 G4LossTableManager* lManager = G4LossTable << 460 << " DEDX(MeV*cm^2/g)= " 461 const std::vector<G4VEnergyLossProcess*> v << 461 << res*gram/(MeV*cm2*mat->GetDensity()) 462 lManager->GetEnergyLossProcessVector(); << 462 << G4endl; 463 std::size_t n = vel.size(); << 463 } 464 << 464 465 //G4cout << "ComputeElectronicDEDX for " < << 465 // emulate smoothing procedure 466 // << " n= " << n << G4endl; << 466 G4double eth = currentModel->LowEnergyLimit(); 467 << 467 // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl; 468 for(std::size_t i=0; i<n; ++i) { << 468 if(loweModel) { 469 if(vel[i]) { << 469 G4double res0 = 0.0; 470 auto p = static_cast<G4VProcess*>(vel[ << 470 G4double res1 = 0.0; 471 if(ActiveForParticle(part, p)) { << 471 if(baseParticle) { 472 //G4cout << "idx= " << i << " " << ( << 472 res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut) 473 // << " " << (vel[i])->Particle()- << 473 * chargeSquare; 474 dedx += ComputeDEDX(kinEnergy,part,( << 474 res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut) 475 } << 475 * chargeSquare; 476 } << 476 } else { 477 } << 477 res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut); 478 } << 478 res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut); 479 return dedx; << 479 } 480 } << 480 if(verbose > 1) { 481 << 481 G4cout << "At boundary energy(MeV)= " << eth/MeV 482 //....oooOO0OOooo........oooOO0OOooo........oo << 482 << " DEDX(MeV/mm)= " << res1*mm/MeV 483 << 483 << G4endl; 484 G4double << 484 } 485 G4EmCalculator::ComputeDEDXForCutInRange(G4dou << 485 486 const << 486 //G4cout << "eth= " << eth << " escaled= " << escaled 487 const << 487 // << " res0= " << res0 << " res1= " 488 G4dou << 488 // << res1 << " q2= " << chargeSquare << G4endl; 489 { << 489 490 SetupMaterial(mat); << 490 if(res1 > 0.0 && escaled > 0.0) { 491 G4double dedx = 0.0; << 491 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled); 492 if(UpdateParticle(part, kinEnergy)) { << 492 } 493 << 493 } 494 G4LossTableManager* lManager = G4LossTable << 494 495 const std::vector<G4VEnergyLossProcess*> v << 495 // low energy correction for ions 496 lManager->GetEnergyLossProcessVector(); << 496 if(isIon) { 497 std::size_t n = vel.size(); << 497 G4double length = CLHEP::nm; 498 << 498 const G4Region* r = 0; 499 if(mat != cutMaterial) { << 499 const G4MaterialCutsCouple* couple = FindCouple(mat, r); 500 cutMaterial = mat; << 500 G4double eloss = res*length; 501 cutenergy[0] = << 501 G4double niel = 0.0; 502 ComputeEnergyCutFromRangeCut(rangecut, << 502 dynParticle.SetKineticEnergy(kinEnergy); 503 cutenergy[1] = << 503 currentModel->GetChargeSquareRatio(p, mat, kinEnergy); 504 ComputeEnergyCutFromRangeCut(rangecut, << 504 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length); 505 cutenergy[2] = << 505 res = eloss/length; 506 ComputeEnergyCutFromRangeCut(rangecut, << 506 507 } << 507 if(verbose > 1) { 508 << 508 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV 509 //G4cout << "ComputeElectronicDEDX for " < << 509 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 510 // << " n= " << n << G4endl; << 510 << G4endl; 511 << 511 } 512 for(std::size_t i=0; i<n; ++i) { << 512 } 513 if(vel[i]) { << 513 } 514 auto p = static_cast<G4VProcess*>(vel[ << 514 515 if(ActiveForParticle(part, p)) { << 515 if(verbose > 0) { 516 //G4cout << "idx= " << i << " " << ( << 516 G4cout << "Sum: E(MeV)= " << kinEnergy/MeV 517 // << " " << (vel[i])->Particle()-> << 517 << " DEDX(MeV/mm)= " << res*mm/MeV 518 const G4ParticleDefinition* sec = (v << 518 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity()) 519 std::size_t idx = 0; << 519 << " cut(MeV)= " << cut/MeV 520 if(sec == G4Electron::Electron()) { << 520 << " " << p->GetParticleName() 521 else if(sec == G4Positron::Positron( << 521 << " in " << currentMaterialName 522 << 522 << " Zi^2= " << chargeSquare 523 dedx += ComputeDEDX(kinEnergy,part,( << 523 << " isIon=" << isIon 524 mat,cutenergy[id << 524 << G4endl; 525 } << 525 } 526 } << 526 } 527 } << 527 return res; 528 } << 528 } 529 return dedx; << 529 530 } << 530 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 531 << 531 532 //....oooOO0OOooo........oooOO0OOooo........oo << 532 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy, 533 << 533 const G4ParticleDefinition* part, 534 G4double G4EmCalculator::ComputeTotalDEDX(G4do << 534 const G4Material* mat, 535 cons << 535 G4double cut) 536 cons << 536 { 537 G4do << 537 currentMaterial = mat; 538 { << 538 currentMaterialName = mat->GetName(); 539 G4double dedx = ComputeElectronicDEDX(kinEne << 539 G4double dedx = 0.0; 540 if(mass > 700.*MeV) { dedx += ComputeNuclear << 540 if(UpdateParticle(part, kinEnergy)) { 541 return dedx; << 541 G4LossTableManager* lManager = G4LossTableManager::Instance(); 542 } << 542 const std::vector<G4VEnergyLossProcess*> vel = 543 << 543 lManager->GetEnergyLossProcessVector(); 544 //....oooOO0OOooo........oooOO0OOooo........oo << 544 G4int n = vel.size(); 545 << 545 for(G4int i=0; i<n; ++i) { 546 G4double G4EmCalculator::ComputeNuclearDEDX(G4 << 546 const G4ParticleDefinition* p = (vel[i])->Particle(); 547 const G4 << 547 if((!isIon && p == part) || (isIon && p == theGenericIon)) { 548 const G4 << 548 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut); 549 { << 549 } 550 G4double res = 0.0; << 550 } 551 G4VEmProcess* nucst = FindDiscreteProcess(p, << 551 } 552 if(nucst) { << 552 return dedx; 553 G4VEmModel* mod = nucst->EmModel(); << 553 } 554 if(mod) { << 554 555 mod->SetFluctuationFlag(false); << 555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 556 res = mod->ComputeDEDXPerVolume(mat, p, << 556 557 } << 557 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy, const G4String& part, 558 } << 558 const G4String& mat, G4double cut) 559 << 559 { 560 if(verbose > 1) { << 560 return ComputeElectronicDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut); 561 G4cout << p->GetParticleName() << " E(MeV << 561 } 562 << " NuclearDEDX(MeV/mm)= " << res* << 562 563 << " NuclearDEDX(MeV*cm^2/g)= " << 563 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 564 << res*gram/(MeV*cm2*mat->GetDensit << 564 565 << G4endl; << 565 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 566 } << 566 const G4ParticleDefinition* part, 567 return res; << 567 const G4Material* mat, 568 } << 568 G4double cut) 569 << 569 { 570 //....oooOO0OOooo........oooOO0OOooo........oo << 570 G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut); 571 << 571 if(mass > 700.*MeV) dedx += ComputeNuclearDEDX(kinEnergy,part,mat); 572 G4double G4EmCalculator::ComputeCrossSectionPe << 572 return dedx; 573 << 573 } 574 c << 574 575 c << 575 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 576 c << 576 577 << 577 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 578 { << 578 const G4String& part, 579 SetupMaterial(mat); << 579 const G4String& mat, 580 G4double res = 0.0; << 580 G4double cut) 581 if(UpdateParticle(p, kinEnergy)) { << 581 { 582 if(FindEmModel(p, processName, kinEnergy)) << 582 return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut); 583 G4double e = kinEnergy; << 583 } 584 G4double aCut = std::max(cut, theParamet << 584 585 if(baseParticle) { << 585 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 586 e *= kinEnergy*massRatio; << 586 587 res = currentModel->CrossSectionPerVol << 587 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy, 588 mat, baseParticle, e, aCut, e) * << 588 const G4String& particle, 589 } else { << 589 const G4String& processName, 590 res = currentModel->CrossSectionPerVol << 590 const G4String& material, 591 } << 591 G4double cut) 592 if(verbose>0) { << 592 { 593 G4cout << "G4EmCalculator::ComputeXSPe << 593 return ComputeDEDX(kinEnergy,FindParticle(particle),processName, 594 << kinEnergy/MeV << 594 FindMaterial(material),cut); 595 << " cross(cm-1)= " << res*cm << 595 } 596 << " cut(keV)= " << aCut/keV << 596 597 << " " << p->GetParticleName( << 597 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 598 << " in " << mat->GetName() << 598 599 << G4endl; << 599 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy, 600 } << 600 const G4ParticleDefinition* p, 601 } << 601 const G4Material* mat) 602 } << 602 { 603 return res; << 603 604 } << 604 G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false); 605 << 605 606 //....oooOO0OOooo........oooOO0OOooo........oo << 606 if(verbose > 1) { 607 << 607 G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV 608 G4double << 608 << " NuclearDEDX(MeV/mm)= " << res*mm/MeV 609 G4EmCalculator::ComputeCrossSectionPerAtom(G4d << 609 << " NuclearDEDX(MeV*cm^2/g)= " 610 con << 610 << res*gram/(MeV*cm2*mat->GetDensity()) 611 con << 611 << G4endl; 612 G4d << 612 } 613 G4d << 613 return res; 614 { << 614 } 615 G4double res = 0.0; << 615 616 if(UpdateParticle(p, kinEnergy)) { << 616 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 617 G4int iz = G4lrint(Z); << 617 618 CheckMaterial(iz); << 618 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy, 619 if(FindEmModel(p, processName, kinEnergy)) << 619 const G4String& particle, 620 G4double e = kinEnergy; << 620 const G4String& material) 621 G4double aCut = std::max(cut, theParamet << 621 { 622 if(baseParticle) { << 622 return ComputeNuclearDEDX(kinEnergy,FindParticle(particle), 623 e *= kinEnergy*massRatio; << 623 FindMaterial(material)); 624 currentModel->InitialiseForElement(bas << 624 } 625 res = currentModel->ComputeCrossSectio << 625 626 baseParticle, e, Z, A, aCut) * c << 626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 627 } else { << 627 628 currentModel->InitialiseForElement(p, << 628 G4double G4EmCalculator::ComputeCrossSectionPerVolume( 629 res = currentModel->ComputeCrossSectio << 629 G4double kinEnergy, 630 } << 630 const G4ParticleDefinition* p, 631 if(verbose>0) { << 631 const G4String& processName, 632 G4cout << "E(MeV)= " << kinEnergy/MeV << 632 const G4Material* mat, 633 << " cross(barn)= " << res/barn << 633 G4double cut) 634 << " " << p->GetParticleName( << 634 { 635 << " Z= " << Z << " A= " << A/ << 635 currentMaterial = mat; 636 << " cut(keV)= " << aCut/keV << 636 currentMaterialName = mat->GetName(); 637 << G4endl; << 637 G4double res = 0.0; 638 } << 638 if(UpdateParticle(p, kinEnergy)) { 639 } << 639 if(FindEmModel(p, processName, kinEnergy)) { 640 } << 640 G4double e = kinEnergy; 641 return res; << 641 if(baseParticle) { 642 } << 642 e *= kinEnergy*massRatio; 643 << 643 res = currentModel->CrossSectionPerVolume( 644 //....oooOO0OOooo........oooOO0OOooo........oo << 644 mat, baseParticle, e, cut, e) * chargeSquare; 645 << 645 } else { 646 G4double << 646 res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e); 647 G4EmCalculator::ComputeCrossSectionPerShell(G4 << 647 } 648 co << 648 if(verbose>0) { 649 co << 649 G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV 650 G4 << 650 << " cross(cm-1)= " << res*cm 651 G4 << 651 << " cut(keV)= " << cut/keV 652 { << 652 << " " << p->GetParticleName() 653 G4double res = 0.0; << 653 << " in " << mat->GetName() 654 if(UpdateParticle(p, kinEnergy)) { << 654 << G4endl; 655 CheckMaterial(Z); << 655 } 656 if(FindEmModel(p, processName, kinEnergy)) << 656 } 657 G4double e = kinEnergy; << 657 } 658 G4double aCut = std::max(cut, theParamet << 658 return res; 659 if(nullptr != baseParticle) { << 659 } 660 e *= kinEnergy*massRatio; << 660 661 currentModel->InitialiseForElement(bas << 661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 662 res = << 662 663 currentModel->ComputeCrossSectionPer << 663 G4double G4EmCalculator::ComputeCrossSectionPerVolume( 664 << 664 G4double kinEnergy, 665 } else { << 665 const G4String& particle, 666 currentModel->InitialiseForElement(p, << 666 const G4String& processName, 667 res = currentModel->ComputeCrossSectio << 667 const G4String& material, 668 } << 668 G4double cut) 669 if(verbose>0) { << 669 { 670 G4cout << "E(MeV)= " << kinEnergy/MeV << 670 return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle), 671 << " cross(barn)= " << res/barn << 671 processName, 672 << " " << p->GetParticleName( << 672 FindMaterial(material),cut); 673 << " Z= " << Z << " shellIdx= << 673 } 674 << " cut(keV)= " << aCut/keV << 674 675 << G4endl; << 675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 676 } << 676 677 } << 677 G4double G4EmCalculator::ComputeCrossSectionPerAtom( 678 } << 678 G4double kinEnergy, 679 return res; << 679 const G4ParticleDefinition* p, 680 } << 680 const G4String& processName, 681 << 681 G4double Z, G4double A, 682 //....oooOO0OOooo........oooOO0OOooo........oo << 682 G4double cut) 683 << 683 { 684 G4double << 684 G4double res = 0.0; 685 G4EmCalculator::ComputeGammaAttenuationLength( << 685 if(UpdateParticle(p, kinEnergy)) { 686 << 686 if(FindEmModel(p, processName, kinEnergy)) { 687 { << 687 G4double e = kinEnergy; 688 G4double res = 0.0; << 688 if(baseParticle) { 689 const G4ParticleDefinition* gamma = G4Gamma: << 689 e *= kinEnergy*massRatio; 690 res += ComputeCrossSectionPerVolume(kinEnerg << 690 res = currentModel->ComputeCrossSectionPerAtom( 691 res += ComputeCrossSectionPerVolume(kinEnerg << 691 baseParticle, e, Z, A, cut) * chargeSquare; 692 res += ComputeCrossSectionPerVolume(kinEnerg << 692 } else { 693 res += ComputeCrossSectionPerVolume(kinEnerg << 693 res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut); 694 if(res > 0.0) { res = 1.0/res; } << 694 } 695 return res; << 695 if(verbose>0) { 696 } << 696 G4cout << "E(MeV)= " << kinEnergy/MeV 697 << 697 << " cross(barn)= " << res/barn 698 //....oooOO0OOooo........oooOO0OOooo........oo << 698 << " " << p->GetParticleName() 699 << 699 << " Z= " << Z << " A= " << A/(g/mole) << " g/mole" 700 G4double G4EmCalculator::ComputeShellIonisatio << 700 << G4endl; 701 const << 701 } 702 G4int << 702 } 703 G4Ato << 703 } 704 G4dou << 704 return res; 705 const << 705 } 706 { << 706 707 G4double res = 0.0; << 707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 708 const G4ParticleDefinition* p = FindParticle << 708 709 G4VAtomDeexcitation* ad = manager->AtomDeexc << 709 G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy, 710 if(p && ad) { << 710 const G4String& particle, 711 res = ad->ComputeShellIonisationCrossSecti << 711 const G4String& processName, 712 << 712 const G4Element* elm, 713 } << 713 G4double cut) 714 return res; << 714 { 715 } << 715 return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle), 716 << 716 processName, 717 //....oooOO0OOooo........oooOO0OOooo........oo << 717 elm->GetZ(),elm->GetN(),cut); 718 << 718 } 719 G4double G4EmCalculator::ComputeMeanFreePath(G << 719 720 c << 720 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 721 c << 721 722 c << 722 G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom( 723 G << 723 const G4String& particle, 724 { << 724 G4int Z, 725 G4double mfp = DBL_MAX; << 725 G4AtomicShellEnumerator shell, 726 G4double x = << 726 G4double kinEnergy, 727 ComputeCrossSectionPerVolume(kinEnergy, p, << 727 const G4Material* mat) 728 if(x > 0.0) { mfp = 1.0/x; } << 728 { 729 if(verbose>1) { << 729 G4double res = 0.0; 730 G4cout << "E(MeV)= " << kinEnergy/MeV << 730 const G4ParticleDefinition* p = FindParticle(particle); 731 << " MFP(mm)= " << mfp/mm << 731 G4VAtomDeexcitation* ad = manager->AtomDeexcitation(); 732 << " " << p->GetParticleName() << 732 if(p && ad) { 733 << " in " << mat->GetName() << 733 res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell, 734 << G4endl; << 734 kinEnergy, mat); 735 } << 735 } 736 return mfp; << 736 return res; 737 } << 737 } 738 << 738 739 //....oooOO0OOooo........oooOO0OOooo........oo << 739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 740 << 740 741 G4double G4EmCalculator::ComputeEnergyCutFromR << 741 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy, 742 G4double range, << 742 const G4ParticleDefinition* p, 743 const G4ParticleDefin << 743 const G4String& processName, 744 const G4Material* mat << 744 const G4Material* mat, 745 { << 745 G4double cut) 746 return G4ProductionCutsTable::GetProductionC << 746 { 747 ConvertRangeToEnergy(part, mat, range); << 747 G4double mfp = DBL_MAX; 748 } << 748 G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut); 749 << 749 if(x > 0.0) mfp = 1.0/x; 750 //....oooOO0OOooo........oooOO0OOooo........oo << 750 if(verbose>1) { 751 << 751 G4cout << "E(MeV)= " << kinEnergy/MeV 752 G4bool G4EmCalculator::UpdateParticle(const G4 << 752 << " MFP(mm)= " << mfp/mm 753 G4double << 753 << " " << p->GetParticleName() 754 { << 754 << " in " << mat->GetName() 755 if(p != currentParticle) { << 755 << G4endl; 756 << 756 } 757 // new particle << 757 return mfp; 758 currentParticle = p; << 758 } 759 dynParticle->SetDefinition(const_cast<G4Pa << 759 760 dynParticle->SetKineticEnergy(kinEnergy); << 760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 761 baseParticle = nullptr; << 761 762 currentParticleName = p->GetParticleName() << 762 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy, 763 massRatio = 1.0; << 763 const G4String& particle, 764 mass = p->GetPDGMass(); << 764 const G4String& processName, 765 chargeSquare = 1.0; << 765 const G4String& material, 766 currentProcess = manager->GetEnergyLossPr << 766 G4double cut) 767 currentProcessName = ""; << 767 { 768 isIon = false; << 768 return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName, 769 << 769 FindMaterial(material),cut); 770 // ionisation process exist << 770 } 771 if(nullptr != currentProcess) { << 771 772 currentProcessName = currentProcess->Get << 772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 773 baseParticle = currentProcess->BaseParti << 773 774 if(currentProcessName == "ionIoni" && p- << 774 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut( 775 baseParticle = theGenericIon; << 775 G4double range, 776 isIon = true; << 776 const G4ParticleDefinition* part, 777 } << 777 const G4Material* mat) 778 << 778 { 779 // base particle is used << 779 return G4ProductionCutsTable::GetProductionCutsTable()-> 780 if(nullptr != baseParticle) { << 780 ConvertRangeToEnergy(part, mat, range); 781 massRatio = baseParticle->GetPDGMass() << 781 } 782 G4double q = p->GetPDGCharge()/basePar << 782 783 chargeSquare = q*q; << 783 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 784 } << 784 785 } << 785 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut( 786 } << 786 G4double range, 787 // Effective charge for ions << 787 const G4String& particle, 788 if(isIon && nullptr != currentProcess) { << 788 const G4String& material) 789 chargeSquare = << 789 { 790 corr->EffectiveChargeSquareRatio(p, curr << 790 return ComputeEnergyCutFromRangeCut(range,FindParticle(particle), 791 currentProcess->SetDynamicMassCharge(massR << 791 FindMaterial(material)); 792 if(verbose>1) { << 792 } 793 G4cout <<"\n NewIon: massR= "<< massRati << 793 794 << chargeSquare << " " << currentProce << 794 795 } << 795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 796 } << 796 797 return true; << 797 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p, 798 } << 798 G4double kinEnergy) 799 << 799 { 800 //....oooOO0OOooo........oooOO0OOooo........oo << 800 if(p != currentParticle) { 801 << 801 802 const G4ParticleDefinition* G4EmCalculator::Fi << 802 // new particle 803 { << 803 currentParticle = p; 804 const G4ParticleDefinition* p = nullptr; << 804 dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p)); 805 if(name != currentParticleName) { << 805 dynParticle.SetKineticEnergy(kinEnergy); 806 p = G4ParticleTable::GetParticleTable()->F << 806 baseParticle = 0; 807 if(nullptr == p) { << 807 currentParticleName = p->GetParticleName(); 808 G4cout << "### WARNING: G4EmCalculator:: << 808 massRatio = 1.0; 809 << name << G4endl; << 809 mass = p->GetPDGMass(); 810 } << 810 chargeSquare = 1.0; 811 } else { << 811 currentProcess = FindEnergyLossProcess(p); 812 p = currentParticle; << 812 currentProcessName = ""; 813 } << 813 isIon = false; 814 return p; << 814 815 } << 815 // ionisation process exist 816 << 816 if(currentProcess) { 817 //....oooOO0OOooo........oooOO0OOooo........oo << 817 currentProcessName = currentProcess->GetProcessName(); 818 << 818 baseParticle = currentProcess->BaseParticle(); 819 const G4ParticleDefinition* G4EmCalculator::Fi << 819 820 { << 820 // base particle is used 821 const G4ParticleDefinition* p = ionTable->Ge << 821 if(baseParticle) { 822 return p; << 822 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass(); 823 } << 823 G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge(); 824 << 824 chargeSquare = q*q; 825 //....oooOO0OOooo........oooOO0OOooo........oo << 825 } 826 << 826 827 const G4Material* G4EmCalculator::FindMaterial << 827 if(p->GetParticleType() == "nucleus" 828 { << 828 && currentParticleName != "deuteron" 829 if(name != currentMaterialName) { << 829 && currentParticleName != "triton" 830 SetupMaterial(G4Material::GetMaterial(name << 830 && currentParticleName != "alpha+" 831 if(nullptr == currentMaterial) { << 831 && currentParticleName != "helium" 832 G4cout << "### WARNING: G4EmCalculator:: << 832 && currentParticleName != "hydrogen" 833 << name << G4endl; << 833 ) { 834 } << 834 isIon = true; 835 } << 835 massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass(); 836 return currentMaterial; << 836 baseParticle = theGenericIon; 837 } << 837 // G4cout << p->GetParticleName() 838 << 838 // << " in " << currentMaterial->GetName() 839 //....oooOO0OOooo........oooOO0OOooo........oo << 839 // << " e= " << kinEnergy << G4endl; 840 << 840 } 841 const G4Region* G4EmCalculator::FindRegion(con << 841 } 842 { << 842 } 843 return G4EmUtility::FindRegion(reg); << 843 844 } << 844 // Effective charge for ions 845 << 845 if(isIon) { 846 //....oooOO0OOooo........oooOO0OOooo........oo << 846 chargeSquare = 847 << 847 corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy) 848 const G4MaterialCutsCouple* G4EmCalculator::Fi << 848 * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy); 849 const G4Material* << 849 if(currentProcess) { 850 const G4Region* re << 850 currentProcess->SetDynamicMassCharge(massRatio,chargeSquare); 851 { << 851 //G4cout << "NewP: massR= " << massRatio << " q2= " << chargeSquare << G4endl; 852 const G4MaterialCutsCouple* couple = nullptr << 852 } 853 SetupMaterial(material); << 853 } 854 if(nullptr != currentMaterial) { << 854 return true; 855 // Access to materials << 855 } 856 const G4ProductionCutsTable* theCoupleTabl << 856 857 G4ProductionCutsTable::GetProductionCuts << 857 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 858 const G4Region* r = region; << 858 859 if(nullptr != r) { << 859 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name) 860 couple = theCoupleTable->GetMaterialCuts << 860 { 861 << 861 const G4ParticleDefinition* p = 0; 862 } else { << 862 if(name != currentParticleName) { 863 G4RegionStore* store = G4RegionStore::Ge << 863 p = G4ParticleTable::GetParticleTable()->FindParticle(name); 864 std::size_t nr = store->size(); << 864 if(!p) { 865 if(0 < nr) { << 865 G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find " 866 for(std::size_t i=0; i<nr; ++i) { << 866 << name << G4endl; 867 couple = theCoupleTable->GetMaterial << 867 } 868 material, ((*store)[i])->GetProduc << 868 } else { 869 if(nullptr != couple) { break; } << 869 p = currentParticle; 870 } << 870 } 871 } << 871 return p; 872 } << 872 } 873 } << 873 874 if(nullptr == couple) { << 874 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 875 G4ExceptionDescription ed; << 875 876 ed << "G4EmCalculator::FindCouple: fail fo << 876 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A) 877 << currentMaterialName << ">"; << 877 { 878 if(region) { ed << " and region " << regio << 878 const G4ParticleDefinition* p = 879 G4Exception("G4EmCalculator::FindCouple", << 879 G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); 880 FatalException, ed); << 880 return p; 881 } << 881 } 882 return couple; << 882 883 } << 883 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 884 << 884 885 //....oooOO0OOooo........oooOO0OOooo........oo << 885 const G4Material* G4EmCalculator::FindMaterial(const G4String& name) 886 << 886 { 887 G4bool G4EmCalculator::UpdateCouple(const G4Ma << 887 if(name != currentMaterialName) { 888 { << 888 currentMaterial = G4Material::GetMaterial(name); 889 SetupMaterial(material); << 889 currentMaterialName = name; 890 if(!currentMaterial) { return false; } << 890 if(!currentMaterial) 891 for (G4int i=0; i<nLocalMaterials; ++i) { << 891 G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find " 892 if(material == localMaterials[i] && cut == << 892 << name << G4endl; 893 currentCouple = localCouples[i]; << 893 } 894 currentCoupleIndex = currentCouple->GetI << 894 return currentMaterial; 895 currentCut = cut; << 895 } 896 return true; << 896 897 } << 897 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 898 } << 898 899 const G4MaterialCutsCouple* cc = new G4Mater << 899 const G4Region* G4EmCalculator::FindRegion(const G4String& reg) 900 localMaterials.push_back(material); << 900 { 901 localCouples.push_back(cc); << 901 const G4Region* r = 0; 902 localCuts.push_back(cut); << 902 if(reg != "" || reg != "world") 903 ++nLocalMaterials; << 903 r = G4RegionStore::GetInstance()->GetRegion(reg); 904 currentCouple = cc; << 904 else 905 currentCoupleIndex = currentCouple->GetIndex << 905 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld"); 906 currentCut = cut; << 906 return r; 907 return true; << 907 } 908 } << 908 909 << 909 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 910 //....oooOO0OOooo........oooOO0OOooo........oo << 910 911 << 911 const G4MaterialCutsCouple* G4EmCalculator::FindCouple( 912 void G4EmCalculator::FindLambdaTable(const G4P << 912 const G4Material* material, 913 const G4S << 913 const G4Region* region) 914 G4double << 914 { 915 { << 915 if(!material) return 0; 916 // Search for the process << 916 currentMaterial = material; 917 if (!currentLambda || p != lambdaParticle || << 917 currentMaterialName = material->GetName(); 918 lambdaName = processName; << 918 // Access to materials 919 currentLambda = nullptr; << 919 const G4ProductionCutsTable* theCoupleTable= 920 lambdaParticle = p; << 920 G4ProductionCutsTable::GetProductionCutsTable(); 921 isApplicable = false; << 921 const G4Region* r = region; 922 << 922 if(!r) 923 const G4ParticleDefinition* part = (isIon) << 923 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld"); 924 << 924 925 // Search for energy loss process << 925 return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts()); 926 currentName = processName; << 926 927 currentModel = nullptr; << 927 } 928 loweModel = nullptr; << 928 929 << 929 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 930 G4VEnergyLossProcess* elproc = FindEnLossP << 930 931 if(nullptr != elproc) { << 931 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut) 932 currentLambda = elproc->LambdaTable(); << 932 { 933 proctype = 0; << 933 if(!material) return false; 934 if(nullptr != currentLambda) { << 934 currentMaterial = material; 935 isApplicable = true; << 935 currentMaterialName = material->GetName(); 936 if(verbose>1) { << 936 for (G4int i=0; i<nLocalMaterials; ++i) { 937 G4cout << "G4VEnergyLossProcess is f << 937 if(material == localMaterials[i] && cut == localCuts[i]) { 938 << G4endl; << 938 currentCouple = localCouples[i]; 939 } << 939 currentCoupleIndex = currentCouple->GetIndex(); 940 } << 940 currentCut = cut; 941 curProcess = elproc; << 941 return true; 942 return; << 942 } 943 } << 943 } 944 << 944 const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material); 945 // Search for discrete process << 945 localMaterials.push_back(material); 946 G4VEmProcess* proc = FindDiscreteProcess(p << 946 localCouples.push_back(cc); 947 if(nullptr != proc) { << 947 localCuts.push_back(cut); 948 currentLambda = proc->LambdaTable(); << 948 nLocalMaterials++; 949 proctype = 1; << 949 currentCouple = cc; 950 if(nullptr != currentLambda) { << 950 currentCoupleIndex = currentCouple->GetIndex(); 951 isApplicable = true; << 951 currentCut = cut; 952 if(verbose>1) { << 952 return true; 953 G4cout << "G4VEmProcess is found out << 953 } 954 } << 954 955 } << 955 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 956 curProcess = proc; << 956 957 return; << 957 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p, 958 } << 958 const G4String& processName) 959 << 959 { 960 // Search for msc process << 960 // Search for the process 961 G4VMultipleScattering* msc = FindMscProces << 961 if (!currentLambda || p != lambdaParticle || processName != lambdaName) { 962 if(nullptr != msc) { << 962 lambdaName = processName; 963 currentModel = msc->SelectModel(kinEnerg << 963 currentLambda = 0; 964 proctype = 2; << 964 lambdaParticle = p; 965 if(nullptr != currentModel) { << 965 966 currentLambda = currentModel->GetCross << 966 G4String partname = p->GetParticleName(); 967 if(nullptr != currentLambda) { << 967 const G4ParticleDefinition* part = p; 968 isApplicable = true; << 968 if(isIon) { part = theGenericIon; } 969 if(verbose>1) { << 969 970 G4cout << "G4VMultipleScattering i << 970 // energy loss process 971 << G4endl; << 971 G4LossTableManager* lManager = G4LossTableManager::Instance(); 972 } << 972 const std::vector<G4VEnergyLossProcess*> vel = 973 } << 973 lManager->GetEnergyLossProcessVector(); 974 } << 974 G4int n = vel.size(); 975 curProcess = msc; << 975 for(G4int i=0; i<n; ++i) { 976 } << 976 if((vel[i])->GetProcessName() == lambdaName && 977 } << 977 (vel[i])->Particle() == part) 978 } << 978 { 979 << 979 currentLambda = (vel[i])->LambdaTable(); 980 //....oooOO0OOooo........oooOO0OOooo........oo << 980 isApplicable = true; 981 << 981 if(verbose>1) { 982 G4bool G4EmCalculator::FindEmModel(const G4Par << 982 G4cout << "G4VEnergyLossProcess is found out: " 983 const G4Str << 983 << currentName << G4endl; 984 G4 << 984 } 985 { << 985 return; 986 isApplicable = false; << 986 } 987 if(nullptr == p || nullptr == currentMateria << 987 } 988 G4cout << "G4EmCalculator::FindEmModel WAR << 988 989 << " or materail defined; particle: << 989 // discrete process 990 return isApplicable; << 990 if(!currentLambda) { 991 } << 991 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector(); 992 G4String partname = p->GetParticleName(); << 992 G4int n = vem.size(); 993 G4double scaledEnergy = kinEnergy*massRatio; << 993 for(G4int i=0; i<n; ++i) { 994 const G4ParticleDefinition* part = (isIon) ? << 994 if((vem[i])->GetProcessName() == lambdaName && 995 << 995 (vem[i])->Particle() == part) 996 if(verbose > 1) { << 996 { 997 G4cout << "## G4EmCalculator::FindEmModel << 997 currentLambda = (vem[i])->LambdaTable(); 998 << " (type= " << p->GetParticleType << 998 isApplicable = true; 999 << ") and " << processName << " at << 999 if(verbose>1) { 1000 << G4endl; << 1000 G4cout << "G4VEmProcess is found out: " 1001 if(p != part) { G4cout << " GenericIon i << 1001 << currentName << G4endl; 1002 } << 1002 } 1003 << 1003 return; 1004 // Search for energy loss process << 1004 } 1005 currentName = processName; << 1005 } 1006 currentModel = nullptr; << 1006 } 1007 loweModel = nullptr; << 1007 1008 std::size_t idx = 0; << 1008 // msc process 1009 << 1009 if(!currentLambda) { 1010 G4VEnergyLossProcess* elproc = FindEnLossPr << 1010 const std::vector<G4VMultipleScattering*> vmsc = 1011 if(nullptr != elproc) { << 1011 lManager->GetMultipleScatteringVector(); 1012 currentModel = elproc->SelectModelForMate << 1012 G4int n = vmsc.size(); 1013 currentModel->InitialiseForMaterial(part, << 1013 for(G4int i=0; i<n; ++i) { 1014 currentModel->SetupForMaterial(part, curr << 1014 if((vmsc[i])->GetProcessName() == lambdaName && 1015 G4double eth = currentModel->LowEnergyLim << 1015 (vmsc[i])->Particle() == part) 1016 if(eth > 0.0) { << 1016 { 1017 loweModel = elproc->SelectModelForMater << 1017 currentLambda = (vmsc[i])->LambdaTable(); 1018 if(loweModel == currentModel) { loweMod << 1018 isApplicable = true; 1019 else { << 1019 if(verbose>1) { 1020 loweModel->InitialiseForMaterial(part << 1020 G4cout << "G4VMultipleScattering is found out: " 1021 loweModel->SetupForMaterial(part, cur << 1021 << currentName << G4endl; 1022 } << 1022 } 1023 } << 1023 return; 1024 } << 1024 } 1025 << 1025 } 1026 // Search for discrete process << 1026 } 1027 if(nullptr == currentModel) { << 1027 } 1028 G4VEmProcess* proc = FindDiscreteProcess( << 1028 } 1029 if(nullptr != proc) { << 1029 1030 currentModel = proc->SelectModelForMate << 1030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1031 currentModel->InitialiseForMaterial(par << 1031 1032 currentModel->SetupForMaterial(part, cu << 1032 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p, 1033 G4double eth = currentModel->LowEnergyL << 1033 const G4String& processName, 1034 if(eth > 0.0) { << 1034 G4double kinEnergy) 1035 loweModel = proc->SelectModelForMater << 1035 { 1036 if(loweModel == currentModel) { loweM << 1036 isApplicable = false; 1037 else { << 1037 if(!p) { 1038 loweModel->InitialiseForMaterial(pa << 1038 G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined" 1039 loweModel->SetupForMaterial(part, c << 1039 << G4endl; 1040 } << 1040 return isApplicable; 1041 } << 1041 } 1042 } << 1042 G4String partname = p->GetParticleName(); 1043 } << 1043 const G4ParticleDefinition* part = p; 1044 << 1044 G4double scaledEnergy = kinEnergy*massRatio; 1045 // Search for msc process << 1045 if(isIon) { part = theGenericIon; } 1046 if(nullptr == currentModel) { << 1046 1047 G4VMultipleScattering* proc = FindMscProc << 1047 if(verbose > 1) { 1048 if(nullptr != proc) { << 1048 G4cout << "G4EmCalculator::FindEmModel for " << partname 1049 currentModel = proc->SelectModel(kinEne << 1049 << " (type= " << p->GetParticleType() 1050 loweModel = nullptr; << 1050 << ") and " << processName << " at E(MeV)= " << scaledEnergy; 1051 } << 1051 if(p != part) G4cout << " GenericIon is the base particle"; 1052 } << 1052 G4cout << G4endl; 1053 if(nullptr != currentModel) { << 1053 } 1054 if(loweModel == currentModel) { loweModel << 1054 1055 isApplicable = true; << 1055 // Search for energy loss process 1056 currentModel->InitialiseForMaterial(part, << 1056 currentName = processName; 1057 if(loweModel) { << 1057 currentModel = 0; 1058 loweModel->InitialiseForMaterial(part, << 1058 loweModel = 0; 1059 } << 1059 size_t idx = 0; 1060 if(verbose > 1) { << 1060 G4LossTableManager* lManager = G4LossTableManager::Instance(); 1061 G4cout << " Model <" << currentModel- << 1061 const std::vector<G4VEnergyLossProcess*> vel = 1062 << "> Emin(MeV)= " << currentMod << 1062 lManager->GetEnergyLossProcessVector(); 1063 << " for " << part->GetParticleN << 1063 G4int n = vel.size(); 1064 if(nullptr != elproc) { << 1064 G4VEnergyLossProcess* elproc = 0; 1065 G4cout << " and " << elproc->GetProce << 1065 for(G4int i=0; i<n; ++i) { 1066 << G4endl; << 1066 // G4cout << "i= " << i << " part= " 1067 } << 1067 // << (vel[i])->Particle()->GetParticleName() 1068 if(nullptr != loweModel) { << 1068 // << " proc= " << (vel[i])->GetProcessName() << G4endl; 1069 G4cout << " LowEnergy model <" << low << 1069 if((vel[i])->GetProcessName() == currentName) { 1070 } << 1070 if(baseParticle) { 1071 G4cout << G4endl; << 1071 if((vel[i])->Particle() == baseParticle) { 1072 } << 1072 elproc = vel[i]; 1073 } << 1073 break; 1074 return isApplicable; << 1074 } 1075 } << 1075 } else { 1076 << 1076 if((vel[i])->Particle() == part) { 1077 //....oooOO0OOooo........oooOO0OOooo........o << 1077 elproc = vel[i]; 1078 << 1078 break; 1079 G4VEnergyLossProcess* << 1079 } 1080 G4EmCalculator::FindEnLossProcess(const G4Par << 1080 } 1081 const G4Str << 1081 } 1082 { << 1082 } 1083 G4VEnergyLossProcess* proc = nullptr; << 1083 if(elproc) { 1084 const std::vector<G4VEnergyLossProcess*> v << 1084 currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx); 1085 manager->GetEnergyLossProcessVector(); << 1085 G4double eth = currentModel->LowEnergyLimit(); 1086 std::size_t n = v.size(); << 1086 if(eth > 0.0) { 1087 for(std::size_t i=0; i<n; ++i) { << 1087 loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx); 1088 if((v[i])->GetProcessName() == processNam << 1088 } 1089 auto p = static_cast<G4VProcess*>(v[i]) << 1089 } 1090 if(ActiveForParticle(part, p)) { << 1090 1091 proc = v[i]; << 1091 // Search for discrete process 1092 break; << 1092 if(!currentModel) { 1093 } << 1093 const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector(); 1094 } << 1094 G4int n = vem.size(); 1095 } << 1095 for(G4int i=0; i<n; ++i) { 1096 return proc; << 1096 if((vem[i])->GetProcessName() == currentName && 1097 } << 1097 (vem[i])->Particle() == part) 1098 << 1098 { 1099 //....oooOO0OOooo........oooOO0OOooo........o << 1099 currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx); 1100 << 1100 G4double eth = currentModel->LowEnergyLimit(); 1101 G4VEmProcess* << 1101 if(eth > 0.0) { 1102 G4EmCalculator::FindDiscreteProcess(const G4P << 1102 loweModel = (vem[i])->SelectModelForMaterial(eth - CLHEP::eV, idx); 1103 const G4S << 1103 } 1104 { << 1104 break; 1105 G4VEmProcess* proc = nullptr; << 1105 } 1106 auto v = manager->GetEmProcessVector(); << 1106 } 1107 std::size_t n = v.size(); << 1107 } 1108 for(std::size_t i=0; i<n; ++i) { << 1108 1109 const G4String& pName = v[i]->GetProcessN << 1109 // Search for msc process 1110 if(pName == "GammaGeneralProc") { << 1110 if(!currentModel) { 1111 proc = v[i]->GetEmProcess(processName); << 1111 const std::vector<G4VMultipleScattering*> vmsc = 1112 break; << 1112 lManager->GetMultipleScatteringVector(); 1113 } else if(pName == processName) { << 1113 G4int n = vmsc.size(); 1114 const auto p = static_cast<G4VProcess*> << 1114 for(G4int i=0; i<n; ++i) { 1115 if(ActiveForParticle(part, p)) { << 1115 if((vmsc[i])->GetProcessName() == currentName && 1116 proc = v[i]; << 1116 (vmsc[i])->Particle() == part) 1117 break; << 1117 { 1118 } << 1118 currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx); 1119 } << 1119 G4double eth = currentModel->LowEnergyLimit(); 1120 } << 1120 if(eth > 0.0) { 1121 return proc; << 1121 loweModel = (vmsc[i])->SelectModelForMaterial(eth - CLHEP::eV, idx); 1122 } << 1122 } 1123 << 1123 break; 1124 //....oooOO0OOooo........oooOO0OOooo........o << 1124 } 1125 << 1125 } 1126 G4VMultipleScattering* << 1126 } 1127 G4EmCalculator::FindMscProcess(const G4Partic << 1127 if(currentModel) { 1128 const G4String << 1128 if(loweModel == currentModel) { loweModel = 0; } 1129 { << 1129 isApplicable = true; 1130 G4VMultipleScattering* proc = nullptr; << 1130 if(verbose > 1) { 1131 const std::vector<G4VMultipleScattering*> v << 1131 G4cout << "Model <" << currentModel->GetName() 1132 manager->GetMultipleScatteringVector(); << 1132 << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV; 1133 std::size_t n = v.size(); << 1133 if(loweModel) { 1134 for(std::size_t i=0; i<n; ++i) { << 1134 G4cout << " LowEnergy model <" << loweModel->GetName() << ">"; 1135 if((v[i])->GetProcessName() == processNam << 1135 } 1136 auto p = static_cast<G4VProcess*>(v[i]) << 1136 G4cout << G4endl; 1137 if(ActiveForParticle(part, p)) { << 1137 } 1138 proc = v[i]; << 1138 } 1139 break; << 1139 return isApplicable; 1140 } << 1140 } 1141 } << 1141 1142 } << 1142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1143 return proc; << 1143 1144 } << 1144 G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess( 1145 << 1145 const G4ParticleDefinition* p) 1146 //....oooOO0OOooo........oooOO0OOooo........o << 1146 { 1147 << 1147 G4VEnergyLossProcess* elp = 0; 1148 G4VProcess* G4EmCalculator::FindProcess(const << 1148 G4String partname = p->GetParticleName(); 1149 const << 1149 const G4ParticleDefinition* part = p; 1150 { << 1150 1151 G4VProcess* proc = nullptr; << 1151 if(p->GetParticleType() == "nucleus" 1152 const G4ProcessManager* procman = part->Get << 1152 && currentParticleName != "deuteron" 1153 G4ProcessVector* pv = procman->GetProcessLi << 1153 && currentParticleName != "triton" 1154 G4int nproc = (G4int)pv->size(); << 1154 && currentParticleName != "alpha+" 1155 for(G4int i=0; i<nproc; ++i) { << 1155 && currentParticleName != "helium" 1156 if(processName == (*pv)[i]->GetProcessNam << 1156 && currentParticleName != "hydrogen" 1157 proc = (*pv)[i]; << 1157 ) { part = theGenericIon; } 1158 break; << 1158 1159 } << 1159 G4LossTableManager* lManager = G4LossTableManager::Instance(); 1160 } << 1160 const std::vector<G4VEnergyLossProcess*> vel = 1161 return proc; << 1161 lManager->GetEnergyLossProcessVector(); 1162 } << 1162 G4int n = vel.size(); 1163 << 1163 for(G4int i=0; i<n; ++i) { 1164 //....oooOO0OOooo........oooOO0OOooo........o << 1164 if( (vel[i])->Particle() == part ) { 1165 << 1165 elp = vel[i]; 1166 G4bool G4EmCalculator::ActiveForParticle(cons << 1166 break; 1167 G4VP << 1167 } 1168 { << 1168 } 1169 G4ProcessManager* pm = part->GetProcessMana << 1169 return elp; 1170 G4ProcessVector* pv = pm->GetProcessList(); << 1170 } 1171 G4int n = (G4int)pv->size(); << 1171 1172 G4bool res = false; << 1172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1173 for(G4int i=0; i<n; ++i) { << 1173 1174 if((*pv)[i] == proc) { << 1174 void G4EmCalculator::SetVerbose(G4int verb) 1175 if(pm->GetProcessActivation(i)) { res = << 1175 { 1176 break; << 1176 verbose = verb; 1177 } << 1177 } 1178 } << 1178 1179 return res; << 1179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1180 } << 1180 1181 << 1182 //....oooOO0OOooo........oooOO0OOooo........o << 1183 << 1184 void G4EmCalculator::SetupMaterial(const G4Ma << 1185 { << 1186 if(mat) { << 1187 currentMaterial = mat; << 1188 currentMaterialName = mat->GetName(); << 1189 } else { << 1190 currentMaterial = nullptr; << 1191 currentMaterialName = ""; << 1192 } << 1193 } << 1194 << 1195 //....oooOO0OOooo........oooOO0OOooo........o << 1196 << 1197 void G4EmCalculator::SetupMaterial(const G4St << 1198 { << 1199 SetupMaterial(nist->FindOrBuildMaterial(mna << 1200 } << 1201 << 1202 //....oooOO0OOooo........oooOO0OOooo........o << 1203 << 1204 void G4EmCalculator::CheckMaterial(G4int Z) << 1205 { << 1206 G4bool isFound = false; << 1207 if(nullptr != currentMaterial) { << 1208 G4int nn = (G4int)currentMaterial->GetNum << 1209 for(G4int i=0; i<nn; ++i) { << 1210 if(Z == currentMaterial->GetElement(i)- << 1211 isFound = true; << 1212 break; << 1213 } << 1214 } << 1215 } << 1216 if(!isFound) { << 1217 SetupMaterial(nist->FindOrBuildSimpleMate << 1218 } << 1219 } << 1220 << 1221 //....oooOO0OOooo........oooOO0OOooo........o << 1222 << 1223 void G4EmCalculator::SetVerbose(G4int verb) << 1224 { << 1225 verbose = verb; << 1226 } << 1227 << 1228 //....oooOO0OOooo........oooOO0OOooo........o << 1229 << 1230 1181