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