Geant4  10.02.p01
G4PhotonEvaporation.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 // $Id: G4PhotonEvaporation.cc 95483 2016-02-12 10:21:53Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 class file
31 //
32 // CERN, Geneva, Switzerland
33 //
34 // File name: G4PhotonEvaporation
35 //
36 // Author: Vladimir Ivantchenko
37 //
38 // Creation date: 20 December 2011
39 //
40 //Modifications:
41 //
42 //
43 // -------------------------------------------------------------------
44 //
45 
46 #include "G4PhotonEvaporation.hh"
47 
48 #include "Randomize.hh"
49 #include "G4Gamma.hh"
50 #include "G4LorentzVector.hh"
51 #include "G4FragmentVector.hh"
52 #include "G4GammaTransition.hh"
54 #include "G4Pow.hh"
55 #include <CLHEP/Units/SystemOfUnits.h>
56 #include <CLHEP/Units/PhysicalConstants.h>
57 
60 const G4float GREfactor = 5.0f;
61 const G4float GRWfactor = 0.3f;
64 const G4double NormC = 2.5*CLHEP::millibarn/(CLHEP::pi2*CLHEP::hbarc*CLHEP::hbarc);
66 
68  : fLevelManager(nullptr), fTransition(p), fVerbose(0), fPoints(0),
69  vShellNumber(-1), fIndex(0), fTimeLimit(DBL_MAX), fMaxLifeTime(DBL_MAX),
70  fICM(false), fRDM(false), fSampleTime(true)
71 {
73  if(!fTransition) {
74  char* en = getenv("G4UseNuclearPolarization");
75  if(en) { fTransition = new G4PolarizedGammaTransition(); }
76  else { fTransition = new G4GammaTransition(); }
77  }
78 
79  char* env = getenv("G4AddTimeLimitToPhotonEvaporation");
80  if(env) { fTimeLimit = 1.e-16*CLHEP::second; }
81 
82  theA = theZ = fCode = 0;
84 
85  for(G4int i=0; i<MAXDEPOINT; ++i) { fCummProbability[i] = 0.0; }
86  if(0.0f == GREnergy[1]) {
87  G4Pow* g4pow = G4Pow::GetInstance();
88  for (G4int A=1; A<MAXGRDATA; ++A) {
89  GREnergy[A] = (G4float)(40.3*CLHEP::MeV/g4pow->powZ(A,0.2));
91  }
92  }
93 }
94 
96 {
97  delete fTransition;
98 }
99 
100 G4Fragment*
102 {
103  if(fRDM) { fSampleTime = false; }
104  else { fSampleTime = true; }
105 
106  G4Fragment* gamma = GenerateGamma(nucleus);
107  if(fVerbose > 0) {
108  G4cout << "G4PhotonEvaporation::EmittedFragment: RDM= " << fRDM << G4endl;
109  if(gamma) { G4cout << *gamma << G4endl; }
110  G4cout << " Residual: " << *nucleus << G4endl;
111  }
112  return gamma;
113 }
114 
117 {
118  //G4cout << "G4PhotonEvaporation::BreakUpFragment" << G4endl;
119  G4FragmentVector* products = new G4FragmentVector();
120  BreakUpChain(products, nucleus);
121  return products;
122 }
123 
126 {
127  //G4cout << "G4PhotonEvaporation::BreakUp" << G4endl;
128  G4Fragment* aNucleus = new G4Fragment(nucleus);
129  G4FragmentVector* products = new G4FragmentVector();
130  BreakUpChain(products, aNucleus);
131  products->push_back(aNucleus);
132  return products;
133 }
134 
137 {
138  //G4cout << "G4PhotonEvaporation::BreakItUp" << G4endl;
139  G4Fragment* aNucleus = new G4Fragment(nucleus);
140  G4FragmentVector* products = new G4FragmentVector();
141  BreakUpChain(products, aNucleus);
142  products->push_back(aNucleus);
143  return products;
144 }
145 
147  G4Fragment* nucleus)
148 {
149  if(fVerbose > 0) {
150  G4cout << "G4PhotonEvaporation::BreakUpChain RDM= " << fRDM << " "
151  << *nucleus << G4endl;
152  }
153  G4Fragment* gamma = 0;
154  if(fRDM) { fSampleTime = false; }
155  else { fSampleTime = true; }
156 
157  do {
158  gamma = GenerateGamma(nucleus);
159  if(gamma) {
160  products->push_back(gamma);
161  if(fVerbose > 0) {
162  G4cout << "G4PhotonEvaporation::BreakUpChain: "
163  << *gamma << G4endl;
164  G4cout << " Residual: " << *nucleus << G4endl;
165  }
166  // for next decays in the chain always sample time
167  fSampleTime = true;
168  }
169  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
170  } while(gamma);
171  return false;
172 }
173 
174 G4double
176 {
177  fProbability = 0.0;
178  fExcEnergy = nucleus->GetExcitationEnergy();
179  G4int Z = nucleus->GetZ_asInt();
180  G4int A = nucleus->GetA_asInt();
181  fCode = 1000*Z + A;
182  if(fVerbose > 1) {
183  G4cout << "G4PhotonEvaporation::GetEmissionProbability: Z="
184  << Z << " A=" << A << " Eexc(MeV)= " << fExcEnergy << G4endl;
185  }
186  // ignore gamma de-excitation for exotic fragments
187  // and for very low excitations
188  if(0 >= Z || 1 >= A || Z == A || Tolerance >= fExcEnergy)
189  { return fProbability; }
190 
191  // ignore gamma de-excitation for highly excited levels
192  if(A >= MAXGRDATA) { A = MAXGRDATA-1; }
193  //G4cout<<" GREnergy= "<< GREnergy[A]<<" GRWidth= "<<GRWidth[A]<<G4endl;
194  if(fExcEnergy >= (G4double)(GREfactor*GRWidth[A] + GREnergy[A])) {
195  return fProbability;
196  }
197 
198  // probability computed assuming continium transitions
199  // VI: continium transition are limited only to final states
200  // below Fermi energy (this approach needs further evaluation)
201  fFermiEnergy = std::max(0.0, nucleus->ComputeGroundStateMass(Z, A-1)
202  + CLHEP::neutron_mass_c2 - nucleus->GetGroundStateMass());
203 
204  // max energy level for continues transition
206  const G4double eexcfac = 0.99;
207  if(0.0 == emax || fExcEnergy*eexcfac <= emax) { emax = fExcEnergy*eexcfac; }
208 
209  fStep = emax;
211  fStep /= ((G4double)(fPoints - 1));
212  if(fVerbose > 1) {
213  G4cout << "Emax= " << emax << " Npoints= " << fPoints
214  <<" Efermi= " << fFermiEnergy << " Eex= " << fExcEnergy << G4endl;
215  }
216  // integrate probabilities
217  G4double eres = (G4double)GREnergy[A];
218  G4double wres = (G4double)GRWidth[A];
219  G4double eres2= eres*eres;
220  G4double wres2= wres*wres;
221  G4double xsqr = std::sqrt(A*LevelDensity*fExcEnergy);
222 
223  G4double egam = fExcEnergy - emax;
224  G4double gammaE2 = egam*egam;
225  G4double gammaR2 = gammaE2*wres*wres;
226  G4double egdp2 = gammaE2 - eres*eres;
227 
228  G4double p0 = G4Exp(2*(std::sqrt(A*LevelDensity*emax) - xsqr))
229  *gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
230  G4double p1(0.0);
231 
232  for(G4int i=1; i<fPoints; ++i) {
233  egam += fStep;
234  gammaE2 = egam*egam;
235  gammaR2 = gammaE2*wres2;
236  egdp2 = gammaE2 - eres2;
237  //G4cout << "Egamma= " << egam << " Eex= " << fExcEnergy
238  //<< " p0= " << p0 << " p1= " << p1 << G4endl;
239  p1 = G4Exp(2*(std::sqrt(A*LevelDensity*std::abs(fExcEnergy - egam)) - xsqr))
240  *gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
241  fProbability += (p1 + p0);
243  p0 = p1;
244  }
245  fProbability *= 0.5*fStep*NormC*A;
246  if(fVerbose > 1) { G4cout << "prob= " << fProbability << G4endl; }
247  return fProbability;
248 }
249 
250 G4double
252 {
253  G4double E = energy;
255  if(fLevelManager) {
257  if(E > fLevelEnergyMax + Tolerance) { E = energy; }
258  }
259  return E;
260 }
261 
263 {
265  return fLevelEnergyMax;
266 }
267 
268 G4Fragment*
270 {
271  G4Fragment* result = 0;
272  G4double eexc = nucleus->GetExcitationEnergy();
273  if(eexc < Tolerance) { return result; }
274 
275  InitialiseLevelManager(nucleus->GetZ_asInt(), nucleus->GetA_asInt());
276 
277  G4double time = nucleus->GetCreationTime();
278 
279  G4double efinal = 0.0;
280  vShellNumber = -1;
281  size_t shell = 0;
282  G4int deltaS = 1;
283  G4bool isGamma = true;
284  G4bool isLongLived = false;
285  G4bool isX = false;
286  G4bool icm = fICM;
287 
288  if(fVerbose > 1) {
289  G4cout << "GenerateGamma: Exc= " << eexc << " Emax= "
290  << fLevelEnergyMax << " fEex= " << fExcEnergy
291  << " fCode= " << fCode << " fPoints= " << fPoints
292  << " fProb= " << fProbability << G4endl;
293  }
294  // continues part
295  if(!fLevelManager || eexc > fLevelEnergyMax + Tolerance) {
296  //G4cout << "Continues fPoints= " << fPoints << " " << fLevelManager << G4endl;
297 
298  // we compare current excitation versus value used for probability
299  // computation and also Z and A used for probability computation
300  if(fCode != 1000*theZ + theA || eexc != fExcEnergy) {
301  GetEmissionProbability(nucleus);
302  }
303  if(fProbability == 0.0) { return result; }
305  for(G4int i=1; i<fPoints; ++i) {
306  //G4cout << "y= " << y << " cummProb= " << fCummProbability[i] << G4endl;
307  if(y <= fCummProbability[i]) {
308  efinal = fStep*((i - 1) + (y - fCummProbability[i-1])
309  /(fCummProbability[i] - fCummProbability[i-1]));
310  break;
311  }
312  }
313  // final discrete level
314  if(fLevelManager && efinal <= fLevelEnergyMax + Tolerance) {
315  //G4cout << "Efinal= " << efinal << " idx= " << fIndex << G4endl;
318  }
319  //discrete part
320  } else {
322  if(fVerbose > 1) {
323  G4cout << "Discrete emission from level Index= " << fIndex
324  << " Elevel= " << fLevelManager->LevelEnergy(fIndex)
325  << " RDM= " << fRDM << " ICM= " << fICM << G4endl;
326  }
327  if(0 == fIndex) { return result; }
328  const G4NucLevel* level = fLevelManager->GetLevel(fIndex);
329 
330  // radioactive decay model call transitions one by one
331  // for X-level corresponding gamma transition is forced
332  // from the previous one; this is needed, because
333  // X-level assumes a competition between beta amd gamma
334  // decay probabilities
335  isX = level->IsXLevel();
336  if(isX && fRDM && fIndex > 0) {
337  --fIndex;
338  level = fLevelManager->GetLevel(fIndex);
339  isX = level->IsXLevel();
340  }
341 
342  G4double ltime = 0.0;
343  if(fICM) { ltime = (G4double)fLevelManager->LifeTime(fIndex); }
344  else { ltime = (G4double)fLevelManager->LifeTimeGamma(fIndex); }
345 
346  if(ltime >= fMaxLifeTime) { return result; }
347  if(ltime > fTimeLimit) {
348  icm = true;
349  isLongLived = true;
350  }
351  size_t ntrans = level->NumberOfTransitions();
352  size_t idx = 0;
353  if(fVerbose > 1) {
354  G4cout << "Ntrans= " << ntrans << " idx= " << idx << " isX= " << isX
355  << " icm= " << icm << G4endl;
356  }
357  if(!isX) {
358  if(1 < ntrans) {
359  G4double rndm = G4UniformRand();
360  if(icm) { idx = level->SampleGammaETransition(rndm); }
361  else { idx = level->SampleGammaTransition(rndm); }
362  //G4cout << "Sampled idx= " << idx << " rndm= " << rndm << G4endl;
363  }
364  if(icm) {
365  G4double rndm = G4UniformRand();
366  G4double prob = level->GammaProbability(idx);
367  if(rndm > prob) {
368  rndm = (rndm - prob)/(1.0 - prob);
369  shell = level->SampleShell(idx, rndm);
370  vShellNumber = shell;
371  isGamma = false;
372  }
373  }
374  }
375  efinal = (G4double)level->FinalExcitationEnergy(idx);
376  fIndex = idx;
377 
378  if(fSampleTime && ltime > 0.0) {
379  time -= ltime*G4Log(G4UniformRand());
380  }
381  }
382  // no gamma emission for X-level
383  if(isX) {
384  G4ThreeVector v = nucleus->GetMomentum().vect().unit();
385  G4double mass = nucleus->GetGroundStateMass() + efinal;
386  G4double e = std::max(mass,nucleus->GetMomentum().e());
387  G4double mom = std::sqrt((e - mass)*(e + mass));
388  v *= mom;
389  nucleus->SetMomentum(G4LorentzVector(v.x(),v.y(),v.z(),efinal+mass));
390  if(fVerbose > 1) {
391  G4cout << "X-level Eexc= " << nucleus->GetExcitationEnergy()
392  << G4endl;
393  }
394  // normal case
395  } else {
396  result = fTransition->SampleTransition(nucleus, efinal,
397  deltaS, shell,
398  isGamma, isLongLived);
399  if(result) { result->SetCreationTime(time); }
400  }
401  nucleus->SetCreationTime(time);
402 
403  if(fVerbose > 1) {
404  G4cout << "Final level E= " << efinal << " time= " << time
405  << " idx= " << fIndex << " isX " << isX
406  << " isGamma: " << isGamma << " isLongLived: " << isLongLived
407  << " deltaS= " << deltaS << " shell= " << shell << G4endl;
408  }
409  return result;
410 }
411 
413 {
414  static const G4double tfact = G4Pow::GetInstance()->logZ(2);
415  fMaxLifeTime = val/tfact;
416 }
417 
419 {
420  if(p != fTransition) {
421  delete fTransition;
422  fTransition = p;
423  }
424 }
425 
427 {
428  fICM = val;
429 }
430 
432 {
433  fRDM = val;
434 }
435 
virtual void SetICM(G4bool)
virtual G4bool BreakUpChain(G4FragmentVector *theResult, G4Fragment *theNucleus)
G4float GammaProbability(size_t idx) const
Definition: G4NucLevel.hh:115
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
void SetMaxHalfLife(G4double)
virtual G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
static const double MeV
Definition: G4SIunits.hh:211
G4GammaTransition * fTransition
const G4int MAXGRDATA
static G4float GREnergy[MAXGRDATA]
virtual G4double GetUpperLevelEnergy(G4int Z, G4int A)
CLHEP::Hep3Vector G4ThreeVector
const G4int MAXDEPOINT
const G4double Tolerance
Definition: G4Pow.hh:56
virtual void RDMForced(G4bool)
G4float LifeTime(size_t i) const
float G4float
Definition: G4Types.hh:77
virtual G4double GetFinalLevelEnergy(G4int Z, G4int A, G4double energy)
static G4float GRWidth[MAXGRDATA]
G4double fCummProbability[MAXDEPOINT]
G4PhotonEvaporation(G4GammaTransition *ptr=0)
const G4LevelManager * fLevelManager
static const double pi2
Definition: G4SIunits.hh:77
int G4int
Definition: G4Types.hh:78
const G4NucLevel * GetLevel(size_t i) const
virtual G4FragmentVector * BreakUp(const G4Fragment &theNucleus)
G4float FinalExcitationEnergy(size_t idx) const
Definition: G4NucLevel.hh:109
G4double logZ(G4int Z) const
Definition: G4Pow.hh:166
G4float LevelEnergy(size_t i) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:251
G4double GetCreationTime() const
Definition: G4Fragment.hh:405
G4bool IsXLevel() const
Definition: G4NucLevel.hh:104
bool G4bool
Definition: G4Types.hh:79
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:284
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:289
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
static const double second
Definition: G4SIunits.hh:156
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
void SetGammaTransition(G4GammaTransition *)
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:273
size_t SampleGammaTransition(G4double rndm) const
Definition: G4NucLevel.hh:127
G4float NearestLevelEnergy(G4double energy, size_t index=0) const
size_t NearestLevelIndex(G4double energy, size_t index=0) const
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4float GRWfactor
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:410
size_t NumberOfTransitions() const
Definition: G4NucLevel.hh:99
virtual G4Fragment * SampleTransition(G4Fragment *nucleus, G4double newExcEnergy, G4int deltaS, size_t shell, G4bool isGamma, G4bool isLongLived)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4Fragment * GenerateGamma(G4Fragment *nucleus)
const G4double MaxDeltaEnergy
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4NuclearLevelData * fNuclearLevelData
size_t SampleShell(size_t idx, G4double rndm) const
Definition: G4NucLevel.hh:147
static const double millibarn
Definition: G4SIunits.hh:105
G4int GetZ_asInt() const
Definition: G4Fragment.hh:256
const G4float GREfactor
#define G4endl
Definition: G4ios.hh:61
const G4double LevelDensity
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
size_t SampleGammaETransition(G4double rndm) const
Definition: G4NucLevel.hh:137
static const double keV
Definition: G4SIunits.hh:213
virtual G4FragmentVector * BreakUpFragment(G4Fragment *theNucleus)
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:254
double G4double
Definition: G4Types.hh:76
G4float LifeTimeGamma(size_t i) const
void InitialiseLevelManager(G4int Z, G4int A)
#define DBL_MAX
Definition: templates.hh:83
G4double ComputeGroundStateMass(G4int Z, G4int A) const
Definition: G4Fragment.hh:241
static G4NuclearLevelData * GetInstance()
const G4double NormC
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:268
CLHEP::HepLorentzVector G4LorentzVector