Geant4  10.02.p02
G4BoldyshevTripletModel.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 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),isInitialised(false),smallEnergy(4.*MeV)
48 {
49  fParticleChange = 0;
50 
51  lowEnergyLimit = 4.0*electron_mass_c2;
52 
53  verboseLevel= 0;
54  // Verbosity scale for debugging purposes:
55  // 0 = nothing
56  // 1 = calculation of cross sections, file openings...
57  // 2 = entering in methods
58 
59  if(verboseLevel > 0)
60  {
61  G4cout << "G4BoldyshevTripletModel is constructed " << G4endl;
62  }
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68 {
69  if(IsMaster()) {
70  for(G4int i=0; i<maxZ; ++i) {
71  if(data[i]) {
72  delete data[i];
73  data[i] = 0;
74  }
75  }
76  }
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82  const G4ParticleDefinition* particle,
83  const G4DataVector& cuts)
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"
92  << G4endl;
93  }
94 
95  if(IsMaster())
96  {
97 
98  // Initialise element selector
99 
100  InitialiseElementSelectors(particle, cuts);
101 
102  // Access to elements
103 
104  char* path = getenv("G4LEDATA");
105 
106  G4ProductionCutsTable* theCoupleTable =
108 
109  G4int numOfCouples = theCoupleTable->GetTableSize();
110 
111  for(G4int i=0; i<numOfCouples; ++i)
112  {
113  const G4Material* material =
114  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
115  const G4ElementVector* theElementVector = material->GetElementVector();
116  G4int nelm = material->GetNumberOfElements();
117 
118  for (G4int j=0; j<nelm; ++j)
119  {
120  G4int Z = (G4int)(*theElementVector)[j]->GetZ();
121  if(Z < 1) { Z = 1; }
122  else if(Z > maxZ) { Z = maxZ; }
123  if(!data[Z]) { ReadData(Z, path); }
124  }
125  }
126  }
127  if(isInitialised) { return; }
129  isInitialised = true;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 
135  const G4ParticleDefinition*, G4VEmModel* masterModel)
136 {
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141 
142 G4double
144  const G4ParticleDefinition*,
145  G4double)
146 {
147  return lowEnergyLimit;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
152 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path)
153 {
154  if (verboseLevel > 1)
155  {
156  G4cout << "Calling ReadData() of G4BoldyshevTripletModel"
157  << G4endl;
158  }
159 
160  if(data[Z]) { return; }
161 
162  const char* datadir = path;
163 
164  if(!datadir)
165  {
166  datadir = getenv("G4LEDATA");
167  if(!datadir)
168  {
169  G4Exception("G4BoldyshevTripletModel::ReadData()",
170  "em0006",FatalException,
171  "Environment variable G4LEDATA not defined");
172  return;
173  }
174  }
175 
176  //
177 
178  data[Z] = new G4LPhysicsFreeVector();
179 
180  //
181 
182  std::ostringstream ost;
183  ost << datadir << "livermore/tripdata/pp-trip-cs-" << Z <<".dat";
184  std::ifstream fin(ost.str().c_str());
185 
186  if( !fin.is_open())
187  {
189  ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str()
190  << "> is not opened!" << G4endl;
191  G4Exception("G4BoldyshevTripletModel::ReadData()",
192  "em0003",FatalException,
193  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
194  return;
195  }
196 
197  else
198  {
199 
200  if(verboseLevel > 3) { G4cout << "File " << ost.str()
201  << " is opened by G4BoldyshevTripletModel" << G4endl;}
202 
203  data[Z]->Retrieve(fin, true);
204  }
205 
206  // Activation of spline interpolation
207  data[Z] ->SetSpline(true);
208 
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
212 
213 G4double
215  G4double GammaEnergy,
216  G4double Z, G4double,
218 {
219  if (verboseLevel > 1)
220  {
221  G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
222  << G4endl;
223  }
224 
225  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
226 
227  G4double xs = 0.0;
228 
229  G4int intZ=G4int(Z);
230 
231  if(intZ < 1 || intZ > maxZ) { return xs; }
232 
233  G4LPhysicsFreeVector* pv = data[intZ];
234 
235  // if element was not initialised
236  // do initialisation safely for MT mode
237  if(!pv)
238  {
239  InitialiseForElement(0, intZ);
240  pv = data[intZ];
241  if(!pv) { return xs; }
242  }
243  // x-section is taken from the table
244  xs = pv->Value(GammaEnergy);
245 
246  if(verboseLevel > 0)
247  {
248  G4int n = pv->GetVectorLength() - 1;
249  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
250  << GammaEnergy/MeV << G4endl;
251  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
252  G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
253  G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
254  G4cout << "*********************************************************" << G4endl;
255  }
256 
257  return xs;
258 
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
262 
264  std::vector<G4DynamicParticle*>* fvect,
265  const G4MaterialCutsCouple* /*couple*/,
266  const G4DynamicParticle* aDynamicGamma,
268 {
269 
270  // The energies of the secondary particles are sampled using // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
271 
272  if (verboseLevel > 1) {
273  G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel"
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 
284  G4double p0 = electron_mass_c2;
285  G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
286  G4double ener_re=0., theta_re, phi_re, phi;
287 
288  // Calculo de theta - elecron de recoil
289 
290  G4double energyThreshold = sqrt(2.)*electron_mass_c2; // -> momentumThreshold_N = 1
291 
292  energyThreshold = 1.1*electron_mass_c2;
293  // G4cout << energyThreshold << G4endl;
294 
295 
296  G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2); // momentun in MeV/c unit
297  G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; // momentun in mc unit
298 
299  // Calculation of recoil electron production
300 
301 
302  G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ;
303  G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
304  G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0);
305  G4double recoilProb = G4UniformRand();
306  //G4cout << "SIGMA TOT " << SigmaTot << " " << "SigmaQ " << SigmaQ << " " << SigmaQ/SigmaTot << " " << recoilProb << G4endl;
307 
308 
309  if (recoilProb >= SigmaQ/SigmaTot) // create electron recoil
310  {
311 
312  G4double cosThetaMax = ( ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
313  ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
314 
315 
316  if (cosThetaMax > 1) G4cout << "ERRORE " << G4endl;
317 
318  G4double r1;
319  G4double r2;
320  G4double are, bre, loga, f1_re, greject, cost;
321 
322  do {
323  r1 = G4UniformRand();
324  r2 = G4UniformRand();
325  // cost = (pow(4./enern,0.5*r1)) ;
326 
327  cost = pow(cosThetaMax,r1);
328  theta_re = acos(cost);
329  are = 1./(14.*cost*cost);
330  bre = (1.-5.*cost*cost)/(2.*cost);
331  loga = log((1.+ cost)/(1.- cost));
332  f1_re = 1. - bre*loga;
333 
334  if ( theta_re >= 4.47*CLHEP::pi/180.)
335  {
336  greject = are*f1_re;
337  } else {
338  greject = 1. ;
339  }
340  } while(greject < r2);
341 
342  // Calculo de phi - elecron de recoil
343 
344  G4double r3, r4, rt;
345 
346  do {
347  r3 = G4UniformRand();
348  r4 = G4UniformRand();
349  phi_re = twopi*r3 ;
350  G4double sint2 = 1. - cost*cost ;
351  G4double fp = 1. - sint2*loga/(2.*cost) ;
352  rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
353 
354  } while(rt < r4);
355 
356  // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
357 
358  G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
359 
360  G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
361  + (S - electron_mass_c2*electron_mass_c2)
362  *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
363  ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
364 
365  // New Recoil energy calculation
366 
367  G4double momentum_recoil = 2* (electron_mass_c2) * (std::cos(theta_re)/(std::sin(phi_re)*std::sin(phi_re)));
368  G4double ener_recoil = sqrt( momentum_recoil*momentum_recoil + electron_mass_c2*electron_mass_c2);
369  ener_re = ener_recoil;
370 
371  // G4cout << "electron de retroceso " << ener_re << " " << theta_re << " " << phi_re << G4endl;
372 
373  // Recoil electron creation
374  G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
375 
376  G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
377 
378  G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
379  electronRDirection.rotateUz(photonDirection);
380 
382  electronRDirection,
383  electronRKineEnergy);
384  fvect->push_back(particle3);
385 
386  }
387  else
388  {
389  // deposito la energia ener_re - electron_mass_c2
390  // G4cout << "electron de retroceso " << ener_re << G4endl;
391 
392  fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
393  }
394 
395 
396  // Depaola (2004) suggested distribution for e+e- energy
397 
398  // G4double t = 0.5*asinh(momentumThreshold_N);
399  G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
400  //G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
401  G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
402  G4double J2 = (-2./3.)*log(2.*sinh(t)) + t*cosh(t)/sinh(t) + (sinh(t)-t*pow(cosh(t),3))/(3.*pow(sinh(t),3));
403  G4double b = 2.*(J1-J2)/J1;
404 
405  G4double n = 1 - b/6.;
406  G4double re=0.;
407  re = G4UniformRand();
408  G4double a = 0.;
409 
410  G4double b1 = 16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) +
411  6.*pow(b,2.)*re*n;
412  a = pow((b1/b),0.5);
413  G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
414  epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
415 
416  G4double photonEnergy1 = photonEnergy - ener_re ; // resto al foton la energia del electron de retro.
417  positronTotEnergy = epsilon*photonEnergy1;
418  electronTotEnergy = photonEnergy1 - positronTotEnergy; // temporarly
419 
420  G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
421  electron_mass_c2*electron_mass_c2) ;
422  G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
423  electron_mass_c2*electron_mass_c2) ;
424 
425  thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
426  thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
427  phi = twopi * G4UniformRand();
428 
429  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
430  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
431 
432  // Kinematics of the created pair:
433 
434  // the electron and positron are assumed to have a symetric angular
435  // distribution with respect to the Z axis along the parent photon
436 
437 
438  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
439 
440 
441  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
442  electronDirection.rotateUz(photonDirection);
443 
445  electronDirection,
446  electronKineEnergy);
447 
448  // The e+ is always created (even with kinetic energy = 0) for further annihilation
449 
450  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
451 
452  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
453  positronDirection.rotateUz(photonDirection);
454 
455  // Create G4DynamicParticle object for the particle2
456 
458  positronDirection, positronKineEnergy);
459  // Fill output vector
460 
461  fvect->push_back(particle1);
462  fvect->push_back(particle2);
463 
464 
465  // kill incident photon
468 
469 
470 
471 }
472 
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 
475 #include "G4AutoLock.hh"
476 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; }
477 
479  const G4ParticleDefinition*,
480  G4int Z)
481 {
482  G4AutoLock l(&BoldyshevTripletModelMutex);
483  // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= "
484  // << Z << G4endl;
485  if(!data[Z]) { ReadData(Z); }
486  l.unlock();
487 }
488 
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static const double MeV
Definition: G4SIunits.hh:211
G4BoldyshevTripletModel(const G4ParticleDefinition *p=0, const G4String &nam="BoldyshevTripletConversion")
std::vector< G4Element * > G4ElementVector
double S(double temp)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:148
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
void ReadData(size_t Z, const char *path=0)
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
G4double a
Definition: TRTMaterials.hh:39
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma * fParticleChange
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetSpline(G4bool)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
Definition: G4SIunits.hh:75
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static const double GeV
Definition: G4SIunits.hh:214
const G4int n
static const G4double c1
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:802
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4LPhysicsFreeVector * data[100]
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
static G4ProductionCutsTable * GetProductionCutsTable()
G4int G4Mutex
Definition: G4Threading.hh:173
static G4Positron * Positron()
Definition: G4Positron.cc:94
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:810
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
static const G4double b1
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ThreeVector G4ParticleMomentum
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
double epsilon(double density, double temperature)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
const G4Material * GetMaterial() const