ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreGammaConversionModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermoreGammaConversionModel.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 G4LivermoreGammaConversionModel (original version)
29 // and G4LivermoreRayleighModel (MT version)
30 
32 #include "G4Electron.hh"
33 #include "G4Positron.hh"
34 #include "G4EmParameters.hh"
36 #include "G4LPhysicsFreeVector.hh"
37 #include "G4PhysicsLogVector.hh"
38 #include "G4ProductionCutsTable.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4Exp.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
54 
56 (const G4ParticleDefinition*, const G4String& nam)
57  : G4VEmModel(nam),fParticleChange(nullptr)
58 {
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 << "G4LivermoreGammaConversionModel is constructed " << G4endl;
67  }
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
73 {
74  if(IsMaster()) {
75  for(G4int i=0; i<maxZ; ++i) {
76  if(data[i]) {
77  delete data[i];
78  data[i] = nullptr;
79  }
80  if(probTriplet[i]) {
81  delete probTriplet[i];
82  probTriplet[i] = nullptr;
83  }
84  }
85  }
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
92  const G4DataVector& cuts)
93 {
94  if (verboseLevel > 1)
95  {
96  G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel."
97  << G4endl
98  << "Energy range: "
99  << LowEnergyLimit() / MeV << " MeV - "
100  << HighEnergyLimit() / GeV << " GeV isMater: " << IsMaster()
101  << G4endl;
102  }
103 
104  if(!fParticleChange) {
106  if(GetTripletModel()) {
108  }
109  }
110  if(GetTripletModel()) { GetTripletModel()->Initialise(particle, cuts); }
111 
112  if(IsMaster())
113  {
114  // Initialise element selector
115  InitialiseElementSelectors(particle, cuts);
116 
117  // Access to elements
118  char* path = std::getenv("G4LEDATA");
119  G4ProductionCutsTable* theCoupleTable =
121 
122  G4int numOfCouples = theCoupleTable->GetTableSize();
123 
124  for(G4int i=0; i<numOfCouples; ++i)
125  {
126  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
127  SetCurrentCouple(couple);
128  const G4Material* mat = couple->GetMaterial();
129  const G4ElementVector* theElementVector = mat->GetElementVector();
130  G4int nelm = mat->GetNumberOfElements();
131 
132  for (G4int j=0; j<nelm; ++j)
133  {
134  G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
135  if(!data[Z]) { ReadData(Z, path); }
136  if(GetTripletModel()) { InitialiseProbability(particle, Z); }
137  }
138  }
139  }
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 
145  const G4ParticleDefinition*, G4VEmModel* masterModel)
146 {
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
152 G4double
154  const G4ParticleDefinition*,
155  G4double)
156 {
157  return lowEnergyLimit;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
162 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path)
163 {
164  if (verboseLevel > 1)
165  {
166  G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
167  << G4endl;
168  }
169 
170  if(data[Z]) { return; }
171 
172  const char* datadir = path;
173 
174  if(!datadir)
175  {
176  datadir = std::getenv("G4LEDATA");
177  if(!datadir)
178  {
179  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
180  "em0006",FatalException,
181  "Environment variable G4LEDATA not defined");
182  return;
183  }
184  }
185  data[Z] = new G4LPhysicsFreeVector();
186  std::ostringstream ost;
187  ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat";
188  std::ifstream fin(ost.str().c_str());
189 
190  if( !fin.is_open())
191  {
193  ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
194  << "> is not opened!" << G4endl;
195  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
196  "em0003",FatalException,
197  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
198  return;
199  }
200  else
201  {
202 
203  if(verboseLevel > 1) { G4cout << "File " << ost.str()
204  << " is opened by G4LivermoreGammaConversionModel" << G4endl;}
205 
206  data[Z]->Retrieve(fin, true);
207  }
208  // Activation of spline interpolation
209  data[Z] ->SetSpline(true);
210 }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
213 
216  G4double GammaEnergy, G4double Z, G4double, G4double, G4double)
217 {
218  if (verboseLevel > 1)
219  {
220  G4cout << "G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom() Z= "
221  << Z << G4endl;
222  }
223 
224  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
225 
226  G4double xs = 0.0;
227 
228  G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ));
229 
230  G4LPhysicsFreeVector* pv = data[intZ];
231 
232  // if element was not initialised
233  // do initialisation safely for MT mode
234  if(!pv)
235  {
236  InitialiseForElement(particle, intZ);
237  pv = data[intZ];
238  if(!pv) { return xs; }
239  }
240  // x-section is taken from the table
241  xs = pv->Value(GammaEnergy);
242 
243  if(verboseLevel > 0)
244  {
245  G4cout << "*** Gamma conversion xs for Z=" << Z << " at energy E(MeV)="
246  << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl;
247  }
248 
249  return xs;
250 
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 
256  std::vector<G4DynamicParticle*>* fvect,
257  const G4MaterialCutsCouple* couple,
258  const G4DynamicParticle* aDynamicGamma,
260 {
261 
262  // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
263  // cross sections with Coulomb correction. A modified version of the random
264  // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
265 
266  // Note 1 : Effects due to the breakdown of the Born approximation at low
267  // energy are ignored.
268  // Note 2 : The differential cross section implicitly takes account of
269  // pair creation in both nuclear and atomic electron fields. However triplet
270  // prodution is not generated.
271 
272  if (verboseLevel > 1) {
273  G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel"
274  << G4endl;
275  }
276 
277  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
278  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
279 
280  G4double epsilon ;
281  G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
282 
283  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
284 
285  // Do it fast if photon energy < 2. MeV
286  static const G4double smallEnergy = 2.*CLHEP::MeV;
287  if (photonEnergy < smallEnergy )
288  {
289  epsilon = epsilon0Local + (0.5 - epsilon0Local) * rndmEngine->flat();
290  }
291  else
292  {
293  // Select randomly one element in the current material
294 
295  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
296  const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
297  G4int Z = element->GetZasInt();
298 
299  // triplet production
300  if(GetTripletModel()) {
301  if(!probTriplet[Z]) { InitialiseForElement(particle, Z); }
302  /*
303  G4cout << "Liv: E= " << photonEnergy
304  << " prob= " << probTriplet[Z]->Value(photonEnergy)
305  << G4endl;
306  */
307  if(probTriplet[Z] &&
308  rndmEngine->flat() < probTriplet[Z]->Value(photonEnergy)) {
309  GetTripletModel()->SampleSecondaries(fvect, couple, aDynamicGamma);
310  return;
311  }
312  }
313 
314  G4IonisParamElm* ionisation = element->GetIonisation();
315 
316  // Extract Coulomb factor for this Element
317  G4double fZ = 8. * (ionisation->GetlogZ3());
318  static const G4double midEnergy = 50.*CLHEP::MeV;
319  if (photonEnergy > midEnergy) { fZ += 8. * (element->GetfCoulomb()); }
320 
321  // Limits of the screening variable
322  G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3());
323  G4double screenMax = G4Exp((42.24 - fZ)/8.368) + 0.952;
324  G4double screenMin = std::min(4.*screenFactor,screenMax);
325 
326  // Limits of the energy sampling
327  G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
328  G4double epsilonMin = std::max(epsilon0Local,epsilon1);
329  G4double epsilonRange = 0.5 - epsilonMin ;
330 
331  // Sample the energy rate of the created electron (or positron)
332  G4double screen;
333  G4double gReject;
334 
335  G4double f10 = ScreenFunction1(screenMin) - fZ;
336  G4double f20 = ScreenFunction2(screenMin) - fZ;
337  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
338  G4double normF2 = std::max(1.5 * f20,0.);
339 
340  do
341  {
342  if (normF1 > (normF1 + normF2)*rndmEngine->flat() )
343  {
344  epsilon = 0.5 - epsilonRange *G4Exp(G4Log(rndmEngine->flat())/3.);
345  screen = screenFactor / (epsilon * (1. - epsilon));
346  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
347  }
348  else
349  {
350  epsilon = epsilonMin + epsilonRange * rndmEngine->flat();
351  screen = screenFactor / (epsilon * (1 - epsilon));
352  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
353  }
354  } while ( gReject < rndmEngine->flat() );
355 
356  } // End of epsilon sampling
357 
358  // Fix charges randomly
359 
360  G4double electronTotEnergy;
361  G4double positronTotEnergy;
362 
363  if (rndmEngine->flat() > 0.5)
364  {
365  electronTotEnergy = (1. - epsilon) * photonEnergy;
366  positronTotEnergy = epsilon * photonEnergy;
367  }
368  else
369  {
370  positronTotEnergy = (1. - epsilon) * photonEnergy;
371  electronTotEnergy = epsilon * photonEnergy;
372  }
373 
374  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
375  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
376  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
377 
378  static const G4double a1 = 1.6;
379  static const G4double a2 = 0.5333333333;
380  G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat());
381  G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2;
382 
383  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
384  G4double sinte = std::sin(thetaEle);
385  G4double coste = std::cos(thetaEle);
386 
387  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
388  G4double sintp = std::sin(thetaPos);
389  G4double costp = std::cos(thetaPos);
390 
391  G4double phi = twopi * rndmEngine->flat();
392  G4double sinp = std::sin(phi);
393  G4double cosp = std::cos(phi);
394 
395  // Kinematics of the created pair:
396  // the electron and positron are assumed to have a symetric angular
397  // distribution with respect to the Z axis along the parent photon
398 
399  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
400 
401  G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
402  electronDirection.rotateUz(photonDirection);
403 
405  electronDirection,
406  electronKineEnergy);
407 
408  // The e+ is always created
409  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
410 
411  G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
412  positronDirection.rotateUz(photonDirection);
413 
414  // Create G4DynamicParticle object for the particle2
416  positronDirection,
417  positronKineEnergy);
418  // Fill output vector
419  fvect->push_back(particle1);
420  fvect->push_back(particle2);
421 
422  // kill incident photon
425 
426 }
427 
428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
429 
430 #include "G4AutoLock.hh"
431 namespace { G4Mutex LivermoreGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
432 
434  const G4ParticleDefinition* part,
435  G4int Z)
436 {
438  G4AutoLock l(&LivermoreGammaConversionModelMutex);
439  // G4cout << "G4LivermoreGammaConversionModel::InitialiseForElement Z= "
440  // << Z << G4endl;
441  if(!data[Z]) { ReadData(Z); }
442  if(GetTripletModel() && !probTriplet[Z]) { InitialiseProbability(part, Z); }
443  l.unlock();
444 }
445 
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
447 
450 {
451  if(!probTriplet[Z]) {
453  : nullptr;
454  if(0 == nbinsTriplet) {
459  nbinsTriplet = std::max(3,
460  (G4int)(nbins*G4Log(tripletHighEnergy/tripletLowEnergy)/(6*G4Log(10.))));
461  }
462  /*
463  G4cout << "G4LivermoreGammaConversionModel::InitialiseProbability Z= "
464  << Z << " Nbin= " << nbinsTriplet
465  << " Emin(MeV)= " << tripletLowEnergy
466  << " Emax(MeV)= " << tripletHighEnergy << G4endl;
467  */
468  probTriplet[Z] =
470  probTriplet[Z]->SetSpline(true);
471  G4double zz = (G4double)Z;
472  // loop over bins
473  for(G4int j=0; j<=nbinsTriplet; ++j) {
474  G4double e = (probTriplet[Z])->Energy(j);
475  SetupForMaterial(part, mat, e);
477  G4double tcross =
479  tcross = (0.0 < cross) ? tcross/cross : 0.0;
480  (probTriplet[Z])->PutValue(j, tcross);
481  //G4cout << j << ". E= " << e << " prob= " << tcross << G4endl;
482  }
483  }
484 }
485 
486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......