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