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