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