| 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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4GammaGeneralProcess
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 19.07.2018
37 //
38 // Modifications:
39 //
40 // Class Description:
41 //
42
43 // -------------------------------------------------------------------
44 //
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48 #include "G4GammaGeneralProcess.hh"
49 #include "G4VDiscreteProcess.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4ProcessManager.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4LossTableBuilder.hh"
55 #include "G4HadronicProcess.hh"
56 #include "G4LossTableManager.hh"
57 #include "G4Step.hh"
58 #include "G4Track.hh"
59 #include "G4ParticleDefinition.hh"
60 #include "G4DataVector.hh"
61 #include "G4PhysicsTable.hh"
62 #include "G4PhysicsLogVector.hh"
63 #include "G4PhysicsLinearVector.hh"
64 #include "G4VParticleChange.hh"
65 #include "G4PhysicsTableHelper.hh"
66 #include "G4EmParameters.hh"
67 #include "G4EmProcessSubType.hh"
68 #include "G4Material.hh"
69 #include "G4MaterialCutsCouple.hh"
70 #include "G4GammaConversionToMuons.hh"
71 #include "G4Gamma.hh"
72
73 #include "G4Log.hh"
74 #include <iostream>
75
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78 G4EmDataHandler* G4GammaGeneralProcess::theHandler = nullptr;
79 G4bool G4GammaGeneralProcess::theT[nTables] =
80 {true,false,true,true,true,false,true,true,true,
81 true,true,true,true,true,true};
82 G4String G4GammaGeneralProcess::nameT[nTables] =
83 {"0","1","2","3","4","5","6","7","8",
84 "9","10","11","12","13","14"};
85
86 G4GammaGeneralProcess::G4GammaGeneralProcess(const G4String& pname):
87 G4VEmProcess(pname, fElectromagnetic),
88 minPEEnergy(150*CLHEP::keV),
89 minEEEnergy(2*CLHEP::electron_mass_c2),
90 minMMEnergy(100*CLHEP::MeV)
91 {
92 SetVerboseLevel(1);
93 SetParticle(G4Gamma::Gamma());
94 SetProcessSubType(fGammaGeneralProcess);
95 if (nullptr == theHandler) {
96 theHandler = new G4EmDataHandler(nTables);
97 }
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102 G4GammaGeneralProcess::~G4GammaGeneralProcess()
103 {
104 if(isTheMaster) {
105 delete theHandler;
106 theHandler = nullptr;
107 }
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
112 G4bool G4GammaGeneralProcess::IsApplicable(const G4ParticleDefinition&)
113 {
114 return true;
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119 void G4GammaGeneralProcess::AddEmProcess(G4VEmProcess* ptr)
120 {
121 if(nullptr == ptr) { return; }
122 G4int stype = ptr->GetProcessSubType();
123 if(stype == fRayleigh) { theRayleigh = ptr; }
124 else if(stype == fPhotoElectricEffect) { thePhotoElectric = ptr; }
125 else if(stype == fComptonScattering) { theCompton = ptr; }
126 else if(stype == fGammaConversion) { theConversionEE = ptr; }
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130
131 void G4GammaGeneralProcess::AddMMProcess(G4GammaConversionToMuons* ptr)
132 {
133 theConversionMM = ptr;
134 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137
138 void G4GammaGeneralProcess::AddHadProcess(G4HadronicProcess* ptr)
139 {
140 theGammaNuclear = ptr;
141 }
142
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144
145 void G4GammaGeneralProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
146 {
147 SetParticle(&part);
148 preStepLambda = 0.0;
149 idxEnergy = 0;
150 currentCouple = nullptr;
151
152 G4LossTableManager* man = G4LossTableManager::Instance();
153
154 // initialise base material for the current run
155 G4LossTableBuilder* bld = man->GetTableBuilder();
156 bld->InitialiseBaseMaterials();
157 baseMat = bld->GetBaseMaterialFlag();
158
159 if(1 < verboseLevel) {
160 G4cout << "G4GammaGeneralProcess::PreparePhysicsTable() for "
161 << GetProcessName()
162 << " and particle " << part.GetParticleName()
163 << " isMaster: " << isTheMaster << G4endl;
164 }
165
166 // 3 sub-processes must be always defined
167 if(thePhotoElectric == nullptr || theCompton == nullptr ||
168 theConversionEE == nullptr) {
169 G4ExceptionDescription ed;
170 ed << "### G4GeneralGammaProcess is initialized incorrectly"
171 << "\n Photoelectric: " << thePhotoElectric
172 << "\n Compton: " << theCompton
173 << "\n Conversion: " << theConversionEE;
174 G4Exception("G4GeneralGammaProcess","em0004",
175 FatalException, ed,"");
176 return;
177 }
178
179 thePhotoElectric->PreparePhysicsTable(part);
180 theCompton->PreparePhysicsTable(part);
181 theConversionEE->PreparePhysicsTable(part);
182 if (nullptr != theRayleigh) { theRayleigh->PreparePhysicsTable(part); }
183 if (nullptr != theGammaNuclear) { theGammaNuclear->PreparePhysicsTable(part); }
184 if (nullptr != theConversionMM) { theConversionMM->PreparePhysicsTable(part); }
185
186 InitialiseProcess(&part);
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
190
191 void G4GammaGeneralProcess::InitialiseProcess(const G4ParticleDefinition*)
192 {
193 if(isTheMaster) {
194
195 G4EmParameters* param = G4EmParameters::Instance();
196 G4LossTableManager* man = G4LossTableManager::Instance();
197
198 // tables are created and its size is defined only once
199 if (nullptr != theRayleigh) { theT[1] = true; }
200
201 theHandler->SetMasterProcess(thePhotoElectric);
202 theHandler->SetMasterProcess(theCompton);
203 theHandler->SetMasterProcess(theConversionEE);
204 theHandler->SetMasterProcess(theRayleigh);
205
206 auto bld = man->GetTableBuilder();
207
208 const G4ProductionCutsTable* theCoupleTable=
209 G4ProductionCutsTable::GetProductionCutsTable();
210 std::size_t numOfCouples = theCoupleTable->GetTableSize();
211
212 G4double mine = param->MinKinEnergy();
213 G4double maxe = param->MaxKinEnergy();
214 G4int nd = param->NumberOfBinsPerDecade();
215 std::size_t nbin1 = std::max(5, nd*G4lrint(std::log10(minPEEnergy/mine)));
216 std::size_t nbin2 = std::max(5, nd*G4lrint(std::log10(maxe/minMMEnergy)));
217
218 G4PhysicsVector* vec = nullptr;
219 G4PhysicsLogVector aVector(mine,minPEEnergy,nbin1,true);
220 G4PhysicsLogVector bVector(minPEEnergy,minEEEnergy,nLowE,false);
221 G4PhysicsLogVector cVector(minEEEnergy,minMMEnergy,nHighE,false);
222 G4PhysicsLogVector dVector(minMMEnergy,maxe,nbin2,true);
223
224 for(std::size_t i=0; i<nTables; ++i) {
225 if(!theT[i]) { continue; }
226 //G4cout << "## PreparePhysTable " << i << "." << G4endl;
227 G4PhysicsTable* table = theHandler->MakeTable(i);
228 //G4cout << " make table " << table << G4endl;
229 for(std::size_t j=0; j<numOfCouples; ++j) {
230 vec = (*table)[j];
231 if (bld->GetFlag(j) && nullptr == vec) {
232 if(i<=1) {
233 vec = new G4PhysicsVector(aVector);
234 } else if(i<=5) {
235 vec = new G4PhysicsVector(bVector);
236 } else if(i<=9) {
237 vec = new G4PhysicsVector(cVector);
238 } else {
239 vec = new G4PhysicsVector(dVector);
240 }
241 G4PhysicsTableHelper::SetPhysicsVector(table, j, vec);
242 }
243 }
244 }
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249
250 void G4GammaGeneralProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
251 {
252 if(1 < verboseLevel) {
253 G4cout << "### G4VEmProcess::BuildPhysicsTable() for "
254 << GetProcessName()
255 << " and particle " << part.GetParticleName()
256 << G4endl;
257 }
258 if(!isTheMaster) {
259 thePhotoElectric->SetEmMasterProcess(theHandler->GetMasterProcess(0));
260 baseMat = theHandler->GetMasterProcess(0)->UseBaseMaterial();
261 }
262 thePhotoElectric->BuildPhysicsTable(part);
263
264 if(!isTheMaster) {
265 theCompton->SetEmMasterProcess(theHandler->GetMasterProcess(1));
266 }
267 theCompton->BuildPhysicsTable(part);
268
269 if(!isTheMaster) {
270 theConversionEE->SetEmMasterProcess(theHandler->GetMasterProcess(2));
271 }
272 theConversionEE->BuildPhysicsTable(part);
273
274 if(theRayleigh != nullptr) {
275 if(!isTheMaster) {
276 theRayleigh->SetEmMasterProcess(theHandler->GetMasterProcess(3));
277 }
278 theRayleigh->BuildPhysicsTable(part);
279 }
280 if(theGammaNuclear != nullptr) { theGammaNuclear->BuildPhysicsTable(part); }
281 if(theConversionMM != nullptr) { theConversionMM->BuildPhysicsTable(part); }
282
283 if(isTheMaster) {
284 const G4ProductionCutsTable* theCoupleTable=
285 G4ProductionCutsTable::GetProductionCutsTable();
286 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
287
288 G4LossTableBuilder* bld = G4LossTableManager::Instance()->GetTableBuilder();
289 const std::vector<G4PhysicsTable*>& tables = theHandler->GetTables();
290
291 G4CrossSectionDataStore* gn = (nullptr != theGammaNuclear)
292 ? theGammaNuclear->GetCrossSectionDataStore() : nullptr;
293 G4DynamicParticle* dynParticle =
294 new G4DynamicParticle(G4Gamma::Gamma(),G4ThreeVector(1,0,0),1.0);
295
296 G4double sigComp(0.), sigPE(0.), sigConv(0.), sigR(0.),
297 sigN(0.), sigM(0.), val(0.);
298
299 for(G4int i=0; i<numOfCouples; ++i) {
300
301 if (bld->GetFlag(i)) {
302 G4int idx = (!baseMat) ? i : DensityIndex(i);
303 const G4MaterialCutsCouple* couple =
304 theCoupleTable->GetMaterialCutsCouple(i);
305 const G4Material* material = couple->GetMaterial();
306
307 // energy interval 0
308 std::size_t nn = (*(tables[0]))[idx]->GetVectorLength();
309 if(1 < verboseLevel) {
310 G4cout << "======= Zone 0 ======= N= " << nn
311 << " for " << material->GetName() << G4endl;
312 }
313 for(std::size_t j=0; j<nn; ++j) {
314 G4double e = (*(tables[0]))[idx]->Energy(j);
315 G4double loge = G4Log(e);
316 sigComp = theCompton->GetLambda(e, couple, loge);
317 sigR = (nullptr != theRayleigh) ?
318 theRayleigh->GetLambda(e, couple, loge) : 0.0;
319 G4double sum = sigComp + sigR;
320 if(1 < verboseLevel) {
321 G4cout << j << ". E= " << e << " xs= " << sum
322 << " compt= " << sigComp << " Rayl= " << sigR << G4endl;
323 }
324 (*(tables[0]))[idx]->PutValue(j, sum);
325 if(theT[1]) {
326 val = sigR/sum;
327 (*(tables[1]))[idx]->PutValue(j, val);
328 }
329 }
330
331 // energy interval 1
332 nn = (*(tables[2]))[idx]->GetVectorLength();
333 if(1 < verboseLevel) {
334 G4cout << "======= Zone 1 ======= N= " << nn << G4endl;
335 }
336 for(std::size_t j=0; j<nn; ++j) {
337 G4double e = (*(tables[2]))[idx]->Energy(j);
338 G4double loge = G4Log(e);
339 sigComp = theCompton->GetLambda(e, couple, loge);
340 sigR = (nullptr != theRayleigh) ?
341 theRayleigh->GetLambda(e, couple, loge) : 0.0;
342 sigPE = thePhotoElectric->GetLambda(e, couple, loge);
343 G4double sum = sigComp + sigR + sigPE;
344 if(1 < verboseLevel) {
345 G4cout << j << ". E= " << e << " xs= " << sum
346 << " compt= " << sigComp << " conv= " << sigConv
347 << " PE= " << sigPE << " Rayl= " << sigR
348 << " GN= " << sigN << G4endl;
349 }
350 (*(tables[2]))[idx]->PutValue(j, sum);
351
352 val = sigPE/sum;
353 (*(tables[3]))[idx]->PutValue(j, val);
354
355 val = (sigR > 0.0) ? (sigComp + sigPE)/sum : 1.0;
356 (*(tables[4]))[idx]->PutValue(j, val);
357 }
358
359 // energy interval 2
360 nn = (*(tables[6]))[idx]->GetVectorLength();
361 if(1 < verboseLevel) {
362 G4cout << "======= Zone 2 ======= N= " << nn << G4endl;
363 }
364 for(std::size_t j=0; j<nn; ++j) {
365 G4double e = (*(tables[6]))[idx]->Energy(j);
366 G4double loge = G4Log(e);
367 sigComp = theCompton->GetLambda(e, couple, loge);
368 sigConv = theConversionEE->GetLambda(e, couple, loge);
369 sigPE = thePhotoElectric->GetLambda(e, couple, loge);
370 sigN = 0.0;
371 if(nullptr != gn) {
372 dynParticle->SetKineticEnergy(e);
373 sigN = gn->ComputeCrossSection(dynParticle, material);
374 }
375 G4double sum = sigComp + sigConv + sigPE + sigN;
376 if(1 < verboseLevel) {
377 G4cout << j << ". E= " << e << " xs= " << sum
378 << " compt= " << sigComp << " conv= " << sigConv
379 << " PE= " << sigPE
380 << " GN= " << sigN << G4endl;
381 }
382 (*(tables[6]))[idx]->PutValue(j, sum);
383
384 val = sigConv/sum;
385 (*(tables[7]))[idx]->PutValue(j, val);
386
387 val = (sigConv + sigComp)/sum;
388 (*(tables[8]))[idx]->PutValue(j, val);
389
390 val = (sigN > 0.0) ? (sigConv + sigComp + sigPE)/sum : 1.0;
391 (*(tables[9]))[idx]->PutValue(j, val);
392 }
393
394 // energy interval 3
395 nn = (*(tables[10]))[idx]->GetVectorLength();
396 if(1 < verboseLevel) {
397 G4cout << "======= Zone 3 ======= N= " << nn
398 << " for " << material->GetName() << G4endl;
399 }
400 for(std::size_t j=0; j<nn; ++j) {
401 G4double e = (*(tables[10]))[idx]->Energy(j);
402 G4double loge = G4Log(e);
403 sigComp = theCompton->GetLambda(e, couple, loge);
404 sigConv = theConversionEE->GetLambda(e, couple, loge);
405 sigPE = thePhotoElectric->GetLambda(e, couple, loge);
406 sigN = 0.0;
407 if(nullptr != gn) {
408 dynParticle->SetKineticEnergy(e);
409 sigN = gn->ComputeCrossSection(dynParticle, material);
410 }
411 sigM = 0.0;
412 if(nullptr != theConversionMM) {
413 val = theConversionMM->ComputeMeanFreePath(e, material);
414 sigM = (val < DBL_MAX) ? 1./val : 0.0;
415 }
416 G4double sum = sigComp + sigConv + sigPE + sigN + sigM;
417 if(1 < verboseLevel) {
418 G4cout << j << ". E= " << e << " xs= " << sum
419 << " compt= " << sigComp << " conv= " << sigConv
420 << " PE= " << sigPE
421 << " GN= " << sigN << G4endl;
422 }
423 (*(tables[10]))[idx]->PutValue(j, sum);
424
425 val = (sigComp + sigPE + sigN + sigM)/sum;
426 (*(tables[11]))[idx]->PutValue(j, val);
427
428 val = (sigPE + sigN + sigM)/sum;
429 (*(tables[12]))[idx]->PutValue(j, val);
430
431 val = (sigN + sigM)/sum;
432 (*(tables[13]))[idx]->PutValue(j, val);
433
434 val = sigM/sum;
435 (*(tables[14]))[idx]->PutValue(j, val);
436 }
437 for(std::size_t k=0; k<nTables; ++k) {
438 if(theT[k] && (k <= 1 || k >= 10)) {
439 //G4cout <<"BuildPhysTable spline iTable="<<k<<" jCouple="<< idx << G4endl;
440 (*(tables[k]))[idx]->FillSecondDerivatives();
441 }
442 }
443 }
444 }
445 delete dynParticle;
446 }
447
448 if(1 < verboseLevel) {
449 G4cout << "### G4VEmProcess::BuildPhysicsTable() done for "
450 << GetProcessName()
451 << " and particle " << part.GetParticleName()
452 << G4endl;
453 }
454 }
455
456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
457
458 void G4GammaGeneralProcess::StartTracking(G4Track*)
459 {
460 theNumberOfInteractionLengthLeft = -1.0;
461 }
462
463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
464
465 G4double G4GammaGeneralProcess::PostStepGetPhysicalInteractionLength(
466 const G4Track& track,
467 G4double previousStepSize,
468 G4ForceCondition* condition)
469 {
470 *condition = NotForced;
471 G4double x = DBL_MAX;
472
473 G4double energy = track.GetKineticEnergy();
474 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
475
476 // compute mean free path
477 G4bool recompute = false;
478 if(couple != currentCouple) {
479 currentCouple = couple;
480 basedCoupleIndex = currentCoupleIndex = couple->GetIndex();
481 currentMaterial = couple->GetMaterial();
482 factor = 1.0;
483 if(baseMat) {
484 basedCoupleIndex = DensityIndex((G4int)currentCoupleIndex);
485 factor = DensityFactor((G4int)currentCoupleIndex);
486 }
487 recompute = true;
488 }
489 if(energy != preStepKinEnergy) {
490 preStepKinEnergy = energy;
491 preStepLogE = track.GetDynamicParticle()->GetLogKineticEnergy();
492 recompute = true;
493 }
494 if(recompute) {
495 preStepLambda = TotalCrossSectionPerVolume();
496
497 // zero cross section
498 if(preStepLambda <= 0.0) {
499 theNumberOfInteractionLengthLeft = -1.0;
500 currentInteractionLength = DBL_MAX;
501 }
502 }
503
504 // non-zero cross section
505 if(preStepLambda > 0.0) {
506
507 if (theNumberOfInteractionLengthLeft < 0.0) {
508
509 // beggining of tracking (or just after DoIt of this process)
510 theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
511 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
512
513 } else if(currentInteractionLength < DBL_MAX) {
514
515 theNumberOfInteractionLengthLeft -=
516 previousStepSize/currentInteractionLength;
517 theNumberOfInteractionLengthLeft =
518 std::max(theNumberOfInteractionLengthLeft, 0.0);
519 }
520
521 // new mean free path and step limit for the next step
522 currentInteractionLength = 1.0/preStepLambda;
523 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
524 }
525 /*
526 G4cout << "PostStepGetPhysicalInteractionLength: e= " << energy
527 << " idxe= " << idxEnergy << " xs= " << preStepLambda
528 << " x= " << x << G4endl;
529 */
530 return x;
531 }
532
533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
534
535 G4double G4GammaGeneralProcess::TotalCrossSectionPerVolume()
536 {
537 G4double cross = 0.0;
538 /*
539 G4cout << "#Total: " << preStepKinEnergy << " " << minPEEnergy << " "
540 << minEEEnergy << " " << minMMEnergy<< G4endl;
541 G4cout << " idxE= " << idxEnergy
542 << " idxC= " << currentCoupleIndex << G4endl;
543 */
544 if(preStepKinEnergy < minPEEnergy) {
545 cross = ComputeGeneralLambda(0, 0);
546 //G4cout << "XS1: " << cross << G4endl;
547 peLambda = thePhotoElectric->GetLambda(preStepKinEnergy, currentCouple, preStepLogE);
548 cross += peLambda;
549 //G4cout << "XS2: " << peLambda << G4endl;
550
551 } else if(preStepKinEnergy < minEEEnergy) {
552 cross = ComputeGeneralLambda(1, 2);
553 //G4cout << "XS3: " << cross << G4endl;
554
555 } else if(preStepKinEnergy < minMMEnergy) {
556 cross = ComputeGeneralLambda(2, 6);
557 //G4cout << "XS4: " << cross << G4endl;
558
559 } else {
560 cross = ComputeGeneralLambda(3, 10);
561 //G4cout << "XS5: " << cross << G4endl;
562 }
563 /*
564 G4cout << "xs= " << cross << " idxE= " << idxEnergy
565 << " idxC= " << currentCoupleIndex
566 << " E= " << preStepKinEnergy << G4endl;
567 */
568 return cross;
569 }
570
571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
572
573 G4VParticleChange* G4GammaGeneralProcess::PostStepDoIt(const G4Track& track,
574 const G4Step& step)
575 {
576 // In all cases clear number of interaction lengths
577 theNumberOfInteractionLengthLeft = -1.0;
578 selectedProc = nullptr;
579 G4double q = G4UniformRand();
580 /*
581 G4cout << "PostStep: preStepLambda= " << preStepLambda
582 << " PE= " << peLambda << " q= " << q << " idxE= " << idxEnergy
583 << G4endl;
584 */
585 switch (idxEnergy) {
586 case 0:
587 if(preStepLambda*q <= peLambda) {
588 SelectEmProcess(step, thePhotoElectric);
589 } else {
590 if(theT[1] && preStepLambda*q < (preStepLambda - peLambda)*GetProbability(1) + peLambda) {
591 SelectEmProcess(step, theRayleigh);
592 } else {
593 SelectEmProcess(step, theCompton);
594 }
595 }
596 break;
597
598 case 1:
599 if(q <= GetProbability(3)) {
600 SelectEmProcess(step, thePhotoElectric);
601 } else if(q <= GetProbability(4)) {
602 SelectEmProcess(step, theCompton);
603 } else if(theRayleigh) {
604 SelectEmProcess(step, theRayleigh);
605 } else {
606 SelectEmProcess(step, thePhotoElectric);
607 }
608 break;
609
610 case 2:
611 if(q <= GetProbability(7)) {
612 SelectEmProcess(step, theConversionEE);
613 } else if(q <= GetProbability(8)) {
614 SelectEmProcess(step, theCompton);
615 } else if(q <= GetProbability(9)) {
616 SelectEmProcess(step, thePhotoElectric);
617 } else if(theGammaNuclear) {
618 SelectHadProcess(track, step, theGammaNuclear);
619 } else {
620 SelectEmProcess(step, theConversionEE);
621 }
622 break;
623
624 case 3:
625 if(q + GetProbability(11) <= 1.0) {
626 SelectEmProcess(step, theConversionEE);
627 } else if(q + GetProbability(12) <= 1.0) {
628 SelectEmProcess(step, theCompton);
629 } else if(q + GetProbability(13) <= 1.0) {
630 SelectEmProcess(step, thePhotoElectric);
631 } else if(theGammaNuclear && q + GetProbability(14) <= 1.0) {
632 SelectHadProcess(track, step, theGammaNuclear);
633 } else if(theConversionMM) {
634 SelectedProcess(step, theConversionMM);
635 } else {
636 SelectEmProcess(step, theConversionEE);
637 }
638 break;
639 }
640 // sample secondaries
641 if(selectedProc != nullptr) {
642 return selectedProc->PostStepDoIt(track, step);
643 }
644 // no interaction - exception case
645 fParticleChange.InitializeForPostStep(track);
646 return &fParticleChange;
647 }
648
649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
650
651 void G4GammaGeneralProcess::SelectHadProcess(const G4Track& track,
652 const G4Step& step, G4HadronicProcess* proc)
653 {
654 SelectedProcess(step, proc);
655 proc->GetCrossSectionDataStore()->ComputeCrossSection(track.GetDynamicParticle(),
656 currentMaterial);
657 }
658
659 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
660
661 G4bool G4GammaGeneralProcess::StorePhysicsTable(const G4ParticleDefinition* part,
662 const G4String& directory,
663 G4bool ascii)
664 {
665 G4bool yes = true;
666 if(!isTheMaster) { return yes; }
667 if(!thePhotoElectric->StorePhysicsTable(part, directory, ascii))
668 { yes = false; }
669 if(!theCompton->StorePhysicsTable(part, directory, ascii))
670 { yes = false; }
671 if(!theConversionEE->StorePhysicsTable(part, directory, ascii))
672 { yes = false; }
673 if(theRayleigh != nullptr &&
674 !theRayleigh->StorePhysicsTable(part, directory, ascii))
675 { yes = false; }
676
677 for(std::size_t i=0; i<nTables; ++i) {
678 if(theT[i]) {
679 G4String nam = (0==i || 2==i || 6==i || 10==i)
680 ? "LambdaGeneral" + nameT[i] : "ProbGeneral" + nameT[i];
681 G4String fnam = GetPhysicsTableFileName(part,directory,nam,ascii);
682 if(!theHandler->StorePhysicsTable(i, part, fnam, ascii)) { yes = false; }
683 }
684 }
685 return yes;
686 }
687
688 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
689
690 G4bool
691 G4GammaGeneralProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
692 const G4String& directory,
693 G4bool ascii)
694 {
695 if (!isTheMaster) { return true; }
696 if(1 < verboseLevel) {
697 G4cout << "G4GammaGeneralProcess::RetrievePhysicsTable() for "
698 << part->GetParticleName() << " and process "
699 << GetProcessName() << G4endl;
700 }
701 G4bool yes = true;
702 if(!thePhotoElectric->RetrievePhysicsTable(part, directory, ascii))
703 { yes = false; }
704 if(!theCompton->RetrievePhysicsTable(part, directory, ascii))
705 { yes = false; }
706 if(!theConversionEE->RetrievePhysicsTable(part, directory, ascii))
707 { yes = false; }
708 if(theRayleigh != nullptr &&
709 !theRayleigh->RetrievePhysicsTable(part, directory, ascii))
710 { yes = false; }
711
712 for(std::size_t i=0; i<nTables; ++i) {
713 if(theT[i]) {
714 G4String nam = (0==i || 2==i || 6==i || 10==i)
715 ? "LambdaGeneral" + nameT[i] : "ProbGeneral" + nameT[i];
716 G4String fnam = GetPhysicsTableFileName(part,directory,nam,ascii);
717 G4bool spline = (i <= 1 || i >= 10);
718 if(!theHandler->RetrievePhysicsTable(i, part, fnam, ascii, spline))
719 { yes = false; }
720 }
721 }
722 return yes;
723 }
724
725 //....Ooooo0ooooo ........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
726
727 G4double G4GammaGeneralProcess::GetMeanFreePath(const G4Track& track,
728 G4double,
729 G4ForceCondition* condition)
730 {
731 *condition = NotForced;
732 return MeanFreePath(track);
733 }
734
735 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
736
737 void G4GammaGeneralProcess::ProcessDescription(std::ostream& out) const
738 {
739 thePhotoElectric->ProcessDescription(out);
740 theCompton->ProcessDescription(out);
741 theConversionEE->ProcessDescription(out);
742 if(theRayleigh) { theRayleigh->ProcessDescription(out); }
743 if(theGammaNuclear) { theGammaNuclear->ProcessDescription(out); }
744 if(theConversionMM) { theConversionMM->ProcessDescription(out); }
745 }
746
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
748
749 const G4String& G4GammaGeneralProcess::GetSubProcessName() const
750 {
751 return (selectedProc) ? selectedProc->GetProcessName()
752 : G4VProcess::GetProcessName();
753 }
754
755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
756
757 G4int G4GammaGeneralProcess::GetSubProcessSubType() const
758 {
759 return (selectedProc) ? selectedProc->GetProcessSubType()
760 : fGammaGeneralProcess;
761 }
762
763 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
764
765 G4VEmProcess* G4GammaGeneralProcess::GetEmProcess(const G4String& name)
766 {
767 G4VEmProcess* proc = nullptr;
768 if(name == thePhotoElectric->GetProcessName()) {
769 proc = thePhotoElectric;
770 } else if(name == theCompton->GetProcessName()) {
771 proc = theCompton;
772 } else if(name == theConversionEE->GetProcessName()) {
773 proc = theConversionEE;
774 } else if(theRayleigh != nullptr && name == theRayleigh->GetProcessName()) {
775 proc = theRayleigh;
776 }
777 return proc;
778 }
779
780 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
781
782 const G4VProcess* G4GammaGeneralProcess::GetCreatorProcess() const
783 {
784 return selectedProc;
785 }
786
787 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788