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