| 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: G4VEmProcess.cc,v 1.88 2010-08-17 17:36:59 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
>> 28 //
26 // ------------------------------------------- 29 // -------------------------------------------------------------------
27 // 30 //
28 // GEANT4 Class file 31 // GEANT4 Class file
29 // 32 //
30 // 33 //
31 // File name: G4VEmProcess 34 // File name: G4VEmProcess
32 // 35 //
33 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 // 37 //
35 // Creation date: 01.10.2003 38 // Creation date: 01.10.2003
36 // 39 //
37 // Modifications: by V.Ivanchenko << 40 // Modifications:
>> 41 // 30-06-04 make it to be pure discrete process (V.Ivanchenko)
>> 42 // 30-09-08 optimise integral option (V.Ivanchenko)
>> 43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
>> 44 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
>> 45 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
>> 46 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
>> 47 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
>> 48 // 25-07-05 Add protection: integral mode only for charged particles (VI)
>> 49 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
>> 50 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
>> 51 // 12-09-06 add SetModel() (mma)
>> 52 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
>> 53 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
>> 54 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
>> 55 // 17-02-10 Added pointer currentParticle (VI)
38 // 56 //
39 // Class Description: based class for discrete << 57 // Class Description:
40 // 58 //
41 59
42 // ------------------------------------------- 60 // -------------------------------------------------------------------
43 // 61 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 64
47 #include "G4VEmProcess.hh" 65 #include "G4VEmProcess.hh"
48 #include "G4PhysicalConstants.hh" <<
49 #include "G4SystemOfUnits.hh" <<
50 #include "G4ProcessManager.hh" <<
51 #include "G4LossTableManager.hh" 66 #include "G4LossTableManager.hh"
52 #include "G4LossTableBuilder.hh" <<
53 #include "G4Step.hh" 67 #include "G4Step.hh"
54 #include "G4ParticleDefinition.hh" 68 #include "G4ParticleDefinition.hh"
55 #include "G4VEmModel.hh" 69 #include "G4VEmModel.hh"
56 #include "G4DataVector.hh" 70 #include "G4DataVector.hh"
57 #include "G4PhysicsTable.hh" 71 #include "G4PhysicsTable.hh"
58 #include "G4EmDataHandler.hh" << 72 #include "G4PhysicsVector.hh"
59 #include "G4PhysicsLogVector.hh" 73 #include "G4PhysicsLogVector.hh"
60 #include "G4VParticleChange.hh" 74 #include "G4VParticleChange.hh"
61 #include "G4ProductionCutsTable.hh" 75 #include "G4ProductionCutsTable.hh"
62 #include "G4Region.hh" 76 #include "G4Region.hh"
>> 77 #include "G4RegionStore.hh"
63 #include "G4Gamma.hh" 78 #include "G4Gamma.hh"
64 #include "G4Electron.hh" 79 #include "G4Electron.hh"
65 #include "G4Positron.hh" 80 #include "G4Positron.hh"
66 #include "G4PhysicsTableHelper.hh" 81 #include "G4PhysicsTableHelper.hh"
67 #include "G4EmBiasingManager.hh" << 82 #include "G4EmConfigurator.hh"
68 #include "G4EmParameters.hh" <<
69 #include "G4EmProcessSubType.hh" <<
70 #include "G4EmTableUtil.hh" <<
71 #include "G4EmUtility.hh" <<
72 #include "G4DNAModelSubType.hh" <<
73 #include "G4GenericIon.hh" <<
74 #include "G4Log.hh" <<
75 #include <iostream> <<
76 83
77 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 85
79 G4VEmProcess::G4VEmProcess(const G4String& nam 86 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
80 G4VDiscreteProcess(name, type) << 87 G4VDiscreteProcess(name, type),
>> 88 secondaryParticle(0),
>> 89 buildLambdaTable(true),
>> 90 theLambdaTable(0),
>> 91 theEnergyOfCrossSectionMax(0),
>> 92 theCrossSectionMax(0),
>> 93 integral(false),
>> 94 applyCuts(false),
>> 95 startFromNull(true),
>> 96 useDeexcitation(false),
>> 97 nDERegions(0),
>> 98 idxDERegions(0),
>> 99 currentModel(0),
>> 100 particle(0),
>> 101 currentParticle(0),
>> 102 currentCouple(0)
81 { 103 {
82 theParameters = G4EmParameters::Instance(); <<
83 SetVerboseLevel(1); 104 SetVerboseLevel(1);
84 105
85 // Size of tables << 106 // Size of tables assuming spline
86 minKinEnergy = 0.1*CLHEP::keV; << 107 minKinEnergy = 0.1*keV;
87 maxKinEnergy = 100.0*CLHEP::TeV; << 108 maxKinEnergy = 10.0*TeV;
>> 109 nLambdaBins = 77;
88 110
89 // default lambda factor 111 // default lambda factor
90 invLambdaFactor = 1.0/lambdaFactor; << 112 lambdaFactor = 0.8;
>> 113
>> 114 // default limit on polar angle
>> 115 polarAngleLimit = 0.0;
91 116
92 // particle types 117 // particle types
93 theGamma = G4Gamma::Gamma(); << 118 theGamma = G4Gamma::Gamma();
94 theElectron = G4Electron::Electron(); << 119 theElectron = G4Electron::Electron();
95 thePositron = G4Positron::Positron(); << 120 thePositron = G4Positron::Positron();
96 121
97 pParticleChange = &fParticleChange; 122 pParticleChange = &fParticleChange;
98 fParticleChange.SetSecondaryWeightByProcess( <<
99 secParticles.reserve(5); 123 secParticles.reserve(5);
100 124
101 modelManager = new G4EmModelManager(); 125 modelManager = new G4EmModelManager();
102 lManager = G4LossTableManager::Instance(); << 126 (G4LossTableManager::Instance())->Register(this);
103 lManager->Register(this); <<
104 isTheMaster = lManager->IsMaster(); <<
105 G4LossTableBuilder* bld = lManager->GetTable <<
106 theDensityFactor = bld->GetDensityFactors(); <<
107 theDensityIdx = bld->GetCoupleIndexes(); <<
108 } 127 }
109 128
110 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 130
112 G4VEmProcess::~G4VEmProcess() 131 G4VEmProcess::~G4VEmProcess()
113 { 132 {
114 if(isTheMaster) { << 133 if(1 < verboseLevel)
115 delete theData; << 134 G4cout << "G4VEmProcess destruct " << GetProcessName()
116 delete theEnergyOfCrossSectionMax; << 135 << G4endl;
>> 136 Clear();
>> 137 if(theLambdaTable) {
>> 138 theLambdaTable->clearAndDestroy();
>> 139 delete theLambdaTable;
117 } 140 }
118 delete modelManager; 141 delete modelManager;
119 delete biasManager; << 142 (G4LossTableManager::Instance())->DeRegister(this);
120 lManager->DeRegister(this); <<
121 } 143 }
122 144
123 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 146
125 void G4VEmProcess::AddEmModel(G4int order, G4V << 147 void G4VEmProcess::Clear()
126 const G4Region* <<
127 { 148 {
128 if(nullptr == ptr) { return; } << 149 delete [] theEnergyOfCrossSectionMax;
129 G4VEmFluctuationModel* fm = nullptr; << 150 delete [] theCrossSectionMax;
130 modelManager->AddEmModel(order, ptr, fm, reg << 151 delete [] idxDERegions;
131 ptr->SetParticleChange(pParticleChange); << 152 theEnergyOfCrossSectionMax = 0;
>> 153 theCrossSectionMax = 0;
>> 154 idxDERegions = 0;
>> 155 currentCouple = 0;
>> 156 preStepLambda = 0.0;
>> 157 mfpKinEnergy = DBL_MAX;
>> 158 deRegions.clear();
>> 159 nDERegions = 0;
132 } 160 }
133 161
134 //....oooOO0OOooo........oooOO0OOooo........oo 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 163
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 164 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
>> 165 const G4Region* region)
137 { 166 {
138 if(nullptr == ptr) { return; } << 167 G4VEmFluctuationModel* fm = 0;
139 if(!emModels.empty()) { << 168 modelManager->AddEmModel(order, p, fm, region);
140 for(auto & em : emModels) { if(em == ptr) << 169 if(p) { p->SetParticleChange(pParticleChange); }
141 } <<
142 emModels.push_back(ptr); <<
143 } 170 }
144 171
145 //....oooOO0OOooo........oooOO0OOooo........oo 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 173
147 void G4VEmProcess::PreparePhysicsTable(const G << 174 void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
148 { 175 {
149 if(nullptr == particle) { SetParticle(&part) << 176 G4int n = emModels.size();
>> 177 if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
>> 178 emModels[index] = p;
>> 179 }
>> 180
>> 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 182
151 if(part.GetParticleType() == "nucleus" && << 183 G4VEmModel* G4VEmProcess::Model(G4int index)
152 part.GetParticleSubType() == "generic") { << 184 {
>> 185 G4VEmModel* p = 0;
>> 186 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
>> 187 return p;
>> 188 }
153 189
154 G4String pname = part.GetParticleName(); << 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 if(pname != "deuteron" && pname != "triton <<
156 pname != "He3" && pname != "alpha" && p <<
157 pname != "helium" && pname != "hydrogen <<
158 191
159 particle = G4GenericIon::GenericIon(); << 192 void G4VEmProcess::UpdateEmModel(const G4String& nam,
160 isIon = true; << 193 G4double emin, G4double emax)
161 } << 194 {
162 } << 195 modelManager->UpdateEmModel(nam, emin, emax);
163 if(particle != &part) { return; } << 196 }
164 197
165 lManager->PreparePhysicsTable(&part, this); << 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 199
167 // for new run << 200 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
168 currentCouple = nullptr; << 201 {
169 preStepLambda = 0.0; << 202 return modelManager->GetModel(idx, ver);
170 fLambdaEnergy = 0.0; << 203 }
171 204
172 InitialiseProcess(particle); << 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 206
174 G4LossTableBuilder* bld = lManager->GetTable << 207 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
175 const G4ProductionCutsTable* theCoupleTable= << 208 {
176 G4ProductionCutsTable::GetProductionCutsTa << 209 if(!particle) { SetParticle(&part); }
177 theCutsGamma = theCoupleTable->GetEnergyC << 210 if(1 < verboseLevel) {
178 theCutsElectron = theCoupleTable->GetEnergyC << 211 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
179 theCutsPositron = theCoupleTable->GetEnergyC << 212 << GetProcessName()
180 << 213 << " and particle " << part.GetParticleName()
181 // initialisation of the process << 214 << " local particle " << particle->GetParticleName()
182 if(!actMinKinEnergy) { minKinEnergy = thePar << 215 << G4endl;
183 if(!actMaxKinEnergy) { maxKinEnergy = thePar << 216 }
184 <<
185 applyCuts = theParameters->ApplyCuts() <<
186 lambdaFactor = theParameters->LambdaFacto <<
187 invLambdaFactor = 1.0/lambdaFactor; <<
188 theParameters->DefineRegParamForEM(this); <<
189 <<
190 // integral option may be disabled <<
191 if(!theParameters->Integral()) { fXSType = f <<
192 <<
193 // prepare tables <<
194 if(isTheMaster) { <<
195 if(nullptr == theData) { theData = new G4E <<
196 217
197 if(buildLambdaTable) { << 218 (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
198 theLambdaTable = theData->MakeTable(0); << 219
199 bld->InitialiseBaseMaterials(theLambdaTa << 220 if(particle == &part) {
>> 221 Clear();
>> 222 InitialiseProcess(particle);
>> 223
>> 224 // initialisation of models
>> 225 G4int nmod = modelManager->NumberOfModels();
>> 226 for(G4int i=0; i<nmod; ++i) {
>> 227 G4VEmModel* mod = modelManager->GetModel(i);
>> 228 mod->SetPolarAngleLimit(polarAngleLimit);
>> 229 if(mod->HighEnergyLimit() > maxKinEnergy) {
>> 230 mod->SetHighEnergyLimit(maxKinEnergy);
>> 231 }
>> 232 }
>> 233
>> 234 theCuts = modelManager->Initialise(particle,secondaryParticle,
>> 235 2.,verboseLevel);
>> 236 const G4ProductionCutsTable* theCoupleTable=
>> 237 G4ProductionCutsTable::GetProductionCutsTable();
>> 238 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
>> 239 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
>> 240 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
>> 241
>> 242 // prepare tables
>> 243 if(buildLambdaTable){
>> 244 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
>> 245 }
>> 246 }
>> 247 // Deexcitation
>> 248 if (nDERegions>0) {
>> 249
>> 250 const G4ProductionCutsTable* theCoupleTable=
>> 251 G4ProductionCutsTable::GetProductionCutsTable();
>> 252 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 253
>> 254 idxDERegions = new G4bool[numOfCouples];
>> 255
>> 256 for (size_t j=0; j<numOfCouples; ++j) {
>> 257
>> 258 const G4MaterialCutsCouple* couple =
>> 259 theCoupleTable->GetMaterialCutsCouple(j);
>> 260 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
>> 261 G4bool reg = false;
>> 262 for(G4int i=0; i<nDERegions; ++i) {
>> 263 if(deRegions[i]) {
>> 264 if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
>> 265 }
>> 266 }
>> 267 idxDERegions[j] = reg;
200 } 268 }
201 // high energy table << 269 }
202 if(minKinEnergyPrim < maxKinEnergy) { << 270 if (1 < verboseLevel && nDERegions>0) {
203 theLambdaTablePrim = theData->MakeTable( << 271 G4cout << " Deexcitation is activated for regions: " << G4endl;
204 bld->InitialiseBaseMaterials(theLambdaTa << 272 for (G4int i=0; i<nDERegions; ++i) {
205 } << 273 const G4Region* r = deRegions[i];
206 } << 274 G4cout << " " << r->GetName() << G4endl;
207 // models << 275 }
208 baseMat = bld->GetBaseMaterialFlag(); << 276 }
209 numberOfModels = modelManager->NumberOfModel <<
210 currentModel = modelManager->GetModel(0); <<
211 if(nullptr != lManager->AtomDeexcitation()) <<
212 modelManager->SetFluoFlag(true); <<
213 } <<
214 // forced biasing <<
215 if(nullptr != biasManager) { <<
216 biasManager->Initialise(part, GetProcessNa <<
217 biasFlag = false; <<
218 } <<
219 <<
220 theCuts = <<
221 G4EmTableUtil::PrepareEmProcess(this, part <<
222 modelManag <<
223 secID, tri <<
224 verboseLev <<
225 } 277 }
226 278
227 //....oooOO0OOooo........oooOO0OOooo........oo 279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 280
229 void G4VEmProcess::BuildPhysicsTable(const G4P 281 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 { 282 {
231 if(nullptr == masterProc) { << 283 G4String partname = part.GetParticleName();
232 if(isTheMaster) { masterProc = this; } << 284 if(1 < verboseLevel) {
233 else { masterProc = static_cast<const G4VE << 285 G4cout << "G4VEmProcess::BuildPhysicsTable() for "
234 } << 286 << GetProcessName()
235 G4int nModels = modelManager->NumberOfModels << 287 << " and particle " << partname
236 G4bool isLocked = theParameters->IsPrintLock << 288 << " buildLambdaTable= " << buildLambdaTable
237 G4bool toBuild = (buildLambdaTable || minKin << 289 << G4endl;
238 << 290 }
239 G4EmTableUtil::BuildEmProcess(this, masterPr << 291
240 nModels, verbo << 292 (G4LossTableManager::Instance())->BuildPhysicsTable(particle);
241 isLocked, toBu << 293 if(buildLambdaTable) {
>> 294 BuildLambdaTable();
>> 295 FindLambdaMax();
>> 296 }
>> 297
>> 298 // reduce printout for nuclear stopping
>> 299 G4bool gproc = true;
>> 300 G4int st = GetProcessSubType();
>> 301 if(st == fCoulombScattering && part.GetParticleType() == "nucleus" &&
>> 302 partname != "GenericIon" && partname != "alpha") { gproc = false; }
>> 303
>> 304 if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); }
>> 305
>> 306 if(1 < verboseLevel) {
>> 307 G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
>> 308 << GetProcessName()
>> 309 << " and particle " << partname
>> 310 << G4endl;
>> 311 }
242 } 312 }
243 313
244 //....oooOO0OOooo........oooOO0OOooo........oo 314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 315
246 void G4VEmProcess::BuildLambdaTable() 316 void G4VEmProcess::BuildLambdaTable()
247 { 317 {
248 G4double scale = theParameters->MaxKinEnergy << 318 if(1 < verboseLevel) {
249 G4int nbin = << 319 G4cout << "G4EmProcess::BuildLambdaTable() for process "
250 theParameters->NumberOfBinsPerDecade()*G4l << 320 << GetProcessName() << " and particle "
251 if(actBinning) { nbin = std::max(nbin, nLamb << 321 << particle->GetParticleName()
252 scale = nbin/G4Log(scale); << 322 << G4endl;
253 << 323 }
254 G4LossTableBuilder* bld = lManager->GetTable <<
255 G4EmTableUtil::BuildLambdaTable(this, partic <<
256 bld, theLamb <<
257 minKinEnergy <<
258 maxKinEnergy <<
259 startFromNul <<
260 } <<
261 324
262 //....oooOO0OOooo........oooOO0OOooo........oo << 325 // Access to materials
>> 326 const G4ProductionCutsTable* theCoupleTable=
>> 327 G4ProductionCutsTable::GetProductionCutsTable();
>> 328 size_t numOfCouples = theCoupleTable->GetTableSize();
263 329
264 void G4VEmProcess::StreamInfo(std::ostream& ou << 330 G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
265 const G4ParticleDefinition& << 331
266 { << 332 G4PhysicsLogVector* aVector = 0;
267 G4String indent = (rst ? " " : ""); << 333 G4PhysicsLogVector* bVector = 0;
268 out << std::setprecision(6); << 334
269 out << G4endl << indent << GetProcessName() << 335 for(size_t i=0; i<numOfCouples; ++i) {
270 if (!rst) { << 336
271 out << " for " << part.GetParticleName(); << 337 if (theLambdaTable->GetFlag(i)) {
272 } << 338
273 if(fXSType != fEmNoIntegral) { out << " XSt << 339 // create physics vector and fill it
274 if(applyCuts) { out << " applyCuts:1 "; } << 340 const G4MaterialCutsCouple* couple =
275 G4int subtype = GetProcessSubType(); << 341 theCoupleTable->GetMaterialCutsCouple(i);
276 out << " SubType=" << subtype; << 342 if(!bVector) {
277 if (subtype == fAnnihilation) { << 343 aVector =
278 G4int mod = theParameters->PositronAtRestM << 344 static_cast<G4PhysicsLogVector*>(LambdaPhysicsVector(couple));
279 const G4String namp[2] = {"Simple", "Allis << 345 bVector = aVector;
280 out << " AtRestModel:" << namp[mod]; << 346 } else {
281 } << 347 aVector = new G4PhysicsLogVector(*bVector);
282 if(biasFactor != 1.0) { out << " BiasingFac <<
283 out << " BuildTable=" << buildLambdaTable << <<
284 if(buildLambdaTable) { <<
285 if(particle == &part) { <<
286 for(auto & v : *theLambdaTable) { <<
287 if(nullptr != v) { <<
288 out << " Lambda table from "; <<
289 G4double emin = v->Energy(0); <<
290 G4double emax = v->GetMaxEnergy(); <<
291 G4int nbin = G4int(v->GetVectorLengt <<
292 if(emin > minKinEnergy) { out << "th <<
293 else { out << G4BestUnit(emin,"Energ <<
294 out << " to " <<
295 << G4BestUnit(emax,"Energy") <<
296 << ", " << G4lrint(nbin/std::log <<
297 << " bins/decade, spline: " <<
298 << splineFlag << G4endl; <<
299 break; <<
300 } <<
301 } <<
302 } else { <<
303 out << " Used Lambda table of " <<
304 << particle->GetParticleName() << G4endl <<
305 } <<
306 } <<
307 if(minKinEnergyPrim < maxKinEnergy) { <<
308 if(particle == &part) { <<
309 for(auto & v : *theLambdaTablePrim) { <<
310 if(nullptr != v) { <<
311 out << " LambdaPrime table from <<
312 << G4BestUnit(v->Energy(0),"Ener <<
313 << " to " <<
314 << G4BestUnit(v->GetMaxEnergy(), <<
315 << " in " << v->GetVectorLength( <<
316 << " bins " << G4endl; <<
317 break; <<
318 } <<
319 } 348 }
320 } else { << 349 // G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
321 out << " Used LambdaPrime table of << 350 aVector->SetSpline(splineFlag);
322 << particle->GetParticleName() << 351 modelManager->FillLambdaVector(aVector, couple, startFromNull);
>> 352 if(splineFlag) aVector->FillSecondDerivatives();
>> 353 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
323 } 354 }
324 } 355 }
325 StreamProcessInfo(out); <<
326 modelManager->DumpModelList(out, verboseLeve <<
327 356
328 if(verboseLevel > 2 && buildLambdaTable) { << 357 if(1 < verboseLevel) {
329 out << " LambdaTable address= " << th << 358 G4cout << "Lambda table is built for "
330 if(theLambdaTable && particle == &part) { << 359 << particle->GetParticleName()
331 out << (*theLambdaTable) << G4endl; << 360 << G4endl;
332 } <<
333 } 361 }
334 } 362 }
335 363
336 //....oooOO0OOooo........oooOO0OOooo........oo 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337 365
338 void G4VEmProcess::StartTracking(G4Track* trac << 366 void G4VEmProcess::PrintInfoDefinition()
339 { 367 {
340 // reset parameters for the new track << 368 if(verboseLevel > 0) {
341 currentParticle = track->GetParticleDefiniti << 369 G4cout << G4endl << GetProcessName() << ": for "
342 theNumberOfInteractionLengthLeft = -1.0; << 370 << particle->GetParticleName();
343 mfpKinEnergy = DBL_MAX; << 371 if(integral) G4cout << ", integral: 1 ";
344 preStepLambda = 0.0; << 372 if(applyCuts) G4cout << ", applyCuts: 1 ";
345 << 373 G4cout << " SubType= " << GetProcessSubType() << G4endl;
346 if(isIon) { massRatio = proton_mass_c2/curre << 374 if(buildLambdaTable) {
347 << 375 G4cout << " Lambda tables from "
348 // forced biasing only for primary particles << 376 << G4BestUnit(minKinEnergy,"Energy")
349 if(biasManager) { << 377 << " to "
350 if(0 == track->GetParentID()) { << 378 << G4BestUnit(maxKinEnergy,"Energy")
351 // primary particle << 379 << " in " << nLambdaBins << " bins, spline: "
352 biasFlag = true; << 380 << (G4LossTableManager::Instance())->SplineFlag()
353 biasManager->ResetForcedInteraction(); << 381 << G4endl;
354 } 382 }
>> 383 PrintInfo();
>> 384 modelManager->DumpModelList(verboseLevel);
>> 385 }
>> 386
>> 387 if(verboseLevel > 2 && buildLambdaTable) {
>> 388 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
>> 389 if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
355 } 390 }
356 } 391 }
357 392
358 //....oooOO0OOooo........oooOO0OOooo........oo 393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 394
360 G4double G4VEmProcess::PostStepGetPhysicalInte 395 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& tr 396 const G4Track& track,
362 G4double previo 397 G4double previousStepSize,
363 G4ForceCondition* 398 G4ForceCondition* condition)
364 { 399 {
>> 400 // condition is set to "Not Forced"
365 *condition = NotForced; 401 *condition = NotForced;
366 G4double x = DBL_MAX; 402 G4double x = DBL_MAX;
367 << 403 if(previousStepSize <= DBL_MIN) { theNumberOfInteractionLengthLeft = -1.0; }
368 DefineMaterial(track.GetMaterialCutsCouple() << 404 InitialiseStep(track);
369 preStepKinEnergy = track.GetKineticEnergy(); << 405 if(!currentModel->IsActive(preStepKinEnergy)) { return x; }
370 const G4double scaledEnergy = preStepKinEner << 406
371 SelectModel(scaledEnergy, currentCoupleIndex << 407 if(preStepKinEnergy < mfpKinEnergy) {
372 /* << 408 if (integral) ComputeIntegralLambda(preStepKinEnergy);
373 G4cout << "PostStepGetPhysicalInteractionLen << 409 else preStepLambda = GetCurrentLambda(preStepKinEnergy);
374 << " couple: " << currentCouple << G << 410 if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
375 */ <<
376 if(!currentModel->IsActive(scaledEnergy)) { <<
377 theNumberOfInteractionLengthLeft = -1.0; <<
378 currentInteractionLength = DBL_MAX; <<
379 mfpKinEnergy = DBL_MAX; <<
380 preStepLambda = 0.0; <<
381 return x; <<
382 } <<
383 <<
384 // forced biasing only for primary particles <<
385 if(biasManager) { <<
386 if(0 == track.GetParentID()) { <<
387 if(biasFlag && <<
388 biasManager->ForcedInteractionRegion( <<
389 return biasManager->GetStepLimit((G4in <<
390 } <<
391 } <<
392 } 411 }
393 412
394 // compute mean free path << 413 // non-zero cross section
395 << 414 if(preStepLambda > DBL_MIN) {
396 ComputeIntegralLambda(preStepKinEnergy, trac <<
397 <<
398 // zero cross section <<
399 if(preStepLambda <= 0.0) { <<
400 theNumberOfInteractionLengthLeft = -1.0; <<
401 currentInteractionLength = DBL_MAX; <<
402 <<
403 } else { <<
404 <<
405 // non-zero cross section <<
406 if (theNumberOfInteractionLengthLeft < 0.0 415 if (theNumberOfInteractionLengthLeft < 0.0) {
407 <<
408 // beggining of tracking (or just after 416 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Lo << 417 ResetNumberOfInteractionLengthLeft();
410 theInitialNumberOfInteractionLength = th << 418 } else if(currentInteractionLength < DBL_MAX) {
411 << 419 // subtract NumberOfInteractionLengthLeft
412 } else { << 420 SubtractNumberOfInteractionLengthLeft(previousStepSize);
413 << 421 if(theNumberOfInteractionLengthLeft < 0.)
414 theNumberOfInteractionLengthLeft -= << 422 theNumberOfInteractionLengthLeft = perMillion;
415 previousStepSize/currentInteractionLen <<
416 theNumberOfInteractionLengthLeft = <<
417 std::max(theNumberOfInteractionLengthL <<
418 } 423 }
419 424
420 // new mean free path and step limit for t << 425 // get mean free path and step limit
421 currentInteractionLength = 1.0/preStepLamb 426 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * cur 427 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 } << 428 #ifdef G4VERBOSE
424 return x; << 429 if (verboseLevel>2){
425 } << 430 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
426 << 431 G4cout << "[ " << GetProcessName() << "]" << G4endl;
427 //....oooOO0OOooo........oooOO0OOooo........oo << 432 G4cout << " for " << currentParticle->GetParticleName()
>> 433 << " in Material " << currentMaterial->GetName()
>> 434 << " Ekin(MeV)= " << preStepKinEnergy/MeV
>> 435 <<G4endl;
>> 436 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
>> 437 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
>> 438 }
>> 439 #endif
428 440
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 441 // zero cross section case
430 { << 442 } else {
431 if (fXSType == fEmNoIntegral) { << 443 if(theNumberOfInteractionLengthLeft > DBL_MIN &&
432 preStepLambda = GetCurrentLambda(e, LogEki << 444 currentInteractionLength < DBL_MAX) {
433 445
434 } else if (fXSType == fEmIncreasing) { << 446 // subtract NumberOfInteractionLengthLeft
435 if(e*invLambdaFactor < mfpKinEnergy) { << 447 SubtractNumberOfInteractionLengthLeft(previousStepSize);
436 preStepLambda = GetCurrentLambda(e, LogE << 448 if(theNumberOfInteractionLengthLeft < 0.)
437 mfpKinEnergy = (preStepLambda > 0.0) ? e << 449 theNumberOfInteractionLengthLeft = perMillion;
438 } <<
439 <<
440 } else if(fXSType == fEmDecreasing) { <<
441 if(e < mfpKinEnergy) { <<
442 const G4double e1 = e*lambdaFactor; <<
443 preStepLambda = GetCurrentLambda(e1); <<
444 mfpKinEnergy = e1; <<
445 } <<
446 <<
447 } else if(fXSType == fEmOnePeak) { <<
448 const G4double epeak = (*theEnergyOfCrossS <<
449 if(e <= epeak) { <<
450 if(e*invLambdaFactor < mfpKinEnergy) { <<
451 preStepLambda = GetCurrentLambda(e, Lo <<
452 mfpKinEnergy = (preStepLambda > 0.0) ? <<
453 } <<
454 } else if(e < mfpKinEnergy) { <<
455 const G4double e1 = std::max(epeak, e*la <<
456 preStepLambda = GetCurrentLambda(e1); <<
457 mfpKinEnergy = e1; <<
458 } 450 }
459 } else { << 451 currentInteractionLength = DBL_MAX;
460 preStepLambda = GetCurrentLambda(e, LogEki <<
461 } 452 }
>> 453 return x;
462 } 454 }
463 455
464 //....oooOO0OOooo........oooOO0OOooo........oo 456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 457
466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 458 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 << 459 const G4Step&)
468 { 460 {
469 // clear number of interaction lengths in an <<
470 theNumberOfInteractionLengthLeft = -1.0; <<
471 mfpKinEnergy = DBL_MAX; <<
472 <<
473 fParticleChange.InitializeForPostStep(track) 461 fParticleChange.InitializeForPostStep(track);
474 462
475 // Do not make anything if particle is stopp 463 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt! 464 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) 465 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478 466
479 const G4double finalT = track.GetKineticEner << 467 G4double finalT = track.GetKineticEnergy();
480 <<
481 // forced process - should happen only once <<
482 if(biasFlag) { <<
483 if(biasManager->ForcedInteractionRegion((G <<
484 biasFlag = false; <<
485 } <<
486 } <<
487 <<
488 // check active and select model <<
489 const G4double scaledEnergy = finalT*massRat <<
490 SelectModel(scaledEnergy, currentCoupleIndex <<
491 if(!currentModel->IsActive(scaledEnergy)) { <<
492 468
493 // Integral approach 469 // Integral approach
494 if (fXSType != fEmNoIntegral) { << 470 if (integral) {
495 const G4double logFinalT = << 471 G4double lx = GetLambda(finalT, currentCouple);
496 track.GetDynamicParticle()->GetLogKineti << 472 if(preStepLambda<lx && 1 < verboseLevel) {
497 const G4double lx = std::max(GetCurrentLam << 473 G4cout << "WARING: for " << currentParticle->GetParticleName()
498 #ifdef G4VERBOSE << 474 << " and " << GetProcessName()
499 if(preStepLambda < lx && 1 < verboseLevel) << 475 << " E(MeV)= " << finalT/MeV
500 G4cout << "WARNING: for " << currentPart << 476 << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
501 << " and " << GetProcessName() << << 477 << G4endl;
502 << " preLambda= " << preStepLambd <<
503 << " < " << lx << " (postLambda) <<
504 } 478 }
505 #endif << 479
506 // if false interaction then use new cross << 480 if(preStepLambda*G4UniformRand() > lx) {
507 // if both values are zero - no interactio << 481 ClearNumberOfInteractionLengthLeft();
508 if(preStepLambda*G4UniformRand() >= lx) { <<
509 return &fParticleChange; 482 return &fParticleChange;
510 } 483 }
511 } 484 }
512 485
513 // define new weight for primary and seconda << 486 SelectModel(finalT, currentCoupleIndex);
514 G4double weight = fParticleChange.GetParentW << 487 if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
515 if(weightFlag) { << 488 if(useDeexcitation) {
516 weight /= biasFactor; << 489 currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]);
517 fParticleChange.ProposeWeight(weight); <<
518 } 490 }
519 << 491 /*
520 #ifdef G4VERBOSE << 492 if(0 < verboseLevel) {
521 if(1 < verboseLevel) { <<
522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 493 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV 494 << finalT/MeV
524 << " MeV; model= (" << currentModel 495 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnerg 496 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl; 497 << G4endl;
527 } 498 }
528 #endif << 499 */
529 500
>> 501
530 // sample secondaries 502 // sample secondaries
531 secParticles.clear(); 503 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticle 504 currentModel->SampleSecondaries(&secParticles,
533 currentCoupl << 505 currentCouple,
534 track.GetDyn << 506 track.GetDynamicParticle(),
535 (*theCuts)[c << 507 (*theCuts)[currentCoupleIndex]);
536 <<
537 G4int num0 = (G4int)secParticles.size(); <<
538 <<
539 // splitting or Russian roulette <<
540 if(biasManager) { <<
541 if(biasManager->SecondaryBiasingRegion((G4 <<
542 G4double eloss = 0.0; <<
543 weight *= biasManager->ApplySecondaryBia <<
544 secParticles, track, currentModel, &fP <<
545 (G4int)currentCoupleIndex, (*theCuts)[ <<
546 step.GetPostStepPoint()->GetSafety()); <<
547 if(eloss > 0.0) { <<
548 eloss += fParticleChange.GetLocalEnerg <<
549 fParticleChange.ProposeLocalEnergyDepo <<
550 } <<
551 } <<
552 } <<
553 508
554 // save secondaries 509 // save secondaries
555 G4int num = (G4int)secParticles.size(); << 510 G4int num = secParticles.size();
556 if(num > 0) { 511 if(num > 0) {
557 512
558 fParticleChange.SetNumberOfSecondaries(num 513 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEn 514 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime(); <<
561 <<
562 G4int n1(0), n2(0); <<
563 if(num0 > mainSecondaries) { <<
564 currentModel->FillNumberOfSecondaries(n1 <<
565 } <<
566 515
567 for (G4int i=0; i<num; ++i) { 516 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i]; 517 G4DynamicParticle* dp = secParticles[i];
569 if (nullptr != dp) { << 518 const G4ParticleDefinition* p = dp->GetParticleDefinition();
570 const G4ParticleDefinition* p = dp->Ge << 519 G4double e = dp->GetKineticEnergy();
571 G4double e = dp->GetKineticEnergy(); << 520 G4bool good = true;
572 G4bool good = true; << 521 if(applyCuts) {
573 if(applyCuts) { << 522 if (p == theGamma) {
574 if (p == theGamma) { << 523 if (e < (*theCutsGamma)[currentCoupleIndex]) good = false;
575 if (e < (*theCutsGamma)[currentCou << 524
576 << 525 } else if (p == theElectron) {
577 } else if (p == theElectron) { << 526 if (e < (*theCutsElectron)[currentCoupleIndex]) good = false;
578 if (e < (*theCutsElectron)[current << 527
579 << 528 } else if (p == thePositron) {
580 } else if (p == thePositron) { << 529 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
581 if (electron_mass_c2 < (*theCutsGa << 530 e < (*theCutsPositron)[currentCoupleIndex]) {
582 e < (*theCutsPositron)[current << 531 good = false;
583 good = false; << 532 e += 2.0*electron_mass_c2;
584 e += 2.0*electron_mass_c2; << 533 }
585 } << 534 }
586 } << 535 if(!good) {
587 // added secondary if it is good << 536 delete dp;
588 } << 537 edep += e;
589 if (good) { << 538 }
590 G4Track* t = new G4Track(dp, time, t << 539 }
591 t->SetTouchableHandle(track.GetTouch << 540 if (good) fParticleChange.AddSecondary(dp);
592 if (biasManager) { << 541 }
593 t->SetWeight(weight * biasManager- <<
594 } else { <<
595 t->SetWeight(weight); <<
596 } <<
597 pParticleChange->AddSecondary(t); <<
598 <<
599 // define type of secondary <<
600 if(i < mainSecondaries) { <<
601 t->SetCreatorModelID(secID); <<
602 if(GetProcessSubType() == fCompton <<
603 t->SetCreatorModelID(_ComptonGam <<
604 } <<
605 } else if(i < mainSecondaries + n1) <<
606 t->SetCreatorModelID(tripletID); <<
607 } else if(i < mainSecondaries + n1 + <<
608 t->SetCreatorModelID(_IonRecoil); <<
609 } else { <<
610 if(i < num0) { <<
611 if(p == theGamma) { <<
612 t->SetCreatorModelID(fluoID); <<
613 } else { <<
614 t->SetCreatorModelID(augerID); <<
615 } <<
616 } else { <<
617 t->SetCreatorModelID(biasID); <<
618 } <<
619 } <<
620 /* <<
621 G4cout << "Secondary(post step) has <<
622 << ", Ekin= " << t->GetKineti <<
623 << GetProcessName() << " fluo <<
624 << " augerID= " << augerID << <<
625 */ <<
626 } else { <<
627 delete dp; <<
628 edep += e; <<
629 } <<
630 } <<
631 } <<
632 fParticleChange.ProposeLocalEnergyDeposit( 542 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 } 543 }
634 544
635 if(0.0 == fParticleChange.GetProposedKinetic << 545 ClearNumberOfInteractionLengthLeft();
636 fAlive == fParticleChange.GetTrackStatus( <<
637 if(particle->GetProcessManager()->GetAtRes <<
638 { fParticleChange.ProposeTrackStatus( <<
639 else { fParticleChange.ProposeTrackStatus( <<
640 } <<
641 <<
642 return &fParticleChange; 546 return &fParticleChange;
643 } 547 }
644 548
645 //....oooOO0OOooo........oooOO0OOooo........oo 549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 550
647 G4bool G4VEmProcess::StorePhysicsTable(const G 551 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G << 552 const G4String& directory,
649 G4bool << 553 G4bool ascii)
650 { 554 {
651 if(!isTheMaster || part != particle) { retur << 555 G4bool yes = true;
652 if(G4EmTableUtil::StoreTable(this, part, the << 556
653 directory, "Lambda", << 557 if ( theLambdaTable && part == particle) {
654 verboseLevel, a << 558 const G4String name =
655 G4EmTableUtil::StoreTable(this, part, the << 559 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
656 directory, "LambdaPrim", << 560 yes = theLambdaTable->StorePhysicsTable(name,ascii);
657 verboseLevel, a << 561
658 return true; << 562 if ( yes ) {
>> 563 G4cout << "Physics tables are stored for " << particle->GetParticleName()
>> 564 << " and process " << GetProcessName()
>> 565 << " in the directory <" << directory
>> 566 << "> " << G4endl;
>> 567 } else {
>> 568 G4cout << "Fail to store Physics Tables for "
>> 569 << particle->GetParticleName()
>> 570 << " and process " << GetProcessName()
>> 571 << " in the directory <" << directory
>> 572 << "> " << G4endl;
>> 573 }
659 } 574 }
660 return false; << 575 return yes;
661 } 576 }
662 577
663 //....oooOO0OOooo........oooOO0OOooo........oo 578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664 579
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 580 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 cons << 581 const G4String& directory,
667 G4bo << 582 G4bool ascii)
668 { 583 {
669 if(!isTheMaster || part != particle) { retur << 584 if(1 < verboseLevel) {
670 G4bool yes = true; << 585 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
671 if(buildLambdaTable) { << 586 << part->GetParticleName() << " and process "
672 yes = G4EmTableUtil::RetrieveTable(this, p << 587 << GetProcessName() << G4endl;
673 "Lambda <<
674 ascii, <<
675 } <<
676 if(yes && minKinEnergyPrim < maxKinEnergy) { <<
677 yes = G4EmTableUtil::RetrieveTable(this, p <<
678 "Lambda <<
679 ascii, <<
680 } 588 }
681 return yes; << 589 G4bool yes = true;
682 } <<
683 590
684 //....oooOO0OOooo........oooOO0OOooo........oo << 591 if(!buildLambdaTable || particle != part) return yes;
685 592
686 G4double G4VEmProcess::GetCrossSection(G4doubl << 593 const G4String particleName = part->GetParticleName();
687 const G << 594 G4String filename;
688 { <<
689 CurrentSetup(couple, kinEnergy); <<
690 return GetCurrentLambda(kinEnergy, G4Log(kin <<
691 } <<
692 595
693 //....oooOO0OOooo........oooOO0OOooo........oo << 596 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
>> 597 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
>> 598 filename,ascii);
>> 599 if ( yes ) {
>> 600 if (0 < verboseLevel) {
>> 601 G4cout << "Lambda table for " << particleName
>> 602 << " is Retrieved from <"
>> 603 << filename << ">"
>> 604 << G4endl;
>> 605 }
>> 606 if((G4LossTableManager::Instance())->SplineFlag()) {
>> 607 size_t n = theLambdaTable->length();
>> 608 for(size_t i=0; i<n; ++i) {
>> 609 if((* theLambdaTable)[i]) {
>> 610 (* theLambdaTable)[i]->SetSpline(true);
>> 611 }
>> 612 }
>> 613 }
>> 614 } else {
>> 615 if (1 < verboseLevel) {
>> 616 G4cout << "Lambda table for " << particleName << " in file <"
>> 617 << filename << "> is not exist"
>> 618 << G4endl;
>> 619 }
>> 620 }
694 621
695 G4double G4VEmProcess::GetMeanFreePath(const G << 622 return yes;
696 G4doubl <<
697 G4Force <<
698 { <<
699 *condition = NotForced; <<
700 return G4VEmProcess::MeanFreePath(track); <<
701 } 623 }
702 624
703 //....oooOO0OOooo........oooOO0OOooo........oo 625 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 626
705 G4double << 627 void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou <<
707 G4dou <<
708 { 628 {
709 SelectModel(kinEnergy, currentCoupleIndex); << 629 G4RegionStore* regionStore = G4RegionStore::GetInstance();
710 return (currentModel) ? << 630 const G4Region* reg = r;
711 currentModel->ComputeCrossSectionPerAtom(c << 631 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
712 Z << 632
>> 633 // the region is in the list
>> 634 if (nDERegions) {
>> 635 for (G4int i=0; i<nDERegions; ++i) {
>> 636 if (reg == deRegions[i]) {
>> 637 if(!val) deRegions[i] = 0;
>> 638 return;
>> 639 }
>> 640 }
>> 641 }
>> 642
>> 643 // new region
>> 644 if(val) {
>> 645 useDeexcitation = true;
>> 646 deRegions.push_back(reg);
>> 647 nDERegions++;
>> 648 } else {
>> 649 useDeexcitation = false;
>> 650 }
713 } 651 }
714 652
715 //....oooOO0OOooo........oooOO0OOooo........oo 653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 654
717 G4PhysicsVector* << 655 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 656 const G4MaterialCutsCouple* couple)
719 { 657 {
>> 658 // Cross section per atom is calculated
720 DefineMaterial(couple); 659 DefineMaterial(couple);
721 G4PhysicsVector* newv = new G4PhysicsLogVect << 660 G4double cross = 0.0;
722 << 661 if(theLambdaTable) {
723 return newv; << 662 cross = (((*theLambdaTable)[currentCoupleIndex])->Value(kineticEnergy));
724 } << 663 } else {
725 << 664 SelectModel(kineticEnergy, currentCoupleIndex);
726 //....oooOO0OOooo........oooOO0OOooo........oo << 665 cross = currentModel->CrossSectionPerVolume(currentMaterial,
>> 666 currentParticle,kineticEnergy);
>> 667 }
727 668
728 const G4Element* G4VEmProcess::GetCurrentEleme << 669 if(cross < 0.0) { cross = 0.0; }
729 { << 670 return cross;
730 return (nullptr != currentModel) ? <<
731 currentModel->GetCurrentElement(currentMat <<
732 } 671 }
733 672
734 //....oooOO0OOooo........oooOO0OOooo........oo 673 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735 674
736 const G4Element* G4VEmProcess::GetTargetElemen << 675 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
>> 676 G4double,
>> 677 G4ForceCondition* condition)
737 { 678 {
738 return (nullptr != currentModel) ? << 679 *condition = NotForced;
739 currentModel->GetCurrentElement(currentMat << 680 return G4VEmProcess::MeanFreePath(track);
740 } 681 }
741 682
742 //....oooOO0OOooo........oooOO0OOooo........oo 683 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743 684
744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 685 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
745 { 686 {
746 return (nullptr != currentModel) ? << 687 DefineMaterial(track.GetMaterialCutsCouple());
747 currentModel->GetCurrentIsotope(GetCurrent << 688 preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
>> 689 G4double x = DBL_MAX;
>> 690 if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
>> 691 return x;
748 } 692 }
749 693
750 //....oooOO0OOooo........oooOO0OOooo........oo 694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 695
752 void G4VEmProcess::SetCrossSectionBiasingFacto << 696 G4double
>> 697 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy,
>> 698 G4double Z, G4double A, G4double cut)
753 { 699 {
754 if(f > 0.0) { << 700 SelectModel(kineticEnergy, currentCoupleIndex);
755 biasFactor = f; << 701 G4double x = 0.0;
756 weightFlag = flag; << 702 if(currentModel) {
757 if(1 < verboseLevel) { << 703 x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
758 G4cout << "### SetCrossSectionBiasingFac << 704 Z,A,cut);
759 << particle->GetParticleName() <<
760 << " and process " << GetProcessN <<
761 << " biasFactor= " << f << " weig <<
762 << G4endl; <<
763 } <<
764 } 705 }
>> 706 return x;
765 } 707 }
766 708
767 //....oooOO0OOooo........oooOO0OOooo........oo 709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768 710
769 void << 711 void G4VEmProcess::FindLambdaMax()
770 G4VEmProcess::ActivateForcedInteraction(G4doub <<
771 G4bool <<
772 { 712 {
773 if(nullptr == biasManager) { biasManager = n <<
774 if(1 < verboseLevel) { 713 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: << 714 G4cout << "### G4VEmProcess::FindLambdaMax: "
776 << particle->GetParticleName() << 715 << particle->GetParticleName()
777 << " and process " << GetProcessNam << 716 << " and process " << GetProcessName() << G4endl;
778 << " length(mm)= " << length/mm << 717 }
779 << " in G4Region <" << r << 718 size_t n = theLambdaTable->length();
780 << "> weightFlag= " << flag << 719 G4PhysicsVector* pv = (*theLambdaTable)[0];
781 << G4endl; << 720 G4double e, s, emax, smax;
782 } << 721 theEnergyOfCrossSectionMax = new G4double [n];
783 weightFlag = flag; << 722 theCrossSectionMax = new G4double [n];
784 biasManager->ActivateForcedInteraction(lengt << 723
785 } << 724 for (size_t i=0; i<n; ++i) {
786 << 725 pv = (*theLambdaTable)[i];
787 //....oooOO0OOooo........oooOO0OOooo........oo << 726 emax = DBL_MAX;
788 << 727 smax = 0.0;
789 void << 728 if(pv) {
790 G4VEmProcess::ActivateSecondaryBiasing(const G << 729 size_t nb = pv->GetVectorLength();
791 G4double factor, << 730 emax = DBL_MAX;
792 G4double energyLimit) << 731 smax = 0.0;
793 { << 732 if(nb > 0) {
794 if (0.0 <= factor) { << 733 for (size_t j=0; j<nb; ++j) {
795 << 734 e = pv->Energy(j);
796 // Range cut can be applied only for e- << 735 s = (*pv)(j);
797 if(0.0 == factor && secondaryParticle != G << 736 if(s > smax) {
798 { return; } << 737 smax = s;
799 << 738 emax = e;
800 if(!biasManager) { biasManager = new G4EmB << 739 }
801 biasManager->ActivateSecondaryBiasing(regi << 740 }
>> 741 }
>> 742 }
>> 743 theEnergyOfCrossSectionMax[i] = emax;
>> 744 theCrossSectionMax[i] = smax;
802 if(1 < verboseLevel) { 745 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: << 746 G4cout << "For " << particle->GetParticleName()
804 << " process " << GetProcessName() << 747 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
805 << " factor= " << factor << 748 << " lambda= " << smax << G4endl;
806 << " in G4Region <" << region <<
807 << "> energyLimit(MeV)= " << energyLimi <<
808 << G4endl; <<
809 } 749 }
810 } 750 }
811 } 751 }
812 752
813 //....oooOO0OOooo........oooOO0OOooo........oo 753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 754
815 void G4VEmProcess::SetLambdaBinning(G4int n) << 755 G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
816 { <<
817 if(5 < n && n < 10000000) { <<
818 nLambdaBins = n; <<
819 actBinning = true; <<
820 } else { <<
821 G4double e = (G4double)n; <<
822 PrintWarning("SetLambdaBinning", e); <<
823 } <<
824 } <<
825 <<
826 //....oooOO0OOooo........oooOO0OOooo........oo <<
827 <<
828 void G4VEmProcess::SetMinKinEnergy(G4double e) <<
829 { 756 {
830 if(1.e-3*eV < e && e < maxKinEnergy) { << 757 G4PhysicsVector* v =
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma << 758 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
832 /G4Log(maxKinEnergy/ << 759 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
833 minKinEnergy = e; << 760 return v;
834 actMinKinEnergy = true; <<
835 } else { PrintWarning("SetMinKinEnergy", e); <<
836 } 761 }
837 762
838 //....oooOO0OOooo........oooOO0OOooo........oo 763 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839 764
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) << 765 const G4Element* G4VEmProcess::GetCurrentElement() const
841 { <<
842 if(minKinEnergy < e && e < 1.e+6*TeV) { <<
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ <<
844 /G4Log(maxKinEnergy/ <<
845 maxKinEnergy = e; <<
846 actMaxKinEnergy = true; <<
847 } else { PrintWarning("SetMaxKinEnergy", e); <<
848 } <<
849 <<
850 //....oooOO0OOooo........oooOO0OOooo........oo <<
851 <<
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl <<
853 { <<
854 if(theParameters->MinKinEnergy() <= e && <<
855 e <= theParameters->MaxKinEnergy()) { min <<
856 else { PrintWarning("SetMinKinEnergyPrim", e <<
857 } <<
858 <<
859 //....oooOO0OOooo........oooOO0OOooo........oo <<
860 <<
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const <<
862 { <<
863 return (nam == GetProcessName()) ? this : nu <<
864 } <<
865 <<
866 //....oooOO0OOooo........oooOO0OOooo........oo <<
867 <<
868 G4double G4VEmProcess::PolarAngleLimit() const <<
869 { <<
870 return theParameters->MscThetaLimit(); <<
871 } <<
872 <<
873 //....oooOO0OOooo........oooOO0OOooo........oo <<
874 <<
875 void G4VEmProcess::PrintWarning(G4String tit, <<
876 { <<
877 G4String ss = "G4VEmProcess::" + tit; <<
878 G4ExceptionDescription ed; <<
879 ed << "Parameter is out of range: " << val <<
880 << " it will have no effect!\n" << " Pro <<
881 << GetProcessName() << " nbins= " << the <<
882 << " Emin(keV)= " << theParameters->MinKi <<
883 << " Emax(GeV)= " << theParameters->MaxKi <<
884 G4Exception(ss, "em0044", JustWarning, ed); <<
885 } <<
886 <<
887 //....oooOO0OOooo........oooOO0OOooo........oo <<
888 <<
889 void G4VEmProcess::ProcessDescription(std::ost <<
890 { 766 {
891 if(nullptr != particle) { << 767 const G4Element* elm = 0;
892 StreamInfo(out, *particle, true); << 768 if(currentModel) {elm = currentModel->GetCurrentElement(); }
893 } << 769 return elm;
894 } 770 }
895 771
896 //....oooOO0OOooo........oooOO0OOooo........oo 772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 773