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