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