Geant4  10.01.p02
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: G4VFermiBreakUp.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 = 0;
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  // look into pair list
455  size_t nz = list2[A].size();
456  if(0 < nz) {
457  for(size_t j=0; j<nz; ++j) {
458  const G4FermiConfiguration* conf = (list2[A])[j];
459  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
460  }
461  }
462  // look into triple list
463  nz = list3[A].size();
464  if(0 < nz) {
465  for(size_t j=0; j<nz; ++j) {
466  const G4FermiConfiguration* conf = (list3[A])[j];
467  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
468  }
469  }
470  // look into quartet list
471  nz = list4[A].size();
472  if(0 < nz) {
473  for(size_t j=0; j<nz; ++j) {
474  const G4FermiConfiguration* conf = (list4[A])[j];
475  if(Z == conf->GetZ() && mass >= conf->GetMass()) { return true; }
476  }
477  }
478 
479  // search in the pool and if found then return vector with one element
480  nz = list1[A].size();
481  if(0 < nz) {
482  for(size_t j=0; j<nz; ++j) {
483  const G4FermiConfiguration* conf = (list1[A])[j];
484  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
485  if(!(*(conf->GetFragmentList()))[0]->IsStable()) { return true; }
486  }
487  }
488  }
489  return false;
490 }
491 
492 const std::vector<const G4FermiConfiguration*>*
493 G4FermiFragmentsPool::GetConfigurationList(G4int Z, G4int A, G4double mass) const
494 {
495  std::vector<const G4FermiConfiguration*>* v =
496  new std::vector<const G4FermiConfiguration*>;
497  if(Z >= maxZ || A >= maxA) { return v; }
498 
499  //G4cout << "G4FermiFragmentsPool::GetConfigurationList:"
500  // << " Z= " << Z << " A= " << A << " Mass(GeV)= " << mass/GeV<< G4endl;
501 
502  // look into pair list
503  size_t nz = list2[A].size();
504  if(0 < nz) {
505  for(size_t j=0; j<nz; ++j) {
506  const G4FermiConfiguration* conf = (list2[A])[j];
507  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
508  v->push_back(conf);
509  }
510  //G4cout << "Pair dM(MeV)= " << mass - conf->GetMass() << G4endl; }
511  }
512  }
513  // look into triple list
514  nz = list3[A].size();
515  if(0 < nz) {
516  for(size_t j=0; j<nz; ++j) {
517  const G4FermiConfiguration* conf = (list3[A])[j];
518  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
519  v->push_back(conf);
520  }
521  //G4cout << "Triple dM(MeV)= " << mass - conf->GetMass() << G4endl; }
522  }
523  }
524  // look into quartet list
525  nz = list4[A].size();
526  if(0 < nz) {
527  for(size_t j=0; j<nz; ++j) {
528  const G4FermiConfiguration* conf = (list4[A])[j];
529  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
530  v->push_back(conf);
531  }
532  // G4cout << "Quartet dM(MeV)= " << mass - conf->GetMass() << G4endl; }
533  }
534  }
535  // return if vector not empty
536  if(0 < v->size()) {
537  if(verbose > 0) {
538  G4double ExEn= mass - G4NucleiProperties::GetNuclearMass(A,Z);
539  G4cout<<"Total number of configurations = "<<v->size()<<" for A= "
540  <<A<<" Z= "<<Z<<" E*= "<< ExEn<<" MeV"<<G4endl;
541  size_t size_vector_conf = v->size();
542  for(size_t jc=0; jc<size_vector_conf; ++jc) {
543  const std::vector<const G4VFermiFragment*>* v_frag =
544  (*v)[jc]->GetFragmentList();
545  size_t size_vector_fragments = v_frag->size();
546  G4cout<<size_vector_fragments<<"-body configuration "<<jc+1<<": ";
547  for(size_t jf=0;jf<size_vector_fragments;++jf){
548  G4int af= (*v_frag)[jf]->GetA();
549  G4int zf= (*v_frag)[jf]->GetZ();
550  G4double ex=(*v_frag)[jf]->GetExcitationEnergy();
551  G4cout<<"(a="<<af<<", z="<<zf<<", ex="<<ex<<") ";
552  }
553  G4cout<<G4endl;
554  G4cout<<"-----------------------------------------------------"<<G4endl;
555  }
556  }
557  return v;
558  }
559 
560  // search in the pool and if found then return vector with one element
561  nz = list1[A].size();
562  if(0 < nz) {
563  for(size_t j=0; j<nz; ++j) {
564  const G4FermiConfiguration* conf = (list1[A])[j];
565 
566  // G4cout << "Single dM(MeV)= " << mass - conf->GetMass() << G4endl; }
567  if(Z == conf->GetZ() && mass >= conf->GetMass()) {
568  if(!(*(conf->GetFragmentList()))[0]->IsStable()) {
569  v->push_back(conf);
570 
571  if(verbose > 0) {
572  G4double ExEn= mass -G4NucleiProperties::GetNuclearMass(A,Z);
573  G4cout<<"Only 1 configurations for A= "
574  <<A<<" Z= "<<Z<<" E*= "<< ExEn<<" MeV"<<G4endl;
575  const std::vector<const G4VFermiFragment*>* v_frag
576  = (*v)[0]->GetFragmentList();
577  size_t size_vector_fragments=v_frag->size();
578  G4cout<<"1 Fragment configuration: ";
579  for(size_t jf=0;jf<size_vector_fragments;++jf){
580  G4int af= (*v_frag)[jf]->GetA();
581  G4int zf= (*v_frag)[jf]->GetZ();
582  G4double ex=(*v_frag)[jf]->GetExcitationEnergy();
583  G4cout<<"(a="<<af<<", z="<<zf<<", ex="<<ex<<") ";
584  }
585  G4cout<<G4endl;
586  G4cout<<"-----------------------------------------------------"<<G4endl;
587  }
588  return v;
589  }
590  }
591  }
592  }
593  //failer
594  if(verbose > 0) {
595  G4cout << "G4FermiFragmentsPool::GetConfigurationList: WARNING: not "
596  << "able decay fragment Z= " << Z << " A= " << A
597  << " Mass(GeV)= " << mass/GeV<< G4endl;
598  }
599  return v;
600 }
601 
602 G4bool G4FermiFragmentsPool::IsExist(G4int Z, G4int A,
603  std::vector<const G4VFermiFragment*>& newconf) const
604 {
605  size_t nn = newconf.size();
606  G4double mass = 0.0;
607  for(size_t i=0; i<nn; ++i) { mass += newconf[i]->GetTotalEnergy(); }
608  // look into pair list
609  if(2 == nn) {
610  size_t nz = list2[A].size();
611  if(0 < nz) {
612  for(size_t j=0; j<nz; ++j) {
613  const G4FermiConfiguration* conf = (list2[A])[j];
614  if(Z == conf->GetZ() && A == conf->GetA() &&
615  std::fabs(mass - conf->GetMass()) < keV) {return true; }
616  }
617  }
618  return false;
619  }
620  // look into triple list
621  if(3 == nn) {
622  size_t nz = list3[A].size();
623  if(0 < nz) {
624  for(size_t j=0; j<nz; ++j) {
625  const G4FermiConfiguration* conf = (list3[A])[j];
626  if(Z == conf->GetZ() && A == conf->GetA() &&
627  std::fabs(mass - conf->GetMass()) < keV) { return true; }
628  }
629  }
630  return false;
631  }
632  // look into quartet list
633  if(4 == nn) {
634  size_t nz = list4[A].size();
635  if(0 < nz) {
636  for(size_t j=0; j<nz; ++j) {
637  const G4FermiConfiguration* conf = (list4[A])[j];
638  if(Z == conf->GetZ() && A == conf->GetA() &&
639  std::fabs(mass - conf->GetMass()) < keV) { return true; }
640  }
641  }
642  return false;
643  }
644  return false;
645 }
646 
647 const G4VFermiFragment*
648 G4FermiFragmentsPool::GetFragment(G4int Z, G4int A) const
649 {
650  const G4VFermiFragment* f = 0;
651  if(Z < maxZ && A < maxA) {
652  size_t nz = list1[A].size();
653  for(size_t j=0; j<nz; ++j) {
654  const G4FermiConfiguration* conf = (list1[A])[j];
655  if(Z == conf->GetZ()) {
656  f = (*(conf->GetFragmentList()))[0];
657  break;
658  }
659  }
660  }
661  return f;
662 }
663 
664 void G4FermiFragmentsPool::DumpFragment(const G4VFermiFragment* f) const
665 {
666  if(f) {
667  G4cout << "Z= " << f->GetZ() << " A= " << f->GetA()
668  << " Mass(GeV)= " << f->GetFragmentMass()/GeV
669  << " Eexc(MeV)= " << f->GetExcitationEnergy() << G4endl;
670  }
671 }
672 
673 void G4FermiFragmentsPool::Dump() const
674 {
675  G4cout << "##### List of Fragments in the Fermi Fragment Pool #####"
676  << G4endl;
677  G4int nfrag = fragment_pool.size();
678  for(G4int i=0; i<nfrag; ++i) {
679  DumpFragment(fragment_pool[i]);
680  }
681  G4cout << G4endl;
682 }
MeV
static const double MeV
Definition: G4SIunits.hh:193
G4B9FermiFragment.hh
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:53
G4Li5FermiFragment.hh
G4FermiFragmentsPool::thePhaseSpace
G4FermiPhaseSpaceDecay thePhaseSpace
Definition: G4FermiFragmentsPool.hh:93
G4FermiFragmentsPool.hh
a1
static const G4double a1
Definition: G4BetheHeitlerModel.cc:72
G4FermiFragmentsPool::list3
std::vector< const G4FermiConfiguration * > list3[17]
Definition: G4FermiFragmentsPool.hh:90
G4Be8FermiFragment.hh
G4FermiConfiguration::GetMass
G4double GetMass() const
Definition: G4FermiConfiguration.hh:79
a4
static const G4double a4
Definition: G4BetheHeitlerModel.cc:73
G4FermiFragmentsPool::list1
std::vector< const G4FermiConfiguration * > list1[17]
Definition: G4FermiFragmentsPool.hh:88
G4FermiFragmentsPool::theInstance
static G4FermiFragmentsPool * theInstance
Definition: G4FermiFragmentsPool.hh:79
G4int
int G4int
Definition: G4Types.hh:78
G4VFermiFragment::GetZ
G4int GetZ(void) const
Definition: G4VFermiFragment.hh:72
G4FermiFragmentsPool::GetFragment
const G4VFermiFragment * GetFragment(G4int Z, G4int A) const
Definition: G4FermiFragmentsPool.cc:648
G4StableFermiFragment.hh
G4FermiFragmentsPool::IsAvailable
G4bool IsAvailable(G4int Z, G4int A) const
Definition: G4FermiFragmentsPool.cc:90
G4VFermiFragment::GetA
G4int GetA(void) const
Definition: G4VFermiFragment.hh:67
G4He5FermiFragment.hh
G4FermiFragmentsPool::IsApplicable
G4bool IsApplicable(G4int Z, G4int A, G4double mass) const
Definition: G4FermiFragmentsPool.cc:452
G4cout
G4GLOB_DLL std::ostream G4cout
G4bool
bool G4bool
Definition: G4Types.hh:79
G4FermiFragmentsPool::fragment_pool
std::vector< const G4VFermiFragment * > fragment_pool
Definition: G4FermiFragmentsPool.hh:81
GeV
static const double GeV
Definition: G4SIunits.hh:196
G4FermiFragmentsPool::DumpFragment
void DumpFragment(const G4VFermiFragment *f) const
Definition: G4FermiFragmentsPool.cc:664
A
static const G4double A[nN]
Definition: G4ElectroNuclearCrossSection.cc:115
a3
static const G4double a3
Definition: G4BetheHeitlerModel.cc:73
G4PhysicalConstants.hh
G4VFermiFragment::GetExcitationEnergy
G4double GetExcitationEnergy(void) const
Definition: G4VFermiFragment.hh:82
G4VFermiFragment::GetFragmentMass
G4double GetFragmentMass(void) const
Definition: G4VFermiFragment.hh:87
G4endl
#define G4endl
Definition: G4ios.hh:61
instance
static MCTruthManager * instance
Definition: MCTruthManager.cc:43
keV
static const double keV
Definition: G4SIunits.hh:195
G4FermiFragmentsPool::list2
std::vector< const G4FermiConfiguration * > list2[17]
Definition: G4FermiFragmentsPool.hh:89
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4FermiFragmentsPool::IsExist
G4bool IsExist(G4int Z, G4int A, std::vector< const G4VFermiFragment * > &) const
Definition: G4FermiFragmentsPool.cc:602
G4FermiFragmentsPool::Instance
static G4FermiFragmentsPool * Instance()
Definition: G4FermiFragmentsPool.cc:56
G4FermiConfiguration::GetFragmentList
const std::vector< const G4VFermiFragment * > * GetFragmentList() const
Definition: G4FermiConfiguration.hh:85
G4FermiFragmentsPool::GetConfigurationList
const std::vector< const G4FermiConfiguration * > * GetConfigurationList(G4int Z, G4int A, G4double mass) const
Definition: G4FermiFragmentsPool.cc:493
G4FermiFragmentsPool::GetFermiPhaseSpaceDecay
const G4FermiPhaseSpaceDecay * GetFermiPhaseSpaceDecay() const
Definition: G4FermiFragmentsPool.cc:111
a2
static const G4double a2
Definition: G4BetheHeitlerModel.cc:72
G4FermiFragmentsPool::list4
std::vector< const G4FermiConfiguration * > list4[17]
Definition: G4FermiFragmentsPool.hh:91