ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BoldyshevTripletModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BoldyshevTripletModel.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 G4BoldyshevTripletModel (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  :G4VEmModel(nam),smallEnergy(4.*MeV)
48 {
49  fParticleChange = nullptr;
50 
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 << "G4BoldyshevTripletModel 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] = nullptr;
75  }
76  }
77  }
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
83  const G4DataVector&)
84 {
85  if (verboseLevel > 1)
86  {
87  G4cout << "Calling Initialise() of G4BoldyshevTripletModel."
88  << G4endl
89  << "Energy range: "
90  << LowEnergyLimit() / MeV << " MeV - "
91  << HighEnergyLimit() / GeV << " GeV isMaster: " << IsMaster()
92  << G4endl;
93  }
94  // compute values only once
98  G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2;
99  G4double t = 0.5*G4Log(momentumThreshold_N +
100  std::sqrt(momentumThreshold_N*momentumThreshold_N + 1.0));
101  //G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
102  G4double sinht = std::sinh(t);
103  G4double cosht = std::cosh(t);
104  G4double logsinht = G4Log(2.*sinht);
105  G4double J1 = 0.5*(t*cosht/sinht - logsinht);
106  G4double J2 = (-2./3.)*logsinht + t*cosht/sinht
107  + (sinht - t*cosht*cosht*cosht)/(3.*sinht*sinht*sinht);
108 
109  xb = 2.*(J1-J2)/J1;
110  xn = 1. - xb/6.;
111 
112  if(IsMaster())
113  {
114  // Access to elements
115  char* path = std::getenv("G4LEDATA");
116 
117  G4ProductionCutsTable* theCoupleTable =
119 
120  G4int numOfCouples = theCoupleTable->GetTableSize();
121 
122  for(G4int i=0; i<numOfCouples; ++i)
123  {
124  const G4Material* material =
125  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
126  const G4ElementVector* theElementVector = material->GetElementVector();
127  G4int nelm = material->GetNumberOfElements();
128 
129  for (G4int j=0; j<nelm; ++j)
130  {
131  G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
132  if(!data[Z]) { ReadData(Z, path); }
133  }
134  }
135  }
136  if(!fParticleChange) {
138  }
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
143 G4double
145  const G4ParticleDefinition*,
146  G4double)
147 {
148  return lowEnergyLimit;
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 
153 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path)
154 {
155  if (verboseLevel > 1)
156  {
157  G4cout << "Calling ReadData() of G4BoldyshevTripletModel"
158  << G4endl;
159  }
160 
161  if(data[Z]) { return; }
162 
163  const char* datadir = path;
164 
165  if(!datadir)
166  {
167  datadir = std::getenv("G4LEDATA");
168  if(!datadir)
169  {
170  G4Exception("G4BoldyshevTripletModel::ReadData()",
171  "em0006",FatalException,
172  "Environment variable G4LEDATA not defined");
173  return;
174  }
175  }
176 
177  data[Z] = new G4LPhysicsFreeVector();
178  std::ostringstream ost;
179  ost << datadir << "/livermore/tripdata/pp-trip-cs-" << Z <<".dat";
180  std::ifstream fin(ost.str().c_str());
181 
182  if( !fin.is_open())
183  {
185  ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str()
186  << "> is not opened!" << G4endl;
187  G4Exception("G4BoldyshevTripletModel::ReadData()",
188  "em0003",FatalException,
189  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
190  return;
191  }
192 
193  else
194  {
195 
196  if(verboseLevel > 3) { G4cout << "File " << ost.str()
197  << " is opened by G4BoldyshevTripletModel" << G4endl;}
198 
199  data[Z]->Retrieve(fin, true);
200  }
201 
202  // Activation of spline interpolation
203  data[Z]->SetSpline(true);
204 }
205 
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 
209  const G4ParticleDefinition* part,
210  G4double GammaEnergy, G4double Z, G4double, G4double, G4double)
211 {
212  if (verboseLevel > 1)
213  {
214  G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
215  << G4endl;
216  }
217 
218  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
219 
220  G4double xs = 0.0;
221  G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ));
222  G4LPhysicsFreeVector* pv = data[intZ];
223 
224  // if element was not initialised
225  // do initialisation safely for MT mode
226  if(!pv)
227  {
228  InitialiseForElement(part, intZ);
229  pv = data[intZ];
230  if(!pv) { return xs; }
231  }
232  // x-section is taken from the table
233  xs = pv->Value(GammaEnergy);
234 
235  if(verboseLevel > 1)
236  {
237  G4cout << "*** Triplet conversion xs for Z=" << Z << " at energy E(MeV)="
238  << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl;
239  }
240  return xs;
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 
246  std::vector<G4DynamicParticle*>* fvect,
247  const G4MaterialCutsCouple* /*couple*/,
248  const G4DynamicParticle* aDynamicGamma,
250 {
251 
252  // The energies of the secondary particles are sampled using
253  // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
254  if (verboseLevel > 1) {
255  G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel"
256  << G4endl;
257  }
258 
259  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
260  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
261 
263 
264  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
265 
266  // recoil electron thould be 3d particle
267  G4DynamicParticle* particle3 = nullptr;
268  static const G4double costlim = std::cos(4.47*CLHEP::pi/180.);
269 
270  G4double loga, f1_re, greject, cost;
271  G4double cosThetaMax = (energyThreshold - electron_mass_c2
274  if (cosThetaMax > 1.) {
275  //G4cout << "G4BoldyshevTripletModel::SampleSecondaries: ERROR cosThetaMax= "
276  // << cosThetaMax << G4endl;
277  cosThetaMax = 1.0;
278  }
279 
280  G4double logcostm = G4Log(cosThetaMax);
281  G4int nn = 0;
282  do {
283  cost = G4Exp(logcostm*rndmEngine->flat());
284  G4double are = 1./(14.*cost*cost);
285  G4double bre = (1.-5.*cost*cost)/(2.*cost);
286  loga = G4Log((1.+ cost)/(1.- cost));
287  f1_re = 1. - bre*loga;
288  greject = (cost < costlim) ? are*f1_re : 1.0;
289  // G4cout << nn << ". step of the 1st loop greject= " << greject << G4endl;
290  ++nn;
291  } while(greject < rndmEngine->flat());
292 
293  // Calculo de phi - elecron de recoil
294  G4double sint2 = (1. - cost)*(1. + cost);
295  G4double fp = 1. - sint2*loga/(2.*cost) ;
296  G4double rt, phi_re;
297  nn = 0;
298  do {
299  phi_re = twopi*rndmEngine->flat();
300  rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
301  //G4cout << nn << ". step of the 2nd loop greject= " << rt << G4endl;
302  ++nn;
303  } while(rt < rndmEngine->flat());
304 
305  // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
306  G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
308 
309  G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
310  G4double ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
311 
312  if(ener_re >= energyThreshold)
313  {
314  G4double electronRKineEnergy = ener_re - electron_mass_c2;
315  G4double sint = std::sqrt(sint2);
316  G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
317  electronRDirection.rotateUz(photonDirection);
318  particle3 = new G4DynamicParticle (G4Electron::Electron(),
319  electronRDirection,
320  electronRKineEnergy);
321  }
322  else
323  {
324  // deposito la energia ener_re - electron_mass_c2
325  // G4cout << "electron de retroceso " << ener_re << G4endl;
326  fParticleChange->ProposeLocalEnergyDeposit(std::max(0.0, ener_re - electron_mass_c2));
327  ener_re = 0.0;
328  }
329 
330  // Depaola (2004) suggested distribution for e+e- energy
331  // VI: very suspect that 1 random number is not enough
332  // and sampling below is not correct - should be fixed
333  G4double re = rndmEngine->flat();
334 
335  G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
336  G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + xb + 2*a)*xb*xb)/3.);
337  epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
338 
339  G4double photonEnergy1 = photonEnergy - ener_re ;
340  // resto al foton la energia del electron de retro.
341  G4double positronTotEnergy = std::max(epsilon*photonEnergy1, electron_mass_c2);
342  G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
343 
344  static const G4double a1 = 1.6;
345  static const G4double a2 = 0.5333333333;
346  G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat());
347  G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2;
348 
349  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
350  G4double sinte = std::sin(thetaEle);
351  G4double coste = std::cos(thetaEle);
352 
353  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
354  G4double sintp = std::sin(thetaPos);
355  G4double costp = std::cos(thetaPos);
356 
357  G4double phi = twopi * rndmEngine->flat();
358  G4double sinp = std::sin(phi);
359  G4double cosp = std::cos(phi);
360 
361  // Kinematics of the created pair:
362  // the electron and positron are assumed to have a symetric angular
363  // distribution with respect to the Z axis along the parent photon
364 
365  G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
366 
367  G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
368  electronDirection.rotateUz(photonDirection);
369 
371  electronDirection,
372  electronKineEnergy);
373 
374  G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
375 
376  G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
377  positronDirection.rotateUz(photonDirection);
378 
379  // Create G4DynamicParticle object for the particle2
381  positronDirection, positronKineEnergy);
382  // Fill output vector
383 
384  fvect->push_back(particle1);
385  fvect->push_back(particle2);
386 
387  if(particle3) { fvect->push_back(particle3); }
388 
389  // kill incident photon
392 }
393 
394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
395 
396 #include "G4AutoLock.hh"
397 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; }
398 
400  const G4ParticleDefinition*, G4int Z)
401 {
402  G4AutoLock l(&BoldyshevTripletModelMutex);
403  // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= "
404  // << Z << G4endl;
405  if(!data[Z]) { ReadData(Z); }
406  l.unlock();
407 }
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......