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