Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 
91  const G4ParticleDefinition* particle,
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 = getenv("G4LEDATA");
119 
120  G4ProductionCutsTable* theCoupleTable =
122 
123  G4int numOfCouples = theCoupleTable->GetTableSize();
124 
125  for(G4int i=0; i<numOfCouples; ++i)
126  {
127  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
128  SetCurrentCouple(couple);
129  const G4Material* mat = couple->GetMaterial();
130  const G4ElementVector* theElementVector = mat->GetElementVector();
131  G4int nelm = mat->GetNumberOfElements();
132 
133  for (G4int j=0; j<nelm; ++j)
134  {
135  G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
136  if(!data[Z]) { ReadData(Z, path); }
137  if(GetTripletModel()) { InitialiseProbability(particle, Z); }
138  }
139  }
140  }
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
146  const G4ParticleDefinition*, G4VEmModel* masterModel)
147 {
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 
153 G4double
155  const G4ParticleDefinition*,
156  G4double)
157 {
158  return lowEnergyLimit;
159 }
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 
163 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path)
164 {
165  if (verboseLevel > 1)
166  {
167  G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
168  << G4endl;
169  }
170 
171  if(data[Z]) { return; }
172 
173  const char* datadir = path;
174 
175  if(!datadir)
176  {
177  datadir = getenv("G4LEDATA");
178  if(!datadir)
179  {
180  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
181  "em0006",FatalException,
182  "Environment variable G4LEDATA not defined");
183  return;
184  }
185  }
186  data[Z] = new G4LPhysicsFreeVector();
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  {
194  ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str()
195  << "> is not opened!" << G4endl;
196  G4Exception("G4LivermoreGammaConversionModel::ReadData()",
197  "em0003",FatalException,
198  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
199  return;
200  }
201  else
202  {
203 
204  if(verboseLevel > 1) { G4cout << "File " << ost.str()
205  << " is opened by G4LivermoreGammaConversionModel" << G4endl;}
206 
207  data[Z]->Retrieve(fin, true);
208  }
209  // Activation of spline interpolation
210  data[Z] ->SetSpline(true);
211 }
212 
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214 
216  const G4ParticleDefinition* particle,
217  G4double GammaEnergy, G4double Z, G4double, G4double, G4double)
218 {
219  if (verboseLevel > 1)
220  {
221  G4cout << "G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom() Z= "
222  << Z << G4endl;
223  }
224 
225  if (GammaEnergy < lowEnergyLimit) { return 0.0; }
226 
227  G4double xs = 0.0;
228 
229  G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ));
230 
231  G4LPhysicsFreeVector* pv = data[intZ];
232 
233  // if element was not initialised
234  // do initialisation safely for MT mode
235  if(!pv)
236  {
237  InitialiseForElement(particle, intZ);
238  pv = data[intZ];
239  if(!pv) { return xs; }
240  }
241  // x-section is taken from the table
242  xs = pv->Value(GammaEnergy);
243 
244  if(verboseLevel > 0)
245  {
246  G4cout << "*** Gamma conversion xs for Z=" << Z << " at energy E(MeV)="
247  << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl;
248  }
249 
250  return xs;
251 
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
257  std::vector<G4DynamicParticle*>* fvect,
258  const G4MaterialCutsCouple* couple,
259  const G4DynamicParticle* aDynamicGamma,
261 {
262 
263  // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
264  // cross sections with Coulomb correction. A modified version of the random
265  // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
266 
267  // Note 1 : Effects due to the breakdown of the Born approximation at low
268  // energy are ignored.
269  // Note 2 : The differential cross section implicitly takes account of
270  // pair creation in both nuclear and atomic electron fields. However triplet
271  // prodution is not generated.
272 
273  if (verboseLevel > 1) {
274  G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel"
275  << G4endl;
276  }
277 
278  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
279  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
280 
281  G4double epsilon ;
282  G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
283 
284  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
285 
286  // Do it fast if photon energy < 2. MeV
287  static const G4double smallEnergy = 2.*CLHEP::MeV;
288  if (photonEnergy < smallEnergy )
289  {
290  epsilon = epsilon0Local + (0.5 - epsilon0Local) * rndmEngine->flat();
291  }
292  else
293  {
294  // Select randomly one element in the current material
295 
296  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
297  const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
298  G4int Z = element->GetZasInt();
299 
300  // triplet production
301  if(GetTripletModel()) {
302  if(!probTriplet[Z]) { InitialiseForElement(particle, Z); }
303  /*
304  G4cout << "Liv: E= " << photonEnergy
305  << " prob= " << probTriplet[Z]->Value(photonEnergy)
306  << G4endl;
307  */
308  if(probTriplet[Z] &&
309  rndmEngine->flat() < probTriplet[Z]->Value(photonEnergy)) {
310  GetTripletModel()->SampleSecondaries(fvect, couple, aDynamicGamma);
311  return;
312  }
313  }
314 
315  G4IonisParamElm* ionisation = element->GetIonisation();
316 
317  // Extract Coulomb factor for this Element
318  G4double fZ = 8. * (ionisation->GetlogZ3());
319  static const G4double midEnergy = 50.*CLHEP::MeV;
320  if (photonEnergy > midEnergy) { fZ += 8. * (element->GetfCoulomb()); }
321 
322  // Limits of the screening variable
323  G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3());
324  G4double screenMax = G4Exp((42.24 - fZ)/8.368) + 0.952;
325  G4double screenMin = std::min(4.*screenFactor,screenMax);
326 
327  // Limits of the energy sampling
328  G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
329  G4double epsilonMin = std::max(epsilon0Local,epsilon1);
330  G4double epsilonRange = 0.5 - epsilonMin ;
331 
332  // Sample the energy rate of the created electron (or positron)
333  G4double screen;
334  G4double gReject;
335 
336  G4double f10 = ScreenFunction1(screenMin) - fZ;
337  G4double f20 = ScreenFunction2(screenMin) - fZ;
338  G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
339  G4double normF2 = std::max(1.5 * f20,0.);
340 
341  do
342  {
343  if (normF1 > (normF1 + normF2)*rndmEngine->flat() )
344  {
345  epsilon = 0.5 - epsilonRange *G4Exp(G4Log(rndmEngine->flat())/3.);
346  screen = screenFactor / (epsilon * (1. - epsilon));
347  gReject = (ScreenFunction1(screen) - fZ) / f10 ;
348  }
349  else
350  {
351  epsilon = epsilonMin + epsilonRange * rndmEngine->flat();
352  screen = screenFactor / (epsilon * (1 - epsilon));
353  gReject = (ScreenFunction2(screen) - fZ) / f20 ;
354  }
355  } while ( gReject < rndmEngine->flat() );
356 
357  } // End of epsilon sampling
358 
359  // Fix charges randomly
360 
361  G4double electronTotEnergy;
362  G4double positronTotEnergy;
363 
364  if (rndmEngine->flat() > 0.5)
365  {
366  electronTotEnergy = (1. - epsilon) * photonEnergy;
367  positronTotEnergy = epsilon * photonEnergy;
368  }
369  else
370  {
371  positronTotEnergy = (1. - epsilon) * photonEnergy;
372  electronTotEnergy = epsilon * photonEnergy;
373  }
374 
375  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
376  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
377  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
378 
379  static const G4double a1 = 1.6;
380  static const G4double a2 = 0.5333333333;
381  G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat());
382  G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2;
383 
384  G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
385  G4double sinte = std::sin(thetaEle);
386  G4double coste = std::cos(thetaEle);
387 
388  G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
389  G4double sintp = std::sin(thetaPos);
390  G4double costp = std::cos(thetaPos);
391 
392  G4double phi = twopi * rndmEngine->flat();
393  G4double sinp = std::sin(phi);
394  G4double cosp = std::cos(phi);
395 
396  // Kinematics of the created pair:
397  // the electron and positron are assumed to have a symetric angular
398  // distribution with respect to the Z axis along the parent photon
399 
400  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
401 
402  G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
403  electronDirection.rotateUz(photonDirection);
404 
406  electronDirection,
407  electronKineEnergy);
408 
409  // The e+ is always created
410  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
411 
412  G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
413  positronDirection.rotateUz(photonDirection);
414 
415  // Create G4DynamicParticle object for the particle2
417  positronDirection,
418  positronKineEnergy);
419  // Fill output vector
420  fvect->push_back(particle1);
421  fvect->push_back(particle2);
422 
423  // kill incident photon
426 
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430 
431 #include "G4AutoLock.hh"
432 namespace { G4Mutex LivermoreGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
433 
435  const G4ParticleDefinition* part,
436  G4int Z)
437 {
439  G4AutoLock l(&LivermoreGammaConversionModelMutex);
440  // G4cout << "G4LivermoreGammaConversionModel::InitialiseForElement Z= "
441  // << Z << G4endl;
442  if(!data[Z]) { ReadData(Z); }
443  if(GetTripletModel() && !probTriplet[Z]) { InitialiseProbability(part, Z); }
444  l.unlock();
445 }
446 
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
448 
451 {
452  if(!probTriplet[Z]) {
454  : nullptr;
455  if(0 == nbinsTriplet) {
460  nbinsTriplet = std::max(3,
461  (G4int)(nbins*G4Log(tripletHighEnergy/tripletLowEnergy)/(6*G4Log(10.))));
462  }
463  /*
464  G4cout << "G4LivermoreGammaConversionModel::InitialiseProbability Z= "
465  << Z << " Nbin= " << nbinsTriplet
466  << " Emin(MeV)= " << tripletLowEnergy
467  << " Emax(MeV)= " << tripletHighEnergy << G4endl;
468  */
469  probTriplet[Z] =
471  probTriplet[Z]->SetSpline(true);
472  G4double zz = (G4double)Z;
473  // loop over bins
474  for(G4int j=0; j<=nbinsTriplet; ++j) {
475  G4double e = (probTriplet[Z])->Energy(j);
476  SetupForMaterial(part, mat, e);
477  G4double cross = ComputeCrossSectionPerAtom(part, e, zz);
478  G4double tcross =
480  tcross = (0.0 < cross) ? tcross/cross : 0.0;
481  (probTriplet[Z])->PutValue(j, tcross);
482  //G4cout << j << ". E= " << e << " prob= " << tcross << G4endl;
483  }
484  }
485 }
486 
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ReadData(size_t Z, const char *path=nullptr)
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:792
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:447
G4double ScreenFunction2(G4double screenVariable)
G4int NumberOfBinsPerDecade() const
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
static constexpr double MeV
Definition: G4SIunits.hh:214
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Double_t zz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:552
G4double Value(G4double theEnergy, size_t &lastidx) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:784
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:455
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0.0, G4double cut=0.0, G4double emax=DBL_MAX)
const XML_Char const XML_Char * data
Definition: expat.h:268
fin
Definition: AddMC.C:2
Float_t Z
void InitialiseProbability(const G4ParticleDefinition *, G4int Z)
double G4double
Definition: G4Types.hh:76
TString part[npart]
Definition: Style.C:32
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
static constexpr double millibarn
Definition: G4SIunits.hh:106
G4ParticleDefinition * GetDefinition() const
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual double flat()=0
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4VEmModel.cc:238
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
G4VEmModel * GetTripletModel()
Definition: G4VEmModel.hh:592
static constexpr double electron_mass_c2
G4double GetZ3() const
static constexpr double MeV
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetSpline(G4bool)
static G4Positron * Positron()
Definition: G4Positron.cc:94
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:357
Float_t mat
static G4Electron * Electron()
Definition: G4Electron.cc:94
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
G4int GetZasInt() const
Definition: G4Element.hh:132
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:462
double epsilon(double density, double temperature)
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4lrint(double ad)
Definition: templates.hh:151
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:454
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4LivermoreGammaConversionModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="LivermoreConversion")
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:199
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4LPhysicsFreeVector * data[100]
const G4Material * GetMaterial() const
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
Definition: G4VEmModel.cc:422
G4double ScreenFunction1(G4double screenVariable)
G4double GetfCoulomb() const
Definition: G4Element.hh:191
G4double GetlogZ3() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
double flat()
Definition: G4AblaRandom.cc:48
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
static G4EmParameters * Instance()
static G4PhysicsLogVector * probTriplet[100]
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::mutex G4Mutex
Definition: G4Threading.hh:84