Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 105162 2017-07-14 10:22: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"
53 #include "G4Pow.hh"
56 
57 G4float G4PhotonEvaporation::GREnergy[] = {0.0f};
58 G4float G4PhotonEvaporation::GRWidth[] = {0.0f};
59 
61  : fLevelManager(nullptr), fTransition(p), fVerbose(0), fPoints(0),
62  vShellNumber(-1), fIndex(0), fMaxLifeTime(DBL_MAX),
63  fICM(false), fRDM(false), fSampleTime(true), isInitialised(false)
64 {
65  //G4cout << "### New G4PhotonEvaporation() " << this << G4endl;
66  fNuclearLevelData = G4NuclearLevelData::GetInstance();
67  LevelDensity = 0.125/CLHEP::MeV;
68  Tolerance = 20*CLHEP::eV;
69 
70  if(!fTransition) { fTransition = new G4GammaTransition(); }
71 
72  theA = theZ = fCode = 0;
73  fLevelEnergyMax = fStep = fExcEnergy = fProbability = 0.0;
74 
75  for(G4int i=0; i<MAXDEPOINT; ++i) { fCummProbability[i] = 0.0; }
76  if(0.0f == GREnergy[1]) {
77  G4Pow* g4calc = G4Pow::GetInstance();
78  static const G4float GRWfactor = 0.3f;
79  for (G4int A=1; A<MAXGRDATA; ++A) {
80  GREnergy[A] = (G4float)(40.3*CLHEP::MeV/g4calc->powZ(A,0.2));
81  GRWidth[A] = GRWfactor*GREnergy[A];
82  }
83  }
84 }
85 
87 {
88  delete fTransition;
89 }
90 
92 {
93  if(isInitialised) { return; }
94  isInitialised = true;
95 
96  if(fVerbose > 0) {
97  G4cout << "### G4PhotonEvaporation is initialized " << this << G4endl;
98  }
99  G4DeexPrecoParameters* param = fNuclearLevelData->GetParameters();
100  LevelDensity = param->GetLevelDensity();
101  Tolerance = param->GetMinExcitation();
102  fMaxLifeTime = param->GetMaxLifeTime();
103 
104  fTransition->SetPolarizationFlag(param->CorrelatedGamma());
105  fTransition->SetVerbose(fVerbose);
106 }
107 
108 G4Fragment*
110 {
111  if(fRDM) { fSampleTime = false; }
112  else { fSampleTime = true; }
113 
114  G4Fragment* gamma = GenerateGamma(nucleus);
115  if(fVerbose > 0) {
116  G4cout << "G4PhotonEvaporation::EmittedFragment: RDM= " << fRDM << G4endl;
117  if(gamma) { G4cout << *gamma << G4endl; }
118  G4cout << " Residual: " << *nucleus << G4endl;
119  }
120  return gamma;
121 }
122 
125 {
126  if(fVerbose > 0) {
127  G4cout << "G4PhotonEvaporation::BreakItUp" << G4endl;
128  }
129  G4Fragment* aNucleus = new G4Fragment(nucleus);
130  G4FragmentVector* products = new G4FragmentVector();
131  BreakUpChain(products, aNucleus);
132  products->push_back(aNucleus);
133  return products;
134 }
135 
137  G4Fragment* nucleus)
138 {
139  if(!isInitialised) { Initialise(); }
140  if(fVerbose > 0) {
141  G4cout << "G4PhotonEvaporation::BreakUpChain RDM= " << fRDM << " "
142  << *nucleus << G4endl;
143  }
144  G4Fragment* gamma = nullptr;
145  if(fRDM) { fSampleTime = false; }
146  else { fSampleTime = true; }
147 
148  do {
149  gamma = GenerateGamma(nucleus);
150  if(gamma) {
151  products->push_back(gamma);
152  if(fVerbose > 0) {
153  G4cout << "G4PhotonEvaporation::BreakUpChain: "
154  << *gamma << G4endl;
155  G4cout << " Residual: " << *nucleus << G4endl;
156  }
157  // for next decays in the chain always sample time
158  fSampleTime = true;
159  }
160  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
161  } while(gamma);
162  return false;
163 }
164 
165 G4double
167 {
168  if(!isInitialised) { Initialise(); }
169  fProbability = 0.0;
170  fExcEnergy = nucleus->GetExcitationEnergy();
171  G4int Z = nucleus->GetZ_asInt();
172  G4int A = nucleus->GetA_asInt();
173  fCode = 1000*Z + A;
174  if(fVerbose > 1) {
175  G4cout << "G4PhotonEvaporation::GetEmissionProbability: Z="
176  << Z << " A=" << A << " Eexc(MeV)= " << fExcEnergy << G4endl;
177  }
178  // ignore gamma de-excitation for exotic fragments
179  // and for very low excitations
180  if(0 >= Z || 1 >= A || Z == A || Tolerance >= fExcEnergy)
181  { return fProbability; }
182 
183  // ignore gamma de-excitation for highly excited levels
184  if(A >= MAXGRDATA) { A = MAXGRDATA-1; }
185  //G4cout<<" GREnergy= "<< GREnergy[A]<<" GRWidth= "<<GRWidth[A]<<G4endl;
186 
187  static const G4float GREfactor = 5.0f;
188  if(fExcEnergy >= (G4double)(GREfactor*GRWidth[A] + GREnergy[A])) {
189  return fProbability;
190  }
191  // probability computed assuming continium transitions
192  // VI: continium transition are limited only to final states
193  // below Fermi energy (this approach needs further evaluation)
194  G4double emax = std::max(0.0, nucleus->ComputeGroundStateMass(Z, A-1)
196 
197  // max energy level for continues transition
198  emax = std::min(emax, fExcEnergy);
199  const G4double eexcfac = 0.99;
200  if(0.0 == emax || fExcEnergy*eexcfac <= emax) { emax = fExcEnergy*eexcfac; }
201 
202  fStep = emax;
203  static const G4double MaxDeltaEnergy = CLHEP::MeV;
204  fPoints = std::min((G4int)(fStep/MaxDeltaEnergy) + 2, MAXDEPOINT);
205  fStep /= ((G4double)(fPoints - 1));
206  if(fVerbose > 1) {
207  G4cout << "Emax= " << emax << " Npoints= " << fPoints
208  << " Eex= " << fExcEnergy << G4endl;
209  }
210  // integrate probabilities
211  G4double eres = (G4double)GREnergy[A];
212  G4double wres = (G4double)GRWidth[A];
213  G4double eres2= eres*eres;
214  G4double wres2= wres*wres;
215  G4double xsqr = std::sqrt(A*LevelDensity*fExcEnergy);
216 
217  G4double egam = fExcEnergy;
218  G4double gammaE2 = egam*egam;
219  G4double gammaR2 = gammaE2*wres2;
220  G4double egdp2 = gammaE2 - eres2;
221 
222  G4double p0 = G4Exp(-2.0*xsqr)*gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
223  G4double p1(0.0);
224 
225  for(G4int i=1; i<fPoints; ++i) {
226  egam -= fStep;
227  gammaE2 = egam*egam;
228  gammaR2 = gammaE2*wres2;
229  egdp2 = gammaE2 - eres2;
230  p1 = G4Exp(2.0*(std::sqrt(A*LevelDensity*std::abs(fExcEnergy - egam)) - xsqr))
231  *gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
232  fProbability += (p1 + p0);
233  fCummProbability[i] = fProbability;
234  //G4cout << "Egamma= " << egam << " Eex= " << fExcEnergy
235  //<< " p0= " << p0 << " p1= " << p1 << " sum= " << fCummProbability[i] <<G4endl;
236  p0 = p1;
237  }
238 
239  static const G4double NormC = 1.25*CLHEP::millibarn
241  fProbability *= fStep*NormC*A;
242  if(fVerbose > 1) { G4cout << "prob= " << fProbability << G4endl; }
243  return fProbability;
244 }
245 
246 G4double
248 {
249  G4double E = energy;
250  InitialiseLevelManager(Z, A);
251  if(fLevelManager) {
252  E = (G4double)fLevelManager->NearestLevelEnergy(energy, fIndex);
253  if(E > fLevelEnergyMax + Tolerance) { E = energy; }
254  }
255  return E;
256 }
257 
259 {
260  InitialiseLevelManager(Z, A);
261  return fLevelEnergyMax;
262 }
263 
264 G4Fragment*
265 G4PhotonEvaporation::GenerateGamma(G4Fragment* nucleus)
266 {
267  if(!isInitialised) { Initialise(); }
268  G4Fragment* result = nullptr;
269  G4double eexc = nucleus->GetExcitationEnergy();
270  if(eexc <= Tolerance) { return result; }
271 
272  InitialiseLevelManager(nucleus->GetZ_asInt(), nucleus->GetA_asInt());
273 
274  G4double time = nucleus->GetCreationTime();
275 
276  G4double efinal = 0.0;
277  G4double ratio = 0.0;
278  vShellNumber = -1;
279  G4int JP1 = 0;
280  G4int JP2 = 0;
281  G4int multiP = 0;
282  G4bool isGamma = true;
283  G4bool isDiscrete = false;
284 
285  const G4NucLevel* level = nullptr;
286  size_t ntrans = 0;
287 
288  if(fLevelManager && eexc <= fLevelEnergyMax + Tolerance) {
289  fIndex = fLevelManager->NearestLevelIndex(eexc, fIndex);
290  if(0 < fIndex) {
291  // for discrete transition
292  level = fLevelManager->GetLevel(fIndex);
293  if(level) {
294  ntrans = level->NumberOfTransitions();
295  JP1 = fLevelManager->SpinTwo(fIndex);
296  if(ntrans > 0) { isDiscrete = true; }
297  else if(fLevelManager->FloatingLevel(fIndex) > 0) {
298  --fIndex;
299  level = fLevelManager->GetLevel(fIndex);
300  ntrans = level->NumberOfTransitions();
301  JP1 = fLevelManager->SpinTwo(fIndex);
302  if(ntrans > 0) { isDiscrete = true; }
303  }
304  }
305  }
306  }
307  if(fVerbose > 1) {
308  G4int prec = G4cout.precision(4);
309  G4cout << "GenerateGamma: Z= " << nucleus->GetZ_asInt()
310  << " A= " << nucleus->GetA_asInt()
311  << " Exc= " << eexc << " Emax= "
312  << fLevelEnergyMax << " idx= " << fIndex
313  << " fCode= " << fCode << " fPoints= " << fPoints
314  << " Ntr= " << ntrans << " discrete: " << isDiscrete
315  << " fProb= " << fProbability << G4endl;
316  G4cout.precision(prec);
317  }
318 
319  // continues part
320  if(!isDiscrete) {
321  // G4cout << "Continues fIndex= " << fIndex << G4endl;
322 
323  // we compare current excitation versus value used for probability
324  // computation and also Z and A used for probability computation
325  if(fCode != 1000*theZ + theA || eexc != fExcEnergy) {
326  GetEmissionProbability(nucleus);
327  }
328  if(fProbability == 0.0) { return result; }
329  G4double y = fCummProbability[fPoints-1]*G4UniformRand();
330  for(G4int i=1; i<fPoints; ++i) {
331  //G4cout << "y= " << y << " cummProb= " << fCummProbability[i] << G4endl;
332  if(y <= fCummProbability[i]) {
333  efinal = fStep*((i - 1) + (y - fCummProbability[i-1])
334  /(fCummProbability[i] - fCummProbability[i-1]));
335  break;
336  }
337  }
338  // final discrete level
339  if(fLevelManager) {
340  if(efinal < fLevelEnergyMax) {
341  //G4cout << "Efinal= " << efinal << " idx= " << fIndex << G4endl;
342  fIndex = fLevelManager->NearestLevelIndex(efinal, fIndex);
343  efinal = fLevelManager->LevelEnergy(fIndex);
344  // protection - take level below
345  if(efinal >= eexc && 0 < fIndex) {
346  --fIndex;
347  efinal = fLevelManager->LevelEnergy(fIndex);
348  }
349  nucleus->SetFloatingLevelNumber(fLevelManager->FloatingLevel(fIndex));
350 
351  // not allowed to have final energy above max energy
352  // if G4LevelManager exist
353  } else {
354  efinal = fLevelEnergyMax;
355  fIndex = fLevelManager->NearestLevelIndex(efinal, fIndex);
356  }
357  }
358  if(fVerbose > 1) {
359  G4cout << "Continues emission efinal(MeV)= " << efinal << G4endl;
360  }
361  //discrete part
362  } else {
363  if(fVerbose > 1) {
364  G4cout << "Discrete emission from level Index= " << fIndex
365  << " Elevel= " << fLevelManager->LevelEnergy(fIndex)
366  << " RDM= " << fRDM << " ICM= " << fICM << G4endl;
367  }
368  if(0 == fIndex || !level) { return result; }
369 
370  // stable fragment has life time -1
371  // if called from radioactive decay the life time is not checked
372  G4double ltime = fLevelManager->LifeTime(fIndex);
373  if(ltime < 0.0 || (!fRDM && ltime > fMaxLifeTime)) { return result; }
374 
375  size_t idx = 0;
376  if(1 < ntrans) {
377  idx = level->SampleGammaTransition(G4UniformRand());
378  }
379  if(fVerbose > 1) {
380  G4cout << "Ntrans= " << ntrans << " idx= " << idx
381  << " ICM= " << fICM << G4endl;
382  }
383  G4double prob = (G4double)level->GammaProbability(idx);
384  // prob = 0 means that there is only internal conversion
385  // if((fICM && prob < 1.0) || prob == 0.0) {
386  if(prob < 1.0) {
387  G4double rndm = G4UniformRand();
388  if(rndm > prob) {
389  isGamma = false;
390  rndm = (rndm - prob)/(1.0 - prob);
391  vShellNumber = level->SampleShell(idx, rndm);
392  }
393  }
394  // it is discrete transition with possible gamma correlation
395  isDiscrete = true;
396  ratio = level->MixingRatio(idx);
397  multiP = level->TransitionType(idx);
398  fIndex = level->FinalExcitationIndex(idx);
399  JP2 = fLevelManager->SpinTwo(fIndex);
400 
401  // final energy and time
402  efinal = (G4double)fLevelManager->LevelEnergy(fIndex);
403  if(fSampleTime && ltime > 0.0) {
404  time -= ltime*G4Log(G4UniformRand());
405  }
406  nucleus->SetFloatingLevelNumber(fLevelManager->FloatingLevel(fIndex));
407  }
408  // protection for floating levels
409  if(std::abs(efinal - eexc) <= Tolerance) { return result; }
410 
411  result = fTransition->SampleTransition(nucleus, efinal, ratio, JP1,
412  JP2, multiP, vShellNumber,
413  isDiscrete, isGamma);
414  if(result) { result->SetCreationTime(time); }
415  nucleus->SetCreationTime(time);
416 
417  // ignore the floating levels with zero energy and create ground state
418  if(efinal == 0.0 && fIndex > 0) {
419  fIndex = 0;
420  nucleus->SetFloatingLevelNumber(fLevelManager->FloatingLevel(fIndex));
421  }
422 
423  if(fVerbose > 1) {
424  G4cout << "Final level E= " << efinal << " time= " << time
425  << " idxFinal= " << fIndex << " isDiscrete: " << isDiscrete
426  << " isGamma: " << isGamma << " multiP= " << multiP
427  << " shell= " << vShellNumber << G4endl;
428  }
429  return result;
430 }
431 
433 {}
434 
436 {
437  if(p != fTransition) {
438  delete fTransition;
439  fTransition = p;
440  }
441 }
442 
444 {
445  fICM = val;
446 }
447 
449 {
450  fRDM = val;
451 }
452 
G4double G4ParticleHPJENDLHEData::G4double result
virtual void SetICM(G4bool)
G4float GammaProbability(size_t idx) const
Definition: G4NucLevel.hh:142
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
void SetMaxHalfLife(G4double)
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
G4int SpinTwo(size_t i) const
const G4int MAXGRDATA
virtual G4bool BreakUpChain(G4FragmentVector *theResult, G4Fragment *theNucleus) final
const G4int MAXDEPOINT
Definition: G4Pow.hh:56
virtual void RDMForced(G4bool)
static constexpr double pi2
Definition: SystemOfUnits.h:57
static constexpr double millibarn
Definition: SystemOfUnits.h:86
static constexpr double hbarc
G4PhotonEvaporation(G4GammaTransition *ptr=nullptr)
G4float LifeTime(size_t i) const
const char * p
Definition: xmltok.h:285
float G4float
Definition: G4Types.hh:77
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus) final
virtual void Initialise() final
G4double GetLevelDensity() const
int G4int
Definition: G4Types.hh:78
const G4NucLevel * GetLevel(size_t i) const
G4float MixingRatio(size_t idx) const
Definition: G4NucLevel.hh:134
G4double GetMaxLifeTime() const
void SetPolarizationFlag(G4bool val)
static const double prec
Definition: RanecuEngine.cc:58
virtual G4double GetEmissionProbability(G4Fragment *theNucleus) final
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:266
G4double GetCreationTime() const
Definition: G4Fragment.hh:448
bool G4bool
Definition: G4Types.hh:79
static constexpr double neutron_mass_c2
static constexpr double MeV
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4DeexPrecoParameters * GetParameters()
void SetGammaTransition(G4GammaTransition *)
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:288
size_t SampleGammaTransition(G4double rndm) const
Definition: G4NucLevel.hh:166
G4float NearestLevelEnergy(G4double energy, size_t index=0) const
size_t NearestLevelIndex(G4double energy, size_t index=0) const
static constexpr double eV
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
virtual G4Fragment * SampleTransition(G4Fragment *nucleus, G4double newExcEnergy, G4double mpRatio, G4int JP1, G4int JP2, G4int MP, G4int shell, G4bool isDiscrete, G4bool isGamma)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:453
size_t NumberOfTransitions() const
Definition: G4NucLevel.hh:108
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4int FloatingLevel(size_t i) const
G4int TransitionType(size_t idx) const
Definition: G4NucLevel.hh:121
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
void SetFloatingLevelNumber(G4int value)
Definition: G4Fragment.hh:432
G4int SampleShell(size_t idx, G4double rndm) const
Definition: G4NucLevel.hh:176
#define G4endl
Definition: G4ios.hh:61
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:254
double G4double
Definition: G4Types.hh:76
size_t FinalExcitationIndex(size_t idx) const
Definition: G4NucLevel.hh:113
#define DBL_MAX
Definition: templates.hh:83
virtual G4double GetUpperLevelEnergy(G4int Z, G4int A) final
G4double ComputeGroundStateMass(G4int Z, G4int A) const
Definition: G4Fragment.hh:256
static G4NuclearLevelData * GetInstance()
G4double GetMinExcitation() const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283
virtual G4double GetFinalLevelEnergy(G4int Z, G4int A, G4double energy) final
void SetVerbose(G4int val)