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