Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4FermiFragmentsPoolVI.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: G4FermiFragmentsPoolVI.cc,v 1.5 2006-06-29 20:13:13 gunter Exp $
27 //
28 // FermiBreakUp de-excitation model
29 // by V. Ivanchenko (July 2016)
30 //
31 
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4NuclearLevelData.hh"
36 #include "G4LevelManager.hh"
37 #include "G4DeexPrecoParameters.hh"
38 #include <iomanip>
39 
41 {
42  // G4cout << "### G4FermiFragmentsPoolVI is constructed" << G4endl;
43  G4DeexPrecoParameters* param =
46  timelim = (G4float)param->GetMaxLifeTime();
47 
48  elim = 10*CLHEP::MeV;
49  elimf= (G4float)elim;
50 
51  fragment_pool.reserve(399);
52  funstable.reserve(80);
53 
54  Initialise();
55 }
56 
58 {
59  size_t nn;
60  for(G4int i=0; i<maxA; ++i) {
61  nn = list_p[i].size();
62  for(size_t j=0; j<nn; ++j) { delete (list_p[i])[j]; }
63  nn = list_c[i].size();
64  for(size_t j=0; j<nn; ++j) { delete (list_c[i])[j]; }
65  nn = list_d[i].size();
66  for(size_t j=0; j<nn; ++j) { delete (list_d[i])[j]; }
67  nn = list_u[i].size();
68  for(size_t j=0; j<nn; ++j) { delete (list_u[i])[j]; }
69  }
70  nn = fragment_pool.size();
71  for(size_t j=0; j<nn; ++j) { delete fragment_pool[j]; }
72  nn = funstable.size();
73  for(size_t j=0; j<nn; ++j) { delete funstable[j]; }
74 }
75 
76 G4bool
78 {
79  G4bool isInList = false;
80  size_t nn = list_f[A].size();
81  for(size_t i=0; i<nn; ++i) {
82  if(Z == (list_f[A])[i]->GetZ()) {
83  isInList = true;
84  if(etot <= (list_f[A])[i]->GetFragmentMass() + elim) { return true; }
85  }
86  }
87  if(isInList) { return false; }
88  nn = list_g[A].size();
89  for(size_t i=0; i<nn; ++i) {
90  if(Z == (list_g[A])[i]->GetZ() &&
91  etot <= (list_g[A])[i]->GetFragmentMass() + elim) { return true; }
92  }
93  return false;
94 }
95 
96 const G4FermiChannels*
98 {
99  const G4FermiChannels* res = nullptr;
100  G4double demax = e;
101 
102  // stable channels;
103  size_t nn = list_c[A].size();
104  for(size_t j=0; j<nn; ++j) {
105  const G4FermiFragment* frag = (list_f[A])[j];
106  if(frag->GetZ() != Z) { continue; }
107  G4double de = e - frag->GetTotalEnergy();
108  //G4cout << " Stab check " << j << " channel de= " << de << " tol= " << tolerance << G4endl;
109  if(std::abs(de) <= tolerance) {
110  res = (list_c[A])[j];
111  //G4cout << " Stab chan: " << j << " N= " << res->GetNumberOfChannels() << G4endl;
112  break;
113  } else if(de + tolerance > 0.0) {
114  if(de < demax) {
115  res = (list_c[A])[j];
116  demax = de;
117  //G4cout << " Stab chan: " << j << " N= " << res->GetNumberOfChannels() << G4endl;
118  }
119  }
120  }
121  // unstable channels
122  if(!res) {
123  nn = list_d[A].size();
124  for(size_t j=0; j<nn; ++j) {
125  const G4FermiFragment* frag = (list_g[A])[j];
126  if(frag->GetZ() != Z) { continue; }
127  G4double de = e - frag->GetTotalEnergy();
128  //G4cout << " Unst check " << j << " channel de= " << de << " tol= " << tolerance << G4endl;
129  if(std::abs(de) <= tolerance || de > 0.0) {
130  res = (list_d[A])[j];
131  //G4cout << " Unst chan No: " << j << " N= " << res->GetNumberOfChannels() << G4endl;
132  break;
133  }
134  }
135  }
136  return res;
137 }
138 
140 {
141  G4bool res = false;
142  G4int nn = list_f[A].size();
143  for(G4int i=0; i<nn; ++i) {
144  if((list_f[A])[i]->GetZ() == Z) {
145  res = true;
146  break;
147  }
148  }
149  return res;
150 }
151 
153  G4double exc) const
154 {
155  G4bool res = false;
156  G4int nfrag = fragment_pool.size();
157  for(G4int i=0; i<nfrag; ++i) {
158  const G4FermiFragment* fr = fragment_pool[i];
159  if(fr->GetZ() == Z && fr->GetA() == A &&
160  std::abs(exc - fr->GetExcitationEnergy()) < tolerance) {
161  res = true;
162  break;
163  }
164  }
165  return res;
166 }
167 
169  const G4FermiFragment* f1, const G4FermiFragment* f2) const
170 {
171  G4bool res = false;
172  G4int A1 = f1->GetA();
173  G4int A2 = f2->GetA();
174  G4int A = A1 + A2;
175  G4int nn = list_p[A].size();
176  for(G4int i=0; i<nn; ++i) {
177  if(f1 == (list_p[A])[i]->GetFragment1() &&
178  f2 == (list_p[A])[i]->GetFragment2()) {
179  res = true;
180  break;
181  }
182  }
183  return res;
184 }
185 
187  const G4FermiFragment* f1, const G4FermiFragment* f2) const
188 {
189  G4bool res = false;
190  G4int A1 = f1->GetA();
191  G4int A2 = f2->GetA();
192  G4int A = A1 + A2;
193  G4int nn = list_u[A].size();
194  for(G4int i=0; i<nn; ++i) {
195  if(f1 == (list_u[A])[i]->GetFragment1() &&
196  f2 == (list_u[A])[i]->GetFragment2()) {
197  res = true;
198  break;
199  }
200  }
201  return res;
202 }
203 
205 {
206  static const G4int nmin = 8;
207 
208  // stable particles
209  fragment_pool.push_back(new G4FermiFragment(1, 0, 1, 0.0, true, true));
210  fragment_pool.push_back(new G4FermiFragment(1, 1, 1, 0.0, true, true));
211  fragment_pool.push_back(new G4FermiFragment(2, 1, 2, 0.0, true, true));
212  fragment_pool.push_back(new G4FermiFragment(3, 1, 1, 0.0, true, true));
213  fragment_pool.push_back(new G4FermiFragment(3, 2, 1, 0.0, true, true));
214  fragment_pool.push_back(new G4FermiFragment(4, 2, 0, 0.0, true, true));
215  fragment_pool.push_back(new G4FermiFragment(5, 2, 3, 0.0, true, true));
216  fragment_pool.push_back(new G4FermiFragment(5, 3, 3, 0.0, true, true));
217 
218  // use level data and construct the pool
220  for(G4int Z=1; Z<maxZ; ++Z) {
221  G4int Amin = ndata->GetMinA(Z);
222  G4int Amax = std::min(maxA, ndata->GetMaxA(Z)+1);
223  for(G4int A=Amin; A<Amax; ++A) {
224  const G4LevelManager* man = ndata->GetLevelManager(Z, A);
225  if(man) {
226  size_t nn = man->NumberOfTransitions();
227  // very unstable state
228  if(ndata->MaxLevelEnergy(Z, A) == 0.0f && man->LifeTime(0) == 0.0f) {
229  continue;
230  }
231  for(size_t i=0; i<=nn; ++i) {
232  G4float exc = man->LevelEnergy(i);
233  // only levels below limit are consided
234  if(exc >= elimf) { continue; }
235  G4double excd = (G4double)exc;
236  // only new are considered
237  if(IsInThePool(Z, A, excd)) { continue; }
238  G4float ltime = man->LifeTime(i);
239  G4bool stable = (ltime < 0.0f || ltime > timelim) ? true : false;
240  //G4cout << "Z= " << Z << " A= " << A << " Eex= " << exc
241  // << " t= " << ltime << " tlim= " << timelim
242  // << " stable: "<< stable << G4endl;
243  fragment_pool.push_back(new G4FermiFragment(A,Z,man->SpinTwo(i),excd,stable,true));
244  }
245  }
246  }
247  }
248  // prepare structures per A for normal fragments
249  static const size_t lfmax[maxA] = {
250  0, 2, 1, 2, 1, 2, 6, 14, 16, 22, 45, 53, 37, 44, 33, 58, 63};
251  for(G4int A=1; A<maxA; ++A) {
252  list_f[A].reserve(lfmax[A]);
253  list_c[A].reserve(lfmax[A]);
254  }
255  static const size_t lfch[maxA] = {
256  0, 0, 0, 0, 0, 1, 2, 5, 6, 3, 12, 8, 4, 10, 1, 8, 6};
257  G4int nfrag = fragment_pool.size();
258  for(G4int i=0; i<nfrag; ++i) {
259  const G4FermiFragment* f = fragment_pool[i];
260  G4int A = f->GetA();
261  G4double exc = f->GetExcitationEnergy();
262  list_f[A].push_back(f);
263  list_c[A].push_back(new G4FermiChannels(lfch[A], exc, f->GetTotalEnergy()));
264  }
265 
266  // list of unphysical fragments
267  static const size_t lfun[maxA] = {
268  0, 0, 2, 2, 4, 4, 6, 6, 6, 6, 6, 6, 7, 6, 7, 6, 6};
269  for(G4int A=1; A<maxA; ++A) {
270  list_g[A].reserve(lfun[A]);
271  list_d[A].reserve(lfun[A]);
272  }
273 
274  static const size_t luch[maxA] = {
275  0, 0, 1, 1, 3, 3, 5, 6, 6, 8, 8, 8, 8, 8, 8, 8, 8};
276  for(G4int Z=0; Z<maxZ; ++Z) {
277  G4int A0 = std::max(Z, 1);
278  for(G4int A=A0; A<maxA; ++A) {
279  if(IsInThePool(Z, A, 0.0)) { continue; }
280  const G4FermiFragment* f = new G4FermiFragment(A, Z, -1, 0.0, false, false);
281  funstable.push_back(f);
282  list_g[A].push_back(f);
283  list_d[A].push_back(new G4FermiChannels(luch[A],0.0,f->GetTotalEnergy()));
284  }
285  }
286  static const size_t pphm[maxA] = {
287  0, 0, 2, 2, 4, 4, 6, 6, 6, 6, 6, 6, 7, 6, 7, 6, 6};
288  static const size_t punm[maxA] = {
289  0, 0, 2, 2, 8, 10, 22, 24, 24, 31, 35, 36, 44, 36, 44, 36, 36};
290  for(G4int A=1; A<maxA; ++A) {
291  list_p[A].reserve(pphm[A]);
292  list_u[A].reserve(punm[A]);
293  }
294  /*
295  G4cout << "G4FermiFragmentsPoolVI::Initialise main loop @@@@@@"
296  << " Nfrag= " << nfrag << " Nuns= " << funstable.size() << G4endl;
297  */
298  // list of fragment pairs ordered by A
299  for(G4int i=0; i<nfrag; ++i) {
300  const G4FermiFragment* f1 = fragment_pool[i];
301  G4int Z1 = f1->GetZ();
302  G4int A1 = f1->GetA();
303  G4double e1 = f1->GetTotalEnergy();
304  for(G4int j=0; j<nfrag; ++j) {
305  const G4FermiFragment* f2 = fragment_pool[j];
306  G4int Z2 = f2->GetZ();
307  G4int A2 = f2->GetA();
308  if(A2 < A1 || (A2 == A1 && Z2 < Z1)) { continue; }
309  G4int Z = Z1 + Z2;
310  G4int A = A1 + A2;
311 
312  if(Z >= maxZ || A >= maxA ||
313  IsInPhysPairs(f1, f2) || IsInUnphysPairs(f1, f2)) { continue; }
314 
315  G4double e2 = f2->GetTotalEnergy();
316  G4double minE = e1 + e2;
317  G4double exc = 0.0;
318  if(IsPhysical(Z, A)) {
319  minE += f1->GetCoulombBarrier(A2, Z2, 0.0);
320  exc = minE - G4NucleiProperties::GetNuclearMass(A, Z);
321  }
322  /*
323  G4cout << "Z= " << Z << " A= " << A
324  << " Z1= " << Z1 << " A1= " << A1
325  << " Z2= " << Z2 << " A2= " << A2 << " Eex= " << exc
326  << " Qb= " << f1->GetCoulombBarrier(A2, Z2, 0.0)
327  << " " << e1
328  << " " << e2
329  << " " << G4NucleiProperties::GetNuclearMass(A, Z)
330  << G4endl;
331  */
332  // ignore very excited case
333  if(exc >= elim) { continue; }
334  G4FermiPair* fpair = nullptr;
335  G4int kmax = list_f[A].size();
336  for(G4int k=0; k<kmax; ++k) {
337  const G4FermiFragment* f3 = (list_f[A])[k];
338  if(Z == f3->GetZ() &&
339  f3->GetTotalEnergy() - minE + tolerance >= 0.0) {
340  if(!fpair) {
341  fpair = new G4FermiPair(f1, f2);
342  list_p[A].push_back(fpair);
343  }
344  (list_c[A])[k]->AddChannel(fpair);
345  }
346  }
347  if(fpair) { continue; }
348 
349  kmax = list_g[A].size();
350  for(G4int k=0; k<kmax; ++k) {
351  if((list_d[A])[k]->GetNumberOfChannels() >= nmin) { continue; }
352  const G4FermiFragment* f3 = (list_g[A])[k];
353  /*
354  if(Z==0) {
355  G4cout << "%%% A= " << A << " isStable: " << f3->IsStable()
356  << " de= " << f3->GetTotalEnergy() - minE << G4endl;
357  }
358  */
359  if(Z == f3->GetZ() &&
360  f3->GetTotalEnergy() - minE + tolerance >= 0.0) {
361  if(!fpair) {
362  fpair = new G4FermiPair(f1, f2);
363  list_u[A].push_back(fpair);
364  }
365  (list_d[A])[k]->AddChannel(fpair);
366  /*
367  if(Z==0) {
368  G4cout << " isAdded Unstable" << G4endl;
369  }
370  */
371  }
372  }
373  }
374  }
375 
376  // G4cout << "@@@@@@ sec loop @@@@@@" << G4endl;
377  // list of fragment pairs (stable+unstable) ordered by A
378  G4int unphys = funstable.size();
379  for(G4int i=0; i<nfrag; ++i) {
380  const G4FermiFragment* f1 = fragment_pool[i];
381  G4int Z1 = f1->GetZ();
382  G4int A1 = f1->GetA();
383  G4double e1 = f1->GetTotalEnergy();
384  for(G4int j=0; j<unphys; ++j) {
385  const G4FermiFragment* f2 = funstable[j];
386  G4int Z2 = f2->GetZ();
387  G4int A2 = f2->GetA();
388  G4int Z = Z1 + Z2;
389  G4int A = A1 + A2;
390  if(Z >= maxZ || A >= maxA || IsInUnphysPairs(f1, f2) || IsPhysical(Z, A))
391  { continue; }
392  G4double e2 = f2->GetTotalEnergy();
393  G4double minE = e1 + e2;
394  /*
395  G4cout << "Z= " << Z << " A= " << A << " Z1= " << Z1 << " A1= " << A1
396  << " Z2= " << Z2 << " A2= " << A2 << G4endl;
397  */
398  // check if this is the list of stable pairs
399  G4FermiPair* fpair = nullptr;
400  /*
401  G4int kmax = list_f[A].size();
402  for(G4int k=0; k<kmax; ++k) {
403  const G4FermiFragment* f3 = (list_f[A])[k];
404  if(Z == f3->GetZ() && f3->IsStable() &&
405  (list_c[A])[k]->GetNumberOfChannels() == 0 &&
406  f3->GetTotalEnergy() >= minE)
407  {
408  if(!fpair) {
409  fpair = new G4FermiPair(f1, f2);
410  list_u[A].push_back(fpair);
411  }
412  (list_c[A])[k]->AddChannel(fpair);
413  }
414  }
415  if(fpair) { continue; }
416  */
417  // check unphysics list
418  G4int kmax = list_g[A].size();
419  for(G4int k=0; k<kmax; ++k) {
420  const G4FermiFragment* f3 = (list_g[A])[k];
421  /*
422  if(Z == f3->GetZ())
423  G4cout << " Unst+ST k= " << k << " Z= " << f3->GetZ()
424  << " isStab " << f3->IsStable()
425  << " Nch= " << (list_d[A])[k]->GetNumberOfChannels()
426  << " de= " << f3->GetTotalEnergy() - minE
427  << G4endl;
428  */
429  if(Z == f3->GetZ() &&
430  (list_d[A])[k]->GetNumberOfChannels() < nmin &&
431  f3->GetTotalEnergy() - minE + tolerance >= 0.0)
432  {
433  fpair = new G4FermiPair(f1, f2);
434  list_u[A].push_back(fpair);
435  (list_d[A])[k]->AddChannel(fpair);
436  // G4cout << " isAdded Unstable" << G4endl;
437  break;
438  }
439  }
440  }
441  }
442 
443  // compute static probabilities
444  for(G4int A=1; A<maxA; ++A) {
445  for(size_t j=0; j<list_c[A].size(); ++j) {
446  G4FermiChannels* ch = (list_c[A])[j];
447  const G4FermiFragment* frag = (list_f[A])[j];
448  size_t nch = ch->GetNumberOfChannels();
449  if(1 < nch) {
450  std::vector<G4double>& prob = ch->GetProbabilities();
451  const std::vector<const G4FermiPair*>& pairs = ch->GetChannels();
452  G4double ptot = 0.0;
453  for(size_t i=0; i<nch; ++i) {
454  ptot += theDecay.ComputeProbability(frag->GetZ(), frag->GetA(),
455  frag->GetSpin(),
456  frag->GetTotalEnergy(),
457  pairs[i]->GetFragment1(),
458  pairs[i]->GetFragment2());
459  prob[i] = ptot;
460  }
461  if(0.0 == ptot) {
462  prob[0] = 1.0;
463  } else {
464  ptot = 1./ptot;
465  for(size_t i=0; i<nch-1; ++i) { prob[i] *= ptot; }
466  prob[nch-1] = 1.0;
467  }
468  }
469  }
470  }
471  for(G4int A=1; A<maxA; ++A) {
472  for(size_t j=0; j<list_d[A].size(); ++j) {
473  G4FermiChannels* ch = (list_d[A])[j];
474  const G4FermiFragment* frag = (list_g[A])[j];
475  size_t nch = ch->GetNumberOfChannels();
476  if(1 < nch) {
477  std::vector<G4double>& prob = ch->GetProbabilities();
478  const std::vector<const G4FermiPair*>& pairs = ch->GetChannels();
479  G4double ptot = 0.0;
480  for(size_t i=0; i<nch; ++i) {
481  ptot += theDecay.ComputeProbability(frag->GetZ(), frag->GetA(),
482  frag->GetSpin(),
483  frag->GetTotalEnergy(),
484  pairs[i]->GetFragment1(),
485  pairs[i]->GetFragment2());
486  prob[i] = ptot;
487  }
488  if(0.0 == ptot) {
489  prob[0] = 1.0;
490  } else {
491  ptot = 1./ptot;
492  for(size_t i=0; i<nch-1; ++i) { prob[i] *= ptot; }
493  prob[nch-1] = 1.0;
494  }
495  }
496  }
497  }
498 }
499 
501 {
502  if(f) {
503  G4int prec = G4cout.precision(6);
504  G4cout << " Z= " << f->GetZ() << " A= " << std::setw(2) << f->GetA()
505  << " Mass(GeV)= " << std::setw(8) << f->GetFragmentMass()/GeV
506  << " Eexc(MeV)= " << std::setw(7) << f->GetExcitationEnergy()
507  << " 2s= " << f->GetSpin() << " IsStable: " << f->IsStable()
508  << " IsPhys: " << f->IsPhysical() << G4endl;
509  G4cout.precision(prec);
510  }
511 }
512 
514 {
515  G4cout <<"----------------------------------------------------------------"
516  <<G4endl;
517  G4cout << "##### List of Fragments in the Fermi Fragment Pool #####"
518  << G4endl;
519  G4int nfrag = fragment_pool.size();
520  G4cout << " For stable " << nfrag << " Elim(MeV) = "
521  << elim/CLHEP::MeV << G4endl;
522  for(G4int i=0; i<nfrag; ++i) {
524  }
525  G4cout << G4endl;
526 
527 
528  G4cout << "----------------------------------------------------------------"
529  << G4endl;
530  G4cout << "### G4FermiFragmentPoolVI: fragments sorted by A" << G4endl;
531 
532  G4int prec = G4cout.precision(6);
533  G4int ama[maxA];
534  ama[0] = 0;
535  for(G4int A=1; A<maxA; ++A) {
536  G4cout << " # A= " << A << G4endl;
537  size_t am(0);
538  for(size_t j=0; j<list_f[A].size(); ++j) {
539  const G4FermiFragment* f = (list_f[A])[j];
540  G4int a1 = f->GetA();
541  G4int z1 = f->GetZ();
542  size_t nch = (list_c[A])[j]->GetNumberOfChannels();
543  am = std::max(am, nch);
544  G4cout << " ("<<a1<<","<<z1<<"); Eex(MeV)= "
545  << f->GetExcitationEnergy()
546  << " 2S= " << f->GetSpin()
547  << "; Nchannels= " << nch
548  << " MassExcess= " << f->GetTotalEnergy() -
549  (z1*proton_mass_c2 + (a1 - z1)*neutron_mass_c2)
550  << G4endl;
551  for(size_t k=0; k<nch; ++k) {
552  const G4FermiPair* fpair = ((list_c[A])[j]->GetChannels())[k];
553  G4cout << " (" << fpair->GetFragment1()->GetZ()
554  << ", " << fpair->GetFragment1()->GetA()
555  << ", " << fpair->GetFragment1()->GetExcitationEnergy()
556  << ") ("<< fpair->GetFragment2()->GetZ()
557  << ", " << std::setw(3)<< fpair->GetFragment2()->GetA()
558  << ", " << std::setw(8)<< fpair->GetFragment2()->GetExcitationEnergy()
559  << ") prob= " << ((list_c[A])[j]->GetProbabilities())[k]
560  << G4endl;
561  }
562  }
563  ama[A] = am;
564  }
565  G4cout.precision(prec);
566  G4cout << G4endl;
567 
568  G4cout << " Number of fragments per A:" << G4endl;
569  for(G4int j=0; j<maxA; ++j) { G4cout << list_f[j].size() << ", "; }
570  G4cout << G4endl;
571 
572  G4cout << " Max number of channels per A:" << G4endl;
573  for (size_t j=0; j<maxA; ++j) { G4cout << ama[j] << ", "; }
574  G4cout << G4endl;
575 
576  G4cout << " Number of fragment pairs per A:" << G4endl;
577  for(G4int j=0; j<maxA; ++j) { G4cout << list_p[j].size() << ", "; }
578  G4cout << G4endl;
579 
580  G4cout << "----------------------------------------------------------------"
581  << G4endl;
582 
583  G4cout << "### G4FermiFragmentPoolVI: " << funstable.size()
584  << " unphysical fragments" << G4endl;
585  G4cout << " Number of unphysical fragments per A:" << G4endl;
586  for(G4int j=0; j<maxA; ++j) { G4cout << list_g[j].size() << ", "; }
587  G4cout << G4endl;
588 
589  G4cout << " Number of unphysical fragment pairs per A:" << G4endl;
590  for(G4int j=0; j<maxA; ++j) { G4cout << list_u[j].size() << ", "; }
591  G4cout << G4endl;
592 
593  prec = G4cout.precision(6);
594  for(G4int A=1; A<maxA; ++A) {
595  G4cout << " # A= " << A << G4endl;
596  size_t am(0);
597  for(size_t j=0; j<list_g[A].size(); ++j) {
598  const G4FermiFragment* f = (list_g[A])[j];
599  G4int a1 = f->GetA();
600  G4int z1 = f->GetZ();
601  size_t nch = (list_d[A])[j]->GetNumberOfChannels();
602  am = std::max(am, nch);
603  G4cout << "("<<a1<<","<<z1<<"); Eex(MeV)= "
604  << std::setw(8) << f->GetExcitationEnergy()
605  << "; Nchannels= " << nch
606  << " MassExcess= " << f->GetTotalEnergy() -
607  (z1*proton_mass_c2 + (a1 - z1)*neutron_mass_c2)
608  << G4endl;
609  for(size_t k=0; k<nch; ++k) {
610  const G4FermiPair* fpair = ((list_d[A])[j]->GetChannels())[k];
611  G4cout << " (" << fpair->GetFragment1()->GetZ()
612  << ", " << fpair->GetFragment1()->GetA()
613  << ", " << std::setw(8)<< fpair->GetFragment1()->GetExcitationEnergy()
614  << ") ("<< fpair->GetFragment2()->GetZ()
615  << ", " << std::setw(3)<< fpair->GetFragment2()->GetA()
616  << ", " << std::setw(8)<< fpair->GetFragment2()->GetExcitationEnergy()
617  << ") prob= " << ((list_d[A])[j]->GetProbabilities())[k]
618  << G4endl;
619  }
620  }
621  ama[A] = am;
622  G4cout << G4endl;
623  }
624  G4cout.precision(prec);
625  G4cout << G4endl;
626  G4cout << " Max number of unphysical channels per A:" << G4endl;
627  for (size_t j=0; j<maxA; ++j) { G4cout << ama[j] << ", "; }
628  G4cout << G4endl;
629  G4cout << "----------------------------------------------------------------"
630  << G4endl;
631  G4cout << G4endl;
632  G4cout << "### Pairs of stable fragments: " << G4endl;
633 
634  prec = G4cout.precision(6);
635  for(G4int A=2; A<maxA; ++A) {
636  G4cout << " A= " << A<<G4endl;
637  for(size_t j=0; j<list_p[A].size(); ++j) {
638  const G4FermiFragment* f1 = (list_p[A])[j]->GetFragment1();
639  const G4FermiFragment* f2 = (list_p[A])[j]->GetFragment2();
640  G4int a1 = f1->GetA();
641  G4int z1 = f1->GetZ();
642  G4int a2 = f2->GetA();
643  G4int z2 = f2->GetZ();
644  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<") % Eex(MeV)= "
645  << std::setw(8)<< (list_p[A])[j]->GetExcitationEnergy()
646  << " Eex1= " << std::setw(8)<< f1->GetExcitationEnergy()
647  << " Eex2= " << std::setw(8)<< f2->GetExcitationEnergy()
648  << G4endl;
649  }
650  G4cout << G4endl;
651  G4cout <<"----------------------------------------------------------------"
652  << G4endl;
653  }
654  G4cout.precision(prec);
655  G4cout << "### Pairs of stable+unstable fragments: " << G4endl;
656 
657  prec = G4cout.precision(6);
658  for(G4int A=2; A<maxA; ++A) {
659  G4cout << " A= " << A << G4endl;
660  for(size_t j=0; j<list_u[A].size(); ++j) {
661  const G4FermiFragment* f1 = (list_u[A])[j]->GetFragment1();
662  const G4FermiFragment* f2 = (list_u[A])[j]->GetFragment2();
663  G4int a1 = f1->GetA();
664  G4int z1 = f1->GetZ();
665  G4int a2 = f2->GetA();
666  G4int z2 = f2->GetZ();
667  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<") % Eex(MeV)= "
668  << std::setw(8)<< (list_u[A])[j]->GetExcitationEnergy()
669  << " Eex1= " << std::setw(8)<< f1->GetExcitationEnergy()
670  << " Eex2= " << std::setw(8)<< f2->GetExcitationEnergy()
671  << G4endl;
672  }
673  G4cout << G4endl;
674  G4cout << "----------------------------------------------------------------"
675  << G4endl;
676  }
677  G4cout.precision(prec);
678 }
679 
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetFragmentMass(void) const
static const double prec
Definition: RanecuEngine.cc:58
std::vector< const G4FermiFragment * > funstable
G4float MaxLevelEnergy(G4int Z, G4int A) const
G4bool IsPhysical() const
G4DeexPrecoParameters * GetParameters()
G4int GetMinA(G4int Z) const
G4bool IsInPhysPairs(const G4FermiFragment *f1, const G4FermiFragment *f2) const
const G4FermiFragment * GetFragment1() const
Definition: G4FermiPair.hh:97
G4bool IsStable() const
#define G4endl
Definition: G4ios.hh:61
static G4NuclearLevelData * GetInstance()
float G4float
Definition: G4Types.hh:77
G4FermiDecayProbability theDecay
G4double ComputeProbability(G4int Z, G4int A, G4int spin, G4double TotalE, const G4FermiFragment *f1, const G4FermiFragment *f2) const
Float_t f2
G4double GetMinExcitation() const
G4bool IsInUnphysPairs(const G4FermiFragment *f1, const G4FermiFragment *f2) const
static const G4int maxZ
G4int GetA(void) const
static constexpr double neutron_mass_c2
static constexpr double proton_mass_c2
const std::vector< const G4FermiPair * > & GetChannels() const
std::vector< G4FermiChannels * > list_c[maxA]
Float_t Z
void DumpFragment(const G4FermiFragment *) const
Float_t f3
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4double GetMaxLifeTime() const
G4bool IsPhysical(G4int Z, G4int A) const
static constexpr double MeV
double A(double temperature)
G4double GetCoulombBarrier(G4int Ares, G4int Zres, G4double Eex) const
static const G4int maxA
static G4double GetNuclearMass(const G4double A, const G4double Z)
Float_t f1
G4double LifeTime(size_t i) const
G4double GetTotalEnergy(void) const
G4int GetMaxA(G4int Z) const
std::vector< const G4FermiPair * > list_u[maxA]
Double_t Z2
size_t GetNumberOfChannels() const
Double_t Z1
int G4int
Definition: G4Types.hh:78
std::vector< const G4FermiFragment * > list_f[maxA]
G4int SpinTwo(size_t i) const
G4int GetSpin(void) const
G4bool IsApplicable(G4int Z, G4int A, G4double etot) const
std::vector< G4FermiChannels * > list_d[maxA]
std::vector< const G4FermiFragment * > fragment_pool
size_t NumberOfTransitions() const
const G4FermiChannels * ClosestChannels(G4int Z, G4int A, G4double mass) const
G4GLOB_DLL std::ostream G4cout
G4int GetZ(void) const
double minE
Definition: plot_hist.C:7
G4double GetExcitationEnergy(void) const
std::vector< G4double > & GetProbabilities()
std::vector< const G4FermiPair * > list_p[maxA]
std::vector< const G4FermiFragment * > list_g[maxA]
G4double LevelEnergy(size_t i) const
static constexpr double GeV
Definition: G4SIunits.hh:217
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool IsInThePool(G4int Z, G4int A, G4double exc) const
const G4FermiFragment * GetFragment2() const
Definition: G4FermiPair.hh:102