Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // ------------------------------------------- 26 // ------------------------------------------------------------------- 27 // 27 // 28 // GEANT4 Class file 28 // GEANT4 Class file 29 // 29 // 30 // 30 // 31 // File name: G4VEnergyLossProcess 31 // File name: G4VEnergyLossProcess 32 // 32 // 33 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko 34 // 34 // 35 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002 36 // 36 // 37 // Modifications: Vladimir Ivanchenko 37 // Modifications: Vladimir Ivanchenko 38 // 38 // 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // 41 // 42 // It is the unified energy loss process it ca 42 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 43 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 44 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 45 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 46 // that information on fly. 47 // ------------------------------------------- 47 // ------------------------------------------------------------------- 48 // 48 // 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 51 52 #include "G4VEnergyLossProcess.hh" 52 #include "G4VEnergyLossProcess.hh" 53 #include "G4PhysicalConstants.hh" 53 #include "G4PhysicalConstants.hh" 54 #include "G4SystemOfUnits.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "G4ProcessManager.hh" 55 #include "G4ProcessManager.hh" 56 #include "G4LossTableManager.hh" 56 #include "G4LossTableManager.hh" 57 #include "G4LossTableBuilder.hh" 57 #include "G4LossTableBuilder.hh" 58 #include "G4Step.hh" 58 #include "G4Step.hh" 59 #include "G4ParticleDefinition.hh" 59 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" 60 #include "G4ParticleTable.hh" 61 #include "G4EmParameters.hh" 61 #include "G4EmParameters.hh" 62 #include "G4EmUtility.hh" 62 #include "G4EmUtility.hh" 63 #include "G4EmTableUtil.hh" 63 #include "G4EmTableUtil.hh" 64 #include "G4VEmModel.hh" 64 #include "G4VEmModel.hh" 65 #include "G4VEmFluctuationModel.hh" 65 #include "G4VEmFluctuationModel.hh" 66 #include "G4DataVector.hh" 66 #include "G4DataVector.hh" 67 #include "G4PhysicsLogVector.hh" 67 #include "G4PhysicsLogVector.hh" 68 #include "G4VParticleChange.hh" 68 #include "G4VParticleChange.hh" 69 #include "G4Electron.hh" 69 #include "G4Electron.hh" 70 #include "G4ProcessManager.hh" 70 #include "G4ProcessManager.hh" 71 #include "G4UnitsTable.hh" 71 #include "G4UnitsTable.hh" 72 #include "G4Region.hh" 72 #include "G4Region.hh" 73 #include "G4RegionStore.hh" 73 #include "G4RegionStore.hh" 74 #include "G4PhysicsTableHelper.hh" 74 #include "G4PhysicsTableHelper.hh" 75 #include "G4SafetyHelper.hh" 75 #include "G4SafetyHelper.hh" 76 #include "G4EmDataHandler.hh" 76 #include "G4EmDataHandler.hh" 77 #include "G4TransportationManager.hh" 77 #include "G4TransportationManager.hh" 78 #include "G4VAtomDeexcitation.hh" 78 #include "G4VAtomDeexcitation.hh" 79 #include "G4VSubCutProducer.hh" 79 #include "G4VSubCutProducer.hh" 80 #include "G4EmBiasingManager.hh" 80 #include "G4EmBiasingManager.hh" 81 #include "G4Log.hh" 81 #include "G4Log.hh" 82 #include <iostream> 82 #include <iostream> 83 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 85 86 namespace 86 namespace 87 { 87 { 88 G4String tnames[7] = 88 G4String tnames[7] = 89 {"DEDX","Ionisation","DEDXnr","CSDARange", 89 {"DEDX","Ionisation","DEDXnr","CSDARange","Lambda","Range","InverseRange"}; 90 } 90 } 91 91 92 92 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con 93 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 94 G4P 94 G4ProcessType type): 95 G4VContinuousDiscreteProcess(name, type) 95 G4VContinuousDiscreteProcess(name, type) 96 { 96 { 97 theParameters = G4EmParameters::Instance(); 97 theParameters = G4EmParameters::Instance(); 98 SetVerboseLevel(1); 98 SetVerboseLevel(1); 99 99 100 // low energy limit 100 // low energy limit 101 lowestKinEnergy = theParameters->LowestElect 101 lowestKinEnergy = theParameters->LowestElectronEnergy(); 102 102 103 // Size of tables 103 // Size of tables 104 minKinEnergy = 0.1*CLHEP::keV; 104 minKinEnergy = 0.1*CLHEP::keV; 105 maxKinEnergy = 100.0*CLHEP::TeV; 105 maxKinEnergy = 100.0*CLHEP::TeV; 106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; 106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; 107 nBins = 84; 107 nBins = 84; 108 nBinsCSDA = 35; 108 nBinsCSDA = 35; 109 109 110 invLambdaFactor = 1.0/lambdaFactor; 110 invLambdaFactor = 1.0/lambdaFactor; 111 111 112 // default linear loss limit 112 // default linear loss limit 113 finalRange = 1.*CLHEP::mm; 113 finalRange = 1.*CLHEP::mm; 114 114 115 // run time objects 115 // run time objects 116 pParticleChange = &fParticleChange; 116 pParticleChange = &fParticleChange; 117 fParticleChange.SetSecondaryWeightByProcess( 117 fParticleChange.SetSecondaryWeightByProcess(true); 118 modelManager = new G4EmModelManager(); 118 modelManager = new G4EmModelManager(); 119 safetyHelper = G4TransportationManager::GetT 119 safetyHelper = G4TransportationManager::GetTransportationManager() 120 ->GetSafetyHelper(); 120 ->GetSafetyHelper(); 121 aGPILSelection = CandidateForSelection; 121 aGPILSelection = CandidateForSelection; 122 122 123 // initialise model 123 // initialise model 124 lManager = G4LossTableManager::Instance(); 124 lManager = G4LossTableManager::Instance(); 125 lManager->Register(this); 125 lManager->Register(this); 126 isMaster = lManager->IsMaster(); 126 isMaster = lManager->IsMaster(); 127 127 128 G4LossTableBuilder* bld = lManager->GetTable 128 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 129 theDensityFactor = bld->GetDensityFactors(); 129 theDensityFactor = bld->GetDensityFactors(); 130 theDensityIdx = bld->GetCoupleIndexes(); 130 theDensityIdx = bld->GetCoupleIndexes(); 131 131 132 scTracks.reserve(10); 132 scTracks.reserve(10); 133 secParticles.reserve(12); 133 secParticles.reserve(12); 134 emModels = new std::vector<G4VEmModel*>; 134 emModels = new std::vector<G4VEmModel*>; 135 } 135 } 136 136 137 //....oooOO0OOooo........oooOO0OOooo........oo 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 138 139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 140 { 140 { 141 if (isMaster) { 141 if (isMaster) { 142 if(nullptr == baseParticle) { delete theDa 142 if(nullptr == baseParticle) { delete theData; } 143 delete theEnergyOfCrossSectionMax; 143 delete theEnergyOfCrossSectionMax; 144 if(nullptr != fXSpeaks) { 144 if(nullptr != fXSpeaks) { 145 for(auto const & v : *fXSpeaks) { delete 145 for(auto const & v : *fXSpeaks) { delete v; } 146 delete fXSpeaks; 146 delete fXSpeaks; 147 } 147 } 148 } 148 } 149 delete modelManager; 149 delete modelManager; 150 delete biasManager; 150 delete biasManager; 151 delete scoffRegions; 151 delete scoffRegions; 152 delete emModels; 152 delete emModels; 153 lManager->DeRegister(this); 153 lManager->DeRegister(this); 154 } 154 } 155 155 156 //....oooOO0OOooo........oooOO0OOooo........oo 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 157 158 G4double G4VEnergyLossProcess::MinPrimaryEnerg 158 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*, 159 159 const G4Material*, 160 160 G4double cut) 161 { 161 { 162 return cut; 162 return cut; 163 } 163 } 164 164 165 //....oooOO0OOooo........oooOO0OOooo........oo 165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 166 166 167 void G4VEnergyLossProcess::AddEmModel(G4int or 167 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* ptr, 168 G4VEmFlu 168 G4VEmFluctuationModel* fluc, 169 const G4 169 const G4Region* region) 170 { 170 { 171 if(nullptr == ptr) { return; } 171 if(nullptr == ptr) { return; } 172 G4VEmFluctuationModel* afluc = (nullptr == f 172 G4VEmFluctuationModel* afluc = (nullptr == fluc) ? fluctModel : fluc; 173 modelManager->AddEmModel(order, ptr, afluc, 173 modelManager->AddEmModel(order, ptr, afluc, region); 174 ptr->SetParticleChange(pParticleChange, aflu 174 ptr->SetParticleChange(pParticleChange, afluc); 175 } 175 } 176 176 177 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 178 179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod 179 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* ptr, G4int) 180 { 180 { 181 if(nullptr == ptr) { return; } 181 if(nullptr == ptr) { return; } 182 if(!emModels->empty()) { 182 if(!emModels->empty()) { 183 for(auto & em : *emModels) { if(em == ptr) 183 for(auto & em : *emModels) { if(em == ptr) { return; } } 184 } 184 } 185 emModels->push_back(ptr); 185 emModels->push_back(ptr); 186 } 186 } 187 187 188 //....oooOO0OOooo........oooOO0OOooo........oo 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 189 189 190 void G4VEnergyLossProcess::SetDynamicMassCharg 190 void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio, 191 191 G4double charge2ratio) 192 { 192 { 193 massRatio = massratio; 193 massRatio = massratio; 194 logMassRatio = G4Log(massRatio); 194 logMassRatio = G4Log(massRatio); 195 fFactor = charge2ratio*biasFactor; 195 fFactor = charge2ratio*biasFactor; 196 if(baseMat) { fFactor *= (*theDensityFactor) 196 if(baseMat) { fFactor *= (*theDensityFactor)[currentCoupleIndex]; } 197 chargeSqRatio = charge2ratio; 197 chargeSqRatio = charge2ratio; 198 reduceFactor = 1.0/(fFactor*massRatio); 198 reduceFactor = 1.0/(fFactor*massRatio); 199 } 199 } 200 200 201 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 202 202 203 void 203 void 204 G4VEnergyLossProcess::PreparePhysicsTable(cons 204 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part) 205 { 205 { 206 particle = G4EmTableUtil::CheckIon(this, &pa 206 particle = G4EmTableUtil::CheckIon(this, &part, particle, 207 verboseLe 207 verboseLevel, isIon); 208 208 209 if( particle != &part ) { 209 if( particle != &part ) { 210 if(!isIon) { lManager->RegisterExtraPartic 210 if(!isIon) { lManager->RegisterExtraParticle(&part, this); } 211 if(1 < verboseLevel) { 211 if(1 < verboseLevel) { 212 G4cout << "### G4VEnergyLossProcess::Pre 212 G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()" 213 << " interrupted for " << GetProc << 213 << " interrupted for " 214 << part.GetParticleName() << " is << 214 << part.GetParticleName() << " isIon=" << isIon << G4endl; 215 << " spline=" << spline << G4endl << 216 } 215 } 217 return; 216 return; 218 } 217 } 219 218 220 tablesAreBuilt = false; 219 tablesAreBuilt = false; 221 if (GetProcessSubType() == fIonisation) { Se << 222 220 223 G4LossTableBuilder* bld = lManager->GetTable 221 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 224 lManager->PreparePhysicsTable(&part, this); << 222 lManager->PreparePhysicsTable(&part, this, isMaster); 225 223 226 // Base particle and set of models can be de 224 // Base particle and set of models can be defined here 227 InitialiseEnergyLossProcess(particle, basePa 225 InitialiseEnergyLossProcess(particle, baseParticle); 228 226 229 // parameters of the process 227 // parameters of the process 230 if(!actLossFluc) { lossFluctuationFlag = the 228 if(!actLossFluc) { lossFluctuationFlag = theParameters->LossFluctuation(); } 231 useCutAsFinalRange = theParameters->UseCutAs << 229 rndmStepFlag = theParameters->UseCutAsFinalRange(); 232 if(!actMinKinEnergy) { minKinEnergy = thePar 230 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); } 233 if(!actMaxKinEnergy) { maxKinEnergy = thePar 231 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); } 234 if(!actBinning) { nBins = theParameters->Num 232 if(!actBinning) { nBins = theParameters->NumberOfBins(); } 235 maxKinEnergyCSDA = theParameters->MaxEnergyF 233 maxKinEnergyCSDA = theParameters->MaxEnergyForCSDARange(); 236 nBinsCSDA = theParameters->NumberOfBinsPerDe 234 nBinsCSDA = theParameters->NumberOfBinsPerDecade() 237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi 235 *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy)); 238 if(!actLinLossLimit) { linLossLimit = thePar 236 if(!actLinLossLimit) { linLossLimit = theParameters->LinearLossLimit(); } 239 lambdaFactor = theParameters->LambdaFactor() 237 lambdaFactor = theParameters->LambdaFactor(); 240 invLambdaFactor = 1.0/lambdaFactor; 238 invLambdaFactor = 1.0/lambdaFactor; 241 if(isMaster) { SetVerboseLevel(theParameters 239 if(isMaster) { SetVerboseLevel(theParameters->Verbose()); } 242 else { SetVerboseLevel(theParameters->Worker 240 else { SetVerboseLevel(theParameters->WorkerVerbose()); } 243 // integral option may be disabled 241 // integral option may be disabled 244 if(!theParameters->Integral()) { fXSType = f 242 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; } 245 243 246 theParameters->DefineRegParamForLoss(this); 244 theParameters->DefineRegParamForLoss(this); 247 245 248 fRangeEnergy = 0.0; 246 fRangeEnergy = 0.0; 249 247 250 G4double initialCharge = particle->GetPDGCha 248 G4double initialCharge = particle->GetPDGCharge(); 251 G4double initialMass = particle->GetPDGMas 249 G4double initialMass = particle->GetPDGMass(); 252 250 253 theParameters->FillStepFunction(particle, th 251 theParameters->FillStepFunction(particle, this); 254 252 255 // parameters for scaling from the base part 253 // parameters for scaling from the base particle 256 if (nullptr != baseParticle) { 254 if (nullptr != baseParticle) { 257 massRatio = (baseParticle->GetPDGMass() 255 massRatio = (baseParticle->GetPDGMass())/initialMass; 258 logMassRatio = G4Log(massRatio); 256 logMassRatio = G4Log(massRatio); 259 G4double q = initialCharge/baseParticle->G 257 G4double q = initialCharge/baseParticle->GetPDGCharge(); 260 chargeSqRatio = q*q; 258 chargeSqRatio = q*q; 261 if(chargeSqRatio > 0.0) { reduceFactor = 1 259 if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); } 262 } 260 } 263 lowestKinEnergy = (initialMass < CLHEP::MeV) 261 lowestKinEnergy = (initialMass < CLHEP::MeV) 264 ? theParameters->LowestElectronEnergy() 262 ? theParameters->LowestElectronEnergy() 265 : theParameters->LowestMuHadEnergy(); 263 : theParameters->LowestMuHadEnergy(); 266 264 267 // Tables preparation 265 // Tables preparation 268 if (isMaster && nullptr == baseParticle) { 266 if (isMaster && nullptr == baseParticle) { 269 if(nullptr == theData) { theData = new G4E 267 if(nullptr == theData) { theData = new G4EmDataHandler(7); } 270 268 271 if(nullptr != theDEDXTable && isIonisation 269 if(nullptr != theDEDXTable && isIonisation) { 272 if(nullptr != theIonisationTable && theD 270 if(nullptr != theIonisationTable && theDEDXTable != theIonisationTable) { 273 theData->CleanTable(0); 271 theData->CleanTable(0); 274 theDEDXTable = theIonisationTable; 272 theDEDXTable = theIonisationTable; 275 theIonisationTable = nullptr; 273 theIonisationTable = nullptr; 276 } 274 } 277 } 275 } 278 276 279 theDEDXTable = theData->MakeTable(theDEDXT 277 theDEDXTable = theData->MakeTable(theDEDXTable, 0); 280 bld->InitialiseBaseMaterials(theDEDXTable) 278 bld->InitialiseBaseMaterials(theDEDXTable); 281 theData->UpdateTable(theIonisationTable, 1 279 theData->UpdateTable(theIonisationTable, 1); 282 280 283 if (theParameters->BuildCSDARange()) { 281 if (theParameters->BuildCSDARange()) { 284 theDEDXunRestrictedTable = theData->Make 282 theDEDXunRestrictedTable = theData->MakeTable(2); 285 if(isIonisation) { theCSDARangeTable = t 283 if(isIonisation) { theCSDARangeTable = theData->MakeTable(3); } 286 } 284 } 287 285 288 theLambdaTable = theData->MakeTable(4); 286 theLambdaTable = theData->MakeTable(4); 289 if(isIonisation) { 287 if(isIonisation) { 290 theRangeTableForLoss = theData->MakeTabl 288 theRangeTableForLoss = theData->MakeTable(5); 291 theInverseRangeTable = theData->MakeTabl 289 theInverseRangeTable = theData->MakeTable(6); 292 } 290 } 293 } 291 } 294 292 295 // forced biasing 293 // forced biasing 296 if(nullptr != biasManager) { 294 if(nullptr != biasManager) { 297 biasManager->Initialise(part,GetProcessNam 295 biasManager->Initialise(part,GetProcessName(),verboseLevel); 298 biasFlag = false; 296 biasFlag = false; 299 } 297 } 300 baseMat = bld->GetBaseMaterialFlag(); 298 baseMat = bld->GetBaseMaterialFlag(); 301 numberOfModels = modelManager->NumberOfModel 299 numberOfModels = modelManager->NumberOfModels(); 302 currentModel = modelManager->GetModel(0); 300 currentModel = modelManager->GetModel(0); 303 G4EmTableUtil::UpdateModels(this, modelManag 301 G4EmTableUtil::UpdateModels(this, modelManager, maxKinEnergy, 304 numberOfModels, 302 numberOfModels, secID, biasID, 305 mainSecondaries, 303 mainSecondaries, baseMat, isMaster, 306 theParameters->U 304 theParameters->UseAngularGeneratorForIonisation()); 307 theCuts = modelManager->Initialise(particle, 305 theCuts = modelManager->Initialise(particle, secondaryParticle, 308 verboseLe 306 verboseLevel); 309 // subcut processor 307 // subcut processor 310 if(isIonisation) { 308 if(isIonisation) { 311 subcutProducer = lManager->SubCutProducer( 309 subcutProducer = lManager->SubCutProducer(); 312 } 310 } 313 if(1 == nSCoffRegions) { 311 if(1 == nSCoffRegions) { 314 if((*scoffRegions)[0]->GetName() == "Defau 312 if((*scoffRegions)[0]->GetName() == "DefaultRegionForTheWorld") { 315 delete scoffRegions; 313 delete scoffRegions; 316 scoffRegions = nullptr; 314 scoffRegions = nullptr; 317 nSCoffRegions = 0; 315 nSCoffRegions = 0; 318 } 316 } 319 } 317 } 320 318 321 if(1 < verboseLevel) { 319 if(1 < verboseLevel) { 322 G4cout << "G4VEnergyLossProcess::PrepearPh 320 G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done " 323 << " for " << GetProcessName() << " << 321 << " for local " << particle->GetParticleName() 324 << " isIon= " << isIon << " spline= << 322 << " isIon= " << isIon; 325 if(baseParticle) { 323 if(baseParticle) { 326 G4cout << "; base: " << baseParticle->Ge 324 G4cout << "; base: " << baseParticle->GetParticleName(); 327 } 325 } 328 G4cout << G4endl; << 329 G4cout << " chargeSqRatio= " << chargeSqRa 326 G4cout << " chargeSqRatio= " << chargeSqRatio 330 << " massRatio= " << massRatio 327 << " massRatio= " << massRatio 331 << " reduceFactor= " << reduceFacto 328 << " reduceFactor= " << reduceFactor << G4endl; 332 if (nSCoffRegions > 0) { 329 if (nSCoffRegions > 0) { 333 G4cout << " SubCut secondary production 330 G4cout << " SubCut secondary production is ON for regions: " << G4endl; 334 for (G4int i=0; i<nSCoffRegions; ++i) { 331 for (G4int i=0; i<nSCoffRegions; ++i) { 335 const G4Region* r = (*scoffRegions)[i] 332 const G4Region* r = (*scoffRegions)[i]; 336 G4cout << " " << r->GetName( 333 G4cout << " " << r->GetName() << G4endl; 337 } 334 } 338 } else if(nullptr != subcutProducer) { 335 } else if(nullptr != subcutProducer) { 339 G4cout << " SubCut secondary production 336 G4cout << " SubCut secondary production is ON for all regions" << G4endl; 340 } 337 } 341 } 338 } 342 } 339 } 343 340 344 //....oooOO0OOooo........oooOO0OOooo........oo 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 345 342 346 void G4VEnergyLossProcess::BuildPhysicsTable(c 343 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 347 { 344 { 348 if(1 < verboseLevel) { 345 if(1 < verboseLevel) { 349 G4cout << "### G4VEnergyLossProcess::Build 346 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for " 350 << GetProcessName() 347 << GetProcessName() 351 << " and particle " << part.GetPart 348 << " and particle " << part.GetParticleName() 352 << "; the first particle " << parti << 349 << "; local: " << particle->GetParticleName(); 353 if(baseParticle) { 350 if(baseParticle) { 354 G4cout << "; base: " << baseParticle->Ge 351 G4cout << "; base: " << baseParticle->GetParticleName(); 355 } 352 } 356 G4cout << G4endl; << 353 G4cout << " TablesAreBuilt= " << tablesAreBuilt 357 G4cout << " TablesAreBuilt= " << tables << 354 << " isIon= " << isIon << " " << this << G4endl; 358 << " spline=" << spline << " ptr: " << 359 } 355 } 360 356 361 if(&part == particle) { 357 if(&part == particle) { >> 358 362 if(isMaster) { 359 if(isMaster) { 363 lManager->BuildPhysicsTable(particle, th 360 lManager->BuildPhysicsTable(particle, this); 364 361 365 } else { 362 } else { 366 const auto masterProcess = 363 const auto masterProcess = 367 static_cast<const G4VEnergyLossProcess 364 static_cast<const G4VEnergyLossProcess*>(GetMasterProcess()); 368 365 369 numberOfModels = modelManager->NumberOfM 366 numberOfModels = modelManager->NumberOfModels(); 370 G4EmTableUtil::BuildLocalElossProcess(th 367 G4EmTableUtil::BuildLocalElossProcess(this, masterProcess, 371 pa 368 particle, numberOfModels); 372 tablesAreBuilt = true; 369 tablesAreBuilt = true; 373 baseMat = masterProcess->UseBaseMaterial 370 baseMat = masterProcess->UseBaseMaterial(); 374 lManager->LocalPhysicsTables(particle, t 371 lManager->LocalPhysicsTables(particle, this); 375 } 372 } 376 373 377 // needs to be done only once 374 // needs to be done only once 378 safetyHelper->InitialiseHelper(); 375 safetyHelper->InitialiseHelper(); 379 } 376 } 380 // Added tracking cut to avoid tracking arti 377 // Added tracking cut to avoid tracking artifacts 381 // and identified deexcitation flag 378 // and identified deexcitation flag 382 if(isIonisation) { 379 if(isIonisation) { 383 atomDeexcitation = lManager->AtomDeexcitat 380 atomDeexcitation = lManager->AtomDeexcitation(); 384 if(nullptr != atomDeexcitation) { 381 if(nullptr != atomDeexcitation) { 385 if(atomDeexcitation->IsPIXEActive()) { u 382 if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; } 386 } 383 } 387 } 384 } 388 385 389 // protection against double printout 386 // protection against double printout 390 if(theParameters->IsPrintLocked()) { return; 387 if(theParameters->IsPrintLocked()) { return; } 391 388 392 // explicitly defined printout by particle n 389 // explicitly defined printout by particle name 393 G4String num = part.GetParticleName(); 390 G4String num = part.GetParticleName(); 394 if(1 < verboseLevel || 391 if(1 < verboseLevel || 395 (0 < verboseLevel && (num == "e-" || 392 (0 < verboseLevel && (num == "e-" || 396 num == "e+" || n 393 num == "e+" || num == "mu+" || 397 num == "mu-" || n 394 num == "mu-" || num == "proton"|| 398 num == "pi+" || n 395 num == "pi+" || num == "pi-" || 399 num == "kaon+" || n 396 num == "kaon+" || num == "kaon-" || 400 num == "alpha" || n 397 num == "alpha" || num == "anti_proton" || 401 num == "GenericIon" 398 num == "GenericIon"|| num == "alpha+" ))) { 402 StreamInfo(G4cout, part); 399 StreamInfo(G4cout, part); 403 } 400 } 404 if(1 < verboseLevel) { 401 if(1 < verboseLevel) { 405 G4cout << "### G4VEnergyLossProcess::Build 402 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for " 406 << GetProcessName() 403 << GetProcessName() 407 << " and particle " << part.GetPart 404 << " and particle " << part.GetParticleName(); 408 if(isIonisation) { G4cout << " isIonisati 405 if(isIonisation) { G4cout << " isIonisation flag=1"; } 409 G4cout << " baseMat=" << baseMat << G4endl 406 G4cout << " baseMat=" << baseMat << G4endl; 410 } 407 } 411 } 408 } 412 409 413 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 414 411 415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED 412 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType) 416 { 413 { 417 G4PhysicsTable* table = nullptr; 414 G4PhysicsTable* table = nullptr; 418 G4double emax = maxKinEnergy; 415 G4double emax = maxKinEnergy; 419 G4int bin = nBins; 416 G4int bin = nBins; 420 417 421 if(fTotal == tType) { 418 if(fTotal == tType) { 422 emax = maxKinEnergyCSDA; 419 emax = maxKinEnergyCSDA; 423 bin = nBinsCSDA; 420 bin = nBinsCSDA; 424 table = theDEDXunRestrictedTable; 421 table = theDEDXunRestrictedTable; 425 } else if(fRestricted == tType) { 422 } else if(fRestricted == tType) { 426 table = theDEDXTable; 423 table = theDEDXTable; 427 } else { 424 } else { 428 G4cout << "G4VEnergyLossProcess::BuildDEDX 425 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type " 429 << tType << G4endl; 426 << tType << G4endl; 430 } 427 } 431 if(1 < verboseLevel) { 428 if(1 < verboseLevel) { 432 G4cout << "G4VEnergyLossProcess::BuildDEDX 429 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType 433 << " for " << GetProcessName() 430 << " for " << GetProcessName() 434 << " and " << particle->GetParticle << 431 << " and " << particle->GetParticleName() << G4endl; 435 << "spline=" << spline << G4endl; << 436 } 432 } 437 if(nullptr == table) { return table; } 433 if(nullptr == table) { return table; } 438 434 439 G4LossTableBuilder* bld = lManager->GetTable 435 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 440 G4EmTableUtil::BuildDEDXTable(this, particle 436 G4EmTableUtil::BuildDEDXTable(this, particle, modelManager, bld, 441 table, minKinE 437 table, minKinEnergy, emax, bin, 442 verboseLevel, 438 verboseLevel, tType, spline); 443 return table; 439 return table; 444 } 440 } 445 441 446 //....oooOO0OOooo........oooOO0OOooo........oo 442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 447 443 448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam 444 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType) 449 { 445 { 450 if(nullptr == theLambdaTable) { return theLa 446 if(nullptr == theLambdaTable) { return theLambdaTable; } 451 447 452 G4double scale = theParameters->MaxKinEnergy 448 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy(); 453 G4int nbin = 449 G4int nbin = 454 theParameters->NumberOfBinsPerDecade()*G4l 450 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale)); 455 scale = nbin/G4Log(scale); 451 scale = nbin/G4Log(scale); 456 452 457 G4LossTableBuilder* bld = lManager->GetTable 453 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 458 G4EmTableUtil::BuildLambdaTable(this, partic 454 G4EmTableUtil::BuildLambdaTable(this, particle, modelManager, 459 bld, theLamb 455 bld, theLambdaTable, theCuts, 460 minKinEnergy 456 minKinEnergy, maxKinEnergy, scale, 461 verboseLevel 457 verboseLevel, spline); 462 return theLambdaTable; 458 return theLambdaTable; 463 } 459 } 464 460 465 //....oooOO0OOooo........oooOO0OOooo........oo 461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 466 462 467 void G4VEnergyLossProcess::StreamInfo(std::ost 463 void G4VEnergyLossProcess::StreamInfo(std::ostream& out, 468 const G4ParticleDefinition& pa 464 const G4ParticleDefinition& part, G4bool rst) const 469 { 465 { 470 G4String indent = (rst ? " " : ""); 466 G4String indent = (rst ? " " : ""); 471 out << std::setprecision(6); 467 out << std::setprecision(6); 472 out << G4endl << indent << GetProcessName() 468 out << G4endl << indent << GetProcessName() << ": "; 473 if (!rst) out << " for " << part.GetParticle 469 if (!rst) out << " for " << part.GetParticleName(); 474 out << " XStype:" << fXSType 470 out << " XStype:" << fXSType 475 << " SubType=" << GetProcessSubType() < 471 << " SubType=" << GetProcessSubType() << G4endl 476 << " dE/dx and range tables from " 472 << " dE/dx and range tables from " 477 << G4BestUnit(minKinEnergy,"Energy") 473 << G4BestUnit(minKinEnergy,"Energy") 478 << " to " << G4BestUnit(maxKinEnergy,"En 474 << " to " << G4BestUnit(maxKinEnergy,"Energy") 479 << " in " << nBins << " bins" << G4endl 475 << " in " << nBins << " bins" << G4endl 480 << " Lambda tables from threshold t 476 << " Lambda tables from threshold to " 481 << G4BestUnit(maxKinEnergy,"Energy") 477 << G4BestUnit(maxKinEnergy,"Energy") 482 << ", " << theParameters->NumberOfBinsPe 478 << ", " << theParameters->NumberOfBinsPerDecade() 483 << " bins/decade, spline: " << spline 479 << " bins/decade, spline: " << spline 484 << G4endl; 480 << G4endl; 485 if(nullptr != theRangeTableForLoss && isIoni 481 if(nullptr != theRangeTableForLoss && isIonisation) { 486 out << " StepFunction=(" << dRoverRan 482 out << " StepFunction=(" << dRoverRange << ", " 487 << finalRange/mm << " mm)" 483 << finalRange/mm << " mm)" 488 << ", integ: " << fXSType 484 << ", integ: " << fXSType 489 << ", fluct: " << lossFluctuationFlag 485 << ", fluct: " << lossFluctuationFlag 490 << ", linLossLim= " << linLossLimit 486 << ", linLossLim= " << linLossLimit 491 << G4endl; 487 << G4endl; 492 } 488 } 493 StreamProcessInfo(out); 489 StreamProcessInfo(out); 494 modelManager->DumpModelList(out, verboseLeve 490 modelManager->DumpModelList(out, verboseLevel); 495 if(nullptr != theCSDARangeTable && isIonisat 491 if(nullptr != theCSDARangeTable && isIonisation) { 496 out << " CSDA range table up" 492 out << " CSDA range table up" 497 << " to " << G4BestUnit(maxKinEnergyCS 493 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy") 498 << " in " << nBinsCSDA << " bins" << G 494 << " in " << nBinsCSDA << " bins" << G4endl; 499 } 495 } 500 if(nSCoffRegions>0 && isIonisation) { 496 if(nSCoffRegions>0 && isIonisation) { 501 out << " Subcutoff sampling in " << n 497 out << " Subcutoff sampling in " << nSCoffRegions 502 << " regions" << G4endl; 498 << " regions" << G4endl; 503 } 499 } 504 if(2 < verboseLevel) { 500 if(2 < verboseLevel) { 505 for(std::size_t i=0; i<7; ++i) { 501 for(std::size_t i=0; i<7; ++i) { 506 auto ta = theData->Table(i); 502 auto ta = theData->Table(i); 507 out << " " << tnames[i] << " addres 503 out << " " << tnames[i] << " address: " << ta << G4endl; 508 if(nullptr != ta) { out << *ta << G4endl 504 if(nullptr != ta) { out << *ta << G4endl; } 509 } 505 } 510 } 506 } 511 } 507 } 512 508 513 //....oooOO0OOooo........oooOO0OOooo........oo 509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 514 510 515 void G4VEnergyLossProcess::ActivateSubCutoff(c 511 void G4VEnergyLossProcess::ActivateSubCutoff(const G4Region* r) 516 { 512 { 517 if(nullptr == scoffRegions) { 513 if(nullptr == scoffRegions) { 518 scoffRegions = new std::vector<const G4Reg 514 scoffRegions = new std::vector<const G4Region*>; 519 } 515 } 520 // the region is in the list 516 // the region is in the list 521 if(!scoffRegions->empty()) { 517 if(!scoffRegions->empty()) { 522 for (auto & reg : *scoffRegions) { 518 for (auto & reg : *scoffRegions) { 523 if (reg == r) { return; } 519 if (reg == r) { return; } 524 } 520 } 525 } 521 } 526 // new region 522 // new region 527 scoffRegions->push_back(r); 523 scoffRegions->push_back(r); 528 ++nSCoffRegions; 524 ++nSCoffRegions; 529 } 525 } 530 526 531 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 532 528 533 G4bool G4VEnergyLossProcess::IsRegionForCubcut 529 G4bool G4VEnergyLossProcess::IsRegionForCubcutProcessor(const G4Track& aTrack) 534 { 530 { 535 if(0 == nSCoffRegions) { return true; } 531 if(0 == nSCoffRegions) { return true; } 536 const G4Region* r = aTrack.GetVolume()->GetL 532 const G4Region* r = aTrack.GetVolume()->GetLogicalVolume()->GetRegion(); 537 for(auto & reg : *scoffRegions) { 533 for(auto & reg : *scoffRegions) { 538 if(r == reg) { return true; } 534 if(r == reg) { return true; } 539 } 535 } 540 return false; 536 return false; 541 } 537 } 542 538 543 //....oooOO0OOooo........oooOO0OOooo........oo 539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 544 540 545 void G4VEnergyLossProcess::StartTracking(G4Tra 541 void G4VEnergyLossProcess::StartTracking(G4Track* track) 546 { 542 { 547 // reset parameters for the new track 543 // reset parameters for the new track 548 theNumberOfInteractionLengthLeft = -1.0; 544 theNumberOfInteractionLengthLeft = -1.0; 549 mfpKinEnergy = DBL_MAX; 545 mfpKinEnergy = DBL_MAX; 550 preStepLambda = 0.0; << 551 currentCouple = nullptr; 546 currentCouple = nullptr; 552 547 553 // reset ion 548 // reset ion 554 if(isIon) { 549 if(isIon) { 555 const G4double newmass = track->GetDefinit 550 const G4double newmass = track->GetDefinition()->GetPDGMass(); 556 massRatio = (nullptr == baseParticle) ? CL << 551 if(nullptr != baseParticle) { 557 : baseParticle->GetPDGMass()/newmass; << 552 massRatio = baseParticle->GetPDGMass()/newmass; 558 logMassRatio = G4Log(massRatio); << 553 logMassRatio = G4Log(massRatio); >> 554 } else { >> 555 massRatio = CLHEP::proton_mass_c2/newmass; >> 556 logMassRatio = G4Log(massRatio); >> 557 } 559 } 558 } 560 // forced biasing only for primary particles 559 // forced biasing only for primary particles 561 if(nullptr != biasManager) { 560 if(nullptr != biasManager) { 562 if(0 == track->GetParentID()) { 561 if(0 == track->GetParentID()) { 563 biasFlag = true; 562 biasFlag = true; 564 biasManager->ResetForcedInteraction(); 563 biasManager->ResetForcedInteraction(); 565 } 564 } 566 } 565 } 567 } 566 } 568 567 569 //....oooOO0OOooo........oooOO0OOooo........oo 568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 570 569 571 G4double G4VEnergyLossProcess::AlongStepGetPhy 570 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength( 572 const G4Track& tr << 571 const G4Track&,G4double,G4double,G4double&, 573 G4GPILSelection* 572 G4GPILSelection* selection) 574 { 573 { 575 G4double x = DBL_MAX; 574 G4double x = DBL_MAX; 576 *selection = aGPILSelection; 575 *selection = aGPILSelection; 577 if(isIonisation && currentModel->IsActive(pr 576 if(isIonisation && currentModel->IsActive(preStepScaledEnergy)) { 578 GetScaledRangeForScaledEnergy(preStepScale << 577 GetScaledRangeForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy); 579 x = (useCutAsFinalRange) ? std::min(finalR << 578 const G4double finR = (rndmStepFlag) ? std::min(finalRange, 580 currentCouple->GetProductionCuts()->GetP 579 currentCouple->GetProductionCuts()->GetProductionCut(1)) : finalRange; 581 x = (fRange > x) ? fRange*dRoverRange + x* << 580 x = (fRange > finR) ? 582 : fRange; << 581 fRange*dRoverRange + finR*(1.0-dRoverRange)*(2.0-finR/fRange) : fRange; 583 /* << 584 G4cout<<"AlongStepGPIL: " << GetProcessN << 585 << " fRange=" << fRange << " finR=" << finR << 586 */ << 587 } 582 } >> 583 //G4cout<<"AlongStepGPIL: " << GetProcessName()<<": e= "<<preStepKinEnergy >> 584 //<<" stepLimit= "<<x<<G4endl; 588 return x; 585 return x; 589 } 586 } 590 587 591 //....oooOO0OOooo........oooOO0OOooo........oo 588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 592 589 593 G4double G4VEnergyLossProcess::PostStepGetPhys 590 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength( 594 const G4Track& tr 591 const G4Track& track, 595 G4double previo 592 G4double previousStepSize, 596 G4ForceCondition* 593 G4ForceCondition* condition) 597 { 594 { 598 // condition is set to "Not Forced" 595 // condition is set to "Not Forced" 599 *condition = NotForced; 596 *condition = NotForced; 600 G4double x = DBL_MAX; 597 G4double x = DBL_MAX; 601 598 602 // initialisation of material, mass, charge, 599 // initialisation of material, mass, charge, model 603 // at the beginning of the step 600 // at the beginning of the step 604 DefineMaterial(track.GetMaterialCutsCouple() 601 DefineMaterial(track.GetMaterialCutsCouple()); 605 preStepKinEnergy = track.GetKineticEne 602 preStepKinEnergy = track.GetKineticEnergy(); >> 603 preStepLogKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy(); 606 preStepScaledEnergy = preStepKinEnergy*ma 604 preStepScaledEnergy = preStepKinEnergy*massRatio; >> 605 preStepLogScaledEnergy = preStepLogKinEnergy + logMassRatio; 607 SelectModel(preStepScaledEnergy); 606 SelectModel(preStepScaledEnergy); 608 607 609 if(!currentModel->IsActive(preStepScaledEner 608 if(!currentModel->IsActive(preStepScaledEnergy)) { 610 theNumberOfInteractionLengthLeft = -1.0; 609 theNumberOfInteractionLengthLeft = -1.0; 611 mfpKinEnergy = DBL_MAX; << 612 preStepLambda = 0.0; << 613 currentInteractionLength = DBL_MAX; 610 currentInteractionLength = DBL_MAX; 614 return x; << 611 return x; 615 } 612 } 616 613 617 // change effective charge of a charged part 614 // change effective charge of a charged particle on fly 618 if(isIon) { 615 if(isIon) { 619 const G4double q2 = currentModel->ChargeSq 616 const G4double q2 = currentModel->ChargeSquareRatio(track); 620 fFactor = q2*biasFactor; << 617 if(q2 != chargeSqRatio) { 621 if(baseMat) { fFactor *= (*theDensityFacto << 618 fFactor *= q2/chargeSqRatio; 622 reduceFactor = 1.0/(fFactor*massRatio); << 619 reduceFactor = 1.0/(fFactor*massRatio); >> 620 chargeSqRatio = q2; >> 621 } 623 if (lossFluctuationFlag) { 622 if (lossFluctuationFlag) { 624 auto fluc = currentModel->GetModelOfFluc 623 auto fluc = currentModel->GetModelOfFluctuations(); 625 fluc->SetParticleAndCharge(track.GetDefi 624 fluc->SetParticleAndCharge(track.GetDefinition(), q2); 626 } 625 } 627 } 626 } 628 627 629 // forced biasing only for primary particles 628 // forced biasing only for primary particles 630 if(biasManager) { 629 if(biasManager) { 631 if(0 == track.GetParentID() && biasFlag && 630 if(0 == track.GetParentID() && biasFlag && 632 biasManager->ForcedInteractionRegion((G 631 biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) { 633 return biasManager->GetStepLimit((G4int) 632 return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize); 634 } 633 } 635 } 634 } 636 635 637 ComputeLambdaForScaledEnergy(preStepScaledEn << 636 // compute mean free path 638 << 637 ComputeLambdaForScaledEnergy(preStepScaledEnergy, preStepLogScaledEnergy); >> 638 639 // zero cross section 639 // zero cross section 640 if(preStepLambda <= 0.0) { 640 if(preStepLambda <= 0.0) { 641 theNumberOfInteractionLengthLeft = -1.0; 641 theNumberOfInteractionLengthLeft = -1.0; 642 currentInteractionLength = DBL_MAX; 642 currentInteractionLength = DBL_MAX; 643 } else { 643 } else { 644 644 645 // non-zero cross section 645 // non-zero cross section 646 if (theNumberOfInteractionLengthLeft < 0.0 646 if (theNumberOfInteractionLengthLeft < 0.0) { 647 647 648 // beggining of tracking (or just after 648 // beggining of tracking (or just after DoIt of this process) 649 theNumberOfInteractionLengthLeft = -G4Lo 649 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() ); 650 theInitialNumberOfInteractionLength = th 650 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 651 651 652 } else if(currentInteractionLength < DBL_M 652 } else if(currentInteractionLength < DBL_MAX) { 653 653 654 // subtract NumberOfInteractionLengthLef 654 // subtract NumberOfInteractionLengthLeft using previous step 655 theNumberOfInteractionLengthLeft -= 655 theNumberOfInteractionLengthLeft -= 656 previousStepSize/currentInteractionLen 656 previousStepSize/currentInteractionLength; 657 657 658 theNumberOfInteractionLengthLeft = 658 theNumberOfInteractionLengthLeft = 659 std::max(theNumberOfInteractionLengthL 659 std::max(theNumberOfInteractionLengthLeft, 0.0); 660 } 660 } 661 661 662 // new mean free path and step limit 662 // new mean free path and step limit 663 currentInteractionLength = 1.0/preStepLamb 663 currentInteractionLength = 1.0/preStepLambda; 664 x = theNumberOfInteractionLengthLeft * cur 664 x = theNumberOfInteractionLengthLeft * currentInteractionLength; 665 } 665 } 666 #ifdef G4VERBOSE 666 #ifdef G4VERBOSE 667 if (verboseLevel>2) { 667 if (verboseLevel>2) { 668 G4cout << "G4VEnergyLossProcess::PostStepG 668 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength "; 669 G4cout << "[ " << GetProcessName() << "]" 669 G4cout << "[ " << GetProcessName() << "]" << G4endl; 670 G4cout << " for " << track.GetDefinition() 670 G4cout << " for " << track.GetDefinition()->GetParticleName() 671 << " in Material " << currentMate 671 << " in Material " << currentMaterial->GetName() 672 << " Ekin(MeV)= " << preStepKinEner 672 << " Ekin(MeV)= " << preStepKinEnergy/MeV 673 << " track material: " << track.Get 673 << " track material: " << track.GetMaterial()->GetName() 674 <<G4endl; 674 <<G4endl; 675 G4cout << "MeanFreePath = " << currentInte 675 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 676 << "InteractionLength= " << x/cm << 676 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl; 677 } 677 } 678 #endif 678 #endif 679 return x; 679 return x; 680 } 680 } 681 681 682 //....oooOO0OOooo........oooOO0OOooo........oo 682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 683 683 684 void 684 void 685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 685 G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e, G4double loge) 686 { 686 { 687 // cross section increased with energy 687 // cross section increased with energy 688 if(fXSType == fEmIncreasing) { 688 if(fXSType == fEmIncreasing) { 689 if(e*invLambdaFactor < mfpKinEnergy) { 689 if(e*invLambdaFactor < mfpKinEnergy) { 690 preStepLambda = GetLambdaForScaledEnergy << 690 mfpKinEnergy = e; 691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 691 preStepLambda = GetLambdaForScaledEnergy(e, loge); 692 } 692 } 693 693 694 // cross section has one peak 694 // cross section has one peak 695 } else if(fXSType == fEmOnePeak) { 695 } else if(fXSType == fEmOnePeak) { 696 const G4double epeak = (*theEnergyOfCrossS 696 const G4double epeak = (*theEnergyOfCrossSectionMax)[basedCoupleIndex]; 697 if(e <= epeak) { 697 if(e <= epeak) { 698 if(e*invLambdaFactor < mfpKinEnergy) { 698 if(e*invLambdaFactor < mfpKinEnergy) { 699 preStepLambda = GetLambdaForScaledEner << 699 mfpKinEnergy = e; 700 mfpKinEnergy = (preStepLambda > 0.0) ? << 700 preStepLambda = GetLambdaForScaledEnergy(e, loge); 701 } 701 } 702 } else if(e < mfpKinEnergy) { 702 } else if(e < mfpKinEnergy) { 703 const G4double e1 = std::max(epeak, e*la 703 const G4double e1 = std::max(epeak, e*lambdaFactor); 704 mfpKinEnergy = e1; 704 mfpKinEnergy = e1; 705 preStepLambda = GetLambdaForScaledEnergy 705 preStepLambda = GetLambdaForScaledEnergy(e1); 706 } 706 } 707 707 708 // cross section has more than one peaks 708 // cross section has more than one peaks 709 } else if(fXSType == fEmTwoPeaks) { 709 } else if(fXSType == fEmTwoPeaks) { 710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple 710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCoupleIndex]; 711 const G4double e1peak = xs->e1peak; 711 const G4double e1peak = xs->e1peak; 712 712 713 // below the 1st peak 713 // below the 1st peak 714 if(e <= e1peak) { 714 if(e <= e1peak) { 715 if(e*invLambdaFactor < mfpKinEnergy) { 715 if(e*invLambdaFactor < mfpKinEnergy) { 716 preStepLambda = GetLambdaForScaledEner << 716 mfpKinEnergy = e; 717 mfpKinEnergy = (preStepLambda > 0.0) ? << 717 preStepLambda = GetLambdaForScaledEnergy(e, loge); 718 } 718 } 719 return; 719 return; 720 } 720 } 721 const G4double e1deep = xs->e1deep; 721 const G4double e1deep = xs->e1deep; 722 // above the 1st peak, below the deep 722 // above the 1st peak, below the deep 723 if(e <= e1deep) { 723 if(e <= e1deep) { 724 if(mfpKinEnergy >= e1deep || e <= mfpKin 724 if(mfpKinEnergy >= e1deep || e <= mfpKinEnergy) { 725 const G4double e1 = std::max(e1peak, e 725 const G4double e1 = std::max(e1peak, e*lambdaFactor); 726 mfpKinEnergy = e1; 726 mfpKinEnergy = e1; 727 preStepLambda = GetLambdaForScaledEner 727 preStepLambda = GetLambdaForScaledEnergy(e1); 728 } 728 } 729 return; 729 return; 730 } 730 } 731 const G4double e2peak = xs->e2peak; 731 const G4double e2peak = xs->e2peak; 732 // above the deep, below 2nd peak 732 // above the deep, below 2nd peak 733 if(e <= e2peak) { 733 if(e <= e2peak) { 734 if(e*invLambdaFactor < mfpKinEnergy) { 734 if(e*invLambdaFactor < mfpKinEnergy) { 735 mfpKinEnergy = e; 735 mfpKinEnergy = e; 736 preStepLambda = GetLambdaForScaledEner << 736 preStepLambda = GetLambdaForScaledEnergy(e, loge); 737 } 737 } 738 return; 738 return; 739 } 739 } 740 const G4double e2deep = xs->e2deep; 740 const G4double e2deep = xs->e2deep; 741 // above the 2nd peak, below the deep 741 // above the 2nd peak, below the deep 742 if(e <= e2deep) { 742 if(e <= e2deep) { 743 if(mfpKinEnergy >= e2deep || e <= mfpKin 743 if(mfpKinEnergy >= e2deep || e <= mfpKinEnergy) { 744 const G4double e1 = std::max(e2peak, e 744 const G4double e1 = std::max(e2peak, e*lambdaFactor); 745 mfpKinEnergy = e1; 745 mfpKinEnergy = e1; 746 preStepLambda = GetLambdaForScaledEner 746 preStepLambda = GetLambdaForScaledEnergy(e1); 747 } 747 } 748 return; 748 return; 749 } 749 } 750 const G4double e3peak = xs->e3peak; 750 const G4double e3peak = xs->e3peak; 751 // above the deep, below 3d peak 751 // above the deep, below 3d peak 752 if(e <= e3peak) { 752 if(e <= e3peak) { 753 if(e*invLambdaFactor < mfpKinEnergy) { 753 if(e*invLambdaFactor < mfpKinEnergy) { 754 mfpKinEnergy = e; 754 mfpKinEnergy = e; 755 preStepLambda = GetLambdaForScaledEner << 755 preStepLambda = GetLambdaForScaledEnergy(e, loge); 756 } 756 } 757 return; 757 return; 758 } 758 } 759 // above 3d peak 759 // above 3d peak 760 if(e <= mfpKinEnergy) { 760 if(e <= mfpKinEnergy) { 761 const G4double e1 = std::max(e3peak, e*l 761 const G4double e1 = std::max(e3peak, e*lambdaFactor); 762 mfpKinEnergy = e1; 762 mfpKinEnergy = e1; 763 preStepLambda = GetLambdaForScaledEnergy 763 preStepLambda = GetLambdaForScaledEnergy(e1); 764 } 764 } 765 // integral method is not used 765 // integral method is not used 766 } else { 766 } else { 767 preStepLambda = GetLambdaForScaledEnergy(e << 767 preStepLambda = GetLambdaForScaledEnergy(e, loge); 768 } 768 } 769 } 769 } 770 770 771 //....oooOO0OOooo........oooOO0OOooo........oo 771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 772 772 773 G4VParticleChange* G4VEnergyLossProcess::Along 773 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track, 774 774 const G4Step& step) 775 { 775 { 776 fParticleChange.InitializeForAlongStep(track 776 fParticleChange.InitializeForAlongStep(track); 777 // The process has range table - calculate e 777 // The process has range table - calculate energy loss 778 if(!isIonisation || !currentModel->IsActive( 778 if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) { 779 return &fParticleChange; 779 return &fParticleChange; 780 } 780 } 781 781 >> 782 // Get the actual (true) Step length 782 G4double length = step.GetStepLength(); 783 G4double length = step.GetStepLength(); >> 784 if(length <= 0.0) { return &fParticleChange; } 783 G4double eloss = 0.0; 785 G4double eloss = 0.0; 784 786 785 /* << 787 /* 786 if(-1 < verboseLevel) { 788 if(-1 < verboseLevel) { 787 const G4ParticleDefinition* d = track.GetP 789 const G4ParticleDefinition* d = track.GetParticleDefinition(); 788 G4cout << "AlongStepDoIt for " 790 G4cout << "AlongStepDoIt for " 789 << GetProcessName() << " and partic 791 << GetProcessName() << " and particle " << d->GetParticleName() 790 << " eScaled(MeV)=" << preStepScal 792 << " eScaled(MeV)=" << preStepScaledEnergy/MeV 791 << " range(mm)=" << fRange/mm << " 793 << " range(mm)=" << fRange/mm << " s(mm)=" << length/mm 792 << " rf=" << reduceFactor << " q^ 794 << " rf=" << reduceFactor << " q^2=" << chargeSqRatio 793 << " md=" << d->GetPDGMass() << " 795 << " md=" << d->GetPDGMass() << " status=" << track.GetTrackStatus() 794 << " " << track.GetMaterial()->Get 796 << " " << track.GetMaterial()->GetName() << G4endl; 795 } 797 } 796 */ 798 */ 797 const G4DynamicParticle* dynParticle = track 799 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); 798 800 799 // define new weight for primary and seconda 801 // define new weight for primary and secondaries 800 G4double weight = fParticleChange.GetParentW 802 G4double weight = fParticleChange.GetParentWeight(); 801 if(weightFlag) { 803 if(weightFlag) { 802 weight /= biasFactor; 804 weight /= biasFactor; 803 fParticleChange.ProposeWeight(weight); 805 fParticleChange.ProposeWeight(weight); 804 } 806 } 805 807 806 // stopping, check actual range and kinetic << 808 // stopping 807 if (length >= fRange || preStepKinEnergy <= 809 if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) { 808 eloss = preStepKinEnergy; 810 eloss = preStepKinEnergy; 809 if (useDeexcitation) { 811 if (useDeexcitation) { 810 atomDeexcitation->AlongStepDeexcitation( 812 atomDeexcitation->AlongStepDeexcitation(scTracks, step, 811 813 eloss, (G4int)currentCoupleIndex); 812 if(scTracks.size() > 0) { FillSecondarie 814 if(scTracks.size() > 0) { FillSecondariesAlongStep(weight); } 813 eloss = std::max(eloss, 0.0); 815 eloss = std::max(eloss, 0.0); 814 } 816 } 815 fParticleChange.SetProposedKineticEnergy(0 817 fParticleChange.SetProposedKineticEnergy(0.0); 816 fParticleChange.ProposeLocalEnergyDeposit( 818 fParticleChange.ProposeLocalEnergyDeposit(eloss); 817 return &fParticleChange; 819 return &fParticleChange; 818 } 820 } 819 // zero step length with non-zero range << 820 if(length <= 0.0) { return &fParticleChange; << 821 << 822 // Short step 821 // Short step 823 eloss = length*GetDEDXForScaledEnergy(preSte 822 eloss = length*GetDEDXForScaledEnergy(preStepScaledEnergy, 824 LogSca << 823 preStepLogScaledEnergy); 825 /* << 824 //G4cout << "Short STEP: eloss= " << eloss << G4endl; 826 G4cout << "##### Short STEP: eloss= " << elo << 825 827 << " Escaled=" << preStepScaledEnergy << 828 << " R=" << fRange << 829 << " L=" << length << 830 << " fFactor=" << fFactor << " minE=" << mi << 831 << " idxBase=" << basedCoupleIndex << G4end << 832 */ << 833 // Long step 826 // Long step 834 if(eloss > preStepKinEnergy*linLossLimit) { 827 if(eloss > preStepKinEnergy*linLossLimit) { 835 828 836 const G4double x = (fRange - length)/reduc << 829 G4double x = (fRange - length)/reduceFactor; 837 const G4double de = preStepKinEnergy - Sca << 830 //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl; 838 if(de > 0.0) { eloss = de; } << 831 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio; >> 832 839 /* 833 /* 840 if(-1 < verboseLevel) 834 if(-1 < verboseLevel) 841 G4cout << " Long STEP: rPre(mm)=" << 835 G4cout << "Long STEP: rPre(mm)= " 842 << GetScaledRangeForScaledEnergy( 836 << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm 843 << " x(mm)=" << x/mm << 837 << " rPost(mm)= " << x/mm 844 << " eloss(MeV)=" << eloss/MeV << 838 << " ePre(MeV)= " << preStepScaledEnergy/MeV 845 << " rFactor=" << reduceFactor << 839 << " eloss(MeV)= " << eloss/MeV << " eloss0(MeV)= " 846 << " massRatio=" << massRatio << 840 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV >> 841 << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV 847 << G4endl; 842 << G4endl; 848 */ 843 */ 849 } 844 } 850 845 851 /* 846 /* 852 if(-1 < verboseLevel ) { 847 if(-1 < verboseLevel ) { 853 G4cout << "Before fluct: eloss(MeV)= " << 848 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV 854 << " e-eloss= " << preStepKinEnergy 849 << " e-eloss= " << preStepKinEnergy-eloss 855 << " step(mm)= " << length/mm << " 850 << " step(mm)= " << length/mm << " range(mm)= " << fRange/mm 856 << " fluct= " << lossFluctuationFla 851 << " fluct= " << lossFluctuationFlag << G4endl; 857 } 852 } 858 */ 853 */ 859 854 860 const G4double cut = (*theCuts)[currentCoupl 855 const G4double cut = (*theCuts)[currentCoupleIndex]; 861 G4double esec = 0.0; 856 G4double esec = 0.0; 862 857 863 // Corrections, which cannot be tabulated 858 // Corrections, which cannot be tabulated 864 if(isIon) { 859 if(isIon) { 865 currentModel->CorrectionsAlongStep(current 860 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 866 length, 861 length, eloss); 867 eloss = std::max(eloss, 0.0); 862 eloss = std::max(eloss, 0.0); 868 } 863 } 869 864 870 // Sample fluctuations if not full energy lo 865 // Sample fluctuations if not full energy loss 871 if(eloss >= preStepKinEnergy) { 866 if(eloss >= preStepKinEnergy) { 872 eloss = preStepKinEnergy; 867 eloss = preStepKinEnergy; 873 868 874 } else if (lossFluctuationFlag) { 869 } else if (lossFluctuationFlag) { 875 const G4double tmax = currentModel->MaxSec 870 const G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle); 876 const G4double tcut = std::min(cut, tmax); 871 const G4double tcut = std::min(cut, tmax); 877 G4VEmFluctuationModel* fluc = currentModel 872 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations(); 878 eloss = fluc->SampleFluctuations(currentCo 873 eloss = fluc->SampleFluctuations(currentCouple,dynParticle, 879 tcut, tma 874 tcut, tmax, length, eloss); 880 /* 875 /* 881 if(-1 < verboseLevel) 876 if(-1 < verboseLevel) 882 G4cout << "After fluct: eloss(MeV)= " << 877 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV 883 << " fluc= " << (eloss-eloss0)/Me 878 << " fluc= " << (eloss-eloss0)/MeV 884 << " ChargeSqRatio= " << chargeSq 879 << " ChargeSqRatio= " << chargeSqRatio 885 << " massRatio= " << massRatio << 880 << " massRatio= " << massRatio << " tmax= " << tmax << G4endl; 886 */ 881 */ 887 } 882 } 888 883 889 // deexcitation 884 // deexcitation 890 if (useDeexcitation) { 885 if (useDeexcitation) { 891 G4double esecfluo = preStepKinEnergy; 886 G4double esecfluo = preStepKinEnergy; 892 G4double de = esecfluo; 887 G4double de = esecfluo; 893 atomDeexcitation->AlongStepDeexcitation(sc 888 atomDeexcitation->AlongStepDeexcitation(scTracks, step, 894 de 889 de, (G4int)currentCoupleIndex); 895 890 896 // sum of de-excitation energies 891 // sum of de-excitation energies 897 esecfluo -= de; 892 esecfluo -= de; 898 893 899 // subtracted from energy loss 894 // subtracted from energy loss 900 if(eloss >= esecfluo) { 895 if(eloss >= esecfluo) { 901 esec += esecfluo; 896 esec += esecfluo; 902 eloss -= esecfluo; 897 eloss -= esecfluo; 903 } else { 898 } else { 904 esec += esecfluo; 899 esec += esecfluo; 905 eloss = 0.0; 900 eloss = 0.0; 906 } 901 } 907 } 902 } 908 if(nullptr != subcutProducer && IsRegionForC 903 if(nullptr != subcutProducer && IsRegionForCubcutProcessor(track)) { 909 subcutProducer->SampleSecondaries(step, sc 904 subcutProducer->SampleSecondaries(step, scTracks, eloss, cut); 910 } 905 } 911 // secondaries from atomic de-excitation and 906 // secondaries from atomic de-excitation and subcut 912 if(!scTracks.empty()) { FillSecondariesAlong 907 if(!scTracks.empty()) { FillSecondariesAlongStep(weight); } 913 908 914 // Energy balance 909 // Energy balance 915 G4double finalT = preStepKinEnergy - eloss - 910 G4double finalT = preStepKinEnergy - eloss - esec; 916 if (finalT <= lowestKinEnergy) { 911 if (finalT <= lowestKinEnergy) { 917 eloss += finalT; 912 eloss += finalT; 918 finalT = 0.0; 913 finalT = 0.0; 919 } else if(isIon) { 914 } else if(isIon) { 920 fParticleChange.SetProposedCharge( 915 fParticleChange.SetProposedCharge( 921 currentModel->GetParticleCharge(track.Ge 916 currentModel->GetParticleCharge(track.GetParticleDefinition(), 922 currentM 917 currentMaterial,finalT)); 923 } 918 } 924 eloss = std::max(eloss, 0.0); 919 eloss = std::max(eloss, 0.0); 925 920 926 fParticleChange.SetProposedKineticEnergy(fin 921 fParticleChange.SetProposedKineticEnergy(finalT); 927 fParticleChange.ProposeLocalEnergyDeposit(el 922 fParticleChange.ProposeLocalEnergyDeposit(eloss); 928 /* 923 /* 929 if(-1 < verboseLevel) { 924 if(-1 < verboseLevel) { 930 G4double del = finalT + eloss + esec - pre 925 G4double del = finalT + eloss + esec - preStepKinEnergy; 931 G4cout << "Final value eloss(MeV)= " << el 926 G4cout << "Final value eloss(MeV)= " << eloss/MeV 932 << " preStepKinEnergy= " << preStep 927 << " preStepKinEnergy= " << preStepKinEnergy 933 << " postStepKinEnergy= " << finalT 928 << " postStepKinEnergy= " << finalT 934 << " de(keV)= " << del/keV 929 << " de(keV)= " << del/keV 935 << " lossFlag= " << lossFluctuation 930 << " lossFlag= " << lossFluctuationFlag 936 << " status= " << track.GetTrackSt 931 << " status= " << track.GetTrackStatus() 937 << G4endl; 932 << G4endl; 938 } 933 } 939 */ 934 */ 940 return &fParticleChange; 935 return &fParticleChange; 941 } 936 } 942 937 943 //....oooOO0OOooo........oooOO0OOooo........oo 938 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 944 939 945 void G4VEnergyLossProcess::FillSecondariesAlon 940 void G4VEnergyLossProcess::FillSecondariesAlongStep(G4double wt) 946 { 941 { 947 const std::size_t n0 = scTracks.size(); 942 const std::size_t n0 = scTracks.size(); 948 G4double weight = wt; 943 G4double weight = wt; 949 // weight may be changed by biasing manager 944 // weight may be changed by biasing manager 950 if(biasManager) { 945 if(biasManager) { 951 if(biasManager->SecondaryBiasingRegion((G4 946 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) { 952 weight *= 947 weight *= 953 biasManager->ApplySecondaryBiasing(scT 948 biasManager->ApplySecondaryBiasing(scTracks, (G4int)currentCoupleIndex); 954 } 949 } 955 } 950 } 956 951 957 // fill secondaries 952 // fill secondaries 958 const std::size_t n = scTracks.size(); 953 const std::size_t n = scTracks.size(); 959 fParticleChange.SetNumberOfSecondaries((G4in 954 fParticleChange.SetNumberOfSecondaries((G4int)n); 960 955 961 for(std::size_t i=0; i<n; ++i) { 956 for(std::size_t i=0; i<n; ++i) { 962 G4Track* t = scTracks[i]; 957 G4Track* t = scTracks[i]; 963 if(nullptr != t) { 958 if(nullptr != t) { 964 t->SetWeight(weight); 959 t->SetWeight(weight); 965 pParticleChange->AddSecondary(t); 960 pParticleChange->AddSecondary(t); 966 G4int pdg = t->GetDefinition()->GetPDGEn << 961 if(i >= n0) { t->SetCreatorModelID(biasID); } 967 if (i < n0) { << 968 if (pdg == 22) { << 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 } 962 } 979 } 963 } 980 scTracks.clear(); 964 scTracks.clear(); 981 } 965 } 982 966 983 //....oooOO0OOooo........oooOO0OOooo........oo 967 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 984 968 985 G4VParticleChange* G4VEnergyLossProcess::PostS 969 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track, 986 970 const G4Step& step) 987 { 971 { 988 // clear number of interaction lengths in an 972 // clear number of interaction lengths in any case 989 theNumberOfInteractionLengthLeft = -1.0; 973 theNumberOfInteractionLengthLeft = -1.0; 990 mfpKinEnergy = DBL_MAX; 974 mfpKinEnergy = DBL_MAX; 991 975 992 fParticleChange.InitializeForPostStep(track) 976 fParticleChange.InitializeForPostStep(track); 993 const G4double finalT = track.GetKineticEner 977 const G4double finalT = track.GetKineticEnergy(); 994 978 995 const G4double postStepScaledEnergy = finalT 979 const G4double postStepScaledEnergy = finalT*massRatio; 996 SelectModel(postStepScaledEnergy); 980 SelectModel(postStepScaledEnergy); 997 981 998 if(!currentModel->IsActive(postStepScaledEne 982 if(!currentModel->IsActive(postStepScaledEnergy)) { 999 return &fParticleChange; 983 return &fParticleChange; 1000 } 984 } 1001 /* 985 /* 1002 if(1 < verboseLevel) { 986 if(1 < verboseLevel) { 1003 G4cout<<GetProcessName()<<" PostStepDoIt: 987 G4cout<<GetProcessName()<<" PostStepDoIt: E(MeV)= "<< finalT/MeV<< G4endl; 1004 } 988 } 1005 */ 989 */ 1006 // forced process - should happen only once 990 // forced process - should happen only once per track 1007 if(biasFlag) { 991 if(biasFlag) { 1008 if(biasManager->ForcedInteractionRegion(( 992 if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) { 1009 biasFlag = false; 993 biasFlag = false; 1010 } 994 } 1011 } 995 } 1012 const G4DynamicParticle* dp = track.GetDyna 996 const G4DynamicParticle* dp = track.GetDynamicParticle(); 1013 997 1014 // Integral approach 998 // Integral approach 1015 if (fXSType != fEmNoIntegral) { 999 if (fXSType != fEmNoIntegral) { 1016 const G4double logFinalT = dp->GetLogKine 1000 const G4double logFinalT = dp->GetLogKineticEnergy(); 1017 G4double lx = GetLambdaForScaledEnergy(po 1001 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy, 1018 lo 1002 logFinalT + logMassRatio); 1019 lx = std::max(lx, 0.0); 1003 lx = std::max(lx, 0.0); 1020 1004 1021 // if both lg and lx are zero then no int 1005 // if both lg and lx are zero then no interaction 1022 if(preStepLambda*G4UniformRand() >= lx) { 1006 if(preStepLambda*G4UniformRand() >= lx) { 1023 return &fParticleChange; 1007 return &fParticleChange; 1024 } 1008 } 1025 } 1009 } 1026 1010 1027 // define new weight for primary and second 1011 // define new weight for primary and secondaries 1028 G4double weight = fParticleChange.GetParent 1012 G4double weight = fParticleChange.GetParentWeight(); 1029 if(weightFlag) { 1013 if(weightFlag) { 1030 weight /= biasFactor; 1014 weight /= biasFactor; 1031 fParticleChange.ProposeWeight(weight); 1015 fParticleChange.ProposeWeight(weight); 1032 } 1016 } 1033 1017 1034 const G4double tcut = (*theCuts)[currentCou 1018 const G4double tcut = (*theCuts)[currentCoupleIndex]; 1035 1019 1036 // sample secondaries 1020 // sample secondaries 1037 secParticles.clear(); 1021 secParticles.clear(); 1038 currentModel->SampleSecondaries(&secParticl 1022 currentModel->SampleSecondaries(&secParticles, currentCouple, dp, tcut); 1039 1023 1040 const G4int num0 = (G4int)secParticles.size 1024 const G4int num0 = (G4int)secParticles.size(); 1041 1025 1042 // bremsstrahlung splitting or Russian roul 1026 // bremsstrahlung splitting or Russian roulette 1043 if(biasManager) { 1027 if(biasManager) { 1044 if(biasManager->SecondaryBiasingRegion((G 1028 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) { 1045 G4double eloss = 0.0; 1029 G4double eloss = 0.0; 1046 weight *= biasManager->ApplySecondaryBi 1030 weight *= biasManager->ApplySecondaryBiasing( 1047 secPart 1031 secParticles, 1048 track, 1032 track, currentModel, 1049 &fParti 1033 &fParticleChange, eloss, 1050 (G4int) 1034 (G4int)currentCoupleIndex, tcut, 1051 step.Ge 1035 step.GetPostStepPoint()->GetSafety()); 1052 if(eloss > 0.0) { 1036 if(eloss > 0.0) { 1053 eloss += fParticleChange.GetLocalEner 1037 eloss += fParticleChange.GetLocalEnergyDeposit(); 1054 fParticleChange.ProposeLocalEnergyDep 1038 fParticleChange.ProposeLocalEnergyDeposit(eloss); 1055 } 1039 } 1056 } 1040 } 1057 } 1041 } 1058 1042 1059 // save secondaries 1043 // save secondaries 1060 const G4int num = (G4int)secParticles.size( 1044 const G4int num = (G4int)secParticles.size(); 1061 if(num > 0) { 1045 if(num > 0) { 1062 1046 1063 fParticleChange.SetNumberOfSecondaries(nu 1047 fParticleChange.SetNumberOfSecondaries(num); 1064 G4double time = track.GetGlobalTime(); 1048 G4double time = track.GetGlobalTime(); 1065 1049 1066 G4int n1(0), n2(0); 1050 G4int n1(0), n2(0); 1067 if(num0 > mainSecondaries) { 1051 if(num0 > mainSecondaries) { 1068 currentModel->FillNumberOfSecondaries(n 1052 currentModel->FillNumberOfSecondaries(n1, n2); 1069 } 1053 } 1070 1054 1071 for (G4int i=0; i<num; ++i) { 1055 for (G4int i=0; i<num; ++i) { 1072 if(nullptr != secParticles[i]) { 1056 if(nullptr != secParticles[i]) { 1073 G4Track* t = new G4Track(secParticles 1057 G4Track* t = new G4Track(secParticles[i], time, track.GetPosition()); 1074 t->SetTouchableHandle(track.GetToucha 1058 t->SetTouchableHandle(track.GetTouchableHandle()); 1075 if (biasManager) { 1059 if (biasManager) { 1076 t->SetWeight(weight * biasManager-> 1060 t->SetWeight(weight * biasManager->GetWeight(i)); 1077 } else { 1061 } else { 1078 t->SetWeight(weight); 1062 t->SetWeight(weight); 1079 } 1063 } 1080 if(i < num0) { 1064 if(i < num0) { 1081 t->SetCreatorModelID(secID); 1065 t->SetCreatorModelID(secID); 1082 } else if(i < num0 + n1) { 1066 } else if(i < num0 + n1) { 1083 t->SetCreatorModelID(tripletID); 1067 t->SetCreatorModelID(tripletID); 1084 } else { 1068 } else { 1085 t->SetCreatorModelID(biasID); 1069 t->SetCreatorModelID(biasID); 1086 } 1070 } 1087 1071 1088 //G4cout << "Secondary(post step) has 1072 //G4cout << "Secondary(post step) has weight " << t->GetWeight() 1089 // << ", kenergy " << t->GetKin 1073 // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" 1090 // << " time= " << time/ns << " 1074 // << " time= " << time/ns << " ns " << G4endl; 1091 pParticleChange->AddSecondary(t); 1075 pParticleChange->AddSecondary(t); 1092 } 1076 } 1093 } 1077 } 1094 } 1078 } 1095 1079 1096 if(0.0 == fParticleChange.GetProposedKineti 1080 if(0.0 == fParticleChange.GetProposedKineticEnergy() && 1097 fAlive == fParticleChange.GetTrackStatus 1081 fAlive == fParticleChange.GetTrackStatus()) { 1098 if(particle->GetProcessManager()->GetAtRe 1082 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0) 1099 { fParticleChange.ProposeTrackStatus 1083 { fParticleChange.ProposeTrackStatus(fStopButAlive); } 1100 else { fParticleChange.ProposeTrackStatus 1084 else { fParticleChange.ProposeTrackStatus(fStopAndKill); } 1101 } 1085 } 1102 1086 1103 /* 1087 /* 1104 if(-1 < verboseLevel) { 1088 if(-1 < verboseLevel) { 1105 G4cout << "::PostStepDoIt: Sample seconda 1089 G4cout << "::PostStepDoIt: Sample secondary; Efin= " 1106 << fParticleChange.GetProposedKineticEner 1090 << fParticleChange.GetProposedKineticEnergy()/MeV 1107 << " MeV; model= (" << currentMode 1091 << " MeV; model= (" << currentModel->LowEnergyLimit() 1108 << ", " << currentModel->HighEner 1092 << ", " << currentModel->HighEnergyLimit() << ")" 1109 << " preStepLambda= " << preStepL 1093 << " preStepLambda= " << preStepLambda 1110 << " dir= " << track.GetMomentumD 1094 << " dir= " << track.GetMomentumDirection() 1111 << " status= " << track.GetTrackS 1095 << " status= " << track.GetTrackStatus() 1112 << G4endl; 1096 << G4endl; 1113 } 1097 } 1114 */ 1098 */ 1115 return &fParticleChange; 1099 return &fParticleChange; 1116 } 1100 } 1117 1101 1118 //....oooOO0OOooo........oooOO0OOooo........o 1102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1119 1103 1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl 1104 G4bool G4VEnergyLossProcess::StorePhysicsTable( 1121 const G4ParticleDefinition* part, cons 1105 const G4ParticleDefinition* part, const G4String& dir, G4bool ascii) 1122 { 1106 { 1123 if (!isMaster || nullptr != baseParticle || 1107 if (!isMaster || nullptr != baseParticle || part != particle ) return true; 1124 for(std::size_t i=0; i<7; ++i) { 1108 for(std::size_t i=0; i<7; ++i) { 1125 // ionisation table only for ionisation p << 1109 if(nullptr != theData->Table(i)) { 1126 if (nullptr == theData->Table(i) || (!isI << 1110 if(1 < verboseLevel) { 1127 continue; << 1111 G4cout << "G4VEnergyLossProcess::StorePhysicsTable i=" << i 1128 } << 1112 << " " << particle->GetParticleName() 1129 if (-1 < verboseLevel) { << 1113 << " " << GetProcessName() 1130 G4cout << "G4VEnergyLossProcess::StoreP << 1114 << " " << tnames[i] << " " << theData->Table(i) << G4endl; 1131 << " " << particle->GetParticleName() << 1115 } 1132 << " " << GetProcessName() << 1116 if(!G4EmTableUtil::StoreTable(this, part, theData->Table(i), 1133 << " " << tnames[i] << " " << theDat << 1117 dir, tnames[i], verboseLevel, ascii)) { 1134 } << 1118 return false; 1135 if (!G4EmTableUtil::StoreTable(this, part << 1119 } 1136 dir, tnames[i], verboseLevel, asci << 1137 return false; << 1138 } 1120 } 1139 } 1121 } 1140 return true; 1122 return true; 1141 } 1123 } 1142 1124 1143 //....oooOO0OOooo........oooOO0OOooo........o 1125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1144 1126 1145 G4bool 1127 G4bool 1146 G4VEnergyLossProcess::RetrievePhysicsTable(co 1128 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 1147 co 1129 const G4String& dir, G4bool ascii) 1148 { 1130 { 1149 if (!isMaster || nullptr != baseParticle || 1131 if (!isMaster || nullptr != baseParticle || part != particle ) return true; 1150 for(std::size_t i=0; i<7; ++i) { 1132 for(std::size_t i=0; i<7; ++i) { 1151 // ionisation table only for ionisation p << 1152 if (!isIonisation && 1 == i) { continue; << 1153 if(!G4EmTableUtil::RetrieveTable(this, pa 1133 if(!G4EmTableUtil::RetrieveTable(this, part, theData->Table(i), dir, tnames[i], 1154 verboseL 1134 verboseLevel, ascii, spline)) { 1155 return false; 1135 return false; 1156 } 1136 } 1157 } 1137 } 1158 return true; 1138 return true; 1159 } 1139 } 1160 1140 1161 //....oooOO0OOooo........oooOO0OOooo........o 1141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1162 1142 1163 G4double G4VEnergyLossProcess::GetDEDXDispers 1143 G4double G4VEnergyLossProcess::GetDEDXDispersion( 1164 const G4Mat 1144 const G4MaterialCutsCouple *couple, 1165 const G4Dyn 1145 const G4DynamicParticle* dp, 1166 G4dou 1146 G4double length) 1167 { 1147 { 1168 DefineMaterial(couple); 1148 DefineMaterial(couple); 1169 G4double ekin = dp->GetKineticEnergy(); 1149 G4double ekin = dp->GetKineticEnergy(); 1170 SelectModel(ekin*massRatio); 1150 SelectModel(ekin*massRatio); 1171 G4double tmax = currentModel->MaxSecondaryK 1151 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp); 1172 G4double tcut = std::min(tmax,(*theCuts)[cu 1152 G4double tcut = std::min(tmax,(*theCuts)[currentCoupleIndex]); 1173 G4double d = 0.0; 1153 G4double d = 0.0; 1174 G4VEmFluctuationModel* fm = currentModel->G 1154 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations(); 1175 if(nullptr != fm) { d = fm->Dispersion(curr 1155 if(nullptr != fm) { d = fm->Dispersion(currentMaterial,dp,tcut,tmax,length); } 1176 return d; 1156 return d; 1177 } 1157 } 1178 1158 1179 //....oooOO0OOooo........oooOO0OOooo........o 1159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1180 1160 1181 G4double 1161 G4double 1182 G4VEnergyLossProcess::CrossSectionPerVolume(G 1162 G4VEnergyLossProcess::CrossSectionPerVolume(G4double kineticEnergy, 1183 c 1163 const G4MaterialCutsCouple* couple, 1184 G 1164 G4double logKineticEnergy) 1185 { 1165 { 1186 // Cross section per volume is calculated 1166 // Cross section per volume is calculated 1187 DefineMaterial(couple); 1167 DefineMaterial(couple); 1188 G4double cross = 0.0; 1168 G4double cross = 0.0; 1189 if (nullptr != theLambdaTable) { 1169 if (nullptr != theLambdaTable) { 1190 cross = GetLambdaForScaledEnergy(kineticE 1170 cross = GetLambdaForScaledEnergy(kineticEnergy * massRatio, 1191 logKinet 1171 logKineticEnergy + logMassRatio); 1192 } else { 1172 } else { 1193 SelectModel(kineticEnergy*massRatio); 1173 SelectModel(kineticEnergy*massRatio); 1194 cross = (!baseMat) ? biasFactor : biasFac 1174 cross = (!baseMat) ? biasFactor : biasFactor*(*theDensityFactor)[currentCoupleIndex]; 1195 cross *= (currentModel->CrossSectionPerVo 1175 cross *= (currentModel->CrossSectionPerVolume(currentMaterial, particle, kineticEnergy, 1196 1176 (*theCuts)[currentCoupleIndex])); 1197 } 1177 } 1198 return std::max(cross, 0.0); 1178 return std::max(cross, 0.0); 1199 } 1179 } 1200 1180 1201 //....oooOO0OOooo........oooOO0OOooo........o 1181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1202 1182 1203 G4double G4VEnergyLossProcess::MeanFreePath(c 1183 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track) 1204 { 1184 { 1205 DefineMaterial(track.GetMaterialCutsCouple( 1185 DefineMaterial(track.GetMaterialCutsCouple()); 1206 const G4double kinEnergy = track.GetKine 1186 const G4double kinEnergy = track.GetKineticEnergy(); 1207 const G4double logKinEnergy = track.GetDyna 1187 const G4double logKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy(); 1208 const G4double cs = GetLambdaForScaledEnerg 1188 const G4double cs = GetLambdaForScaledEnergy(kinEnergy * massRatio, 1209 1189 logKinEnergy + logMassRatio); 1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1190 return (0.0 < cs) ? 1.0/cs : DBL_MAX; 1211 } 1191 } 1212 1192 1213 //....oooOO0OOooo........oooOO0OOooo........o 1193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1214 1194 1215 G4double G4VEnergyLossProcess::ContinuousStep 1195 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, 1216 1196 G4double x, G4double y, 1217 1197 G4double& z) 1218 { 1198 { 1219 return AlongStepGetPhysicalInteractionLengt 1199 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &aGPILSelection); 1220 } 1200 } 1221 1201 1222 //....oooOO0OOooo........oooOO0OOooo........o 1202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1223 1203 1224 G4double G4VEnergyLossProcess::GetMeanFreePat 1204 G4double G4VEnergyLossProcess::GetMeanFreePath( 1225 const G4Track& t 1205 const G4Track& track, 1226 G4double, 1206 G4double, 1227 G4ForceCondition 1207 G4ForceCondition* condition) 1228 1208 1229 { 1209 { 1230 *condition = NotForced; 1210 *condition = NotForced; 1231 return MeanFreePath(track); 1211 return MeanFreePath(track); 1232 } 1212 } 1233 1213 1234 //....oooOO0OOooo........oooOO0OOooo........o 1214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1235 1215 1236 G4double G4VEnergyLossProcess::GetContinuousS 1216 G4double G4VEnergyLossProcess::GetContinuousStepLimit( 1237 const G4Track&, 1217 const G4Track&, 1238 G4double, G4double, G4double& 1218 G4double, G4double, G4double&) 1239 { 1219 { 1240 return DBL_MAX; 1220 return DBL_MAX; 1241 } 1221 } 1242 1222 1243 //....oooOO0OOooo........oooOO0OOooo........o 1223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1244 1224 1245 G4PhysicsVector* 1225 G4PhysicsVector* 1246 G4VEnergyLossProcess::LambdaPhysicsVector(con 1226 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple, 1247 G4d 1227 G4double) 1248 { 1228 { 1249 DefineMaterial(couple); 1229 DefineMaterial(couple); 1250 G4PhysicsVector* v = (*theLambdaTable)[base 1230 G4PhysicsVector* v = (*theLambdaTable)[basedCoupleIndex]; 1251 return new G4PhysicsVector(*v); 1231 return new G4PhysicsVector(*v); 1252 } 1232 } 1253 1233 1254 //....oooOO0OOooo........oooOO0OOooo........o 1234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1255 1235 1256 void 1236 void 1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT 1237 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType) 1258 { 1238 { 1259 if(1 < verboseLevel) { 1239 if(1 < verboseLevel) { 1260 G4cout << "### Set DEDX table " << p << " 1240 G4cout << "### Set DEDX table " << p << " " << theDEDXTable 1261 << " " << theDEDXunRestrictedTable << 1241 << " " << theDEDXunRestrictedTable << " " << theIonisationTable 1262 << " for " << particle->GetParticl 1242 << " for " << particle->GetParticleName() 1263 << " and process " << GetProcessNa 1243 << " and process " << GetProcessName() 1264 << " type=" << tType << " isIonisation:" 1244 << " type=" << tType << " isIonisation:" << isIonisation << G4endl; 1265 } 1245 } 1266 if(fTotal == tType) { 1246 if(fTotal == tType) { 1267 theDEDXunRestrictedTable = p; 1247 theDEDXunRestrictedTable = p; 1268 } else if(fRestricted == tType) { 1248 } else if(fRestricted == tType) { 1269 theDEDXTable = p; 1249 theDEDXTable = p; 1270 if(isMaster && nullptr == baseParticle) { 1250 if(isMaster && nullptr == baseParticle) { 1271 theData->UpdateTable(theDEDXTable, 0); 1251 theData->UpdateTable(theDEDXTable, 0); 1272 } 1252 } 1273 } else if(fIsIonisation == tType) { 1253 } else if(fIsIonisation == tType) { 1274 theIonisationTable = p; 1254 theIonisationTable = p; 1275 if(isMaster && nullptr == baseParticle) { 1255 if(isMaster && nullptr == baseParticle) { 1276 theData->UpdateTable(theIonisationTable 1256 theData->UpdateTable(theIonisationTable, 1); 1277 } 1257 } 1278 } 1258 } 1279 } 1259 } 1280 1260 1281 //....oooOO0OOooo........oooOO0OOooo........o 1261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1282 1262 1283 void G4VEnergyLossProcess::SetCSDARangeTable( 1263 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p) 1284 { 1264 { 1285 theCSDARangeTable = p; 1265 theCSDARangeTable = p; 1286 } 1266 } 1287 1267 1288 //....oooOO0OOooo........oooOO0OOooo........o 1268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1289 1269 1290 void G4VEnergyLossProcess::SetRangeTableForLo 1270 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p) 1291 { 1271 { 1292 theRangeTableForLoss = p; 1272 theRangeTableForLoss = p; 1293 } 1273 } 1294 1274 1295 //....oooOO0OOooo........oooOO0OOooo........o 1275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1296 1276 1297 void G4VEnergyLossProcess::SetInverseRangeTab 1277 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p) 1298 { 1278 { 1299 theInverseRangeTable = p; 1279 theInverseRangeTable = p; 1300 } 1280 } 1301 1281 1302 //....oooOO0OOooo........oooOO0OOooo........o 1282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1303 1283 1304 void G4VEnergyLossProcess::SetLambdaTable(G4P 1284 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p) 1305 { 1285 { 1306 if(1 < verboseLevel) { 1286 if(1 < verboseLevel) { 1307 G4cout << "### Set Lambda table " << p << 1287 G4cout << "### Set Lambda table " << p << " " << theLambdaTable 1308 << " for " << particle->GetParticl 1288 << " for " << particle->GetParticleName() 1309 << " and process " << GetProcessNa 1289 << " and process " << GetProcessName() << G4endl; 1310 } 1290 } 1311 theLambdaTable = p; 1291 theLambdaTable = p; 1312 tablesAreBuilt = true; 1292 tablesAreBuilt = true; 1313 1293 1314 if(isMaster && nullptr != p) { 1294 if(isMaster && nullptr != p) { 1315 delete theEnergyOfCrossSectionMax; 1295 delete theEnergyOfCrossSectionMax; 1316 theEnergyOfCrossSectionMax = nullptr; 1296 theEnergyOfCrossSectionMax = nullptr; 1317 if(fEmTwoPeaks == fXSType) { 1297 if(fEmTwoPeaks == fXSType) { 1318 if(nullptr != fXSpeaks) { 1298 if(nullptr != fXSpeaks) { 1319 for(auto & ptr : *fXSpeaks) { delete ptr; } 1299 for(auto & ptr : *fXSpeaks) { delete ptr; } 1320 delete fXSpeaks; 1300 delete fXSpeaks; 1321 } 1301 } 1322 G4LossTableBuilder* bld = lManager->Get 1302 G4LossTableBuilder* bld = lManager->GetTableBuilder(); 1323 fXSpeaks = G4EmUtility::FillPeaksStruct 1303 fXSpeaks = G4EmUtility::FillPeaksStructure(p, bld); 1324 if(nullptr == fXSpeaks) { fXSType = fEm 1304 if(nullptr == fXSpeaks) { fXSType = fEmOnePeak; } 1325 } 1305 } 1326 if(fXSType == fEmOnePeak) { 1306 if(fXSType == fEmOnePeak) { 1327 theEnergyOfCrossSectionMax = G4EmUtilit 1307 theEnergyOfCrossSectionMax = G4EmUtility::FindCrossSectionMax(p); 1328 if(nullptr == theEnergyOfCrossSectionMa 1308 if(nullptr == theEnergyOfCrossSectionMax) { fXSType = fEmIncreasing; } 1329 } 1309 } 1330 } 1310 } 1331 } 1311 } 1332 1312 1333 //....oooOO0OOooo........oooOO0OOooo........o 1313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1334 1314 1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe 1315 void G4VEnergyLossProcess::SetEnergyOfCrossSectionMax(std::vector<G4double>* p) 1336 { 1316 { 1337 theEnergyOfCrossSectionMax = p; 1317 theEnergyOfCrossSectionMax = p; 1338 } 1318 } 1339 1319 1340 //....oooOO0OOooo........oooOO0OOooo........o 1320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1341 1321 1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: 1322 void G4VEnergyLossProcess::SetTwoPeaksXS(std::vector<G4TwoPeaksXS*>* ptr) 1343 { 1323 { 1344 fXSpeaks = ptr; 1324 fXSpeaks = ptr; 1345 } 1325 } 1346 1326 1347 //....oooOO0OOooo........oooOO0OOooo........o 1327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1348 1328 1349 const G4Element* G4VEnergyLossProcess::GetCur 1329 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const 1350 { 1330 { 1351 return (nullptr != currentModel) 1331 return (nullptr != currentModel) 1352 ? currentModel->GetCurrentElement(current 1332 ? currentModel->GetCurrentElement(currentMaterial) : nullptr; 1353 } 1333 } 1354 1334 1355 //....oooOO0OOooo........oooOO0OOooo........o 1335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1356 1336 1357 void G4VEnergyLossProcess::SetCrossSectionBia 1337 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f, 1358 1338 G4bool flag) 1359 { 1339 { 1360 if(f > 0.0) { 1340 if(f > 0.0) { 1361 biasFactor = f; 1341 biasFactor = f; 1362 weightFlag = flag; 1342 weightFlag = flag; 1363 if(1 < verboseLevel) { 1343 if(1 < verboseLevel) { 1364 G4cout << "### SetCrossSectionBiasingFa 1344 G4cout << "### SetCrossSectionBiasingFactor: for " 1365 << " process " << GetProcessName 1345 << " process " << GetProcessName() 1366 << " biasFactor= " << f << " wei 1346 << " biasFactor= " << f << " weightFlag= " << flag 1367 << G4endl; 1347 << G4endl; 1368 } 1348 } 1369 } 1349 } 1370 } 1350 } 1371 1351 1372 //....oooOO0OOooo........oooOO0OOooo........o 1352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1373 1353 1374 void G4VEnergyLossProcess::ActivateForcedInte 1354 void G4VEnergyLossProcess::ActivateForcedInteraction(G4double length, 1375 1355 const G4String& region, 1376 1356 G4bool flag) 1377 { 1357 { 1378 if(nullptr == biasManager) { biasManager = 1358 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); } 1379 if(1 < verboseLevel) { 1359 if(1 < verboseLevel) { 1380 G4cout << "### ActivateForcedInteraction: 1360 G4cout << "### ActivateForcedInteraction: for " 1381 << " process " << GetProcessName() 1361 << " process " << GetProcessName() 1382 << " length(mm)= " << length/mm 1362 << " length(mm)= " << length/mm 1383 << " in G4Region <" << region 1363 << " in G4Region <" << region 1384 << "> weightFlag= " << flag 1364 << "> weightFlag= " << flag 1385 << G4endl; 1365 << G4endl; 1386 } 1366 } 1387 weightFlag = flag; 1367 weightFlag = flag; 1388 biasManager->ActivateForcedInteraction(leng 1368 biasManager->ActivateForcedInteraction(length, region); 1389 } 1369 } 1390 1370 1391 //....oooOO0OOooo........oooOO0OOooo........o 1371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1392 1372 1393 void 1373 void 1394 G4VEnergyLossProcess::ActivateSecondaryBiasin 1374 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region, 1395 1375 G4double factor, 1396 1376 G4double energyLimit) 1397 { 1377 { 1398 if (0.0 <= factor) { 1378 if (0.0 <= factor) { 1399 // Range cut can be applied only for e- 1379 // Range cut can be applied only for e- 1400 if(0.0 == factor && secondaryParticle != 1380 if(0.0 == factor && secondaryParticle != G4Electron::Electron()) 1401 { return; } 1381 { return; } 1402 1382 1403 if(nullptr == biasManager) { biasManager 1383 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); } 1404 biasManager->ActivateSecondaryBiasing(reg 1384 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit); 1405 if(1 < verboseLevel) { 1385 if(1 < verboseLevel) { 1406 G4cout << "### ActivateSecondaryBiasing 1386 G4cout << "### ActivateSecondaryBiasing: for " 1407 << " process " << GetProcessName 1387 << " process " << GetProcessName() 1408 << " factor= " << factor 1388 << " factor= " << factor 1409 << " in G4Region <" << region 1389 << " in G4Region <" << region 1410 << "> energyLimit(MeV)= " << ene 1390 << "> energyLimit(MeV)= " << energyLimit/MeV 1411 << G4endl; 1391 << G4endl; 1412 } 1392 } 1413 } 1393 } 1414 } 1394 } 1415 1395 1416 //....oooOO0OOooo........oooOO0OOooo........o 1396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1417 1397 1418 void G4VEnergyLossProcess::SetIonisation(G4bo 1398 void G4VEnergyLossProcess::SetIonisation(G4bool val) 1419 { 1399 { 1420 isIonisation = val; 1400 isIonisation = val; 1421 aGPILSelection = (val) ? CandidateForSelect 1401 aGPILSelection = (val) ? CandidateForSelection : NotCandidateForSelection; 1422 } 1402 } 1423 1403 1424 //....oooOO0OOooo........oooOO0OOooo........o 1404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1425 1405 1426 void G4VEnergyLossProcess::SetLinearLossLimi 1406 void G4VEnergyLossProcess::SetLinearLossLimit(G4double val) 1427 { 1407 { 1428 if(0.0 < val && val < 1.0) { 1408 if(0.0 < val && val < 1.0) { 1429 linLossLimit = val; 1409 linLossLimit = val; 1430 actLinLossLimit = true; 1410 actLinLossLimit = true; 1431 } else { PrintWarning("SetLinearLossLimit", 1411 } else { PrintWarning("SetLinearLossLimit", val); } 1432 } 1412 } 1433 1413 1434 //....oooOO0OOooo........oooOO0OOooo........o 1414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1435 1415 1436 void G4VEnergyLossProcess::SetStepFunction(G4 1416 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2) 1437 { 1417 { 1438 if(0.0 < v1 && 0.0 < v2) { 1418 if(0.0 < v1 && 0.0 < v2) { 1439 dRoverRange = std::min(1.0, v1); 1419 dRoverRange = std::min(1.0, v1); 1440 finalRange = std::min(v2, 1.e+50); 1420 finalRange = std::min(v2, 1.e+50); 1441 } else { 1421 } else { 1442 PrintWarning("SetStepFunctionV1", v1); 1422 PrintWarning("SetStepFunctionV1", v1); 1443 PrintWarning("SetStepFunctionV2", v2); 1423 PrintWarning("SetStepFunctionV2", v2); 1444 } 1424 } 1445 } 1425 } 1446 1426 1447 //....oooOO0OOooo........oooOO0OOooo........o 1427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1448 1428 1449 void G4VEnergyLossProcess::SetLowestEnergyLim 1429 void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val) 1450 { 1430 { 1451 if(1.e-18 < val && val < 1.e+50) { lowestKi 1431 if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; } 1452 else { PrintWarning("SetLowestEnergyLimit", 1432 else { PrintWarning("SetLowestEnergyLimit", val); } 1453 } 1433 } 1454 1434 1455 //....oooOO0OOooo........oooOO0OOooo........o 1435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1456 1436 1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i 1437 void G4VEnergyLossProcess::SetDEDXBinning(G4int n) 1458 { 1438 { 1459 if(2 < n && n < 1000000000) { 1439 if(2 < n && n < 1000000000) { 1460 nBins = n; 1440 nBins = n; 1461 actBinning = true; 1441 actBinning = true; 1462 } else { 1442 } else { 1463 G4double e = (G4double)n; 1443 G4double e = (G4double)n; 1464 PrintWarning("SetDEDXBinning", e); 1444 PrintWarning("SetDEDXBinning", e); 1465 } 1445 } 1466 } 1446 } 1467 1447 1468 //....oooOO0OOooo........oooOO0OOooo........o 1448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1469 1449 1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 1450 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e) 1471 { 1451 { 1472 if(1.e-18 < e && e < maxKinEnergy) { 1452 if(1.e-18 < e && e < maxKinEnergy) { 1473 minKinEnergy = e; 1453 minKinEnergy = e; 1474 actMinKinEnergy = true; 1454 actMinKinEnergy = true; 1475 } else { PrintWarning("SetMinKinEnergy", e) 1455 } else { PrintWarning("SetMinKinEnergy", e); } 1476 } 1456 } 1477 1457 1478 //....oooOO0OOooo........oooOO0OOooo........o 1458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1479 1459 1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 1460 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e) 1481 { 1461 { 1482 if(minKinEnergy < e && e < 1.e+50) { 1462 if(minKinEnergy < e && e < 1.e+50) { 1483 maxKinEnergy = e; 1463 maxKinEnergy = e; 1484 actMaxKinEnergy = true; 1464 actMaxKinEnergy = true; 1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS 1465 if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; } 1486 } else { PrintWarning("SetMaxKinEnergy", e) 1466 } else { PrintWarning("SetMaxKinEnergy", e); } 1487 } 1467 } 1488 1468 1489 //....oooOO0OOooo........oooOO0OOooo........o 1469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1490 1470 1491 void G4VEnergyLossProcess::PrintWarning(const 1471 void G4VEnergyLossProcess::PrintWarning(const G4String& tit, G4double val) const 1492 { 1472 { 1493 G4String ss = "G4VEnergyLossProcess::" + ti 1473 G4String ss = "G4VEnergyLossProcess::" + tit; 1494 G4ExceptionDescription ed; 1474 G4ExceptionDescription ed; 1495 ed << "Parameter is out of range: " << val 1475 ed << "Parameter is out of range: " << val 1496 << " it will have no effect!\n" << " Pr 1476 << " it will have no effect!\n" << " Process " 1497 << GetProcessName() << " nbins= " << nB 1477 << GetProcessName() << " nbins= " << nBins 1498 << " Emin(keV)= " << minKinEnergy/keV 1478 << " Emin(keV)= " << minKinEnergy/keV 1499 << " Emax(GeV)= " << maxKinEnergy/GeV; 1479 << " Emax(GeV)= " << maxKinEnergy/GeV; 1500 G4Exception(ss, "em0044", JustWarning, ed); 1480 G4Exception(ss, "em0044", JustWarning, ed); 1501 } 1481 } 1502 1482 1503 //....oooOO0OOooo........oooOO0OOooo........o 1483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1504 1484 1505 void G4VEnergyLossProcess::ProcessDescription 1485 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const 1506 { 1486 { 1507 if(nullptr != particle) { StreamInfo(out, * 1487 if(nullptr != particle) { StreamInfo(out, *particle, true); } 1508 } 1488 } 1509 1489 1510 //....oooOO0OOooo........oooOO0OOooo........o 1490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1511 1491