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