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