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.35 2006/06/29 19:55:19 gunter Exp $ >> 27 // GEANT4 tag $Name: geant4-08-01 $ >> 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 // 38 // 52 // 39 // Class Description: based class for discrete << 53 // Class Description: 40 // 54 // 41 55 42 // ------------------------------------------- 56 // ------------------------------------------------------------------- 43 // 57 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 60 47 #include "G4VEmProcess.hh" 61 #include "G4VEmProcess.hh" 48 #include "G4PhysicalConstants.hh" << 49 #include "G4SystemOfUnits.hh" << 50 #include "G4ProcessManager.hh" << 51 #include "G4LossTableManager.hh" 62 #include "G4LossTableManager.hh" 52 #include "G4LossTableBuilder.hh" << 53 #include "G4Step.hh" 63 #include "G4Step.hh" 54 #include "G4ParticleDefinition.hh" 64 #include "G4ParticleDefinition.hh" 55 #include "G4VEmModel.hh" 65 #include "G4VEmModel.hh" 56 #include "G4DataVector.hh" 66 #include "G4DataVector.hh" 57 #include "G4PhysicsTable.hh" 67 #include "G4PhysicsTable.hh" 58 #include "G4EmDataHandler.hh" << 68 #include "G4PhysicsVector.hh" 59 #include "G4PhysicsLogVector.hh" 69 #include "G4PhysicsLogVector.hh" 60 #include "G4VParticleChange.hh" 70 #include "G4VParticleChange.hh" 61 #include "G4ProductionCutsTable.hh" 71 #include "G4ProductionCutsTable.hh" 62 #include "G4Region.hh" 72 #include "G4Region.hh" >> 73 #include "G4RegionStore.hh" 63 #include "G4Gamma.hh" 74 #include "G4Gamma.hh" 64 #include "G4Electron.hh" 75 #include "G4Electron.hh" 65 #include "G4Positron.hh" 76 #include "G4Positron.hh" 66 #include "G4PhysicsTableHelper.hh" 77 #include "G4PhysicsTableHelper.hh" 67 #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 78 77 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 80 79 G4VEmProcess::G4VEmProcess(const G4String& nam 81 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type): 80 G4VDiscreteProcess(name, type) << 82 G4VDiscreteProcess(name, type), >> 83 theLambdaTable(0), >> 84 theEnergyOfCrossSectionMax(0), >> 85 theCrossSectionMax(0), >> 86 particle(0), >> 87 secondaryParticle(0), >> 88 nLambdaBins(90), >> 89 lambdaFactor(0.8), >> 90 currentCouple(0), >> 91 integral(false), >> 92 meanFreePath(true), >> 93 aboveCSmax(true), >> 94 buildLambdaTable(true), >> 95 applyCuts(false), >> 96 startFromNull(true), >> 97 nRegions(0) 81 { 98 { 82 theParameters = G4EmParameters::Instance(); << 83 SetVerboseLevel(1); 99 SetVerboseLevel(1); 84 << 100 minKinEnergy = 0.1*keV; 85 // Size of tables << 101 maxKinEnergy = 100.0*GeV; 86 minKinEnergy = 0.1*CLHEP::keV; << 102 theGamma = G4Gamma::Gamma(); 87 maxKinEnergy = 100.0*CLHEP::TeV; << 103 theElectron = G4Electron::Electron(); 88 << 104 thePositron = G4Positron::Positron(); 89 // default lambda factor << 90 invLambdaFactor = 1.0/lambdaFactor; << 91 << 92 // particle types << 93 theGamma = G4Gamma::Gamma(); << 94 theElectron = G4Electron::Electron(); << 95 thePositron = G4Positron::Positron(); << 96 105 97 pParticleChange = &fParticleChange; 106 pParticleChange = &fParticleChange; 98 fParticleChange.SetSecondaryWeightByProcess( << 99 secParticles.reserve(5); << 100 107 101 modelManager = new G4EmModelManager(); 108 modelManager = new G4EmModelManager(); 102 lManager = G4LossTableManager::Instance(); << 109 (G4LossTableManager::Instance())->Register(this); 103 lManager->Register(this); << 104 isTheMaster = lManager->IsMaster(); << 105 G4LossTableBuilder* bld = lManager->GetTable << 106 theDensityFactor = bld->GetDensityFactors(); << 107 theDensityIdx = bld->GetCoupleIndexes(); << 108 } 110 } 109 111 110 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 111 113 112 G4VEmProcess::~G4VEmProcess() 114 G4VEmProcess::~G4VEmProcess() 113 { 115 { 114 if(isTheMaster) { << 116 Clear(); 115 delete theData; << 117 if(theLambdaTable) theLambdaTable->clearAndDestroy(); 116 delete theEnergyOfCrossSectionMax; << 117 } << 118 delete modelManager; 118 delete modelManager; 119 delete biasManager; << 119 (G4LossTableManager::Instance())->DeRegister(this); 120 lManager->DeRegister(this); << 121 } 120 } 122 121 123 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 123 125 void G4VEmProcess::AddEmModel(G4int order, G4V << 124 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part) 126 const G4Region* << 127 { 125 { 128 if(nullptr == ptr) { return; } << 126 if(!particle) particle = ∂ 129 G4VEmFluctuationModel* fm = nullptr; << 127 if(1 < verboseLevel) { 130 modelManager->AddEmModel(order, ptr, fm, reg << 128 G4cout << "G4VEmProcess::PreparePhysicsTable() for " 131 ptr->SetParticleChange(pParticleChange); << 129 << GetProcessName() 132 } << 130 << " and particle " << part.GetParticleName() 133 << 131 << " local particle " << particle->GetParticleName() 134 //....oooOO0OOooo........oooOO0OOooo........oo << 132 << G4endl; >> 133 } 135 134 136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 135 if(particle == &part) { 137 { << 136 Clear(); 138 if(nullptr == ptr) { return; } << 137 InitialiseProcess(particle); 139 if(!emModels.empty()) { << 138 theCutsGamma = 140 for(auto & em : emModels) { if(em == ptr) << 139 modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel); >> 140 const G4ProductionCutsTable* theCoupleTable= >> 141 G4ProductionCutsTable::GetProductionCutsTable(); >> 142 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut); >> 143 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut); >> 144 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut); >> 145 if(buildLambdaTable) >> 146 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); 141 } 147 } 142 emModels.push_back(ptr); << 143 } 148 } 144 149 145 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 146 151 147 void G4VEmProcess::PreparePhysicsTable(const G << 152 void G4VEmProcess::Clear() 148 { << 153 { 149 if(nullptr == particle) { SetParticle(&part) << 154 if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax; 150 << 155 if(theCrossSectionMax) delete [] theCrossSectionMax; 151 if(part.GetParticleType() == "nucleus" && << 156 theEnergyOfCrossSectionMax = 0; 152 part.GetParticleSubType() == "generic") { << 157 theCrossSectionMax = 0; 153 << 158 modelManager->Clear(); 154 G4String pname = part.GetParticleName(); << 159 currentCouple = 0; 155 if(pname != "deuteron" && pname != "triton << 156 pname != "He3" && pname != "alpha" && p << 157 pname != "helium" && pname != "hydrogen << 158 << 159 particle = G4GenericIon::GenericIon(); << 160 isIon = true; << 161 } << 162 } << 163 if(particle != &part) { return; } << 164 << 165 lManager->PreparePhysicsTable(&part, this); << 166 << 167 // for new run << 168 currentCouple = nullptr; << 169 preStepLambda = 0.0; 160 preStepLambda = 0.0; 170 fLambdaEnergy = 0.0; << 161 mfpKinEnergy = DBL_MAX; 171 << 162 preStepMFP = DBL_MAX; 172 InitialiseProcess(particle); << 173 << 174 G4LossTableBuilder* bld = lManager->GetTable << 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 << 197 if(buildLambdaTable) { << 198 theLambdaTable = theData->MakeTable(0); << 199 bld->InitialiseBaseMaterials(theLambdaTa << 200 } << 201 // high energy table << 202 if(minKinEnergyPrim < maxKinEnergy) { << 203 theLambdaTablePrim = theData->MakeTable( << 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 } << 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 } 163 } 226 164 227 //....oooOO0OOooo........oooOO0OOooo........oo 165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 228 166 229 void G4VEmProcess::BuildPhysicsTable(const G4P 167 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 230 { 168 { 231 if(nullptr == masterProc) { << 169 if(1 < verboseLevel) { 232 if(isTheMaster) { masterProc = this; } << 170 G4cout << "G4VEmProcess::BuildPhysicsTable() for " 233 else { masterProc = static_cast<const G4VE << 171 << GetProcessName() 234 } << 172 << " and particle " << part.GetParticleName() 235 G4int nModels = modelManager->NumberOfModels << 173 << " buildLambdaTable= " << buildLambdaTable 236 G4bool isLocked = theParameters->IsPrintLock << 174 << G4endl; 237 G4bool toBuild = (buildLambdaTable || minKin << 175 } 238 << 239 G4EmTableUtil::BuildEmProcess(this, masterPr << 240 nModels, verbo << 241 isLocked, toBu << 242 } << 243 176 244 //....oooOO0OOooo........oooOO0OOooo........oo << 177 if(buildLambdaTable) { >> 178 BuildLambdaTable(); >> 179 FindLambdaMax(); >> 180 } >> 181 if(0 < verboseLevel) PrintInfoDefinition(); 245 182 246 void G4VEmProcess::BuildLambdaTable() << 183 if(1 < verboseLevel) { 247 { << 184 G4cout << "G4VEmProcess::BuildPhysicsTable() done for " 248 G4double scale = theParameters->MaxKinEnergy << 185 << GetProcessName() 249 G4int nbin = << 186 << " and particle " << part.GetParticleName() 250 theParameters->NumberOfBinsPerDecade()*G4l << 187 << G4endl; 251 if(actBinning) { nbin = std::max(nbin, nLamb << 188 } 252 scale = nbin/G4Log(scale); << 253 << 254 G4LossTableBuilder* bld = lManager->GetTable << 255 G4EmTableUtil::BuildLambdaTable(this, partic << 256 bld, theLamb << 257 minKinEnergy << 258 maxKinEnergy << 259 startFromNul << 260 } 189 } 261 190 262 //....oooOO0OOooo........oooOO0OOooo........oo 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 263 192 264 void G4VEmProcess::StreamInfo(std::ostream& ou << 193 void G4VEmProcess::BuildLambdaTable() 265 const G4ParticleDefinition& << 194 { 266 { << 195 if(1 < verboseLevel) { 267 G4String indent = (rst ? " " : ""); << 196 G4cout << "G4EmProcess::BuildLambdaTable() for process " 268 out << std::setprecision(6); << 197 << GetProcessName() << " and particle " 269 out << G4endl << indent << GetProcessName() << 198 << particle->GetParticleName() 270 if (!rst) { << 199 << G4endl; 271 out << " for " << part.GetParticleName(); << 272 } << 273 if(fXSType != fEmNoIntegral) { out << " XSt << 274 if(applyCuts) { out << " applyCuts:1 "; } << 275 G4int subtype = GetProcessSubType(); << 276 out << " SubType=" << subtype; << 277 if (subtype == fAnnihilation) { << 278 G4int mod = theParameters->PositronAtRestM << 279 const G4String namp[2] = {"Simple", "Allis << 280 out << " AtRestModel:" << namp[mod]; << 281 } << 282 if(biasFactor != 1.0) { out << " BiasingFac << 283 out << " BuildTable=" << buildLambdaTable << << 284 if(buildLambdaTable) { << 285 if(particle == &part) { << 286 for(auto & v : *theLambdaTable) { << 287 if(nullptr != v) { << 288 out << " Lambda table from "; << 289 G4double emin = v->Energy(0); << 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 } << 302 } else { << 303 out << " Used Lambda table of " << 304 << particle->GetParticleName() << G4endl << 305 } << 306 } 200 } 307 if(minKinEnergyPrim < maxKinEnergy) { << 201 308 if(particle == &part) { << 202 // Access to materials 309 for(auto & v : *theLambdaTablePrim) { << 203 const G4ProductionCutsTable* theCoupleTable= 310 if(nullptr != v) { << 204 G4ProductionCutsTable::GetProductionCutsTable(); 311 out << " LambdaPrime table from << 205 size_t numOfCouples = theCoupleTable->GetTableSize(); 312 << G4BestUnit(v->Energy(0),"Ener << 206 for(size_t i=0; i<numOfCouples; i++) { 313 << " to " << 207 314 << G4BestUnit(v->GetMaxEnergy(), << 208 if (theLambdaTable->GetFlag(i)) { 315 << " in " << v->GetVectorLength( << 209 316 << " bins " << G4endl; << 210 // create physics vector and fill it 317 break; << 211 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); 318 } << 212 G4PhysicsVector* aVector = LambdaPhysicsVector(couple); 319 } << 213 modelManager->FillLambdaVector(aVector, couple, startFromNull); 320 } else { << 214 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector); 321 out << " Used LambdaPrime table of << 322 << particle->GetParticleName() << 323 } 215 } 324 } 216 } 325 StreamProcessInfo(out); << 326 modelManager->DumpModelList(out, verboseLeve << 327 217 328 if(verboseLevel > 2 && buildLambdaTable) { << 218 if(1 < verboseLevel) { 329 out << " LambdaTable address= " << th << 219 G4cout << "Lambda table is built for " 330 if(theLambdaTable && particle == &part) { << 220 << particle->GetParticleName() 331 out << (*theLambdaTable) << G4endl; << 221 << G4endl; >> 222 if(2 < verboseLevel) { >> 223 G4cout << *theLambdaTable << G4endl; 332 } 224 } 333 } 225 } 334 } 226 } 335 227 336 //....oooOO0OOooo........oooOO0OOooo........oo 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 337 229 338 void G4VEmProcess::StartTracking(G4Track* trac << 230 void G4VEmProcess::SetParticle(const G4ParticleDefinition* p) 339 { 231 { 340 // reset parameters for the new track << 232 particle = p; 341 currentParticle = track->GetParticleDefiniti << 342 theNumberOfInteractionLengthLeft = -1.0; << 343 mfpKinEnergy = DBL_MAX; << 344 preStepLambda = 0.0; << 345 << 346 if(isIon) { massRatio = proton_mass_c2/curre << 347 << 348 // forced biasing only for primary particles << 349 if(biasManager) { << 350 if(0 == track->GetParentID()) { << 351 // primary particle << 352 biasFlag = true; << 353 biasManager->ResetForcedInteraction(); << 354 } << 355 } << 356 } 233 } 357 234 358 //....oooOO0OOooo........oooOO0OOooo........oo 235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 359 236 360 G4double G4VEmProcess::PostStepGetPhysicalInte << 237 void G4VEmProcess::SetSecondaryParticle(const G4ParticleDefinition* p) 361 const G4Track& tr << 238 { 362 G4double previo << 239 secondaryParticle = p; 363 G4ForceCondition* << 364 { << 365 *condition = NotForced; << 366 G4double x = DBL_MAX; << 367 << 368 DefineMaterial(track.GetMaterialCutsCouple() << 369 preStepKinEnergy = track.GetKineticEnergy(); << 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 << 384 // forced biasing only for primary particles << 385 if(biasManager) { << 386 if(0 == track.GetParentID()) { << 387 if(biasFlag && << 388 biasManager->ForcedInteractionRegion( << 389 return biasManager->GetStepLimit((G4in << 390 } << 391 } << 392 } << 393 << 394 // compute mean free path << 395 << 396 ComputeIntegralLambda(preStepKinEnergy, trac << 397 << 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 << 407 << 408 // beggining of tracking (or just after << 409 theNumberOfInteractionLengthLeft = -G4Lo << 410 theInitialNumberOfInteractionLength = th << 411 << 412 } else { << 413 << 414 theNumberOfInteractionLengthLeft -= << 415 previousStepSize/currentInteractionLen << 416 theNumberOfInteractionLengthLeft = << 417 std::max(theNumberOfInteractionLengthL << 418 } << 419 << 420 // new mean free path and step limit for t << 421 currentInteractionLength = 1.0/preStepLamb << 422 x = theNumberOfInteractionLengthLeft * cur << 423 } << 424 return x; << 425 } 240 } 426 241 427 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 428 243 429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 244 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p, >> 245 const G4Region* region) 430 { 246 { 431 if (fXSType == fEmNoIntegral) { << 247 modelManager->AddEmModel(order, p, 0, region); 432 preStepLambda = GetCurrentLambda(e, LogEki << 248 if(p) p->SetParticleChange(pParticleChange); 433 << 249 } 434 } else if (fXSType == fEmIncreasing) { << 435 if(e*invLambdaFactor < mfpKinEnergy) { << 436 preStepLambda = GetCurrentLambda(e, LogE << 437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 438 } << 439 250 440 } else if(fXSType == fEmDecreasing) { << 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 441 if(e < mfpKinEnergy) { << 442 const G4double e1 = e*lambdaFactor; << 443 preStepLambda = GetCurrentLambda(e1); << 444 mfpKinEnergy = e1; << 445 } << 446 252 447 } else if(fXSType == fEmOnePeak) { << 253 void G4VEmProcess::UpdateEmModel(const G4String& nam, 448 const G4double epeak = (*theEnergyOfCrossS << 254 G4double emin, G4double emax) 449 if(e <= epeak) { << 255 { 450 if(e*invLambdaFactor < mfpKinEnergy) { << 256 modelManager->UpdateEmModel(nam, emin, emax); 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 } << 459 } else { << 460 preStepLambda = GetCurrentLambda(e, LogEki << 461 } << 462 } 257 } 463 258 464 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 465 260 466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 261 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track, 467 262 const G4Step& step) 468 { 263 { 469 // clear number of interaction lengths in an << 470 theNumberOfInteractionLengthLeft = -1.0; << 471 mfpKinEnergy = DBL_MAX; << 472 << 473 fParticleChange.InitializeForPostStep(track) 264 fParticleChange.InitializeForPostStep(track); 474 265 475 // Do not make anything if particle is stopp 266 // Do not make anything if particle is stopped, the annihilation then 476 // should be performed by the AtRestDoIt! 267 // should be performed by the AtRestDoIt! 477 if (track.GetTrackStatus() == fStopButAlive) << 268 if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange; 478 << 479 const G4double finalT = track.GetKineticEner << 480 269 481 // forced process - should happen only once << 270 G4double finalT = track.GetKineticEnergy(); 482 if(biasFlag) { << 483 if(biasManager->ForcedInteractionRegion((G << 484 biasFlag = false; << 485 } << 486 } << 487 << 488 // check active and select model << 489 const G4double scaledEnergy = finalT*massRat << 490 SelectModel(scaledEnergy, currentCoupleIndex << 491 if(!currentModel->IsActive(scaledEnergy)) { << 492 271 493 // Integral approach 272 // Integral approach 494 if (fXSType != fEmNoIntegral) { << 273 if (integral) { 495 const G4double logFinalT = << 274 G4double lx = GetLambda(finalT, currentCouple); 496 track.GetDynamicParticle()->GetLogKineti << 275 if(preStepLambda<lx && 1 < verboseLevel) { 497 const G4double lx = std::max(GetCurrentLam << 276 G4cout << "WARING: for " << particle->GetParticleName() 498 #ifdef G4VERBOSE << 277 << " and " << GetProcessName() 499 if(preStepLambda < lx && 1 < verboseLevel) << 278 << " E(MeV)= " << finalT/MeV 500 G4cout << "WARNING: for " << currentPart << 279 << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) " 501 << " and " << GetProcessName() << << 280 << G4endl; 502 << " preLambda= " << preStepLambd << 503 << " < " << lx << " (postLambda) << 504 } 281 } 505 #endif << 506 // if false interaction then use new cross << 507 // if both values are zero - no interactio << 508 if(preStepLambda*G4UniformRand() >= lx) { << 509 return &fParticleChange; << 510 } << 511 } << 512 282 513 // define new weight for primary and seconda << 283 if(preStepLambda*G4UniformRand() > lx) 514 G4double weight = fParticleChange.GetParentW << 284 return G4VDiscreteProcess::PostStepDoIt(track,step); 515 if(weightFlag) { << 516 weight /= biasFactor; << 517 fParticleChange.ProposeWeight(weight); << 518 } 285 } 519 << 286 520 #ifdef G4VERBOSE << 287 G4VEmModel* currentModel = SelectModel(finalT); 521 if(1 < verboseLevel) { << 288 >> 289 /* >> 290 if(0 < verboseLevel) { 522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 291 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= " 523 << finalT/MeV 292 << finalT/MeV 524 << " MeV; model= (" << currentModel 293 << " MeV; model= (" << currentModel->LowEnergyLimit() 525 << ", " << currentModel->HighEnerg 294 << ", " << currentModel->HighEnergyLimit() << ")" 526 << G4endl; 295 << G4endl; 527 } 296 } 528 #endif << 297 */ 529 << 530 // sample secondaries << 531 secParticles.clear(); << 532 currentModel->SampleSecondaries(&secParticle << 533 currentCoupl << 534 track.GetDyn << 535 (*theCuts)[c << 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 298 554 // save secondaries << 299 555 G4int num = (G4int)secParticles.size(); << 300 std::vector<G4DynamicParticle*>* newp = 556 if(num > 0) { << 301 SecondariesPostStep(currentModel,currentCouple,track.GetDynamicParticle()); 557 302 >> 303 if (newp) { >> 304 G4int num = newp->size(); 558 fParticleChange.SetNumberOfSecondaries(num 305 fParticleChange.SetNumberOfSecondaries(num); 559 G4double edep = fParticleChange.GetLocalEn 306 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 307 567 for (G4int i=0; i<num; ++i) { << 308 for (G4int i=0; i<num; i++) { 568 G4DynamicParticle* dp = secParticles[i]; << 309 G4DynamicParticle* dp = (*newp)[i]; 569 if (nullptr != dp) { << 310 const G4ParticleDefinition* p = dp->GetDefinition(); 570 const G4ParticleDefinition* p = dp->Ge << 311 G4double e = dp->GetKineticEnergy(); 571 G4double e = dp->GetKineticEnergy(); << 312 G4bool good = true; 572 G4bool good = true; << 313 if (p == theGamma) { 573 if(applyCuts) { << 314 if (e < (*theCutsGamma)[currentMaterialIndex]) good = false; 574 if (p == theGamma) { << 315 575 if (e < (*theCutsGamma)[currentCou << 316 } else if (p == theElectron) { 576 << 317 if (e < (*theCutsElectron)[currentMaterialIndex]) good = false; 577 } else if (p == theElectron) { << 318 578 if (e < (*theCutsElectron)[current << 319 } else if (p == thePositron) { 579 << 320 if (e < (*theCutsPositron)[currentMaterialIndex]) { 580 } else if (p == thePositron) { << 321 good = false; 581 if (electron_mass_c2 < (*theCutsGa << 322 e += 2.0*electron_mass_c2; 582 e < (*theCutsPositron)[current << 323 } 583 good = false; << 324 } 584 e += 2.0*electron_mass_c2; << 585 } << 586 } << 587 // added secondary if it is good << 588 } << 589 if (good) { << 590 G4Track* t = new G4Track(dp, time, t << 591 t->SetTouchableHandle(track.GetTouch << 592 if (biasManager) { << 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( << 633 } << 634 325 635 if(0.0 == fParticleChange.GetProposedKinetic << 326 if (good || !applyCuts) { 636 fAlive == fParticleChange.GetTrackStatus( << 327 fParticleChange.AddSecondary(dp); 637 if(particle->GetProcessManager()->GetAtRes << 328 638 { fParticleChange.ProposeTrackStatus( << 329 } else { 639 else { fParticleChange.ProposeTrackStatus( << 330 delete dp; >> 331 edep += e; >> 332 } >> 333 } >> 334 fParticleChange.ProposeLocalEnergyDeposit(edep); >> 335 delete newp; 640 } 336 } 641 337 642 return &fParticleChange; << 338 return G4VDiscreteProcess::PostStepDoIt(track,step); 643 } 339 } 644 340 645 //....oooOO0OOooo........oooOO0OOooo........oo 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 646 342 647 G4bool G4VEmProcess::StorePhysicsTable(const G << 343 void G4VEmProcess::PrintInfoDefinition() 648 const G << 649 G4bool << 650 { 344 { 651 if(!isTheMaster || part != particle) { retur << 345 if(verboseLevel > 0) { 652 if(G4EmTableUtil::StoreTable(this, part, the << 346 G4cout << G4endl << GetProcessName() << ": " ; 653 directory, "Lambda", << 347 PrintInfo(); 654 verboseLevel, a << 348 if(integral) { 655 G4EmTableUtil::StoreTable(this, part, the << 349 G4cout << " Integral mode is used "<< G4endl; 656 directory, "LambdaPrim", << 350 } 657 verboseLevel, a << 351 } 658 return true; << 352 >> 353 if (!buildLambdaTable) return; >> 354 >> 355 if(verboseLevel > 0) { >> 356 G4cout << " tables are built for " >> 357 << particle->GetParticleName() >> 358 << G4endl >> 359 << " Lambda tables from " >> 360 << G4BestUnit(minKinEnergy,"Energy") >> 361 << " to " >> 362 << G4BestUnit(maxKinEnergy,"Energy") >> 363 << " in " << nLambdaBins << " bins." >> 364 << G4endl; >> 365 } >> 366 >> 367 if(verboseLevel > 1) { >> 368 G4cout << "Tables are built for " << particle->GetParticleName() >> 369 << G4endl; >> 370 >> 371 if(verboseLevel > 2) { >> 372 G4cout << "LambdaTable address= " << theLambdaTable << G4endl; >> 373 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl; >> 374 } 659 } 375 } 660 return false; << 661 } 376 } 662 377 663 //....oooOO0OOooo........oooOO0OOooo........oo << 378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 664 379 665 G4bool G4VEmProcess::RetrievePhysicsTable(cons << 380 G4double G4VEmProcess::MicroscopicCrossSection(G4double kineticEnergy, 666 cons << 381 const G4MaterialCutsCouple* couple) 667 G4bo << 668 { 382 { 669 if(!isTheMaster || part != particle) { retur << 383 // Cross section per atom is calculated 670 G4bool yes = true; << 384 DefineMaterial(couple); 671 if(buildLambdaTable) { << 385 G4double cross = 0.0; 672 yes = G4EmTableUtil::RetrieveTable(this, p << 386 G4bool b; 673 "Lambda << 387 if(theLambdaTable) { 674 ascii, << 388 cross = (((*theLambdaTable)[currentMaterialIndex])-> 675 } << 389 GetValue(kineticEnergy, b)); 676 if(yes && minKinEnergyPrim < maxKinEnergy) { << 390 677 yes = G4EmTableUtil::RetrieveTable(this, p << 391 cross /= currentMaterial->GetTotNbOfAtomsPerVolume(); 678 "Lambda << 392 } else { 679 ascii, << 393 G4VEmModel* model = SelectModel(kineticEnergy); >> 394 cross = >> 395 model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy); 680 } 396 } 681 return yes; << 397 >> 398 return cross; 682 } 399 } 683 400 684 //....oooOO0OOooo........oooOO0OOooo........oo 401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 685 402 686 G4double G4VEmProcess::GetCrossSection(G4doubl << 403 G4double G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy, 687 const G << 404 G4double Z, G4double A) 688 { 405 { 689 CurrentSetup(couple, kinEnergy); << 406 G4VEmModel* model = SelectModel(kineticEnergy); 690 return GetCurrentLambda(kinEnergy, G4Log(kin << 407 return model->ComputeCrossSectionPerAtom(particle,kineticEnergy,Z,A); 691 } 408 } 692 409 693 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 694 411 695 G4double G4VEmProcess::GetMeanFreePath(const G << 412 G4double G4VEmProcess::MeanFreePath(const G4Track& track, 696 G4doubl << 413 G4double s, 697 G4Force << 414 G4ForceCondition* cond) 698 { 415 { 699 *condition = NotForced; << 416 return GetMeanFreePath(track, s, cond); 700 return G4VEmProcess::MeanFreePath(track); << 701 } 417 } 702 418 703 //....oooOO0OOooo........oooOO0OOooo........oo 419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 704 420 705 G4double << 421 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part, 706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou << 422 const G4String& directory, 707 G4dou << 423 G4bool ascii) 708 { 424 { 709 SelectModel(kinEnergy, currentCoupleIndex); << 425 G4bool yes = true; 710 return (currentModel) ? << 426 711 currentModel->ComputeCrossSectionPerAtom(c << 427 if ( theLambdaTable && part == particle) { 712 Z << 428 const G4String name = >> 429 GetPhysicsTableFileName(part,directory,"Lambda",ascii); >> 430 yes = theLambdaTable->StorePhysicsTable(name,ascii); >> 431 >> 432 if ( yes ) { >> 433 G4cout << "Physics tables are stored for " << particle->GetParticleName() >> 434 << " and process " << GetProcessName() >> 435 << " in the directory <" << directory >> 436 << "> " << G4endl; >> 437 } else { >> 438 G4cout << "Fail to store Physics Tables for " >> 439 << particle->GetParticleName() >> 440 << " and process " << GetProcessName() >> 441 << " in the directory <" << directory >> 442 << "> " << G4endl; >> 443 } >> 444 } >> 445 return yes; 713 } 446 } 714 447 715 //....oooOO0OOooo........oooOO0OOooo........oo << 448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 716 449 717 G4PhysicsVector* << 450 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 451 const G4String& directory, >> 452 G4bool ascii) 719 { 453 { 720 DefineMaterial(couple); << 454 if(1 < verboseLevel) { 721 G4PhysicsVector* newv = new G4PhysicsLogVect << 455 G4cout << "G4VEmProcess::RetrievePhysicsTable() for " 722 << 456 << part->GetParticleName() << " and process " 723 return newv; << 457 << GetProcessName() << G4endl; >> 458 } >> 459 G4bool yes = true; >> 460 >> 461 if(!buildLambdaTable || particle != part) return yes; >> 462 >> 463 const G4String particleName = part->GetParticleName(); >> 464 G4String filename; >> 465 >> 466 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii); >> 467 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable, >> 468 filename,ascii); >> 469 if ( yes ) { >> 470 if (0 < verboseLevel) { >> 471 G4cout << "Lambda table for " << particleName << " is Retrieved from <" >> 472 << filename << ">" >> 473 << G4endl; >> 474 } >> 475 } else { >> 476 if (1 < verboseLevel) { >> 477 G4cout << "Lambda table for " << particleName << " in file <" >> 478 << filename << "> is not exist" >> 479 << G4endl; >> 480 } >> 481 } >> 482 >> 483 return yes; 724 } 484 } 725 485 726 //....oooOO0OOooo........oooOO0OOooo........oo 486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 727 487 728 const G4Element* G4VEmProcess::GetCurrentEleme << 488 void G4VEmProcess::FindLambdaMax() 729 { 489 { 730 return (nullptr != currentModel) ? << 490 if(1 < verboseLevel) { 731 currentModel->GetCurrentElement(currentMat << 491 G4cout << "### G4VEmProcess::FindLambdaMax: " >> 492 << particle->GetParticleName() >> 493 << " and process " << GetProcessName() << G4endl; >> 494 } >> 495 size_t n = theLambdaTable->length(); >> 496 G4PhysicsVector* pv = (*theLambdaTable)[0]; >> 497 G4double e, s, emax, smax; >> 498 theEnergyOfCrossSectionMax = new G4double [n]; >> 499 theCrossSectionMax = new G4double [n]; >> 500 G4bool b; >> 501 >> 502 for (size_t i=0; i<n; i++) { >> 503 pv = (*theLambdaTable)[i]; >> 504 emax = DBL_MAX; >> 505 smax = 0.0; >> 506 if(pv) { >> 507 size_t nb = pv->GetVectorLength(); >> 508 emax = pv->GetLowEdgeEnergy(nb); >> 509 smax = 0.0; >> 510 for (size_t j=0; j<nb; j++) { >> 511 e = pv->GetLowEdgeEnergy(j); >> 512 s = pv->GetValue(e,b); >> 513 if(s > smax) { >> 514 smax = s; >> 515 emax = e; >> 516 } >> 517 } >> 518 } >> 519 theEnergyOfCrossSectionMax[i] = emax; >> 520 theCrossSectionMax[i] = smax; >> 521 if(2 < verboseLevel) { >> 522 G4cout << "For " << particle->GetParticleName() >> 523 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV >> 524 << " lambda= " << smax << G4endl; >> 525 } >> 526 } 732 } 527 } 733 528 734 //....oooOO0OOooo........oooOO0OOooo........oo 529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 735 530 736 const G4Element* G4VEmProcess::GetTargetElemen << 531 void G4VEmProcess::SetLambdaBinning(G4int nbins) 737 { 532 { 738 return (nullptr != currentModel) ? << 533 nLambdaBins = nbins; 739 currentModel->GetCurrentElement(currentMat << 740 } 534 } 741 535 742 //....oooOO0OOooo........oooOO0OOooo........oo 536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 743 537 744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 538 G4int G4VEmProcess::LambdaBinning() const 745 { 539 { 746 return (nullptr != currentModel) ? << 540 return nLambdaBins; 747 currentModel->GetCurrentIsotope(GetCurrent << 748 } 541 } 749 542 750 //....oooOO0OOooo........oooOO0OOooo........oo 543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 751 544 752 void G4VEmProcess::SetCrossSectionBiasingFacto << 545 void G4VEmProcess::SetMinKinEnergy(G4double e) 753 { 546 { 754 if(f > 0.0) { << 547 minKinEnergy = e; 755 biasFactor = f; << 756 weightFlag = flag; << 757 if(1 < verboseLevel) { << 758 G4cout << "### SetCrossSectionBiasingFac << 759 << particle->GetParticleName() << 760 << " and process " << GetProcessN << 761 << " biasFactor= " << f << " weig << 762 << G4endl; << 763 } << 764 } << 765 } 548 } 766 549 767 //....oooOO0OOooo........oooOO0OOooo........oo 550 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 768 551 769 void << 552 G4double G4VEmProcess::MinKinEnergy() const 770 G4VEmProcess::ActivateForcedInteraction(G4doub << 771 G4bool << 772 { 553 { 773 if(nullptr == biasManager) { biasManager = n << 554 return minKinEnergy; 774 if(1 < verboseLevel) { << 775 G4cout << "### ActivateForcedInteraction: << 776 << particle->GetParticleName() << 777 << " and process " << GetProcessNam << 778 << " length(mm)= " << length/mm << 779 << " in G4Region <" << r << 780 << "> weightFlag= " << flag << 781 << G4endl; << 782 } << 783 weightFlag = flag; << 784 biasManager->ActivateForcedInteraction(lengt << 785 } 555 } 786 556 787 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 788 558 789 void << 559 void G4VEmProcess::SetMaxKinEnergy(G4double e) 790 G4VEmProcess::ActivateSecondaryBiasing(const G << 560 { 791 G4double factor, << 561 maxKinEnergy = e; 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 } 562 } 812 563 813 //....oooOO0OOooo........oooOO0OOooo........oo 564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 814 565 815 void G4VEmProcess::SetLambdaBinning(G4int n) << 566 G4double G4VEmProcess::MaxKinEnergy() const 816 { 567 { 817 if(5 < n && n < 10000000) { << 568 return maxKinEnergy; 818 nLambdaBins = n; << 819 actBinning = true; << 820 } else { << 821 G4double e = (G4double)n; << 822 PrintWarning("SetLambdaBinning", e); << 823 } << 824 } 569 } 825 570 826 //....oooOO0OOooo........oooOO0OOooo........oo 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 827 572 828 void G4VEmProcess::SetMinKinEnergy(G4double e) << 573 void G4VEmProcess::ActivateDeexcitation(G4bool, const G4Region*) >> 574 {} >> 575 >> 576 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 577 >> 578 const G4PhysicsTable* G4VEmProcess::LambdaTable() const 829 { 579 { 830 if(1.e-3*eV < e && e < maxKinEnergy) { << 580 return theLambdaTable; 831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma << 832 /G4Log(maxKinEnergy/ << 833 minKinEnergy = e; << 834 actMinKinEnergy = true; << 835 } else { PrintWarning("SetMinKinEnergy", e); << 836 } 581 } 837 582 838 //....oooOO0OOooo........oooOO0OOooo........oo 583 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 839 584 840 void G4VEmProcess::SetMaxKinEnergy(G4double e) << 585 void G4VEmProcess::SetIntegral(G4bool val) 841 { 586 { 842 if(minKinEnergy < e && e < 1.e+6*TeV) { << 587 if(particle && particle != theGamma) integral = val; 843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ << 588 if(integral) buildLambdaTable = true; 844 /G4Log(maxKinEnergy/ << 845 maxKinEnergy = e; << 846 actMaxKinEnergy = true; << 847 } else { PrintWarning("SetMaxKinEnergy", e); << 848 } 589 } 849 590 850 //....oooOO0OOooo........oooOO0OOooo........oo 591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 851 592 852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl << 593 G4bool G4VEmProcess::IsIntegral() const 853 { 594 { 854 if(theParameters->MinKinEnergy() <= e && << 595 return integral; 855 e <= theParameters->MaxKinEnergy()) { min << 856 else { PrintWarning("SetMinKinEnergyPrim", e << 857 } 596 } 858 597 859 //....oooOO0OOooo........oooOO0OOooo........oo 598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 860 599 861 G4VEmProcess* G4VEmProcess::GetEmProcess(const << 600 G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*) 862 { 601 { 863 return (nam == GetProcessName()) ? this : nu << 602 G4PhysicsVector* v = >> 603 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins); >> 604 return v; 864 } 605 } 865 606 866 //....oooOO0OOooo........oooOO0OOooo........oo 607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 867 608 868 G4double G4VEmProcess::PolarAngleLimit() const << 609 void G4VEmProcess::SetBuildTableFlag(G4bool val) 869 { 610 { 870 return theParameters->MscThetaLimit(); << 611 buildLambdaTable = val; 871 } 612 } 872 613 873 //....oooOO0OOooo........oooOO0OOooo........oo 614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 874 615 875 void G4VEmProcess::PrintWarning(G4String tit, << 616 void G4VEmProcess::SetStartFromNullFlag(G4bool val) 876 { 617 { 877 G4String ss = "G4VEmProcess::" + tit; << 618 startFromNull = val; 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 } 619 } 886 620 887 //....oooOO0OOooo........oooOO0OOooo........oo 621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 888 622 889 void G4VEmProcess::ProcessDescription(std::ost << 623 void G4VEmProcess::SetApplyCuts(G4bool val) 890 { 624 { 891 if(nullptr != particle) { << 625 applyCuts = val; 892 StreamInfo(out, *particle, true); << 893 } << 894 } 626 } 895 627 896 //....oooOO0OOooo........oooOO0OOooo........oo 628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 897 629