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