Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PhotonEvaporation.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 // $Id: G4PhotonEvaporation.cc 106723 2017-10-20 09:50:34Z 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 
49 #include "Randomize.hh"
50 #include "G4Gamma.hh"
51 #include "G4LorentzVector.hh"
52 #include "G4FragmentVector.hh"
53 #include "G4GammaTransition.hh"
54 #include "G4Pow.hh"
57 
60 
61 #ifdef G4MULTITHREADED
62 G4Mutex G4PhotonEvaporation::PhotonEvaporationMutex = G4MUTEX_INITIALIZER;
63 #endif
64 
66  : fLevelManager(nullptr), fTransition(p), fPolarization(nullptr),
67  fVerbose(0), fPoints(0), vShellNumber(-1), fIndex(0),
68  fMaxLifeTime(DBL_MAX),
69  fICM(true), fRDM(false), fSampleTime(true),
70  fCorrelatedGamma(false), isInitialised(false)
71 {
72  //G4cout << "### New G4PhotonEvaporation() " << this << G4endl;
75  LevelDensity = 0.125/CLHEP::MeV;
76  Tolerance = 20*CLHEP::eV;
77 
78  if(!fTransition) { fTransition = new G4GammaTransition(); }
79 
80  theA = theZ = fCode = 0;
82 
83  for(G4int i=0; i<MAXDEPOINT; ++i) { fCummProbability[i] = 0.0; }
84  if(0.0f == GREnergy[1]) { InitialiseGRData(); }
85 }
86 
88 {
89  delete fTransition;
90 }
91 
93 {
94  if(isInitialised) { return; }
95  isInitialised = true;
96 
97  if(fVerbose > 0) {
98  G4cout << "### G4PhotonEvaporation is initialized " << this << G4endl;
99  }
101  LevelDensity = param->GetLevelDensity();
102  Tolerance = param->GetMinExcitation();
103  fMaxLifeTime = param->GetMaxLifeTime();
105  fICM = param->GetInternalConversionFlag();
106 
108  fTransition->SetTwoJMAX(param->GetTwoJMAX());
110 }
111 
113 {
114 #ifdef G4MULTITHREADED
115  G4MUTEXLOCK(&G4PhotonEvaporation::PhotonEvaporationMutex);
116 #endif
117  if(0.0f == GREnergy[1]) {
118  G4Pow* g4calc = G4Pow::GetInstance();
119  static const G4float GRWfactor = 0.3f;
120  for (G4int A=1; A<MAXGRDATA; ++A) {
121  GREnergy[A] = (G4float)(40.3*CLHEP::MeV/g4calc->powZ(A,0.2));
122  GRWidth[A] = GRWfactor*GREnergy[A];
123  }
124  }
125 #ifdef G4MULTITHREADED
126  G4MUTEXUNLOCK(&G4PhotonEvaporation::PhotonEvaporationMutex);
127 #endif
128 }
129 
130 G4Fragment*
132 {
133  if(!isInitialised) { Initialise(); }
134  fSampleTime = (fRDM) ? false : true;
135 
136  // potentially external code may set initial polarization
137  // but only for radioactive decay nuclear polarization is considered
138  if(fCorrelatedGamma && fRDM) {
139  if(nucleus->GetNuclearPolarization()) {
141  delete nucleus->GetNuclearPolarization();
142  }
144  nucleus->GetA_asInt(),
145  nucleus->GetExcitationEnergy());
147  }
148  if(fVerbose > 1) {
149  G4cout << "G4PhotonEvaporation::EmittedFragment: "
150  << *nucleus << G4endl;
151  if(fPolarization) { G4cout << "NucPolar: " << fPolarization << G4endl; }
152  G4cout << " CorrGamma: " << fCorrelatedGamma << " RDM: " << fRDM
153  << " fPolarization: " << fPolarization << G4endl;
154  }
155  G4Fragment* gamma = GenerateGamma(nucleus);
156 
157  // remove G4NuclearPolarizaton when reach ground state
158  if(fPolarization && 0 == fIndex) {
159  if(fVerbose > 1) {
160  G4cout << "G4PhotonEvaporation::EmittedFragment: remove "
161  << fPolarization << G4endl;
162  }
164  fPolarization = nullptr;
166  }
167 
168  if(fVerbose > 1) {
169  G4cout << "G4PhotonEvaporation::EmittedFragment: RDM= "
170  << fRDM << " done:" << G4endl;
171  if(gamma) { G4cout << *gamma << G4endl; }
172  G4cout << " Residual: " << *nucleus << G4endl;
173  }
174  return gamma;
175 }
176 
179 {
180  if(fVerbose > 0) {
181  G4cout << "G4PhotonEvaporation::BreakItUp" << G4endl;
182  }
183  G4Fragment* aNucleus = new G4Fragment(nucleus);
184  G4FragmentVector* products = new G4FragmentVector();
185  BreakUpChain(products, aNucleus);
186  products->push_back(aNucleus);
187  return products;
188 }
189 
191  G4Fragment* nucleus)
192 {
193  if(!isInitialised) { Initialise(); }
194  if(fVerbose > 0) {
195  G4cout << "G4PhotonEvaporation::BreakUpChain RDM= " << fRDM << " "
196  << *nucleus << G4endl;
197  }
198  G4Fragment* gamma = nullptr;
199  fSampleTime = (fRDM) ? false : true;
200 
201  // start decay chain from unpolarized state
202  if(fCorrelatedGamma) {
204  nucleus->GetA_asInt(),
205  nucleus->GetExcitationEnergy());
207  }
208 
209  do {
210  gamma = GenerateGamma(nucleus);
211  if(gamma) {
212  products->push_back(gamma);
213  if(fVerbose > 0) {
214  G4cout << "G4PhotonEvaporation::BreakUpChain: "
215  << *gamma << G4endl;
216  G4cout << " Residual: " << *nucleus << G4endl;
217  }
218  // for next decays in the chain always sample time
219  fSampleTime = true;
220  }
221  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
222  } while(gamma);
223 
224  // clear nuclear polarization end of chain
225  if(fPolarization) {
226  delete fPolarization;
227  fPolarization = nullptr;
229  }
230  return false;
231 }
232 
233 G4double
235 {
236  if(!isInitialised) { Initialise(); }
237  fProbability = 0.0;
238  fExcEnergy = nucleus->GetExcitationEnergy();
239  G4int Z = nucleus->GetZ_asInt();
240  G4int A = nucleus->GetA_asInt();
241  fCode = 1000*Z + A;
242  if(fVerbose > 1) {
243  G4cout << "G4PhotonEvaporation::GetEmissionProbability: Z="
244  << Z << " A=" << A << " Eexc(MeV)= " << fExcEnergy << G4endl;
245  }
246  // ignore gamma de-excitation for exotic fragments
247  // and for very low excitations
248  if(0 >= Z || 1 >= A || Z == A || Tolerance >= fExcEnergy)
249  { return fProbability; }
250 
251  // ignore gamma de-excitation for highly excited levels
252  if(A >= MAXGRDATA) { A = MAXGRDATA-1; }
253  //G4cout<<" GREnergy= "<< GREnergy[A]<<" GRWidth= "<<GRWidth[A]<<G4endl;
254 
255  static const G4float GREfactor = 5.0f;
256  if(fExcEnergy >= (G4double)(GREfactor*GRWidth[A] + GREnergy[A])) {
257  return fProbability;
258  }
259  // probability computed assuming continium transitions
260  // VI: continium transition are limited only to final states
261  // below Fermi energy (this approach needs further evaluation)
262  G4double emax = std::max(0.0, nucleus->ComputeGroundStateMass(Z, A-1)
264 
265  // max energy level for continues transition
266  emax = std::min(emax, fExcEnergy);
267  const G4double eexcfac = 0.99;
268  if(0.0 == emax || fExcEnergy*eexcfac <= emax) { emax = fExcEnergy*eexcfac; }
269 
270  fStep = emax;
271  static const G4double MaxDeltaEnergy = CLHEP::MeV;
272  fPoints = std::min((G4int)(fStep/MaxDeltaEnergy) + 2, MAXDEPOINT);
273  fStep /= ((G4double)(fPoints - 1));
274  if(fVerbose > 1) {
275  G4cout << "Emax= " << emax << " Npoints= " << fPoints
276  << " Eex= " << fExcEnergy << G4endl;
277  }
278  // integrate probabilities
279  G4double eres = (G4double)GREnergy[A];
280  G4double wres = (G4double)GRWidth[A];
281  G4double eres2= eres*eres;
282  G4double wres2= wres*wres;
283  G4double xsqr = std::sqrt(A*LevelDensity*fExcEnergy);
284 
285  G4double egam = fExcEnergy;
286  G4double gammaE2 = egam*egam;
287  G4double gammaR2 = gammaE2*wres2;
288  G4double egdp2 = gammaE2 - eres2;
289 
290  G4double p0 = G4Exp(-2.0*xsqr)*gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
291  G4double p1(0.0);
292 
293  for(G4int i=1; i<fPoints; ++i) {
294  egam -= fStep;
295  gammaE2 = egam*egam;
296  gammaR2 = gammaE2*wres2;
297  egdp2 = gammaE2 - eres2;
298  p1 = G4Exp(2.0*(std::sqrt(A*LevelDensity*std::abs(fExcEnergy - egam)) - xsqr))
299  *gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
300  fProbability += (p1 + p0);
302  //G4cout << "Egamma= " << egam << " Eex= " << fExcEnergy
303  //<< " p0= " << p0 << " p1= " << p1 << " sum= " << fCummProbability[i] <<G4endl;
304  p0 = p1;
305  }
306 
307  static const G4double NormC = 1.25*CLHEP::millibarn
309  fProbability *= fStep*NormC*A;
310  if(fVerbose > 1) { G4cout << "prob= " << fProbability << G4endl; }
311  return fProbability;
312 }
313 
314 G4double
316 {
317  G4double E = energy;
319  if(fLevelManager) {
320  E = fLevelManager->NearestLevelEnergy(energy, fIndex);
321  if(E > fLevelEnergyMax + Tolerance) { E = energy; }
322  }
323  return E;
324 }
325 
327 {
329  return fLevelEnergyMax;
330 }
331 
332 G4Fragment*
334 {
335  if(!isInitialised) { Initialise(); }
336  G4Fragment* result = nullptr;
337  G4double eexc = nucleus->GetExcitationEnergy();
338  if(eexc <= Tolerance) { return result; }
339 
340  InitialiseLevelManager(nucleus->GetZ_asInt(), nucleus->GetA_asInt());
341 
342  G4double time = nucleus->GetCreationTime();
343 
344  G4double efinal = 0.0;
345  G4double ratio = 0.0;
346  vShellNumber = -1;
347  G4int JP1 = 0;
348  G4int JP2 = 0;
349  G4int multiP = 0;
350  G4bool isGamma = true;
351  G4bool isDiscrete = false;
352 
353  const G4NucLevel* level = nullptr;
354  size_t ntrans = 0;
355 
356  // initial discrete state
357  if(fLevelManager && eexc <= fLevelEnergyMax + Tolerance) {
359  if(0 < fIndex) {
360  // for discrete transition
361  level = fLevelManager->GetLevel(fIndex);
362  if(level) {
363  ntrans = level->NumberOfTransitions();
364  JP1 = fLevelManager->SpinTwo(fIndex);
365  if(ntrans > 0) {
366  isDiscrete = true;
367  } else {
368  // if no transition available nothing is done for RDM
369  if(fRDM) {return result; }
371  --fIndex;
372  level = fLevelManager->GetLevel(fIndex);
373  ntrans = level->NumberOfTransitions();
374  JP1 = fLevelManager->SpinTwo(fIndex);
375  if(ntrans > 0) { isDiscrete = true; }
376  }
377  }
378  }
379  }
380  }
381  if(fVerbose > 1) {
382  G4int prec = G4cout.precision(4);
383  G4cout << "GenerateGamma: Z= " << nucleus->GetZ_asInt()
384  << " A= " << nucleus->GetA_asInt()
385  << " Exc= " << eexc << " Emax= "
386  << fLevelEnergyMax << " idx= " << fIndex
387  << " fCode= " << fCode << " fPoints= " << fPoints
388  << " Ntr= " << ntrans << " discrete: " << isDiscrete
389  << " fProb= " << fProbability << G4endl;
390  G4cout.precision(prec);
391  }
392 
393  // continues part
394  if(!isDiscrete) {
395  // G4cout << "Continues fIndex= " << fIndex << G4endl;
396 
397  // we compare current excitation versus value used for probability
398  // computation and also Z and A used for probability computation
399  if(fCode != 1000*theZ + theA || eexc != fExcEnergy) {
400  GetEmissionProbability(nucleus);
401  }
402  if(fProbability == 0.0) { return result; }
404  for(G4int i=1; i<fPoints; ++i) {
405  //G4cout << "y= " << y << " cummProb= " << fCummProbability[i] << G4endl;
406  if(y <= fCummProbability[i]) {
407  efinal = fStep*((i - 1) + (y - fCummProbability[i-1])
408  /(fCummProbability[i] - fCummProbability[i-1]));
409  break;
410  }
411  }
412  // final discrete level
413  if(fLevelManager) {
414  if(efinal < fLevelEnergyMax) {
415  //G4cout << "Efinal= " << efinal << " idx= " << fIndex << G4endl;
417  efinal = fLevelManager->LevelEnergy(fIndex);
418  // protection - take level below
419  if(efinal >= eexc && 0 < fIndex) {
420  --fIndex;
421  efinal = fLevelManager->LevelEnergy(fIndex);
422  }
424 
425  // not allowed to have final energy above max energy
426  // if G4LevelManager exist
427  } else {
428  efinal = fLevelEnergyMax;
430  }
431  }
432  if(fVerbose > 1) {
433  G4cout << "Continues emission efinal(MeV)= " << efinal << G4endl;
434  }
435  //discrete part
436  } else {
437  if(fVerbose > 1) {
438  G4cout << "Discrete emission from level Index= " << fIndex
439  << " Elevel= " << fLevelManager->LevelEnergy(fIndex)
440  << " RDM= " << fRDM << " ICM= " << fICM << G4endl;
441  }
442  if(0 == fIndex || !level) { return result; }
443 
444  // stable fragment has life time -1
445  // if called from radioactive decay the life time is not checked
447  if(ltime < 0.0 || (!fRDM && ltime > fMaxLifeTime)) { return result; }
448 
449  size_t idx = 0;
450  if(1 < ntrans) {
451  idx = level->SampleGammaTransition(G4UniformRand());
452  }
453  if(fVerbose > 1) {
454  G4cout << "Ntrans= " << ntrans << " idx= " << idx
455  << " ICM= " << fICM << " JP1= " << JP1 << G4endl;
456  }
457  G4double prob = (G4double)level->GammaProbability(idx);
458  // prob = 0 means that there is only internal conversion
459  if(fICM && prob < 1.0) {
460  G4double rndm = G4UniformRand();
461  if(rndm > prob) {
462  isGamma = false;
463  rndm = (rndm - prob)/(1.0 - prob);
464  vShellNumber = level->SampleShell(idx, rndm);
465  }
466  }
467  // it is discrete transition with possible gamma correlation
468  isDiscrete = true;
469  ratio = level->MultipolarityRatio(idx);
470  multiP = level->TransitionType(idx);
471  fIndex = level->FinalExcitationIndex(idx);
472  JP2 = fLevelManager->SpinTwo(fIndex);
473 
474  // final energy and time
475  efinal = fLevelManager->LevelEnergy(fIndex);
476  if(fSampleTime && ltime > 0.0) {
477  time -= ltime*G4Log(G4UniformRand());
478  }
480  }
481  // protection for floating levels
482  if(std::abs(efinal - eexc) <= Tolerance) { return result; }
483 
484  result = fTransition->SampleTransition(nucleus, efinal, ratio, JP1,
485  JP2, multiP, vShellNumber,
486  isDiscrete, isGamma);
487  if(result) { result->SetCreationTime(time); }
488 
489  // updated residual nucleus
490  nucleus->SetCreationTime(time);
491  nucleus->SetSpin(0.5*JP2);
493 
494  // ignore the floating levels with zero energy and create ground state
495  if(efinal == 0.0 && fIndex > 0) {
496  fIndex = 0;
498  }
499 
500  if(fVerbose > 1) {
501  G4cout << "Final level E= " << efinal << " time= " << time
502  << " idxFinal= " << fIndex << " isDiscrete: " << isDiscrete
503  << " isGamma: " << isGamma << " multiP= " << multiP
504  << " shell= " << vShellNumber
505  << " JP1= " << JP1 << " JP2= " << JP2 << G4endl;
506  }
507  return result;
508 }
509 
511 {
512  if(p != fTransition) {
513  delete fTransition;
514  fTransition = p;
515  }
516 }
517 
519 {
520  fICM = val;
521 }
522 
524 {
525  fRDM = val;
526 }
527 
static G4NuclearPolarizationStore * GetInstance()
G4int SampleShell(size_t idx, G4double rndm) const
Definition: G4NucLevel.hh:169
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
G4double GetCreationTime() const
Definition: G4Fragment.hh:441
static const double prec
Definition: RanecuEngine.cc:58
void SetNuclearPolarization(G4NuclearPolarization *)
Definition: G4Fragment.hh:461
virtual G4bool BreakUpChain(G4FragmentVector *theResult, G4Fragment *theNucleus) final
virtual void SetICM(G4bool)
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus) final
void SetTwoJMAX(G4int val)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:446
G4DeexPrecoParameters * GetParameters()
G4bool GetInternalConversionFlag() const
static G4float GRWidth[MAXGRDATA]
G4NuclearPolarization * GetNuclearPolarization() const
Definition: G4Fragment.hh:456
size_t SampleGammaTransition(G4double rndm) const
Definition: G4NucLevel.hh:159
#define G4endl
Definition: G4ios.hh:61
const G4int MAXGRDATA
const G4int MAXDEPOINT
Float_t y
Definition: compare.C:6
const char * p
Definition: xmltok.h:285
static G4NuclearLevelData * GetInstance()
static constexpr double hbarc
float G4float
Definition: G4Types.hh:77
static G4float GREnergy[MAXGRDATA]
G4PhotonEvaporation(G4GammaTransition *ptr=nullptr)
G4int GetA_asInt() const
Definition: G4Fragment.hh:259
static constexpr double pi2
Definition: SystemOfUnits.h:57
void SetExcitationEnergy(G4double val)
static const G4double emax
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:281
G4double GetMinExcitation() const
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
G4float MultipolarityRatio(size_t idx) const
Definition: G4NucLevel.hh:151
G4double G4Log(G4double x)
Definition: G4Log.hh:230
const G4NucLevel * GetLevel(size_t i) const
void SetPolarizationFlag(G4bool val)
static constexpr double neutron_mass_c2
Float_t Z
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:237
G4int FloatingLevel(size_t i) const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double GetMaxLifeTime() const
size_t FinalExcitationIndex(size_t idx) const
Definition: G4NucLevel.hh:114
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4double NearestLevelEnergy(G4double energy, size_t index=0) const
G4float GammaProbability(size_t idx) const
Definition: G4NucLevel.hh:135
G4NuclearPolarizationStore * fNucPStore
void SetSpin(G4double value)
Definition: G4Fragment.hh:415
G4double ComputeGroundStateMass(G4int Z, G4int A) const
Definition: G4Fragment.hh:249
void SetGammaTransition(G4GammaTransition *)
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
double energy
Definition: plottest35.C:25
double A(double temperature)
virtual void Initialise() final
static constexpr double eV
G4double LifeTime(size_t i) const
virtual G4double GetEmissionProbability(G4Fragment *theNucleus) final
G4NuclearLevelData * fNuclearLevelData
G4int GetZ_asInt() const
Definition: G4Fragment.hh:264
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
G4Fragment * GenerateGamma(G4Fragment *nucleus)
G4double G4ParticleHPJENDLHEData::G4double result
void SetFloatingLevelNumber(G4int value)
Definition: G4Fragment.hh:425
virtual G4Fragment * SampleTransition(G4Fragment *nucleus, G4double newExcEnergy, G4double mpRatio, G4int JP1, G4int JP2, G4int MP, G4int shell, G4bool isDiscrete, G4bool isGamma)
G4double GetLevelDensity() const
void RemoveMe(G4NuclearPolarization *ptr)
int G4int
Definition: G4Types.hh:78
G4int SpinTwo(size_t i) const
void SetVerbose(G4int val)
G4NuclearPolarization * FindOrBuild(G4int Z, G4int A, G4double Eexc)
virtual G4double GetFinalLevelEnergy(G4int Z, G4int A, G4double energy) final
Definition: G4Pow.hh:56
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4GammaTransition * fTransition
G4GLOB_DLL std::ostream G4cout
G4double fCummProbability[MAXDEPOINT]
const G4LevelManager * fLevelManager
virtual void RDMForced(G4bool)
size_t NumberOfTransitions() const
Definition: G4NucLevel.hh:109
G4NuclearPolarization * fPolarization
G4double LevelEnergy(size_t i) const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:276
#define DBL_MAX
Definition: templates.hh:83
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
void InitialiseLevelManager(G4int Z, G4int A)
G4int TransitionType(size_t idx) const
Definition: G4NucLevel.hh:122
static constexpr double millibarn
Definition: SystemOfUnits.h:86
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
size_t NearestLevelIndex(G4double energy, size_t index=0) const
virtual G4double GetUpperLevelEnergy(G4int Z, G4int A) final
std::mutex G4Mutex
Definition: G4Threading.hh:84