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