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,v 1.88 2010-08-17 17:36:59 vnivanch Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ >> 28 // 26 // ------------------------------------------- 29 // ------------------------------------------------------------------- 27 // 30 // 28 // GEANT4 Class file 31 // GEANT4 Class file 29 // 32 // 30 // 33 // 31 // File name: G4VEmProcess 34 // File name: G4VEmProcess 32 // 35 // 33 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of Laszlo Urban code 34 // 37 // 35 // Creation date: 01.10.2003 38 // Creation date: 01.10.2003 36 // 39 // 37 // Modifications: by V.Ivanchenko << 40 // Modifications: >> 41 // 30-06-04 make it to be pure discrete process (V.Ivanchenko) >> 42 // 30-09-08 optimise integral option (V.Ivanchenko) >> 43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko) >> 44 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI) >> 45 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko) >> 46 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko) >> 47 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko) >> 48 // 25-07-05 Add protection: integral mode only for charged particles (VI) >> 49 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko) >> 50 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI) >> 51 // 12-09-06 add SetModel() (mma) >> 52 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko) >> 53 // 27-10-07 Virtual functions moved to source (V.Ivanchenko) >> 54 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) >> 55 // 17-02-10 Added pointer currentParticle (VI) 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" << 49 #include "G4SystemOfUnits.hh" << 50 #include "G4ProcessManager.hh" 66 #include "G4ProcessManager.hh" 51 #include "G4LossTableManager.hh" 67 #include "G4LossTableManager.hh" 52 #include "G4LossTableBuilder.hh" 68 #include "G4LossTableBuilder.hh" 53 #include "G4Step.hh" 69 #include "G4Step.hh" 54 #include "G4ParticleDefinition.hh" 70 #include "G4ParticleDefinition.hh" 55 #include "G4VEmModel.hh" 71 #include "G4VEmModel.hh" 56 #include "G4DataVector.hh" 72 #include "G4DataVector.hh" 57 #include "G4PhysicsTable.hh" 73 #include "G4PhysicsTable.hh" 58 #include "G4EmDataHandler.hh" << 59 #include "G4PhysicsLogVector.hh" 74 #include "G4PhysicsLogVector.hh" 60 #include "G4VParticleChange.hh" 75 #include "G4VParticleChange.hh" 61 #include "G4ProductionCutsTable.hh" 76 #include "G4ProductionCutsTable.hh" 62 #include "G4Region.hh" 77 #include "G4Region.hh" 63 #include "G4Gamma.hh" 78 #include "G4Gamma.hh" 64 #include "G4Electron.hh" 79 #include "G4Electron.hh" 65 #include "G4Positron.hh" 80 #include "G4Positron.hh" 66 #include "G4PhysicsTableHelper.hh" 81 #include "G4PhysicsTableHelper.hh" 67 #include "G4EmBiasingManager.hh" 82 #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" << 74 #include "G4Log.hh" << 75 #include <iostream> << 76 83 77 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 85 79 G4VEmProcess::G4VEmProcess(const G4String& nam 86 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type): 80 G4VDiscreteProcess(name, type) << 87 G4VDiscreteProcess(name, type), >> 88 secondaryParticle(0), >> 89 buildLambdaTable(true), >> 90 theLambdaTable(0), >> 91 theDensityFactor(0), >> 92 theDensityIdx(0), >> 93 integral(false), >> 94 applyCuts(false), >> 95 startFromNull(false), >> 96 currentModel(0), >> 97 particle(0), >> 98 currentParticle(0), >> 99 currentCouple(0) 81 { 100 { 82 theParameters = G4EmParameters::Instance(); << 83 SetVerboseLevel(1); 101 SetVerboseLevel(1); 84 102 85 // Size of tables << 103 // Size of tables assuming spline 86 minKinEnergy = 0.1*CLHEP::keV; << 104 minKinEnergy = 0.1*keV; 87 maxKinEnergy = 100.0*CLHEP::TeV; << 105 maxKinEnergy = 10.0*TeV; >> 106 nLambdaBins = 77; 88 107 89 // default lambda factor 108 // default lambda factor 90 invLambdaFactor = 1.0/lambdaFactor; << 109 lambdaFactor = 0.8; >> 110 >> 111 // default limit on polar angle >> 112 polarAngleLimit = 0.0; >> 113 biasFactor = 1.0; 91 114 92 // particle types 115 // particle types 93 theGamma = G4Gamma::Gamma(); << 116 theGamma = G4Gamma::Gamma(); 94 theElectron = G4Electron::Electron(); << 117 theElectron = G4Electron::Electron(); 95 thePositron = G4Positron::Positron(); << 118 thePositron = G4Positron::Positron(); 96 119 97 pParticleChange = &fParticleChange; 120 pParticleChange = &fParticleChange; 98 fParticleChange.SetSecondaryWeightByProcess( << 99 secParticles.reserve(5); 121 secParticles.reserve(5); 100 122 >> 123 preStepLambda = 0.0; >> 124 mfpKinEnergy = DBL_MAX; >> 125 101 modelManager = new G4EmModelManager(); 126 modelManager = new G4EmModelManager(); 102 lManager = G4LossTableManager::Instance(); << 127 biasManager = 0; 103 lManager->Register(this); << 128 biasFlag = false; 104 isTheMaster = lManager->IsMaster(); << 129 weightFlag = false; 105 G4LossTableBuilder* bld = lManager->GetTable << 130 (G4LossTableManager::Instance())->Register(this); 106 theDensityFactor = bld->GetDensityFactors(); << 107 theDensityIdx = bld->GetCoupleIndexes(); << 108 } 131 } 109 132 110 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 111 134 112 G4VEmProcess::~G4VEmProcess() 135 G4VEmProcess::~G4VEmProcess() 113 { 136 { 114 if(isTheMaster) { << 137 if(1 < verboseLevel) { 115 delete theData; << 138 G4cout << "G4VEmProcess destruct " << GetProcessName() 116 delete theEnergyOfCrossSectionMax; << 139 << G4endl; >> 140 } >> 141 Clear(); >> 142 if(theLambdaTable) { >> 143 theLambdaTable->clearAndDestroy(); >> 144 delete theLambdaTable; 117 } 145 } 118 delete modelManager; 146 delete modelManager; 119 delete biasManager; 147 delete biasManager; 120 lManager->DeRegister(this); << 148 (G4LossTableManager::Instance())->DeRegister(this); 121 } 149 } 122 150 123 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 152 125 void G4VEmProcess::AddEmModel(G4int order, G4V << 153 void G4VEmProcess::Clear() 126 const G4Region* << 127 { 154 { 128 if(nullptr == ptr) { return; } << 155 currentCouple = 0; 129 G4VEmFluctuationModel* fm = nullptr; << 156 preStepLambda = 0.0; 130 modelManager->AddEmModel(order, ptr, fm, reg << 157 mfpKinEnergy = DBL_MAX; 131 ptr->SetParticleChange(pParticleChange); << 132 } 158 } 133 159 134 //....oooOO0OOooo........oooOO0OOooo........oo 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 135 161 136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 162 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*, >> 163 const G4Material*) 137 { 164 { 138 if(nullptr == ptr) { return; } << 165 return 0.0; 139 if(!emModels.empty()) { << 140 for(auto & em : emModels) { if(em == ptr) << 141 } << 142 emModels.push_back(ptr); << 143 } 166 } 144 167 145 //....oooOO0OOooo........oooOO0OOooo........oo 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 146 169 147 void G4VEmProcess::PreparePhysicsTable(const G << 170 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p, >> 171 const G4Region* region) 148 { 172 { 149 if(nullptr == particle) { SetParticle(&part) << 173 G4VEmFluctuationModel* fm = 0; >> 174 modelManager->AddEmModel(order, p, fm, region); >> 175 if(p) { p->SetParticleChange(pParticleChange); } >> 176 } 150 177 151 if(part.GetParticleType() == "nucleus" && << 178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 152 part.GetParticleSubType() == "generic") { << 153 179 154 G4String pname = part.GetParticleName(); << 180 void G4VEmProcess::SetModel(G4VEmModel* p, G4int index) 155 if(pname != "deuteron" && pname != "triton << 181 { 156 pname != "He3" && pname != "alpha" && p << 182 G4int n = emModels.size(); 157 pname != "helium" && pname != "hydrogen << 183 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} } >> 184 emModels[index] = p; >> 185 } 158 186 159 particle = G4GenericIon::GenericIon(); << 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 160 isIon = true; << 161 } << 162 } << 163 if(particle != &part) { return; } << 164 188 165 lManager->PreparePhysicsTable(&part, this); << 189 G4VEmModel* G4VEmProcess::Model(G4int index) >> 190 { >> 191 G4VEmModel* p = 0; >> 192 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; } >> 193 return p; >> 194 } 166 195 167 // for new run << 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 168 currentCouple = nullptr; << 169 preStepLambda = 0.0; << 170 fLambdaEnergy = 0.0; << 171 197 172 InitialiseProcess(particle); << 198 void G4VEmProcess::UpdateEmModel(const G4String& nam, >> 199 G4double emin, G4double emax) >> 200 { >> 201 modelManager->UpdateEmModel(nam, emin, emax); >> 202 } 173 203 174 G4LossTableBuilder* bld = lManager->GetTable << 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 175 const G4ProductionCutsTable* theCoupleTable= << 176 G4ProductionCutsTable::GetProductionCutsTa << 177 theCutsGamma = theCoupleTable->GetEnergyC << 178 theCutsElectron = theCoupleTable->GetEnergyC << 179 theCutsPositron = theCoupleTable->GetEnergyC << 180 << 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 << 194 if(isTheMaster) { << 195 if(nullptr == theData) { theData = new G4E << 196 205 197 if(buildLambdaTable) { << 206 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver) 198 theLambdaTable = theData->MakeTable(0); << 207 { >> 208 return modelManager->GetModel(idx, ver); >> 209 } >> 210 >> 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 212 >> 213 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part) >> 214 { >> 215 if(!particle) { SetParticle(&part); } >> 216 if(1 < verboseLevel) { >> 217 G4cout << "G4VEmProcess::PreparePhysicsTable() for " >> 218 << GetProcessName() >> 219 << " and particle " << part.GetParticleName() >> 220 << " local particle " << particle->GetParticleName() >> 221 << G4endl; >> 222 } >> 223 >> 224 G4LossTableManager* man = G4LossTableManager::Instance(); >> 225 G4LossTableBuilder* bld = man->GetTableBuilder(); >> 226 >> 227 man->PreparePhysicsTable(&part, this); >> 228 >> 229 if(particle == &part) { >> 230 Clear(); >> 231 InitialiseProcess(particle); >> 232 >> 233 const G4ProductionCutsTable* theCoupleTable= >> 234 G4ProductionCutsTable::GetProductionCutsTable(); >> 235 size_t n = theCoupleTable->GetTableSize(); >> 236 >> 237 theEnergyOfCrossSectionMax.resize(n, 0.0); >> 238 theCrossSectionMax.resize(n, DBL_MAX); >> 239 >> 240 // initialisation of models >> 241 G4int nmod = modelManager->NumberOfModels(); >> 242 for(G4int i=0; i<nmod; ++i) { >> 243 G4VEmModel* mod = modelManager->GetModel(i); >> 244 mod->SetPolarAngleLimit(polarAngleLimit); >> 245 if(mod->HighEnergyLimit() > maxKinEnergy) { >> 246 mod->SetHighEnergyLimit(maxKinEnergy); >> 247 } >> 248 } >> 249 >> 250 if(man->AtomDeexcitation()) { modelManager->SetFluoFlag(true); } >> 251 theCuts = modelManager->Initialise(particle,secondaryParticle, >> 252 2.,verboseLevel); >> 253 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut); >> 254 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut); >> 255 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut); >> 256 >> 257 // prepare tables >> 258 if(buildLambdaTable){ >> 259 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); 199 bld->InitialiseBaseMaterials(theLambdaTa 260 bld->InitialiseBaseMaterials(theLambdaTable); 200 } 261 } 201 // high energy table << 262 // forced biasing 202 if(minKinEnergyPrim < maxKinEnergy) { << 263 if(biasManager) { 203 theLambdaTablePrim = theData->MakeTable( << 264 biasManager->Initialise(part,GetProcessName(),verboseLevel); 204 bld->InitialiseBaseMaterials(theLambdaTa << 265 biasFlag = false; 205 } << 266 } 206 } << 267 } 207 // models << 268 theDensityFactor = bld->GetDensityFactors(); 208 baseMat = bld->GetBaseMaterialFlag(); << 269 theDensityIdx = bld->GetCoupleIndexes(); 209 numberOfModels = modelManager->NumberOfModel << 210 currentModel = modelManager->GetModel(0); << 211 if(nullptr != lManager->AtomDeexcitation()) << 212 modelManager->SetFluoFlag(true); << 213 } << 214 // forced biasing << 215 if(nullptr != biasManager) { << 216 biasManager->Initialise(part, GetProcessNa << 217 biasFlag = false; << 218 } << 219 << 220 theCuts = << 221 G4EmTableUtil::PrepareEmProcess(this, part << 222 modelManag << 223 secID, tri << 224 verboseLev << 225 } 270 } 226 271 227 //....oooOO0OOooo........oooOO0OOooo........oo 272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 228 273 229 void G4VEmProcess::BuildPhysicsTable(const G4P 274 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 230 { 275 { 231 if(nullptr == masterProc) { << 276 G4String partname = part.GetParticleName(); 232 if(isTheMaster) { masterProc = this; } << 277 if(1 < verboseLevel) { 233 else { masterProc = static_cast<const G4VE << 278 G4cout << "G4VEmProcess::BuildPhysicsTable() for " 234 } << 279 << GetProcessName() 235 G4int nModels = modelManager->NumberOfModels << 280 << " and particle " << partname 236 G4bool isLocked = theParameters->IsPrintLock << 281 << " buildLambdaTable= " << buildLambdaTable 237 G4bool toBuild = (buildLambdaTable || minKin << 282 << G4endl; 238 << 283 } 239 G4EmTableUtil::BuildEmProcess(this, masterPr << 284 240 nModels, verbo << 285 (G4LossTableManager::Instance())->BuildPhysicsTable(particle); 241 isLocked, toBu << 286 >> 287 if(buildLambdaTable) { >> 288 BuildLambdaTable(); >> 289 FindLambdaMax(); >> 290 } >> 291 >> 292 // reduce printout for nuclear stopping >> 293 G4bool gproc = true; >> 294 G4int st = GetProcessSubType(); >> 295 if((st == fCoulombScattering || st == fNuclearStopping) && >> 296 part.GetParticleType() == "nucleus" && >> 297 partname != "GenericIon" && partname != "alpha") { gproc = false; } >> 298 >> 299 if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); } >> 300 >> 301 if(1 < verboseLevel) { >> 302 G4cout << "G4VEmProcess::BuildPhysicsTable() done for " >> 303 << GetProcessName() >> 304 << " and particle " << partname >> 305 << G4endl; >> 306 } 242 } 307 } 243 308 244 //....oooOO0OOooo........oooOO0OOooo........oo 309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 245 310 246 void G4VEmProcess::BuildLambdaTable() 311 void G4VEmProcess::BuildLambdaTable() 247 { 312 { 248 G4double scale = theParameters->MaxKinEnergy << 313 if(1 < verboseLevel) { 249 G4int nbin = << 314 G4cout << "G4EmProcess::BuildLambdaTable() for process " 250 theParameters->NumberOfBinsPerDecade()*G4l << 315 << GetProcessName() << " and particle " 251 if(actBinning) { nbin = std::max(nbin, nLamb << 316 << particle->GetParticleName() 252 scale = nbin/G4Log(scale); << 317 << G4endl; 253 << 318 } 254 G4LossTableBuilder* bld = lManager->GetTable << 255 G4EmTableUtil::BuildLambdaTable(this, partic << 256 bld, theLamb << 257 minKinEnergy << 258 maxKinEnergy << 259 startFromNul << 260 } << 261 319 262 //....oooOO0OOooo........oooOO0OOooo........oo << 320 // Access to materials >> 321 const G4ProductionCutsTable* theCoupleTable= >> 322 G4ProductionCutsTable::GetProductionCutsTable(); >> 323 size_t numOfCouples = theCoupleTable->GetTableSize(); 263 324 264 void G4VEmProcess::StreamInfo(std::ostream& ou << 325 G4LossTableBuilder* bld = (G4LossTableManager::Instance())->GetTableBuilder(); 265 const G4ParticleDefinition& << 326 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag(); 266 { << 327 267 G4String indent = (rst ? " " : ""); << 328 G4PhysicsLogVector* aVector = 0; 268 out << std::setprecision(6); << 329 G4PhysicsLogVector* bVector = 0; 269 out << G4endl << indent << GetProcessName() << 330 270 if (!rst) { << 331 G4double scale = 0.0; 271 out << " for " << part.GetParticleName(); << 332 if(startFromNull) { scale = std::log(maxKinEnergy/minKinEnergy); } 272 } << 333 273 if(fXSType != fEmNoIntegral) { out << " XSt << 334 for(size_t i=0; i<numOfCouples; ++i) { 274 if(applyCuts) { out << " applyCuts:1 "; } << 335 275 G4int subtype = GetProcessSubType(); << 336 if (/*theLambdaTable->GetFlag(i)*/bld->GetFlag(i)) { 276 out << " SubType=" << subtype; << 337 277 if (subtype == fAnnihilation) { << 338 // create physics vector and fill it 278 G4int mod = theParameters->PositronAtRestM << 339 const G4MaterialCutsCouple* couple = 279 const G4String namp[2] = {"Simple", "Allis << 340 theCoupleTable->GetMaterialCutsCouple(i); 280 out << " AtRestModel:" << namp[mod]; << 341 281 } << 342 // if start from zero then change the scale 282 if(biasFactor != 1.0) { out << " BiasingFac << 343 if(startFromNull) { 283 out << " BuildTable=" << buildLambdaTable << << 344 G4double emin = MinPrimaryEnergy(particle,couple->GetMaterial()); 284 if(buildLambdaTable) { << 345 if(0.0 >= emin) { emin = eV; } 285 if(particle == &part) { << 346 else if(maxKinEnergy <= emin) { emin = 0.5*maxKinEnergy; } 286 for(auto & v : *theLambdaTable) { << 347 G4int bin = 287 if(nullptr != v) { << 348 G4int(nLambdaBins*std::log(maxKinEnergy/emin)/scale + 0.5); 288 out << " Lambda table from "; << 349 if(bin < 3) { bin = 3; } 289 G4double emin = v->Energy(0); << 350 aVector = new G4PhysicsLogVector(emin, maxKinEnergy, bin); 290 G4double emax = v->GetMaxEnergy(); << 351 291 G4int nbin = G4int(v->GetVectorLengt << 352 // start not from zero 292 if(emin > minKinEnergy) { out << "th << 353 } else if(!bVector) { 293 else { out << G4BestUnit(emin,"Energ << 354 aVector = 294 out << " to " << 355 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins); 295 << G4BestUnit(emax,"Energy") << 356 bVector = aVector; 296 << ", " << G4lrint(nbin/std::log << 357 } else { 297 << " bins/decade, spline: " << 358 aVector = new G4PhysicsLogVector(*bVector); 298 << splineFlag << G4endl; << 299 break; << 300 } << 301 } << 302 } else { << 303 out << " Used Lambda table of " << 304 << particle->GetParticleName() << G4endl << 305 } << 306 } << 307 if(minKinEnergyPrim < maxKinEnergy) { << 308 if(particle == &part) { << 309 for(auto & v : *theLambdaTablePrim) { << 310 if(nullptr != v) { << 311 out << " LambdaPrime table from << 312 << G4BestUnit(v->Energy(0),"Ener << 313 << " to " << 314 << G4BestUnit(v->GetMaxEnergy(), << 315 << " in " << v->GetVectorLength( << 316 << " bins " << G4endl; << 317 break; << 318 } << 319 } 359 } 320 } else { << 360 aVector->SetSpline(splineFlag); 321 out << " Used LambdaPrime table of << 361 modelManager->FillLambdaVector(aVector, couple, startFromNull); 322 << particle->GetParticleName() << 362 if(splineFlag) { aVector->FillSecondDerivatives(); } >> 363 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector); 323 } 364 } 324 } 365 } 325 StreamProcessInfo(out); << 326 modelManager->DumpModelList(out, verboseLeve << 327 366 328 if(verboseLevel > 2 && buildLambdaTable) { << 367 if(1 < verboseLevel) { 329 out << " LambdaTable address= " << th << 368 G4cout << "Lambda table is built for " 330 if(theLambdaTable && particle == &part) { << 369 << particle->GetParticleName() 331 out << (*theLambdaTable) << G4endl; << 370 << G4endl; 332 } << 333 } 371 } 334 } 372 } 335 373 336 //....oooOO0OOooo........oooOO0OOooo........oo 374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 337 375 338 void G4VEmProcess::StartTracking(G4Track* trac << 376 void G4VEmProcess::PrintInfoDefinition() 339 { 377 { 340 // reset parameters for the new track << 378 if(verboseLevel > 0) { 341 currentParticle = track->GetParticleDefiniti << 379 G4cout << G4endl << GetProcessName() << ": for " 342 theNumberOfInteractionLengthLeft = -1.0; << 380 << particle->GetParticleName(); 343 mfpKinEnergy = DBL_MAX; << 381 if(integral) { G4cout << ", integral: 1 "; } 344 preStepLambda = 0.0; << 382 if(applyCuts) { G4cout << ", applyCuts: 1 "; } 345 << 383 G4cout << " SubType= " << GetProcessSubType();; 346 if(isIon) { massRatio = proton_mass_c2/curre << 384 if(biasFactor != 1.0) { G4cout << " BiasingFactor= " << biasFactor; } 347 << 385 G4cout << G4endl; 348 // forced biasing only for primary particles << 386 if(buildLambdaTable) { 349 if(biasManager) { << 387 G4cout << " Lambda tables from " 350 if(0 == track->GetParentID()) { << 388 << G4BestUnit(minKinEnergy,"Energy") 351 // primary particle << 389 << " to " 352 biasFlag = true; << 390 << G4BestUnit(maxKinEnergy,"Energy") 353 biasManager->ResetForcedInteraction(); << 391 << " in " << nLambdaBins << " bins, spline: " >> 392 << (G4LossTableManager::Instance())->SplineFlag() >> 393 << G4endl; 354 } 394 } >> 395 PrintInfo(); >> 396 modelManager->DumpModelList(verboseLevel); >> 397 } >> 398 >> 399 if(verboseLevel > 2 && buildLambdaTable) { >> 400 G4cout << " LambdaTable address= " << theLambdaTable << G4endl; >> 401 if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; } 355 } 402 } 356 } 403 } 357 404 358 //....oooOO0OOooo........oooOO0OOooo........oo 405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 359 406 360 G4double G4VEmProcess::PostStepGetPhysicalInte 407 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength( 361 const G4Track& tr 408 const G4Track& track, 362 G4double previo 409 G4double previousStepSize, 363 G4ForceCondition* 410 G4ForceCondition* condition) 364 { 411 { >> 412 // condition is set to "Not Forced" 365 *condition = NotForced; 413 *condition = NotForced; 366 G4double x = DBL_MAX; 414 G4double x = DBL_MAX; 367 << 415 if(previousStepSize <= 0.0) { theNumberOfInteractionLengthLeft = -1.0; } 368 DefineMaterial(track.GetMaterialCutsCouple() << 416 InitialiseStep(track); 369 preStepKinEnergy = track.GetKineticEnergy(); << 417 if(!currentModel->IsActive(preStepKinEnergy)) { return x; } 370 const G4double scaledEnergy = preStepKinEner << 371 SelectModel(scaledEnergy, currentCoupleIndex << 372 /* << 373 G4cout << "PostStepGetPhysicalInteractionLen << 374 << " couple: " << currentCouple << G << 375 */ << 376 if(!currentModel->IsActive(scaledEnergy)) { << 377 theNumberOfInteractionLengthLeft = -1.0; << 378 currentInteractionLength = DBL_MAX; << 379 mfpKinEnergy = DBL_MAX; << 380 preStepLambda = 0.0; << 381 return x; << 382 } << 383 418 384 // forced biasing only for primary particles 419 // forced biasing only for primary particles 385 if(biasManager) { 420 if(biasManager) { 386 if(0 == track.GetParentID()) { 421 if(0 == track.GetParentID()) { 387 if(biasFlag && << 422 if(0 == track.GetCurrentStepNumber()) { 388 biasManager->ForcedInteractionRegion( << 423 biasFlag = true; 389 return biasManager->GetStepLimit((G4in << 424 biasManager->ResetForcedInteraction(); >> 425 } >> 426 if(biasFlag && biasManager->ForcedInteractionRegion(currentCoupleIndex)) { >> 427 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize); 390 } 428 } 391 } 429 } 392 } 430 } 393 431 394 // compute mean free path 432 // compute mean free path >> 433 if(preStepKinEnergy < mfpKinEnergy) { >> 434 if (integral) { ComputeIntegralLambda(preStepKinEnergy); } >> 435 else { preStepLambda = GetCurrentLambda(preStepKinEnergy); } >> 436 if(preStepLambda <= 0.0) { mfpKinEnergy = 0.0; } >> 437 } 395 438 396 ComputeIntegralLambda(preStepKinEnergy, trac << 439 // non-zero cross sect}ion 397 << 440 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 << 405 // non-zero cross section << 406 if (theNumberOfInteractionLengthLeft < 0.0 441 if (theNumberOfInteractionLengthLeft < 0.0) { 407 << 408 // beggining of tracking (or just after 442 // beggining of tracking (or just after DoIt of this process) 409 theNumberOfInteractionLengthLeft = -G4Lo << 443 ResetNumberOfInteractionLengthLeft(); 410 theInitialNumberOfInteractionLength = th << 444 } else if(currentInteractionLength < DBL_MAX) { 411 << 445 // subtract NumberOfInteractionLengthLeft 412 } else { << 446 SubtractNumberOfInteractionLengthLeft(previousStepSize); 413 << 447 if(theNumberOfInteractionLengthLeft < 0.) 414 theNumberOfInteractionLengthLeft -= << 448 theNumberOfInteractionLengthLeft = perMillion; 415 previousStepSize/currentInteractionLen << 416 theNumberOfInteractionLengthLeft = << 417 std::max(theNumberOfInteractionLengthL << 418 } 449 } 419 450 420 // new mean free path and step limit for t << 451 // get mean free path and step limit 421 currentInteractionLength = 1.0/preStepLamb 452 currentInteractionLength = 1.0/preStepLambda; 422 x = theNumberOfInteractionLengthLeft * cur 453 x = theNumberOfInteractionLengthLeft * currentInteractionLength; 423 } << 454 #ifdef G4VERBOSE 424 return x; << 455 if (verboseLevel>2){ 425 } << 456 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength "; 426 << 457 G4cout << "[ " << GetProcessName() << "]" << G4endl; 427 //....oooOO0OOooo........oooOO0OOooo........oo << 458 G4cout << " for " << currentParticle->GetParticleName() >> 459 << " in Material " << currentMaterial->GetName() >> 460 << " Ekin(MeV)= " << preStepKinEnergy/MeV >> 461 <<G4endl; >> 462 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" >> 463 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl; >> 464 } >> 465 #endif 428 466 429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 467 // zero cross section case 430 { << 468 } else { 431 if (fXSType == fEmNoIntegral) { << 469 if(theNumberOfInteractionLengthLeft > DBL_MIN && 432 preStepLambda = GetCurrentLambda(e, LogEki << 470 currentInteractionLength < DBL_MAX) { 433 471 434 } else if (fXSType == fEmIncreasing) { << 472 // subtract NumberOfInteractionLengthLeft 435 if(e*invLambdaFactor < mfpKinEnergy) { << 473 SubtractNumberOfInteractionLengthLeft(previousStepSize); 436 preStepLambda = GetCurrentLambda(e, LogE << 474 if(theNumberOfInteractionLengthLeft < 0.) 437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 475 theNumberOfInteractionLengthLeft = perMillion; 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 } 476 } 459 } else { << 477 currentInteractionLength = DBL_MAX; 460 preStepLambda = GetCurrentLambda(e, LogEki << 461 } 478 } >> 479 return x; 462 } 480 } 463 481 464 //....oooOO0OOooo........oooOO0OOooo........oo 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 465 483 466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 484 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track, 467 << 485 const G4Step&) 468 { 486 { 469 // clear number of interaction lengths in an << 470 theNumberOfInteractionLengthLeft = -1.0; << 471 mfpKinEnergy = DBL_MAX; << 472 << 473 fParticleChange.InitializeForPostStep(track) 487 fParticleChange.InitializeForPostStep(track); 474 488 475 // Do not make anything if particle is stopp 489 // Do not make anything if particle is stopped, the annihilation then 476 // should be performed by the AtRestDoIt! 490 // should be performed by the AtRestDoIt! 477 if (track.GetTrackStatus() == fStopButAlive) 491 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; } 478 492 479 const G4double finalT = track.GetKineticEner << 493 G4double finalT = track.GetKineticEnergy(); 480 494 481 // forced process - should happen only once 495 // forced process - should happen only once per track 482 if(biasFlag) { 496 if(biasFlag) { 483 if(biasManager->ForcedInteractionRegion((G << 497 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) { 484 biasFlag = false; 498 biasFlag = false; 485 } 499 } 486 } 500 } 487 501 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 502 // Integral approach 494 if (fXSType != fEmNoIntegral) { << 503 if (integral) { 495 const G4double logFinalT = << 504 G4double lx = GetLambda(finalT, currentCouple); 496 track.GetDynamicParticle()->GetLogKineti << 505 if(preStepLambda<lx && 1 < verboseLevel) { 497 const G4double lx = std::max(GetCurrentLam << 506 G4cout << "WARING: for " << currentParticle->GetParticleName() 498 #ifdef G4VERBOSE << 507 << " and " << GetProcessName() 499 if(preStepLambda < lx && 1 < verboseLevel) << 508 << " E(MeV)= " << finalT/MeV 500 G4cout << "WARNING: for " << currentPart << 509 << " preLambda= " << preStepLambda << " < " 501 << " and " << GetProcessName() << << 510 << lx << " (postLambda) " 502 << " preLambda= " << preStepLambd << 511 << G4endl; 503 << " < " << lx << " (postLambda) << 504 } 512 } 505 #endif << 513 506 // if false interaction then use new cross << 514 if(preStepLambda*G4UniformRand() > lx) { 507 // if both values are zero - no interactio << 515 ClearNumberOfInteractionLengthLeft(); 508 if(preStepLambda*G4UniformRand() >= lx) { << 509 return &fParticleChange; 516 return &fParticleChange; 510 } 517 } 511 } 518 } 512 519 >> 520 SelectModel(finalT, currentCoupleIndex); >> 521 if(!currentModel->IsActive(finalT)) { return &fParticleChange; } >> 522 513 // define new weight for primary and seconda 523 // define new weight for primary and secondaries 514 G4double weight = fParticleChange.GetParentW << 524 if(weightFlag) { 515 if(weightFlag) { << 525 G4double weight = fParticleChange.GetParentWeight()/biasFactor; 516 weight /= biasFactor; << 526 fParticleChange.ProposeParentWeight(weight); 517 fParticleChange.ProposeWeight(weight); << 518 } 527 } 519 << 528 /* 520 #ifdef G4VERBOSE << 529 if(0 < verboseLevel) { 521 if(1 < verboseLevel) { << 522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 530 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= " 523 << finalT/MeV 531 << finalT/MeV 524 << " MeV; model= (" << currentModel 532 << " MeV; model= (" << currentModel->LowEnergyLimit() 525 << ", " << currentModel->HighEnerg 533 << ", " << currentModel->HighEnergyLimit() << ")" 526 << G4endl; 534 << G4endl; 527 } 535 } 528 #endif << 536 */ 529 537 530 // sample secondaries 538 // sample secondaries 531 secParticles.clear(); 539 secParticles.clear(); 532 currentModel->SampleSecondaries(&secParticle 540 currentModel->SampleSecondaries(&secParticles, 533 currentCoupl << 541 currentCouple, 534 track.GetDyn << 542 track.GetDynamicParticle(), 535 (*theCuts)[c << 543 (*theCuts)[currentCoupleIndex]); 536 << 537 G4int num0 = (G4int)secParticles.size(); << 538 << 539 // splitting or Russian roulette << 540 if(biasManager) { << 541 if(biasManager->SecondaryBiasingRegion((G4 << 542 G4double eloss = 0.0; << 543 weight *= biasManager->ApplySecondaryBia << 544 secParticles, track, currentModel, &fP << 545 (G4int)currentCoupleIndex, (*theCuts)[ << 546 step.GetPostStepPoint()->GetSafety()); << 547 if(eloss > 0.0) { << 548 eloss += fParticleChange.GetLocalEnerg << 549 fParticleChange.ProposeLocalEnergyDepo << 550 } << 551 } << 552 } << 553 544 554 // save secondaries 545 // save secondaries 555 G4int num = (G4int)secParticles.size(); << 546 G4int num = secParticles.size(); 556 if(num > 0) { 547 if(num > 0) { 557 548 558 fParticleChange.SetNumberOfSecondaries(num 549 fParticleChange.SetNumberOfSecondaries(num); 559 G4double edep = fParticleChange.GetLocalEn 550 G4double edep = fParticleChange.GetLocalEnergyDeposit(); 560 G4double time = track.GetGlobalTime(); << 561 << 562 G4int n1(0), n2(0); << 563 if(num0 > mainSecondaries) { << 564 currentModel->FillNumberOfSecondaries(n1 << 565 } << 566 551 567 for (G4int i=0; i<num; ++i) { 552 for (G4int i=0; i<num; ++i) { 568 G4DynamicParticle* dp = secParticles[i]; 553 G4DynamicParticle* dp = secParticles[i]; 569 if (nullptr != dp) { << 554 const G4ParticleDefinition* p = dp->GetParticleDefinition(); 570 const G4ParticleDefinition* p = dp->Ge << 555 G4double e = dp->GetKineticEnergy(); 571 G4double e = dp->GetKineticEnergy(); << 556 G4bool good = true; 572 G4bool good = true; << 557 if(applyCuts) { 573 if(applyCuts) { << 558 if (p == theGamma) { 574 if (p == theGamma) { << 559 if (e < (*theCutsGamma)[currentCoupleIndex]) good = false; 575 if (e < (*theCutsGamma)[currentCou << 560 576 << 561 } else if (p == theElectron) { 577 } else if (p == theElectron) { << 562 if (e < (*theCutsElectron)[currentCoupleIndex]) good = false; 578 if (e < (*theCutsElectron)[current << 563 579 << 564 } else if (p == thePositron) { 580 } else if (p == thePositron) { << 565 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] && 581 if (electron_mass_c2 < (*theCutsGa << 566 e < (*theCutsPositron)[currentCoupleIndex]) { 582 e < (*theCutsPositron)[current << 567 good = false; 583 good = false; << 568 e += 2.0*electron_mass_c2; 584 e += 2.0*electron_mass_c2; << 569 } 585 } << 570 } 586 } << 571 if(!good) { 587 // added secondary if it is good << 572 delete dp; 588 } << 573 edep += e; 589 if (good) { << 574 } 590 G4Track* t = new G4Track(dp, time, t << 575 } 591 t->SetTouchableHandle(track.GetTouch << 576 if (good) { fParticleChange.AddSecondary(dp); } 592 if (biasManager) { << 577 } 593 t->SetWeight(weight * biasManager- << 594 } else { << 595 t->SetWeight(weight); << 596 } << 597 pParticleChange->AddSecondary(t); << 598 << 599 // define type of secondary << 600 if(i < mainSecondaries) { << 601 t->SetCreatorModelID(secID); << 602 if(GetProcessSubType() == fCompton << 603 t->SetCreatorModelID(_ComptonGam << 604 } << 605 } else if(i < mainSecondaries + n1) << 606 t->SetCreatorModelID(tripletID); << 607 } else if(i < mainSecondaries + n1 + << 608 t->SetCreatorModelID(_IonRecoil); << 609 } else { << 610 if(i < num0) { << 611 if(p == theGamma) { << 612 t->SetCreatorModelID(fluoID); << 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 { << 627 delete dp; << 628 edep += e; << 629 } << 630 } << 631 } << 632 fParticleChange.ProposeLocalEnergyDeposit( 578 fParticleChange.ProposeLocalEnergyDeposit(edep); 633 } 579 } 634 580 635 if(0.0 == fParticleChange.GetProposedKinetic 581 if(0.0 == fParticleChange.GetProposedKineticEnergy() && 636 fAlive == fParticleChange.GetTrackStatus( 582 fAlive == fParticleChange.GetTrackStatus()) { 637 if(particle->GetProcessManager()->GetAtRes 583 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0) 638 { fParticleChange.ProposeTrackStatus( 584 { fParticleChange.ProposeTrackStatus(fStopButAlive); } 639 else { fParticleChange.ProposeTrackStatus( 585 else { fParticleChange.ProposeTrackStatus(fStopAndKill); } 640 } 586 } 641 587 >> 588 ClearNumberOfInteractionLengthLeft(); 642 return &fParticleChange; 589 return &fParticleChange; 643 } 590 } 644 591 645 //....oooOO0OOooo........oooOO0OOooo........oo 592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 646 593 647 G4bool G4VEmProcess::StorePhysicsTable(const G 594 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part, 648 const G << 595 const G4String& directory, 649 G4bool << 596 G4bool ascii) 650 { 597 { 651 if(!isTheMaster || part != particle) { retur << 598 G4bool yes = true; 652 if(G4EmTableUtil::StoreTable(this, part, the << 599 653 directory, "Lambda", << 600 if ( theLambdaTable && part == particle) { 654 verboseLevel, a << 601 const G4String name = 655 G4EmTableUtil::StoreTable(this, part, the << 602 GetPhysicsTableFileName(part,directory,"Lambda",ascii); 656 directory, "LambdaPrim", << 603 yes = theLambdaTable->StorePhysicsTable(name,ascii); 657 verboseLevel, a << 604 658 return true; << 605 if ( yes ) { >> 606 G4cout << "Physics tables are stored for " << particle->GetParticleName() >> 607 << " and process " << GetProcessName() >> 608 << " in the directory <" << directory >> 609 << "> " << G4endl; >> 610 } else { >> 611 G4cout << "Fail to store Physics Tables for " >> 612 << particle->GetParticleName() >> 613 << " and process " << GetProcessName() >> 614 << " in the directory <" << directory >> 615 << "> " << G4endl; >> 616 } 659 } 617 } 660 return false; << 618 return yes; 661 } 619 } 662 620 663 //....oooOO0OOooo........oooOO0OOooo........oo 621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 664 622 665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 623 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 666 cons << 624 const G4String& directory, 667 G4bo << 625 G4bool ascii) 668 { 626 { 669 if(!isTheMaster || part != particle) { retur << 627 if(1 < verboseLevel) { >> 628 G4cout << "G4VEmProcess::RetrievePhysicsTable() for " >> 629 << part->GetParticleName() << " and process " >> 630 << GetProcessName() << G4endl; >> 631 } 670 G4bool yes = true; 632 G4bool yes = true; 671 if(buildLambdaTable) { << 633 672 yes = G4EmTableUtil::RetrieveTable(this, p << 634 if(!buildLambdaTable || particle != part) return yes; 673 "Lambda << 635 674 ascii, << 636 const G4String particleName = part->GetParticleName(); 675 } << 637 G4String filename; 676 if(yes && minKinEnergyPrim < maxKinEnergy) { << 638 677 yes = G4EmTableUtil::RetrieveTable(this, p << 639 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii); 678 "Lambda << 640 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable, 679 ascii, << 641 filename,ascii); >> 642 if ( yes ) { >> 643 if (0 < verboseLevel) { >> 644 G4cout << "Lambda table for " << particleName >> 645 << " is Retrieved from <" >> 646 << filename << ">" >> 647 << G4endl; >> 648 } >> 649 if((G4LossTableManager::Instance())->SplineFlag()) { >> 650 size_t n = theLambdaTable->length(); >> 651 for(size_t i=0; i<n; ++i) { >> 652 if((* theLambdaTable)[i]) { >> 653 (* theLambdaTable)[i]->SetSpline(true); >> 654 } >> 655 } >> 656 } >> 657 } else { >> 658 if (1 < verboseLevel) { >> 659 G4cout << "Lambda table for " << particleName << " in file <" >> 660 << filename << "> is not exist" >> 661 << G4endl; >> 662 } 680 } 663 } >> 664 681 return yes; 665 return yes; 682 } 666 } 683 667 684 //....oooOO0OOooo........oooOO0OOooo........oo 668 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 685 669 686 G4double G4VEmProcess::GetCrossSection(G4doubl << 670 G4double 687 const G << 671 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy, >> 672 const G4MaterialCutsCouple* couple) 688 { 673 { 689 CurrentSetup(couple, kinEnergy); << 674 // Cross section per atom is calculated 690 return GetCurrentLambda(kinEnergy, G4Log(kin << 675 DefineMaterial(couple); >> 676 G4double cross = 0.0; >> 677 if(theLambdaTable) { >> 678 cross = (*theDensityFactor)[currentCoupleIndex]* >> 679 (((*theLambdaTable)[basedCoupleIndex])->Value(kineticEnergy)); >> 680 } else { >> 681 SelectModel(kineticEnergy, currentCoupleIndex); >> 682 cross = currentModel->CrossSectionPerVolume(currentMaterial, >> 683 currentParticle,kineticEnergy); >> 684 } >> 685 >> 686 if(cross < 0.0) { cross = 0.0; } >> 687 return cross; 691 } 688 } 692 689 693 //....oooOO0OOooo........oooOO0OOooo........oo 690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 694 691 695 G4double G4VEmProcess::GetMeanFreePath(const G 692 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track, 696 G4doubl << 693 G4double, 697 G4Force << 694 G4ForceCondition* condition) 698 { 695 { 699 *condition = NotForced; 696 *condition = NotForced; 700 return G4VEmProcess::MeanFreePath(track); 697 return G4VEmProcess::MeanFreePath(track); 701 } 698 } 702 699 703 //....oooOO0OOooo........oooOO0OOooo........oo 700 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 704 701 705 G4double << 702 G4double G4VEmProcess::MeanFreePath(const G4Track& track) 706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou << 707 G4dou << 708 { 703 { 709 SelectModel(kinEnergy, currentCoupleIndex); << 704 DefineMaterial(track.GetMaterialCutsCouple()); 710 return (currentModel) ? << 705 preStepLambda = GetCurrentLambda(track.GetKineticEnergy()); 711 currentModel->ComputeCrossSectionPerAtom(c << 706 G4double x = DBL_MAX; 712 Z << 707 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; } >> 708 return x; 713 } 709 } 714 710 715 //....oooOO0OOooo........oooOO0OOooo........oo 711 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 716 712 717 G4PhysicsVector* << 713 G4double 718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 714 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy, >> 715 G4double Z, G4double A, G4double cut) 719 { 716 { 720 DefineMaterial(couple); << 717 SelectModel(kineticEnergy, currentCoupleIndex); 721 G4PhysicsVector* newv = new G4PhysicsLogVect << 718 G4double x = 0.0; 722 << 719 if(currentModel) { 723 return newv; << 720 x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy, >> 721 Z,A,cut); >> 722 } >> 723 return x; 724 } 724 } 725 725 726 //....oooOO0OOooo........oooOO0OOooo........oo 726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 727 727 728 const G4Element* G4VEmProcess::GetCurrentEleme << 728 void G4VEmProcess::FindLambdaMax() 729 { 729 { 730 return (nullptr != currentModel) ? << 730 if(1 < verboseLevel) { 731 currentModel->GetCurrentElement(currentMat << 731 G4cout << "### G4VEmProcess::FindLambdaMax: " >> 732 << particle->GetParticleName() >> 733 << " and process " << GetProcessName() << G4endl; >> 734 } >> 735 size_t n = theLambdaTable->length(); >> 736 G4PhysicsVector* pv; >> 737 G4double e, s, emax, smax; >> 738 >> 739 size_t i; >> 740 >> 741 // first loop on existing vectors >> 742 for (i=0; i<n; ++i) { >> 743 pv = (*theLambdaTable)[i]; >> 744 if(pv) { >> 745 size_t nb = pv->GetVectorLength(); >> 746 emax = DBL_MAX; >> 747 smax = 0.0; >> 748 if(nb > 0) { >> 749 for (size_t j=0; j<nb; ++j) { >> 750 e = pv->Energy(j); >> 751 s = (*pv)(j); >> 752 if(s > smax) { >> 753 smax = s; >> 754 emax = e; >> 755 } >> 756 } >> 757 } >> 758 theEnergyOfCrossSectionMax[i] = emax; >> 759 theCrossSectionMax[i] = smax; >> 760 if(1 < verboseLevel) { >> 761 G4cout << "For " << particle->GetParticleName() >> 762 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV >> 763 << " lambda= " << smax << G4endl; >> 764 } >> 765 } >> 766 } >> 767 // second loop using base materials >> 768 for (size_t i=0; i<n; ++i) { >> 769 pv = (*theLambdaTable)[i]; >> 770 if(!pv){ >> 771 G4int j = (*theDensityIdx)[i]; >> 772 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j]; >> 773 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j]; >> 774 } >> 775 } 732 } 776 } 733 777 734 //....oooOO0OOooo........oooOO0OOooo........oo 778 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 735 779 736 const G4Element* G4VEmProcess::GetTargetElemen << 780 G4PhysicsVector* >> 781 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*) 737 { 782 { 738 return (nullptr != currentModel) ? << 783 G4PhysicsVector* v = 739 currentModel->GetCurrentElement(currentMat << 784 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins); >> 785 v->SetSpline((G4LossTableManager::Instance())->SplineFlag()); >> 786 return v; 740 } 787 } 741 788 742 //....oooOO0OOooo........oooOO0OOooo........oo 789 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 743 790 744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 791 const G4Element* G4VEmProcess::GetCurrentElement() const 745 { 792 { 746 return (nullptr != currentModel) ? << 793 const G4Element* elm = 0; 747 currentModel->GetCurrentIsotope(GetCurrent << 794 if(currentModel) {elm = currentModel->GetCurrentElement(); } >> 795 return elm; 748 } 796 } 749 797 750 //....oooOO0OOooo........oooOO0OOooo........oo 798 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 751 799 752 void G4VEmProcess::SetCrossSectionBiasingFacto 800 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag) 753 { 801 { 754 if(f > 0.0) { 802 if(f > 0.0) { 755 biasFactor = f; 803 biasFactor = f; 756 weightFlag = flag; 804 weightFlag = flag; 757 if(1 < verboseLevel) { 805 if(1 < verboseLevel) { 758 G4cout << "### SetCrossSectionBiasingFac 806 G4cout << "### SetCrossSectionBiasingFactor: for " 759 << particle->GetParticleName() << 807 << particle->GetParticleName() 760 << " and process " << GetProcessN << 808 << " and process " << GetProcessName() 761 << " biasFactor= " << f << " weig << 809 << " biasFactor= " << f << " weightFlag= " << flag 762 << G4endl; << 810 << G4endl; 763 } 811 } 764 } 812 } 765 } 813 } 766 814 767 //....oooOO0OOooo........oooOO0OOooo........oo 815 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 768 816 769 void 817 void 770 G4VEmProcess::ActivateForcedInteraction(G4doub 818 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r, 771 G4bool << 819 G4bool flag) 772 { 820 { 773 if(nullptr == biasManager) { biasManager = n << 821 if(!biasManager) { biasManager = new G4EmBiasingManager(); } 774 if(1 < verboseLevel) { 822 if(1 < verboseLevel) { 775 G4cout << "### ActivateForcedInteraction: 823 G4cout << "### ActivateForcedInteraction: for " 776 << particle->GetParticleName() << 824 << particle->GetParticleName() 777 << " and process " << GetProcessNam << 825 << " and process " << GetProcessName() 778 << " length(mm)= " << length/mm << 826 << " length(mm)= " << length/mm 779 << " in G4Region <" << r << 827 << " in G4Region <" << r 780 << "> weightFlag= " << flag << 828 << "> weightFlag= " << flag 781 << G4endl; << 829 << G4endl; 782 } 830 } 783 weightFlag = flag; 831 weightFlag = flag; 784 biasManager->ActivateForcedInteraction(lengt 832 biasManager->ActivateForcedInteraction(length, r); 785 } << 786 << 787 //....oooOO0OOooo........oooOO0OOooo........oo << 788 << 789 void << 790 G4VEmProcess::ActivateSecondaryBiasing(const G << 791 G4double factor, << 792 G4double energyLimit) << 793 { << 794 if (0.0 <= factor) { << 795 << 796 // Range cut can be applied only for e- << 797 if(0.0 == factor && secondaryParticle != G << 798 { return; } << 799 << 800 if(!biasManager) { biasManager = new G4EmB << 801 biasManager->ActivateSecondaryBiasing(regi << 802 if(1 < verboseLevel) { << 803 G4cout << "### ActivateSecondaryBiasing: << 804 << " process " << GetProcessName() << 805 << " factor= " << factor << 806 << " in G4Region <" << region << 807 << "> energyLimit(MeV)= " << energyLimi << 808 << G4endl; << 809 } << 810 } << 811 } << 812 << 813 //....oooOO0OOooo........oooOO0OOooo........oo << 814 << 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) << 829 { << 830 if(1.e-3*eV < e && e < maxKinEnergy) { << 831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma << 832 /G4Log(maxKinEnergy/ << 833 minKinEnergy = e; << 834 actMinKinEnergy = true; << 835 } else { PrintWarning("SetMinKinEnergy", e); << 836 } << 837 << 838 //....oooOO0OOooo........oooOO0OOooo........oo << 839 << 840 void G4VEmProcess::SetMaxKinEnergy(G4double e) << 841 { << 842 if(minKinEnergy < e && e < 1.e+6*TeV) { << 843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ << 844 /G4Log(maxKinEnergy/ << 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 } 833 } 895 834 896 //....oooOO0OOooo........oooOO0OOooo........oo 835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 897 836