ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FermiFragmentsPoolVI.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 // FermiBreakUp de-excitation model
28 // by V. Ivanchenko (July 2016)
29 //
30 
31 #include "G4FermiFragmentsPoolVI.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4NuclearLevelData.hh"
35 #include "G4LevelManager.hh"
36 #include "G4DeexPrecoParameters.hh"
37 #include <iomanip>
38 
39 G4FermiFragmentsPoolVI::G4FermiFragmentsPoolVI()
40 {
41  // G4cout << "### G4FermiFragmentsPoolVI is constructed" << G4endl;
42  G4DeexPrecoParameters* param =
43  G4NuclearLevelData::GetInstance()->GetParameters();
44  tolerance = param->GetMinExcitation();
45  timelim = (G4float)param->GetMaxLifeTime();
46 
47  elim = param->GetFBUEnergyLimit();
48  elimf= (G4float)elim;
49  /*
50  G4cout << "G4FermiFragmentsPoolVI: tolerance= " << tolerance
51  << " timelim= " << timelim << " elim= " << elim << G4endl;
52  */
53  fragment_pool.reserve(991);
54  Initialise();
55 }
56 
57 G4FermiFragmentsPoolVI::~G4FermiFragmentsPoolVI()
58 {
59  for(G4int i=0; i<maxA; ++i) {
60  for(auto & ptr : list_p[i]) { delete ptr; ptr = nullptr; }
61  for(auto & ptr : list_c[i]) { delete ptr; ptr = nullptr; }
62  }
63  for(auto & ptr : fragment_pool) { delete ptr; ptr = nullptr; }
64 }
65 
66 const G4FermiChannels*
67 G4FermiFragmentsPoolVI::ClosestChannels(G4int Z, G4int A, G4double e) const
68 {
69  const G4FermiChannels* res = nullptr;
70  G4double demax = 1.e+9;
71 
72  // stable channels
73  for(size_t j=0; j<(list_c[A]).size(); ++j) {
74  const G4FermiFragment* frag = (list_f[A])[j];
75  if(frag->GetZ() != Z) { continue; }
76  G4double de = e - frag->GetTotalEnergy();
77  //G4cout << " Stab check " << j << " channel de= " << de
78  // << " tol= " << tolerance << G4endl;
79  // an excitation coincide with a level
80  if(std::abs(de) <= tolerance) {
81  res = (list_c[A])[j];
82  break;
83  } else {
84  // closest level selected
85  de += tolerance;
86  if(de >= 0.0 && de <= demax) {
87  res = (list_c[A])[j];
88  demax = de;
89  }
90  //G4cout << " Stab chan: " << j << " N= "
91  //<< res->GetNumberOfChannels() << G4endl;
92  }
93  }
94  return res;
95 }
96 
97 G4bool G4FermiFragmentsPoolVI::IsPhysical(G4int Z, G4int A) const
98 {
99  for(auto const& ptr : list_f[A]) {
100  if(ptr->GetZ() == Z) { return true; }
101  }
102  return false;
103 }
104 
105 G4bool G4FermiFragmentsPoolVI::IsInThePool(G4int Z, G4int A,
106  G4double exc) const
107 {
108  for(auto const& fr : fragment_pool) {
109  if(fr->GetZ() == Z && fr->GetA() == A &&
110  std::abs(exc - fr->GetExcitationEnergy()) < tolerance)
111  { return true; }
112  }
113  return false;
114 }
115 
116 G4bool
117 G4FermiFragmentsPoolVI::HasChannels(G4int Z, G4int A, G4double exc) const
118 {
119  // stable fragment
120  for(size_t j=0; j<(list_f[A]).size(); ++j) {
121  const G4FermiFragment* frag = (list_f[A])[j];
122  if(frag->GetZ() == Z) {
123  if(exc > frag->GetExcitationEnergy() &&
124  (list_c[A])[j]->GetNumberOfChannels() > 0) { return true; }
125  }
126  }
127  return false;
128 }
129 
130 G4bool G4FermiFragmentsPoolVI::IsInPhysPairs(
131  const G4FermiFragment* f1, const G4FermiFragment* f2) const
132 {
133  const G4int A = f1->GetA() + f2->GetA();
134  for(auto const& ptr : list_p[A]) {
135  if(f1 == ptr->GetFragment1() && f2 == ptr->GetFragment2()) {
136  return true;
137  }
138  }
139  return false;
140 }
141 
142 void G4FermiFragmentsPoolVI::Initialise()
143 {
144  //G4cout << "G4FermiFragmentsPoolVI::Initialise main loop @@@@@@" << G4endl;
145 
146  // stable particles
147  fragment_pool.push_back(new G4FermiFragment(1, 0, 1, 0.0));
148  fragment_pool.push_back(new G4FermiFragment(1, 1, 1, 0.0));
149  fragment_pool.push_back(new G4FermiFragment(2, 1, 2, 0.0));
150  fragment_pool.push_back(new G4FermiFragment(3, 1, 1, 0.0));
151  fragment_pool.push_back(new G4FermiFragment(3, 2, 1, 0.0));
152  fragment_pool.push_back(new G4FermiFragment(4, 2, 0, 0.0));
153  fragment_pool.push_back(new G4FermiFragment(5, 2, 3, 0.0));
154  fragment_pool.push_back(new G4FermiFragment(5, 3, 3, 0.0));
155 
156  // use level data and construct the pool
157  G4NuclearLevelData* ndata = G4NuclearLevelData::GetInstance();
158  for(G4int Z=1; Z<maxZ; ++Z) {
159  G4int Amin = ndata->GetMinA(Z);
160  G4int Amax = std::min(maxA, ndata->GetMaxA(Z)+1);
161  for(G4int A=Amin; A<Amax; ++A) {
162  const G4LevelManager* man = ndata->GetLevelManager(Z, A);
163  if(man) {
164  size_t nn = man->NumberOfTransitions();
165  // very unstable state
166  if(ndata->MaxLevelEnergy(Z, A) == 0.0f && man->LifeTime(0) == 0.0f) {
167  continue;
168  }
169  for(size_t i=0; i<=nn; ++i) {
170  G4float exc = man->LevelEnergy(i);
171  /*
172  G4cout << "Z= " << Z << " A= " << A << " Eex= " << exc
173  << " elimf= " << elimf << " toler= " << tolerance
174  << " time= " << man->LifeTime(i) << " i= " << i << G4endl;
175  */
176  // only levels below limit are consided
177  if(exc >= elimf) { continue; }
178  G4double excd = (G4double)exc;
179  // only new are considered
180  if(IsInThePool(Z, A, excd)) { continue; }
181  fragment_pool.push_back(new G4FermiFragment(A,Z,man->SpinTwo(i),excd));
182  }
183  }
184  }
185  }
186  G4int nfrag = fragment_pool.size();
187  // prepare structures per A for normal fragments
188  const size_t lfmax[maxA] = {
189  0, 2, 1, 2, 1, 2, 8, 19, 28, 56, 70, 104, 74, 109, 143, 212, 160};
190  for(G4int A=1; A<maxA; ++A) {
191  list_f[A].reserve(lfmax[A]);
192  list_c[A].reserve(lfmax[A]);
193  }
194  const size_t lfch[maxA] = {
195  0, 0, 0, 0, 0, 1, 4, 8, 6, 13, 27, 40, 29, 21, 31, 32, 30};
196 
197  for(auto const& f : fragment_pool) {
198  G4int A = f->GetA();
199  G4double exc = f->GetExcitationEnergy();
200  list_f[A].push_back(f);
201  list_c[A].push_back(new G4FermiChannels(lfch[A], exc, f->GetTotalEnergy()));
202  }
203  /*
204  G4cout << "Defined fragments @@@@@@"
205  << " PhysicalFrag= " << nfrag
206  << " UnphysicalFrag= " << funstable.size() << G4endl;
207  */
208  // list of fragment pairs ordered by A
209  for(G4int i=0; i<nfrag; ++i) {
210  const G4FermiFragment* f1 = fragment_pool[i];
211  G4int Z1 = f1->GetZ();
212  G4int A1 = f1->GetA();
213  G4double e1 = f1->GetTotalEnergy();
214  for(G4int j=0; j<nfrag; ++j) {
215  const G4FermiFragment* f2 = fragment_pool[j];
216  G4int Z2 = f2->GetZ();
217  G4int A2 = f2->GetA();
218  if(A2 < A1 || (A2 == A1 && Z2 < Z1)) { continue; }
219  G4int Z = Z1 + Z2;
220  G4int A = A1 + A2;
221 
222  if(Z >= maxZ || A >= maxA || IsInPhysPairs(f1, f2)) { continue; }
223 
224  G4double e2 = f2->GetTotalEnergy();
225  G4double minE = e1 + e2;
226  G4double exc = 0.0;
227  if(IsPhysical(Z, A)) {
228  minE += f1->GetCoulombBarrier(A2, Z2, 0.0);
229  exc = minE - G4NucleiProperties::GetNuclearMass(A, Z);
230  }
231  /*
232  G4cout << "Z= " << Z << " A= " << A
233  << " Z1= " << Z1 << " A1= " << A1
234  << " Z2= " << Z2 << " A2= " << A2 << " Eex= " << exc
235  << " Qb= " << f1->GetCoulombBarrier(A2, Z2, 0.0)
236  << " " << e1
237  << " " << e2
238  << " " << G4NucleiProperties::GetNuclearMass(A, Z)
239  << G4endl;
240  */
241  // ignore very excited case
242  if(exc >= elim) { continue; }
243  G4FermiPair* fpair = nullptr;
244  G4int kmax = list_f[A].size();
245  for(G4int k=0; k<kmax; ++k) {
246  const G4FermiFragment* f3 = (list_f[A])[k];
247  if(Z == f3->GetZ() &&
248  f3->GetTotalEnergy() - minE + tolerance >= 0.0) {
249  if(!fpair) {
250  fpair = new G4FermiPair(f1, f2);
251  list_p[A].push_back(fpair);
252  }
253  (list_c[A])[k]->AddChannel(fpair);
254  }
255  }
256  }
257  }
258  // compute static probabilities
259  for(G4int A=1; A<maxA; ++A) {
260  for(size_t j=0; j<list_c[A].size(); ++j) {
261  G4FermiChannels* ch = (list_c[A])[j];
262  const G4FermiFragment* frag = (list_f[A])[j];
263  size_t nch = ch->GetNumberOfChannels();
264  if(1 < nch) {
265  std::vector<G4double>& prob = ch->GetProbabilities();
266  const std::vector<const G4FermiPair*>& pairs = ch->GetChannels();
267  G4double ptot = 0.0;
268  for(size_t i=0; i<nch; ++i) {
269  ptot += theDecay.ComputeProbability(frag->GetZ(), frag->GetA(),
270  frag->GetSpin(),
271  frag->GetTotalEnergy(),
272  pairs[i]->GetFragment1(),
273  pairs[i]->GetFragment2());
274  prob[i] = ptot;
275  }
276  if(0.0 == ptot) {
277  prob[0] = 1.0;
278  } else {
279  ptot = 1./ptot;
280  for(size_t i=0; i<nch-1; ++i) { prob[i] *= ptot; }
281  prob[nch-1] = 1.0;
282  }
283  }
284  }
285  }
286 }
287 
288 void G4FermiFragmentsPoolVI::DumpFragment(const G4FermiFragment* f) const
289 {
290  if(f) {
291  G4int prec = G4cout.precision(6);
292  G4cout << " Z= " << f->GetZ() << " A= " << std::setw(2) << f->GetA()
293  << " Mass(GeV)= " << std::setw(8) << f->GetFragmentMass()/GeV
294  << " Eexc(MeV)= " << std::setw(7) << f->GetExcitationEnergy()
295  << " 2s= " << f->GetSpin() << " IsStable: "
296  << HasChannels(f->GetZ(), f->GetA(), f->GetExcitationEnergy())
297  << G4endl;
298  G4cout.precision(prec);
299  }
300 }
301 
302 void G4FermiFragmentsPoolVI::Dump() const
303 {
304  G4cout <<"----------------------------------------------------------------"
305  <<G4endl;
306  G4cout << "##### List of Fragments in the Fermi Fragment Pool #####"
307  << G4endl;
308  G4int nfrag = fragment_pool.size();
309  G4cout << " For stable " << nfrag << " Elim(MeV) = "
310  << elim/CLHEP::MeV << G4endl;
311  for(G4int i=0; i<nfrag; ++i) {
312  DumpFragment(fragment_pool[i]);
313  }
314  G4cout << G4endl;
315 
316 
317  G4cout << "----------------------------------------------------------------"
318  << G4endl;
319  G4cout << "### G4FermiFragmentPoolVI: fragments sorted by A" << G4endl;
320 
321  G4int prec = G4cout.precision(6);
322  G4int ama[maxA];
323  ama[0] = 0;
324  for(G4int A=1; A<maxA; ++A) {
325  G4cout << " # A= " << A << G4endl;
326  size_t am(0);
327  for(size_t j=0; j<list_f[A].size(); ++j) {
328  const G4FermiFragment* f = (list_f[A])[j];
329  G4int a1 = f->GetA();
330  G4int z1 = f->GetZ();
331  size_t nch = (list_c[A])[j]->GetNumberOfChannels();
332  am = std::max(am, nch);
333  G4cout << " ("<<a1<<","<<z1<<"); Eex(MeV)= "
334  << f->GetExcitationEnergy()
335  << " 2S= " << f->GetSpin()
336  << "; Nchannels= " << nch
337  << " MassExcess= " << f->GetTotalEnergy() -
338  (z1*proton_mass_c2 + (a1 - z1)*neutron_mass_c2)
339  << G4endl;
340  for(size_t k=0; k<nch; ++k) {
341  const G4FermiPair* fpair = ((list_c[A])[j]->GetChannels())[k];
342  G4cout << " (" << fpair->GetFragment1()->GetZ()
343  << ", " << fpair->GetFragment1()->GetA()
344  << ", " << fpair->GetFragment1()->GetExcitationEnergy()
345  << ") ("<< fpair->GetFragment2()->GetZ()
346  << ", " << std::setw(3)<< fpair->GetFragment2()->GetA()
347  << ", " << std::setw(8)<< fpair->GetFragment2()->GetExcitationEnergy()
348  << ") prob= " << ((list_c[A])[j]->GetProbabilities())[k]
349  << G4endl;
350  }
351  }
352  ama[A] = am;
353  }
354  G4cout.precision(prec);
355  G4cout << G4endl;
356 
357  G4cout << " Number of fragments per A:" << G4endl;
358  for(G4int j=0; j<maxA; ++j) { G4cout << list_f[j].size() << ", "; }
359  G4cout << G4endl;
360 
361  G4cout << " Max number of channels per A:" << G4endl;
362  for (size_t j=0; j<maxA; ++j) { G4cout << ama[j] << ", "; }
363  G4cout << G4endl;
364 
365  G4cout << " Number of fragment pairs per A:" << G4endl;
366  for(G4int j=0; j<maxA; ++j) { G4cout << list_p[j].size() << ", "; }
367  G4cout << G4endl;
368 
369  G4cout << "----------------------------------------------------------------"
370  << G4endl;
371  G4cout << "### Pairs of stable fragments: " << G4endl;
372 
373  prec = G4cout.precision(6);
374  for(G4int A=2; A<maxA; ++A) {
375  G4cout << " A= " << A<<G4endl;
376  for(size_t j=0; j<list_p[A].size(); ++j) {
377  const G4FermiFragment* f1 = (list_p[A])[j]->GetFragment1();
378  const G4FermiFragment* f2 = (list_p[A])[j]->GetFragment2();
379  G4int a1 = f1->GetA();
380  G4int z1 = f1->GetZ();
381  G4int a2 = f2->GetA();
382  G4int z2 = f2->GetZ();
383  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<") % Eex(MeV)= "
384  << std::setw(8)<< (list_p[A])[j]->GetExcitationEnergy()
385  << " Eex1= " << std::setw(8)<< f1->GetExcitationEnergy()
386  << " Eex2= " << std::setw(8)<< f2->GetExcitationEnergy()
387  << G4endl;
388  }
389  G4cout << G4endl;
390  G4cout <<"----------------------------------------------------------------"
391  << G4endl;
392  }
393  G4cout.precision(prec);
394 }
395