Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreGammaConversionModel.cc
Go to the documentation of this file.
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 // Author: Sebastien Incerti
27 // 22 January 2012
28 // on base of G4LivermoreGammaConversionModel (original version)
29 // and G4LivermoreRayleighModel (MT version)
30 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Log.hh"
35 #include "G4Exp.hh"
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 
39 using namespace std;
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 G4int G4LivermoreGammaConversionModel::maxZ = 99;
44 G4LPhysicsFreeVector* G4LivermoreGammaConversionModel::data[] = {nullptr};
45 
47 (const G4ParticleDefinition*, const G4String& nam)
48 :G4VEmModel(nam),isInitialised(false),smallEnergy(2.*MeV)
49 {
50  fParticleChange = nullptr;
51 
52  lowEnergyLimit = 2.0*electron_mass_c2;
53 
54  verboseLevel= 0;
55  // Verbosity scale for debugging purposes:
56  // 0 = nothing
57  // 1 = calculation of cross sections, file openings...
58  // 2 = entering in methods
59 
60  if(verboseLevel > 0)
61  {
62  G4cout << "G4LivermoreGammaConversionModel is constructed " << G4endl;
63  }
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
69 {
70  if(IsMaster()) {
71  for(G4int i=0; i<maxZ; ++i) {
72  if(data[i]) {
73  delete data[i];
74  data[i] = 0;
75  }
76  }
77  }
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
83  const G4ParticleDefinition* particle,
84  const G4DataVector& cuts)
85 {
86  if (verboseLevel > 1)
87  {
88  G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel."
89  << G4endl
90  << "Energy range: "
91  << LowEnergyLimit() / MeV << " MeV - "
92  << HighEnergyLimit() / GeV << " GeV"
93  << G4endl;
94  }
95 
96  if(IsMaster())
97  {
98  // Initialise element selector
99  InitialiseElementSelectors(particle, cuts);
100 
101  // Access to elements
102  char* path = getenv("G4LEDATA");
103 
104  G4ProductionCutsTable* theCoupleTable =
106 
107  G4int numOfCouples = theCoupleTable->GetTableSize();
108 
109  for(G4int i=0; i<numOfCouples; ++i)
110  {
111  const G4Material* material =
112  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
113  const G4ElementVector* theElementVector = material->GetElementVector();
114  G4int nelm = material->GetNumberOfElements();
115 
116  for (G4int j=0; j<nelm; ++j)
117  {
118  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
119  if(Z < 1) { Z = 1; }
120  else if(Z > maxZ) { Z = maxZ; }
121  if(!data[Z]) { ReadData(Z, path); }
122  }
123  }
124  }
125  if(isInitialised) { return; }
126  fParticleChange = GetParticleChangeForGamma();
127  isInitialised = true;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131 
133  const G4ParticleDefinition*, G4VEmModel* masterModel)
134 {
135  SetElementSelectors(masterModel->GetElementSelectors());
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 
140 G4double
142  const G4ParticleDefinition*,
143  G4double)
144 {
145  return lowEnergyLimit;
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 
150 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path)
151 {
152  if (verboseLevel > 1)
153  {
154  G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
155  << G4endl;
156  }
157 
158  if(data[Z]) { return; }
159 
160  const char* datadir = path;
161 
162  if(!datadir)
163  {
164  datadir = getenv("G4LEDATA");
165  if(!datadir)
166  {
167  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
168  "em0006",FatalException,
169  "Environment variable G4LEDATA not defined");
170  return;
171  }
172  }
173 
174  //
175 
176  data[Z] = new G4LPhysicsFreeVector();
177 
178  //
179 
180  std::ostringstream ost;
181  ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
182  std::ifstream fin(ost.str().c_str());
183 
184  if( !fin.is_open())
185  {
187  ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
188  << "> is not opened!" << G4endl;
189  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
190  "em0003",FatalException,
191  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
192  return;
193  }
194 
195  else
196  {
197 
198  if(verboseLevel > 3) { G4cout << "File " << ost.str()
199  << " is opened by G4LivermoreGammaConversionModel" << G4endl;}
200 
201  data[Z]->Retrieve(fin, true);
202  }
203 
204  // Activation of spline interpolation
205  data[Z] ->SetSpline(true);
206 
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
211 G4double
213  G4double GammaEnergy,
214  G4double Z, G4double,
216 {
217  if (verboseLevel > 1)
218  {
219  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
220  << G4endl;
221  }
222 
223  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
224 
225  G4double xs = 0.0;
226 
227  G4int intZ=G4int(Z);
228 
229  if(intZ < 1 || intZ > maxZ) { return xs; }
230 
231  G4LPhysicsFreeVector* pv = data[intZ];
232 
233  // if element was not initialised
234  // do initialisation safely for MT mode
235  if(!pv)
236  {
237  InitialiseForElement(0, intZ);
238  pv = data[intZ];
239  if(!pv) { return xs; }
240  }
241  // x-section is taken from the table
242  xs = pv->Value(GammaEnergy);
243 
244  if(verboseLevel > 0)
245  {
246  G4int n = pv->GetVectorLength() - 1;
247  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
248  << GammaEnergy/MeV << G4endl;
249  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
250  G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
251  G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
252  G4cout << "*********************************************************" << G4endl;
253  }
254 
255  return xs;
256 
257 }
258 
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260 
262  std::vector<G4DynamicParticle*>* fvect,
263  const G4MaterialCutsCouple* couple,
264  const G4DynamicParticle* aDynamicGamma,
266 {
267 
268 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
269 // cross sections with Coulomb correction. A modified version of the random
270 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
271 
272 // Note 1 : Effects due to the breakdown of the Born approximation at low
273 // energy are ignored.
274 // Note 2 : The differential cross section implicitly takes account of
275 // pair creation in both nuclear and atomic electron fields. However triplet
276 // prodution is not generated.
277 
278  if (verboseLevel > 1) {
279  G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel"
280  << G4endl;
281  }
282 
283  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
284  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
285 
286  G4double epsilon ;
287  G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
288 
289  // Do it fast if photon energy < 2. MeV
290  if (photonEnergy < smallEnergy )
291  {
292  epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
293  }
294  else
295  {
296  // Select randomly one element in the current material
297 
298  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
299  const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
300 
301  if (element == 0)
302  {
303  G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
304  << G4endl;
305  return;
306  }
307  G4IonisParamElm* ionisation = element->GetIonisation();
308  if (ionisation == 0)
309  {
310  G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
311  << G4endl;
312  return;
313  }
314 
315  // Extract Coulomb factor for this Element
316  G4double fZ = 8. * (ionisation->GetlogZ3());
317  if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
318 
319  // Limits of the screening variable
320  G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
321  G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ;
322  G4double screenMin = std::min(4.*screenFactor,screenMax) ;
323 
324  // Limits of the energy sampling
325  G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
326  G4double epsilonMin = std::max(epsilon0Local,epsilon1);
327  G4double epsilonRange = 0.5 - epsilonMin ;
328 
329  // Sample the energy rate of the created electron (or positron)
330  G4double screen;
331  G4double gReject ;
332 
333  G4double f10 = ScreenFunction1(screenMin) - fZ;
334  G4double f20 = ScreenFunction2(screenMin) - fZ;
335  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
336  G4double normF2 = std::max(1.5 * f20,0.);
337 
338  do
339  {
340  if (normF1 / (normF1 + normF2) > G4UniformRand() )
341  {
342  epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
343  screen = screenFactor / (epsilon * (1. - epsilon));
344  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
345  }
346  else
347  {
348  epsilon = epsilonMin + epsilonRange * G4UniformRand();
349  screen = screenFactor / (epsilon * (1 - epsilon));
350  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
351  }
352  } while ( gReject < G4UniformRand() );
353 
354  } // End of epsilon sampling
355 
356  // Fix charges randomly
357 
358  G4double electronTotEnergy;
359  G4double positronTotEnergy;
360 
361  if (G4UniformRand() > 0.5)
362  {
363  electronTotEnergy = (1. - epsilon) * photonEnergy;
364  positronTotEnergy = epsilon * photonEnergy;
365  }
366  else
367  {
368  positronTotEnergy = (1. - epsilon) * photonEnergy;
369  electronTotEnergy = epsilon * photonEnergy;
370  }
371 
372  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
373  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
374  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
375 
376  G4double u;
377  const G4double a1 = 0.625;
378  G4double a2 = 3. * a1;
379  // G4double d = 27. ;
380 
381  // if (9. / (9. + d) > G4UniformRand())
382  if (0.25 > G4UniformRand())
383  {
384  u = - G4Log(G4UniformRand() * G4UniformRand()) / a1 ;
385  }
386  else
387  {
388  u = - G4Log(G4UniformRand() * G4UniformRand()) / a2 ;
389  }
390 
391  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
392  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
393  G4double phi = twopi * G4UniformRand();
394 
395  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
396  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
397 
398 
399  // Kinematics of the created pair:
400  // the electron and positron are assumed to have a symetric angular
401  // distribution with respect to the Z axis along the parent photon
402 
403  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
404 
405  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
406  electronDirection.rotateUz(photonDirection);
407 
409  electronDirection,
410  electronKineEnergy);
411 
412  // The e+ is always created
413  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
414 
415  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
416  positronDirection.rotateUz(photonDirection);
417 
418  // Create G4DynamicParticle object for the particle2
420  positronDirection,
421  positronKineEnergy);
422  // Fill output vector
423  fvect->push_back(particle1);
424  fvect->push_back(particle2);
425 
426  // kill incident photon
427  fParticleChange->SetProposedKineticEnergy(0.);
428  fParticleChange->ProposeTrackStatus(fStopAndKill);
429 
430 }
431 
432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
433 
434 G4double
435 G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable)
436 {
437  // Compute the value of the screening function 3*phi1 - phi2
438 
439  G4double value;
440 
441  if (screenVariable > 1.)
442  value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
443  else
444  value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
445 
446  return value;
447 }
448 
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450 
451 G4double
452 G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable)
453 {
454  // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
455 
456  G4double value;
457 
458  if (screenVariable > 1.)
459  value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
460  else
461  value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
462 
463  return value;
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
467 
468 #include "G4AutoLock.hh"
469 namespace { G4Mutex LivermoreGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
470 
472  const G4ParticleDefinition*,
473  G4int Z)
474 {
475  G4AutoLock l(&LivermoreGammaConversionModelMutex);
476  // G4cout << "G4LivermoreGammaConversionModel::InitialiseForElement Z= "
477  // << Z << G4endl;
478  if(!data[Z]) { ReadData(Z); }
479  l.unlock();
480 }
481 
482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4double GetfCoulomb() const
Definition: G4Element.hh:191
G4ParticleDefinition * GetDefinition() const
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
static constexpr double twopi
Definition: G4SIunits.hh:76
const XML_Char const XML_Char * data
Definition: expat.h:268
string material
Definition: eplot.py:19
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4LivermoreGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreConversion")
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
float electron_mass_c2
Definition: hepunit.py:274
const G4int n
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double GetlogZ3() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4ProductionCutsTable * GetProductionCutsTable()
G4int G4Mutex
Definition: G4Threading.hh:173
static G4Positron * Positron()
Definition: G4Positron.cc:94
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:199
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
double G4double
Definition: G4Types.hh:76
G4double GetZ3() const
double epsilon(double density, double temperature)
const G4Material * GetMaterial() const