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