| 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 98778 2016-08-09 14:41:08Z gcosmo $
>> 27 //
26 // ------------------------------------------- 28 // -------------------------------------------------------------------
27 // 29 //
28 // GEANT4 Class file 30 // GEANT4 Class file
29 // 31 //
30 // 32 //
31 // File name: G4VEmProcess 33 // File name: G4VEmProcess
32 // 34 //
33 // Author: Vladimir Ivanchenko on base 35 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 // 36 //
35 // Creation date: 01.10.2003 37 // Creation date: 01.10.2003
36 // 38 //
37 // Modifications: by V.Ivanchenko << 39 // Modifications:
>> 40 // 30-06-04 make it to be pure discrete process (V.Ivanchenko)
>> 41 // 30-09-08 optimise integral option (V.Ivanchenko)
>> 42 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
>> 43 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
>> 44 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
>> 45 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
>> 46 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
>> 47 // 25-07-05 Add protection: integral mode only for charged particles (VI)
>> 48 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
>> 49 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
>> 50 // 12-09-06 add SetModel() (mma)
>> 51 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
>> 52 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
>> 53 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
>> 54 // 17-02-10 Added pointer currentParticle (VI)
>> 55 // 30-05-12 allow Russian roulette, brem splitting (D. Sawkey)
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" 66 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 67 #include "G4SystemOfUnits.hh"
50 #include "G4ProcessManager.hh" 68 #include "G4ProcessManager.hh"
51 #include "G4LossTableManager.hh" 69 #include "G4LossTableManager.hh"
52 #include "G4LossTableBuilder.hh" 70 #include "G4LossTableBuilder.hh"
53 #include "G4Step.hh" 71 #include "G4Step.hh"
54 #include "G4ParticleDefinition.hh" 72 #include "G4ParticleDefinition.hh"
55 #include "G4VEmModel.hh" 73 #include "G4VEmModel.hh"
56 #include "G4DataVector.hh" 74 #include "G4DataVector.hh"
57 #include "G4PhysicsTable.hh" 75 #include "G4PhysicsTable.hh"
58 #include "G4EmDataHandler.hh" <<
59 #include "G4PhysicsLogVector.hh" 76 #include "G4PhysicsLogVector.hh"
60 #include "G4VParticleChange.hh" 77 #include "G4VParticleChange.hh"
61 #include "G4ProductionCutsTable.hh" 78 #include "G4ProductionCutsTable.hh"
62 #include "G4Region.hh" 79 #include "G4Region.hh"
63 #include "G4Gamma.hh" 80 #include "G4Gamma.hh"
64 #include "G4Electron.hh" 81 #include "G4Electron.hh"
65 #include "G4Positron.hh" 82 #include "G4Positron.hh"
66 #include "G4PhysicsTableHelper.hh" 83 #include "G4PhysicsTableHelper.hh"
67 #include "G4EmBiasingManager.hh" 84 #include "G4EmBiasingManager.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" 85 #include "G4GenericIon.hh"
74 #include "G4Log.hh" 86 #include "G4Log.hh"
75 #include <iostream> 87 #include <iostream>
76 88
77 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 90
79 G4VEmProcess::G4VEmProcess(const G4String& nam 91 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
80 G4VDiscreteProcess(name, type) << 92 G4VDiscreteProcess(name, type),
>> 93 secondaryParticle(nullptr),
>> 94 buildLambdaTable(true),
>> 95 numberOfModels(0),
>> 96 theLambdaTable(nullptr),
>> 97 theLambdaTablePrim(nullptr),
>> 98 theDensityFactor(nullptr),
>> 99 theDensityIdx(nullptr),
>> 100 integral(false),
>> 101 applyCuts(false),
>> 102 startFromNull(false),
>> 103 splineFlag(true),
>> 104 currentModel(nullptr),
>> 105 particle(nullptr),
>> 106 currentParticle(nullptr),
>> 107 currentCouple(nullptr)
81 { 108 {
82 theParameters = G4EmParameters::Instance(); 109 theParameters = G4EmParameters::Instance();
83 SetVerboseLevel(1); 110 SetVerboseLevel(1);
84 111
85 // Size of tables << 112 // Size of tables assuming spline
86 minKinEnergy = 0.1*CLHEP::keV; << 113 minKinEnergy = 0.1*keV;
87 maxKinEnergy = 100.0*CLHEP::TeV; << 114 maxKinEnergy = 100.0*TeV;
>> 115 nLambdaBins = 77;
>> 116 minKinEnergyPrim = DBL_MAX;
>> 117 actBinning = actSpline = actMinKinEnergy = actMaxKinEnergy = false;
88 118
89 // default lambda factor 119 // default lambda factor
90 invLambdaFactor = 1.0/lambdaFactor; << 120 lambdaFactor = 0.8;
>> 121
>> 122 // default limit on polar angle
>> 123 biasFactor = fFactor = 1.0;
91 124
92 // particle types 125 // particle types
93 theGamma = G4Gamma::Gamma(); << 126 theGamma = G4Gamma::Gamma();
94 theElectron = G4Electron::Electron(); << 127 theElectron = G4Electron::Electron();
95 thePositron = G4Positron::Positron(); << 128 thePositron = G4Positron::Positron();
>> 129
>> 130 theCuts = theCutsGamma = theCutsElectron = theCutsPositron = nullptr;
96 131
97 pParticleChange = &fParticleChange; 132 pParticleChange = &fParticleChange;
98 fParticleChange.SetSecondaryWeightByProcess( 133 fParticleChange.SetSecondaryWeightByProcess(true);
99 secParticles.reserve(5); 134 secParticles.reserve(5);
100 135
>> 136 baseMaterial = currentMaterial = nullptr;
>> 137
>> 138 preStepLambda = preStepKinEnergy = 0.0;
>> 139 mfpKinEnergy = DBL_MAX;
>> 140
>> 141 idxLambda = idxLambdaPrim = currentCoupleIndex
>> 142 = basedCoupleIndex = 0;
>> 143
101 modelManager = new G4EmModelManager(); 144 modelManager = new G4EmModelManager();
>> 145 biasManager = nullptr;
>> 146 biasFlag = false;
>> 147 weightFlag = false;
102 lManager = G4LossTableManager::Instance(); 148 lManager = G4LossTableManager::Instance();
103 lManager->Register(this); 149 lManager->Register(this);
104 isTheMaster = lManager->IsMaster(); << 150 secID = fluoID = augerID = biasID = -1;
105 G4LossTableBuilder* bld = lManager->GetTable << 151 mainSecondaries = 100;
106 theDensityFactor = bld->GetDensityFactors(); << 152 if("phot" == GetProcessName() || "compt" == GetProcessName()) {
107 theDensityIdx = bld->GetCoupleIndexes(); << 153 mainSecondaries = 1;
>> 154 }
108 } 155 }
109 156
110 //....oooOO0OOooo........oooOO0OOooo........oo 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 158
112 G4VEmProcess::~G4VEmProcess() 159 G4VEmProcess::~G4VEmProcess()
113 { 160 {
114 if(isTheMaster) { << 161 /*
115 delete theData; << 162 if(1 < verboseLevel) {
116 delete theEnergyOfCrossSectionMax; << 163 G4cout << "G4VEmProcess destruct " << GetProcessName()
>> 164 << " " << this << " " << theLambdaTable <<G4endl;
>> 165 }
>> 166 */
>> 167 if(lManager->IsMaster()) {
>> 168 if(theLambdaTable) {
>> 169 theLambdaTable->clearAndDestroy();
>> 170 delete theLambdaTable;
>> 171 }
>> 172 if(theLambdaTablePrim) {
>> 173 theLambdaTablePrim->clearAndDestroy();
>> 174 delete theLambdaTablePrim;
>> 175 }
117 } 176 }
118 delete modelManager; 177 delete modelManager;
119 delete biasManager; 178 delete biasManager;
120 lManager->DeRegister(this); 179 lManager->DeRegister(this);
>> 180 //G4cout << "G4VEmProcess removed " << G4endl;
121 } 181 }
122 182
123 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 184
125 void G4VEmProcess::AddEmModel(G4int order, G4V << 185 void G4VEmProcess::Clear()
>> 186 {
>> 187 currentCouple = nullptr;
>> 188 preStepLambda = 0.0;
>> 189 mfpKinEnergy = DBL_MAX;
>> 190 idxLambda = idxLambdaPrim = 0;
>> 191 }
>> 192
>> 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 194
>> 195 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
>> 196 const G4Material*)
>> 197 {
>> 198 return 0.0;
>> 199 }
>> 200
>> 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 202
>> 203 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p,
126 const G4Region* 204 const G4Region* region)
127 { 205 {
128 if(nullptr == ptr) { return; } <<
129 G4VEmFluctuationModel* fm = nullptr; 206 G4VEmFluctuationModel* fm = nullptr;
130 modelManager->AddEmModel(order, ptr, fm, reg << 207 modelManager->AddEmModel(order, p, fm, region);
131 ptr->SetParticleChange(pParticleChange); << 208 if(p) { p->SetParticleChange(pParticleChange); }
132 } 209 }
133 210
134 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 212
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 213 void G4VEmProcess::SetEmModel(G4VEmModel* p, G4int index)
137 { 214 {
138 if(nullptr == ptr) { return; } << 215 G4int n = emModels.size();
139 if(!emModels.empty()) { << 216 if(index >= n) {
140 for(auto & em : emModels) { if(em == ptr) << 217 for(G4int i=n; i<=index; ++i) {emModels.push_back(nullptr);}
141 } 218 }
142 emModels.push_back(ptr); << 219 emModels[index] = p;
>> 220 }
>> 221
>> 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 223
>> 224 G4VEmModel* G4VEmProcess::EmModel(G4int index) const
>> 225 {
>> 226 G4VEmModel* p = nullptr;
>> 227 if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
>> 228 return p;
>> 229 }
>> 230
>> 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 232
>> 233 void G4VEmProcess::UpdateEmModel(const G4String& nam,
>> 234 G4double emin, G4double emax)
>> 235 {
>> 236 modelManager->UpdateEmModel(nam, emin, emax);
>> 237 }
>> 238
>> 239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 240
>> 241 G4int G4VEmProcess::GetNumberOfModels() const
>> 242 {
>> 243 return modelManager->NumberOfModels();
>> 244 }
>> 245
>> 246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 247
>> 248 G4int G4VEmProcess::GetNumberOfRegionModels(size_t couple_index) const
>> 249 {
>> 250 return modelManager->NumberOfRegionModels(couple_index);
>> 251 }
>> 252
>> 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 254
>> 255 G4VEmModel* G4VEmProcess::GetRegionModel(G4int idx, size_t couple_index) const
>> 256 {
>> 257 return modelManager->GetRegionModel(idx, couple_index);
>> 258 }
>> 259
>> 260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 261
>> 262 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver) const
>> 263 {
>> 264 return modelManager->GetModel(idx, ver);
143 } 265 }
144 266
145 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 268
147 void G4VEmProcess::PreparePhysicsTable(const G 269 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
148 { 270 {
149 if(nullptr == particle) { SetParticle(&part) << 271 G4bool isMaster = true;
>> 272 const G4VEmProcess* masterProcess =
>> 273 static_cast<const G4VEmProcess*>(GetMasterProcess());
>> 274 if(masterProcess && masterProcess != this) { isMaster = false; }
>> 275
>> 276 if(!particle) { SetParticle(&part); }
150 277
151 if(part.GetParticleType() == "nucleus" && 278 if(part.GetParticleType() == "nucleus" &&
152 part.GetParticleSubType() == "generic") { 279 part.GetParticleSubType() == "generic") {
153 280
154 G4String pname = part.GetParticleName(); 281 G4String pname = part.GetParticleName();
155 if(pname != "deuteron" && pname != "triton 282 if(pname != "deuteron" && pname != "triton" &&
156 pname != "He3" && pname != "alpha" && p << 283 pname != "alpha" && pname != "He3" &&
157 pname != "helium" && pname != "hydrogen << 284 pname != "alpha+" && pname != "helium" &&
>> 285 pname != "hydrogen") {
158 286
159 particle = G4GenericIon::GenericIon(); 287 particle = G4GenericIon::GenericIon();
160 isIon = true; <<
161 } 288 }
162 } 289 }
>> 290
>> 291 if(1 < verboseLevel) {
>> 292 G4cout << "G4VEmProcess::PreparePhysicsTable() for "
>> 293 << GetProcessName()
>> 294 << " and particle " << part.GetParticleName()
>> 295 << " local particle " << particle->GetParticleName()
>> 296 << G4endl;
>> 297 }
>> 298
163 if(particle != &part) { return; } 299 if(particle != &part) { return; }
164 300
165 lManager->PreparePhysicsTable(&part, this); << 301 G4LossTableBuilder* bld = lManager->GetTableBuilder();
166 302
167 // for new run << 303 lManager->PreparePhysicsTable(&part, this, isMaster);
168 currentCouple = nullptr; <<
169 preStepLambda = 0.0; <<
170 fLambdaEnergy = 0.0; <<
171 304
>> 305 Clear();
172 InitialiseProcess(particle); 306 InitialiseProcess(particle);
173 307
174 G4LossTableBuilder* bld = lManager->GetTable <<
175 const G4ProductionCutsTable* theCoupleTable= 308 const G4ProductionCutsTable* theCoupleTable=
176 G4ProductionCutsTable::GetProductionCutsTa 309 G4ProductionCutsTable::GetProductionCutsTable();
177 theCutsGamma = theCoupleTable->GetEnergyC << 310 size_t n = theCoupleTable->GetTableSize();
178 theCutsElectron = theCoupleTable->GetEnergyC << 311
179 theCutsPositron = theCoupleTable->GetEnergyC << 312 theEnergyOfCrossSectionMax.resize(n, 0.0);
>> 313 theCrossSectionMax.resize(n, DBL_MAX);
180 314
181 // initialisation of the process 315 // initialisation of the process
182 if(!actMinKinEnergy) { minKinEnergy = thePar 316 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
183 if(!actMaxKinEnergy) { maxKinEnergy = thePar 317 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
>> 318 if(!actSpline) { splineFlag = theParameters->Spline(); }
184 319
185 applyCuts = theParameters->ApplyCuts() << 320 if(isMaster) { SetVerboseLevel(theParameters->Verbose()); }
186 lambdaFactor = theParameters->LambdaFacto << 321 else { SetVerboseLevel(theParameters->WorkerVerbose()); }
187 invLambdaFactor = 1.0/lambdaFactor; << 322 applyCuts = theParameters->ApplyCuts();
>> 323 lambdaFactor = theParameters->LambdaFactor();
188 theParameters->DefineRegParamForEM(this); 324 theParameters->DefineRegParamForEM(this);
189 325
190 // integral option may be disabled << 326 // initialisation of models
191 if(!theParameters->Integral()) { fXSType = f << 327 numberOfModels = modelManager->NumberOfModels();
192 << 328 for(G4int i=0; i<numberOfModels; ++i) {
193 // prepare tables << 329 G4VEmModel* mod = modelManager->GetModel(i);
194 if(isTheMaster) { << 330 if(0 == i) { currentModel = mod; }
195 if(nullptr == theData) { theData = new G4E << 331 mod->SetPolarAngleLimit(theParameters->MscThetaLimit());
196 << 332 mod->SetMasterThread(isMaster);
197 if(buildLambdaTable) { << 333 if(mod->HighEnergyLimit() > maxKinEnergy) {
198 theLambdaTable = theData->MakeTable(0); << 334 mod->SetHighEnergyLimit(maxKinEnergy);
199 bld->InitialiseBaseMaterials(theLambdaTa <<
200 } <<
201 // high energy table <<
202 if(minKinEnergyPrim < maxKinEnergy) { <<
203 theLambdaTablePrim = theData->MakeTable( <<
204 bld->InitialiseBaseMaterials(theLambdaTa <<
205 } 335 }
206 } 336 }
207 // models << 337
208 baseMat = bld->GetBaseMaterialFlag(); << 338 if(lManager->AtomDeexcitation()) { modelManager->SetFluoFlag(true); }
209 numberOfModels = modelManager->NumberOfModel << 339 theCuts = modelManager->Initialise(particle,secondaryParticle,
210 currentModel = modelManager->GetModel(0); << 340 2.,verboseLevel);
211 if(nullptr != lManager->AtomDeexcitation()) << 341 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
212 modelManager->SetFluoFlag(true); << 342 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
>> 343 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
>> 344
>> 345 // prepare tables
>> 346 if(buildLambdaTable && isMaster){
>> 347 theLambdaTable =
>> 348 G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
>> 349 bld->InitialiseBaseMaterials(theLambdaTable);
>> 350 }
>> 351 // high energy table
>> 352 if(isMaster && minKinEnergyPrim < maxKinEnergy){
>> 353 theLambdaTablePrim =
>> 354 G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTablePrim);
>> 355 bld->InitialiseBaseMaterials(theLambdaTablePrim);
213 } 356 }
214 // forced biasing 357 // forced biasing
215 if(nullptr != biasManager) { << 358 if(biasManager) {
216 biasManager->Initialise(part, GetProcessNa << 359 biasManager->Initialise(part,GetProcessName(),verboseLevel);
217 biasFlag = false; << 360 biasFlag = false;
218 } << 361 }
219 << 362 // defined ID of secondary particles
220 theCuts = << 363 G4String nam1 = GetProcessName();
221 G4EmTableUtil::PrepareEmProcess(this, part << 364 G4String nam2 = nam1 + "_fluo" ;
222 modelManag << 365 G4String nam3 = nam1 + "_auger";
223 secID, tri << 366 G4String nam4 = nam1 + "_split";
224 verboseLev << 367 secID = G4PhysicsModelCatalog::Register(nam1);
>> 368 fluoID = G4PhysicsModelCatalog::Register(nam2);
>> 369 augerID = G4PhysicsModelCatalog::Register(nam3);
>> 370 biasID = G4PhysicsModelCatalog::Register(nam4);
225 } 371 }
226 372
227 //....oooOO0OOooo........oooOO0OOooo........oo 373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228 374
229 void G4VEmProcess::BuildPhysicsTable(const G4P 375 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 { 376 {
231 if(nullptr == masterProc) { << 377 G4bool isMaster = true;
232 if(isTheMaster) { masterProc = this; } << 378 const G4VEmProcess* masterProc =
233 else { masterProc = static_cast<const G4VE << 379 static_cast<const G4VEmProcess*>(GetMasterProcess());
234 } << 380 if(masterProc && masterProc != this) { isMaster = false; }
235 G4int nModels = modelManager->NumberOfModels << 381
236 G4bool isLocked = theParameters->IsPrintLock << 382 G4String num = part.GetParticleName();
237 G4bool toBuild = (buildLambdaTable || minKin << 383 if(1 < verboseLevel) {
238 << 384 G4cout << "### G4VEmProcess::BuildPhysicsTable() for "
239 G4EmTableUtil::BuildEmProcess(this, masterPr << 385 << GetProcessName()
240 nModels, verbo << 386 << " and particle " << num
241 isLocked, toBu << 387 << " buildLambdaTable= " << buildLambdaTable
>> 388 << " isMaster= " << isMaster
>> 389 << G4endl;
>> 390 }
>> 391
>> 392 if(particle == &part) {
>> 393
>> 394 G4LossTableBuilder* bld = lManager->GetTableBuilder();
>> 395
>> 396 // worker initialisation
>> 397 if(!isMaster) {
>> 398 theLambdaTable = masterProc->LambdaTable();
>> 399 theLambdaTablePrim = masterProc->LambdaTablePrim();
>> 400
>> 401 if(theLambdaTable) {
>> 402 bld->InitialiseBaseMaterials(theLambdaTable);
>> 403 } else if(theLambdaTablePrim) {
>> 404 bld->InitialiseBaseMaterials(theLambdaTablePrim);
>> 405 }
>> 406 theDensityFactor = bld->GetDensityFactors();
>> 407 theDensityIdx = bld->GetCoupleIndexes();
>> 408 if(theLambdaTable) { FindLambdaMax(); }
>> 409
>> 410 // local initialisation of models
>> 411 G4bool printing = true;
>> 412 numberOfModels = modelManager->NumberOfModels();
>> 413 for(G4int i=0; i<numberOfModels; ++i) {
>> 414 G4VEmModel* mod = GetModelByIndex(i, printing);
>> 415 G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing);
>> 416 mod->InitialiseLocal(particle, mod0);
>> 417 }
>> 418 // master thread
>> 419 } else {
>> 420 theDensityFactor = bld->GetDensityFactors();
>> 421 theDensityIdx = bld->GetCoupleIndexes();
>> 422 if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
>> 423 BuildLambdaTable();
>> 424 }
>> 425 }
>> 426 }
>> 427
>> 428 // explicitly defined printout by particle name
>> 429 if(1 < verboseLevel ||
>> 430 (0 < verboseLevel && (num == "gamma" || num == "e-" ||
>> 431 num == "e+" || num == "mu+" ||
>> 432 num == "mu-" || num == "proton"||
>> 433 num == "pi+" || num == "pi-" ||
>> 434 num == "kaon+" || num == "kaon-" ||
>> 435 num == "alpha" || num == "anti_proton" ||
>> 436 num == "GenericIon")))
>> 437 {
>> 438 PrintInfoProcess(part);
>> 439 }
>> 440
>> 441 if(1 < verboseLevel) {
>> 442 G4cout << "### G4VEmProcess::BuildPhysicsTable() done for "
>> 443 << GetProcessName()
>> 444 << " and particle " << num
>> 445 << G4endl;
>> 446 }
242 } 447 }
243 448
244 //....oooOO0OOooo........oooOO0OOooo........oo 449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 450
246 void G4VEmProcess::BuildLambdaTable() 451 void G4VEmProcess::BuildLambdaTable()
247 { 452 {
>> 453 if(1 < verboseLevel) {
>> 454 G4cout << "G4EmProcess::BuildLambdaTable() for process "
>> 455 << GetProcessName() << " and particle "
>> 456 << particle->GetParticleName() << " " << this
>> 457 << G4endl;
>> 458 }
>> 459
>> 460 // Access to materials
>> 461 const G4ProductionCutsTable* theCoupleTable=
>> 462 G4ProductionCutsTable::GetProductionCutsTable();
>> 463 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 464
>> 465 G4LossTableBuilder* bld = lManager->GetTableBuilder();
>> 466
>> 467 G4PhysicsLogVector* aVector = nullptr;
>> 468 G4PhysicsLogVector* aVectorPrim = nullptr;
>> 469 G4PhysicsLogVector* bVectorPrim = nullptr;
>> 470
248 G4double scale = theParameters->MaxKinEnergy 471 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
249 G4int nbin = << 472 G4int nbin = theParameters->NumberOfBinsPerDecade()
250 theParameters->NumberOfBinsPerDecade()*G4l << 473 *G4lrint(std::log10(scale));
>> 474 scale = G4Log(scale);
251 if(actBinning) { nbin = std::max(nbin, nLamb 475 if(actBinning) { nbin = std::max(nbin, nLambdaBins); }
252 scale = nbin/G4Log(scale); << 476 G4double emax1 = std::min(maxKinEnergy, minKinEnergyPrim);
253 << 477
254 G4LossTableBuilder* bld = lManager->GetTable << 478 for(size_t i=0; i<numOfCouples; ++i) {
255 G4EmTableUtil::BuildLambdaTable(this, partic << 479
256 bld, theLamb << 480 if (bld->GetFlag(i)) {
257 minKinEnergy << 481
258 maxKinEnergy << 482 // create physics vector and fill it
259 startFromNul << 483 const G4MaterialCutsCouple* couple =
>> 484 theCoupleTable->GetMaterialCutsCouple(i);
>> 485
>> 486 // build main table
>> 487 if(buildLambdaTable) {
>> 488 delete (*theLambdaTable)[i];
>> 489
>> 490 // if start from zero then change the scale
>> 491 G4double emin = minKinEnergy;
>> 492 G4bool startNull = false;
>> 493 if(startFromNull) {
>> 494 G4double e = MinPrimaryEnergy(particle,couple->GetMaterial());
>> 495 if(e >= emin) {
>> 496 emin = e;
>> 497 startNull = true;
>> 498 }
>> 499 }
>> 500 G4double emax = emax1;
>> 501 if(emax <= emin) { emax = 2*emin; }
>> 502 G4int bin = G4lrint(nbin*G4Log(emax/emin)/scale);
>> 503 if(bin < 3) { bin = 3; }
>> 504 aVector = new G4PhysicsLogVector(emin, emax, bin);
>> 505 aVector->SetSpline(splineFlag);
>> 506 modelManager->FillLambdaVector(aVector, couple, startNull);
>> 507 if(splineFlag) { aVector->FillSecondDerivatives(); }
>> 508 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
>> 509 }
>> 510 // build high energy table
>> 511 if(minKinEnergyPrim < maxKinEnergy) {
>> 512 delete (*theLambdaTablePrim)[i];
>> 513
>> 514 // start not from zero
>> 515 if(!bVectorPrim) {
>> 516 G4int bin = G4lrint(nbin*G4Log(maxKinEnergy/minKinEnergyPrim)/scale);
>> 517 if(bin < 3) { bin = 3; }
>> 518 aVectorPrim =
>> 519 new G4PhysicsLogVector(minKinEnergyPrim, maxKinEnergy, bin);
>> 520 bVectorPrim = aVectorPrim;
>> 521 } else {
>> 522 aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
>> 523 }
>> 524 // always use spline
>> 525 aVectorPrim->SetSpline(splineFlag);
>> 526 modelManager->FillLambdaVector(aVectorPrim, couple, false,
>> 527 fIsCrossSectionPrim);
>> 528 aVectorPrim->FillSecondDerivatives();
>> 529 G4PhysicsTableHelper::SetPhysicsVector(theLambdaTablePrim, i,
>> 530 aVectorPrim);
>> 531 }
>> 532 }
>> 533 }
>> 534
>> 535 if(buildLambdaTable) { FindLambdaMax(); }
>> 536
>> 537 if(1 < verboseLevel) {
>> 538 G4cout << "Lambda table is built for "
>> 539 << particle->GetParticleName()
>> 540 << G4endl;
>> 541 }
260 } 542 }
261 543
262 //....oooOO0OOooo........oooOO0OOooo........oo 544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 545
264 void G4VEmProcess::StreamInfo(std::ostream& ou << 546 void G4VEmProcess::PrintInfoProcess(const G4ParticleDefinition& part)
265 const G4ParticleDefinition& << 547 {
266 { << 548 if(verboseLevel > 0) {
267 G4String indent = (rst ? " " : ""); << 549 G4cout << std::setprecision(6);
268 out << std::setprecision(6); << 550 G4cout << G4endl << GetProcessName() << ": for "
269 out << G4endl << indent << GetProcessName() << 551 << part.GetParticleName();
270 if (!rst) { << 552 if(integral) { G4cout << ", integral: 1 "; }
271 out << " for " << part.GetParticleName(); << 553 if(applyCuts) { G4cout << ", applyCuts: 1 "; }
272 } << 554 G4cout << " SubType= " << GetProcessSubType();;
273 if(fXSType != fEmNoIntegral) { out << " XSt << 555 if(biasFactor != 1.0) { G4cout << " BiasingFactor= " << biasFactor; }
274 if(applyCuts) { out << " applyCuts:1 "; } << 556 G4cout << " BuildTable= " << buildLambdaTable;
275 G4int subtype = GetProcessSubType(); << 557 G4cout << G4endl;
276 out << " SubType=" << subtype; << 558 if(buildLambdaTable) {
277 if (subtype == fAnnihilation) { << 559 if(particle == &part) {
278 G4int mod = theParameters->PositronAtRestM << 560 size_t length = theLambdaTable->length();
279 const G4String namp[2] = {"Simple", "Allis << 561 for(size_t i=0; i<length; ++i) {
280 out << " AtRestModel:" << namp[mod]; << 562 G4PhysicsVector* v = (*theLambdaTable)[i];
281 } << 563 if(v) {
282 if(biasFactor != 1.0) { out << " BiasingFac << 564 G4cout << " Lambda table from ";
283 out << " BuildTable=" << buildLambdaTable << << 565 G4double emin = v->Energy(0);
284 if(buildLambdaTable) { << 566 G4double emax = v->GetMaxEnergy();
285 if(particle == &part) { << 567 G4int nbin = v->GetVectorLength() - 1;
286 for(auto & v : *theLambdaTable) { << 568 if(emin > minKinEnergy) { G4cout << "threshold "; }
287 if(nullptr != v) { << 569 else { G4cout << G4BestUnit(emin,"Energy"); }
288 out << " Lambda table from "; << 570 G4cout << " to "
289 G4double emin = v->Energy(0); << 571 << G4BestUnit(emax,"Energy")
290 G4double emax = v->GetMaxEnergy(); << 572 << ", " << G4lrint(nbin/std::log10(emax/emin))
291 G4int nbin = G4int(v->GetVectorLengt << 573 << " bins per decade, spline: "
292 if(emin > minKinEnergy) { out << "th << 574 << splineFlag
293 else { out << G4BestUnit(emin,"Energ << 575 << G4endl;
294 out << " to " << 576 break;
295 << G4BestUnit(emax,"Energy") << 577 }
296 << ", " << G4lrint(nbin/std::log <<
297 << " bins/decade, spline: " <<
298 << splineFlag << G4endl; <<
299 break; <<
300 } 578 }
>> 579 } else {
>> 580 G4cout << " Used Lambda table of "
>> 581 << particle->GetParticleName() << G4endl;;
301 } 582 }
302 } else { <<
303 out << " Used Lambda table of " <<
304 << particle->GetParticleName() << G4endl <<
305 } 583 }
306 } << 584 if(minKinEnergyPrim < maxKinEnergy) {
307 if(minKinEnergyPrim < maxKinEnergy) { << 585 if(particle == &part) {
308 if(particle == &part) { << 586 size_t length = theLambdaTablePrim->length();
309 for(auto & v : *theLambdaTablePrim) { << 587 for(size_t i=0; i<length; ++i) {
310 if(nullptr != v) { << 588 G4PhysicsVector* v = (*theLambdaTablePrim)[i];
311 out << " LambdaPrime table from << 589 if(v) {
312 << G4BestUnit(v->Energy(0),"Ener << 590 G4cout << " LambdaPrime table from "
313 << " to " << 591 << G4BestUnit(v->Energy(0),"Energy")
314 << G4BestUnit(v->GetMaxEnergy(), << 592 << " to "
315 << " in " << v->GetVectorLength( << 593 << G4BestUnit(v->GetMaxEnergy(),"Energy")
316 << " bins " << G4endl; << 594 << " in " << v->GetVectorLength()-1
317 break; << 595 << " bins "
>> 596 << G4endl;
>> 597 break;
>> 598 }
318 } 599 }
>> 600 } else {
>> 601 G4cout << " Used LambdaPrime table of "
>> 602 << particle->GetParticleName() << G4endl;;
319 } 603 }
320 } else { <<
321 out << " Used LambdaPrime table of <<
322 << particle->GetParticleName() <<
323 } 604 }
>> 605 PrintInfo();
>> 606 modelManager->DumpModelList(verboseLevel);
324 } 607 }
325 StreamProcessInfo(out); <<
326 modelManager->DumpModelList(out, verboseLeve <<
327 608
328 if(verboseLevel > 2 && buildLambdaTable) { 609 if(verboseLevel > 2 && buildLambdaTable) {
329 out << " LambdaTable address= " << th << 610 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
330 if(theLambdaTable && particle == &part) { 611 if(theLambdaTable && particle == &part) {
331 out << (*theLambdaTable) << G4endl; << 612 G4cout << (*theLambdaTable) << G4endl;
332 } 613 }
333 } 614 }
334 } 615 }
335 616
336 //....oooOO0OOooo........oooOO0OOooo........oo 617 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337 618
338 void G4VEmProcess::StartTracking(G4Track* trac 619 void G4VEmProcess::StartTracking(G4Track* track)
339 { 620 {
340 // reset parameters for the new track 621 // reset parameters for the new track
341 currentParticle = track->GetParticleDefiniti 622 currentParticle = track->GetParticleDefinition();
342 theNumberOfInteractionLengthLeft = -1.0; 623 theNumberOfInteractionLengthLeft = -1.0;
343 mfpKinEnergy = DBL_MAX; << 624 mfpKinEnergy = DBL_MAX;
344 preStepLambda = 0.0; <<
345 <<
346 if(isIon) { massRatio = proton_mass_c2/curre <<
347 625
348 // forced biasing only for primary particles 626 // forced biasing only for primary particles
349 if(biasManager) { 627 if(biasManager) {
350 if(0 == track->GetParentID()) { 628 if(0 == track->GetParentID()) {
351 // primary particle 629 // primary particle
352 biasFlag = true; 630 biasFlag = true;
353 biasManager->ResetForcedInteraction(); 631 biasManager->ResetForcedInteraction();
354 } 632 }
355 } 633 }
356 } 634 }
357 635
358 //....oooOO0OOooo........oooOO0OOooo........oo 636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 637
360 G4double G4VEmProcess::PostStepGetPhysicalInte 638 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& tr 639 const G4Track& track,
362 G4double previo 640 G4double previousStepSize,
363 G4ForceCondition* 641 G4ForceCondition* condition)
364 { 642 {
365 *condition = NotForced; 643 *condition = NotForced;
366 G4double x = DBL_MAX; 644 G4double x = DBL_MAX;
367 645
368 DefineMaterial(track.GetMaterialCutsCouple() <<
369 preStepKinEnergy = track.GetKineticEnergy(); 646 preStepKinEnergy = track.GetKineticEnergy();
370 const G4double scaledEnergy = preStepKinEner << 647 DefineMaterial(track.GetMaterialCutsCouple());
371 SelectModel(scaledEnergy, currentCoupleIndex << 648 SelectModel(preStepKinEnergy, currentCoupleIndex);
372 /* << 649
373 G4cout << "PostStepGetPhysicalInteractionLen << 650 if(!currentModel->IsActive(preStepKinEnergy)) {
374 << " couple: " << currentCouple << G <<
375 */ <<
376 if(!currentModel->IsActive(scaledEnergy)) { <<
377 theNumberOfInteractionLengthLeft = -1.0; 651 theNumberOfInteractionLengthLeft = -1.0;
378 currentInteractionLength = DBL_MAX; 652 currentInteractionLength = DBL_MAX;
379 mfpKinEnergy = DBL_MAX; <<
380 preStepLambda = 0.0; <<
381 return x; 653 return x;
382 } 654 }
383 655
384 // forced biasing only for primary particles 656 // forced biasing only for primary particles
385 if(biasManager) { 657 if(biasManager) {
386 if(0 == track.GetParentID()) { 658 if(0 == track.GetParentID()) {
387 if(biasFlag && 659 if(biasFlag &&
388 biasManager->ForcedInteractionRegion( << 660 biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
389 return biasManager->GetStepLimit((G4in << 661 return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
390 } 662 }
391 } 663 }
392 } 664 }
393 665
394 // compute mean free path 666 // compute mean free path
>> 667 if(preStepKinEnergy < mfpKinEnergy) {
>> 668 if (integral) { ComputeIntegralLambda(preStepKinEnergy); }
>> 669 else { preStepLambda = GetCurrentLambda(preStepKinEnergy); }
>> 670
>> 671 // zero cross section
>> 672 if(preStepLambda <= 0.0) {
>> 673 theNumberOfInteractionLengthLeft = -1.0;
>> 674 currentInteractionLength = DBL_MAX;
>> 675 }
>> 676 }
395 677
396 ComputeIntegralLambda(preStepKinEnergy, trac << 678 // non-zero cross section
397 << 679 if(preStepLambda > 0.0) {
398 // zero cross section <<
399 if(preStepLambda <= 0.0) { <<
400 theNumberOfInteractionLengthLeft = -1.0; <<
401 currentInteractionLength = DBL_MAX; <<
402 <<
403 } else { <<
404 680
405 // non-zero cross section <<
406 if (theNumberOfInteractionLengthLeft < 0.0 681 if (theNumberOfInteractionLengthLeft < 0.0) {
407 682
408 // beggining of tracking (or just after 683 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Lo << 684 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
410 theInitialNumberOfInteractionLength = th 685 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
411 686
412 } else { << 687 } else if(currentInteractionLength < DBL_MAX) {
413 688
414 theNumberOfInteractionLengthLeft -= 689 theNumberOfInteractionLengthLeft -=
415 previousStepSize/currentInteractionLen 690 previousStepSize/currentInteractionLength;
416 theNumberOfInteractionLengthLeft = 691 theNumberOfInteractionLengthLeft =
417 std::max(theNumberOfInteractionLengthL 692 std::max(theNumberOfInteractionLengthLeft, 0.0);
418 } 693 }
419 694
420 // new mean free path and step limit for t 695 // new mean free path and step limit for the next step
421 currentInteractionLength = 1.0/preStepLamb 696 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * cur 697 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 } << 698 /*
424 return x; << 699 #ifdef G4VERBOSE
425 } << 700 if (verboseLevel>2){
426 << 701 G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
427 //....oooOO0OOooo........oooOO0OOooo........oo << 702 G4cout << "[ " << GetProcessName() << "]" << G4endl;
428 << 703 G4cout << " for " << currentParticle->GetParticleName()
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 704 << " in Material " << currentMaterial->GetName()
430 { << 705 << " Ekin(MeV)= " << preStepKinEnergy/MeV
431 if (fXSType == fEmNoIntegral) { << 706 <<G4endl;
432 preStepLambda = GetCurrentLambda(e, LogEki << 707 G4cout << " MeanFreePath = " << currentInteractionLength/cm << "[cm]"
433 << 708 << " InteractionLength= " << x/cm <<"[cm] " <<G4endl;
434 } else if (fXSType == fEmIncreasing) { <<
435 if(e*invLambdaFactor < mfpKinEnergy) { <<
436 preStepLambda = GetCurrentLambda(e, LogE <<
437 mfpKinEnergy = (preStepLambda > 0.0) ? e <<
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 } 709 }
459 } else { << 710 #endif
460 preStepLambda = GetCurrentLambda(e, LogEki << 711 */
461 } 712 }
>> 713 return x;
462 } 714 }
463 715
464 //....oooOO0OOooo........oooOO0OOooo........oo 716 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 717
466 G4VParticleChange* G4VEmProcess::PostStepDoIt( 718 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 719 const G4Step& step)
468 { 720 {
469 // clear number of interaction lengths in an << 721 // In all cases clear number of interaction lengths
470 theNumberOfInteractionLengthLeft = -1.0; 722 theNumberOfInteractionLengthLeft = -1.0;
471 mfpKinEnergy = DBL_MAX; << 723 mfpKinEnergy = DBL_MAX;
472 724
473 fParticleChange.InitializeForPostStep(track) 725 fParticleChange.InitializeForPostStep(track);
474 726
475 // Do not make anything if particle is stopp 727 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt! 728 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) 729 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478 730
479 const G4double finalT = track.GetKineticEner << 731 G4double finalT = track.GetKineticEnergy();
480 732
481 // forced process - should happen only once 733 // forced process - should happen only once per track
482 if(biasFlag) { 734 if(biasFlag) {
483 if(biasManager->ForcedInteractionRegion((G << 735 if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
484 biasFlag = false; 736 biasFlag = false;
485 } 737 }
486 } 738 }
487 739
488 // check active and select model <<
489 const G4double scaledEnergy = finalT*massRat <<
490 SelectModel(scaledEnergy, currentCoupleIndex <<
491 if(!currentModel->IsActive(scaledEnergy)) { <<
492 <<
493 // Integral approach 740 // Integral approach
494 if (fXSType != fEmNoIntegral) { << 741 if (integral) {
495 const G4double logFinalT = << 742 G4double lx = GetLambda(finalT, currentCouple);
496 track.GetDynamicParticle()->GetLogKineti << 743 if(preStepLambda<lx && 1 < verboseLevel) {
497 const G4double lx = std::max(GetCurrentLam <<
498 #ifdef G4VERBOSE <<
499 if(preStepLambda < lx && 1 < verboseLevel) <<
500 G4cout << "WARNING: for " << currentPart 744 G4cout << "WARNING: for " << currentParticle->GetParticleName()
501 << " and " << GetProcessName() << << 745 << " and " << GetProcessName()
502 << " preLambda= " << preStepLambd << 746 << " E(MeV)= " << finalT/MeV
503 << " < " << lx << " (postLambda) << 747 << " preLambda= " << preStepLambda << " < "
>> 748 << lx << " (postLambda) "
>> 749 << G4endl;
504 } 750 }
505 #endif << 751
506 // if false interaction then use new cross << 752 if(preStepLambda*G4UniformRand() > lx) {
507 // if both values are zero - no interactio << 753 ClearNumberOfInteractionLengthLeft();
508 if(preStepLambda*G4UniformRand() >= lx) { <<
509 return &fParticleChange; 754 return &fParticleChange;
510 } 755 }
511 } 756 }
512 757
>> 758 SelectModel(finalT, currentCoupleIndex);
>> 759 if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
>> 760
513 // define new weight for primary and seconda 761 // define new weight for primary and secondaries
514 G4double weight = fParticleChange.GetParentW 762 G4double weight = fParticleChange.GetParentWeight();
515 if(weightFlag) { 763 if(weightFlag) {
516 weight /= biasFactor; 764 weight /= biasFactor;
517 fParticleChange.ProposeWeight(weight); 765 fParticleChange.ProposeWeight(weight);
518 } 766 }
519 767
520 #ifdef G4VERBOSE << 768 /*
521 if(1 < verboseLevel) { << 769 if(0 < verboseLevel) {
522 G4cout << "G4VEmProcess::PostStepDoIt: Sam 770 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV 771 << finalT/MeV
524 << " MeV; model= (" << currentModel 772 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnerg 773 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl; 774 << G4endl;
527 } 775 }
528 #endif << 776 */
529 777
530 // sample secondaries 778 // sample secondaries
531 secParticles.clear(); 779 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticle 780 currentModel->SampleSecondaries(&secParticles,
533 currentCoupl 781 currentCouple,
534 track.GetDyn 782 track.GetDynamicParticle(),
535 (*theCuts)[c 783 (*theCuts)[currentCoupleIndex]);
536 784
537 G4int num0 = (G4int)secParticles.size(); << 785 G4int num0 = secParticles.size();
538 786
539 // splitting or Russian roulette 787 // splitting or Russian roulette
540 if(biasManager) { 788 if(biasManager) {
541 if(biasManager->SecondaryBiasingRegion((G4 << 789 if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
542 G4double eloss = 0.0; 790 G4double eloss = 0.0;
543 weight *= biasManager->ApplySecondaryBia 791 weight *= biasManager->ApplySecondaryBiasing(
544 secParticles, track, currentModel, &fP 792 secParticles, track, currentModel, &fParticleChange, eloss,
545 (G4int)currentCoupleIndex, (*theCuts)[ << 793 currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546 step.GetPostStepPoint()->GetSafety()); 794 step.GetPostStepPoint()->GetSafety());
547 if(eloss > 0.0) { 795 if(eloss > 0.0) {
548 eloss += fParticleChange.GetLocalEnerg 796 eloss += fParticleChange.GetLocalEnergyDeposit();
549 fParticleChange.ProposeLocalEnergyDepo 797 fParticleChange.ProposeLocalEnergyDeposit(eloss);
550 } 798 }
551 } 799 }
552 } 800 }
553 801
554 // save secondaries 802 // save secondaries
555 G4int num = (G4int)secParticles.size(); << 803 G4int num = secParticles.size();
556 if(num > 0) { 804 if(num > 0) {
557 805
558 fParticleChange.SetNumberOfSecondaries(num 806 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEn 807 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime(); 808 G4double time = track.GetGlobalTime();
561 <<
562 G4int n1(0), n2(0); <<
563 if(num0 > mainSecondaries) { <<
564 currentModel->FillNumberOfSecondaries(n1 <<
565 } <<
566 809
567 for (G4int i=0; i<num; ++i) { 810 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i]; << 811 if (secParticles[i]) {
569 if (nullptr != dp) { << 812 G4DynamicParticle* dp = secParticles[i];
570 const G4ParticleDefinition* p = dp->Ge 813 const G4ParticleDefinition* p = dp->GetParticleDefinition();
571 G4double e = dp->GetKineticEnergy(); 814 G4double e = dp->GetKineticEnergy();
572 G4bool good = true; 815 G4bool good = true;
573 if(applyCuts) { 816 if(applyCuts) {
574 if (p == theGamma) { 817 if (p == theGamma) {
575 if (e < (*theCutsGamma)[currentCou 818 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576 819
577 } else if (p == theElectron) { 820 } else if (p == theElectron) {
578 if (e < (*theCutsElectron)[current 821 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579 822
580 } else if (p == thePositron) { 823 } else if (p == thePositron) {
581 if (electron_mass_c2 < (*theCutsGa 824 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582 e < (*theCutsPositron)[current 825 e < (*theCutsPositron)[currentCoupleIndex]) {
583 good = false; 826 good = false;
584 e += 2.0*electron_mass_c2; 827 e += 2.0*electron_mass_c2;
585 } 828 }
586 } 829 }
587 // added secondary if it is good 830 // added secondary if it is good
588 } 831 }
589 if (good) { 832 if (good) {
590 G4Track* t = new G4Track(dp, time, t 833 G4Track* t = new G4Track(dp, time, track.GetPosition());
591 t->SetTouchableHandle(track.GetTouch 834 t->SetTouchableHandle(track.GetTouchableHandle());
592 if (biasManager) { << 835 t->SetWeight(weight);
593 t->SetWeight(weight * biasManager- <<
594 } else { <<
595 t->SetWeight(weight); <<
596 } <<
597 pParticleChange->AddSecondary(t); 836 pParticleChange->AddSecondary(t);
598 837
599 // define type of secondary 838 // define type of secondary
600 if(i < mainSecondaries) { << 839 if(i < mainSecondaries) { t->SetCreatorModelIndex(secID); }
601 t->SetCreatorModelID(secID); << 840 else if(i < num0) {
602 if(GetProcessSubType() == fCompton << 841 if(p == theGamma) {
603 t->SetCreatorModelID(_ComptonGam << 842 t->SetCreatorModelIndex(fluoID);
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 { 843 } else {
617 t->SetCreatorModelID(biasID); << 844 t->SetCreatorModelIndex(augerID);
618 } 845 }
>> 846 } else {
>> 847 t->SetCreatorModelIndex(biasID);
619 } 848 }
620 /* << 849
621 G4cout << "Secondary(post step) has << 850 //G4cout << "Secondary(post step) has weight " << t->GetWeight()
622 << ", Ekin= " << t->GetKineti << 851 // << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
623 << GetProcessName() << " fluo <<
624 << " augerID= " << augerID << <<
625 */ <<
626 } else { 852 } else {
627 delete dp; 853 delete dp;
628 edep += e; 854 edep += e;
629 } 855 }
630 } 856 }
631 } 857 }
632 fParticleChange.ProposeLocalEnergyDeposit( 858 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 } 859 }
634 860
635 if(0.0 == fParticleChange.GetProposedKinetic 861 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636 fAlive == fParticleChange.GetTrackStatus( 862 fAlive == fParticleChange.GetTrackStatus()) {
637 if(particle->GetProcessManager()->GetAtRes 863 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638 { fParticleChange.ProposeTrackStatus( 864 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639 else { fParticleChange.ProposeTrackStatus( 865 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640 } 866 }
641 867
642 return &fParticleChange; 868 return &fParticleChange;
643 } 869 }
644 870
645 //....oooOO0OOooo........oooOO0OOooo........oo 871 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646 872
647 G4bool G4VEmProcess::StorePhysicsTable(const G 873 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G << 874 const G4String& directory,
649 G4bool << 875 G4bool ascii)
650 { 876 {
651 if(!isTheMaster || part != particle) { retur << 877 G4bool yes = true;
652 if(G4EmTableUtil::StoreTable(this, part, the << 878 const G4VEmProcess* masterProc =
653 directory, "Lambda", << 879 static_cast<const G4VEmProcess*>(GetMasterProcess());
654 verboseLevel, a << 880 if(masterProc && masterProc != this) { return yes; }
655 G4EmTableUtil::StoreTable(this, part, the << 881
656 directory, "LambdaPrim", << 882 if ( theLambdaTable && part == particle) {
657 verboseLevel, a << 883 const G4String name =
658 return true; << 884 GetPhysicsTableFileName(part,directory,"Lambda",ascii);
>> 885 yes = theLambdaTable->StorePhysicsTable(name,ascii);
>> 886
>> 887 if ( yes ) {
>> 888 G4cout << "Physics table is stored for " << particle->GetParticleName()
>> 889 << " and process " << GetProcessName()
>> 890 << " in the directory <" << directory
>> 891 << "> " << G4endl;
>> 892 } else {
>> 893 G4cout << "Fail to store Physics Table for "
>> 894 << particle->GetParticleName()
>> 895 << " and process " << GetProcessName()
>> 896 << " in the directory <" << directory
>> 897 << "> " << G4endl;
>> 898 }
>> 899 }
>> 900 if ( theLambdaTablePrim && part == particle) {
>> 901 const G4String name =
>> 902 GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
>> 903 yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
>> 904
>> 905 if ( yes ) {
>> 906 G4cout << "Physics table prim is stored for "
>> 907 << particle->GetParticleName()
>> 908 << " and process " << GetProcessName()
>> 909 << " in the directory <" << directory
>> 910 << "> " << G4endl;
>> 911 } else {
>> 912 G4cout << "Fail to store Physics Table Prim for "
>> 913 << particle->GetParticleName()
>> 914 << " and process " << GetProcessName()
>> 915 << " in the directory <" << directory
>> 916 << "> " << G4endl;
>> 917 }
659 } 918 }
660 return false; << 919 return yes;
661 } 920 }
662 921
663 //....oooOO0OOooo........oooOO0OOooo........oo 922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664 923
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons 924 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 cons << 925 const G4String& directory,
667 G4bo 926 G4bool ascii)
668 { 927 {
669 if(!isTheMaster || part != particle) { retur << 928 if(1 < verboseLevel) {
>> 929 G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
>> 930 << part->GetParticleName() << " and process "
>> 931 << GetProcessName() << G4endl;
>> 932 }
670 G4bool yes = true; 933 G4bool yes = true;
>> 934
>> 935 if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy)
>> 936 || particle != part) { return yes; }
>> 937
>> 938 const G4String particleName = part->GetParticleName();
>> 939 G4String filename;
>> 940
671 if(buildLambdaTable) { 941 if(buildLambdaTable) {
672 yes = G4EmTableUtil::RetrieveTable(this, p << 942 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
673 "Lambda << 943 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
674 ascii, << 944 filename,ascii);
675 } << 945 if ( yes ) {
676 if(yes && minKinEnergyPrim < maxKinEnergy) { << 946 if (0 < verboseLevel) {
677 yes = G4EmTableUtil::RetrieveTable(this, p << 947 G4cout << "Lambda table for " << particleName
678 "Lambda << 948 << " is Retrieved from <"
679 ascii, << 949 << filename << ">"
>> 950 << G4endl;
>> 951 }
>> 952 if(theParameters->Spline()) {
>> 953 size_t n = theLambdaTable->length();
>> 954 for(size_t i=0; i<n; ++i) {
>> 955 if((* theLambdaTable)[i]) {
>> 956 (* theLambdaTable)[i]->SetSpline(true);
>> 957 }
>> 958 }
>> 959 }
>> 960 } else {
>> 961 if (1 < verboseLevel) {
>> 962 G4cout << "Lambda table for " << particleName << " in file <"
>> 963 << filename << "> is not exist"
>> 964 << G4endl;
>> 965 }
>> 966 }
>> 967 }
>> 968 if(minKinEnergyPrim < maxKinEnergy) {
>> 969 filename = GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
>> 970 yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTablePrim,
>> 971 filename,ascii);
>> 972 if ( yes ) {
>> 973 if (0 < verboseLevel) {
>> 974 G4cout << "Lambda table prim for " << particleName
>> 975 << " is Retrieved from <"
>> 976 << filename << ">"
>> 977 << G4endl;
>> 978 }
>> 979 if(theParameters->Spline()) {
>> 980 size_t n = theLambdaTablePrim->length();
>> 981 for(size_t i=0; i<n; ++i) {
>> 982 if((* theLambdaTablePrim)[i]) {
>> 983 (* theLambdaTablePrim)[i]->SetSpline(true);
>> 984 }
>> 985 }
>> 986 }
>> 987 } else {
>> 988 if (1 < verboseLevel) {
>> 989 G4cout << "Lambda table prim for " << particleName << " in file <"
>> 990 << filename << "> is not exist"
>> 991 << G4endl;
>> 992 }
>> 993 }
680 } 994 }
>> 995
681 return yes; 996 return yes;
682 } 997 }
683 998
684 //....oooOO0OOooo........oooOO0OOooo........oo 999 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685 1000
686 G4double G4VEmProcess::GetCrossSection(G4doubl << 1001 G4double
687 const G << 1002 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
>> 1003 const G4MaterialCutsCouple* couple)
688 { 1004 {
689 CurrentSetup(couple, kinEnergy); << 1005 // Cross section per atom is calculated
690 return GetCurrentLambda(kinEnergy, G4Log(kin << 1006 DefineMaterial(couple);
>> 1007 G4double cross = 0.0;
>> 1008 if(buildLambdaTable && theLambdaTable) {
>> 1009 cross = GetCurrentLambda(kineticEnergy);
>> 1010 } else {
>> 1011 SelectModel(kineticEnergy, currentCoupleIndex);
>> 1012 cross = fFactor*currentModel->CrossSectionPerVolume(currentMaterial,
>> 1013 currentParticle,
>> 1014 kineticEnergy);
>> 1015 }
>> 1016
>> 1017 if(cross < 0.0) { cross = 0.0; }
>> 1018 return cross;
691 } 1019 }
692 1020
693 //....oooOO0OOooo........oooOO0OOooo........oo 1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694 1022
695 G4double G4VEmProcess::GetMeanFreePath(const G 1023 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696 G4doubl 1024 G4double,
697 G4Force 1025 G4ForceCondition* condition)
698 { 1026 {
699 *condition = NotForced; 1027 *condition = NotForced;
700 return G4VEmProcess::MeanFreePath(track); 1028 return G4VEmProcess::MeanFreePath(track);
701 } 1029 }
702 1030
703 //....oooOO0OOooo........oooOO0OOooo........oo 1031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 1032
705 G4double << 1033 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou <<
707 G4dou <<
708 { 1034 {
709 SelectModel(kinEnergy, currentCoupleIndex); << 1035 DefineMaterial(track.GetMaterialCutsCouple());
710 return (currentModel) ? << 1036 preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
711 currentModel->ComputeCrossSectionPerAtom(c << 1037 G4double x = DBL_MAX;
712 Z << 1038 if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
>> 1039 return x;
713 } 1040 }
714 1041
715 //....oooOO0OOooo........oooOO0OOooo........oo 1042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 1043
717 G4PhysicsVector* << 1044 G4double
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 1045 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy,
>> 1046 G4double Z, G4double A, G4double cut)
719 { 1047 {
720 DefineMaterial(couple); << 1048 SelectModel(kineticEnergy, currentCoupleIndex);
721 G4PhysicsVector* newv = new G4PhysicsLogVect << 1049 G4double x = 0.0;
722 << 1050 if(currentModel) {
723 return newv; << 1051 x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
>> 1052 Z,A,cut);
>> 1053 }
>> 1054 return x;
724 } 1055 }
725 1056
726 //....oooOO0OOooo........oooOO0OOooo........oo 1057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
727 1058
728 const G4Element* G4VEmProcess::GetCurrentEleme << 1059 void G4VEmProcess::FindLambdaMax()
729 { 1060 {
730 return (nullptr != currentModel) ? << 1061 if(1 < verboseLevel) {
731 currentModel->GetCurrentElement(currentMat << 1062 G4cout << "### G4VEmProcess::FindLambdaMax: "
>> 1063 << particle->GetParticleName()
>> 1064 << " and process " << GetProcessName() << " " << G4endl;
>> 1065 }
>> 1066 size_t n = theLambdaTable->length();
>> 1067 G4PhysicsVector* pv;
>> 1068 G4double e, ss, emax, smax;
>> 1069
>> 1070 size_t i;
>> 1071
>> 1072 // first loop on existing vectors
>> 1073 for (i=0; i<n; ++i) {
>> 1074 pv = (*theLambdaTable)[i];
>> 1075 if(pv) {
>> 1076 size_t nb = pv->GetVectorLength();
>> 1077 emax = DBL_MAX;
>> 1078 smax = 0.0;
>> 1079 if(nb > 0) {
>> 1080 for (size_t j=0; j<nb; ++j) {
>> 1081 e = pv->Energy(j);
>> 1082 ss = (*pv)(j);
>> 1083 if(ss > smax) {
>> 1084 smax = ss;
>> 1085 emax = e;
>> 1086 }
>> 1087 }
>> 1088 }
>> 1089 theEnergyOfCrossSectionMax[i] = emax;
>> 1090 theCrossSectionMax[i] = smax;
>> 1091 if(1 < verboseLevel) {
>> 1092 G4cout << "For " << particle->GetParticleName()
>> 1093 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
>> 1094 << " lambda= " << smax << G4endl;
>> 1095 }
>> 1096 }
>> 1097 }
>> 1098 // second loop using base materials
>> 1099 for (i=0; i<n; ++i) {
>> 1100 pv = (*theLambdaTable)[i];
>> 1101 if(!pv){
>> 1102 G4int j = (*theDensityIdx)[i];
>> 1103 theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
>> 1104 theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
>> 1105 }
>> 1106 }
732 } 1107 }
733 1108
734 //....oooOO0OOooo........oooOO0OOooo........oo 1109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735 1110
736 const G4Element* G4VEmProcess::GetTargetElemen << 1111 G4PhysicsVector*
>> 1112 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
737 { 1113 {
738 return (nullptr != currentModel) ? << 1114 G4PhysicsVector* v =
739 currentModel->GetCurrentElement(currentMat << 1115 new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
>> 1116 v->SetSpline(theParameters->Spline());
>> 1117 return v;
740 } 1118 }
741 1119
742 //....oooOO0OOooo........oooOO0OOooo........oo 1120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743 1121
744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 1122 const G4Element* G4VEmProcess::GetCurrentElement() const
745 { 1123 {
746 return (nullptr != currentModel) ? << 1124 const G4Element* elm = nullptr;
747 currentModel->GetCurrentIsotope(GetCurrent << 1125 if(currentModel) {elm = currentModel->GetCurrentElement(); }
>> 1126 return elm;
748 } 1127 }
749 1128
750 //....oooOO0OOooo........oooOO0OOooo........oo 1129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751 1130
752 void G4VEmProcess::SetCrossSectionBiasingFacto 1131 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 { 1132 {
754 if(f > 0.0) { 1133 if(f > 0.0) {
755 biasFactor = f; 1134 biasFactor = f;
756 weightFlag = flag; 1135 weightFlag = flag;
757 if(1 < verboseLevel) { 1136 if(1 < verboseLevel) {
758 G4cout << "### SetCrossSectionBiasingFac 1137 G4cout << "### SetCrossSectionBiasingFactor: for "
759 << particle->GetParticleName() 1138 << particle->GetParticleName()
760 << " and process " << GetProcessN 1139 << " and process " << GetProcessName()
761 << " biasFactor= " << f << " weig 1140 << " biasFactor= " << f << " weightFlag= " << flag
762 << G4endl; 1141 << G4endl;
763 } 1142 }
764 } 1143 }
765 } 1144 }
766 1145
767 //....oooOO0OOooo........oooOO0OOooo........oo 1146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768 1147
769 void 1148 void
770 G4VEmProcess::ActivateForcedInteraction(G4doub 1149 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771 G4bool 1150 G4bool flag)
772 { 1151 {
773 if(nullptr == biasManager) { biasManager = n << 1152 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
774 if(1 < verboseLevel) { 1153 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: 1154 G4cout << "### ActivateForcedInteraction: for "
776 << particle->GetParticleName() 1155 << particle->GetParticleName()
777 << " and process " << GetProcessNam 1156 << " and process " << GetProcessName()
778 << " length(mm)= " << length/mm 1157 << " length(mm)= " << length/mm
779 << " in G4Region <" << r 1158 << " in G4Region <" << r
780 << "> weightFlag= " << flag 1159 << "> weightFlag= " << flag
781 << G4endl; 1160 << G4endl;
782 } 1161 }
783 weightFlag = flag; 1162 weightFlag = flag;
784 biasManager->ActivateForcedInteraction(lengt 1163 biasManager->ActivateForcedInteraction(length, r);
785 } 1164 }
786 1165
787 //....oooOO0OOooo........oooOO0OOooo........oo 1166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788 1167
789 void 1168 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G 1169 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791 G4double factor, 1170 G4double factor,
792 G4double energyLimit) 1171 G4double energyLimit)
793 { 1172 {
794 if (0.0 <= factor) { 1173 if (0.0 <= factor) {
795 1174
796 // Range cut can be applied only for e- 1175 // Range cut can be applied only for e-
797 if(0.0 == factor && secondaryParticle != G 1176 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798 { return; } 1177 { return; }
799 1178
800 if(!biasManager) { biasManager = new G4EmB 1179 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801 biasManager->ActivateSecondaryBiasing(regi 1180 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802 if(1 < verboseLevel) { 1181 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: 1182 G4cout << "### ActivateSecondaryBiasing: for "
804 << " process " << GetProcessName() 1183 << " process " << GetProcessName()
805 << " factor= " << factor 1184 << " factor= " << factor
806 << " in G4Region <" << region 1185 << " in G4Region <" << region
807 << "> energyLimit(MeV)= " << energyLimi 1186 << "> energyLimit(MeV)= " << energyLimit/MeV
808 << G4endl; 1187 << G4endl;
809 } 1188 }
810 } 1189 }
811 } 1190 }
812 1191
813 //....oooOO0OOooo........oooOO0OOooo........oo 1192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 1193
815 void G4VEmProcess::SetLambdaBinning(G4int n) 1194 void G4VEmProcess::SetLambdaBinning(G4int n)
816 { 1195 {
817 if(5 < n && n < 10000000) { 1196 if(5 < n && n < 10000000) {
818 nLambdaBins = n; 1197 nLambdaBins = n;
819 actBinning = true; 1198 actBinning = true;
820 } else { 1199 } else {
821 G4double e = (G4double)n; 1200 G4double e = (G4double)n;
822 PrintWarning("SetLambdaBinning", e); 1201 PrintWarning("SetLambdaBinning", e);
823 } 1202 }
824 } 1203 }
825 1204
826 //....oooOO0OOooo........oooOO0OOooo........oo 1205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827 1206
828 void G4VEmProcess::SetMinKinEnergy(G4double e) 1207 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 { 1208 {
830 if(1.e-3*eV < e && e < maxKinEnergy) { 1209 if(1.e-3*eV < e && e < maxKinEnergy) {
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(ma 1210 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e)
832 /G4Log(maxKinEnergy/ 1211 /G4Log(maxKinEnergy/minKinEnergy));
833 minKinEnergy = e; 1212 minKinEnergy = e;
834 actMinKinEnergy = true; 1213 actMinKinEnergy = true;
835 } else { PrintWarning("SetMinKinEnergy", e); 1214 } else { PrintWarning("SetMinKinEnergy", e); }
836 } 1215 }
837 1216
838 //....oooOO0OOooo........oooOO0OOooo........oo 1217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839 1218
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) 1219 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 { 1220 {
842 if(minKinEnergy < e && e < 1.e+6*TeV) { 1221 if(minKinEnergy < e && e < 1.e+6*TeV) {
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ 1222 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy)
844 /G4Log(maxKinEnergy/ 1223 /G4Log(maxKinEnergy/minKinEnergy));
845 maxKinEnergy = e; 1224 maxKinEnergy = e;
846 actMaxKinEnergy = true; 1225 actMaxKinEnergy = true;
847 } else { PrintWarning("SetMaxKinEnergy", e); 1226 } else { PrintWarning("SetMaxKinEnergy", e); }
848 } 1227 }
849 1228
850 //....oooOO0OOooo........oooOO0OOooo........oo 1229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851 1230
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl 1231 void G4VEmProcess::SetMinKinEnergyPrim(G4double e)
853 { 1232 {
854 if(theParameters->MinKinEnergy() <= e && 1233 if(theParameters->MinKinEnergy() <= e &&
855 e <= theParameters->MaxKinEnergy()) { min 1234 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; }
856 else { PrintWarning("SetMinKinEnergyPrim", e 1235 else { PrintWarning("SetMinKinEnergyPrim", e); }
857 } 1236 }
858 1237
859 //....oooOO0OOooo........oooOO0OOooo........oo 1238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
860 1239
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, 1240 void G4VEmProcess::PrintWarning(G4String tit, G4double val)
876 { 1241 {
877 G4String ss = "G4VEmProcess::" + tit; 1242 G4String ss = "G4VEmProcess::" + tit;
878 G4ExceptionDescription ed; 1243 G4ExceptionDescription ed;
879 ed << "Parameter is out of range: " << val 1244 ed << "Parameter is out of range: " << val
880 << " it will have no effect!\n" << " Pro 1245 << " it will have no effect!\n" << " Process "
881 << GetProcessName() << " nbins= " << the 1246 << GetProcessName() << " nbins= " << theParameters->NumberOfBins()
882 << " Emin(keV)= " << theParameters->MinKi 1247 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV
883 << " Emax(GeV)= " << theParameters->MaxKi 1248 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV;
884 G4Exception(ss, "em0044", JustWarning, ed); 1249 G4Exception(ss, "em0044", JustWarning, ed);
885 } 1250 }
886 1251
887 //....oooOO0OOooo........oooOO0OOooo........oo 1252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 1253
889 void G4VEmProcess::ProcessDescription(std::ost << 1254 void G4VEmProcess::ProcessDescription(std::ostream& outFile) const
890 { 1255 {
891 if(nullptr != particle) { << 1256 outFile << "EM process <" << GetProcessName() << ">" << G4endl;
892 StreamInfo(out, *particle, true); <<
893 } <<
894 } 1257 }
895 1258
896 //....oooOO0OOooo........oooOO0OOooo........oo 1259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 1260