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