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