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