ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreGammaConversionModelRC.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermoreGammaConversionModelRC.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 // Author: Sebastien Incerti
27 // 22 January 2012
28 // on base of G4LivermoreGammaConversionModelRC (original version)
29 // and G4LivermoreRayleighModel (MT version)
30 
31 #include "G4LivermoreGammaConversionModelRC.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Log.hh"
35 #include "G4Electron.hh"
36 #include "G4Positron.hh"
37 #include "G4Gamma.hh"
38 #include "G4ParticleChangeForGamma.hh"
39 #include "G4Exp.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 using namespace std;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
47 G4int G4LivermoreGammaConversionModelRC::maxZ = 99;
48 G4LPhysicsFreeVector* G4LivermoreGammaConversionModelRC::data[] = {nullptr};
49 
50 G4LivermoreGammaConversionModelRC::G4LivermoreGammaConversionModelRC
51 (const G4ParticleDefinition*, const G4String& nam)
52 :G4VEmModel(nam),isInitialised(false),smallEnergy(2.*MeV)
53 {
54  fParticleChange = nullptr;
55 
56  lowEnergyLimit = 2.0*electron_mass_c2;
57 
58  verboseLevel= 0;
59  // Verbosity scale for debugging purposes:
60  // 0 = nothing
61  // 1 = calculation of cross sections, file openings...
62  // 2 = entering in methods
63 
64  if(verboseLevel > 0)
65  {
66  G4cout << "G4LivermoreGammaConversionModelRC is constructed " << G4endl;
67  }
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
72 G4LivermoreGammaConversionModelRC::~G4LivermoreGammaConversionModelRC()
73 {
74  if(IsMaster()) {
75  for(G4int i=0; i<maxZ; ++i) {
76  if(data[i]) {
77  delete data[i];
78  data[i] = 0;
79  }
80  }
81  }
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 void G4LivermoreGammaConversionModelRC::Initialise(
87  const G4ParticleDefinition* particle,
88  const G4DataVector& cuts)
89 {
90  if (verboseLevel > 1)
91  {
92  G4cout << "Calling Initialise() of G4LivermoreGammaConversionModelRC."
93  << G4endl
94  << "Energy range: "
95  << LowEnergyLimit() / MeV << " MeV - "
96  << HighEnergyLimit() / GeV << " GeV"
97  << G4endl;
98  }
99 
100  if(IsMaster())
101  {
102 
103  // Initialise element selector
104 
105  InitialiseElementSelectors(particle, cuts);
106 
107  // Access to elements
108 
109  char* path = std::getenv("G4LEDATA");
110 
111  G4ProductionCutsTable* theCoupleTable =
112  G4ProductionCutsTable::GetProductionCutsTable();
113 
114  G4int numOfCouples = theCoupleTable->GetTableSize();
115 
116  for(G4int i=0; i<numOfCouples; ++i)
117  {
118  const G4Material* material =
119  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
120  const G4ElementVector* theElementVector = material->GetElementVector();
121  G4int nelm = material->GetNumberOfElements();
122 
123  for (G4int j=0; j<nelm; ++j)
124  {
125  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
126  if(Z < 1) { Z = 1; }
127  else if(Z > maxZ) { Z = maxZ; }
128  if(!data[Z]) { ReadData(Z, path); }
129  }
130  }
131  }
132  if(isInitialised) { return; }
133  fParticleChange = GetParticleChangeForGamma();
134  isInitialised = true;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
139 void G4LivermoreGammaConversionModelRC::InitialiseLocal(
140  const G4ParticleDefinition*, G4VEmModel* masterModel)
141 {
142  SetElementSelectors(masterModel->GetElementSelectors());
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 
147 G4double
148 G4LivermoreGammaConversionModelRC::MinPrimaryEnergy(const G4Material*,
149  const G4ParticleDefinition*,
150  G4double)
151 {
152  return lowEnergyLimit;
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
157 void G4LivermoreGammaConversionModelRC::ReadData(size_t Z, const char* path)
158 {
159  if (verboseLevel > 1)
160  {
161  G4cout << "Calling ReadData() of G4LivermoreGammaConversionModelRC"
162  << G4endl;
163  }
164 
165  if(data[Z]) { return; }
166 
167  const char* datadir = path;
168 
169  if(!datadir)
170  {
171  datadir = std::getenv("G4LEDATA");
172  if(!datadir)
173  {
174  G4Exception("G4LivermoreGammaConversionModelRC::ReadData()",
175  "em0006",FatalException,
176  "Environment variable G4LEDATA not defined");
177  return;
178  }
179  }
180 
181  //
182 
183  data[Z] = new G4LPhysicsFreeVector();
184 
185  //
186 
187  std::ostringstream ost;
188  ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
189  std::ifstream fin(ost.str().c_str());
190 
191  if( !fin.is_open())
192  {
193  G4ExceptionDescription ed;
194  ed << "G4LivermoreGammaConversionModelRC data file <" << ost.str().c_str()
195  << "> is not opened!" << G4endl;
196  G4Exception("G4LivermoreGammaConversionModelRC::ReadData()",
197  "em0003",FatalException,
198  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
199  return;
200  }
201 
202  else
203  {
204 
205  if(verboseLevel > 3) { G4cout << "File " << ost.str()
206  << " is opened by G4LivermoreGammaConversionModelRC" << G4endl;}
207 
208  data[Z]->Retrieve(fin, true);
209  }
210 
211  // Activation of spline interpolation
212  data[Z] ->SetSpline(true);
213 
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217 
218 G4double
219 G4LivermoreGammaConversionModelRC::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
220  G4double GammaEnergy,
221  G4double Z, G4double,
222  G4double, G4double)
223 {
224  if (verboseLevel > 1)
225  {
226  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC"
227  << G4endl;
228  }
229 
230  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
231 
232  G4double xs = 0.0;
233 
234  G4int intZ=G4int(Z);
235 
236  if(intZ < 1 || intZ > maxZ) { return xs; }
237 
238  G4LPhysicsFreeVector* pv = data[intZ];
239 
240  // if element was not initialised
241  // do initialisation safely for MT mode
242  if(!pv)
243  {
244  InitialiseForElement(0, intZ);
245  pv = data[intZ];
246  if(!pv) { return xs; }
247  }
248  // x-section is taken from the table
249  xs = pv->Value(GammaEnergy);
250 
251  if(verboseLevel > 0)
252  {
253  G4int n = pv->GetVectorLength() - 1;
254  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
255  << GammaEnergy/MeV << G4endl;
256  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
257  G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
258  G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
259  G4cout << "*********************************************************" << G4endl;
260  }
261 
262  return xs;
263 
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
267 
268 void G4LivermoreGammaConversionModelRC::SampleSecondaries(
269  std::vector<G4DynamicParticle*>* fvect,
270  const G4MaterialCutsCouple* couple,
271  const G4DynamicParticle* aDynamicGamma,
272  G4double, G4double)
273 {
274 
275 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
276 // cross sections with Coulomb correction. A modified version of the random
277 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
278 
279 // Note 1 : Effects due to the breakdown of the Born approximation at low
280 // energy are ignored.
281 // Note 2 : The differential cross section implicitly takes account of
282 // pair creation in both nuclear and atomic electron fields. However triplet
283 // prodution is not generated.
284 
285  if (verboseLevel > 1) {
286  G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC"
287  << G4endl;
288  }
289 
290  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
291  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
292 
293  G4double epsilon ;
294  G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
295 
296  G4double electronTotEnergy = 0.0;
297  G4double positronTotEnergy = 0.0;
298  G4double HardPhotonEnergy = 0.0;
299 
300 
301  // Do it fast if photon energy < 2. MeV
302  if (photonEnergy < smallEnergy )
303  {
304  epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
305  if (G4UniformRand() > 0.5)
306  {
307  electronTotEnergy = (1. - epsilon) * photonEnergy;
308  positronTotEnergy = epsilon * photonEnergy;
309  }
310  else
311  {
312  positronTotEnergy = (1. - epsilon) * photonEnergy;
313  electronTotEnergy = epsilon * photonEnergy;
314  }
315  }
316  else
317  {
318  // Select randomly one element in the current material
319 
320  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
321  const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
322 
323  if (element == 0)
324  {
325  G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
326  << G4endl;
327  return;
328  }
329  G4IonisParamElm* ionisation = element->GetIonisation();
330  if (ionisation == 0)
331  {
332  G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
333  << G4endl;
334  return;
335  }
336 
337  // Extract Coulomb factor for this Element
338  G4double fZ = 8. * (ionisation->GetlogZ3());
339  if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
340 
341  // Limits of the screening variable
342  G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
343  G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ;
344  G4double screenMin = std::min(4.*screenFactor,screenMax) ;
345 
346  // Limits of the energy sampling
347  G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
348  G4double epsilonMin = std::max(epsilon0Local,epsilon1);
349  G4double epsilonRange = 0.5 - epsilonMin ;
350 
351  // Sample the energy rate of the created electron (or positron)
352  G4double screen;
353  G4double gReject ;
354 
355  G4double f10 = ScreenFunction1(screenMin) - fZ;
356  G4double f20 = ScreenFunction2(screenMin) - fZ;
357  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
358  G4double normF2 = std::max(1.5 * f20,0.);
359 
360  // Method for Radiative corrections
361 
362  G4double a=393.3750918, b=115.3070201, c=810.6428451, d=19.96497475, e=1016.874592, f=1.936685510,
363  gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
364  G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287, ff=-117602.840,
365  gg=-721455.467, hh=693957.8635, ii=156266.1085, jj=533209.9347;
366  G4double Rechazo = 0.;
367  G4double logepsMin = log(epsilonMin);
368  G4double NormaRC = a + b*logepsMin + c/logepsMin + d*pow(logepsMin,2.) + e/pow(logepsMin,2.) + f*pow(logepsMin,3.) +
369  gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
370  k/pow(logepsMin,5.);
371 
372  G4double HardPhotonThreshold = 0.08;
373  G4double r1, r2, r3, beta=0, gbeta, sigt = 582.068, sigh, rejet;
374  // , Pi = 2.*acos(0.);
375  G4double cg = (11./2.)/(G4Exp(-11.*HardPhotonThreshold/2.)-G4Exp(-11./2.));
376 
377  r1 = G4UniformRand();
378  sigh = 1028.58*G4Exp(-HardPhotonThreshold/0.09033) + 136.63; // sigma hard
379 
380 
381  if (r1 > 1.- sigh/sigt) {
382  r2 = G4UniformRand();
383  rejet = 0.;
384  while (r2 > rejet) {
385  r3 = G4UniformRand();
386  beta = (-2./11.)*log(G4Exp(-0.08*11./2.)-r3*11./(2.*cg));
387  gbeta = G4Exp(-11.*beta/2.);
388  rejet = fbeta(beta)/(8000.*gbeta);
389  }
390  HardPhotonEnergy = beta * photonEnergy;
391  }
392  else{
393  HardPhotonEnergy = 0.;
394  }
395 
396  photonEnergy -= HardPhotonEnergy;
397 
398  do
399  {
400  do
401  {
402  if (normF1 / (normF1 + normF2) > G4UniformRand() )
403  {
404  epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
405  screen = screenFactor / (epsilon * (1. - epsilon));
406  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
407  }
408  else
409  {
410  epsilon = epsilonMin + epsilonRange * G4UniformRand();
411  screen = screenFactor / (epsilon * (1 - epsilon));
412  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
413  }
414  } while ( gReject < G4UniformRand() );
415 
416 
417 
418  if (G4UniformRand()>0.5) epsilon = (1. - epsilon); // Extención de Epsilon hasta 1.
419 
420  G4double logepsilon = log(epsilon);
421  G4double deltaP_R1 = 1. + (a + b*logepsilon + c/logepsilon + d*pow(logepsilon,2.) + e/pow(logepsilon,2.) +
422  f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
423  j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
424  G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
425  / (1. + bb*logepsilon + dd*pow(logepsilon,2.) + ff*pow(logepsilon,3.) + hh*pow(logepsilon,4.)
426  + jj*pow(logepsilon,5.) ))/100.;
427 
428  if (epsilon <= 0.5)
429  {
430  Rechazo = deltaP_R1/NormaRC;
431  }
432  else
433  {
434  Rechazo = deltaP_R2/NormaRC;
435  }
436  //G4cout << Rechazo << " " << NormaRC << " " << epsilon << G4endl;
437 
438  } while (Rechazo < G4UniformRand() );
439 
440  electronTotEnergy = (1. - epsilon) * photonEnergy;
441  positronTotEnergy = epsilon * photonEnergy;
442 
443  } // End of epsilon sampling
444 
445 
446  // Fix charges randomly
447 
448  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
449  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
450  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
451 
452  G4double u;
453  const G4double a1 = 0.625;
454  G4double a2 = 3. * a1;
455  // G4double d = 27. ;
456 
457  // if (9. / (9. + d) > G4UniformRand())
458  if (0.25 > G4UniformRand())
459  {
460  u = - G4Log(G4UniformRand() * G4UniformRand()) / a1 ;
461  }
462  else
463  {
464  u = - G4Log(G4UniformRand() * G4UniformRand()) / a2 ;
465  }
466 
467  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
468  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
469  G4double phi = twopi * G4UniformRand();
470 
471  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
472  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
473 
474 
475  // Kinematics of the created pair:
476  // the electron and positron are assumed to have a symetric angular
477  // distribution with respect to the Z axis along the parent photon
478 
479  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
480 
481  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
482  electronDirection.rotateUz(photonDirection);
483 
484  G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
485  electronDirection,
486  electronKineEnergy);
487 
488  // The e+ is always created
489  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
490 
491  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
492  positronDirection.rotateUz(photonDirection);
493 
494  // Create G4DynamicParticle object for the particle2
495  G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
496  positronDirection,
497  positronKineEnergy);
498  // Fill output vector
499  fvect->push_back(particle1);
500  fvect->push_back(particle2);
501 
502  if (HardPhotonEnergy > 0.)
503  {
504  G4double thetaHardPhoton = u*electron_mass_c2/HardPhotonEnergy;
505  phi = twopi * G4UniformRand();
506  G4double dxHardP= std::sin(thetaHardPhoton)*std::cos(phi);
507  G4double dyHardP= std::sin(thetaHardPhoton)*std::sin(phi);
508  G4double dzHardP =std::cos(thetaHardPhoton);
509 
510  G4ThreeVector hardPhotonDirection (dxHardP, dyHardP, dzHardP);
511  hardPhotonDirection.rotateUz(photonDirection);
512  G4DynamicParticle* particle3 = new G4DynamicParticle (G4Gamma::Gamma(),
513  hardPhotonDirection,
514  HardPhotonEnergy);
515  fvect->push_back(particle3);
516  }
517 
518  // kill incident photon
519  fParticleChange->SetProposedKineticEnergy(0.);
520  fParticleChange->ProposeTrackStatus(fStopAndKill);
521 
522 }
523 
524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
525 
526 G4double G4LivermoreGammaConversionModelRC::ScreenFunction1(G4double screenVariable)
527 {
528  // Compute the value of the screening function 3*phi1 - phi2
529 
530  G4double value;
531 
532  if (screenVariable > 1.)
533  value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
534  else
535  value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
536 
537  return value;
538 }
539 
540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
541 
542 G4double G4LivermoreGammaConversionModelRC::ScreenFunction2(G4double screenVariable)
543 {
544  // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
545 
546  G4double value;
547 
548  if (screenVariable > 1.)
549  value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
550  else
551  value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
552 
553  return value;
554 }
555 
556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
557 
558 G4double G4LivermoreGammaConversionModelRC::fbeta(G4double x)
559 {
560  // compute the probabililty distribution for hard photon
561  G4double Pi, gamma, eta, d, p1, p2, p3, p4, p5, p6, p7, ffbeta;
562  gamma = (1.-x)*(1.-x)/x;
563  eta = (1.-x)/(1.+x);
564  d = Dilog(1./x)-Dilog(x);
565  Pi = 2.*acos(0.);
566  p1 = -1.*(25528.*pow(gamma,2) + 116044.* gamma +151556.)/105.;
567  p2 = 256.* pow(gamma,3) + 1092.* pow(gamma,2) +1260.*gamma + 420.;
568  p3 = (676.*pow(gamma,3) + 9877.*pow(gamma,2) + 58415.*gamma + 62160.)/105.;
569  p4 = 64.*pow(gamma,3) + 305.*pow(gamma,2) + 475.*gamma + 269. - 276./gamma;
570  p5 = (676.*pow(gamma,3) + 38109.*pow(gamma,2) + 211637.*gamma + 266660. - 53632./gamma)/105.;
571  p6 = 32.*pow(gamma,2) + 416.*gamma + 1310. +1184./gamma;
572  p7 = 128.*pow(gamma,3) + 802.*pow(gamma,2) + 1028.*gamma - 470. - 1184./gamma;
573  ffbeta = (1.-x) * (p1 + p2*Pi*Pi/6. + p3*log(gamma) +
574  p4*pow(log(x),2) + (p5 + p6*log(gamma))*eta*log(x) + p7*d*eta);
575  return ffbeta;
576 }
577 
578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
579 
580 G4double G4LivermoreGammaConversionModelRC::Dilog(G4double y)
581 {
582  G4double fdilog = 0.0;
583  G4double Pi = 2.*acos(0.); // serve?
584  if (y <= 0.5) {
585  fdilog = pow(Pi,2)/6. + (1.-y)*(log(1-y)-1.)+pow((1.-y),2)*((1./2.)*log(1.-y)-1./4.)
586  +pow((1.-y),3)*((1./3.)*log(1.-y)-1./9.)+pow((1.-y),4)*((1./4.)*log(1.-y)-1./16.);
587  }
588  if (0.5 < y && y < 2.) {
589  fdilog = 1.-y+pow((1.-y),2)/4.+pow((1.-y),3)/9.+pow((1.-y),4)/16.+
590  pow((1.-y),5)/25.+pow((1.-y),6)/36.+pow((1.-y),7)/49.;
591  }
592  if (y >= 2.) {
593  fdilog = -pow(log(y),2)/2. - pow(Pi,2)/6. + (log(y)+1.)/y +
594  (log(y)/2.+1./4.)/pow(y,2) + (log(y)/3.+1./9.)/pow(y,3);
595  }
596  return fdilog;
597 }
598 
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
600 
601 #include "G4AutoLock.hh"
602 namespace { G4Mutex LivermoreGammaConversionModelRCMutex = G4MUTEX_INITIALIZER; }
603 
604 void G4LivermoreGammaConversionModelRC::InitialiseForElement(
605  const G4ParticleDefinition*,
606  G4int Z)
607 {
608  G4AutoLock l(&LivermoreGammaConversionModelRCMutex);
609  // G4cout << "G4LivermoreGammaConversionModelRC::InitialiseForElement Z= "
610  // << Z << G4endl;
611  if(!data[Z]) { ReadData(Z); }
612  l.unlock();
613 }
614 
615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......