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