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