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