Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FermiFragmentsPool.cc
Go to the documentation of this file.
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: G4FermiFragmentsPool.cc,v 1.5 2006-06-29 20:13:13 gunter Exp $
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara
30 //
31 // Modifications:
32 // J.M.Quesada, July 2009, bug fixed in excitation energies:
33 // ALL of them are in MeV instead of keV (as they were expressed previously)
34 // source: http://www.nndc.bnl.gov/chart
35 // Unknown excitation energies in He5 and Li5 have been suppressed
36 // Long lived levels (half-lives of the order ps-fs have been included)
37 //
38 // J. M. Quesada, April 2010: excitation energies according to tabulated values
39 // in PhotonEvaporatoion2.0. Fake photons eliminated.
40 //
41 // 01.04.2011 General cleanup by V.Ivanchenko - more clean usage of static
42 //
43 // 04.05.2011 J. M. Quesada: added detailed printout for testing
44 
45 #include "G4FermiFragmentsPool.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4StableFermiFragment.hh"
49 #include "G4B9FermiFragment.hh"
50 #include "G4Be8FermiFragment.hh"
51 #include "G4He5FermiFragment.hh"
52 #include "G4Li5FermiFragment.hh"
53 
54 G4FermiFragmentsPool* G4FermiFragmentsPool::theInstance = nullptr;
55 
56 G4FermiFragmentsPool* G4FermiFragmentsPool::Instance()
57 {
58  if(0 == theInstance) {
59  static G4FermiFragmentsPool instance;
60  theInstance = &instance;
61  }
62  return theInstance;
63 }
64 
65 G4FermiFragmentsPool::G4FermiFragmentsPool()
66 {
67  maxZ = 9;
68  maxA = 17;
69  verbose = 0;
70  Initialise();
71 }
72 
73 G4FermiFragmentsPool::~G4FermiFragmentsPool()
74 {
75  size_t nn;
76  for(size_t i=0; i<17; ++i) {
77  nn = list1[i].size();
78  if(0 < nn) { for(size_t j=0; j<nn; ++j) { delete (list1[i])[j]; }}
79  nn = list2[i].size();
80  if(0 < nn) { for(size_t j=0; j<nn; ++j) { delete (list2[i])[j]; }}
81  nn = list3[i].size();
82  if(0 < nn) { for(size_t j=0; j<nn; ++j) { delete (list3[i])[j]; }}
83  nn = list4[i].size();
84  if(0 < nn) { for(size_t j=0; j<nn; ++j) { delete (list4[i])[j]; }}
85  }
86  nn = fragment_pool.size();
87  if(0 < nn) { for(size_t j=0; j<nn; ++j) { delete fragment_pool[j]; }}
88 }
89 
90 G4bool G4FermiFragmentsPool::IsAvailable(G4int Z, G4int A) const
91 {
92  G4bool res = true;
93  if (2 == Z && 5 == A) { res = false; }
94  else if(3 == Z && 5 == A) { res = false; }
95  else if(4 == Z && 8 == A) { res = false; }
96  else if(5 == Z && 9 == A) { res = false; }
97  return res;
98 }
99 
100 G4int G4FermiFragmentsPool::GetMaxZ() const
101 {
102  return maxZ;
103 }
104 
105 G4int G4FermiFragmentsPool::GetMaxA() const
106 {
107  return maxA;
108 }
109 
110 const G4FermiPhaseSpaceDecay*
111 G4FermiFragmentsPool::GetFermiPhaseSpaceDecay() const
112 {
113  return &thePhaseSpace;
114 }
115 
116 void G4FermiFragmentsPool::Initialise()
117 {
118  // JMQ 30/06/09 unknown levels have been supressed
119  // JMQ 01/07/09 corrected excitation energies for 64-66, according to
120  // http://www.nndc.bnl.gov/chart
121  // JMQ 19/04/10 new level, fragment numbering shifted accordingly from here onwards
122  // A Z Pol ExcitE
123  fragment_pool.push_back(new G4StableFermiFragment( 1, 0, 2, 0.00*MeV ));
124  fragment_pool.push_back(new G4StableFermiFragment( 1, 1, 2, 0.00*MeV ));
125  fragment_pool.push_back(new G4StableFermiFragment( 2, 1, 3, 0.00*MeV ));
126  fragment_pool.push_back(new G4StableFermiFragment( 3, 1, 2, 0.00*MeV ));
127  fragment_pool.push_back(new G4StableFermiFragment( 3, 2, 2, 0.00*MeV ));
128  fragment_pool.push_back(new G4StableFermiFragment( 4, 2, 1, 0.00*MeV ));
129  fragment_pool.push_back(new G4He5FermiFragment ( 5, 2, 4, 0.00*MeV ));
130  fragment_pool.push_back(new G4Li5FermiFragment ( 5, 3, 4, 0.00*MeV ));
131  fragment_pool.push_back(new G4StableFermiFragment( 6, 2, 1, 0.00*MeV ));
132  fragment_pool.push_back(new G4StableFermiFragment( 6, 3, 3, 0.00*MeV ));
133  fragment_pool.push_back(new G4StableFermiFragment( 6, 3, 1, 3.562880*MeV ));
134  fragment_pool.push_back(new G4StableFermiFragment( 7, 3, 4, 0.00*MeV ));
135  fragment_pool.push_back(new G4StableFermiFragment( 7, 3, 2, 0.4776120*MeV ));
136  fragment_pool.push_back(new G4StableFermiFragment( 7, 4, 4, 0.00*MeV ));
137  fragment_pool.push_back(new G4StableFermiFragment( 7, 4, 2, 0.4290800*MeV ));
138  fragment_pool.push_back(new G4StableFermiFragment( 8, 3, 5, 0.00*MeV ));
139  fragment_pool.push_back(new G4StableFermiFragment( 8, 3, 3, 0.9808000*MeV ));
140  fragment_pool.push_back(new G4Be8FermiFragment ( 8, 4, 1, 0.00*MeV ));
141  fragment_pool.push_back(new G4StableFermiFragment( 9, 4, 4, 0.00*MeV ));
142  fragment_pool.push_back(new G4B9FermiFragment ( 9, 5, 4, 0.00*MeV ));
143  fragment_pool.push_back(new G4StableFermiFragment( 10, 4, 1, 0.00*MeV ));
144  fragment_pool.push_back(new G4StableFermiFragment( 10, 4, 5, 3.368030*MeV ));
145  fragment_pool.push_back(new G4StableFermiFragment( 10, 4, 8, 5.958390*MeV ));
146  fragment_pool.push_back(new G4StableFermiFragment( 10, 4, 1, 6.179300*MeV ));
147  fragment_pool.push_back(new G4StableFermiFragment( 10, 4, 5, 6.263300*MeV ));
148  fragment_pool.push_back(new G4StableFermiFragment( 10, 5, 7, 0.00*MeV ));
149  fragment_pool.push_back(new G4StableFermiFragment( 10, 5, 3, 0.7183800*MeV ));
150  fragment_pool.push_back(new G4StableFermiFragment( 10, 5, 1, 1.740050*MeV ));
151  fragment_pool.push_back(new G4StableFermiFragment( 10, 5, 3, 2.154270*MeV ));
152  fragment_pool.push_back(new G4StableFermiFragment( 10, 5, 5, 3.587130*MeV ));
153  fragment_pool.push_back(new G4StableFermiFragment( 10, 6, 3, 0.00*MeV ));
154  fragment_pool.push_back(new G4StableFermiFragment( 10, 6, 5, 3.353700*MeV ));
155  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 4, 0.00*MeV ));
156  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 2, 2.124693*MeV ));
157  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 6, 4.444980*MeV ));
158  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 4, 5.020300*MeV ));
159  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 8, 6.741850*MeV ));
160  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 2, 6.791800*MeV ));
161  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 6, 7.285510*MeV ));
162  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 4, 7.977840*MeV ));
163  fragment_pool.push_back(new G4StableFermiFragment( 11, 5, 6, 8.560100*MeV ));
164  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 4, 0.00*MeV ));
165  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 2, 2.00*MeV ));
166  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 6, 4.318800*MeV ));
167  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 4, 4.804200*MeV ));
168  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 2, 6.339200*MeV ));
169  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 8, 6.478200*MeV ));
170  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 6, 6.904800*MeV ));
171  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 4, 7.499700*MeV ));
172  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 4, 8.104500*MeV ));
173  fragment_pool.push_back(new G4StableFermiFragment( 11, 6, 6, 8.420000*MeV ));
174  fragment_pool.push_back(new G4StableFermiFragment( 12, 5, 3, 0.00*MeV ));
175  fragment_pool.push_back(new G4StableFermiFragment( 12, 5, 5, 0.9531400*MeV ));
176  fragment_pool.push_back(new G4StableFermiFragment( 12, 5, 5, 1.673650*MeV ));
177  fragment_pool.push_back(new G4StableFermiFragment( 12, 5, 3, 2.620800*MeV ));
178  fragment_pool.push_back(new G4StableFermiFragment( 12, 6, 1, 0.00*MeV ));
179  fragment_pool.push_back(new G4StableFermiFragment( 12, 6, 5, 4.438910*MeV ));
180  fragment_pool.push_back(new G4StableFermiFragment( 13, 6, 2, 0.00*MeV ));
181  fragment_pool.push_back(new G4StableFermiFragment( 13, 6, 2, 3.089443*MeV ));
182  fragment_pool.push_back(new G4StableFermiFragment( 13, 6, 4, 3.684507*MeV ));
183  fragment_pool.push_back(new G4StableFermiFragment( 13, 6, 6, 3.853807*MeV ));
184  fragment_pool.push_back(new G4StableFermiFragment( 13, 7, 2, 0.00*MeV ));
185  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 1, 0.00*MeV ));
186  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 3, 6.093800*MeV ));
187  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 1, 6.589400*MeV ));
188  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 7, 6.728200*MeV ));
189  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 1, 6.902600*MeV ));
190  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 5, 7.012000*MeV ));
191  fragment_pool.push_back(new G4StableFermiFragment( 14, 6, 5, 7.341000*MeV ));
192  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 3, 0.00*MeV ));
193  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 1, 2.312798*MeV ));
194  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 3, 3.948100*MeV ));
195  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 1, 4.915100*MeV ));
196  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 5, 5.105890*MeV ));
197  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 3, 5.691440*MeV ));
198  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 7, 5.834250*MeV ));
199  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 3, 6.203500*MeV ));
200  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 7, 6.446170*MeV ));
201  fragment_pool.push_back(new G4StableFermiFragment( 14, 7, 5, 7.029120*MeV ));
202  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 2, 0.00*MeV ));
203  // JMQ 010709 two very close levels instead of only one, with their own spins
204  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 6, 5.270155*MeV ));
205  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 2, 5.298822*MeV ));
206  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 6.323780*MeV ));
207  //JMQ 010709 new level and corrected energy and spins
208  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 6, 7.155050*MeV ));
209  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 7.300830*MeV ));
210  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 8, 7.567100*MeV ));
211  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 2, 8.312620*MeV ));
212  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 8.571400*MeV ));
213  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 2, 9.049710*MeV ));
214  //JMQ 010709 new levels for N15
215  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 9.151900*MeV ));
216  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 6, 9.154900*MeV ));
217  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 2, 9.222100*MeV ));
218  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 6, 9.760000*MeV ));
219  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 8, 9.829000*MeV ));
220  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 9.925000*MeV ));
221  fragment_pool.push_back(new G4StableFermiFragment( 15, 7, 4, 10.06600*MeV ));
222  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 2, 0.00*MeV ));
223  //JMQ 010709 new level and spins
224  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 2, 5.183000*MeV ));
225  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 6, 5.240900*MeV ));
226  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 4, 6.176300*MeV ));
227  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 4, 6.793100*MeV ));
228  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 6, 6.859400*MeV ));
229  fragment_pool.push_back(new G4StableFermiFragment( 15, 8, 8, 7.275900*MeV ));
230  fragment_pool.push_back(new G4StableFermiFragment( 16, 7, 5, 0.00*MeV ));
231  fragment_pool.push_back(new G4StableFermiFragment( 16, 7, 1, 0.1204200*MeV ));
232  fragment_pool.push_back(new G4StableFermiFragment( 16, 7, 7, 0.2982200*MeV ));
233  fragment_pool.push_back(new G4StableFermiFragment( 16, 7, 3, 0.3972700*MeV ));
234  //JMQ 010709 some energies and spins have been changed
235  fragment_pool.push_back(new G4StableFermiFragment( 16, 8, 1, 0.00*MeV ));
236  fragment_pool.push_back(new G4StableFermiFragment( 16, 8, 1, 6.049400*MeV ));
237  fragment_pool.push_back(new G4StableFermiFragment( 16, 8, 7, 6.129890*MeV ));
238  fragment_pool.push_back(new G4StableFermiFragment( 16, 8, 5, 6.917100*MeV ));
239  //JMQ 180510 fixed fragment 111
240  fragment_pool.push_back(new G4StableFermiFragment( 16, 8, 3, 7.116850*MeV ));
241 
242  G4int nfrag = fragment_pool.size();
243  std::vector<const G4VFermiFragment*> newvec;
244  newvec.reserve(4);
245 
246  // list of fragments ordered by A
247  newvec.resize(1);
248  for(G4int i=0; i<nfrag; ++i) {
249  newvec[0] = fragment_pool[i];
250  const G4FermiConfiguration* conf = new G4FermiConfiguration(newvec);
251  G4int A = fragment_pool[i]->GetA();
252  list1[A].push_back(conf);
253  }
254  if(verbose > 0) {
255  G4cout << "### G4FermiFragmentPool: " << nfrag
256  << " fragments" << G4endl;
257  for(G4int A=1; A<maxA; ++A) {
258  G4cout << " A= " << A << " : Z= ";
259  for(size_t j=0; j<list1[A].size(); ++j) {
260  G4cout << (list1[A])[j]->GetZ() << " ";
261  }
262  G4cout << G4endl;
263  }
264  }
265 
266  // list of fragment pairs ordered by A
267  G4int counter = 0;
268  G4int tot = 0;
269  newvec.resize(2);
270  for(G4int i=0; i<nfrag; ++i) {
271  G4int Z1 = fragment_pool[i]->GetZ();
272  G4int A1 = fragment_pool[i]->GetA();
273  for(G4int j=0; j<nfrag; ++j) {
274  G4int Z2 = fragment_pool[j]->GetZ();
275  G4int A2 = fragment_pool[j]->GetA();
276  G4int Z = Z1 + Z2;
277  G4int A = A1 + A2;
278  if(Z < maxZ && A < maxA && IsAvailable(Z, A)) {
279  newvec[0] = fragment_pool[i];
280  newvec[1] = fragment_pool[j];
281  if(!IsExist(Z, A, newvec)) {
282  const G4FermiConfiguration* conf = new G4FermiConfiguration(newvec);
283  list2[A].push_back(conf);
284  ++counter;
285  }
286  }
287  }
288  }
289  if(verbose > 0) {
290  G4cout << G4endl;
291  G4cout << "### Pairs of fragments: " << counter << G4endl;
292  for(G4int A=2; A<maxA; ++A) {
293  G4cout << " A= " << A<<G4endl;
294  for(size_t j=0; j<list2[A].size(); ++j) {
295  const std::vector<const G4VFermiFragment*>* vec
296  = (list2[A])[j]->GetFragmentList();
297  G4int a1=(*vec)[0]->GetA();
298  G4int z1=(*vec)[0]->GetZ();
299  G4int a2=(*vec)[1]->GetA();
300  G4int z2=(*vec)[1]->GetZ();
301  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<") % ";
302  }
303  G4cout<<G4endl;
304  G4cout<<"----------------------------------------------------------------"
305  << G4endl;
306  }
307  }
308 
309  // list of fragment triples ordered by A
310  tot += counter;
311  counter = 0;
312  newvec.resize(3);
313  for(G4int A1=2; A1<maxA; ++A1) {
314  size_t nz = list2[A1].size();
315  for(size_t idx=0; idx<nz; ++idx) {
316  const G4FermiConfiguration* conf2 = (list2[A1])[idx];
317  G4int Z1 = conf2->GetZ();
318  const std::vector<const G4VFermiFragment*>* vec2 =
319  conf2->GetFragmentList();
320  for(G4int j=0; j<nfrag; ++j) {
321  G4int Z2 = fragment_pool[j]->GetZ();
322  G4int A2 = fragment_pool[j]->GetA();
323  G4int Z = Z1 + Z2;
324  G4int A = A1 + A2;
325  if(Z < maxZ && A < maxA && IsAvailable(Z, A)) {
326  newvec[0] = (*vec2)[0];
327  newvec[1] = (*vec2)[1];
328  newvec[2] = fragment_pool[j];
329  if(!IsExist(Z, A, newvec)) {
330  const G4FermiConfiguration* conf3 = new G4FermiConfiguration(newvec);
331  list3[A].push_back(conf3);
332  ++counter;
333  }
334  }
335  }
336  }
337  }
338  if(verbose > 0) {
339  G4cout << G4endl;
340  G4cout << "### Triples of fragments: " << counter << G4endl;
341  for(G4int A=3; A<maxA; ++A) {
342  G4cout << " A= " << A<<G4endl;
343  for(size_t j=0; j<list3[A].size(); ++j) {
344  const std::vector<const G4VFermiFragment*>* vec
345  = (list3[A])[j]->GetFragmentList();
346  G4int a1=(*vec)[0]->GetA();
347  G4int z1=(*vec)[0]->GetZ();
348  G4int a2=(*vec)[1]->GetA();
349  G4int z2=(*vec)[1]->GetZ();
350  G4int a3=(*vec)[2]->GetA();
351  G4int z3=(*vec)[2]->GetZ();
352  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<")("<<a3<<","<<z3<<") % ";
353  }
354  G4cout<<G4endl;
355  G4cout<<"----------------------------------------------------------------"
356  << G4endl;
357  }
358  }
359 
360  // list of fragment quartets (3 + 1) ordered by A
361  tot += counter;
362  counter = 0;
363  newvec.resize(4);
364  for(G4int A1=3; A1<maxA; ++A1) {
365  size_t nz = list3[A1].size();
366  for(size_t idx=0; idx<nz; ++idx) {
367  const G4FermiConfiguration* conf3 = (list3[A1])[idx];
368  G4int Z1 = conf3->GetZ();
369  const std::vector<const G4VFermiFragment*>* vec3 =
370  conf3->GetFragmentList();
371  for(G4int j=0; j<nfrag; ++j) {
372  G4int Z2 = fragment_pool[j]->GetZ();
373  G4int A2 = fragment_pool[j]->GetA();
374  G4int Z = Z1 + Z2;
375  G4int A = A1 + A2;
376  if(Z < maxZ && A < maxA && IsAvailable(Z, A)) {
377  newvec[0] = (*vec3)[0];
378  newvec[1] = (*vec3)[1];
379  newvec[2] = (*vec3)[2];
380  newvec[3] = fragment_pool[j];
381  if(!IsExist(Z, A, newvec)) {
382  const G4FermiConfiguration* conf4 = new G4FermiConfiguration(newvec);
383  list4[A].push_back(conf4);
384  ++counter;
385  }
386  }
387  }
388  }
389  }
390  // list of fragment quartets (2 + 2) ordered by A
391  for(G4int A1=2; A1<maxA; ++A1) {
392  size_t nz1 = list2[A1].size();
393  for(size_t id1=0; id1<nz1; ++id1) {
394  const G4FermiConfiguration* conf1 = (list2[A1])[id1];
395  G4int Z1 = conf1->GetZ();
396  const std::vector<const G4VFermiFragment*>* vec1 =
397  conf1->GetFragmentList();
398  for(G4int A2=2; A2<maxA; ++A2) {
399  size_t nz2 = list2[A2].size();
400  for(size_t id2=0; id2<nz2; ++id2) {
401  const G4FermiConfiguration* conf2 = (list2[A2])[id2];
402  G4int Z2 = conf2->GetZ();
403  const std::vector<const G4VFermiFragment*>* vec2 =
404  conf2->GetFragmentList();
405  G4int Z = Z1 + Z2;
406  G4int A = A1 + A2;
407  if(Z < maxZ && A < maxA && IsAvailable(Z, A)) {
408  newvec[0] = (*vec1)[0];
409  newvec[1] = (*vec1)[1];
410  newvec[2] = (*vec2)[0];
411  newvec[3] = (*vec2)[1];
412  if(!IsExist(Z, A, newvec)) {
413  const G4FermiConfiguration* conf4 =
414  new G4FermiConfiguration(newvec);
415  list4[A].push_back(conf4);
416  ++counter;
417  }
418  }
419  }
420  }
421  }
422  }
423  if(verbose > 0) {
424  tot += counter;
425  G4cout << G4endl;
426  G4cout << "### Quartets of fragments: " << counter << G4endl;
427  for(G4int A=4; A<maxA; ++A) {
428  G4cout << " A= " << A<<G4endl;
429  for(size_t j=0; j<list4[A].size(); ++j) {
430  const std::vector<const G4VFermiFragment*>* vec
431  = (list4[A])[j]->GetFragmentList();
432  G4int a1=(*vec)[0]->GetA();
433  G4int z1=(*vec)[0]->GetZ();
434  G4int a2=(*vec)[1]->GetA();
435  G4int z2=(*vec)[1]->GetZ();
436  G4int a3=(*vec)[2]->GetA();
437  G4int z3=(*vec)[2]->GetZ();
438  G4int a4=(*vec)[3]->GetA();
439  G4int z4=(*vec)[3]->GetZ();
440 
441  G4cout << "("<<a1<<","<<z1<<")("<<a2<<","<<z2<<")("<<a3<<","<<z3<<")("
442  <<a4<<","<<z4<<") % ";
443  }
444  G4cout<<G4endl;
445  G4cout<<"----------------------------------------------------------------"
446  << G4endl;
447  }
448  G4cout << "Total number: " << tot << G4endl;
449  }
450 }
451 
452 G4bool G4FermiFragmentsPool::IsApplicable(G4int Z, G4int A, G4double mass) const
453 {
454  if(Z >= maxZ || A >= maxA || A <= 0) { return false; }
455  // look into pair list
456  size_t nz = list2[A].size();
457  if(0 < nz) {
458  for(size_t j=0; j<nz; ++j) {
459  const G4FermiConfiguration* conf = (list2[A])[j];
460  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
461  }
462  }
463  // look into triple list
464  nz = list3[A].size();
465  if(0 < nz) {
466  for(size_t j=0; j<nz; ++j) {
467  const G4FermiConfiguration* conf = (list3[A])[j];
468  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
469  }
470  }
471  // look into quartet list
472  nz = list4[A].size();
473  if(0 < nz) {
474  for(size_t j=0; j<nz; ++j) {
475  const G4FermiConfiguration* conf = (list4[A])[j];
476  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
477  }
478  }
479 
480  // search in the pool and if found then return vector with one element
481  nz = list1[A].size();
482  if(0 < nz) {
483  for(size_t j=0; j<nz; ++j) {
484  const G4FermiConfiguration* conf = (list1[A])[j];
485  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
486  if(!(*(conf->GetFragmentList()))[0]->IsStable()) { return true; }
487  }
488  }
489  }
490  return false;
491 }
492 
493 const std::vector<const G4FermiConfiguration*>*
494 G4FermiFragmentsPool::GetConfigurationList(G4int Z, G4int A, G4double mass) const
495 {
496  std::vector<const G4FermiConfiguration*>* v =
497  new std::vector<const G4FermiConfiguration*>;
498  if(Z >= maxZ || A >= maxA) { return v; }
499 
500  //G4cout << "G4FermiFragmentsPool::GetConfigurationList:"
501  // << " Z= " << Z << " A= " << A << " Mass(GeV)= " << mass/GeV<< G4endl;
502 
503  // look into pair list
504  size_t nz = list2[A].size();
505  if(0 < nz) {
506  for(size_t j=0; j<nz; ++j) {
507  const G4FermiConfiguration* conf = (list2[A])[j];
508  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
509  v->push_back(conf);
510  }
511  //G4cout << "Pair dM(MeV)= " << mass - conf->GetMass() << G4endl; }
512  }
513  }
514  // look into triple list
515  nz = list3[A].size();
516  if(0 < nz) {
517  for(size_t j=0; j<nz; ++j) {
518  const G4FermiConfiguration* conf = (list3[A])[j];
519  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
520  v->push_back(conf);
521  }
522  //G4cout << "Triple dM(MeV)= " << mass - conf->GetMass() << G4endl; }
523  }
524  }
525  // look into quartet list
526  nz = list4[A].size();
527  if(0 < nz) {
528  for(size_t j=0; j<nz; ++j) {
529  const G4FermiConfiguration* conf = (list4[A])[j];
530  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
531  v->push_back(conf);
532  }
533  // G4cout << "Quartet dM(MeV)= " << mass - conf->GetMass() << G4endl; }
534  }
535  }
536  // return if vector not empty
537  if(0 < v->size()) {
538  if(verbose > 0) {
539  G4double ExEn= mass - G4NucleiProperties::GetNuclearMass(A,Z);
540  G4cout<<"Total number of configurations = "<<v->size()<<" for A= "
541  <<A<<" Z= "<<Z<<" E*= "<< ExEn<<" MeV"<<G4endl;
542  size_t size_vector_conf = v->size();
543  for(size_t jc=0; jc<size_vector_conf; ++jc) {
544  const std::vector<const G4VFermiFragment*>* v_frag =
545  (*v)[jc]->GetFragmentList();
546  size_t size_vector_fragments = v_frag->size();
547  G4cout<<size_vector_fragments<<"-body configuration "<<jc+1<<": ";
548  for(size_t jf=0;jf<size_vector_fragments;++jf){
549  G4int af= (*v_frag)[jf]->GetA();
550  G4int zf= (*v_frag)[jf]->GetZ();
551  G4double ex=(*v_frag)[jf]->GetExcitationEnergy();
552  G4cout<<"(a="<<af<<", z="<<zf<<", ex="<<ex<<") ";
553  }
554  G4cout<<G4endl;
555  G4cout<<"-----------------------------------------------------"<<G4endl;
556  }
557  }
558  return v;
559  }
560 
561  // search in the pool and if found then return vector with one element
562  nz = list1[A].size();
563  if(0 < nz) {
564  for(size_t j=0; j<nz; ++j) {
565  const G4FermiConfiguration* conf = (list1[A])[j];
566 
567  // G4cout << "Single dM(MeV)= " << mass - conf->GetMass() << G4endl; }
568  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
569  if(!(*(conf->GetFragmentList()))[0]->IsStable()) {
570  v->push_back(conf);
571 
572  if(verbose > 0) {
573  G4double ExEn= mass -G4NucleiProperties::GetNuclearMass(A,Z);
574  G4cout<<"Only 1 configurations for A= "
575  <<A<<" Z= "<<Z<<" E*= "<< ExEn<<" MeV"<<G4endl;
576  const std::vector<const G4VFermiFragment*>* v_frag
577  = (*v)[0]->GetFragmentList();
578  size_t size_vector_fragments=v_frag->size();
579  G4cout<<"1 Fragment configuration: ";
580  for(size_t jf=0;jf<size_vector_fragments;++jf){
581  G4int af= (*v_frag)[jf]->GetA();
582  G4int zf= (*v_frag)[jf]->GetZ();
583  G4double ex=(*v_frag)[jf]->GetExcitationEnergy();
584  G4cout<<"(a="<<af<<", z="<<zf<<", ex="<<ex<<") ";
585  }
586  G4cout<<G4endl;
587  G4cout<<"-----------------------------------------------------"<<G4endl;
588  }
589  return v;
590  }
591  }
592  }
593  }
594  //failer
595  if(verbose > 0) {
596  G4cout << "G4FermiFragmentsPool::GetConfigurationList: WARNING: not "
597  << "able decay fragment Z= " << Z << " A= " << A
598  << " Mass(GeV)= " << mass/GeV<< G4endl;
599  }
600  return v;
601 }
602 
603 G4bool G4FermiFragmentsPool::IsExist(G4int Z, G4int A,
604  std::vector<const G4VFermiFragment*>& newconf) const
605 {
606  size_t nn = newconf.size();
607  G4double mass = 0.0;
608  for(size_t i=0; i<nn; ++i) { mass += newconf[i]->GetTotalEnergy(); }
609  // look into pair list
610  if(2 == nn) {
611  size_t nz = list2[A].size();
612  if(0 < nz) {
613  for(size_t j=0; j<nz; ++j) {
614  const G4FermiConfiguration* conf = (list2[A])[j];
615  if(Z == conf->GetZ() && A == conf->GetA() &&
616  std::fabs(mass - conf->GetMass()) < keV) {return true; }
617  }
618  }
619  return false;
620  }
621  // look into triple list
622  if(3 == nn) {
623  size_t nz = list3[A].size();
624  if(0 < nz) {
625  for(size_t j=0; j<nz; ++j) {
626  const G4FermiConfiguration* conf = (list3[A])[j];
627  if(Z == conf->GetZ() && A == conf->GetA() &&
628  std::fabs(mass - conf->GetMass()) < keV) { return true; }
629  }
630  }
631  return false;
632  }
633  // look into quartet list
634  if(4 == nn) {
635  size_t nz = list4[A].size();
636  if(0 < nz) {
637  for(size_t j=0; j<nz; ++j) {
638  const G4FermiConfiguration* conf = (list4[A])[j];
639  if(Z == conf->GetZ() && A == conf->GetA() &&
640  std::fabs(mass - conf->GetMass()) < keV) { return true; }
641  }
642  }
643  return false;
644  }
645  return false;
646 }
647 
648 const G4VFermiFragment*
649 G4FermiFragmentsPool::GetFragment(G4int Z, G4int A) const
650 {
651  const G4VFermiFragment* f = 0;
652  if(Z < maxZ && A < maxA) {
653  size_t nz = list1[A].size();
654  for(size_t j=0; j<nz; ++j) {
655  const G4FermiConfiguration* conf = (list1[A])[j];
656  if(Z == conf->GetZ()) {
657  f = (*(conf->GetFragmentList()))[0];
658  break;
659  }
660  }
661  }
662  return f;
663 }
664 
665 void G4FermiFragmentsPool::DumpFragment(const G4VFermiFragment* f) const
666 {
667  if(f) {
668  G4cout << "Z= " << f->GetZ() << " A= " << f->GetA()
669  << " Mass(GeV)= " << f->GetFragmentMass()/GeV
670  << " Eexc(MeV)= " << f->GetExcitationEnergy() << G4endl;
671  }
672 }
673 
674 void G4FermiFragmentsPool::Dump() const
675 {
676  G4cout << "##### List of Fragments in the Fermi Fragment Pool #####"
677  << G4endl;
678  G4int nfrag = fragment_pool.size();
679  for(G4int i=0; i<nfrag; ++i) {
680  DumpFragment(fragment_pool[i]);
681  }
682  G4cout << G4endl;
683 }
G4B9FermiFragment.hh
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:64
G4Li5FermiFragment.hh
G4FermiFragmentsPool.hh
G4Be8FermiFragment.hh
G4FermiConfiguration::GetMass
G4double GetMass() const
Definition: G4FermiConfiguration.hh:79
G4int
int G4int
Definition: G4Types.hh:78
G4VFermiFragment::GetZ
G4int GetZ(void) const
Definition: G4VFermiFragment.hh:64
G4FermiFragmentsPool::GetFragment
const G4VFermiFragment * GetFragment(G4int Z, G4int A) const
Definition: G4FermiFragmentsPool.cc:649
G4StableFermiFragment.hh
G4VFermiFragment::GetA
G4int GetA(void) const
Definition: G4VFermiFragment.hh:59
G4He5FermiFragment.hh
G4FermiFragmentsPool::IsApplicable
G4bool IsApplicable(G4int Z, G4int A, G4double mass) const
Definition: G4FermiFragmentsPool.cc:452
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4bool
bool G4bool
Definition: G4Types.hh:79
G4FermiFragmentsPool::DumpFragment
void DumpFragment(const G4VFermiFragment *f) const
Definition: G4FermiFragmentsPool.cc:665
test.v
tuple v
Definition: test.py:18
G4PhysicalConstants.hh
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4VFermiFragment::GetExcitationEnergy
G4double GetExcitationEnergy(void) const
Definition: G4VFermiFragment.hh:74
G4VFermiFragment::GetFragmentMass
G4double GetFragmentMass(void) const
Definition: G4VFermiFragment.hh:79
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
instance
static MCTruthManager * instance
Definition: MCTruthManager.cc:47
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4FermiFragmentsPool::Instance
static G4FermiFragmentsPool * Instance()
Definition: G4FermiFragmentsPool.cc:56
G4FermiConfiguration::GetFragmentList
const std::vector< const G4VFermiFragment * > * GetFragmentList() const
Definition: G4FermiConfiguration.hh:85
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4FermiFragmentsPool::GetConfigurationList
const std::vector< const G4FermiConfiguration * > * GetConfigurationList(G4int Z, G4int A, G4double mass) const
Definition: G4FermiFragmentsPool.cc:494
G4FermiFragmentsPool::GetFermiPhaseSpaceDecay
const G4FermiPhaseSpaceDecay * GetFermiPhaseSpaceDecay() const
Definition: G4FermiFragmentsPool.cc:111
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262