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 // GEANT4 Class file 28 // GEANT4 Class file 29 // 29 // 30 // 30 // 31 // File name: G4VEmProcess 31 // File name: G4VEmProcess 32 // 32 // 33 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code 34 // 34 // 35 // Creation date: 01.10.2003 35 // Creation date: 01.10.2003 36 // 36 // 37 // Modifications: by V.Ivanchenko 37 // Modifications: by V.Ivanchenko 38 // 38 // 39 // Class Description: based class for discrete 39 // Class Description: based class for discrete and rest/discrete EM processes 40 // 40 // 41 41 42 // ------------------------------------------- 42 // ------------------------------------------------------------------- 43 // 43 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 46 47 #include "G4VEmProcess.hh" 47 #include "G4VEmProcess.hh" 48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4ProcessManager.hh" 50 #include "G4ProcessManager.hh" 51 #include "G4LossTableManager.hh" 51 #include "G4LossTableManager.hh" 52 #include "G4LossTableBuilder.hh" 52 #include "G4LossTableBuilder.hh" 53 #include "G4Step.hh" 53 #include "G4Step.hh" 54 #include "G4ParticleDefinition.hh" 54 #include "G4ParticleDefinition.hh" 55 #include "G4VEmModel.hh" 55 #include "G4VEmModel.hh" 56 #include "G4DataVector.hh" 56 #include "G4DataVector.hh" 57 #include "G4PhysicsTable.hh" 57 #include "G4PhysicsTable.hh" 58 #include "G4EmDataHandler.hh" 58 #include "G4EmDataHandler.hh" 59 #include "G4PhysicsLogVector.hh" 59 #include "G4PhysicsLogVector.hh" 60 #include "G4VParticleChange.hh" 60 #include "G4VParticleChange.hh" >> 61 #include "G4PhysicsModelCatalog.hh" 61 #include "G4ProductionCutsTable.hh" 62 #include "G4ProductionCutsTable.hh" 62 #include "G4Region.hh" 63 #include "G4Region.hh" 63 #include "G4Gamma.hh" 64 #include "G4Gamma.hh" 64 #include "G4Electron.hh" 65 #include "G4Electron.hh" 65 #include "G4Positron.hh" 66 #include "G4Positron.hh" 66 #include "G4PhysicsTableHelper.hh" 67 #include "G4PhysicsTableHelper.hh" 67 #include "G4EmBiasingManager.hh" 68 #include "G4EmBiasingManager.hh" 68 #include "G4EmParameters.hh" 69 #include "G4EmParameters.hh" 69 #include "G4EmProcessSubType.hh" 70 #include "G4EmProcessSubType.hh" 70 #include "G4EmTableUtil.hh" << 71 #include "G4LowEnergyEmProcessSubType.hh" 71 #include "G4EmUtility.hh" << 72 #include "G4DNAModelSubType.hh" 72 #include "G4DNAModelSubType.hh" 73 #include "G4GenericIon.hh" 73 #include "G4GenericIon.hh" 74 #include "G4Log.hh" 74 #include "G4Log.hh" 75 #include <iostream> 75 #include <iostream> 76 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 78 79 G4VEmProcess::G4VEmProcess(const G4String& nam 79 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type): 80 G4VDiscreteProcess(name, type) 80 G4VDiscreteProcess(name, type) 81 { 81 { 82 theParameters = G4EmParameters::Instance(); 82 theParameters = G4EmParameters::Instance(); 83 SetVerboseLevel(1); 83 SetVerboseLevel(1); 84 84 85 // Size of tables 85 // Size of tables 86 minKinEnergy = 0.1*CLHEP::keV; 86 minKinEnergy = 0.1*CLHEP::keV; 87 maxKinEnergy = 100.0*CLHEP::TeV; 87 maxKinEnergy = 100.0*CLHEP::TeV; 88 88 89 // default lambda factor 89 // default lambda factor 90 invLambdaFactor = 1.0/lambdaFactor; << 90 logLambdaFactor = G4Log(lambdaFactor); 91 91 92 // particle types 92 // particle types 93 theGamma = G4Gamma::Gamma(); << 93 theGamma = G4Gamma::Gamma(); 94 theElectron = G4Electron::Electron(); << 94 theElectron = G4Electron::Electron(); 95 thePositron = G4Positron::Positron(); << 95 thePositron = G4Positron::Positron(); 96 96 97 pParticleChange = &fParticleChange; 97 pParticleChange = &fParticleChange; 98 fParticleChange.SetSecondaryWeightByProcess( 98 fParticleChange.SetSecondaryWeightByProcess(true); 99 secParticles.reserve(5); 99 secParticles.reserve(5); 100 100 101 modelManager = new G4EmModelManager(); 101 modelManager = new G4EmModelManager(); 102 lManager = G4LossTableManager::Instance(); 102 lManager = G4LossTableManager::Instance(); 103 lManager->Register(this); 103 lManager->Register(this); 104 isTheMaster = lManager->IsMaster(); << 105 G4LossTableBuilder* bld = lManager->GetTable 104 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 106 theDensityFactor = bld->GetDensityFactors(); 105 theDensityFactor = bld->GetDensityFactors(); 107 theDensityIdx = bld->GetCoupleIndexes(); 106 theDensityIdx = bld->GetCoupleIndexes(); 108 } 107 } 109 108 110 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 111 110 112 G4VEmProcess::~G4VEmProcess() 111 G4VEmProcess::~G4VEmProcess() 113 { 112 { >> 113 /* >> 114 if(1 < verboseLevel) { >> 115 G4cout << "G4VEmProcess destruct " << GetProcessName() >> 116 << " " << this << " " << theLambdaTable <<G4endl; >> 117 } >> 118 */ 114 if(isTheMaster) { 119 if(isTheMaster) { 115 delete theData; 120 delete theData; 116 delete theEnergyOfCrossSectionMax; 121 delete theEnergyOfCrossSectionMax; 117 } 122 } 118 delete modelManager; 123 delete modelManager; 119 delete biasManager; 124 delete biasManager; 120 lManager->DeRegister(this); 125 lManager->DeRegister(this); 121 } 126 } 122 127 123 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 129 >> 130 void G4VEmProcess::Clear() >> 131 { >> 132 currentCouple = nullptr; >> 133 preStepLambda = 0.0; >> 134 } >> 135 >> 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 137 >> 138 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*, >> 139 const G4Material*) >> 140 { >> 141 return 0.0; >> 142 } >> 143 >> 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 145 125 void G4VEmProcess::AddEmModel(G4int order, G4V 146 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* ptr, 126 const G4Region* 147 const G4Region* region) 127 { 148 { 128 if(nullptr == ptr) { return; } 149 if(nullptr == ptr) { return; } 129 G4VEmFluctuationModel* fm = nullptr; 150 G4VEmFluctuationModel* fm = nullptr; 130 modelManager->AddEmModel(order, ptr, fm, reg 151 modelManager->AddEmModel(order, ptr, fm, region); 131 ptr->SetParticleChange(pParticleChange); 152 ptr->SetParticleChange(pParticleChange); 132 } 153 } 133 154 134 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 135 156 136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, 157 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, G4int) 137 { 158 { 138 if(nullptr == ptr) { return; } 159 if(nullptr == ptr) { return; } 139 if(!emModels.empty()) { 160 if(!emModels.empty()) { 140 for(auto & em : emModels) { if(em == ptr) 161 for(auto & em : emModels) { if(em == ptr) { return; } } 141 } 162 } 142 emModels.push_back(ptr); 163 emModels.push_back(ptr); 143 } 164 } 144 165 145 //....oooOO0OOooo........oooOO0OOooo........oo 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 146 167 147 void G4VEmProcess::PreparePhysicsTable(const G 168 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part) 148 { 169 { >> 170 isTheMaster = lManager->IsMaster(); 149 if(nullptr == particle) { SetParticle(&part) 171 if(nullptr == particle) { SetParticle(&part); } 150 172 151 if(part.GetParticleType() == "nucleus" && 173 if(part.GetParticleType() == "nucleus" && 152 part.GetParticleSubType() == "generic") { 174 part.GetParticleSubType() == "generic") { 153 175 154 G4String pname = part.GetParticleName(); 176 G4String pname = part.GetParticleName(); 155 if(pname != "deuteron" && pname != "triton 177 if(pname != "deuteron" && pname != "triton" && 156 pname != "He3" && pname != "alpha" && p << 178 pname != "alpha" && pname != "He3" && 157 pname != "helium" && pname != "hydrogen << 179 pname != "alpha+" && pname != "helium" && >> 180 pname != "hydrogen") { 158 181 159 particle = G4GenericIon::GenericIon(); 182 particle = G4GenericIon::GenericIon(); 160 isIon = true; 183 isIon = true; 161 } 184 } 162 } 185 } 163 if(particle != &part) { return; } << 164 186 165 lManager->PreparePhysicsTable(&part, this); << 187 if(1 < verboseLevel) { >> 188 G4cout << "G4VEmProcess::PreparePhysicsTable() for " >> 189 << GetProcessName() >> 190 << " and particle " << part.GetParticleName() >> 191 << " local particle " << particle->GetParticleName() >> 192 << G4endl; >> 193 } 166 194 167 // for new run << 195 if(particle != &part) { return; } 168 currentCouple = nullptr; << 169 preStepLambda = 0.0; << 170 fLambdaEnergy = 0.0; << 171 196 >> 197 lManager->PreparePhysicsTable(&part, this, isTheMaster); >> 198 >> 199 Clear(); 172 InitialiseProcess(particle); 200 InitialiseProcess(particle); 173 201 174 G4LossTableBuilder* bld = lManager->GetTable 202 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 175 const G4ProductionCutsTable* theCoupleTable= 203 const G4ProductionCutsTable* theCoupleTable= 176 G4ProductionCutsTable::GetProductionCutsTa 204 G4ProductionCutsTable::GetProductionCutsTable(); 177 theCutsGamma = theCoupleTable->GetEnergyC << 178 theCutsElectron = theCoupleTable->GetEnergyC << 179 theCutsPositron = theCoupleTable->GetEnergyC << 180 205 181 // initialisation of the process 206 // initialisation of the process 182 if(!actMinKinEnergy) { minKinEnergy = thePar 207 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); } 183 if(!actMaxKinEnergy) { maxKinEnergy = thePar 208 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); } 184 209 >> 210 if(isTheMaster) { >> 211 SetVerboseLevel(theParameters->Verbose()); >> 212 if(nullptr == theData) { theData = new G4EmDataHandler(2); } >> 213 if(fEmOnePeak == fXSType) { >> 214 if(nullptr == theEnergyOfCrossSectionMax) { >> 215 theEnergyOfCrossSectionMax = new std::vector<G4double>; >> 216 } >> 217 size_t n = theCoupleTable->GetTableSize(); >> 218 theEnergyOfCrossSectionMax->resize(n, DBL_MAX); >> 219 } >> 220 } else { >> 221 SetVerboseLevel(theParameters->WorkerVerbose()); >> 222 } 185 applyCuts = theParameters->ApplyCuts() 223 applyCuts = theParameters->ApplyCuts(); 186 lambdaFactor = theParameters->LambdaFacto 224 lambdaFactor = theParameters->LambdaFactor(); 187 invLambdaFactor = 1.0/lambdaFactor; << 225 logLambdaFactor = G4Log(lambdaFactor); 188 theParameters->DefineRegParamForEM(this); 226 theParameters->DefineRegParamForEM(this); 189 227 190 // integral option may be disabled 228 // integral option may be disabled 191 if(!theParameters->Integral()) { fXSType = f 229 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; } 192 230 193 // prepare tables 231 // prepare tables 194 if(isTheMaster) { << 232 if(buildLambdaTable && isTheMaster){ 195 if(nullptr == theData) { theData = new G4E << 233 theLambdaTable = theData->MakeTable(0); 196 << 234 bld->InitialiseBaseMaterials(theLambdaTable); 197 if(buildLambdaTable) { << 235 } 198 theLambdaTable = theData->MakeTable(0); << 236 // high energy table 199 bld->InitialiseBaseMaterials(theLambdaTa << 237 if(isTheMaster && minKinEnergyPrim < maxKinEnergy){ 200 } << 238 theLambdaTablePrim = theData->MakeTable(1); 201 // high energy table << 239 bld->InitialiseBaseMaterials(theLambdaTablePrim); 202 if(minKinEnergyPrim < maxKinEnergy) { << 203 theLambdaTablePrim = theData->MakeTable( << 204 bld->InitialiseBaseMaterials(theLambdaTa << 205 } << 206 } 240 } 207 // models << 208 baseMat = bld->GetBaseMaterialFlag(); 241 baseMat = bld->GetBaseMaterialFlag(); >> 242 >> 243 // initialisation of models 209 numberOfModels = modelManager->NumberOfModel 244 numberOfModels = modelManager->NumberOfModels(); 210 currentModel = modelManager->GetModel(0); << 245 for(G4int i=0; i<numberOfModels; ++i) { >> 246 G4VEmModel* mod = modelManager->GetModel(i); >> 247 if(nullptr == mod) { continue; } >> 248 if(nullptr == currentModel) { currentModel = mod; } >> 249 mod->SetPolarAngleLimit(theParameters->MscThetaLimit()); >> 250 mod->SetMasterThread(isTheMaster); >> 251 if(mod->HighEnergyLimit() > maxKinEnergy) { >> 252 mod->SetHighEnergyLimit(maxKinEnergy); >> 253 } >> 254 SetEmModel(mod); >> 255 mod->SetUseBaseMaterials(baseMat); >> 256 } >> 257 211 if(nullptr != lManager->AtomDeexcitation()) 258 if(nullptr != lManager->AtomDeexcitation()) { 212 modelManager->SetFluoFlag(true); 259 modelManager->SetFluoFlag(true); 213 } 260 } >> 261 fLambdaEnergy = 0.0; >> 262 >> 263 theCuts = >> 264 modelManager->Initialise(particle,secondaryParticle,1.0,verboseLevel); >> 265 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut); >> 266 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut); >> 267 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut); >> 268 214 // forced biasing 269 // forced biasing 215 if(nullptr != biasManager) { << 270 if(biasManager) { 216 biasManager->Initialise(part, GetProcessNa << 271 biasManager->Initialise(part,GetProcessName(),verboseLevel); 217 biasFlag = false; 272 biasFlag = false; 218 } 273 } 219 274 220 theCuts = << 275 // defined ID of secondary particles 221 G4EmTableUtil::PrepareEmProcess(this, part << 276 G4int stype = GetProcessSubType(); 222 modelManag << 277 if(stype == fAnnihilation) { 223 secID, tri << 278 secID = _Annihilation; 224 verboseLev << 279 tripletID = _TripletGamma; >> 280 } else if(stype == fGammaConversion) { >> 281 secID = _PairProduction; >> 282 mainSecondaries = 2; >> 283 } else if(stype == fPhotoElectricEffect) { >> 284 secID = _PhotoElectron; >> 285 } else if(stype == fComptonScattering) { >> 286 secID = _ComptonElectron; >> 287 } else if(stype >= fLowEnergyElastic) { >> 288 secID = fDNAUnknownModel; >> 289 } >> 290 if(1 < verboseLevel) { >> 291 G4cout << "### G4VEmProcess::PreparePhysicsTable() done for " >> 292 << GetProcessName() >> 293 << " and particle " << part.GetParticleName() >> 294 << " baseMat=" << baseMat << G4endl; >> 295 } 225 } 296 } 226 297 227 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 228 299 229 void G4VEmProcess::BuildPhysicsTable(const G4P 300 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 230 { 301 { 231 if(nullptr == masterProc) { 302 if(nullptr == masterProc) { 232 if(isTheMaster) { masterProc = this; } 303 if(isTheMaster) { masterProc = this; } 233 else { masterProc = static_cast<const G4VE 304 else { masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess());} 234 } 305 } 235 G4int nModels = modelManager->NumberOfModels << 306 236 G4bool isLocked = theParameters->IsPrintLock << 307 G4String num = part.GetParticleName(); 237 G4bool toBuild = (buildLambdaTable || minKin << 308 if(1 < verboseLevel) { 238 << 309 G4cout << "### G4VEmProcess::BuildPhysicsTable() for " 239 G4EmTableUtil::BuildEmProcess(this, masterPr << 310 << GetProcessName() 240 nModels, verbo << 311 << " and particle " << num 241 isLocked, toBu << 312 << " buildLambdaTable= " << buildLambdaTable >> 313 << " isTheMaster= " << isTheMaster >> 314 << " " << masterProc >> 315 << G4endl; >> 316 } >> 317 >> 318 if(particle == &part) { >> 319 >> 320 // worker initialisation >> 321 if(!isTheMaster) { >> 322 theLambdaTable = masterProc->LambdaTable(); >> 323 theLambdaTablePrim = masterProc->LambdaTablePrim(); >> 324 theEnergyOfCrossSectionMax = masterProc->EnergyOfCrossSectionMax(); >> 325 baseMat = masterProc->UseBaseMaterial(); >> 326 >> 327 // local initialisation of models >> 328 G4bool printing = true; >> 329 for(G4int i=0; i<numberOfModels; ++i) { >> 330 G4VEmModel* mod = GetModelByIndex(i, printing); >> 331 G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing); >> 332 //G4cout << i << ". " << mod << " " << mod0 << " " >> 333 // << particle->GetParticleName() << G4endl; >> 334 mod->SetUseBaseMaterials(baseMat); >> 335 mod->InitialiseLocal(particle, mod0); >> 336 } >> 337 // master thread >> 338 } else { >> 339 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) { >> 340 BuildLambdaTable(); >> 341 } >> 342 } >> 343 } >> 344 // protection against double printout >> 345 if(theParameters->IsPrintLocked()) { return; } >> 346 >> 347 // explicitly defined printout by particle name >> 348 if(1 < verboseLevel || >> 349 (0 < verboseLevel && (num == "gamma" || num == "e-" || >> 350 num == "e+" || num == "mu+" || >> 351 num == "mu-" || num == "proton"|| >> 352 num == "pi+" || num == "pi-" || >> 353 num == "kaon+" || num == "kaon-" || >> 354 num == "alpha" || num == "anti_proton" || >> 355 num == "GenericIon"|| num == "alpha++" || >> 356 num == "alpha+" || num == "helium" || >> 357 num == "hydrogen"))) >> 358 { >> 359 StreamInfo(G4cout, part); >> 360 } >> 361 >> 362 if(1 < verboseLevel) { >> 363 G4cout << "### G4VEmProcess::BuildPhysicsTable() done for " >> 364 << GetProcessName() >> 365 << " and particle " << num >> 366 << " baseMat=" << baseMat >> 367 << G4endl; >> 368 } 242 } 369 } 243 370 244 //....oooOO0OOooo........oooOO0OOooo........oo 371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 245 372 246 void G4VEmProcess::BuildLambdaTable() 373 void G4VEmProcess::BuildLambdaTable() 247 { 374 { >> 375 if(1 < verboseLevel) { >> 376 G4cout << "G4EmProcess::BuildLambdaTable() for process " >> 377 << GetProcessName() << " and particle " >> 378 << particle->GetParticleName() << " " << this >> 379 << G4endl; >> 380 } >> 381 >> 382 // Access to materials >> 383 const G4ProductionCutsTable* theCoupleTable= >> 384 G4ProductionCutsTable::GetProductionCutsTable(); >> 385 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 386 >> 387 G4LossTableBuilder* bld = lManager->GetTableBuilder(); >> 388 >> 389 G4PhysicsLogVector* aVector = nullptr; >> 390 G4PhysicsLogVector* aVectorPrim = nullptr; >> 391 G4PhysicsLogVector* bVectorPrim = nullptr; >> 392 248 G4double scale = theParameters->MaxKinEnergy 393 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy(); 249 G4int nbin = 394 G4int nbin = 250 theParameters->NumberOfBinsPerDecade()*G4l 395 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale)); >> 396 scale = G4Log(scale); 251 if(actBinning) { nbin = std::max(nbin, nLamb 397 if(actBinning) { nbin = std::max(nbin, nLambdaBins); } 252 scale = nbin/G4Log(scale); << 398 G4double emax1 = std::min(maxKinEnergy, minKinEnergyPrim); 253 << 399 254 G4LossTableBuilder* bld = lManager->GetTable << 400 for(size_t i=0; i<numOfCouples; ++i) { 255 G4EmTableUtil::BuildLambdaTable(this, partic << 401 256 bld, theLamb << 402 if (bld->GetFlag(i)) { 257 minKinEnergy << 403 258 maxKinEnergy << 404 // create physics vector and fill it 259 startFromNul << 405 const G4MaterialCutsCouple* couple = >> 406 theCoupleTable->GetMaterialCutsCouple(i); >> 407 >> 408 // build main table >> 409 if(buildLambdaTable) { >> 410 delete (*theLambdaTable)[i]; >> 411 >> 412 // if start from zero then change the scale >> 413 G4double emin = minKinEnergy; >> 414 G4bool startNull = false; >> 415 if(startFromNull) { >> 416 G4double e = MinPrimaryEnergy(particle,couple->GetMaterial()); >> 417 if(e >= emin) { >> 418 emin = e; >> 419 startNull = true; >> 420 } >> 421 } >> 422 G4double emax = emax1; >> 423 if(emax <= emin) { emax = 2*emin; } >> 424 G4int bin = G4lrint(nbin*G4Log(emax/emin)/scale); >> 425 if(bin < 3) { bin = 3; } >> 426 aVector = new G4PhysicsLogVector(emin, emax, bin, splineFlag); >> 427 modelManager->FillLambdaVector(aVector, couple, startNull); >> 428 if(splineFlag) { aVector->FillSecondDerivatives(); } >> 429 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector); >> 430 } >> 431 // build high energy table >> 432 if(minKinEnergyPrim < maxKinEnergy) { >> 433 delete (*theLambdaTablePrim)[i]; >> 434 >> 435 // start not from zero and always use spline >> 436 if(!bVectorPrim) { >> 437 G4int bin = G4lrint(nbin*G4Log(maxKinEnergy/minKinEnergyPrim)/scale); >> 438 if(bin < 3) { bin = 3; } >> 439 aVectorPrim = >> 440 new G4PhysicsLogVector(minKinEnergyPrim, maxKinEnergy, bin, true); >> 441 bVectorPrim = aVectorPrim; >> 442 } else { >> 443 aVectorPrim = new G4PhysicsLogVector(*bVectorPrim); >> 444 } >> 445 modelManager->FillLambdaVector(aVectorPrim, couple, false, >> 446 fIsCrossSectionPrim); >> 447 aVectorPrim->FillSecondDerivatives(); >> 448 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTablePrim, i, >> 449 aVectorPrim); >> 450 } >> 451 } >> 452 } >> 453 >> 454 if(buildLambdaTable && fXSType == fEmOnePeak) { FindLambdaMax(); } >> 455 >> 456 if(1 < verboseLevel) { >> 457 G4cout << "Lambda table is built for " >> 458 << particle->GetParticleName() >> 459 << G4endl; >> 460 } 260 } 461 } 261 462 262 //....oooOO0OOooo........oooOO0OOooo........oo 463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 263 464 264 void G4VEmProcess::StreamInfo(std::ostream& ou 465 void G4VEmProcess::StreamInfo(std::ostream& out, 265 const G4ParticleDefinition& 466 const G4ParticleDefinition& part, G4bool rst) const 266 { 467 { 267 G4String indent = (rst ? " " : ""); 468 G4String indent = (rst ? " " : ""); 268 out << std::setprecision(6); 469 out << std::setprecision(6); 269 out << G4endl << indent << GetProcessName() 470 out << G4endl << indent << GetProcessName() << ": "; 270 if (!rst) { 471 if (!rst) { 271 out << " for " << part.GetParticleName(); 472 out << " for " << part.GetParticleName(); 272 } 473 } 273 if(fXSType != fEmNoIntegral) { out << " XSt 474 if(fXSType != fEmNoIntegral) { out << " XStype:" << fXSType; } 274 if(applyCuts) { out << " applyCuts:1 "; } 475 if(applyCuts) { out << " applyCuts:1 "; } 275 G4int subtype = GetProcessSubType(); << 476 out << " SubType=" << GetProcessSubType(); 276 out << " SubType=" << subtype; << 477 if(biasFactor != 1.0) { out << " BiasingFactor= " << biasFactor; } 277 if (subtype == fAnnihilation) { << 278 G4int mod = theParameters->PositronAtRestM << 279 const G4String namp[2] = {"Simple", "Allis << 280 out << " AtRestModel:" << namp[mod]; << 281 } << 282 if(biasFactor != 1.0) { out << " BiasingFac << 283 out << " BuildTable=" << buildLambdaTable << 478 out << " BuildTable=" << buildLambdaTable << G4endl; 284 if(buildLambdaTable) { 479 if(buildLambdaTable) { 285 if(particle == &part) { 480 if(particle == &part) { 286 for(auto & v : *theLambdaTable) { << 481 size_t length = theLambdaTable->length(); 287 if(nullptr != v) { << 482 for(size_t i=0; i<length; ++i) { >> 483 G4PhysicsVector* v = (*theLambdaTable)[i]; >> 484 if(v) { 288 out << " Lambda table from "; 485 out << " Lambda table from "; 289 G4double emin = v->Energy(0); 486 G4double emin = v->Energy(0); 290 G4double emax = v->GetMaxEnergy(); 487 G4double emax = v->GetMaxEnergy(); 291 G4int nbin = G4int(v->GetVectorLengt << 488 G4int nbin = v->GetVectorLength() - 1; 292 if(emin > minKinEnergy) { out << "th 489 if(emin > minKinEnergy) { out << "threshold "; } 293 else { out << G4BestUnit(emin,"Energ 490 else { out << G4BestUnit(emin,"Energy"); } 294 out << " to " 491 out << " to " 295 << G4BestUnit(emax,"Energy") 492 << G4BestUnit(emax,"Energy") 296 << ", " << G4lrint(nbin/std::log 493 << ", " << G4lrint(nbin/std::log10(emax/emin)) 297 << " bins/decade, spline: " 494 << " bins/decade, spline: " 298 << splineFlag << G4endl; 495 << splineFlag << G4endl; 299 break; 496 break; 300 } 497 } 301 } 498 } 302 } else { 499 } else { 303 out << " Used Lambda table of " 500 out << " Used Lambda table of " 304 << particle->GetParticleName() << G4endl 501 << particle->GetParticleName() << G4endl; 305 } 502 } 306 } 503 } 307 if(minKinEnergyPrim < maxKinEnergy) { 504 if(minKinEnergyPrim < maxKinEnergy) { 308 if(particle == &part) { << 505 if(particle == &part) { 309 for(auto & v : *theLambdaTablePrim) { << 506 size_t length = theLambdaTablePrim->length(); 310 if(nullptr != v) { << 507 for(size_t i=0; i<length; ++i) { >> 508 G4PhysicsVector* v = (*theLambdaTablePrim)[i]; >> 509 if(v) { 311 out << " LambdaPrime table from 510 out << " LambdaPrime table from " 312 << G4BestUnit(v->Energy(0),"Ener 511 << G4BestUnit(v->Energy(0),"Energy") 313 << " to " 512 << " to " 314 << G4BestUnit(v->GetMaxEnergy(), 513 << G4BestUnit(v->GetMaxEnergy(),"Energy") 315 << " in " << v->GetVectorLength( 514 << " in " << v->GetVectorLength()-1 316 << " bins " << G4endl; 515 << " bins " << G4endl; 317 break; 516 break; 318 } 517 } 319 } 518 } 320 } else { 519 } else { 321 out << " Used LambdaPrime table of 520 out << " Used LambdaPrime table of " 322 << particle->GetParticleName() 521 << particle->GetParticleName() << G4endl; 323 } 522 } 324 } 523 } 325 StreamProcessInfo(out); 524 StreamProcessInfo(out); 326 modelManager->DumpModelList(out, verboseLeve 525 modelManager->DumpModelList(out, verboseLevel); 327 526 328 if(verboseLevel > 2 && buildLambdaTable) { 527 if(verboseLevel > 2 && buildLambdaTable) { 329 out << " LambdaTable address= " << th 528 out << " LambdaTable address= " << theLambdaTable << G4endl; 330 if(theLambdaTable && particle == &part) { 529 if(theLambdaTable && particle == &part) { 331 out << (*theLambdaTable) << G4endl; 530 out << (*theLambdaTable) << G4endl; 332 } 531 } 333 } 532 } 334 } 533 } 335 534 336 //....oooOO0OOooo........oooOO0OOooo........oo 535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 337 536 338 void G4VEmProcess::StartTracking(G4Track* trac 537 void G4VEmProcess::StartTracking(G4Track* track) 339 { 538 { 340 // reset parameters for the new track 539 // reset parameters for the new track 341 currentParticle = track->GetParticleDefiniti 540 currentParticle = track->GetParticleDefinition(); 342 theNumberOfInteractionLengthLeft = -1.0; 541 theNumberOfInteractionLengthLeft = -1.0; 343 mfpKinEnergy = DBL_MAX; 542 mfpKinEnergy = DBL_MAX; 344 preStepLambda = 0.0; << 345 543 346 if(isIon) { massRatio = proton_mass_c2/curre 544 if(isIon) { massRatio = proton_mass_c2/currentParticle->GetPDGMass(); } 347 545 348 // forced biasing only for primary particles 546 // forced biasing only for primary particles 349 if(biasManager) { 547 if(biasManager) { 350 if(0 == track->GetParentID()) { 548 if(0 == track->GetParentID()) { 351 // primary particle 549 // primary particle 352 biasFlag = true; 550 biasFlag = true; 353 biasManager->ResetForcedInteraction(); 551 biasManager->ResetForcedInteraction(); 354 } 552 } 355 } 553 } 356 } 554 } 357 555 358 //....oooOO0OOooo........oooOO0OOooo........oo 556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 359 557 360 G4double G4VEmProcess::PostStepGetPhysicalInte 558 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength( 361 const G4Track& tr 559 const G4Track& track, 362 G4double previo 560 G4double previousStepSize, 363 G4ForceCondition* 561 G4ForceCondition* condition) 364 { 562 { 365 *condition = NotForced; 563 *condition = NotForced; 366 G4double x = DBL_MAX; 564 G4double x = DBL_MAX; 367 565 368 DefineMaterial(track.GetMaterialCutsCouple() 566 DefineMaterial(track.GetMaterialCutsCouple()); 369 preStepKinEnergy = track.GetKineticEnergy(); << 567 preStepKinEnergy = track.GetKineticEnergy(); >> 568 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy(); 370 const G4double scaledEnergy = preStepKinEner 569 const G4double scaledEnergy = preStepKinEnergy*massRatio; 371 SelectModel(scaledEnergy, currentCoupleIndex 570 SelectModel(scaledEnergy, currentCoupleIndex); 372 /* 571 /* 373 G4cout << "PostStepGetPhysicalInteractionLen 572 G4cout << "PostStepGetPhysicalInteractionLength: idx= " << currentCoupleIndex 374 << " couple: " << currentCouple << G 573 << " couple: " << currentCouple << G4endl; 375 */ 574 */ 376 if(!currentModel->IsActive(scaledEnergy)) { 575 if(!currentModel->IsActive(scaledEnergy)) { 377 theNumberOfInteractionLengthLeft = -1.0; 576 theNumberOfInteractionLengthLeft = -1.0; 378 currentInteractionLength = DBL_MAX; 577 currentInteractionLength = DBL_MAX; 379 mfpKinEnergy = DBL_MAX; << 380 preStepLambda = 0.0; << 381 return x; 578 return x; 382 } 579 } 383 580 384 // forced biasing only for primary particles 581 // forced biasing only for primary particles 385 if(biasManager) { 582 if(biasManager) { 386 if(0 == track.GetParentID()) { 583 if(0 == track.GetParentID()) { 387 if(biasFlag && 584 if(biasFlag && 388 biasManager->ForcedInteractionRegion( << 585 biasManager->ForcedInteractionRegion(currentCoupleIndex)) { 389 return biasManager->GetStepLimit((G4in << 586 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize); 390 } 587 } 391 } 588 } 392 } 589 } 393 590 394 // compute mean free path 591 // compute mean free path 395 << 592 ComputeIntegralLambda(preStepKinEnergy, preStepLogKinEnergy); 396 ComputeIntegralLambda(preStepKinEnergy, trac << 397 593 398 // zero cross section 594 // zero cross section 399 if(preStepLambda <= 0.0) { 595 if(preStepLambda <= 0.0) { 400 theNumberOfInteractionLengthLeft = -1.0; 596 theNumberOfInteractionLengthLeft = -1.0; 401 currentInteractionLength = DBL_MAX; 597 currentInteractionLength = DBL_MAX; 402 598 403 } else { 599 } else { 404 600 405 // non-zero cross section 601 // non-zero cross section 406 if (theNumberOfInteractionLengthLeft < 0.0 602 if (theNumberOfInteractionLengthLeft < 0.0) { 407 603 408 // beggining of tracking (or just after 604 // beggining of tracking (or just after DoIt of this process) 409 theNumberOfInteractionLengthLeft = -G4Lo 605 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() ); 410 theInitialNumberOfInteractionLength = th 606 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 411 607 412 } else { << 608 } else if(currentInteractionLength < DBL_MAX) { 413 609 414 theNumberOfInteractionLengthLeft -= 610 theNumberOfInteractionLengthLeft -= 415 previousStepSize/currentInteractionLen 611 previousStepSize/currentInteractionLength; 416 theNumberOfInteractionLengthLeft = 612 theNumberOfInteractionLengthLeft = 417 std::max(theNumberOfInteractionLengthL 613 std::max(theNumberOfInteractionLengthLeft, 0.0); 418 } 614 } 419 615 420 // new mean free path and step limit for t 616 // new mean free path and step limit for the next step 421 currentInteractionLength = 1.0/preStepLamb 617 currentInteractionLength = 1.0/preStepLambda; 422 x = theNumberOfInteractionLengthLeft * cur 618 x = theNumberOfInteractionLengthLeft * currentInteractionLength; 423 } 619 } 424 return x; 620 return x; 425 } 621 } 426 622 427 //....oooOO0OOooo........oooOO0OOooo........oo 623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 428 624 429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 625 void G4VEmProcess::ComputeIntegralLambda(G4double e, G4double loge) 430 { 626 { 431 if (fXSType == fEmNoIntegral) { << 627 if(fXSType == fEmNoIntegral) { 432 preStepLambda = GetCurrentLambda(e, LogEki << 628 preStepLambda = GetCurrentLambda(e, loge); 433 629 434 } else if (fXSType == fEmIncreasing) { << 630 } else if(fXSType == fEmIncreasing) { 435 if(e*invLambdaFactor < mfpKinEnergy) { << 631 if(e/lambdaFactor < mfpKinEnergy) { 436 preStepLambda = GetCurrentLambda(e, LogE << 632 mfpKinEnergy = e; 437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 633 preStepLambda = GetCurrentLambda(e, loge); 438 } 634 } 439 635 440 } else if(fXSType == fEmDecreasing) { 636 } else if(fXSType == fEmDecreasing) { 441 if(e < mfpKinEnergy) { 637 if(e < mfpKinEnergy) { 442 const G4double e1 = e*lambdaFactor; 638 const G4double e1 = e*lambdaFactor; 443 preStepLambda = GetCurrentLambda(e1); 639 preStepLambda = GetCurrentLambda(e1); 444 mfpKinEnergy = e1; 640 mfpKinEnergy = e1; 445 } 641 } 446 642 447 } else if(fXSType == fEmOnePeak) { 643 } else if(fXSType == fEmOnePeak) { 448 const G4double epeak = (*theEnergyOfCrossS 644 const G4double epeak = (*theEnergyOfCrossSectionMax)[currentCoupleIndex]; 449 if(e <= epeak) { 645 if(e <= epeak) { 450 if(e*invLambdaFactor < mfpKinEnergy) { << 646 if(e/lambdaFactor < mfpKinEnergy) { 451 preStepLambda = GetCurrentLambda(e, Lo << 647 mfpKinEnergy = e; 452 mfpKinEnergy = (preStepLambda > 0.0) ? << 648 preStepLambda = GetCurrentLambda(e, loge); 453 } 649 } 454 } else if(e < mfpKinEnergy) { 650 } else if(e < mfpKinEnergy) { 455 const G4double e1 = std::max(epeak, e*la 651 const G4double e1 = std::max(epeak, e*lambdaFactor); 456 preStepLambda = GetCurrentLambda(e1); << 652 preStepLambda = GetCurrentLambda(e1); 457 mfpKinEnergy = e1; 653 mfpKinEnergy = e1; 458 } 654 } >> 655 459 } else { 656 } else { 460 preStepLambda = GetCurrentLambda(e, LogEki << 657 preStepLambda = GetCurrentLambda(e, loge); 461 } 658 } 462 } 659 } 463 660 464 //....oooOO0OOooo........oooOO0OOooo........oo 661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 465 662 466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 663 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track, 467 664 const G4Step& step) 468 { 665 { 469 // clear number of interaction lengths in an << 666 // In all cases clear number of interaction lengths 470 theNumberOfInteractionLengthLeft = -1.0; 667 theNumberOfInteractionLengthLeft = -1.0; 471 mfpKinEnergy = DBL_MAX; 668 mfpKinEnergy = DBL_MAX; 472 669 473 fParticleChange.InitializeForPostStep(track) 670 fParticleChange.InitializeForPostStep(track); 474 671 475 // Do not make anything if particle is stopp 672 // Do not make anything if particle is stopped, the annihilation then 476 // should be performed by the AtRestDoIt! 673 // should be performed by the AtRestDoIt! 477 if (track.GetTrackStatus() == fStopButAlive) 674 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; } 478 675 479 const G4double finalT = track.GetKineticEner << 676 const G4double finalT = track.GetKineticEnergy(); 480 677 481 // forced process - should happen only once 678 // forced process - should happen only once per track 482 if(biasFlag) { 679 if(biasFlag) { 483 if(biasManager->ForcedInteractionRegion((G << 680 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) { 484 biasFlag = false; 681 biasFlag = false; 485 } 682 } 486 } 683 } 487 684 488 // check active and select model 685 // check active and select model 489 const G4double scaledEnergy = finalT*massRat 686 const G4double scaledEnergy = finalT*massRatio; 490 SelectModel(scaledEnergy, currentCoupleIndex 687 SelectModel(scaledEnergy, currentCoupleIndex); 491 if(!currentModel->IsActive(scaledEnergy)) { 688 if(!currentModel->IsActive(scaledEnergy)) { return &fParticleChange; } 492 689 493 // Integral approach 690 // Integral approach 494 if (fXSType != fEmNoIntegral) { 691 if (fXSType != fEmNoIntegral) { 495 const G4double logFinalT = << 692 const G4double logFinalT = track.GetDynamicParticle()->GetLogKineticEnergy(); 496 track.GetDynamicParticle()->GetLogKineti << 497 const G4double lx = std::max(GetCurrentLam 693 const G4double lx = std::max(GetCurrentLambda(finalT, logFinalT), 0.0); >> 694 const G4double lg = preStepLambda; >> 695 if(finalT < mfpKinEnergy) { >> 696 mfpKinEnergy = finalT; >> 697 preStepLambda = lx; >> 698 } 498 #ifdef G4VERBOSE 699 #ifdef G4VERBOSE 499 if(preStepLambda < lx && 1 < verboseLevel) << 700 if(lg < lx && 1 < verboseLevel) { 500 G4cout << "WARNING: for " << currentPart 701 G4cout << "WARNING: for " << currentParticle->GetParticleName() 501 << " and " << GetProcessName() << << 702 << " and " << GetProcessName() 502 << " preLambda= " << preStepLambd << 703 << " E(MeV)= " << finalT/MeV 503 << " < " << lx << " (postLambda) << 704 << " preLambda= " << lg << " < " << lx << " (postLambda) " >> 705 << G4endl; 504 } 706 } 505 #endif 707 #endif 506 // if false interaction then use new cross << 708 if(lg*G4UniformRand() >= lx) { 507 // if both values are zero - no interactio << 508 if(preStepLambda*G4UniformRand() >= lx) { << 509 return &fParticleChange; 709 return &fParticleChange; 510 } 710 } 511 } 711 } 512 712 513 // define new weight for primary and seconda 713 // define new weight for primary and secondaries 514 G4double weight = fParticleChange.GetParentW 714 G4double weight = fParticleChange.GetParentWeight(); 515 if(weightFlag) { 715 if(weightFlag) { 516 weight /= biasFactor; 716 weight /= biasFactor; 517 fParticleChange.ProposeWeight(weight); 717 fParticleChange.ProposeWeight(weight); 518 } 718 } 519 719 520 #ifdef G4VERBOSE 720 #ifdef G4VERBOSE 521 if(1 < verboseLevel) { 721 if(1 < verboseLevel) { 522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 722 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= " 523 << finalT/MeV 723 << finalT/MeV 524 << " MeV; model= (" << currentModel 724 << " MeV; model= (" << currentModel->LowEnergyLimit() 525 << ", " << currentModel->HighEnerg 725 << ", " << currentModel->HighEnergyLimit() << ")" 526 << G4endl; 726 << G4endl; 527 } 727 } 528 #endif 728 #endif 529 729 530 // sample secondaries 730 // sample secondaries 531 secParticles.clear(); 731 secParticles.clear(); 532 currentModel->SampleSecondaries(&secParticle 732 currentModel->SampleSecondaries(&secParticles, 533 currentCoupl 733 currentCouple, 534 track.GetDyn 734 track.GetDynamicParticle(), 535 (*theCuts)[c 735 (*theCuts)[currentCoupleIndex]); 536 736 537 G4int num0 = (G4int)secParticles.size(); << 737 G4int num0 = secParticles.size(); 538 738 539 // splitting or Russian roulette 739 // splitting or Russian roulette 540 if(biasManager) { 740 if(biasManager) { 541 if(biasManager->SecondaryBiasingRegion((G4 << 741 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) { 542 G4double eloss = 0.0; 742 G4double eloss = 0.0; 543 weight *= biasManager->ApplySecondaryBia 743 weight *= biasManager->ApplySecondaryBiasing( 544 secParticles, track, currentModel, &fP 744 secParticles, track, currentModel, &fParticleChange, eloss, 545 (G4int)currentCoupleIndex, (*theCuts)[ << 745 currentCoupleIndex, (*theCuts)[currentCoupleIndex], 546 step.GetPostStepPoint()->GetSafety()); 746 step.GetPostStepPoint()->GetSafety()); 547 if(eloss > 0.0) { 747 if(eloss > 0.0) { 548 eloss += fParticleChange.GetLocalEnerg 748 eloss += fParticleChange.GetLocalEnergyDeposit(); 549 fParticleChange.ProposeLocalEnergyDepo 749 fParticleChange.ProposeLocalEnergyDeposit(eloss); 550 } 750 } 551 } 751 } 552 } 752 } 553 753 554 // save secondaries 754 // save secondaries 555 G4int num = (G4int)secParticles.size(); << 755 G4int num = secParticles.size(); 556 if(num > 0) { 756 if(num > 0) { 557 757 558 fParticleChange.SetNumberOfSecondaries(num 758 fParticleChange.SetNumberOfSecondaries(num); 559 G4double edep = fParticleChange.GetLocalEn 759 G4double edep = fParticleChange.GetLocalEnergyDeposit(); 560 G4double time = track.GetGlobalTime(); 760 G4double time = track.GetGlobalTime(); 561 761 562 G4int n1(0), n2(0); 762 G4int n1(0), n2(0); 563 if(num0 > mainSecondaries) { << 763 if(num > mainSecondaries) { 564 currentModel->FillNumberOfSecondaries(n1 764 currentModel->FillNumberOfSecondaries(n1, n2); 565 } 765 } 566 766 567 for (G4int i=0; i<num; ++i) { 767 for (G4int i=0; i<num; ++i) { 568 G4DynamicParticle* dp = secParticles[i]; 768 G4DynamicParticle* dp = secParticles[i]; 569 if (nullptr != dp) { 769 if (nullptr != dp) { 570 const G4ParticleDefinition* p = dp->Ge 770 const G4ParticleDefinition* p = dp->GetParticleDefinition(); 571 G4double e = dp->GetKineticEnergy(); 771 G4double e = dp->GetKineticEnergy(); 572 G4bool good = true; 772 G4bool good = true; 573 if(applyCuts) { 773 if(applyCuts) { 574 if (p == theGamma) { 774 if (p == theGamma) { 575 if (e < (*theCutsGamma)[currentCou 775 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; } 576 776 577 } else if (p == theElectron) { 777 } else if (p == theElectron) { 578 if (e < (*theCutsElectron)[current 778 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; } 579 779 580 } else if (p == thePositron) { 780 } else if (p == thePositron) { 581 if (electron_mass_c2 < (*theCutsGa 781 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] && 582 e < (*theCutsPositron)[current 782 e < (*theCutsPositron)[currentCoupleIndex]) { 583 good = false; 783 good = false; 584 e += 2.0*electron_mass_c2; 784 e += 2.0*electron_mass_c2; 585 } 785 } 586 } 786 } 587 // added secondary if it is good 787 // added secondary if it is good 588 } 788 } 589 if (good) { 789 if (good) { 590 G4Track* t = new G4Track(dp, time, t 790 G4Track* t = new G4Track(dp, time, track.GetPosition()); 591 t->SetTouchableHandle(track.GetTouch 791 t->SetTouchableHandle(track.GetTouchableHandle()); 592 if (biasManager) { 792 if (biasManager) { 593 t->SetWeight(weight * biasManager- 793 t->SetWeight(weight * biasManager->GetWeight(i)); 594 } else { 794 } else { 595 t->SetWeight(weight); 795 t->SetWeight(weight); 596 } 796 } 597 pParticleChange->AddSecondary(t); 797 pParticleChange->AddSecondary(t); 598 798 599 // define type of secondary 799 // define type of secondary 600 if(i < mainSecondaries) { 800 if(i < mainSecondaries) { 601 t->SetCreatorModelID(secID); 801 t->SetCreatorModelID(secID); 602 if(GetProcessSubType() == fCompton 802 if(GetProcessSubType() == fComptonScattering && p == theGamma) { 603 t->SetCreatorModelID(_ComptonGam 803 t->SetCreatorModelID(_ComptonGamma); 604 } 804 } 605 } else if(i < mainSecondaries + n1) 805 } else if(i < mainSecondaries + n1) { 606 t->SetCreatorModelID(tripletID); 806 t->SetCreatorModelID(tripletID); 607 } else if(i < mainSecondaries + n1 + 807 } else if(i < mainSecondaries + n1 + n2) { 608 t->SetCreatorModelID(_IonRecoil); 808 t->SetCreatorModelID(_IonRecoil); 609 } else { 809 } else { 610 if(i < num0) { 810 if(i < num0) { 611 if(p == theGamma) { 811 if(p == theGamma) { 612 t->SetCreatorModelID(fluoID); 812 t->SetCreatorModelID(fluoID); 613 } else { 813 } else { 614 t->SetCreatorModelID(augerID); 814 t->SetCreatorModelID(augerID); 615 } 815 } 616 } else { 816 } else { 617 t->SetCreatorModelID(biasID); << 817 t->SetCreatorModelID(secID); 618 } 818 } 619 } 819 } 620 /* 820 /* 621 G4cout << "Secondary(post step) has 821 G4cout << "Secondary(post step) has weight " << t->GetWeight() 622 << ", Ekin= " << t->GetKineti 822 << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV " 623 << GetProcessName() << " fluo 823 << GetProcessName() << " fluoID= " << fluoID 624 << " augerID= " << augerID << 824 << " augerID= " << augerID <<G4endl; 625 */ 825 */ 626 } else { 826 } else { 627 delete dp; 827 delete dp; 628 edep += e; 828 edep += e; 629 } 829 } 630 } 830 } 631 } 831 } 632 fParticleChange.ProposeLocalEnergyDeposit( 832 fParticleChange.ProposeLocalEnergyDeposit(edep); 633 } 833 } 634 834 635 if(0.0 == fParticleChange.GetProposedKinetic 835 if(0.0 == fParticleChange.GetProposedKineticEnergy() && 636 fAlive == fParticleChange.GetTrackStatus( 836 fAlive == fParticleChange.GetTrackStatus()) { 637 if(particle->GetProcessManager()->GetAtRes 837 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0) 638 { fParticleChange.ProposeTrackStatus( 838 { fParticleChange.ProposeTrackStatus(fStopButAlive); } 639 else { fParticleChange.ProposeTrackStatus( 839 else { fParticleChange.ProposeTrackStatus(fStopAndKill); } 640 } 840 } 641 841 642 return &fParticleChange; 842 return &fParticleChange; 643 } 843 } 644 844 645 //....oooOO0OOooo........oooOO0OOooo........oo 845 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 646 846 647 G4bool G4VEmProcess::StorePhysicsTable(const G 847 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part, 648 const G 848 const G4String& directory, 649 G4bool 849 G4bool ascii) 650 { 850 { 651 if(!isTheMaster || part != particle) { retur << 851 G4bool yes = true; 652 if(G4EmTableUtil::StoreTable(this, part, the << 852 if(!isTheMaster) { return yes; } 653 directory, "Lambda", << 853 654 verboseLevel, a << 854 if ( theLambdaTable && part == particle) { 655 G4EmTableUtil::StoreTable(this, part, the << 855 const G4String& nam = 656 directory, "LambdaPrim", << 856 GetPhysicsTableFileName(part,directory,"Lambda",ascii); 657 verboseLevel, a << 857 yes = theLambdaTable->StorePhysicsTable(nam,ascii); 658 return true; << 858 >> 859 if ( yes ) { >> 860 if(0 < verboseLevel) G4cout << "Stored: " << nam << G4endl; >> 861 } else { >> 862 G4cout << "Fail to store Physics Table for " >> 863 << particle->GetParticleName() >> 864 << " and process " << GetProcessName() >> 865 << " in the directory <" << directory >> 866 << "> " << G4endl; >> 867 } 659 } 868 } 660 return false; << 869 if ( theLambdaTablePrim && part == particle) { >> 870 const G4String& name = >> 871 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii); >> 872 yes = theLambdaTablePrim->StorePhysicsTable(name,ascii); >> 873 >> 874 if ( yes ) { >> 875 if(0 < verboseLevel) { >> 876 G4cout << "Physics table prim is stored for " >> 877 << particle->GetParticleName() >> 878 << " and process " << GetProcessName() >> 879 << " in the directory <" << directory >> 880 << "> " << G4endl; >> 881 } >> 882 } else { >> 883 G4cout << "Fail to store Physics Table Prim for " >> 884 << particle->GetParticleName() >> 885 << " and process " << GetProcessName() >> 886 << " in the directory <" << directory >> 887 << "> " << G4endl; >> 888 } >> 889 } >> 890 return yes; 661 } 891 } 662 892 663 //....oooOO0OOooo........oooOO0OOooo........oo 893 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 664 894 665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 895 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 666 cons << 896 const G4String& directory, 667 G4bo 897 G4bool ascii) 668 { 898 { 669 if(!isTheMaster || part != particle) { retur << 899 if(1 < verboseLevel) { >> 900 G4cout << "G4VEmProcess::RetrievePhysicsTable() for " >> 901 << part->GetParticleName() << " and process " >> 902 << GetProcessName() << G4endl; >> 903 } 670 G4bool yes = true; 904 G4bool yes = true; >> 905 >> 906 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy) >> 907 || particle != part) { return yes; } >> 908 >> 909 const G4String particleName = part->GetParticleName(); >> 910 671 if(buildLambdaTable) { 911 if(buildLambdaTable) { 672 yes = G4EmTableUtil::RetrieveTable(this, p << 912 const G4String& filename = 673 "Lambda << 913 GetPhysicsTableFileName(part,directory,"Lambda",ascii); 674 ascii, << 914 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable, 675 } << 915 filename,ascii, 676 if(yes && minKinEnergyPrim < maxKinEnergy) { << 916 splineFlag); 677 yes = G4EmTableUtil::RetrieveTable(this, p << 917 if ( yes ) { 678 "Lambda << 918 if (0 < verboseLevel) { 679 ascii, << 919 G4cout << "Lambda table for " << particleName >> 920 << " is Retrieved from <" >> 921 << filename << ">" >> 922 << G4endl; >> 923 } >> 924 if(splineFlag) { >> 925 for(auto & v : *theLambdaTable) { >> 926 if(nullptr != v) { v->FillSecondDerivatives(); } >> 927 } >> 928 } >> 929 >> 930 } else { >> 931 if (1 < verboseLevel) { >> 932 G4cout << "Lambda table for " << particleName << " in file <" >> 933 << filename << "> is not exist" >> 934 << G4endl; >> 935 } >> 936 } >> 937 } >> 938 if(minKinEnergyPrim < maxKinEnergy) { >> 939 const G4String& filename = >> 940 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii); >> 941 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTablePrim, >> 942 filename,ascii,true); >> 943 if ( yes ) { >> 944 if (0 < verboseLevel) { >> 945 G4cout << "Lambda table prim for " << particleName >> 946 << " is Retrieved from <" >> 947 << filename << ">" >> 948 << G4endl; >> 949 } >> 950 for(auto & v : *theLambdaTablePrim) { >> 951 if(nullptr != v) { v->FillSecondDerivatives(); } >> 952 } >> 953 } else { >> 954 if (1 < verboseLevel) { >> 955 G4cout << "Lambda table prim for " << particleName << " in file <" >> 956 << filename << "> is not exist" >> 957 << G4endl; >> 958 } >> 959 } 680 } 960 } 681 return yes; 961 return yes; 682 } 962 } 683 963 684 //....oooOO0OOooo........oooOO0OOooo........oo 964 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 685 965 686 G4double G4VEmProcess::GetCrossSection(G4doubl << 966 G4double 687 const G << 967 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy, >> 968 const G4MaterialCutsCouple* couple, >> 969 G4double logKinEnergy) 688 { 970 { 689 CurrentSetup(couple, kinEnergy); << 971 // Cross section per atom is calculated 690 return GetCurrentLambda(kinEnergy, G4Log(kin << 972 DefineMaterial(couple); >> 973 G4double cross = 0.0; >> 974 if(buildLambdaTable) { >> 975 cross = GetCurrentLambda(kineticEnergy, >> 976 (logKinEnergy < DBL_MAX) ? logKinEnergy : G4Log(kineticEnergy)); >> 977 } else { >> 978 SelectModel(kineticEnergy, currentCoupleIndex); >> 979 if(currentModel) { >> 980 cross = fFactor*currentModel->CrossSectionPerVolume(currentMaterial, >> 981 currentParticle, >> 982 kineticEnergy); >> 983 } >> 984 } >> 985 return std::max(cross, 0.0); 691 } 986 } 692 987 693 //....oooOO0OOooo........oooOO0OOooo........oo 988 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 694 989 695 G4double G4VEmProcess::GetMeanFreePath(const G 990 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track, 696 G4doubl 991 G4double, 697 G4Force 992 G4ForceCondition* condition) 698 { 993 { 699 *condition = NotForced; 994 *condition = NotForced; 700 return G4VEmProcess::MeanFreePath(track); 995 return G4VEmProcess::MeanFreePath(track); 701 } 996 } 702 997 703 //....oooOO0OOooo........oooOO0OOooo........oo 998 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 704 999 >> 1000 G4double G4VEmProcess::MeanFreePath(const G4Track& track) >> 1001 { >> 1002 const G4double kinEnergy = track.GetKineticEnergy(); >> 1003 CurrentSetup(track.GetMaterialCutsCouple(), kinEnergy); >> 1004 const G4double xs = GetCurrentLambda(kinEnergy, >> 1005 track.GetDynamicParticle()->GetLogKineticEnergy()); >> 1006 return (0.0 < xs) ? 1.0/xs : DBL_MAX; >> 1007 } >> 1008 >> 1009 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1010 705 G4double 1011 G4double 706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou 1012 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kinEnergy, 707 G4dou 1013 G4double Z, G4double A, G4double cut) 708 { 1014 { 709 SelectModel(kinEnergy, currentCoupleIndex); 1015 SelectModel(kinEnergy, currentCoupleIndex); 710 return (currentModel) ? 1016 return (currentModel) ? 711 currentModel->ComputeCrossSectionPerAtom(c 1017 currentModel->ComputeCrossSectionPerAtom(currentParticle, kinEnergy, 712 Z 1018 Z, A, cut) : 0.0; 713 } 1019 } 714 1020 715 //....oooOO0OOooo........oooOO0OOooo........oo 1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 716 1022 717 G4PhysicsVector* << 1023 void G4VEmProcess::FindLambdaMax() 718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 719 { 1024 { 720 DefineMaterial(couple); << 1025 if(1 < verboseLevel) { 721 G4PhysicsVector* newv = new G4PhysicsLogVect << 1026 G4cout << "### G4VEmProcess::FindLambdaMax: " 722 << 1027 << particle->GetParticleName() 723 return newv; << 1028 << " and process " << GetProcessName() << " " << G4endl; 724 } << 1029 } >> 1030 size_t n = theLambdaTable->length(); >> 1031 >> 1032 G4PhysicsVector* pv; >> 1033 G4double e, ss, emax, smax; 725 1034 726 //....oooOO0OOooo........oooOO0OOooo........oo << 1035 size_t i; 727 1036 728 const G4Element* G4VEmProcess::GetCurrentEleme << 1037 // first loop on existing vectors 729 { << 1038 for (i=0; i<n; ++i) { 730 return (nullptr != currentModel) ? << 1039 pv = (*theLambdaTable)[i]; 731 currentModel->GetCurrentElement(currentMat << 1040 if(nullptr != pv) { >> 1041 size_t nb = pv->GetVectorLength(); >> 1042 emax = DBL_MAX; >> 1043 smax = 0.0; >> 1044 if(nb > 0) { >> 1045 for (size_t j=0; j<nb; ++j) { >> 1046 e = pv->Energy(j); >> 1047 ss = (*pv)(j); >> 1048 if(ss > smax) { >> 1049 smax = ss; >> 1050 emax = e; >> 1051 } else { >> 1052 break; >> 1053 } >> 1054 } >> 1055 } >> 1056 (*theEnergyOfCrossSectionMax)[i] = emax; >> 1057 if(1 < verboseLevel) { >> 1058 G4cout << "For " << particle->GetParticleName() >> 1059 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV >> 1060 << " lambda= " << smax << G4endl; >> 1061 } >> 1062 } >> 1063 } >> 1064 // second loop using base materials >> 1065 for (i=0; i<n; ++i) { >> 1066 pv = (*theLambdaTable)[i]; >> 1067 if(nullptr == pv) { >> 1068 G4int j = (*theDensityIdx)[i]; >> 1069 (*theEnergyOfCrossSectionMax)[i] = (*theEnergyOfCrossSectionMax)[j]; >> 1070 } >> 1071 } 732 } 1072 } 733 1073 734 //....oooOO0OOooo........oooOO0OOooo........oo 1074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 735 1075 736 const G4Element* G4VEmProcess::GetTargetElemen << 1076 G4PhysicsVector* >> 1077 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple) 737 { 1078 { 738 return (nullptr != currentModel) ? << 1079 DefineMaterial(couple); 739 currentModel->GetCurrentElement(currentMat << 1080 G4PhysicsVector* newv = nullptr; >> 1081 if(nullptr == theLambdaTable) { >> 1082 newv = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, >> 1083 nLambdaBins, splineFlag); >> 1084 } else { >> 1085 newv = new G4PhysicsVector(*((*theLambdaTable)[basedCoupleIndex])); >> 1086 } >> 1087 return newv; 740 } 1088 } 741 1089 742 //....oooOO0OOooo........oooOO0OOooo........oo 1090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 743 1091 744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 1092 const G4Element* G4VEmProcess::GetCurrentElement() const 745 { 1093 { 746 return (nullptr != currentModel) ? << 1094 return (nullptr != currentModel) ? currentModel->GetCurrentElement() : nullptr; 747 currentModel->GetCurrentIsotope(GetCurrent << 748 } 1095 } 749 1096 750 //....oooOO0OOooo........oooOO0OOooo........oo 1097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 751 1098 752 void G4VEmProcess::SetCrossSectionBiasingFacto 1099 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag) 753 { 1100 { 754 if(f > 0.0) { 1101 if(f > 0.0) { 755 biasFactor = f; 1102 biasFactor = f; 756 weightFlag = flag; 1103 weightFlag = flag; 757 if(1 < verboseLevel) { 1104 if(1 < verboseLevel) { 758 G4cout << "### SetCrossSectionBiasingFac 1105 G4cout << "### SetCrossSectionBiasingFactor: for " 759 << particle->GetParticleName() 1106 << particle->GetParticleName() 760 << " and process " << GetProcessN 1107 << " and process " << GetProcessName() 761 << " biasFactor= " << f << " weig 1108 << " biasFactor= " << f << " weightFlag= " << flag 762 << G4endl; 1109 << G4endl; 763 } 1110 } 764 } 1111 } 765 } 1112 } 766 1113 767 //....oooOO0OOooo........oooOO0OOooo........oo 1114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 768 1115 769 void 1116 void 770 G4VEmProcess::ActivateForcedInteraction(G4doub 1117 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r, 771 G4bool 1118 G4bool flag) 772 { 1119 { 773 if(nullptr == biasManager) { biasManager = n 1120 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); } 774 if(1 < verboseLevel) { 1121 if(1 < verboseLevel) { 775 G4cout << "### ActivateForcedInteraction: 1122 G4cout << "### ActivateForcedInteraction: for " 776 << particle->GetParticleName() 1123 << particle->GetParticleName() 777 << " and process " << GetProcessNam 1124 << " and process " << GetProcessName() 778 << " length(mm)= " << length/mm 1125 << " length(mm)= " << length/mm 779 << " in G4Region <" << r 1126 << " in G4Region <" << r 780 << "> weightFlag= " << flag 1127 << "> weightFlag= " << flag 781 << G4endl; 1128 << G4endl; 782 } 1129 } 783 weightFlag = flag; 1130 weightFlag = flag; 784 biasManager->ActivateForcedInteraction(lengt 1131 biasManager->ActivateForcedInteraction(length, r); 785 } 1132 } 786 1133 787 //....oooOO0OOooo........oooOO0OOooo........oo 1134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 788 1135 789 void 1136 void 790 G4VEmProcess::ActivateSecondaryBiasing(const G 1137 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region, 791 G4double factor, 1138 G4double factor, 792 G4double energyLimit) 1139 G4double energyLimit) 793 { 1140 { 794 if (0.0 <= factor) { 1141 if (0.0 <= factor) { 795 1142 796 // Range cut can be applied only for e- 1143 // Range cut can be applied only for e- 797 if(0.0 == factor && secondaryParticle != G 1144 if(0.0 == factor && secondaryParticle != G4Electron::Electron()) 798 { return; } 1145 { return; } 799 1146 800 if(!biasManager) { biasManager = new G4EmB 1147 if(!biasManager) { biasManager = new G4EmBiasingManager(); } 801 biasManager->ActivateSecondaryBiasing(regi 1148 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit); 802 if(1 < verboseLevel) { 1149 if(1 < verboseLevel) { 803 G4cout << "### ActivateSecondaryBiasing: 1150 G4cout << "### ActivateSecondaryBiasing: for " 804 << " process " << GetProcessName() 1151 << " process " << GetProcessName() 805 << " factor= " << factor 1152 << " factor= " << factor 806 << " in G4Region <" << region 1153 << " in G4Region <" << region 807 << "> energyLimit(MeV)= " << energyLimi 1154 << "> energyLimit(MeV)= " << energyLimit/MeV 808 << G4endl; 1155 << G4endl; 809 } 1156 } 810 } 1157 } 811 } 1158 } 812 1159 813 //....oooOO0OOooo........oooOO0OOooo........oo 1160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 814 1161 815 void G4VEmProcess::SetLambdaBinning(G4int n) 1162 void G4VEmProcess::SetLambdaBinning(G4int n) 816 { 1163 { 817 if(5 < n && n < 10000000) { 1164 if(5 < n && n < 10000000) { 818 nLambdaBins = n; 1165 nLambdaBins = n; 819 actBinning = true; 1166 actBinning = true; 820 } else { 1167 } else { 821 G4double e = (G4double)n; 1168 G4double e = (G4double)n; 822 PrintWarning("SetLambdaBinning", e); 1169 PrintWarning("SetLambdaBinning", e); 823 } 1170 } 824 } 1171 } 825 1172 826 //....oooOO0OOooo........oooOO0OOooo........oo 1173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 827 1174 828 void G4VEmProcess::SetMinKinEnergy(G4double e) 1175 void G4VEmProcess::SetMinKinEnergy(G4double e) 829 { 1176 { 830 if(1.e-3*eV < e && e < maxKinEnergy) { 1177 if(1.e-3*eV < e && e < maxKinEnergy) { 831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma 1178 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e) 832 /G4Log(maxKinEnergy/ 1179 /G4Log(maxKinEnergy/minKinEnergy)); 833 minKinEnergy = e; 1180 minKinEnergy = e; 834 actMinKinEnergy = true; 1181 actMinKinEnergy = true; 835 } else { PrintWarning("SetMinKinEnergy", e); 1182 } else { PrintWarning("SetMinKinEnergy", e); } 836 } 1183 } 837 1184 838 //....oooOO0OOooo........oooOO0OOooo........oo 1185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 839 1186 840 void G4VEmProcess::SetMaxKinEnergy(G4double e) 1187 void G4VEmProcess::SetMaxKinEnergy(G4double e) 841 { 1188 { 842 if(minKinEnergy < e && e < 1.e+6*TeV) { 1189 if(minKinEnergy < e && e < 1.e+6*TeV) { 843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ 1190 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy) 844 /G4Log(maxKinEnergy/ 1191 /G4Log(maxKinEnergy/minKinEnergy)); 845 maxKinEnergy = e; 1192 maxKinEnergy = e; 846 actMaxKinEnergy = true; 1193 actMaxKinEnergy = true; 847 } else { PrintWarning("SetMaxKinEnergy", e); 1194 } else { PrintWarning("SetMaxKinEnergy", e); } 848 } 1195 } 849 1196 850 //....oooOO0OOooo........oooOO0OOooo........oo 1197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 851 1198 852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl 1199 void G4VEmProcess::SetMinKinEnergyPrim(G4double e) 853 { 1200 { 854 if(theParameters->MinKinEnergy() <= e && 1201 if(theParameters->MinKinEnergy() <= e && 855 e <= theParameters->MaxKinEnergy()) { min 1202 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; } 856 else { PrintWarning("SetMinKinEnergyPrim", e 1203 else { PrintWarning("SetMinKinEnergyPrim", e); } 857 } 1204 } 858 1205 859 //....oooOO0OOooo........oooOO0OOooo........oo 1206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 860 1207 861 G4VEmProcess* G4VEmProcess::GetEmProcess(const 1208 G4VEmProcess* G4VEmProcess::GetEmProcess(const G4String& nam) 862 { 1209 { 863 return (nam == GetProcessName()) ? this : nu 1210 return (nam == GetProcessName()) ? this : nullptr; 864 } 1211 } 865 1212 866 //....oooOO0OOooo........oooOO0OOooo........oo 1213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 867 1214 >> 1215 G4double >> 1216 G4VEmProcess::GetLambda(G4double kinEnergy, const G4MaterialCutsCouple* couple) >> 1217 { >> 1218 CurrentSetup(couple, kinEnergy); >> 1219 return GetCurrentLambda(kinEnergy, G4Log(kinEnergy)); >> 1220 } >> 1221 >> 1222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1223 868 G4double G4VEmProcess::PolarAngleLimit() const 1224 G4double G4VEmProcess::PolarAngleLimit() const 869 { 1225 { 870 return theParameters->MscThetaLimit(); 1226 return theParameters->MscThetaLimit(); 871 } 1227 } 872 1228 873 //....oooOO0OOooo........oooOO0OOooo........oo 1229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 874 1230 875 void G4VEmProcess::PrintWarning(G4String tit, 1231 void G4VEmProcess::PrintWarning(G4String tit, G4double val) 876 { 1232 { 877 G4String ss = "G4VEmProcess::" + tit; 1233 G4String ss = "G4VEmProcess::" + tit; 878 G4ExceptionDescription ed; 1234 G4ExceptionDescription ed; 879 ed << "Parameter is out of range: " << val 1235 ed << "Parameter is out of range: " << val 880 << " it will have no effect!\n" << " Pro 1236 << " it will have no effect!\n" << " Process " 881 << GetProcessName() << " nbins= " << the 1237 << GetProcessName() << " nbins= " << theParameters->NumberOfBins() 882 << " Emin(keV)= " << theParameters->MinKi 1238 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV 883 << " Emax(GeV)= " << theParameters->MaxKi 1239 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV; 884 G4Exception(ss, "em0044", JustWarning, ed); 1240 G4Exception(ss, "em0044", JustWarning, ed); 885 } 1241 } 886 1242 887 //....oooOO0OOooo........oooOO0OOooo........oo 1243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 888 1244 889 void G4VEmProcess::ProcessDescription(std::ost 1245 void G4VEmProcess::ProcessDescription(std::ostream& out) const 890 { 1246 { 891 if(nullptr != particle) { << 1247 if(particle) { 892 StreamInfo(out, *particle, true); 1248 StreamInfo(out, *particle, true); 893 } 1249 } 894 } 1250 } 895 1251 896 //....oooOO0OOooo........oooOO0OOooo........oo 1252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 897 1253