Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 = 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 = 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......
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
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double MeV
Definition: G4SIunits.hh:214
G4ParticleChangeForGamma * fParticleChange
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
static const G4double * P2[nN]
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double S(double temp)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
const XML_Char const XML_Char * data
Definition: expat.h:268
fin
Definition: AddMC.C:2
Float_t 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
virtual double flat()=0
static constexpr double electron_mass_c2
static constexpr double twopi
Definition: G4SIunits.hh:76
static G4Positron * Positron()
Definition: G4Positron.cc:94
FILE * fp
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4LPhysicsFreeVector * data[100]
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
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4Material * GetMaterial() const
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
G4BoldyshevTripletModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="BoldyshevTripletConversion")
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void ProposeTrackStatus(G4TrackStatus status)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
static constexpr double GeV
Definition: G4SIunits.hh:217
double flat()
Definition: G4AblaRandom.cc:48
static constexpr double pi
Definition: SystemOfUnits.h:54
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
void ReadData(size_t Z, const char *path=nullptr)