| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4VEmProcess
32 //
33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 //
35 // Creation date: 01.10.2003
36 //
37 // Modifications: by V.Ivanchenko
38 //
39 // Class Description: based class for discrete and rest/discrete EM processes
40 //
41
42 // -------------------------------------------------------------------
43 //
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
47 #include "G4VEmProcess.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4ProcessManager.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4LossTableBuilder.hh"
53 #include "G4Step.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4VEmModel.hh"
56 #include "G4DataVector.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4EmDataHandler.hh"
59 #include "G4PhysicsLogVector.hh"
60 #include "G4VParticleChange.hh"
61 #include "G4ProductionCutsTable.hh"
62 #include "G4Region.hh"
63 #include "G4Gamma.hh"
64 #include "G4Electron.hh"
65 #include "G4Positron.hh"
66 #include "G4PhysicsTableHelper.hh"
67 #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"
74 #include "G4Log.hh"
75 #include <iostream>
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
80 G4VDiscreteProcess(name, type)
81 {
82 theParameters = G4EmParameters::Instance();
83 SetVerboseLevel(1);
84
85 // Size of tables
86 minKinEnergy = 0.1*CLHEP::keV;
87 maxKinEnergy = 100.0*CLHEP::TeV;
88
89 // default lambda factor
90 invLambdaFactor = 1.0/lambdaFactor;
91
92 // particle types
93 theGamma = G4Gamma::Gamma();
94 theElectron = G4Electron::Electron();
95 thePositron = G4Positron::Positron();
96
97 pParticleChange = &fParticleChange;
98 fParticleChange.SetSecondaryWeightByProcess(true);
99 secParticles.reserve(5);
100
101 modelManager = new G4EmModelManager();
102 lManager = G4LossTableManager::Instance();
103 lManager->Register(this);
104 isTheMaster = lManager->IsMaster();
105 G4LossTableBuilder* bld = lManager->GetTableBuilder();
106 theDensityFactor = bld->GetDensityFactors();
107 theDensityIdx = bld->GetCoupleIndexes();
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
112 G4VEmProcess::~G4VEmProcess()
113 {
114 if(isTheMaster) {
115 delete theData;
116 delete theEnergyOfCrossSectionMax;
117 }
118 delete modelManager;
119 delete biasManager;
120 lManager->DeRegister(this);
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* ptr,
126 const G4Region* region)
127 {
128 if(nullptr == ptr) { return; }
129 G4VEmFluctuationModel* fm = nullptr;
130 modelManager->AddEmModel(order, ptr, fm, region);
131 ptr->SetParticleChange(pParticleChange);
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, G4int)
137 {
138 if(nullptr == ptr) { return; }
139 if(!emModels.empty()) {
140 for(auto & em : emModels) { if(em == ptr) { return; } }
141 }
142 emModels.push_back(ptr);
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146
147 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
148 {
149 if(nullptr == particle) { SetParticle(&part); }
150
151 if(part.GetParticleType() == "nucleus" &&
152 part.GetParticleSubType() == "generic") {
153
154 G4String pname = part.GetParticleName();
155 if(pname != "deuteron" && pname != "triton" &&
156 pname != "He3" && pname != "alpha" && pname != "alpha+" &&
157 pname != "helium" && pname != "hydrogen") {
158
159 particle = G4GenericIon::GenericIon();
160 isIon = true;
161 }
162 }
163 if(particle != &part) { return; }
164
165 lManager->PreparePhysicsTable(&part, this);
166
167 // for new run
168 currentCouple = nullptr;
169 preStepLambda = 0.0;
170 fLambdaEnergy = 0.0;
171
172 InitialiseProcess(particle);
173
174 G4LossTableBuilder* bld = lManager->GetTableBuilder();
175 const G4ProductionCutsTable* theCoupleTable=
176 G4ProductionCutsTable::GetProductionCutsTable();
177 theCutsGamma = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
178 theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
179 theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
180
181 // initialisation of the process
182 if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
183 if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
184
185 applyCuts = theParameters->ApplyCuts();
186 lambdaFactor = theParameters->LambdaFactor();
187 invLambdaFactor = 1.0/lambdaFactor;
188 theParameters->DefineRegParamForEM(this);
189
190 // integral option may be disabled
191 if(!theParameters->Integral()) { fXSType = fEmNoIntegral; }
192
193 // prepare tables
194 if(isTheMaster) {
195 if(nullptr == theData) { theData = new G4EmDataHandler(2); }
196
197 if(buildLambdaTable) {
198 theLambdaTable = theData->MakeTable(0);
199 bld->InitialiseBaseMaterials(theLambdaTable);
200 }
201 // high energy table
202 if(minKinEnergyPrim < maxKinEnergy) {
203 theLambdaTablePrim = theData->MakeTable(1);
204 bld->InitialiseBaseMaterials(theLambdaTablePrim);
205 }
206 }
207 // models
208 baseMat = bld->GetBaseMaterialFlag();
209 numberOfModels = modelManager->NumberOfModels();
210 currentModel = modelManager->GetModel(0);
211 if(nullptr != lManager->AtomDeexcitation()) {
212 modelManager->SetFluoFlag(true);
213 }
214 // forced biasing
215 if(nullptr != biasManager) {
216 biasManager->Initialise(part, GetProcessName(), verboseLevel);
217 biasFlag = false;
218 }
219
220 theCuts =
221 G4EmTableUtil::PrepareEmProcess(this, particle, secondaryParticle,
222 modelManager, maxKinEnergy,
223 secID, tripletID, mainSecondaries,
224 verboseLevel, isTheMaster);
225 }
226
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228
229 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 {
231 if(nullptr == masterProc) {
232 if(isTheMaster) { masterProc = this; }
233 else { masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess());}
234 }
235 G4int nModels = modelManager->NumberOfModels();
236 G4bool isLocked = theParameters->IsPrintLocked();
237 G4bool toBuild = (buildLambdaTable || minKinEnergyPrim < maxKinEnergy);
238
239 G4EmTableUtil::BuildEmProcess(this, masterProc, particle, &part,
240 nModels, verboseLevel, isTheMaster,
241 isLocked, toBuild, baseMat);
242 }
243
244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245
246 void G4VEmProcess::BuildLambdaTable()
247 {
248 G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
249 G4int nbin =
250 theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale));
251 if(actBinning) { nbin = std::max(nbin, nLambdaBins); }
252 scale = nbin/G4Log(scale);
253
254 G4LossTableBuilder* bld = lManager->GetTableBuilder();
255 G4EmTableUtil::BuildLambdaTable(this, particle, modelManager,
256 bld, theLambdaTable, theLambdaTablePrim,
257 minKinEnergy, minKinEnergyPrim,
258 maxKinEnergy, scale, verboseLevel,
259 startFromNull, splineFlag);
260 }
261
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263
264 void G4VEmProcess::StreamInfo(std::ostream& out,
265 const G4ParticleDefinition& part, G4bool rst) const
266 {
267 G4String indent = (rst ? " " : "");
268 out << std::setprecision(6);
269 out << G4endl << indent << GetProcessName() << ": ";
270 if (!rst) {
271 out << " for " << part.GetParticleName();
272 }
273 if(fXSType != fEmNoIntegral) { out << " XStype:" << fXSType; }
274 if(applyCuts) { out << " applyCuts:1 "; }
275 G4int subtype = GetProcessSubType();
276 out << " SubType=" << subtype;
277 if (subtype == fAnnihilation) {
278 G4int mod = theParameters->PositronAtRestModelType();
279 const G4String namp[2] = {"Simple", "Allison"};
280 out << " AtRestModel:" << namp[mod];
281 }
282 if(biasFactor != 1.0) { out << " BiasingFactor=" << biasFactor; }
283 out << " BuildTable=" << buildLambdaTable << G4endl;
284 if(buildLambdaTable) {
285 if(particle == &part) {
286 for(auto & v : *theLambdaTable) {
287 if(nullptr != v) {
288 out << " Lambda table from ";
289 G4double emin = v->Energy(0);
290 G4double emax = v->GetMaxEnergy();
291 G4int nbin = G4int(v->GetVectorLength() - 1);
292 if(emin > minKinEnergy) { out << "threshold "; }
293 else { out << G4BestUnit(emin,"Energy"); }
294 out << " to "
295 << G4BestUnit(emax,"Energy")
296 << ", " << G4lrint(nbin/std::log10(emax/emin))
297 << " bins/decade, spline: "
298 << splineFlag << G4endl;
299 break;
300 }
301 }
302 } else {
303 out << " Used Lambda table of "
304 << particle->GetParticleName() << G4endl;
305 }
306 }
307 if(minKinEnergyPrim < maxKinEnergy) {
308 if(particle == &part) {
309 for(auto & v : *theLambdaTablePrim) {
310 if(nullptr != v) {
311 out << " LambdaPrime table from "
312 << G4BestUnit(v->Energy(0),"Energy")
313 << " to "
314 << G4BestUnit(v->GetMaxEnergy(),"Energy")
315 << " in " << v->GetVectorLength()-1
316 << " bins " << G4endl;
317 break;
318 }
319 }
320 } else {
321 out << " Used LambdaPrime table of "
322 << particle->GetParticleName() << G4endl;
323 }
324 }
325 StreamProcessInfo(out);
326 modelManager->DumpModelList(out, verboseLevel);
327
328 if(verboseLevel > 2 && buildLambdaTable) {
329 out << " LambdaTable address= " << theLambdaTable << G4endl;
330 if(theLambdaTable && particle == &part) {
331 out << (*theLambdaTable) << G4endl;
332 }
333 }
334 }
335
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337
338 void G4VEmProcess::StartTracking(G4Track* track)
339 {
340 // reset parameters for the new track
341 currentParticle = track->GetParticleDefinition();
342 theNumberOfInteractionLengthLeft = -1.0;
343 mfpKinEnergy = DBL_MAX;
344 preStepLambda = 0.0;
345
346 if(isIon) { massRatio = proton_mass_c2/currentParticle->GetPDGMass(); }
347
348 // forced biasing only for primary particles
349 if(biasManager) {
350 if(0 == track->GetParentID()) {
351 // primary particle
352 biasFlag = true;
353 biasManager->ResetForcedInteraction();
354 }
355 }
356 }
357
358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359
360 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361 const G4Track& track,
362 G4double previousStepSize,
363 G4ForceCondition* condition)
364 {
365 *condition = NotForced;
366 G4double x = DBL_MAX;
367
368 DefineMaterial(track.GetMaterialCutsCouple());
369 preStepKinEnergy = track.GetKineticEnergy();
370 const G4double scaledEnergy = preStepKinEnergy*massRatio;
371 SelectModel(scaledEnergy, currentCoupleIndex);
372 /*
373 G4cout << "PostStepGetPhysicalInteractionLength: idx= " << currentCoupleIndex
374 << " couple: " << currentCouple << G4endl;
375 */
376 if(!currentModel->IsActive(scaledEnergy)) {
377 theNumberOfInteractionLengthLeft = -1.0;
378 currentInteractionLength = DBL_MAX;
379 mfpKinEnergy = DBL_MAX;
380 preStepLambda = 0.0;
381 return x;
382 }
383
384 // forced biasing only for primary particles
385 if(biasManager) {
386 if(0 == track.GetParentID()) {
387 if(biasFlag &&
388 biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
389 return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize);
390 }
391 }
392 }
393
394 // compute mean free path
395
396 ComputeIntegralLambda(preStepKinEnergy, track);
397
398 // zero cross section
399 if(preStepLambda <= 0.0) {
400 theNumberOfInteractionLengthLeft = -1.0;
401 currentInteractionLength = DBL_MAX;
402
403 } else {
404
405 // non-zero cross section
406 if (theNumberOfInteractionLengthLeft < 0.0) {
407
408 // beggining of tracking (or just after DoIt of this process)
409 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
410 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
411
412 } else {
413
414 theNumberOfInteractionLengthLeft -=
415 previousStepSize/currentInteractionLength;
416 theNumberOfInteractionLengthLeft =
417 std::max(theNumberOfInteractionLengthLeft, 0.0);
418 }
419
420 // new mean free path and step limit for the next step
421 currentInteractionLength = 1.0/preStepLambda;
422 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423 }
424 return x;
425 }
426
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
428
429 void G4VEmProcess::ComputeIntegralLambda(G4double e, const G4Track& track)
430 {
431 if (fXSType == fEmNoIntegral) {
432 preStepLambda = GetCurrentLambda(e, LogEkin(track));
433
434 } else if (fXSType == fEmIncreasing) {
435 if(e*invLambdaFactor < mfpKinEnergy) {
436 preStepLambda = GetCurrentLambda(e, LogEkin(track));
437 mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
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 = (*theEnergyOfCrossSectionMax)[currentCoupleIndex];
449 if(e <= epeak) {
450 if(e*invLambdaFactor < mfpKinEnergy) {
451 preStepLambda = GetCurrentLambda(e, LogEkin(track));
452 mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
453 }
454 } else if(e < mfpKinEnergy) {
455 const G4double e1 = std::max(epeak, e*lambdaFactor);
456 preStepLambda = GetCurrentLambda(e1);
457 mfpKinEnergy = e1;
458 }
459 } else {
460 preStepLambda = GetCurrentLambda(e, LogEkin(track));
461 }
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465
466 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467 const G4Step& step)
468 {
469 // clear number of interaction lengths in any case
470 theNumberOfInteractionLengthLeft = -1.0;
471 mfpKinEnergy = DBL_MAX;
472
473 fParticleChange.InitializeForPostStep(track);
474
475 // Do not make anything if particle is stopped, the annihilation then
476 // should be performed by the AtRestDoIt!
477 if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478
479 const G4double finalT = track.GetKineticEnergy();
480
481 // forced process - should happen only once per track
482 if(biasFlag) {
483 if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
484 biasFlag = false;
485 }
486 }
487
488 // check active and select model
489 const G4double scaledEnergy = finalT*massRatio;
490 SelectModel(scaledEnergy, currentCoupleIndex);
491 if(!currentModel->IsActive(scaledEnergy)) { return &fParticleChange; }
492
493 // Integral approach
494 if (fXSType != fEmNoIntegral) {
495 const G4double logFinalT =
496 track.GetDynamicParticle()->GetLogKineticEnergy();
497 const G4double lx = std::max(GetCurrentLambda(finalT, logFinalT), 0.0);
498 #ifdef G4VERBOSE
499 if(preStepLambda < lx && 1 < verboseLevel) {
500 G4cout << "WARNING: for " << currentParticle->GetParticleName()
501 << " and " << GetProcessName() << " E(MeV)= " << finalT/MeV
502 << " preLambda= " << preStepLambda
503 << " < " << lx << " (postLambda) " << G4endl;
504 }
505 #endif
506 // if false interaction then use new cross section value
507 // if both values are zero - no interaction
508 if(preStepLambda*G4UniformRand() >= lx) {
509 return &fParticleChange;
510 }
511 }
512
513 // define new weight for primary and secondaries
514 G4double weight = fParticleChange.GetParentWeight();
515 if(weightFlag) {
516 weight /= biasFactor;
517 fParticleChange.ProposeWeight(weight);
518 }
519
520 #ifdef G4VERBOSE
521 if(1 < verboseLevel) {
522 G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523 << finalT/MeV
524 << " MeV; model= (" << currentModel->LowEnergyLimit()
525 << ", " << currentModel->HighEnergyLimit() << ")"
526 << G4endl;
527 }
528 #endif
529
530 // sample secondaries
531 secParticles.clear();
532 currentModel->SampleSecondaries(&secParticles,
533 currentCouple,
534 track.GetDynamicParticle(),
535 (*theCuts)[currentCoupleIndex]);
536
537 G4int num0 = (G4int)secParticles.size();
538
539 // splitting or Russian roulette
540 if(biasManager) {
541 if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
542 G4double eloss = 0.0;
543 weight *= biasManager->ApplySecondaryBiasing(
544 secParticles, track, currentModel, &fParticleChange, eloss,
545 (G4int)currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546 step.GetPostStepPoint()->GetSafety());
547 if(eloss > 0.0) {
548 eloss += fParticleChange.GetLocalEnergyDeposit();
549 fParticleChange.ProposeLocalEnergyDeposit(eloss);
550 }
551 }
552 }
553
554 // save secondaries
555 G4int num = (G4int)secParticles.size();
556 if(num > 0) {
557
558 fParticleChange.SetNumberOfSecondaries(num);
559 G4double edep = fParticleChange.GetLocalEnergyDeposit();
560 G4double time = track.GetGlobalTime();
561
562 G4int n1(0), n2(0);
563 if(num0 > mainSecondaries) {
564 currentModel->FillNumberOfSecondaries(n1, n2);
565 }
566
567 for (G4int i=0; i<num; ++i) {
568 G4DynamicParticle* dp = secParticles[i];
569 if (nullptr != dp) {
570 const G4ParticleDefinition* p = dp->GetParticleDefinition();
571 G4double e = dp->GetKineticEnergy();
572 G4bool good = true;
573 if(applyCuts) {
574 if (p == theGamma) {
575 if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576
577 } else if (p == theElectron) {
578 if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579
580 } else if (p == thePositron) {
581 if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582 e < (*theCutsPositron)[currentCoupleIndex]) {
583 good = false;
584 e += 2.0*electron_mass_c2;
585 }
586 }
587 // added secondary if it is good
588 }
589 if (good) {
590 G4Track* t = new G4Track(dp, time, track.GetPosition());
591 t->SetTouchableHandle(track.GetTouchableHandle());
592 if (biasManager) {
593 t->SetWeight(weight * biasManager->GetWeight(i));
594 } else {
595 t->SetWeight(weight);
596 }
597 pParticleChange->AddSecondary(t);
598
599 // define type of secondary
600 if(i < mainSecondaries) {
601 t->SetCreatorModelID(secID);
602 if(GetProcessSubType() == fComptonScattering && p == theGamma) {
603 t->SetCreatorModelID(_ComptonGamma);
604 }
605 } else if(i < mainSecondaries + n1) {
606 t->SetCreatorModelID(tripletID);
607 } else if(i < mainSecondaries + n1 + n2) {
608 t->SetCreatorModelID(_IonRecoil);
609 } else {
610 if(i < num0) {
611 if(p == theGamma) {
612 t->SetCreatorModelID(fluoID);
613 } else {
614 t->SetCreatorModelID(augerID);
615 }
616 } else {
617 t->SetCreatorModelID(biasID);
618 }
619 }
620 /*
621 G4cout << "Secondary(post step) has weight " << t->GetWeight()
622 << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV "
623 << GetProcessName() << " fluoID= " << fluoID
624 << " augerID= " << augerID <<G4endl;
625 */
626 } else {
627 delete dp;
628 edep += e;
629 }
630 }
631 }
632 fParticleChange.ProposeLocalEnergyDeposit(edep);
633 }
634
635 if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636 fAlive == fParticleChange.GetTrackStatus()) {
637 if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638 { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639 else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640 }
641
642 return &fParticleChange;
643 }
644
645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646
647 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648 const G4String& directory,
649 G4bool ascii)
650 {
651 if(!isTheMaster || part != particle) { return true; }
652 if(G4EmTableUtil::StoreTable(this, part, theLambdaTable,
653 directory, "Lambda",
654 verboseLevel, ascii) &&
655 G4EmTableUtil::StoreTable(this, part, theLambdaTablePrim,
656 directory, "LambdaPrim",
657 verboseLevel, ascii)) {
658 return true;
659 }
660 return false;
661 }
662
663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664
665 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666 const G4String& dir,
667 G4bool ascii)
668 {
669 if(!isTheMaster || part != particle) { return true; }
670 G4bool yes = true;
671 if(buildLambdaTable) {
672 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTable, dir,
673 "Lambda", verboseLevel,
674 ascii, splineFlag);
675 }
676 if(yes && minKinEnergyPrim < maxKinEnergy) {
677 yes = G4EmTableUtil::RetrieveTable(this, part, theLambdaTablePrim, dir,
678 "LambdaPrim", verboseLevel,
679 ascii, splineFlag);
680 }
681 return yes;
682 }
683
684 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685
686 G4double G4VEmProcess::GetCrossSection(G4double kinEnergy,
687 const G4MaterialCutsCouple* couple)
688 {
689 CurrentSetup(couple, kinEnergy);
690 return GetCurrentLambda(kinEnergy, G4Log(kinEnergy));
691 }
692
693 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694
695 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696 G4double,
697 G4ForceCondition* condition)
698 {
699 *condition = NotForced;
700 return G4VEmProcess::MeanFreePath(track);
701 }
702
703 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704
705 G4double
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kinEnergy,
707 G4double Z, G4double A, G4double cut)
708 {
709 SelectModel(kinEnergy, currentCoupleIndex);
710 return (currentModel) ?
711 currentModel->ComputeCrossSectionPerAtom(currentParticle, kinEnergy,
712 Z, A, cut) : 0.0;
713 }
714
715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716
717 G4PhysicsVector*
718 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple)
719 {
720 DefineMaterial(couple);
721 G4PhysicsVector* newv = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy,
722 nLambdaBins, splineFlag);
723 return newv;
724 }
725
726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
727
728 const G4Element* G4VEmProcess::GetCurrentElement() const
729 {
730 return (nullptr != currentModel) ?
731 currentModel->GetCurrentElement(currentMaterial) : nullptr;
732 }
733
734 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735
736 const G4Element* G4VEmProcess::GetTargetElement() const
737 {
738 return (nullptr != currentModel) ?
739 currentModel->GetCurrentElement(currentMaterial) : nullptr;
740 }
741
742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743
744 const G4Isotope* G4VEmProcess::GetTargetIsotope() const
745 {
746 return (nullptr != currentModel) ?
747 currentModel->GetCurrentIsotope(GetCurrentElement()) : nullptr;
748 }
749
750 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751
752 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 {
754 if(f > 0.0) {
755 biasFactor = f;
756 weightFlag = flag;
757 if(1 < verboseLevel) {
758 G4cout << "### SetCrossSectionBiasingFactor: for "
759 << particle->GetParticleName()
760 << " and process " << GetProcessName()
761 << " biasFactor= " << f << " weightFlag= " << flag
762 << G4endl;
763 }
764 }
765 }
766
767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768
769 void
770 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771 G4bool flag)
772 {
773 if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
774 if(1 < verboseLevel) {
775 G4cout << "### ActivateForcedInteraction: for "
776 << particle->GetParticleName()
777 << " and process " << GetProcessName()
778 << " length(mm)= " << length/mm
779 << " in G4Region <" << r
780 << "> weightFlag= " << flag
781 << G4endl;
782 }
783 weightFlag = flag;
784 biasManager->ActivateForcedInteraction(length, r);
785 }
786
787 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788
789 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791 G4double factor,
792 G4double energyLimit)
793 {
794 if (0.0 <= factor) {
795
796 // Range cut can be applied only for e-
797 if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798 { return; }
799
800 if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801 biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802 if(1 < verboseLevel) {
803 G4cout << "### ActivateSecondaryBiasing: for "
804 << " process " << GetProcessName()
805 << " factor= " << factor
806 << " in G4Region <" << region
807 << "> energyLimit(MeV)= " << energyLimit/MeV
808 << G4endl;
809 }
810 }
811 }
812
813 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814
815 void G4VEmProcess::SetLambdaBinning(G4int n)
816 {
817 if(5 < n && n < 10000000) {
818 nLambdaBins = n;
819 actBinning = true;
820 } else {
821 G4double e = (G4double)n;
822 PrintWarning("SetLambdaBinning", e);
823 }
824 }
825
826 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
827
828 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 {
830 if(1.e-3*eV < e && e < maxKinEnergy) {
831 nLambdaBins = G4lrint(nLambdaBins*G4Log(maxKinEnergy/e)
832 /G4Log(maxKinEnergy/minKinEnergy));
833 minKinEnergy = e;
834 actMinKinEnergy = true;
835 } else { PrintWarning("SetMinKinEnergy", e); }
836 }
837
838 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839
840 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 {
842 if(minKinEnergy < e && e < 1.e+6*TeV) {
843 nLambdaBins = G4lrint(nLambdaBins*G4Log(e/minKinEnergy)
844 /G4Log(maxKinEnergy/minKinEnergy));
845 maxKinEnergy = e;
846 actMaxKinEnergy = true;
847 } else { PrintWarning("SetMaxKinEnergy", e); }
848 }
849
850 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851
852 void G4VEmProcess::SetMinKinEnergyPrim(G4double e)
853 {
854 if(theParameters->MinKinEnergy() <= e &&
855 e <= theParameters->MaxKinEnergy()) { minKinEnergyPrim = e; }
856 else { PrintWarning("SetMinKinEnergyPrim", e); }
857 }
858
859 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
860
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const G4String& nam)
862 {
863 return (nam == GetProcessName()) ? this : nullptr;
864 }
865
866 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
867
868 G4double G4VEmProcess::PolarAngleLimit() const
869 {
870 return theParameters->MscThetaLimit();
871 }
872
873 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
874
875 void G4VEmProcess::PrintWarning(G4String tit, G4double val)
876 {
877 G4String ss = "G4VEmProcess::" + tit;
878 G4ExceptionDescription ed;
879 ed << "Parameter is out of range: " << val
880 << " it will have no effect!\n" << " Process "
881 << GetProcessName() << " nbins= " << theParameters->NumberOfBins()
882 << " Emin(keV)= " << theParameters->MinKinEnergy()/keV
883 << " Emax(GeV)= " << theParameters->MaxKinEnergy()/GeV;
884 G4Exception(ss, "em0044", JustWarning, ed);
885 }
886
887 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888
889 void G4VEmProcess::ProcessDescription(std::ostream& out) const
890 {
891 if(nullptr != particle) {
892 StreamInfo(out, *particle, true);
893 }
894 }
895
896 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897