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