Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // 26 // ------------------------------------------- 27 // ------------------------------------------------------------------- 27 // 28 // 28 // GEANT4 Class file 29 // GEANT4 Class file 29 // 30 // 30 // 31 // 31 // File name: G4VEnergyLossProcess 32 // File name: G4VEnergyLossProcess 32 // 33 // 33 // Author: Vladimir Ivanchenko 34 // Author: Vladimir Ivanchenko 34 // 35 // 35 // Creation date: 03.01.2002 36 // Creation date: 03.01.2002 36 // 37 // 37 // Modifications: Vladimir Ivanchenko 38 // Modifications: Vladimir Ivanchenko 38 // 39 // 39 // 40 // 40 // Class Description: 41 // Class Description: 41 // 42 // 42 // It is the unified energy loss process it ca 43 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 44 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 45 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 46 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 47 // that information on fly. 47 // ------------------------------------------- 48 // ------------------------------------------------------------------- 48 // 49 // 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 52 52 #include "G4VEnergyLossProcess.hh" 53 #include "G4VEnergyLossProcess.hh" 53 #include "G4PhysicalConstants.hh" 54 #include "G4PhysicalConstants.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "G4SystemOfUnits.hh" 55 #include "G4ProcessManager.hh" 56 #include "G4ProcessManager.hh" 56 #include "G4LossTableManager.hh" 57 #include "G4LossTableManager.hh" 57 #include "G4LossTableBuilder.hh" 58 #include "G4LossTableBuilder.hh" 58 #include "G4Step.hh" 59 #include "G4Step.hh" 59 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" 61 #include "G4ParticleTable.hh" 61 #include "G4EmParameters.hh" << 62 #include "G4EmUtility.hh" << 63 #include "G4EmTableUtil.hh" << 64 #include "G4VEmModel.hh" 62 #include "G4VEmModel.hh" 65 #include "G4VEmFluctuationModel.hh" 63 #include "G4VEmFluctuationModel.hh" 66 #include "G4DataVector.hh" 64 #include "G4DataVector.hh" 67 #include "G4PhysicsLogVector.hh" 65 #include "G4PhysicsLogVector.hh" 68 #include "G4VParticleChange.hh" 66 #include "G4VParticleChange.hh" >> 67 #include "G4Gamma.hh" 69 #include "G4Electron.hh" 68 #include "G4Electron.hh" >> 69 #include "G4Positron.hh" 70 #include "G4ProcessManager.hh" 70 #include "G4ProcessManager.hh" 71 #include "G4UnitsTable.hh" 71 #include "G4UnitsTable.hh" >> 72 #include "G4ProductionCutsTable.hh" 72 #include "G4Region.hh" 73 #include "G4Region.hh" 73 #include "G4RegionStore.hh" 74 #include "G4RegionStore.hh" 74 #include "G4PhysicsTableHelper.hh" 75 #include "G4PhysicsTableHelper.hh" 75 #include "G4SafetyHelper.hh" 76 #include "G4SafetyHelper.hh" 76 #include "G4EmDataHandler.hh" << 77 #include "G4TransportationManager.hh" 77 #include "G4TransportationManager.hh" >> 78 #include "G4EmConfigurator.hh" 78 #include "G4VAtomDeexcitation.hh" 79 #include "G4VAtomDeexcitation.hh" 79 #include "G4VSubCutProducer.hh" 80 #include "G4VSubCutProducer.hh" 80 #include "G4EmBiasingManager.hh" 81 #include "G4EmBiasingManager.hh" 81 #include "G4Log.hh" 82 #include "G4Log.hh" 82 #include <iostream> 83 #include <iostream> 83 84 84 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 86 86 namespace << 87 { << 88 G4String tnames[7] = << 89 {"DEDX","Ionisation","DEDXnr","CSDARange", << 90 } << 91 << 92 << 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con 87 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 94 G4P 88 G4ProcessType type): 95 G4VContinuousDiscreteProcess(name, type) << 89 G4VContinuousDiscreteProcess(name, type), >> 90 secondaryParticle(nullptr), >> 91 nSCoffRegions(0), >> 92 idxSCoffRegions(nullptr), >> 93 nProcesses(0), >> 94 theDEDXTable(nullptr), >> 95 theDEDXSubTable(nullptr), >> 96 theDEDXunRestrictedTable(nullptr), >> 97 theIonisationTable(nullptr), >> 98 theIonisationSubTable(nullptr), >> 99 theRangeTableForLoss(nullptr), >> 100 theCSDARangeTable(nullptr), >> 101 theSecondaryRangeTable(nullptr), >> 102 theInverseRangeTable(nullptr), >> 103 theLambdaTable(nullptr), >> 104 theSubLambdaTable(nullptr), >> 105 theDensityFactor(nullptr), >> 106 theDensityIdx(nullptr), >> 107 baseParticle(nullptr), >> 108 lossFluctuationFlag(true), >> 109 rndmStepFlag(false), >> 110 tablesAreBuilt(false), >> 111 integral(true), >> 112 isIon(false), >> 113 isIonisation(true), >> 114 useSubCutoff(false), >> 115 useDeexcitation(false), >> 116 currentCouple(nullptr), >> 117 mfpKinEnergy(0.0), >> 118 particle(nullptr) 96 { 119 { 97 theParameters = G4EmParameters::Instance(); 120 theParameters = G4EmParameters::Instance(); 98 SetVerboseLevel(1); 121 SetVerboseLevel(1); 99 122 100 // low energy limit 123 // low energy limit 101 lowestKinEnergy = theParameters->LowestElect << 124 lowestKinEnergy = theParameters->LowestElectronEnergy(); 102 << 125 preStepKinEnergy = 0.0; 103 // Size of tables << 126 preStepRangeEnergy = 0.0; 104 minKinEnergy = 0.1*CLHEP::keV; << 127 computedRange = DBL_MAX; 105 maxKinEnergy = 100.0*CLHEP::TeV; << 128 106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; << 129 // Size of tables assuming spline >> 130 minKinEnergy = 0.1*keV; >> 131 maxKinEnergy = 100.0*TeV; 107 nBins = 84; 132 nBins = 84; >> 133 maxKinEnergyCSDA = 1.0*GeV; 108 nBinsCSDA = 35; 134 nBinsCSDA = 35; >> 135 actMinKinEnergy = actMaxKinEnergy = actBinning = actLinLossLimit >> 136 = actLossFluc = actIntegral = actStepFunc = false; 109 137 110 invLambdaFactor = 1.0/lambdaFactor; << 138 // default linear loss limit for spline 111 << 139 linLossLimit = 0.01; 112 // default linear loss limit << 140 dRoverRange = 0.2; 113 finalRange = 1.*CLHEP::mm; << 141 finalRange = CLHEP::mm; >> 142 >> 143 // default lambda factor >> 144 lambdaFactor = 0.8; >> 145 >> 146 // cross section biasing >> 147 biasFactor = 1.0; >> 148 >> 149 // particle types >> 150 theElectron = G4Electron::Electron(); >> 151 thePositron = G4Positron::Positron(); >> 152 theGamma = G4Gamma::Gamma(); >> 153 theGenericIon = nullptr; 114 154 115 // run time objects 155 // run time objects 116 pParticleChange = &fParticleChange; 156 pParticleChange = &fParticleChange; 117 fParticleChange.SetSecondaryWeightByProcess( 157 fParticleChange.SetSecondaryWeightByProcess(true); 118 modelManager = new G4EmModelManager(); 158 modelManager = new G4EmModelManager(); 119 safetyHelper = G4TransportationManager::GetT 159 safetyHelper = G4TransportationManager::GetTransportationManager() 120 ->GetSafetyHelper(); 160 ->GetSafetyHelper(); 121 aGPILSelection = CandidateForSelection; 161 aGPILSelection = CandidateForSelection; 122 162 123 // initialise model 163 // initialise model 124 lManager = G4LossTableManager::Instance(); 164 lManager = G4LossTableManager::Instance(); 125 lManager->Register(this); 165 lManager->Register(this); 126 isMaster = lManager->IsMaster(); << 166 fluctModel = nullptr; 127 << 167 currentModel = nullptr; 128 G4LossTableBuilder* bld = lManager->GetTable << 168 atomDeexcitation = nullptr; 129 theDensityFactor = bld->GetDensityFactors(); << 169 subcutProducer = nullptr; 130 theDensityIdx = bld->GetCoupleIndexes(); << 170 >> 171 biasManager = nullptr; >> 172 biasFlag = false; >> 173 weightFlag = false; >> 174 isMaster = true; >> 175 lastIdx = 0; >> 176 >> 177 idxDEDX = idxDEDXSub = idxDEDXunRestricted = idxIonisation = >> 178 idxIonisationSub = idxRange = idxCSDA = idxSecRange = >> 179 idxInverseRange = idxLambda = idxSubLambda = 0; >> 180 >> 181 scTracks.reserve(5); >> 182 secParticles.reserve(5); >> 183 >> 184 theCuts = theSubCuts = nullptr; >> 185 currentMaterial = nullptr; >> 186 currentCoupleIndex = basedCoupleIndex = 0; >> 187 massRatio = fFactor = reduceFactor = chargeSqRatio = 1.0; >> 188 preStepLambda = preStepScaledEnergy = fRange = 0.0; 131 189 132 scTracks.reserve(10); << 190 secID = biasID = subsecID = -1; 133 secParticles.reserve(12); << 134 emModels = new std::vector<G4VEmModel*>; << 135 } 191 } 136 192 137 //....oooOO0OOooo........oooOO0OOooo........oo 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 194 139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 195 G4VEnergyLossProcess::~G4VEnergyLossProcess() 140 { 196 { 141 if (isMaster) { << 197 /* 142 if(nullptr == baseParticle) { delete theDa << 198 G4cout << "** G4VEnergyLossProcess::~G4VEnergyLossProcess() for " 143 delete theEnergyOfCrossSectionMax; << 199 << GetProcessName() << " isMaster: " << isMaster 144 if(nullptr != fXSpeaks) { << 200 << " basePart: " << baseParticle 145 for(auto const & v : *fXSpeaks) { delete << 201 << G4endl; 146 delete fXSpeaks; << 202 */ >> 203 Clean(); >> 204 >> 205 // G4cout << " isIonisation " << isIonisation << " " >> 206 // << theDEDXTable << " " << theIonisationTable << G4endl; >> 207 >> 208 if (isMaster && !baseParticle) { >> 209 if(theDEDXTable) { >> 210 >> 211 //G4cout << " theIonisationTable " << theIonisationTable << G4endl; >> 212 if(theIonisationTable == theDEDXTable) { theIonisationTable = nullptr; } >> 213 //G4cout << " delete theDEDXTable " << theDEDXTable << G4endl; >> 214 theDEDXTable->clearAndDestroy(); >> 215 delete theDEDXTable; >> 216 theDEDXTable = nullptr; >> 217 if(theDEDXSubTable) { >> 218 if(theIonisationSubTable == theDEDXSubTable) >> 219 { theIonisationSubTable = nullptr; } >> 220 theDEDXSubTable->clearAndDestroy(); >> 221 delete theDEDXSubTable; >> 222 theDEDXSubTable = nullptr; >> 223 } >> 224 } >> 225 //G4cout << " theIonisationTable " << theIonisationTable << G4endl; >> 226 if(theIonisationTable) { >> 227 //G4cout << " delete theIonisationTable " << theIonisationTable << G4endl; >> 228 theIonisationTable->clearAndDestroy(); >> 229 delete theIonisationTable; >> 230 theIonisationTable = nullptr; >> 231 } >> 232 if(theIonisationSubTable) { >> 233 theIonisationSubTable->clearAndDestroy(); >> 234 delete theIonisationSubTable; >> 235 theIonisationSubTable = nullptr; >> 236 } >> 237 if(theDEDXunRestrictedTable && isIonisation) { >> 238 theDEDXunRestrictedTable->clearAndDestroy(); >> 239 delete theDEDXunRestrictedTable; >> 240 theDEDXunRestrictedTable = nullptr; >> 241 } >> 242 if(theCSDARangeTable && isIonisation) { >> 243 theCSDARangeTable->clearAndDestroy(); >> 244 delete theCSDARangeTable; >> 245 theCSDARangeTable = nullptr; >> 246 } >> 247 //G4cout << "delete RangeTable: " << theRangeTableForLoss << G4endl; >> 248 if(theRangeTableForLoss && isIonisation) { >> 249 theRangeTableForLoss->clearAndDestroy(); >> 250 delete theRangeTableForLoss; >> 251 theRangeTableForLoss = nullptr; >> 252 } >> 253 //G4cout << "delete InvRangeTable: " << theInverseRangeTable << G4endl; >> 254 if(theInverseRangeTable && isIonisation /*&& !isIon*/) { >> 255 theInverseRangeTable->clearAndDestroy(); >> 256 delete theInverseRangeTable; >> 257 theInverseRangeTable = nullptr; >> 258 } >> 259 //G4cout << "delete LambdaTable: " << theLambdaTable << G4endl; >> 260 if(theLambdaTable) { >> 261 theLambdaTable->clearAndDestroy(); >> 262 delete theLambdaTable; >> 263 theLambdaTable = nullptr; >> 264 } >> 265 if(theSubLambdaTable) { >> 266 theSubLambdaTable->clearAndDestroy(); >> 267 delete theSubLambdaTable; >> 268 theSubLambdaTable = nullptr; 147 } 269 } 148 } 270 } >> 271 149 delete modelManager; 272 delete modelManager; 150 delete biasManager; 273 delete biasManager; 151 delete scoffRegions; << 152 delete emModels; << 153 lManager->DeRegister(this); 274 lManager->DeRegister(this); >> 275 //G4cout << "** all removed" << G4endl; >> 276 } >> 277 >> 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 279 >> 280 void G4VEnergyLossProcess::Clean() >> 281 { >> 282 /* >> 283 if(1 < verboseLevel) { >> 284 G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName() >> 285 << G4endl; >> 286 } >> 287 */ >> 288 delete [] idxSCoffRegions; >> 289 >> 290 tablesAreBuilt = false; >> 291 >> 292 scProcesses.clear(); >> 293 nProcesses = 0; >> 294 >> 295 idxDEDX = idxDEDXSub = idxDEDXunRestricted = idxIonisation = >> 296 idxIonisationSub = idxRange = idxCSDA = idxSecRange = >> 297 idxInverseRange = idxLambda = idxSubLambda = 0; 154 } 298 } 155 299 156 //....oooOO0OOooo........oooOO0OOooo........oo 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 301 158 G4double G4VEnergyLossProcess::MinPrimaryEnerg 302 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*, 159 303 const G4Material*, 160 304 G4double cut) 161 { 305 { 162 return cut; 306 return cut; 163 } 307 } 164 308 165 //....oooOO0OOooo........oooOO0OOooo........oo 309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 166 310 167 void G4VEnergyLossProcess::AddEmModel(G4int or << 311 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p, 168 G4VEmFlu 312 G4VEmFluctuationModel* fluc, 169 const G4 313 const G4Region* region) 170 { 314 { 171 if(nullptr == ptr) { return; } << 315 modelManager->AddEmModel(order, p, fluc, region); 172 G4VEmFluctuationModel* afluc = (nullptr == f << 316 if(p) { p->SetParticleChange(pParticleChange, fluc); } 173 modelManager->AddEmModel(order, ptr, afluc, << 317 } 174 ptr->SetParticleChange(pParticleChange, aflu << 318 >> 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 320 >> 321 void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam, >> 322 G4double emin, G4double emax) >> 323 { >> 324 modelManager->UpdateEmModel(nam, emin, emax); 175 } 325 } 176 326 177 //....oooOO0OOooo........oooOO0OOooo........oo 327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 328 179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod 329 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* ptr, G4int) 180 { 330 { 181 if(nullptr == ptr) { return; } << 331 for(auto & em : emModels) { if(em == ptr) { return; } } 182 if(!emModels->empty()) { << 332 emModels.push_back(ptr); 183 for(auto & em : *emModels) { if(em == ptr) << 333 } 184 } << 334 185 emModels->push_back(ptr); << 335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 336 >> 337 G4VEmModel* G4VEnergyLossProcess::EmModel(size_t index) const >> 338 { >> 339 return (index < emModels.size()) ? emModels[index] : nullptr; 186 } 340 } 187 341 188 //....oooOO0OOooo........oooOO0OOooo........oo 342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 189 343 190 void G4VEnergyLossProcess::SetDynamicMassCharg << 344 G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver) const 191 << 192 { 345 { 193 massRatio = massratio; << 346 return modelManager->GetModel(idx, ver); 194 logMassRatio = G4Log(massRatio); << 347 } 195 fFactor = charge2ratio*biasFactor; << 348 196 if(baseMat) { fFactor *= (*theDensityFactor) << 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 197 chargeSqRatio = charge2ratio; << 350 198 reduceFactor = 1.0/(fFactor*massRatio); << 351 G4int G4VEnergyLossProcess::NumberOfModels() const >> 352 { >> 353 return modelManager->NumberOfModels(); 199 } 354 } 200 355 201 //....oooOO0OOooo........oooOO0OOooo........oo 356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 202 357 203 void 358 void 204 G4VEnergyLossProcess::PreparePhysicsTable(cons 359 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part) 205 { 360 { 206 particle = G4EmTableUtil::CheckIon(this, &pa << 361 if(1 < verboseLevel) { 207 verboseLe << 362 G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for " >> 363 << GetProcessName() << " for " << part.GetParticleName() >> 364 << " " << this << G4endl; >> 365 } >> 366 isMaster = lManager->IsMaster(); >> 367 >> 368 currentCouple = nullptr; >> 369 preStepLambda = 0.0; >> 370 mfpKinEnergy = DBL_MAX; >> 371 fRange = DBL_MAX; >> 372 preStepKinEnergy = 0.0; >> 373 preStepRangeEnergy = 0.0; >> 374 chargeSqRatio = 1.0; >> 375 massRatio = 1.0; >> 376 reduceFactor = 1.0; >> 377 fFactor = 1.0; >> 378 lastIdx = 0; >> 379 >> 380 // Are particle defined? >> 381 if( !particle ) { particle = ∂ } >> 382 >> 383 if(part.GetParticleType() == "nucleus") { >> 384 >> 385 G4String pname = part.GetParticleName(); >> 386 if(pname != "deuteron" && pname != "triton" && >> 387 pname != "alpha+" && pname != "helium" && >> 388 pname != "hydrogen") { >> 389 >> 390 if(!theGenericIon) { >> 391 theGenericIon = >> 392 G4ParticleTable::GetParticleTable()->FindParticle("GenericIon"); >> 393 } >> 394 isIon = true; >> 395 if(theGenericIon && particle != theGenericIon) { >> 396 G4ProcessManager* pm = theGenericIon->GetProcessManager(); >> 397 G4ProcessVector* v = pm->GetAlongStepProcessVector(); >> 398 size_t n = v->size(); >> 399 for(size_t j=0; j<n; ++j) { >> 400 if((*v)[j] == this) { >> 401 particle = theGenericIon; >> 402 break; >> 403 } >> 404 } >> 405 } >> 406 } >> 407 } 208 408 209 if( particle != &part ) { 409 if( particle != &part ) { 210 if(!isIon) { lManager->RegisterExtraPartic << 410 if(!isIon) { >> 411 lManager->RegisterExtraParticle(&part, this); >> 412 } 211 if(1 < verboseLevel) { 413 if(1 < verboseLevel) { 212 G4cout << "### G4VEnergyLossProcess::Pre 414 G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()" 213 << " interrupted for " << GetProc << 415 << " interrupted for " 214 << part.GetParticleName() << " is << 416 << part.GetParticleName() << " isIon= " << isIon 215 << " spline=" << spline << G4endl << 417 << " particle " << particle << " GenericIon " << theGenericIon >> 418 << G4endl; 216 } 419 } 217 return; 420 return; 218 } 421 } 219 422 220 tablesAreBuilt = false; << 423 Clean(); 221 if (GetProcessSubType() == fIonisation) { Se << 424 lManager->PreparePhysicsTable(&part, this, isMaster); 222 << 223 G4LossTableBuilder* bld = lManager->GetTable 425 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 224 lManager->PreparePhysicsTable(&part, this); << 225 426 226 // Base particle and set of models can be de 427 // Base particle and set of models can be defined here 227 InitialiseEnergyLossProcess(particle, basePa 428 InitialiseEnergyLossProcess(particle, baseParticle); 228 429 >> 430 const G4ProductionCutsTable* theCoupleTable= >> 431 G4ProductionCutsTable::GetProductionCutsTable(); >> 432 size_t n = theCoupleTable->GetTableSize(); >> 433 >> 434 theDEDXAtMaxEnergy.resize(n, 0.0); >> 435 theRangeAtMaxEnergy.resize(n, 0.0); >> 436 theEnergyOfCrossSectionMax.resize(n, 0.0); >> 437 theCrossSectionMax.resize(n, DBL_MAX); >> 438 229 // parameters of the process 439 // parameters of the process >> 440 if(!actIntegral) { integral = theParameters->Integral(); } 230 if(!actLossFluc) { lossFluctuationFlag = the 441 if(!actLossFluc) { lossFluctuationFlag = theParameters->LossFluctuation(); } 231 useCutAsFinalRange = theParameters->UseCutAs << 442 rndmStepFlag = theParameters->UseCutAsFinalRange(); 232 if(!actMinKinEnergy) { minKinEnergy = thePar 443 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); } 233 if(!actMaxKinEnergy) { maxKinEnergy = thePar 444 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); } 234 if(!actBinning) { nBins = theParameters->Num << 445 if(!actBinning) { >> 446 nBins = theParameters->NumberOfBinsPerDecade() >> 447 *G4lrint(std::log10(maxKinEnergy/minKinEnergy)); >> 448 } 235 maxKinEnergyCSDA = theParameters->MaxEnergyF 449 maxKinEnergyCSDA = theParameters->MaxEnergyForCSDARange(); 236 nBinsCSDA = theParameters->NumberOfBinsPerDe 450 nBinsCSDA = theParameters->NumberOfBinsPerDecade() 237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi 451 *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy)); 238 if(!actLinLossLimit) { linLossLimit = thePar 452 if(!actLinLossLimit) { linLossLimit = theParameters->LinearLossLimit(); } 239 lambdaFactor = theParameters->LambdaFactor() 453 lambdaFactor = theParameters->LambdaFactor(); 240 invLambdaFactor = 1.0/lambdaFactor; << 241 if(isMaster) { SetVerboseLevel(theParameters 454 if(isMaster) { SetVerboseLevel(theParameters->Verbose()); } 242 else { SetVerboseLevel(theParameters->Worker << 455 else { SetVerboseLevel(theParameters->WorkerVerbose()); } 243 // integral option may be disabled << 244 if(!theParameters->Integral()) { fXSType = f << 245 456 246 theParameters->DefineRegParamForLoss(this); << 457 G4bool isElec = true; 247 << 458 if(particle->GetPDGMass() > CLHEP::MeV) { isElec = false; } 248 fRangeEnergy = 0.0; << 459 theParameters->DefineRegParamForLoss(this, isElec); 249 460 250 G4double initialCharge = particle->GetPDGCha 461 G4double initialCharge = particle->GetPDGCharge(); 251 G4double initialMass = particle->GetPDGMas 462 G4double initialMass = particle->GetPDGMass(); 252 463 253 theParameters->FillStepFunction(particle, th << 464 if (baseParticle) { 254 << 465 massRatio = (baseParticle->GetPDGMass())/initialMass; 255 // parameters for scaling from the base part << 256 if (nullptr != baseParticle) { << 257 massRatio = (baseParticle->GetPDGMass() << 258 logMassRatio = G4Log(massRatio); << 259 G4double q = initialCharge/baseParticle->G 466 G4double q = initialCharge/baseParticle->GetPDGCharge(); 260 chargeSqRatio = q*q; 467 chargeSqRatio = q*q; 261 if(chargeSqRatio > 0.0) { reduceFactor = 1 468 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); } 262 } 469 } 263 lowestKinEnergy = (initialMass < CLHEP::MeV) << 470 lowestKinEnergy = (initialMass < MeV) ? theParameters->LowestElectronEnergy() 264 ? theParameters->LowestElectronEnergy() << 265 : theParameters->LowestMuHadEnergy(); 471 : theParameters->LowestMuHadEnergy(); 266 472 267 // Tables preparation 473 // Tables preparation 268 if (isMaster && nullptr == baseParticle) { << 474 if (isMaster && !baseParticle) { 269 if(nullptr == theData) { theData = new G4E << 270 475 271 if(nullptr != theDEDXTable && isIonisation << 476 if(theDEDXTable && isIonisation) { 272 if(nullptr != theIonisationTable && theD << 477 if(theIonisationTable && theDEDXTable != theIonisationTable) { 273 theData->CleanTable(0); << 478 theDEDXTable->clearAndDestroy(); 274 theDEDXTable = theIonisationTable; << 479 delete theDEDXTable; 275 theIonisationTable = nullptr; << 480 theDEDXTable = theIonisationTable; 276 } << 481 } >> 482 if(theDEDXSubTable && theIonisationSubTable && >> 483 theDEDXSubTable != theIonisationSubTable) { >> 484 theDEDXSubTable->clearAndDestroy(); >> 485 delete theDEDXSubTable; >> 486 theDEDXSubTable = theIonisationSubTable; >> 487 } 277 } 488 } 278 489 279 theDEDXTable = theData->MakeTable(theDEDXT << 490 theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable); 280 bld->InitialiseBaseMaterials(theDEDXTable) 491 bld->InitialiseBaseMaterials(theDEDXTable); 281 theData->UpdateTable(theIonisationTable, 1 << 492 >> 493 if(theDEDXSubTable) { >> 494 theDEDXSubTable = >> 495 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable); >> 496 } 282 497 283 if (theParameters->BuildCSDARange()) { 498 if (theParameters->BuildCSDARange()) { 284 theDEDXunRestrictedTable = theData->Make << 499 theDEDXunRestrictedTable = 285 if(isIonisation) { theCSDARangeTable = t << 500 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable); >> 501 theCSDARangeTable = >> 502 G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable); 286 } 503 } 287 504 288 theLambdaTable = theData->MakeTable(4); << 505 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); >> 506 289 if(isIonisation) { 507 if(isIonisation) { 290 theRangeTableForLoss = theData->MakeTabl << 508 theRangeTableForLoss = 291 theInverseRangeTable = theData->MakeTabl << 509 G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss); >> 510 theInverseRangeTable = >> 511 G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable); 292 } 512 } 293 } << 294 513 >> 514 if (nSCoffRegions && !lManager->SubCutProducer()) { >> 515 theDEDXSubTable = >> 516 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable); >> 517 theSubLambdaTable = >> 518 G4PhysicsTableHelper::PreparePhysicsTable(theSubLambdaTable); >> 519 } >> 520 } >> 521 /* >> 522 G4cout << "** G4VEnergyLossProcess::PreparePhysicsTable() for " >> 523 << GetProcessName() << " and " << particle->GetParticleName() >> 524 << " isMaster: " << isMaster << " isIonisation: " >> 525 << isIonisation << G4endl; >> 526 G4cout << " theDEDX: " << theDEDXTable >> 527 << " theRange: " << theRangeTableForLoss >> 528 << " theInverse: " << theInverseRangeTable >> 529 << " theLambda: " << theLambdaTable << G4endl; >> 530 */ 295 // forced biasing 531 // forced biasing 296 if(nullptr != biasManager) { << 532 if(biasManager) { 297 biasManager->Initialise(part,GetProcessNam 533 biasManager->Initialise(part,GetProcessName(),verboseLevel); 298 biasFlag = false; 534 biasFlag = false; 299 } 535 } 300 baseMat = bld->GetBaseMaterialFlag(); << 536 301 numberOfModels = modelManager->NumberOfModel << 537 // defined ID of secondary particles 302 currentModel = modelManager->GetModel(0); << 538 if(isMaster) { 303 G4EmTableUtil::UpdateModels(this, modelManag << 539 G4String nam1 = GetProcessName(); 304 numberOfModels, << 540 G4String nam4 = nam1 + "_split"; 305 mainSecondaries, << 541 G4String nam5 = nam1 + "_subcut"; 306 theParameters->U << 542 secID = G4PhysicsModelCatalog::Register(nam1); 307 theCuts = modelManager->Initialise(particle, << 543 biasID = G4PhysicsModelCatalog::Register(nam4); 308 verboseLe << 544 subsecID= G4PhysicsModelCatalog::Register(nam5); 309 // subcut processor << 545 } 310 if(isIonisation) { << 546 311 subcutProducer = lManager->SubCutProducer( << 547 // initialisation of models >> 548 G4int nmod = modelManager->NumberOfModels(); >> 549 for(G4int i=0; i<nmod; ++i) { >> 550 G4VEmModel* mod = modelManager->GetModel(i); >> 551 mod->SetMasterThread(isMaster); >> 552 mod->SetAngularGeneratorFlag( >> 553 theParameters->UseAngularGeneratorForIonisation()); >> 554 if(mod->HighEnergyLimit() > maxKinEnergy) { >> 555 mod->SetHighEnergyLimit(maxKinEnergy); >> 556 } 312 } 557 } 313 if(1 == nSCoffRegions) { << 558 theCuts = modelManager->Initialise(particle, secondaryParticle, 314 if((*scoffRegions)[0]->GetName() == "Defau << 559 theParameters->MinSubRange(), 315 delete scoffRegions; << 560 verboseLevel); 316 scoffRegions = nullptr; << 561 317 nSCoffRegions = 0; << 562 // Sub Cutoff >> 563 if(nSCoffRegions > 0) { >> 564 if(theParameters->MinSubRange() < 1.0) { useSubCutoff = true; } >> 565 >> 566 theSubCuts = modelManager->SubCutoff(); >> 567 >> 568 idxSCoffRegions = new G4bool[n]; >> 569 for (size_t j=0; j<n; ++j) { >> 570 >> 571 const G4MaterialCutsCouple* couple = >> 572 theCoupleTable->GetMaterialCutsCouple(j); >> 573 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); >> 574 >> 575 G4bool reg = false; >> 576 for(G4int i=0; i<nSCoffRegions; ++i) { >> 577 if( pcuts == scoffRegions[i]->GetProductionCuts()) { >> 578 reg = true; >> 579 break; >> 580 } >> 581 } >> 582 idxSCoffRegions[j] = reg; 318 } 583 } 319 } 584 } 320 585 321 if(1 < verboseLevel) { 586 if(1 < verboseLevel) { 322 G4cout << "G4VEnergyLossProcess::PrepearPh 587 G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done " 323 << " for " << GetProcessName() << " << 588 << " for local " << particle->GetParticleName() 324 << " isIon= " << isIon << " spline= << 589 << " isIon= " << isIon; 325 if(baseParticle) { 590 if(baseParticle) { 326 G4cout << "; base: " << baseParticle->Ge 591 G4cout << "; base: " << baseParticle->GetParticleName(); 327 } 592 } 328 G4cout << G4endl; << 329 G4cout << " chargeSqRatio= " << chargeSqRa 593 G4cout << " chargeSqRatio= " << chargeSqRatio 330 << " massRatio= " << massRatio 594 << " massRatio= " << massRatio 331 << " reduceFactor= " << reduceFacto 595 << " reduceFactor= " << reduceFactor << G4endl; 332 if (nSCoffRegions > 0) { << 596 if (nSCoffRegions) { 333 G4cout << " SubCut secondary production << 597 G4cout << " SubCutoff Regime is ON for regions: " << G4endl; 334 for (G4int i=0; i<nSCoffRegions; ++i) { 598 for (G4int i=0; i<nSCoffRegions; ++i) { 335 const G4Region* r = (*scoffRegions)[i] << 599 const G4Region* r = scoffRegions[i]; 336 G4cout << " " << r->GetName( 600 G4cout << " " << r->GetName() << G4endl; 337 } 601 } 338 } else if(nullptr != subcutProducer) { << 339 G4cout << " SubCut secondary production << 340 } 602 } 341 } 603 } 342 } 604 } 343 605 344 //....oooOO0OOooo........oooOO0OOooo........oo 606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 345 607 346 void G4VEnergyLossProcess::BuildPhysicsTable(c 608 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 347 { 609 { 348 if(1 < verboseLevel) { 610 if(1 < verboseLevel) { 349 G4cout << "### G4VEnergyLossProcess::Build 611 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for " 350 << GetProcessName() 612 << GetProcessName() 351 << " and particle " << part.GetPart 613 << " and particle " << part.GetParticleName() 352 << "; the first particle " << parti << 614 << "; local: " << particle->GetParticleName(); 353 if(baseParticle) { 615 if(baseParticle) { 354 G4cout << "; base: " << baseParticle->Ge 616 G4cout << "; base: " << baseParticle->GetParticleName(); 355 } 617 } 356 G4cout << G4endl; << 618 G4cout << " TablesAreBuilt= " << tablesAreBuilt 357 G4cout << " TablesAreBuilt= " << tables << 619 << " isIon= " << isIon << " " << this << G4endl; 358 << " spline=" << spline << " ptr: " << 359 } 620 } 360 621 361 if(&part == particle) { 622 if(&part == particle) { >> 623 >> 624 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 362 if(isMaster) { 625 if(isMaster) { >> 626 theDensityFactor = bld->GetDensityFactors(); >> 627 theDensityIdx = bld->GetCoupleIndexes(); 363 lManager->BuildPhysicsTable(particle, th 628 lManager->BuildPhysicsTable(particle, this); 364 629 365 } else { 630 } else { 366 const auto masterProcess = << 631 >> 632 const G4VEnergyLossProcess* masterProcess = 367 static_cast<const G4VEnergyLossProcess 633 static_cast<const G4VEnergyLossProcess*>(GetMasterProcess()); 368 634 369 numberOfModels = modelManager->NumberOfM << 635 // define density factors for worker thread 370 G4EmTableUtil::BuildLocalElossProcess(th << 636 bld->InitialiseBaseMaterials(masterProcess->DEDXTable()); 371 pa << 637 theDensityFactor = bld->GetDensityFactors(); >> 638 theDensityIdx = bld->GetCoupleIndexes(); >> 639 >> 640 // copy table pointers from master thread >> 641 SetDEDXTable(masterProcess->DEDXTable(),fRestricted); >> 642 SetDEDXTable(masterProcess->DEDXTableForSubsec(),fSubRestricted); >> 643 SetDEDXTable(masterProcess->DEDXunRestrictedTable(),fTotal); >> 644 SetDEDXTable(masterProcess->IonisationTable(),fIsIonisation); >> 645 SetDEDXTable(masterProcess->IonisationTableForSubsec(),fIsSubIonisation); >> 646 SetRangeTableForLoss(masterProcess->RangeTableForLoss()); >> 647 SetCSDARangeTable(masterProcess->CSDARangeTable()); >> 648 SetSecondaryRangeTable(masterProcess->SecondaryRangeTable()); >> 649 SetInverseRangeTable(masterProcess->InverseRangeTable()); >> 650 SetLambdaTable(masterProcess->LambdaTable()); >> 651 SetSubLambdaTable(masterProcess->SubLambdaTable()); >> 652 isIonisation = masterProcess->IsIonisationProcess(); >> 653 372 tablesAreBuilt = true; 654 tablesAreBuilt = true; 373 baseMat = masterProcess->UseBaseMaterial << 655 // local initialisation of models >> 656 G4bool printing = true; >> 657 G4int numberOfModels = modelManager->NumberOfModels(); >> 658 for(G4int i=0; i<numberOfModels; ++i) { >> 659 G4VEmModel* mod = GetModelByIndex(i, printing); >> 660 G4VEmModel* mod0= masterProcess->GetModelByIndex(i,printing); >> 661 mod->InitialiseLocal(particle, mod0); >> 662 } >> 663 374 lManager->LocalPhysicsTables(particle, t 664 lManager->LocalPhysicsTables(particle, this); 375 } 665 } 376 666 377 // needs to be done only once 667 // needs to be done only once 378 safetyHelper->InitialiseHelper(); 668 safetyHelper->InitialiseHelper(); 379 } 669 } 380 // Added tracking cut to avoid tracking arti << 381 // and identified deexcitation flag << 382 if(isIonisation) { << 383 atomDeexcitation = lManager->AtomDeexcitat << 384 if(nullptr != atomDeexcitation) { << 385 if(atomDeexcitation->IsPIXEActive()) { u << 386 } << 387 } << 388 << 389 // protection against double printout << 390 if(theParameters->IsPrintLocked()) { return; << 391 << 392 // explicitly defined printout by particle n 670 // explicitly defined printout by particle name 393 G4String num = part.GetParticleName(); 671 G4String num = part.GetParticleName(); 394 if(1 < verboseLevel || 672 if(1 < verboseLevel || 395 (0 < verboseLevel && (num == "e-" || 673 (0 < verboseLevel && (num == "e-" || 396 num == "e+" || n 674 num == "e+" || num == "mu+" || 397 num == "mu-" || n 675 num == "mu-" || num == "proton"|| 398 num == "pi+" || n 676 num == "pi+" || num == "pi-" || 399 num == "kaon+" || n 677 num == "kaon+" || num == "kaon-" || 400 num == "alpha" || n 678 num == "alpha" || num == "anti_proton" || 401 num == "GenericIon" << 679 num == "GenericIon"|| num == "alpha++" || 402 StreamInfo(G4cout, part); << 680 num == "alpha+" ))) >> 681 { >> 682 StreamInfo(G4cout, part); >> 683 } >> 684 >> 685 // Added tracking cut to avoid tracking artifacts >> 686 // identify deexcitation flag >> 687 if(isIonisation) { >> 688 atomDeexcitation = lManager->AtomDeexcitation(); >> 689 if(nSCoffRegions > 0) { subcutProducer = lManager->SubCutProducer(); } >> 690 if(atomDeexcitation) { >> 691 if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; } >> 692 } 403 } 693 } >> 694 /* >> 695 G4cout << "** G4VEnergyLossProcess::BuildPhysicsTable() for " >> 696 << GetProcessName() << " and " << particle->GetParticleName() >> 697 << " isMaster: " << isMaster << " isIonisation: " >> 698 << isIonisation << G4endl; >> 699 G4cout << " theDEDX: " << theDEDXTable >> 700 << " theRange: " << theRangeTableForLoss >> 701 << " theInverse: " << theInverseRangeTable >> 702 << " theLambda: " << theLambdaTable << G4endl; >> 703 */ >> 704 //if(1 < verboseLevel || verb) { 404 if(1 < verboseLevel) { 705 if(1 < verboseLevel) { 405 G4cout << "### G4VEnergyLossProcess::Build 706 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for " 406 << GetProcessName() 707 << GetProcessName() 407 << " and particle " << part.GetPart 708 << " and particle " << part.GetParticleName(); 408 if(isIonisation) { G4cout << " isIonisati << 709 if(isIonisation) { G4cout << " isIonisation flag = 1"; } 409 G4cout << " baseMat=" << baseMat << G4endl << 710 G4cout << G4endl; 410 } 711 } 411 } 712 } 412 713 413 //....oooOO0OOooo........oooOO0OOooo........oo 714 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 414 715 415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED 716 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType) 416 { 717 { >> 718 if(1 < verboseLevel ) { >> 719 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType >> 720 << " for " << GetProcessName() >> 721 << " and particle " << particle->GetParticleName() >> 722 << G4endl; >> 723 } 417 G4PhysicsTable* table = nullptr; 724 G4PhysicsTable* table = nullptr; 418 G4double emax = maxKinEnergy; 725 G4double emax = maxKinEnergy; 419 G4int bin = nBins; 726 G4int bin = nBins; 420 727 421 if(fTotal == tType) { 728 if(fTotal == tType) { 422 emax = maxKinEnergyCSDA; 729 emax = maxKinEnergyCSDA; 423 bin = nBinsCSDA; 730 bin = nBinsCSDA; 424 table = theDEDXunRestrictedTable; 731 table = theDEDXunRestrictedTable; 425 } else if(fRestricted == tType) { 732 } else if(fRestricted == tType) { 426 table = theDEDXTable; 733 table = theDEDXTable; >> 734 } else if(fSubRestricted == tType) { >> 735 table = theDEDXSubTable; 427 } else { 736 } else { 428 G4cout << "G4VEnergyLossProcess::BuildDEDX 737 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type " 429 << tType << G4endl; 738 << tType << G4endl; 430 } 739 } >> 740 >> 741 // Access to materials >> 742 const G4ProductionCutsTable* theCoupleTable= >> 743 G4ProductionCutsTable::GetProductionCutsTable(); >> 744 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 745 431 if(1 < verboseLevel) { 746 if(1 < verboseLevel) { 432 G4cout << "G4VEnergyLossProcess::BuildDEDX << 747 G4cout << numOfCouples << " materials" 433 << " for " << GetProcessName() << 748 << " minKinEnergy= " << minKinEnergy 434 << " and " << particle->GetParticle << 749 << " maxKinEnergy= " << emax 435 << "spline=" << spline << G4endl; << 750 << " nbin= " << bin >> 751 << " EmTableType= " << tType >> 752 << " table= " << table << " " << this >> 753 << G4endl; 436 } 754 } 437 if(nullptr == table) { return table; } << 755 if(!table) { return table; } 438 756 439 G4LossTableBuilder* bld = lManager->GetTable 757 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 440 G4EmTableUtil::BuildDEDXTable(this, particle << 758 G4bool splineFlag = theParameters->Spline(); 441 table, minKinE << 759 G4PhysicsLogVector* aVector = nullptr; 442 verboseLevel, << 760 G4PhysicsLogVector* bVector = nullptr; >> 761 >> 762 for(size_t i=0; i<numOfCouples; ++i) { >> 763 >> 764 if(1 < verboseLevel) { >> 765 G4cout << "G4VEnergyLossProcess::BuildDEDXVector Idx= " << i >> 766 << " flagTable= " << table->GetFlag(i) >> 767 << " Flag= " << bld->GetFlag(i) << G4endl; >> 768 } >> 769 if(bld->GetFlag(i)) { >> 770 >> 771 // create physics vector and fill it >> 772 const G4MaterialCutsCouple* couple = >> 773 theCoupleTable->GetMaterialCutsCouple(i); >> 774 if((*table)[i]) { delete (*table)[i]; } >> 775 if(bVector) { >> 776 aVector = new G4PhysicsLogVector(*bVector); >> 777 } else { >> 778 bVector = new G4PhysicsLogVector(minKinEnergy, emax, bin); >> 779 aVector = bVector; >> 780 } >> 781 aVector->SetSpline(splineFlag); >> 782 >> 783 modelManager->FillDEDXVector(aVector, couple, tType); >> 784 if(splineFlag) { aVector->FillSecondDerivatives(); } >> 785 >> 786 // Insert vector for this material into the table >> 787 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); >> 788 } >> 789 } >> 790 >> 791 if(1 < verboseLevel) { >> 792 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for " >> 793 << particle->GetParticleName() >> 794 << " and process " << GetProcessName() >> 795 << G4endl; >> 796 if(2 < verboseLevel) G4cout << (*table) << G4endl; >> 797 } >> 798 443 return table; 799 return table; 444 } 800 } 445 801 446 //....oooOO0OOooo........oooOO0OOooo........oo 802 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 447 803 448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam << 804 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType) 449 { 805 { 450 if(nullptr == theLambdaTable) { return theLa << 806 G4PhysicsTable* table = nullptr; >> 807 >> 808 if(fRestricted == tType) { >> 809 table = theLambdaTable; >> 810 } else if(fSubRestricted == tType) { >> 811 table = theSubLambdaTable; >> 812 } else { >> 813 G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type " >> 814 << tType << G4endl; >> 815 } >> 816 >> 817 if(1 < verboseLevel) { >> 818 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type " >> 819 << tType << " for process " >> 820 << GetProcessName() << " and particle " >> 821 << particle->GetParticleName() >> 822 << " EmTableType= " << tType >> 823 << " table= " << table >> 824 << G4endl; >> 825 } >> 826 if(!table) {return table;} >> 827 >> 828 // Access to materials >> 829 const G4ProductionCutsTable* theCoupleTable= >> 830 G4ProductionCutsTable::GetProductionCutsTable(); >> 831 size_t numOfCouples = theCoupleTable->GetTableSize(); 451 832 452 G4double scale = theParameters->MaxKinEnergy << 453 G4int nbin = << 454 theParameters->NumberOfBinsPerDecade()*G4l << 455 scale = nbin/G4Log(scale); << 456 << 457 G4LossTableBuilder* bld = lManager->GetTable 833 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 458 G4EmTableUtil::BuildLambdaTable(this, partic << 834 theDensityFactor = bld->GetDensityFactors(); 459 bld, theLamb << 835 theDensityIdx = bld->GetCoupleIndexes(); 460 minKinEnergy << 836 461 verboseLevel << 837 G4bool splineFlag = theParameters->Spline(); 462 return theLambdaTable; << 838 G4PhysicsLogVector* aVector = nullptr; >> 839 G4double scale = G4Log(maxKinEnergy/minKinEnergy); >> 840 >> 841 for(size_t i=0; i<numOfCouples; ++i) { >> 842 >> 843 if (bld->GetFlag(i)) { >> 844 >> 845 // create physics vector and fill it >> 846 const G4MaterialCutsCouple* couple = >> 847 theCoupleTable->GetMaterialCutsCouple(i); >> 848 delete (*table)[i]; >> 849 >> 850 G4bool startNull = true; >> 851 G4double emin = >> 852 MinPrimaryEnergy(particle,couple->GetMaterial(),(*theCuts)[i]); >> 853 if(minKinEnergy > emin) { >> 854 emin = minKinEnergy; >> 855 startNull = false; >> 856 } >> 857 >> 858 G4double emax = maxKinEnergy; >> 859 if(emax <= emin) { emax = 2*emin; } >> 860 G4int bin = G4lrint(nBins*G4Log(emax/emin)/scale); >> 861 bin = std::max(bin, 3); >> 862 aVector = new G4PhysicsLogVector(emin, emax, bin); >> 863 aVector->SetSpline(splineFlag); >> 864 >> 865 modelManager->FillLambdaVector(aVector, couple, startNull, tType); >> 866 if(splineFlag) { aVector->FillSecondDerivatives(); } >> 867 >> 868 // Insert vector for this material into the table >> 869 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); >> 870 } >> 871 } >> 872 >> 873 if(1 < verboseLevel) { >> 874 G4cout << "Lambda table is built for " >> 875 << particle->GetParticleName() >> 876 << G4endl; >> 877 } >> 878 >> 879 return table; 463 } 880 } 464 881 465 //....oooOO0OOooo........oooOO0OOooo........oo 882 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 466 883 467 void G4VEnergyLossProcess::StreamInfo(std::ost 884 void G4VEnergyLossProcess::StreamInfo(std::ostream& out, 468 const G4ParticleDefinition& pa 885 const G4ParticleDefinition& part, G4bool rst) const 469 { 886 { 470 G4String indent = (rst ? " " : ""); 887 G4String indent = (rst ? " " : ""); 471 out << std::setprecision(6); 888 out << std::setprecision(6); 472 out << G4endl << indent << GetProcessName() 889 out << G4endl << indent << GetProcessName() << ": "; 473 if (!rst) out << " for " << part.GetParticle << 890 if (!rst) out << " for " << part.GetParticleName(); 474 out << " XStype:" << fXSType << 891 out << " SubType= " << GetProcessSubType() << G4endl 475 << " SubType=" << GetProcessSubType() < << 476 << " dE/dx and range tables from " 892 << " dE/dx and range tables from " 477 << G4BestUnit(minKinEnergy,"Energy") 893 << G4BestUnit(minKinEnergy,"Energy") 478 << " to " << G4BestUnit(maxKinEnergy,"En 894 << " to " << G4BestUnit(maxKinEnergy,"Energy") 479 << " in " << nBins << " bins" << G4endl 895 << " in " << nBins << " bins" << G4endl 480 << " Lambda tables from threshold t 896 << " Lambda tables from threshold to " 481 << G4BestUnit(maxKinEnergy,"Energy") 897 << G4BestUnit(maxKinEnergy,"Energy") 482 << ", " << theParameters->NumberOfBinsPe 898 << ", " << theParameters->NumberOfBinsPerDecade() 483 << " bins/decade, spline: " << spline << 899 << " bins per decade, spline: " >> 900 << theParameters->Spline() 484 << G4endl; 901 << G4endl; 485 if(nullptr != theRangeTableForLoss && isIoni << 902 if(theRangeTableForLoss && isIonisation) { 486 out << " StepFunction=(" << dRoverRan << 903 out << " finalRange(mm)= " << finalRange/mm 487 << finalRange/mm << " mm)" << 904 << ", dRoverRange= " << dRoverRange 488 << ", integ: " << fXSType << 905 << ", integral: " << integral 489 << ", fluct: " << lossFluctuationFlag << 906 << ", fluct: " << lossFluctuationFlag 490 << ", linLossLim= " << linLossLimit << 907 << ", linLossLimit= " << linLossLimit 491 << G4endl; << 908 << G4endl; 492 } 909 } 493 StreamProcessInfo(out); 910 StreamProcessInfo(out); 494 modelManager->DumpModelList(out, verboseLeve 911 modelManager->DumpModelList(out, verboseLevel); 495 if(nullptr != theCSDARangeTable && isIonisat << 912 if(theCSDARangeTable && isIonisation) { 496 out << " CSDA range table up" 913 out << " CSDA range table up" 497 << " to " << G4BestUnit(maxKinEnergyCS << 914 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy") 498 << " in " << nBinsCSDA << " bins" << G << 915 << " in " << nBinsCSDA << " bins" << G4endl; 499 } 916 } 500 if(nSCoffRegions>0 && isIonisation) { 917 if(nSCoffRegions>0 && isIonisation) { 501 out << " Subcutoff sampling in " << n 918 out << " Subcutoff sampling in " << nSCoffRegions 502 << " regions" << G4endl; << 919 << " regions" << G4endl; 503 } 920 } 504 if(2 < verboseLevel) { 921 if(2 < verboseLevel) { 505 for(std::size_t i=0; i<7; ++i) { << 922 out << " DEDXTable address= " << theDEDXTable << G4endl; 506 auto ta = theData->Table(i); << 923 if(theDEDXTable && isIonisation) out << (*theDEDXTable) << G4endl; 507 out << " " << tnames[i] << " addres << 924 out << "non restricted DEDXTable address= " 508 if(nullptr != ta) { out << *ta << G4endl << 925 << theDEDXunRestrictedTable << G4endl; >> 926 if(theDEDXunRestrictedTable && isIonisation) { >> 927 out << (*theDEDXunRestrictedTable) << G4endl; >> 928 } >> 929 if(theDEDXSubTable && isIonisation) { >> 930 out << (*theDEDXSubTable) << G4endl; >> 931 } >> 932 out << " CSDARangeTable address= " << theCSDARangeTable << G4endl; >> 933 if(theCSDARangeTable && isIonisation) { >> 934 out << (*theCSDARangeTable) << G4endl; >> 935 } >> 936 out << " RangeTableForLoss address= " << theRangeTableForLoss << G4endl; >> 937 if(theRangeTableForLoss && isIonisation) { >> 938 out << (*theRangeTableForLoss) << G4endl; >> 939 } >> 940 out << " InverseRangeTable address= " << theInverseRangeTable << G4endl; >> 941 if(theInverseRangeTable && isIonisation) { >> 942 out << (*theInverseRangeTable) << G4endl; >> 943 } >> 944 out << " LambdaTable address= " << theLambdaTable << G4endl; >> 945 if(theLambdaTable && isIonisation) { >> 946 out << (*theLambdaTable) << G4endl; >> 947 } >> 948 out << " SubLambdaTable address= " << theSubLambdaTable << G4endl; >> 949 if(theSubLambdaTable && isIonisation) { >> 950 out << (*theSubLambdaTable) << G4endl; 509 } 951 } 510 } 952 } 511 } 953 } 512 954 513 //....oooOO0OOooo........oooOO0OOooo........oo 955 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 514 956 515 void G4VEnergyLossProcess::ActivateSubCutoff(c << 957 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r) 516 { 958 { 517 if(nullptr == scoffRegions) { << 959 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 518 scoffRegions = new std::vector<const G4Reg << 960 const G4Region* reg = r; >> 961 if (!reg) { >> 962 reg = regionStore->GetRegion("DefaultRegionForTheWorld", false); 519 } 963 } >> 964 520 // the region is in the list 965 // the region is in the list 521 if(!scoffRegions->empty()) { << 966 if (nSCoffRegions > 0) { 522 for (auto & reg : *scoffRegions) { << 967 for (G4int i=0; i<nSCoffRegions; ++i) { 523 if (reg == r) { return; } << 968 if (reg == scoffRegions[i]) { >> 969 return; >> 970 } 524 } 971 } 525 } 972 } 526 // new region 973 // new region 527 scoffRegions->push_back(r); << 974 if(val) { 528 ++nSCoffRegions; << 975 scoffRegions.push_back(reg); 529 } << 976 ++nSCoffRegions; 530 << 531 //....oooOO0OOooo........oooOO0OOooo........oo << 532 << 533 G4bool G4VEnergyLossProcess::IsRegionForCubcut << 534 { << 535 if(0 == nSCoffRegions) { return true; } << 536 const G4Region* r = aTrack.GetVolume()->GetL << 537 for(auto & reg : *scoffRegions) { << 538 if(r == reg) { return true; } << 539 } 977 } 540 return false; << 541 } 978 } 542 979 543 //....oooOO0OOooo........oooOO0OOooo........oo 980 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 544 981 545 void G4VEnergyLossProcess::StartTracking(G4Tra 982 void G4VEnergyLossProcess::StartTracking(G4Track* track) 546 { 983 { >> 984 /* >> 985 G4cout << track->GetDefinition()->GetParticleName() >> 986 << " e(MeV)= " << track->GetKineticEnergy() >> 987 << " baseParticle " << baseParticle << " proc " << this; >> 988 if(particle) G4cout << " " << particle->GetParticleName(); >> 989 G4cout << " isIon= " << isIon << " dedx " << theDEDXTable <<G4endl; >> 990 */ 547 // reset parameters for the new track 991 // reset parameters for the new track 548 theNumberOfInteractionLengthLeft = -1.0; 992 theNumberOfInteractionLengthLeft = -1.0; 549 mfpKinEnergy = DBL_MAX; << 993 currentInteractionLength = mfpKinEnergy = DBL_MAX; 550 preStepLambda = 0.0; << 994 preStepRangeEnergy = 0.0; 551 currentCouple = nullptr; << 552 995 553 // reset ion 996 // reset ion 554 if(isIon) { 997 if(isIon) { 555 const G4double newmass = track->GetDefinit << 998 chargeSqRatio = 0.5; 556 massRatio = (nullptr == baseParticle) ? CL << 999 557 : baseParticle->GetPDGMass()/newmass; << 1000 G4double newmass = track->GetDefinition()->GetPDGMass(); 558 logMassRatio = G4Log(massRatio); << 1001 if(baseParticle) { >> 1002 massRatio = baseParticle->GetPDGMass()/newmass; >> 1003 } else if(theGenericIon) { >> 1004 massRatio = proton_mass_c2/newmass; >> 1005 } else { >> 1006 massRatio = 1.0; >> 1007 } 559 } 1008 } 560 // forced biasing only for primary particles 1009 // forced biasing only for primary particles 561 if(nullptr != biasManager) { << 1010 if(biasManager) { 562 if(0 == track->GetParentID()) { 1011 if(0 == track->GetParentID()) { >> 1012 // primary particle 563 biasFlag = true; 1013 biasFlag = true; 564 biasManager->ResetForcedInteraction(); 1014 biasManager->ResetForcedInteraction(); 565 } 1015 } 566 } 1016 } 567 } 1017 } 568 1018 569 //....oooOO0OOooo........oooOO0OOooo........oo 1019 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 570 1020 571 G4double G4VEnergyLossProcess::AlongStepGetPhy 1021 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength( 572 const G4Track& tr << 1022 const G4Track&,G4double,G4double,G4double&, 573 G4GPILSelection* 1023 G4GPILSelection* selection) 574 { 1024 { 575 G4double x = DBL_MAX; 1025 G4double x = DBL_MAX; 576 *selection = aGPILSelection; 1026 *selection = aGPILSelection; 577 if(isIonisation && currentModel->IsActive(pr 1027 if(isIonisation && currentModel->IsActive(preStepScaledEnergy)) { 578 GetScaledRangeForScaledEnergy(preStepScale << 1028 fRange = GetScaledRangeForScaledEnergy(preStepScaledEnergy)*reduceFactor; 579 x = (useCutAsFinalRange) ? std::min(finalR << 1029 G4double finR = (rndmStepFlag) ? std::min(finalRange, 580 currentCouple->GetProductionCuts()->GetP 1030 currentCouple->GetProductionCuts()->GetProductionCut(1)) : finalRange; 581 x = (fRange > x) ? fRange*dRoverRange + x* << 1031 x = (fRange > finR) ? 582 : fRange; << 1032 fRange*dRoverRange + finR*(1.0 - dRoverRange)*(2.0 - finR/fRange) : fRange; 583 /* << 1033 // if(particle->GetPDGMass() > 0.9*GeV) 584 G4cout<<"AlongStepGPIL: " << GetProcessN << 1034 /* 585 << " fRange=" << fRange << " finR=" << finR << 1035 G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy >> 1036 <<" range= "<<fRange << " idx= " << basedCoupleIndex >> 1037 << " finR= " << finR >> 1038 << " limit= " << x <<G4endl; >> 1039 G4cout << "massRatio= " << massRatio << " Q^2= " << chargeSqRatio >> 1040 << " finR= " << finR << " dRoverRange= " << dRoverRange >> 1041 << " finalRange= " << finalRange << G4endl; 586 */ 1042 */ 587 } 1043 } >> 1044 //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy >> 1045 //<<" stepLimit= "<<x<<G4endl; 588 return x; 1046 return x; 589 } 1047 } 590 1048 591 //....oooOO0OOooo........oooOO0OOooo........oo 1049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 592 1050 593 G4double G4VEnergyLossProcess::PostStepGetPhys 1051 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength( 594 const G4Track& tr 1052 const G4Track& track, 595 G4double previo 1053 G4double previousStepSize, 596 G4ForceCondition* 1054 G4ForceCondition* condition) 597 { 1055 { 598 // condition is set to "Not Forced" 1056 // condition is set to "Not Forced" 599 *condition = NotForced; 1057 *condition = NotForced; 600 G4double x = DBL_MAX; 1058 G4double x = DBL_MAX; 601 1059 602 // initialisation of material, mass, charge, 1060 // initialisation of material, mass, charge, model 603 // at the beginning of the step 1061 // at the beginning of the step 604 DefineMaterial(track.GetMaterialCutsCouple() 1062 DefineMaterial(track.GetMaterialCutsCouple()); 605 preStepKinEnergy = track.GetKineticEne << 1063 preStepKinEnergy = track.GetKineticEnergy(); 606 preStepScaledEnergy = preStepKinEnergy*ma << 1064 preStepScaledEnergy = preStepKinEnergy*massRatio; 607 SelectModel(preStepScaledEnergy); 1065 SelectModel(preStepScaledEnergy); 608 1066 609 if(!currentModel->IsActive(preStepScaledEner 1067 if(!currentModel->IsActive(preStepScaledEnergy)) { 610 theNumberOfInteractionLengthLeft = -1.0; 1068 theNumberOfInteractionLengthLeft = -1.0; 611 mfpKinEnergy = DBL_MAX; << 612 preStepLambda = 0.0; << 613 currentInteractionLength = DBL_MAX; 1069 currentInteractionLength = DBL_MAX; 614 return x; << 1070 return x; 615 } 1071 } 616 1072 617 // change effective charge of a charged part << 1073 // change effective charge of an ion on fly 618 if(isIon) { 1074 if(isIon) { 619 const G4double q2 = currentModel->ChargeSq << 1075 G4double q2 = currentModel->ChargeSquareRatio(track); 620 fFactor = q2*biasFactor; << 1076 if(q2 != chargeSqRatio && q2 > 0.0) { 621 if(baseMat) { fFactor *= (*theDensityFacto << 1077 chargeSqRatio = q2; 622 reduceFactor = 1.0/(fFactor*massRatio); << 1078 fFactor = q2*biasFactor*(*theDensityFactor)[currentCoupleIndex]; 623 if (lossFluctuationFlag) { << 1079 reduceFactor = 1.0/(fFactor*massRatio); 624 auto fluc = currentModel->GetModelOfFluc << 625 fluc->SetParticleAndCharge(track.GetDefi << 626 } 1080 } 627 } 1081 } >> 1082 // if(particle->GetPDGMass() > 0.9*GeV) >> 1083 //G4cout << "q2= "<<chargeSqRatio << " massRatio= " << massRatio << G4endl; 628 1084 629 // forced biasing only for primary particles 1085 // forced biasing only for primary particles 630 if(biasManager) { 1086 if(biasManager) { 631 if(0 == track.GetParentID() && biasFlag && 1087 if(0 == track.GetParentID() && biasFlag && 632 biasManager->ForcedInteractionRegion((G << 1088 biasManager->ForcedInteractionRegion(currentCoupleIndex)) { 633 return biasManager->GetStepLimit((G4int) << 1089 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize); 634 } 1090 } 635 } 1091 } 636 1092 637 ComputeLambdaForScaledEnergy(preStepScaledEn << 1093 // compute mean free path 638 << 1094 if(preStepScaledEnergy < mfpKinEnergy) { 639 // zero cross section << 1095 if (integral) { ComputeLambdaForScaledEnergy(preStepScaledEnergy); } 640 if(preStepLambda <= 0.0) { << 1096 else { preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy); } 641 theNumberOfInteractionLengthLeft = -1.0; << 1097 642 currentInteractionLength = DBL_MAX; << 1098 // zero cross section 643 } else { << 1099 if(preStepLambda <= 0.0) { >> 1100 theNumberOfInteractionLengthLeft = -1.0; >> 1101 currentInteractionLength = DBL_MAX; >> 1102 } >> 1103 } 644 1104 645 // non-zero cross section << 1105 // non-zero cross section >> 1106 if(preStepLambda > 0.0) { 646 if (theNumberOfInteractionLengthLeft < 0.0 1107 if (theNumberOfInteractionLengthLeft < 0.0) { 647 1108 648 // beggining of tracking (or just after 1109 // beggining of tracking (or just after DoIt of this process) 649 theNumberOfInteractionLengthLeft = -G4Lo << 1110 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() ); 650 theInitialNumberOfInteractionLength = th 1111 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 651 1112 652 } else if(currentInteractionLength < DBL_M 1113 } else if(currentInteractionLength < DBL_MAX) { 653 1114 654 // subtract NumberOfInteractionLengthLef 1115 // subtract NumberOfInteractionLengthLeft using previous step 655 theNumberOfInteractionLengthLeft -= 1116 theNumberOfInteractionLengthLeft -= 656 previousStepSize/currentInteractionLen 1117 previousStepSize/currentInteractionLength; 657 1118 658 theNumberOfInteractionLengthLeft = 1119 theNumberOfInteractionLengthLeft = 659 std::max(theNumberOfInteractionLengthL 1120 std::max(theNumberOfInteractionLengthLeft, 0.0); 660 } 1121 } 661 1122 662 // new mean free path and step limit 1123 // new mean free path and step limit 663 currentInteractionLength = 1.0/preStepLamb 1124 currentInteractionLength = 1.0/preStepLambda; 664 x = theNumberOfInteractionLengthLeft * cur 1125 x = theNumberOfInteractionLengthLeft * currentInteractionLength; 665 } 1126 } 666 #ifdef G4VERBOSE 1127 #ifdef G4VERBOSE 667 if (verboseLevel>2) { << 1128 if (verboseLevel>2){ >> 1129 // if(particle->GetPDGMass() > 0.9*GeV){ 668 G4cout << "G4VEnergyLossProcess::PostStepG 1130 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength "; 669 G4cout << "[ " << GetProcessName() << "]" 1131 G4cout << "[ " << GetProcessName() << "]" << G4endl; 670 G4cout << " for " << track.GetDefinition() 1132 G4cout << " for " << track.GetDefinition()->GetParticleName() 671 << " in Material " << currentMate 1133 << " in Material " << currentMaterial->GetName() 672 << " Ekin(MeV)= " << preStepKinEner 1134 << " Ekin(MeV)= " << preStepKinEnergy/MeV 673 << " track material: " << track.Get << 1135 << " " << track.GetMaterial()->GetName() 674 <<G4endl; 1136 <<G4endl; 675 G4cout << "MeanFreePath = " << currentInte 1137 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 676 << "InteractionLength= " << x/cm << 1138 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl; 677 } 1139 } 678 #endif 1140 #endif 679 return x; 1141 return x; 680 } 1142 } 681 1143 682 //....oooOO0OOooo........oooOO0OOooo........oo 1144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 683 1145 684 void << 1146 void G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e) 685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 686 { 1147 { 687 // cross section increased with energy << 1148 // condition to skip recomputation of cross section 688 if(fXSType == fEmIncreasing) { << 1149 G4double epeak = theEnergyOfCrossSectionMax[currentCoupleIndex]; 689 if(e*invLambdaFactor < mfpKinEnergy) { << 1150 if(e <= epeak && e/lambdaFactor >= mfpKinEnergy) { return; } 690 preStepLambda = GetLambdaForScaledEnergy << 1151 691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 1152 // recomputation is needed 692 } << 1153 if (e <= epeak) { 693 << 1154 preStepLambda = GetLambdaForScaledEnergy(e); 694 // cross section has one peak << 1155 mfpKinEnergy = e; 695 } else if(fXSType == fEmOnePeak) { << 1156 696 const G4double epeak = (*theEnergyOfCrossS << 1157 } else { 697 if(e <= epeak) { << 1158 G4double e1 = e*lambdaFactor; 698 if(e*invLambdaFactor < mfpKinEnergy) { << 1159 if(e1 > epeak) { 699 preStepLambda = GetLambdaForScaledEner << 1160 preStepLambda = GetLambdaForScaledEnergy(e); 700 mfpKinEnergy = (preStepLambda > 0.0) ? << 1161 mfpKinEnergy = e; 701 } << 1162 G4double preStepLambda1 = GetLambdaForScaledEnergy(e1); 702 } else if(e < mfpKinEnergy) { << 1163 if(preStepLambda1 > preStepLambda) { 703 const G4double e1 = std::max(epeak, e*la << 704 mfpKinEnergy = e1; << 705 preStepLambda = GetLambdaForScaledEnergy << 706 } << 707 << 708 // cross section has more than one peaks << 709 } else if(fXSType == fEmTwoPeaks) { << 710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple << 711 const G4double e1peak = xs->e1peak; << 712 << 713 // below the 1st peak << 714 if(e <= e1peak) { << 715 if(e*invLambdaFactor < mfpKinEnergy) { << 716 preStepLambda = GetLambdaForScaledEner << 717 mfpKinEnergy = (preStepLambda > 0.0) ? << 718 } << 719 return; << 720 } << 721 const G4double e1deep = xs->e1deep; << 722 // above the 1st peak, below the deep << 723 if(e <= e1deep) { << 724 if(mfpKinEnergy >= e1deep || e <= mfpKin << 725 const G4double e1 = std::max(e1peak, e << 726 mfpKinEnergy = e1; << 727 preStepLambda = GetLambdaForScaledEner << 728 } << 729 return; << 730 } << 731 const G4double e2peak = xs->e2peak; << 732 // above the deep, below 2nd peak << 733 if(e <= e2peak) { << 734 if(e*invLambdaFactor < mfpKinEnergy) { << 735 mfpKinEnergy = e; << 736 preStepLambda = GetLambdaForScaledEner << 737 } << 738 return; << 739 } << 740 const G4double e2deep = xs->e2deep; << 741 // above the 2nd peak, below the deep << 742 if(e <= e2deep) { << 743 if(mfpKinEnergy >= e2deep || e <= mfpKin << 744 const G4double e1 = std::max(e2peak, e << 745 mfpKinEnergy = e1; 1164 mfpKinEnergy = e1; 746 preStepLambda = GetLambdaForScaledEner << 1165 preStepLambda = preStepLambda1; 747 } << 748 return; << 749 } << 750 const G4double e3peak = xs->e3peak; << 751 // above the deep, below 3d peak << 752 if(e <= e3peak) { << 753 if(e*invLambdaFactor < mfpKinEnergy) { << 754 mfpKinEnergy = e; << 755 preStepLambda = GetLambdaForScaledEner << 756 } 1166 } 757 return; << 1167 } else { 758 } << 1168 preStepLambda = fFactor*theCrossSectionMax[currentCoupleIndex]; 759 // above 3d peak << 1169 mfpKinEnergy = epeak; 760 if(e <= mfpKinEnergy) { << 761 const G4double e1 = std::max(e3peak, e*l << 762 mfpKinEnergy = e1; << 763 preStepLambda = GetLambdaForScaledEnergy << 764 } 1170 } 765 // integral method is not used << 766 } else { << 767 preStepLambda = GetLambdaForScaledEnergy(e << 768 } 1171 } 769 } 1172 } 770 1173 771 //....oooOO0OOooo........oooOO0OOooo........oo 1174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 772 1175 773 G4VParticleChange* G4VEnergyLossProcess::Along 1176 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track, 774 1177 const G4Step& step) 775 { 1178 { 776 fParticleChange.InitializeForAlongStep(track 1179 fParticleChange.InitializeForAlongStep(track); 777 // The process has range table - calculate e 1180 // The process has range table - calculate energy loss 778 if(!isIonisation || !currentModel->IsActive( 1181 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) { 779 return &fParticleChange; 1182 return &fParticleChange; 780 } 1183 } 781 1184 >> 1185 // Get the actual (true) Step length 782 G4double length = step.GetStepLength(); 1186 G4double length = step.GetStepLength(); >> 1187 if(length <= 0.0) { return &fParticleChange; } 783 G4double eloss = 0.0; 1188 G4double eloss = 0.0; 784 1189 785 /* << 1190 /* 786 if(-1 < verboseLevel) { 1191 if(-1 < verboseLevel) { 787 const G4ParticleDefinition* d = track.GetP 1192 const G4ParticleDefinition* d = track.GetParticleDefinition(); 788 G4cout << "AlongStepDoIt for " 1193 G4cout << "AlongStepDoIt for " 789 << GetProcessName() << " and partic << 1194 << GetProcessName() << " and particle " 790 << " eScaled(MeV)=" << preStepScal << 1195 << d->GetParticleName() 791 << " range(mm)=" << fRange/mm << " << 1196 << " eScaled(MeV)= " << preStepScaledEnergy/MeV 792 << " rf=" << reduceFactor << " q^ << 1197 << " range(mm)= " << fRange/mm 793 << " md=" << d->GetPDGMass() << " << 1198 << " s(mm)= " << length/mm 794 << " " << track.GetMaterial()->Get << 1199 << " rf= " << reduceFactor >> 1200 << " q^2= " << chargeSqRatio >> 1201 << " md= " << d->GetPDGMass() >> 1202 << " status= " << track.GetTrackStatus() >> 1203 << " " << track.GetMaterial()->GetName() >> 1204 << G4endl; 795 } 1205 } 796 */ 1206 */ >> 1207 797 const G4DynamicParticle* dynParticle = track 1208 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); 798 1209 799 // define new weight for primary and seconda 1210 // define new weight for primary and secondaries 800 G4double weight = fParticleChange.GetParentW 1211 G4double weight = fParticleChange.GetParentWeight(); 801 if(weightFlag) { 1212 if(weightFlag) { 802 weight /= biasFactor; 1213 weight /= biasFactor; 803 fParticleChange.ProposeWeight(weight); 1214 fParticleChange.ProposeWeight(weight); 804 } 1215 } 805 1216 806 // stopping, check actual range and kinetic << 1217 // stopping 807 if (length >= fRange || preStepKinEnergy <= 1218 if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) { 808 eloss = preStepKinEnergy; 1219 eloss = preStepKinEnergy; 809 if (useDeexcitation) { 1220 if (useDeexcitation) { 810 atomDeexcitation->AlongStepDeexcitation( 1221 atomDeexcitation->AlongStepDeexcitation(scTracks, step, 811 << 1222 eloss, currentCoupleIndex); 812 if(scTracks.size() > 0) { FillSecondarie << 1223 if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); } 813 eloss = std::max(eloss, 0.0); 1224 eloss = std::max(eloss, 0.0); 814 } 1225 } 815 fParticleChange.SetProposedKineticEnergy(0 1226 fParticleChange.SetProposedKineticEnergy(0.0); 816 fParticleChange.ProposeLocalEnergyDeposit( 1227 fParticleChange.ProposeLocalEnergyDeposit(eloss); 817 return &fParticleChange; 1228 return &fParticleChange; 818 } 1229 } 819 // zero step length with non-zero range << 1230 //G4cout << theDEDXTable << " idx= " << basedCoupleIndex 820 if(length <= 0.0) { return &fParticleChange; << 1231 // << " " << GetProcessName() << " "<< currentMaterial->GetName()<<G4endl; 821 << 1232 //if(particle->GetParticleName() == "e-")G4cout << (*theDEDXTable) <<G4endl; 822 // Short step 1233 // Short step 823 eloss = length*GetDEDXForScaledEnergy(preSte << 1234 eloss = GetDEDXForScaledEnergy(preStepScaledEnergy)*length; 824 LogSca << 1235 825 /* << 1236 //G4cout << "eloss= " << eloss << G4endl; 826 G4cout << "##### Short STEP: eloss= " << elo << 1237 827 << " Escaled=" << preStepScaledEnergy << 828 << " R=" << fRange << 829 << " L=" << length << 830 << " fFactor=" << fFactor << " minE=" << mi << 831 << " idxBase=" << basedCoupleIndex << G4end << 832 */ << 833 // Long step 1238 // Long step 834 if(eloss > preStepKinEnergy*linLossLimit) { 1239 if(eloss > preStepKinEnergy*linLossLimit) { 835 1240 836 const G4double x = (fRange - length)/reduc << 1241 G4double x = (fRange - length)/reduceFactor; 837 const G4double de = preStepKinEnergy - Sca << 1242 //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl; 838 if(de > 0.0) { eloss = de; } << 1243 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio; >> 1244 839 /* 1245 /* 840 if(-1 < verboseLevel) 1246 if(-1 < verboseLevel) 841 G4cout << " Long STEP: rPre(mm)=" << 1247 G4cout << "Long STEP: rPre(mm)= " 842 << GetScaledRangeForScaledEnergy( 1248 << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm 843 << " x(mm)=" << x/mm << 1249 << " rPost(mm)= " << x/mm 844 << " eloss(MeV)=" << eloss/MeV << 1250 << " ePre(MeV)= " << preStepScaledEnergy/MeV 845 << " rFactor=" << reduceFactor << 1251 << " eloss(MeV)= " << eloss/MeV 846 << " massRatio=" << massRatio << 1252 << " eloss0(MeV)= " >> 1253 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV >> 1254 << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV 847 << G4endl; 1255 << G4endl; 848 */ 1256 */ 849 } 1257 } 850 1258 851 /* << 1259 /* >> 1260 G4double eloss0 = eloss; 852 if(-1 < verboseLevel ) { 1261 if(-1 < verboseLevel ) { 853 G4cout << "Before fluct: eloss(MeV)= " << 1262 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV 854 << " e-eloss= " << preStepKinEnergy 1263 << " e-eloss= " << preStepKinEnergy-eloss 855 << " step(mm)= " << length/mm << " << 1264 << " step(mm)= " << length/mm 856 << " fluct= " << lossFluctuationFla << 1265 << " range(mm)= " << fRange/mm >> 1266 << " fluct= " << lossFluctuationFlag >> 1267 << G4endl; 857 } 1268 } 858 */ 1269 */ 859 1270 860 const G4double cut = (*theCuts)[currentCoupl << 1271 G4double cut = (*theCuts)[currentCoupleIndex]; 861 G4double esec = 0.0; 1272 G4double esec = 0.0; 862 1273 >> 1274 //G4cout << "cut= " << cut << " useSubCut= " << useSubCutoff << G4endl; >> 1275 >> 1276 // SubCutOff >> 1277 if(useSubCutoff && !subcutProducer) { >> 1278 if(idxSCoffRegions[currentCoupleIndex]) { >> 1279 >> 1280 G4bool yes = false; >> 1281 const G4StepPoint* prePoint = step.GetPreStepPoint(); >> 1282 >> 1283 // Check boundary >> 1284 if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; } >> 1285 >> 1286 // Check PrePoint >> 1287 else { >> 1288 G4double preSafety = prePoint->GetSafety(); >> 1289 G4double rcut = >> 1290 currentCouple->GetProductionCuts()->GetProductionCut(1); >> 1291 >> 1292 // recompute presafety >> 1293 if(preSafety < rcut) { >> 1294 preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition(), >> 1295 rcut); >> 1296 } >> 1297 >> 1298 if(preSafety < rcut) { yes = true; } >> 1299 >> 1300 // Check PostPoint >> 1301 else { >> 1302 G4double postSafety = preSafety - length; >> 1303 if(postSafety < rcut) { >> 1304 postSafety = safetyHelper->ComputeSafety( >> 1305 step.GetPostStepPoint()->GetPosition(), rcut); >> 1306 if(postSafety < rcut) { yes = true; } >> 1307 } >> 1308 } >> 1309 } >> 1310 >> 1311 // Decided to start subcut sampling >> 1312 if(yes) { >> 1313 >> 1314 cut = (*theSubCuts)[currentCoupleIndex]; >> 1315 eloss -= GetSubDEDXForScaledEnergy(preStepScaledEnergy)*length; >> 1316 esec = SampleSubCutSecondaries(scTracks, step, >> 1317 currentModel,currentCoupleIndex); >> 1318 // add bremsstrahlung sampling >> 1319 /* >> 1320 if(nProcesses > 0) { >> 1321 for(G4int i=0; i<nProcesses; ++i) { >> 1322 (scProcesses[i])->SampleSubCutSecondaries( >> 1323 scTracks, step, (scProcesses[i])-> >> 1324 SelectModelForMaterial(preStepKinEnergy, currentCoupleIndex), >> 1325 currentCoupleIndex); >> 1326 } >> 1327 } >> 1328 */ >> 1329 } >> 1330 } >> 1331 } >> 1332 863 // Corrections, which cannot be tabulated 1333 // Corrections, which cannot be tabulated 864 if(isIon) { 1334 if(isIon) { >> 1335 G4double eadd = 0.0; >> 1336 G4double eloss_before = eloss; 865 currentModel->CorrectionsAlongStep(current 1337 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 866 length, << 1338 eloss, eadd, length); 867 eloss = std::max(eloss, 0.0); << 1339 if(eloss < 0.0) { eloss = 0.5*eloss_before; } 868 } 1340 } 869 1341 870 // Sample fluctuations if not full energy lo << 1342 // Sample fluctuations 871 if(eloss >= preStepKinEnergy) { << 1343 if (lossFluctuationFlag) { 872 eloss = preStepKinEnergy; << 873 << 874 } else if (lossFluctuationFlag) { << 875 const G4double tmax = currentModel->MaxSec << 876 const G4double tcut = std::min(cut, tmax); << 877 G4VEmFluctuationModel* fluc = currentModel 1344 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations(); 878 eloss = fluc->SampleFluctuations(currentCo << 1345 if(eloss + esec < preStepKinEnergy) { 879 tcut, tma << 1346 880 /* << 1347 G4double tmax = 881 if(-1 < verboseLevel) << 1348 std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut); >> 1349 eloss = fluc->SampleFluctuations(currentCouple,dynParticle, >> 1350 tmax,length,eloss); >> 1351 /* >> 1352 if(-1 < verboseLevel) 882 G4cout << "After fluct: eloss(MeV)= " << 1353 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV 883 << " fluc= " << (eloss-eloss0)/Me 1354 << " fluc= " << (eloss-eloss0)/MeV 884 << " ChargeSqRatio= " << chargeSq 1355 << " ChargeSqRatio= " << chargeSqRatio 885 << " massRatio= " << massRatio << << 1356 << " massRatio= " << massRatio 886 */ << 1357 << " tmax= " << tmax >> 1358 << G4endl; >> 1359 */ >> 1360 } 887 } 1361 } 888 1362 889 // deexcitation 1363 // deexcitation 890 if (useDeexcitation) { 1364 if (useDeexcitation) { 891 G4double esecfluo = preStepKinEnergy; << 1365 G4double esecfluo = preStepKinEnergy - esec; 892 G4double de = esecfluo; 1366 G4double de = esecfluo; >> 1367 //G4double eloss0 = eloss; >> 1368 /* >> 1369 G4cout << "### 1: E(keV)= " << preStepKinEnergy/keV >> 1370 << " Efluomax(keV)= " << de/keV >> 1371 << " Eloss(keV)= " << eloss/keV << G4endl; >> 1372 */ 893 atomDeexcitation->AlongStepDeexcitation(sc 1373 atomDeexcitation->AlongStepDeexcitation(scTracks, step, 894 de << 1374 de, currentCoupleIndex); 895 1375 896 // sum of de-excitation energies 1376 // sum of de-excitation energies 897 esecfluo -= de; 1377 esecfluo -= de; 898 1378 899 // subtracted from energy loss 1379 // subtracted from energy loss 900 if(eloss >= esecfluo) { 1380 if(eloss >= esecfluo) { 901 esec += esecfluo; 1381 esec += esecfluo; 902 eloss -= esecfluo; 1382 eloss -= esecfluo; 903 } else { 1383 } else { 904 esec += esecfluo; 1384 esec += esecfluo; 905 eloss = 0.0; 1385 eloss = 0.0; 906 } 1386 } >> 1387 /* >> 1388 if(esecfluo > 0.0) { >> 1389 G4cout << "### 2: E(keV)= " << preStepKinEnergy/keV >> 1390 << " Esec(keV)= " << esec/keV >> 1391 << " Esecf(kV)= " << esecfluo/keV >> 1392 << " Eloss0(kV)= " << eloss0/keV >> 1393 << " Eloss(keV)= " << eloss/keV >> 1394 << G4endl; >> 1395 } >> 1396 */ 907 } 1397 } 908 if(nullptr != subcutProducer && IsRegionForC << 1398 if(subcutProducer && idxSCoffRegions[currentCoupleIndex]) { 909 subcutProducer->SampleSecondaries(step, sc 1399 subcutProducer->SampleSecondaries(step, scTracks, eloss, cut); 910 } 1400 } 911 // secondaries from atomic de-excitation and << 1401 if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); } 912 if(!scTracks.empty()) { FillSecondariesAlong << 913 1402 914 // Energy balance 1403 // Energy balance 915 G4double finalT = preStepKinEnergy - eloss - 1404 G4double finalT = preStepKinEnergy - eloss - esec; 916 if (finalT <= lowestKinEnergy) { 1405 if (finalT <= lowestKinEnergy) { 917 eloss += finalT; 1406 eloss += finalT; 918 finalT = 0.0; 1407 finalT = 0.0; 919 } else if(isIon) { 1408 } else if(isIon) { 920 fParticleChange.SetProposedCharge( 1409 fParticleChange.SetProposedCharge( 921 currentModel->GetParticleCharge(track.Ge 1410 currentModel->GetParticleCharge(track.GetParticleDefinition(), 922 currentM 1411 currentMaterial,finalT)); 923 } 1412 } >> 1413 924 eloss = std::max(eloss, 0.0); 1414 eloss = std::max(eloss, 0.0); 925 1415 926 fParticleChange.SetProposedKineticEnergy(fin 1416 fParticleChange.SetProposedKineticEnergy(finalT); 927 fParticleChange.ProposeLocalEnergyDeposit(el 1417 fParticleChange.ProposeLocalEnergyDeposit(eloss); 928 /* 1418 /* 929 if(-1 < verboseLevel) { 1419 if(-1 < verboseLevel) { 930 G4double del = finalT + eloss + esec - pre 1420 G4double del = finalT + eloss + esec - preStepKinEnergy; 931 G4cout << "Final value eloss(MeV)= " << el 1421 G4cout << "Final value eloss(MeV)= " << eloss/MeV 932 << " preStepKinEnergy= " << preStep 1422 << " preStepKinEnergy= " << preStepKinEnergy 933 << " postStepKinEnergy= " << finalT 1423 << " postStepKinEnergy= " << finalT 934 << " de(keV)= " << del/keV 1424 << " de(keV)= " << del/keV 935 << " lossFlag= " << lossFluctuation 1425 << " lossFlag= " << lossFluctuationFlag 936 << " status= " << track.GetTrackSt 1426 << " status= " << track.GetTrackStatus() 937 << G4endl; 1427 << G4endl; 938 } 1428 } 939 */ 1429 */ 940 return &fParticleChange; 1430 return &fParticleChange; 941 } 1431 } 942 1432 943 //....oooOO0OOooo........oooOO0OOooo........oo 1433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 944 1434 945 void G4VEnergyLossProcess::FillSecondariesAlon << 1435 void >> 1436 G4VEnergyLossProcess::FillSecondariesAlongStep(G4double&, G4double& weight) 946 { 1437 { 947 const std::size_t n0 = scTracks.size(); << 1438 G4int n0 = scTracks.size(); 948 G4double weight = wt; << 1439 949 // weight may be changed by biasing manager 1440 // weight may be changed by biasing manager 950 if(biasManager) { 1441 if(biasManager) { 951 if(biasManager->SecondaryBiasingRegion((G4 << 1442 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) { 952 weight *= 1443 weight *= 953 biasManager->ApplySecondaryBiasing(scT << 1444 biasManager->ApplySecondaryBiasing(scTracks, currentCoupleIndex); 954 } 1445 } 955 } 1446 } 956 1447 957 // fill secondaries 1448 // fill secondaries 958 const std::size_t n = scTracks.size(); << 1449 G4int n = scTracks.size(); 959 fParticleChange.SetNumberOfSecondaries((G4in << 1450 fParticleChange.SetNumberOfSecondaries(n); 960 1451 961 for(std::size_t i=0; i<n; ++i) { << 1452 for(G4int i=0; i<n; ++i) { 962 G4Track* t = scTracks[i]; 1453 G4Track* t = scTracks[i]; 963 if(nullptr != t) { << 1454 if(t) { 964 t->SetWeight(weight); 1455 t->SetWeight(weight); 965 pParticleChange->AddSecondary(t); 1456 pParticleChange->AddSecondary(t); 966 G4int pdg = t->GetDefinition()->GetPDGEn << 1457 if(i >= n0) { t->SetCreatorModelIndex(biasID); } 967 if (i < n0) { << 1458 //G4cout << "Secondary(along step) has weight " << t->GetWeight() 968 if (pdg == 22) { << 1459 //<< ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl; 969 t->SetCreatorModelID(gpixeID); << 970 } else if (pdg == 11) { << 971 t->SetCreatorModelID(epixeID); << 972 } else { << 973 t->SetCreatorModelID(biasID); << 974 } << 975 } else { << 976 t->SetCreatorModelID(biasID); << 977 } << 978 } 1460 } 979 } 1461 } 980 scTracks.clear(); 1462 scTracks.clear(); 981 } 1463 } 982 1464 983 //....oooOO0OOooo........oooOO0OOooo........oo 1465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 984 1466 >> 1467 G4double >> 1468 G4VEnergyLossProcess::SampleSubCutSecondaries(std::vector<G4Track*>& tracks, >> 1469 const G4Step& step, >> 1470 G4VEmModel* model, >> 1471 G4int idx) >> 1472 { >> 1473 // Fast check weather subcutoff can work >> 1474 G4double esec = 0.0; >> 1475 G4double subcut = (*theSubCuts)[idx]; >> 1476 G4double cut = (*theCuts)[idx]; >> 1477 if(cut <= subcut) { return esec; } >> 1478 >> 1479 const G4Track* track = step.GetTrack(); >> 1480 const G4DynamicParticle* dp = track->GetDynamicParticle(); >> 1481 G4double e = dp->GetKineticEnergy()*massRatio; >> 1482 G4double cross = (*theDensityFactor)[idx]*chargeSqRatio >> 1483 *(((*theSubLambdaTable)[(*theDensityIdx)[idx]])->Value(e, idxSubLambda)); >> 1484 G4double length = step.GetStepLength(); >> 1485 >> 1486 // negligible probability to get any interaction >> 1487 if(length*cross < perMillion) { return esec; } >> 1488 /* >> 1489 if(-1 < verboseLevel) >> 1490 G4cout << "<<< Subcutoff for " << GetProcessName() >> 1491 << " cross(1/mm)= " << cross*mm << ">>>" >> 1492 << " e(MeV)= " << preStepScaledEnergy >> 1493 << " matIdx= " << currentCoupleIndex >> 1494 << G4endl; >> 1495 */ >> 1496 >> 1497 // Sample subcutoff secondaries >> 1498 G4StepPoint* preStepPoint = step.GetPreStepPoint(); >> 1499 G4StepPoint* postStepPoint = step.GetPostStepPoint(); >> 1500 G4ThreeVector prepoint = preStepPoint->GetPosition(); >> 1501 G4ThreeVector dr = postStepPoint->GetPosition() - prepoint; >> 1502 G4double pretime = preStepPoint->GetGlobalTime(); >> 1503 G4double dt = postStepPoint->GetGlobalTime() - pretime; >> 1504 G4double fragment = 0.0; >> 1505 >> 1506 do { >> 1507 G4double del = -G4Log(G4UniformRand())/cross; >> 1508 fragment += del/length; >> 1509 if (fragment > 1.0) { break; } >> 1510 >> 1511 // sample secondaries >> 1512 secParticles.clear(); >> 1513 model->SampleSecondaries(&secParticles,track->GetMaterialCutsCouple(), >> 1514 dp,subcut,cut); >> 1515 >> 1516 // position of subcutoff particles >> 1517 G4ThreeVector r = prepoint + fragment*dr; >> 1518 std::vector<G4DynamicParticle*>::iterator it; >> 1519 for(it=secParticles.begin(); it!=secParticles.end(); ++it) { >> 1520 >> 1521 G4Track* t = new G4Track((*it), pretime + fragment*dt, r); >> 1522 t->SetTouchableHandle(track->GetTouchableHandle()); >> 1523 t->SetCreatorModelIndex(subsecID); >> 1524 tracks.push_back(t); >> 1525 esec += t->GetKineticEnergy(); >> 1526 if (t->GetParticleDefinition() == thePositron) { >> 1527 esec += 2.0*electron_mass_c2; >> 1528 } >> 1529 >> 1530 /* >> 1531 if(-1 < verboseLevel) >> 1532 G4cout << "New track " >> 1533 << t->GetParticleDefinition()->GetParticleName() >> 1534 << " e(keV)= " << t->GetKineticEnergy()/keV >> 1535 << " fragment= " << fragment >> 1536 << G4endl; >> 1537 */ >> 1538 } >> 1539 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko >> 1540 } while (fragment <= 1.0); >> 1541 return esec; >> 1542 } >> 1543 >> 1544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1545 985 G4VParticleChange* G4VEnergyLossProcess::PostS 1546 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track, 986 1547 const G4Step& step) 987 { 1548 { 988 // clear number of interaction lengths in an << 1549 // In all cases clear number of interaction lengths 989 theNumberOfInteractionLengthLeft = -1.0; 1550 theNumberOfInteractionLengthLeft = -1.0; 990 mfpKinEnergy = DBL_MAX; << 1551 mfpKinEnergy = currentInteractionLength = DBL_MAX; 991 1552 992 fParticleChange.InitializeForPostStep(track) 1553 fParticleChange.InitializeForPostStep(track); 993 const G4double finalT = track.GetKineticEner << 1554 G4double finalT = track.GetKineticEnergy(); >> 1555 if(finalT <= lowestKinEnergy) { return &fParticleChange; } 994 1556 995 const G4double postStepScaledEnergy = finalT << 1557 G4double postStepScaledEnergy = finalT*massRatio; 996 SelectModel(postStepScaledEnergy); 1558 SelectModel(postStepScaledEnergy); 997 1559 998 if(!currentModel->IsActive(postStepScaledEne 1560 if(!currentModel->IsActive(postStepScaledEnergy)) { 999 return &fParticleChange; 1561 return &fParticleChange; 1000 } 1562 } 1001 /* 1563 /* 1002 if(1 < verboseLevel) { << 1564 if(-1 < verboseLevel) { 1003 G4cout<<GetProcessName()<<" PostStepDoIt: << 1565 G4cout << GetProcessName() >> 1566 << "::PostStepDoIt: E(MeV)= " << finalT/MeV >> 1567 << G4endl; 1004 } 1568 } 1005 */ 1569 */ >> 1570 1006 // forced process - should happen only once 1571 // forced process - should happen only once per track 1007 if(biasFlag) { 1572 if(biasFlag) { 1008 if(biasManager->ForcedInteractionRegion(( << 1573 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) { 1009 biasFlag = false; 1574 biasFlag = false; 1010 } 1575 } 1011 } 1576 } 1012 const G4DynamicParticle* dp = track.GetDyna << 1013 1577 1014 // Integral approach 1578 // Integral approach 1015 if (fXSType != fEmNoIntegral) { << 1579 if (integral) { 1016 const G4double logFinalT = dp->GetLogKine << 1580 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy); 1017 G4double lx = GetLambdaForScaledEnergy(po << 1581 /* 1018 lo << 1582 if(preStepLambda<lx && 1 < verboseLevel) { 1019 lx = std::max(lx, 0.0); << 1583 G4cout << "WARNING: for " << particle->GetParticleName() 1020 << 1584 << " and " << GetProcessName() 1021 // if both lg and lx are zero then no int << 1585 << " E(MeV)= " << finalT/MeV 1022 if(preStepLambda*G4UniformRand() >= lx) { << 1586 << " preLambda= " << preStepLambda >> 1587 << " < " << lx << " (postLambda) " >> 1588 << G4endl; >> 1589 } >> 1590 */ >> 1591 if(lx <= 0.0 || preStepLambda*G4UniformRand() > lx) { 1023 return &fParticleChange; 1592 return &fParticleChange; 1024 } 1593 } 1025 } 1594 } 1026 1595 >> 1596 SelectModel(postStepScaledEnergy); >> 1597 1027 // define new weight for primary and second 1598 // define new weight for primary and secondaries 1028 G4double weight = fParticleChange.GetParent 1599 G4double weight = fParticleChange.GetParentWeight(); 1029 if(weightFlag) { 1600 if(weightFlag) { 1030 weight /= biasFactor; 1601 weight /= biasFactor; 1031 fParticleChange.ProposeWeight(weight); 1602 fParticleChange.ProposeWeight(weight); 1032 } 1603 } 1033 1604 1034 const G4double tcut = (*theCuts)[currentCou << 1605 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); >> 1606 G4double tcut = (*theCuts)[currentCoupleIndex]; 1035 1607 1036 // sample secondaries 1608 // sample secondaries 1037 secParticles.clear(); 1609 secParticles.clear(); 1038 currentModel->SampleSecondaries(&secParticl << 1610 //G4cout<< "@@@ Eprimary= "<<dynParticle->GetKineticEnergy()/MeV >> 1611 // << " cut= " << tcut/MeV << G4endl; >> 1612 currentModel->SampleSecondaries(&secParticles, currentCouple, >> 1613 dynParticle, tcut); 1039 1614 1040 const G4int num0 = (G4int)secParticles.size << 1615 G4int num0 = secParticles.size(); 1041 1616 1042 // bremsstrahlung splitting or Russian roul 1617 // bremsstrahlung splitting or Russian roulette 1043 if(biasManager) { 1618 if(biasManager) { 1044 if(biasManager->SecondaryBiasingRegion((G << 1619 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) { 1045 G4double eloss = 0.0; 1620 G4double eloss = 0.0; 1046 weight *= biasManager->ApplySecondaryBi 1621 weight *= biasManager->ApplySecondaryBiasing( 1047 secPart 1622 secParticles, 1048 track, 1623 track, currentModel, 1049 &fParti 1624 &fParticleChange, eloss, 1050 (G4int) << 1625 currentCoupleIndex, tcut, 1051 step.Ge 1626 step.GetPostStepPoint()->GetSafety()); 1052 if(eloss > 0.0) { 1627 if(eloss > 0.0) { 1053 eloss += fParticleChange.GetLocalEner 1628 eloss += fParticleChange.GetLocalEnergyDeposit(); 1054 fParticleChange.ProposeLocalEnergyDep 1629 fParticleChange.ProposeLocalEnergyDeposit(eloss); 1055 } 1630 } 1056 } 1631 } 1057 } 1632 } 1058 1633 1059 // save secondaries 1634 // save secondaries 1060 const G4int num = (G4int)secParticles.size( << 1635 G4int num = secParticles.size(); 1061 if(num > 0) { 1636 if(num > 0) { 1062 1637 1063 fParticleChange.SetNumberOfSecondaries(nu 1638 fParticleChange.SetNumberOfSecondaries(num); 1064 G4double time = track.GetGlobalTime(); 1639 G4double time = track.GetGlobalTime(); 1065 1640 1066 G4int n1(0), n2(0); << 1067 if(num0 > mainSecondaries) { << 1068 currentModel->FillNumberOfSecondaries(n << 1069 } << 1070 << 1071 for (G4int i=0; i<num; ++i) { 1641 for (G4int i=0; i<num; ++i) { 1072 if(nullptr != secParticles[i]) { << 1642 if(secParticles[i]) { 1073 G4Track* t = new G4Track(secParticles 1643 G4Track* t = new G4Track(secParticles[i], time, track.GetPosition()); 1074 t->SetTouchableHandle(track.GetToucha 1644 t->SetTouchableHandle(track.GetTouchableHandle()); 1075 if (biasManager) { 1645 if (biasManager) { 1076 t->SetWeight(weight * biasManager-> << 1646 t->SetWeight(biasManager->GetWeight(i)); 1077 } else { 1647 } else { 1078 t->SetWeight(weight); 1648 t->SetWeight(weight); 1079 } 1649 } 1080 if(i < num0) { << 1650 if(i < num0) { t->SetCreatorModelIndex(secID); } 1081 t->SetCreatorModelID(secID); << 1651 else { t->SetCreatorModelIndex(biasID); } 1082 } else if(i < num0 + n1) { << 1083 t->SetCreatorModelID(tripletID); << 1084 } else { << 1085 t->SetCreatorModelID(biasID); << 1086 } << 1087 1652 1088 //G4cout << "Secondary(post step) has 1653 //G4cout << "Secondary(post step) has weight " << t->GetWeight() 1089 // << ", kenergy " << t->GetKin 1654 // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" 1090 // << " time= " << time/ns << " 1655 // << " time= " << time/ns << " ns " << G4endl; 1091 pParticleChange->AddSecondary(t); 1656 pParticleChange->AddSecondary(t); 1092 } 1657 } 1093 } 1658 } 1094 } 1659 } 1095 1660 1096 if(0.0 == fParticleChange.GetProposedKineti 1661 if(0.0 == fParticleChange.GetProposedKineticEnergy() && 1097 fAlive == fParticleChange.GetTrackStatus 1662 fAlive == fParticleChange.GetTrackStatus()) { 1098 if(particle->GetProcessManager()->GetAtRe 1663 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0) 1099 { fParticleChange.ProposeTrackStatus 1664 { fParticleChange.ProposeTrackStatus(fStopButAlive); } 1100 else { fParticleChange.ProposeTrackStatus 1665 else { fParticleChange.ProposeTrackStatus(fStopAndKill); } 1101 } 1666 } 1102 1667 1103 /* 1668 /* 1104 if(-1 < verboseLevel) { 1669 if(-1 < verboseLevel) { 1105 G4cout << "::PostStepDoIt: Sample seconda 1670 G4cout << "::PostStepDoIt: Sample secondary; Efin= " 1106 << fParticleChange.GetProposedKineticEner 1671 << fParticleChange.GetProposedKineticEnergy()/MeV 1107 << " MeV; model= (" << currentMode 1672 << " MeV; model= (" << currentModel->LowEnergyLimit() 1108 << ", " << currentModel->HighEner 1673 << ", " << currentModel->HighEnergyLimit() << ")" 1109 << " preStepLambda= " << preStepL 1674 << " preStepLambda= " << preStepLambda 1110 << " dir= " << track.GetMomentumD 1675 << " dir= " << track.GetMomentumDirection() 1111 << " status= " << track.GetTrackS 1676 << " status= " << track.GetTrackStatus() 1112 << G4endl; 1677 << G4endl; 1113 } 1678 } 1114 */ 1679 */ 1115 return &fParticleChange; 1680 return &fParticleChange; 1116 } 1681 } 1117 1682 1118 //....oooOO0OOooo........oooOO0OOooo........o 1683 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1119 1684 1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl 1685 G4bool G4VEnergyLossProcess::StorePhysicsTable( 1121 const G4ParticleDefinition* part, cons << 1686 const G4ParticleDefinition* part, const G4String& directory, >> 1687 G4bool ascii) 1122 { 1688 { 1123 if (!isMaster || nullptr != baseParticle || << 1689 G4bool res = true; 1124 for(std::size_t i=0; i<7; ++i) { << 1690 //G4cout << "G4VEnergyLossProcess::StorePhysicsTable: " << part->GetParticleName() 1125 // ionisation table only for ionisation p << 1691 // << " " << directory << " " << ascii << G4endl; 1126 if (nullptr == theData->Table(i) || (!isI << 1692 if (!isMaster || baseParticle || part != particle ) return res; 1127 continue; << 1693 1128 } << 1694 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX")) 1129 if (-1 < verboseLevel) { << 1695 {res = false;} 1130 G4cout << "G4VEnergyLossProcess::StoreP << 1696 1131 << " " << particle->GetParticleName() << 1697 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr")) 1132 << " " << GetProcessName() << 1698 {res = false;} 1133 << " " << tnames[i] << " " << theDat << 1699 1134 } << 1700 if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX")) 1135 if (!G4EmTableUtil::StoreTable(this, part << 1701 {res = false;} 1136 dir, tnames[i], verboseLevel, asci << 1702 1137 return false; << 1703 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation")) >> 1704 {res = false;} >> 1705 >> 1706 if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation")) >> 1707 {res = false;} >> 1708 >> 1709 if(isIonisation && >> 1710 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange")) >> 1711 {res = false;} >> 1712 >> 1713 if(isIonisation && >> 1714 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range")) >> 1715 {res = false;} >> 1716 >> 1717 if(isIonisation && >> 1718 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange")) >> 1719 {res = false;} >> 1720 >> 1721 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda")) >> 1722 {res = false;} >> 1723 >> 1724 if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda")) >> 1725 {res = false;} >> 1726 >> 1727 if ( !res ) { >> 1728 if(1 < verboseLevel) { >> 1729 G4cout << "Physics tables are stored for " >> 1730 << particle->GetParticleName() >> 1731 << " and process " << GetProcessName() >> 1732 << " in the directory <" << directory >> 1733 << "> " << G4endl; 1138 } 1734 } >> 1735 } else { >> 1736 G4cout << "Fail to store Physics Tables for " >> 1737 << particle->GetParticleName() >> 1738 << " and process " << GetProcessName() >> 1739 << " in the directory <" << directory >> 1740 << "> " << G4endl; 1139 } 1741 } 1140 return true; << 1742 return res; 1141 } 1743 } 1142 1744 1143 //....oooOO0OOooo........oooOO0OOooo........o 1745 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1144 1746 1145 G4bool 1747 G4bool 1146 G4VEnergyLossProcess::RetrievePhysicsTable(co 1748 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 1147 co << 1749 const G4String& directory, >> 1750 G4bool ascii) 1148 { 1751 { 1149 if (!isMaster || nullptr != baseParticle || << 1752 G4bool res = true; 1150 for(std::size_t i=0; i<7; ++i) { << 1753 if (!isMaster) return res; 1151 // ionisation table only for ionisation p << 1754 const G4String particleName = part->GetParticleName(); 1152 if (!isIonisation && 1 == i) { continue; << 1755 1153 if(!G4EmTableUtil::RetrieveTable(this, pa << 1756 if(1 < verboseLevel) { 1154 verboseL << 1757 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for " 1155 return false; << 1758 << particleName << " and process " << GetProcessName() >> 1759 << "; tables_are_built= " << tablesAreBuilt >> 1760 << G4endl; >> 1761 } >> 1762 if(particle == part) { >> 1763 >> 1764 if ( !baseParticle ) { >> 1765 >> 1766 G4bool fpi = true; >> 1767 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi)) >> 1768 {fpi = false;} >> 1769 >> 1770 // ionisation table keeps individual dEdx and not sum of sub-processes >> 1771 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false)) >> 1772 {fpi = false;} >> 1773 >> 1774 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi)) >> 1775 {res = false;} >> 1776 >> 1777 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory, >> 1778 "DEDXnr",false)) >> 1779 {res = false;} >> 1780 >> 1781 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory, >> 1782 "CSDARange",false)) >> 1783 {res = false;} >> 1784 >> 1785 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory, >> 1786 "InverseRange",fpi)) >> 1787 {res = false;} >> 1788 >> 1789 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true)) >> 1790 {res = false;} >> 1791 >> 1792 G4bool yes = false; >> 1793 if(nSCoffRegions > 0) {yes = true;} >> 1794 >> 1795 if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes)) >> 1796 {res = false;} >> 1797 >> 1798 if(!RetrieveTable(part,theSubLambdaTable,ascii,directory, >> 1799 "SubLambda",yes)) >> 1800 {res = false;} >> 1801 >> 1802 if(!fpi) yes = false; >> 1803 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory, >> 1804 "SubIonisation",yes)) >> 1805 {res = false;} 1156 } 1806 } 1157 } 1807 } >> 1808 >> 1809 return res; >> 1810 } >> 1811 >> 1812 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... >> 1813 >> 1814 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part, >> 1815 G4PhysicsTable* aTable, G4bool ascii, >> 1816 const G4String& directory, >> 1817 const G4String& tname) >> 1818 { >> 1819 //G4cout << "G4VEnergyLossProcess::StoreTable: " << aTable >> 1820 // << " " << directory << " " << tname << G4endl; >> 1821 G4bool res = true; >> 1822 if ( aTable ) { >> 1823 const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii); >> 1824 G4cout << name << G4endl; >> 1825 //G4cout << *aTable << G4endl; >> 1826 if( !aTable->StorePhysicsTable(name,ascii)) res = false; >> 1827 } >> 1828 return res; >> 1829 } >> 1830 >> 1831 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... >> 1832 >> 1833 G4bool >> 1834 G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, >> 1835 G4PhysicsTable* aTable, >> 1836 G4bool ascii, >> 1837 const G4String& directory, >> 1838 const G4String& tname, >> 1839 G4bool mandatory) >> 1840 { >> 1841 G4bool isRetrieved = false; >> 1842 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii); >> 1843 if(aTable) { >> 1844 if(aTable->ExistPhysicsTable(filename)) { >> 1845 if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) { >> 1846 isRetrieved = true; >> 1847 if(theParameters->Spline()) { >> 1848 size_t n = aTable->length(); >> 1849 for(size_t i=0; i<n; ++i) { >> 1850 if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); } >> 1851 } >> 1852 } >> 1853 if (0 < verboseLevel) { >> 1854 G4cout << tname << " table for " << part->GetParticleName() >> 1855 << " is Retrieved from <" << filename << ">" >> 1856 << G4endl; >> 1857 } >> 1858 } >> 1859 } >> 1860 } >> 1861 if(mandatory && !isRetrieved) { >> 1862 if(0 < verboseLevel) { >> 1863 G4cout << tname << " table for " << part->GetParticleName() >> 1864 << " from file <" >> 1865 << filename << "> is not Retrieved" >> 1866 << G4endl; >> 1867 } >> 1868 return false; >> 1869 } 1158 return true; 1870 return true; 1159 } 1871 } 1160 1872 1161 //....oooOO0OOooo........oooOO0OOooo........o 1873 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1162 1874 1163 G4double G4VEnergyLossProcess::GetDEDXDispers 1875 G4double G4VEnergyLossProcess::GetDEDXDispersion( 1164 const G4Mat 1876 const G4MaterialCutsCouple *couple, 1165 const G4Dyn 1877 const G4DynamicParticle* dp, 1166 G4dou 1878 G4double length) 1167 { 1879 { 1168 DefineMaterial(couple); 1880 DefineMaterial(couple); 1169 G4double ekin = dp->GetKineticEnergy(); 1881 G4double ekin = dp->GetKineticEnergy(); 1170 SelectModel(ekin*massRatio); 1882 SelectModel(ekin*massRatio); 1171 G4double tmax = currentModel->MaxSecondaryK 1883 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp); 1172 G4double tcut = std::min(tmax,(*theCuts)[cu << 1884 tmax = std::min(tmax,(*theCuts)[currentCoupleIndex]); 1173 G4double d = 0.0; 1885 G4double d = 0.0; 1174 G4VEmFluctuationModel* fm = currentModel->G 1886 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations(); 1175 if(nullptr != fm) { d = fm->Dispersion(curr << 1887 if(fm) { d = fm->Dispersion(currentMaterial,dp,tmax,length); } 1176 return d; 1888 return d; 1177 } 1889 } 1178 1890 1179 //....oooOO0OOooo........oooOO0OOooo........o 1891 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1180 1892 1181 G4double << 1893 G4double G4VEnergyLossProcess::CrossSectionPerVolume( 1182 G4VEnergyLossProcess::CrossSectionPerVolume(G << 1894 G4double kineticEnergy, const G4MaterialCutsCouple* couple) 1183 c << 1184 G << 1185 { 1895 { 1186 // Cross section per volume is calculated 1896 // Cross section per volume is calculated 1187 DefineMaterial(couple); 1897 DefineMaterial(couple); 1188 G4double cross = 0.0; 1898 G4double cross = 0.0; 1189 if (nullptr != theLambdaTable) { << 1899 if(theLambdaTable) { 1190 cross = GetLambdaForScaledEnergy(kineticE << 1900 cross = GetLambdaForScaledEnergy(kineticEnergy*massRatio); 1191 logKinet << 1192 } else { 1901 } else { 1193 SelectModel(kineticEnergy*massRatio); 1902 SelectModel(kineticEnergy*massRatio); 1194 cross = (!baseMat) ? biasFactor : biasFac << 1903 cross = biasFactor*(*theDensityFactor)[currentCoupleIndex] 1195 cross *= (currentModel->CrossSectionPerVo << 1904 *(currentModel->CrossSectionPerVolume(currentMaterial, 1196 << 1905 particle, kineticEnergy, >> 1906 (*theCuts)[currentCoupleIndex])); 1197 } 1907 } 1198 return std::max(cross, 0.0); 1908 return std::max(cross, 0.0); 1199 } 1909 } 1200 1910 1201 //....oooOO0OOooo........oooOO0OOooo........o 1911 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1202 1912 1203 G4double G4VEnergyLossProcess::MeanFreePath(c 1913 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track) 1204 { 1914 { 1205 DefineMaterial(track.GetMaterialCutsCouple( 1915 DefineMaterial(track.GetMaterialCutsCouple()); 1206 const G4double kinEnergy = track.GetKine << 1916 G4double cs = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio); 1207 const G4double logKinEnergy = track.GetDyna << 1208 const G4double cs = GetLambdaForScaledEnerg << 1209 << 1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1917 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1211 } 1918 } 1212 1919 1213 //....oooOO0OOooo........oooOO0OOooo........o 1920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1214 1921 1215 G4double G4VEnergyLossProcess::ContinuousStep 1922 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, 1216 1923 G4double x, G4double y, 1217 1924 G4double& z) 1218 { 1925 { 1219 return AlongStepGetPhysicalInteractionLengt << 1926 G4GPILSelection sel; >> 1927 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel); 1220 } 1928 } 1221 1929 1222 //....oooOO0OOooo........oooOO0OOooo........o 1930 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1223 1931 1224 G4double G4VEnergyLossProcess::GetMeanFreePat 1932 G4double G4VEnergyLossProcess::GetMeanFreePath( 1225 const G4Track& t 1933 const G4Track& track, 1226 G4double, 1934 G4double, 1227 G4ForceCondition 1935 G4ForceCondition* condition) 1228 1936 1229 { 1937 { 1230 *condition = NotForced; 1938 *condition = NotForced; 1231 return MeanFreePath(track); 1939 return MeanFreePath(track); 1232 } 1940 } 1233 1941 1234 //....oooOO0OOooo........oooOO0OOooo........o 1942 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1235 1943 1236 G4double G4VEnergyLossProcess::GetContinuousS 1944 G4double G4VEnergyLossProcess::GetContinuousStepLimit( 1237 const G4Track&, 1945 const G4Track&, 1238 G4double, G4double, G4double& 1946 G4double, G4double, G4double&) 1239 { 1947 { 1240 return DBL_MAX; 1948 return DBL_MAX; 1241 } 1949 } 1242 1950 1243 //....oooOO0OOooo........oooOO0OOooo........o 1951 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1244 1952 1245 G4PhysicsVector* 1953 G4PhysicsVector* 1246 G4VEnergyLossProcess::LambdaPhysicsVector(con << 1954 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*, 1247 G4d 1955 G4double) 1248 { 1956 { 1249 DefineMaterial(couple); << 1957 G4PhysicsVector* v = 1250 G4PhysicsVector* v = (*theLambdaTable)[base << 1958 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins); 1251 return new G4PhysicsVector(*v); << 1959 v->SetSpline(theParameters->Spline()); >> 1960 return v; >> 1961 } >> 1962 >> 1963 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1964 >> 1965 void G4VEnergyLossProcess::AddCollaborativeProcess( >> 1966 G4VEnergyLossProcess* p) >> 1967 { >> 1968 G4bool add = true; >> 1969 if(p->GetProcessName() != "eBrem") { add = false; } >> 1970 if(add && nProcesses > 0) { >> 1971 for(G4int i=0; i<nProcesses; ++i) { >> 1972 if(p == scProcesses[i]) { >> 1973 add = false; >> 1974 break; >> 1975 } >> 1976 } >> 1977 } >> 1978 if(add) { >> 1979 scProcesses.push_back(p); >> 1980 ++nProcesses; >> 1981 if (1 < verboseLevel) { >> 1982 G4cout << "### The process " << p->GetProcessName() >> 1983 << " is added to the list of collaborative processes of " >> 1984 << GetProcessName() << G4endl; >> 1985 } >> 1986 } 1252 } 1987 } 1253 1988 1254 //....oooOO0OOooo........oooOO0OOooo........o 1989 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1255 1990 1256 void 1991 void 1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT 1992 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType) 1258 { 1993 { 1259 if(1 < verboseLevel) { << 1260 G4cout << "### Set DEDX table " << p << " << 1261 << " " << theDEDXunRestrictedTable << << 1262 << " for " << particle->GetParticl << 1263 << " and process " << GetProcessNa << 1264 << " type=" << tType << " isIonisation:" << 1265 } << 1266 if(fTotal == tType) { 1994 if(fTotal == tType) { 1267 theDEDXunRestrictedTable = p; 1995 theDEDXunRestrictedTable = p; >> 1996 if(p) { >> 1997 size_t n = p->length(); >> 1998 G4PhysicsVector* pv = (*p)[0]; >> 1999 G4double emax = maxKinEnergyCSDA; >> 2000 >> 2001 G4LossTableBuilder* bld = lManager->GetTableBuilder(); >> 2002 theDensityFactor = bld->GetDensityFactors(); >> 2003 theDensityIdx = bld->GetCoupleIndexes(); >> 2004 >> 2005 for (size_t i=0; i<n; ++i) { >> 2006 G4double dedx = 0.0; >> 2007 pv = (*p)[i]; >> 2008 if(pv) { >> 2009 dedx = pv->Value(emax, idxDEDXunRestricted); >> 2010 } else { >> 2011 pv = (*p)[(*theDensityIdx)[i]]; >> 2012 if(pv) { >> 2013 dedx = >> 2014 pv->Value(emax, idxDEDXunRestricted)*(*theDensityFactor)[i]; >> 2015 } >> 2016 } >> 2017 theDEDXAtMaxEnergy[i] = dedx; >> 2018 //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= " >> 2019 // << dedx << G4endl; >> 2020 } >> 2021 } >> 2022 1268 } else if(fRestricted == tType) { 2023 } else if(fRestricted == tType) { >> 2024 /* >> 2025 G4cout<< "G4VEnergyLossProcess::SetDEDXTable " >> 2026 << particle->GetParticleName() >> 2027 << " oldTable " << theDEDXTable << " newTable " << p >> 2028 << " ion " << theIonisationTable >> 2029 << " IsMaster " << isMaster >> 2030 << " " << GetProcessName() << G4endl; >> 2031 G4cout << (*p) << G4endl; >> 2032 */ 1269 theDEDXTable = p; 2033 theDEDXTable = p; 1270 if(isMaster && nullptr == baseParticle) { << 2034 } else if(fSubRestricted == tType) { 1271 theData->UpdateTable(theDEDXTable, 0); << 2035 theDEDXSubTable = p; 1272 } << 1273 } else if(fIsIonisation == tType) { 2036 } else if(fIsIonisation == tType) { >> 2037 /* >> 2038 G4cout<< "G4VEnergyLossProcess::SetIonisationTable " >> 2039 << particle->GetParticleName() >> 2040 << " oldTable " << theDEDXTable << " newTable " << p >> 2041 << " ion " << theIonisationTable >> 2042 << " IsMaster " << isMaster >> 2043 << " " << GetProcessName() << G4endl; >> 2044 */ 1274 theIonisationTable = p; 2045 theIonisationTable = p; 1275 if(isMaster && nullptr == baseParticle) { << 2046 } else if(fIsSubIonisation == tType) { 1276 theData->UpdateTable(theIonisationTable << 2047 theIonisationSubTable = p; 1277 } << 1278 } 2048 } 1279 } 2049 } 1280 2050 1281 //....oooOO0OOooo........oooOO0OOooo........o 2051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1282 2052 1283 void G4VEnergyLossProcess::SetCSDARangeTable( 2053 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p) 1284 { 2054 { 1285 theCSDARangeTable = p; << 2055 theCSDARangeTable = p; >> 2056 >> 2057 if(p) { >> 2058 size_t n = p->length(); >> 2059 G4PhysicsVector* pv; >> 2060 G4double emax = maxKinEnergyCSDA; >> 2061 >> 2062 for (size_t i=0; i<n; ++i) { >> 2063 pv = (*p)[i]; >> 2064 G4double rmax = 0.0; >> 2065 if(pv) { rmax = pv->Value(emax, idxCSDA); } >> 2066 else { >> 2067 pv = (*p)[(*theDensityIdx)[i]]; >> 2068 if(pv) { rmax = pv->Value(emax, idxCSDA)/(*theDensityFactor)[i]; } >> 2069 } >> 2070 theRangeAtMaxEnergy[i] = rmax; >> 2071 //G4cout << "i= " << i << " Emax(MeV)= " << emax/MeV << " Rmax= " >> 2072 //<< rmax<< G4endl; >> 2073 } >> 2074 } 1286 } 2075 } 1287 2076 1288 //....oooOO0OOooo........oooOO0OOooo........o 2077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1289 2078 1290 void G4VEnergyLossProcess::SetRangeTableForLo 2079 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p) 1291 { 2080 { 1292 theRangeTableForLoss = p; 2081 theRangeTableForLoss = p; >> 2082 if(1 < verboseLevel) { >> 2083 G4cout << "### Set Range table " << p >> 2084 << " for " << particle->GetParticleName() >> 2085 << " and process " << GetProcessName() << G4endl; >> 2086 } >> 2087 } >> 2088 >> 2089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 2090 >> 2091 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p) >> 2092 { >> 2093 theSecondaryRangeTable = p; >> 2094 if(1 < verboseLevel) { >> 2095 G4cout << "### Set SecondaryRange table " << p >> 2096 << " for " << particle->GetParticleName() >> 2097 << " and process " << GetProcessName() << G4endl; >> 2098 } 1293 } 2099 } 1294 2100 1295 //....oooOO0OOooo........oooOO0OOooo........o 2101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1296 2102 1297 void G4VEnergyLossProcess::SetInverseRangeTab 2103 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p) 1298 { 2104 { 1299 theInverseRangeTable = p; 2105 theInverseRangeTable = p; >> 2106 if(1 < verboseLevel) { >> 2107 G4cout << "### Set InverseRange table " << p >> 2108 << " for " << particle->GetParticleName() >> 2109 << " and process " << GetProcessName() << G4endl; >> 2110 } 1300 } 2111 } 1301 2112 1302 //....oooOO0OOooo........oooOO0OOooo........o 2113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1303 2114 1304 void G4VEnergyLossProcess::SetLambdaTable(G4P 2115 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p) 1305 { 2116 { 1306 if(1 < verboseLevel) { 2117 if(1 < verboseLevel) { 1307 G4cout << "### Set Lambda table " << p << << 2118 G4cout << "### Set Lambda table " << p 1308 << " for " << particle->GetParticl 2119 << " for " << particle->GetParticleName() 1309 << " and process " << GetProcessNa 2120 << " and process " << GetProcessName() << G4endl; >> 2121 //G4cout << *p << G4endl; 1310 } 2122 } 1311 theLambdaTable = p; << 2123 theLambdaTable = p; 1312 tablesAreBuilt = true; 2124 tablesAreBuilt = true; 1313 2125 1314 if(isMaster && nullptr != p) { << 2126 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 1315 delete theEnergyOfCrossSectionMax; << 2127 theDensityFactor = bld->GetDensityFactors(); 1316 theEnergyOfCrossSectionMax = nullptr; << 2128 theDensityIdx = bld->GetCoupleIndexes(); 1317 if(fEmTwoPeaks == fXSType) { << 2129 1318 if(nullptr != fXSpeaks) { << 2130 if(theLambdaTable) { 1319 for(auto & ptr : *fXSpeaks) { delete ptr; } << 2131 size_t n = theLambdaTable->length(); 1320 delete fXSpeaks; << 2132 G4PhysicsVector* pv = (*theLambdaTable)[0]; >> 2133 G4double e, ss, smax, emax; >> 2134 >> 2135 size_t i; >> 2136 >> 2137 // first loop on existing vectors >> 2138 for (i=0; i<n; ++i) { >> 2139 pv = (*theLambdaTable)[i]; >> 2140 if(pv) { >> 2141 size_t nb = pv->GetVectorLength(); >> 2142 emax = DBL_MAX; >> 2143 smax = 0.0; >> 2144 if(nb > 0) { >> 2145 for (size_t j=0; j<nb; ++j) { >> 2146 e = pv->Energy(j); >> 2147 ss = (*pv)(j); >> 2148 if(ss > smax) { >> 2149 smax = ss; >> 2150 emax = e; >> 2151 } >> 2152 } >> 2153 } >> 2154 theEnergyOfCrossSectionMax[i] = emax; >> 2155 theCrossSectionMax[i] = smax; >> 2156 if(1 < verboseLevel) { >> 2157 G4cout << "For " << particle->GetParticleName() >> 2158 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV >> 2159 << " lambda= " << smax << G4endl; >> 2160 } 1321 } 2161 } 1322 G4LossTableBuilder* bld = lManager->Get << 1323 fXSpeaks = G4EmUtility::FillPeaksStruct << 1324 if(nullptr == fXSpeaks) { fXSType = fEm << 1325 } 2162 } 1326 if(fXSType == fEmOnePeak) { << 2163 // second loop using base materials 1327 theEnergyOfCrossSectionMax = G4EmUtilit << 2164 for (i=0; i<n; ++i) { 1328 if(nullptr == theEnergyOfCrossSectionMa << 2165 pv = (*theLambdaTable)[i]; >> 2166 if(!pv){ >> 2167 G4int j = (*theDensityIdx)[i]; >> 2168 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j]; >> 2169 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j]; >> 2170 } 1329 } 2171 } 1330 } 2172 } 1331 } 2173 } 1332 2174 1333 //....oooOO0OOooo........oooOO0OOooo........o 2175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1334 2176 1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe << 2177 void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p) 1336 { 2178 { 1337 theEnergyOfCrossSectionMax = p; << 2179 theSubLambdaTable = p; 1338 } << 2180 if(1 < verboseLevel) { 1339 << 2181 G4cout << "### Set SebLambda table " << p 1340 //....oooOO0OOooo........oooOO0OOooo........o << 2182 << " for " << particle->GetParticleName() 1341 << 2183 << " and process " << GetProcessName() << G4endl; 1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: << 2184 } 1343 { << 1344 fXSpeaks = ptr; << 1345 } 2185 } 1346 2186 1347 //....oooOO0OOooo........oooOO0OOooo........o 2187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1348 2188 1349 const G4Element* G4VEnergyLossProcess::GetCur 2189 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const 1350 { 2190 { 1351 return (nullptr != currentModel) << 2191 const G4Element* elm = nullptr; 1352 ? currentModel->GetCurrentElement(current << 2192 if(currentModel) { elm = currentModel->GetCurrentElement(); } >> 2193 return elm; 1353 } 2194 } 1354 2195 1355 //....oooOO0OOooo........oooOO0OOooo........o 2196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1356 2197 1357 void G4VEnergyLossProcess::SetCrossSectionBia 2198 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f, 1358 2199 G4bool flag) 1359 { 2200 { 1360 if(f > 0.0) { 2201 if(f > 0.0) { 1361 biasFactor = f; 2202 biasFactor = f; 1362 weightFlag = flag; 2203 weightFlag = flag; 1363 if(1 < verboseLevel) { 2204 if(1 < verboseLevel) { 1364 G4cout << "### SetCrossSectionBiasingFa 2205 G4cout << "### SetCrossSectionBiasingFactor: for " 1365 << " process " << GetProcessName 2206 << " process " << GetProcessName() 1366 << " biasFactor= " << f << " wei 2207 << " biasFactor= " << f << " weightFlag= " << flag 1367 << G4endl; 2208 << G4endl; 1368 } 2209 } 1369 } 2210 } 1370 } 2211 } 1371 2212 1372 //....oooOO0OOooo........oooOO0OOooo........o 2213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1373 2214 1374 void G4VEnergyLossProcess::ActivateForcedInte << 2215 void 1375 << 2216 G4VEnergyLossProcess::ActivateForcedInteraction(G4double length, 1376 << 2217 const G4String& region, >> 2218 G4bool flag) 1377 { 2219 { 1378 if(nullptr == biasManager) { biasManager = << 2220 if(!biasManager) { biasManager = new G4EmBiasingManager(); } 1379 if(1 < verboseLevel) { 2221 if(1 < verboseLevel) { 1380 G4cout << "### ActivateForcedInteraction: 2222 G4cout << "### ActivateForcedInteraction: for " 1381 << " process " << GetProcessName() 2223 << " process " << GetProcessName() 1382 << " length(mm)= " << length/mm 2224 << " length(mm)= " << length/mm 1383 << " in G4Region <" << region 2225 << " in G4Region <" << region 1384 << "> weightFlag= " << flag 2226 << "> weightFlag= " << flag 1385 << G4endl; 2227 << G4endl; 1386 } 2228 } 1387 weightFlag = flag; 2229 weightFlag = flag; 1388 biasManager->ActivateForcedInteraction(leng 2230 biasManager->ActivateForcedInteraction(length, region); 1389 } 2231 } 1390 2232 1391 //....oooOO0OOooo........oooOO0OOooo........o 2233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1392 2234 1393 void 2235 void 1394 G4VEnergyLossProcess::ActivateSecondaryBiasin 2236 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region, 1395 2237 G4double factor, 1396 2238 G4double energyLimit) 1397 { 2239 { 1398 if (0.0 <= factor) { 2240 if (0.0 <= factor) { >> 2241 1399 // Range cut can be applied only for e- 2242 // Range cut can be applied only for e- 1400 if(0.0 == factor && secondaryParticle != 2243 if(0.0 == factor && secondaryParticle != G4Electron::Electron()) 1401 { return; } 2244 { return; } 1402 2245 1403 if(nullptr == biasManager) { biasManager << 2246 if(!biasManager) { biasManager = new G4EmBiasingManager(); } 1404 biasManager->ActivateSecondaryBiasing(reg 2247 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit); 1405 if(1 < verboseLevel) { 2248 if(1 < verboseLevel) { 1406 G4cout << "### ActivateSecondaryBiasing 2249 G4cout << "### ActivateSecondaryBiasing: for " 1407 << " process " << GetProcessName 2250 << " process " << GetProcessName() 1408 << " factor= " << factor 2251 << " factor= " << factor 1409 << " in G4Region <" << region 2252 << " in G4Region <" << region 1410 << "> energyLimit(MeV)= " << ene 2253 << "> energyLimit(MeV)= " << energyLimit/MeV 1411 << G4endl; 2254 << G4endl; 1412 } 2255 } 1413 } 2256 } 1414 } 2257 } 1415 2258 1416 //....oooOO0OOooo........oooOO0OOooo........o 2259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1417 2260 1418 void G4VEnergyLossProcess::SetIonisation(G4bo 2261 void G4VEnergyLossProcess::SetIonisation(G4bool val) 1419 { 2262 { 1420 isIonisation = val; 2263 isIonisation = val; 1421 aGPILSelection = (val) ? CandidateForSelect 2264 aGPILSelection = (val) ? CandidateForSelection : NotCandidateForSelection; 1422 } 2265 } 1423 2266 1424 //....oooOO0OOooo........oooOO0OOooo........o 2267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1425 2268 1426 void G4VEnergyLossProcess::SetLinearLossLimi 2269 void G4VEnergyLossProcess::SetLinearLossLimit(G4double val) 1427 { 2270 { 1428 if(0.0 < val && val < 1.0) { 2271 if(0.0 < val && val < 1.0) { 1429 linLossLimit = val; 2272 linLossLimit = val; 1430 actLinLossLimit = true; 2273 actLinLossLimit = true; 1431 } else { PrintWarning("SetLinearLossLimit", 2274 } else { PrintWarning("SetLinearLossLimit", val); } 1432 } 2275 } 1433 2276 1434 //....oooOO0OOooo........oooOO0OOooo........o 2277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1435 2278 1436 void G4VEnergyLossProcess::SetStepFunction(G4 << 2279 void >> 2280 G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2, G4bool lock) 1437 { 2281 { 1438 if(0.0 < v1 && 0.0 < v2) { << 2282 if(actStepFunc) { return; } >> 2283 actStepFunc = lock; >> 2284 if(0.0 < v1 && 0.0 < v2 && v2 < 1.e+50) { 1439 dRoverRange = std::min(1.0, v1); 2285 dRoverRange = std::min(1.0, v1); 1440 finalRange = std::min(v2, 1.e+50); << 2286 finalRange = v2; >> 2287 } else if(v1 <= 0.0) { >> 2288 PrintWarning("SetStepFunction", v1); 1441 } else { 2289 } else { 1442 PrintWarning("SetStepFunctionV1", v1); << 2290 PrintWarning("SetStepFunction", v2); 1443 PrintWarning("SetStepFunctionV2", v2); << 1444 } 2291 } 1445 } 2292 } 1446 2293 1447 //....oooOO0OOooo........oooOO0OOooo........o 2294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1448 2295 1449 void G4VEnergyLossProcess::SetLowestEnergyLim 2296 void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val) 1450 { 2297 { 1451 if(1.e-18 < val && val < 1.e+50) { lowestKi 2298 if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; } 1452 else { PrintWarning("SetLowestEnergyLimit", 2299 else { PrintWarning("SetLowestEnergyLimit", val); } 1453 } 2300 } 1454 2301 1455 //....oooOO0OOooo........oooOO0OOooo........o 2302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1456 2303 1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i 2304 void G4VEnergyLossProcess::SetDEDXBinning(G4int n) 1458 { 2305 { 1459 if(2 < n && n < 1000000000) { 2306 if(2 < n && n < 1000000000) { 1460 nBins = n; 2307 nBins = n; 1461 actBinning = true; 2308 actBinning = true; 1462 } else { 2309 } else { 1463 G4double e = (G4double)n; 2310 G4double e = (G4double)n; 1464 PrintWarning("SetDEDXBinning", e); 2311 PrintWarning("SetDEDXBinning", e); 1465 } 2312 } 1466 } 2313 } 1467 2314 1468 //....oooOO0OOooo........oooOO0OOooo........o 2315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1469 2316 1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 2317 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e) 1471 { 2318 { 1472 if(1.e-18 < e && e < maxKinEnergy) { 2319 if(1.e-18 < e && e < maxKinEnergy) { 1473 minKinEnergy = e; 2320 minKinEnergy = e; 1474 actMinKinEnergy = true; 2321 actMinKinEnergy = true; 1475 } else { PrintWarning("SetMinKinEnergy", e) 2322 } else { PrintWarning("SetMinKinEnergy", e); } 1476 } 2323 } 1477 2324 1478 //....oooOO0OOooo........oooOO0OOooo........o 2325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1479 2326 1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 2327 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e) 1481 { 2328 { 1482 if(minKinEnergy < e && e < 1.e+50) { 2329 if(minKinEnergy < e && e < 1.e+50) { 1483 maxKinEnergy = e; 2330 maxKinEnergy = e; 1484 actMaxKinEnergy = true; 2331 actMaxKinEnergy = true; 1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS 2332 if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; } 1486 } else { PrintWarning("SetMaxKinEnergy", e) 2333 } else { PrintWarning("SetMaxKinEnergy", e); } 1487 } 2334 } 1488 2335 1489 //....oooOO0OOooo........oooOO0OOooo........o 2336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1490 2337 1491 void G4VEnergyLossProcess::PrintWarning(const << 2338 void G4VEnergyLossProcess::PrintWarning(G4String tit, G4double val) 1492 { 2339 { 1493 G4String ss = "G4VEnergyLossProcess::" + ti 2340 G4String ss = "G4VEnergyLossProcess::" + tit; 1494 G4ExceptionDescription ed; 2341 G4ExceptionDescription ed; 1495 ed << "Parameter is out of range: " << val 2342 ed << "Parameter is out of range: " << val 1496 << " it will have no effect!\n" << " Pr 2343 << " it will have no effect!\n" << " Process " 1497 << GetProcessName() << " nbins= " << nB 2344 << GetProcessName() << " nbins= " << nBins 1498 << " Emin(keV)= " << minKinEnergy/keV 2345 << " Emin(keV)= " << minKinEnergy/keV 1499 << " Emax(GeV)= " << maxKinEnergy/GeV; 2346 << " Emax(GeV)= " << maxKinEnergy/GeV; 1500 G4Exception(ss, "em0044", JustWarning, ed); 2347 G4Exception(ss, "em0044", JustWarning, ed); 1501 } 2348 } 1502 2349 1503 //....oooOO0OOooo........oooOO0OOooo........o 2350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1504 2351 1505 void G4VEnergyLossProcess::ProcessDescription 2352 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const 1506 { 2353 { 1507 if(nullptr != particle) { StreamInfo(out, * << 2354 if(particle) { StreamInfo(out, *particle, true); } 1508 } 2355 } 1509 2356 1510 //....oooOO0OOooo........oooOO0OOooo........o 2357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1511 2358