ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhotonEvaporation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 class file
29 //
30 // CERN, Geneva, Switzerland
31 //
32 // File name: G4PhotonEvaporation
33 //
34 // Author: Vladimir Ivantchenko
35 //
36 // Creation date: 20 December 2011
37 //
38 //Modifications:
39 //
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 #include "G4PhotonEvaporation.hh"
45 
47 #include "Randomize.hh"
48 #include "G4Gamma.hh"
49 #include "G4LorentzVector.hh"
50 #include "G4FragmentVector.hh"
51 #include "G4GammaTransition.hh"
52 #include "G4Pow.hh"
55 
58 
59 #ifdef G4MULTITHREADED
60 G4Mutex G4PhotonEvaporation::PhotonEvaporationMutex = G4MUTEX_INITIALIZER;
61 #endif
62 
64  : fLevelManager(nullptr), fTransition(p), fPolarization(nullptr),
65  fVerbose(1), fPoints(0), vShellNumber(-1), fIndex(0),
66  fMaxLifeTime(DBL_MAX),
67  fICM(true), fRDM(false), fSampleTime(true),
68  fCorrelatedGamma(false), isInitialised(false)
69 {
70  //G4cout << "### New G4PhotonEvaporation() " << this << G4endl;
72  Tolerance = 20*CLHEP::eV;
73 
74  if(!fTransition) { fTransition = new G4GammaTransition(); }
75 
76  theA = theZ = fCode = 0;
78 
79  for(G4int i=0; i<MAXDEPOINT; ++i) { fCummProbability[i] = 0.0; }
80  if(0.0f == GREnergy[1]) { InitialiseGRData(); }
81 }
82 
84 {
85  delete fTransition;
86 }
87 
89 {
90  if(isInitialised) { return; }
91  isInitialised = true;
92 
94  Tolerance = param->GetMinExcitation();
95  fMaxLifeTime = param->GetMaxLifeTime();
98  fVerbose = param->GetVerbose();
99 
101  fTransition->SetTwoJMAX(param->GetTwoJMAX());
103  if(fVerbose > 1) {
104  G4cout << "### G4PhotonEvaporation is initialized " << this << G4endl;
105  }
106 }
107 
109 {
110 #ifdef G4MULTITHREADED
111  G4MUTEXLOCK(&G4PhotonEvaporation::PhotonEvaporationMutex);
112 #endif
113  if(0.0f == GREnergy[1]) {
114  G4Pow* g4calc = G4Pow::GetInstance();
115  const G4float GRWfactor = 0.3f;
116  for (G4int A=1; A<MAXGRDATA; ++A) {
117  GREnergy[A] = (G4float)(40.3*CLHEP::MeV/g4calc->powZ(A,0.2));
118  GRWidth[A] = GRWfactor*GREnergy[A];
119  }
120  }
121 #ifdef G4MULTITHREADED
122  G4MUTEXUNLOCK(&G4PhotonEvaporation::PhotonEvaporationMutex);
123 #endif
124 }
125 
126 G4Fragment*
128 {
129  if(!isInitialised) { Initialise(); }
130  fSampleTime = (fRDM) ? false : true;
131 
132  // potentially external code may set initial polarization
133  // but only for radioactive decay nuclear polarization is considered
134  G4NuclearPolarizationStore* fNucPStore = nullptr;
135  if(fCorrelatedGamma && fRDM) {
137  if(nucleus->GetNuclearPolarization()) {
138  fNucPStore->RemoveMe(nucleus->GetNuclearPolarization());
139  delete nucleus->GetNuclearPolarization();
140  }
141  fPolarization = fNucPStore->FindOrBuild(nucleus->GetZ_asInt(),
142  nucleus->GetA_asInt(),
143  nucleus->GetExcitationEnergy());
145  }
146  if(fVerbose > 2) {
147  G4cout << "G4PhotonEvaporation::EmittedFragment: "
148  << *nucleus << G4endl;
149  if(fPolarization) { G4cout << "NucPolar: " << fPolarization << G4endl; }
150  G4cout << " CorrGamma: " << fCorrelatedGamma << " RDM: " << fRDM
151  << " fPolarization: " << fPolarization << G4endl;
152  }
153  G4Fragment* gamma = GenerateGamma(nucleus);
154 
155  // remove G4NuclearPolarizaton when reach ground state
156  if(fNucPStore && fPolarization && 0 == fIndex) {
157  if(fVerbose > 3) {
158  G4cout << "G4PhotonEvaporation::EmittedFragment: remove "
159  << fPolarization << G4endl;
160  }
161  fNucPStore->RemoveMe(fPolarization);
162  fPolarization = nullptr;
164  }
165 
166  if(fVerbose > 2) {
167  G4cout << "G4PhotonEvaporation::EmittedFragment: RDM= "
168  << fRDM << " done:" << G4endl;
169  if(gamma) { G4cout << *gamma << G4endl; }
170  G4cout << " Residual: " << *nucleus << G4endl;
171  }
172  return gamma;
173 }
174 
177 {
178  if(fVerbose > 1) {
179  G4cout << "G4PhotonEvaporation::BreakItUp" << G4endl;
180  }
181  G4Fragment* aNucleus = new G4Fragment(nucleus);
182  G4FragmentVector* products = new G4FragmentVector();
183  BreakUpChain(products, aNucleus);
184  products->push_back(aNucleus);
185  return products;
186 }
187 
189  G4Fragment* nucleus)
190 {
191  if(!isInitialised) { Initialise(); }
192  if(fVerbose > 1) {
193  G4cout << "G4PhotonEvaporation::BreakUpChain RDM= " << fRDM << " "
194  << *nucleus << G4endl;
195  }
196  G4Fragment* gamma = nullptr;
197  fSampleTime = (fRDM) ? false : true;
198 
199  // start decay chain from unpolarized state
200  if(fCorrelatedGamma) {
202  nucleus->GetA_asInt(),
203  nucleus->GetExcitationEnergy());
205  }
206 
207  do {
208  gamma = GenerateGamma(nucleus);
209  if(gamma) {
210  products->push_back(gamma);
211  if(fVerbose > 2) {
212  G4cout << "G4PhotonEvaporation::BreakUpChain: "
213  << *gamma << G4endl;
214  G4cout << " Residual: " << *nucleus << G4endl;
215  }
216  // for next decays in the chain always sample time
217  fSampleTime = true;
218  }
219  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
220  } while(gamma);
221 
222  // clear nuclear polarization end of chain
223  if(fPolarization) {
224  delete fPolarization;
225  fPolarization = nullptr;
227  }
228  return false;
229 }
230 
231 G4double
233 {
234  if(!isInitialised) { Initialise(); }
235  fProbability = 0.0;
236  fExcEnergy = nucleus->GetExcitationEnergy();
237  G4int Z = nucleus->GetZ_asInt();
238  G4int A = nucleus->GetA_asInt();
239  fCode = 1000*Z + A;
240  if(fVerbose > 2) {
241  G4cout << "G4PhotonEvaporation::GetEmissionProbability: Z="
242  << Z << " A=" << A << " Eexc(MeV)= " << fExcEnergy << G4endl;
243  }
244  // ignore gamma de-excitation for exotic fragments
245  // and for very low excitations
246  if(0 >= Z || 1 >= A || Z == A || Tolerance >= fExcEnergy)
247  { return fProbability; }
248 
249  // ignore gamma de-excitation for highly excited levels
250  if(A >= MAXGRDATA) { A = MAXGRDATA-1; }
251  //G4cout<<" GREnergy= "<< GREnergy[A]<<" GRWidth= "<<GRWidth[A]<<G4endl;
252 
253  static const G4float GREfactor = 5.0f;
254  if(fExcEnergy >= (G4double)(GREfactor*GRWidth[A] + GREnergy[A])) {
255  return fProbability;
256  }
257  // probability computed assuming continium transitions
258  // VI: continium transition are limited only to final states
259  // below Fermi energy (this approach needs further evaluation)
260  G4double emax = std::max(0.0, nucleus->ComputeGroundStateMass(Z, A-1)
262 
263  // max energy level for continues transition
264  emax = std::min(emax, fExcEnergy);
265  const G4double eexcfac = 0.99;
266  if(0.0 == emax || fExcEnergy*eexcfac <= emax) { emax = fExcEnergy*eexcfac; }
267 
268  fStep = emax;
269  const G4double MaxDeltaEnergy = CLHEP::MeV;
270  fPoints = std::min((G4int)(fStep/MaxDeltaEnergy) + 2, MAXDEPOINT);
271  fStep /= ((G4double)(fPoints - 1));
272  if(fVerbose > 2) {
273  G4cout << "Emax= " << emax << " Npoints= " << fPoints
274  << " Eex= " << fExcEnergy << G4endl;
275  }
276  // integrate probabilities
277  G4double eres = (G4double)GREnergy[A];
278  G4double wres = (G4double)GRWidth[A];
279  G4double eres2= eres*eres;
280  G4double wres2= wres*wres;
282  G4double xsqr = std::sqrt(levelDensity*fExcEnergy);
283 
284  G4double egam = fExcEnergy;
285  G4double gammaE2 = egam*egam;
286  G4double gammaR2 = gammaE2*wres2;
287  G4double egdp2 = gammaE2 - eres2;
288 
289  G4double p0 = G4Exp(-2.0*xsqr)*gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
290  G4double p1(0.0);
291 
292  for(G4int i=1; i<fPoints; ++i) {
293  egam -= fStep;
294  gammaE2 = egam*egam;
295  gammaR2 = gammaE2*wres2;
296  egdp2 = gammaE2 - eres2;
297  p1 = G4Exp(2.0*(std::sqrt(levelDensity*std::abs(fExcEnergy - egam)) - xsqr))
298  *gammaR2*gammaE2/(egdp2*egdp2 + gammaR2);
299  fProbability += (p1 + p0);
301  if(fVerbose > 3) {
302  G4cout << "Egamma= " << egam << " Eex= " << fExcEnergy
303  << " p0= " << p0 << " p1= " << p1 << " sum= "
304  << fCummProbability[i] <<G4endl;
305  }
306  p0 = p1;
307  }
308 
309  static const G4double NormC = 1.25*CLHEP::millibarn
311  fProbability *= fStep*NormC*A;
312  if(fVerbose > 1) { G4cout << "prob= " << fProbability << G4endl; }
313  return fProbability;
314 }
315 
316 G4double
318 {
319  G4double E = energy;
321  if(fLevelManager) {
322  E = fLevelManager->NearestLevelEnergy(energy, fIndex);
323  if(E > fLevelEnergyMax + Tolerance) { E = energy; }
324  }
325  return E;
326 }
327 
329 {
331  return fLevelEnergyMax;
332 }
333 
334 G4Fragment*
336 {
337  if(!isInitialised) { Initialise(); }
338  G4Fragment* result = nullptr;
339  G4double eexc = nucleus->GetExcitationEnergy();
340  if(eexc <= Tolerance) { return result; }
341 
342  InitialiseLevelManager(nucleus->GetZ_asInt(), nucleus->GetA_asInt());
343 
344  G4double time = nucleus->GetCreationTime();
345 
346  G4double efinal = 0.0;
347  G4double ratio = 0.0;
348  vShellNumber = -1;
349  G4int JP1 = 0;
350  G4int JP2 = 0;
351  G4int multiP = 0;
352  G4bool isGamma = true;
353  G4bool isDiscrete = false;
354 
355  const G4NucLevel* level = nullptr;
356  size_t ntrans = 0;
357 
358  if(fVerbose > 2) {
359  G4cout << "GenerateGamma: " << " Eex= " << eexc
360  << " Eexmax= " << fLevelEnergyMax << G4endl;
361  }
362  // initial discrete state
363  if(fLevelManager && eexc <= fLevelEnergyMax + Tolerance) {
365  // initial state below 1st level
366  if(0 == fIndex && eexc >= Tolerance
367  && fLevelManager->NumberOfTransitions() > 0) { fIndex = 1; }
368  isDiscrete = true;
369  if(fVerbose > 2) {
370  G4cout << " index= " << fIndex
371  << " lTime= " << fLevelManager->LifeTime(fIndex) << G4endl;
372  }
373  if(0 < fIndex) {
374  // for discrete transition
375  level = fLevelManager->GetLevel(fIndex);
376  if(level) {
377  ntrans = level->NumberOfTransitions();
378  JP1 = fLevelManager->SpinTwo(fIndex);
379  if(fVerbose > 2) {
380  G4cout << " ntrans= " << ntrans << " JP= " << JP1
381  << " RDM: " << fRDM << G4endl;
382  }
383  if(0 == ntrans && fLevelManager->FloatingLevel(fIndex) > 0) {
384  --fIndex;
385  level = fLevelManager->GetLevel(fIndex);
386  ntrans = level->NumberOfTransitions();
387  JP1 = fLevelManager->SpinTwo(fIndex);
388  }
389  }
390  }
391  }
392  if(fVerbose > 2) {
393  G4int prec = G4cout.precision(4);
394  G4cout << "GenerateGamma: Z= " << nucleus->GetZ_asInt()
395  << " A= " << nucleus->GetA_asInt()
396  << " Exc= " << eexc << " Emax= "
397  << fLevelEnergyMax << " idx= " << fIndex
398  << " fCode= " << fCode << " fPoints= " << fPoints
399  << " Ntr= " << ntrans << " discrete: " << isDiscrete
400  << " fProb= " << fProbability << G4endl;
401  G4cout.precision(prec);
402  }
403 
404  // continues part
405  if(!isDiscrete) {
406  // we compare current excitation versus value used for probability
407  // computation and also Z and A used for probability computation
408  if(fCode != 1000*theZ + theA || eexc != fExcEnergy) {
409  GetEmissionProbability(nucleus);
410  }
411  if(fProbability == 0.0) {
412  fPoints = 1;
413  efinal = 0.0;
414  } else {
416  for(G4int i=1; i<fPoints; ++i) {
417  if(fVerbose > 3) {
418  G4cout << "y= " << y << " cummProb= " << fCummProbability[i]
419  << " fPoints= " << fPoints << " fStep= " << fStep << G4endl;
420  }
421  if(y <= fCummProbability[i]) {
422  efinal = fStep*((i - 1) + (y - fCummProbability[i-1])
423  /(fCummProbability[i] - fCummProbability[i-1]));
424  break;
425  }
426  }
427  }
428  // final discrete level
429  if(fVerbose > 2) {
430  G4cout << "Continues proposes Efinal= " << efinal << G4endl;
431  }
432  if(fLevelManager) {
433  if(efinal < fLevelEnergyMax) {
435  efinal = fLevelManager->LevelEnergy(fIndex);
436  // protection - take level below
437  if(efinal >= eexc && 0 < fIndex) {
438  --fIndex;
439  efinal = fLevelManager->LevelEnergy(fIndex);
440  }
442 
443  // not allowed to have final energy above max energy
444  // if G4LevelManager exist
445  } else {
446  efinal = fLevelEnergyMax;
448  }
449  }
450  if(fVerbose > 2) {
451  G4cout << "Continues emission efinal(MeV)= " << efinal << G4endl;
452  }
453  //discrete part ground state
454  } else if(0 == fIndex) {
455  return result;
456 
457  //discrete part
458  } else {
459 
460  if(fVerbose > 2) {
461  G4cout << "Discrete emission from level Index= " << fIndex
462  << " Elevel= " << fLevelManager->LevelEnergy(fIndex)
463  << " Ltime= " << fLevelManager->LifeTime(fIndex)
464  << " LtimeMax= " << fMaxLifeTime
465  << " RDM= " << fRDM << " ICM= " << fICM << G4endl;
466  }
467 
468  // stable fragment has life time -1 or above the limit
469  // if is called from the radioactive decay the life time is not checked
471  if(ltime < 0.0 || (!fRDM && ltime > fMaxLifeTime)) { return result; }
472 
473  // no transitions: force transition to the closest level
474  if(0 == ntrans) {
475  G4int ii = fIndex - 1;
476  for(; ii>0; --ii) {
477  const G4NucLevel* fl = fLevelManager->GetLevel(ii);
478  if(fl && 0 < fl->NumberOfTransitions()) { break; }
479  }
480  fIndex = ii;
481  } else {
482  size_t idx = 0;
483  if(1 < ntrans) {
484  idx = level->SampleGammaTransition(G4UniformRand());
485  }
486  if(fVerbose > 2) {
487  G4cout << "Ntrans= " << ntrans << " idx= " << idx
488  << " ICM= " << fICM << " JP1= " << JP1 << G4endl;
489  }
490  G4double prob = (G4double)level->GammaProbability(idx);
491  // prob = 0 means that there is only internal conversion
492  if(fICM && prob < 1.0) {
493  G4double rndm = G4UniformRand();
494  if(rndm > prob) {
495  isGamma = false;
496  rndm = (rndm - prob)/(1.0 - prob);
497  vShellNumber = level->SampleShell(idx, rndm);
498  }
499  }
500  // it is discrete transition with possible gamma correlation
501  ratio = level->MultipolarityRatio(idx);
502  multiP = level->TransitionType(idx);
503  fIndex = level->FinalExcitationIndex(idx);
504  }
505  JP2 = fLevelManager->SpinTwo(fIndex);
506 
507  // final energy and time
508  efinal = fLevelManager->LevelEnergy(fIndex);
509  if(fSampleTime && ltime > 0.0) {
510  time -= ltime*G4Log(G4UniformRand());
511  }
513  }
514  // protection for floating levels
515  if(std::abs(efinal - eexc) <= Tolerance) { return result; }
516 
517  result = fTransition->SampleTransition(nucleus, efinal, ratio, JP1,
518  JP2, multiP, vShellNumber,
519  isDiscrete, isGamma);
520  if(result) { result->SetCreationTime(time); }
521 
522  // updated residual nucleus
523  nucleus->SetCreationTime(time);
524  nucleus->SetSpin(0.5*JP2);
526 
527  // ignore the floating levels with zero energy and create ground state
528  if(efinal == 0.0 && fIndex > 0) {
529  fIndex = 0;
531  }
532 
533  if(fVerbose > 2) {
534  G4cout << "Final level E= " << efinal << " time= " << time
535  << " idxFinal= " << fIndex << " isDiscrete: " << isDiscrete
536  << " isGamma: " << isGamma << " multiP= " << multiP
537  << " shell= " << vShellNumber
538  << " JP1= " << JP1 << " JP2= " << JP2 << G4endl;
539  }
540  return result;
541 }
542 
544 {
545  if(p != fTransition) {
546  delete fTransition;
547  fTransition = p;
548  }
549 }
550 
552 {
553  fICM = val;
554 }
555 
557 {
558  fRDM = val;
559 }
560