source: trunk/source/processes/hadronic/models/cascade/evaporation/test/ratioTest.cc @ 1334

Last change on this file since 1334 was 1199, checked in by garnier, 15 years ago

nvx fichiers dans CVS
File size: 27.7 KB
Line 
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//
27// $Id: ratioTest.cc,v 1.4 2007/06/21 15:04:38 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Test program for G4 Bertini Evaporation.
31
32#define G4VERBOSE 1
33
34#include "G4LayeredNucleus.hh"
35#include "G4BertiniEvaporation.hh"
36
37int main(int argc, char *argv[])
38{
39  G4cout << "%*********************************" << endl
40         << "%*  Ratio test                   *" << endl
41         << "%*  Output in Latex table format *" << endl
42         << "%*********************************" << endl << endl;
43
44  G4LayeredNucleus nucl;
45  G4BertiniEvaporation bert;
46  G4VParticleChange * pc;
47  G4int A, Z;
48  G4double E;
49
50  G4int rounds=10000;
51  bert.setVerboseLevel(0);
52
53  G4int ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
54  G4int n=0, p=0, d=0, t=0, h3=0, h4=0, g=0, pn=0, _2n=0, _2p=0, alfan=0, p2n=0, _2pn=0;
55  G4int j;
56
57  G4cout << "% * Rounds : "<< rounds << endl;
58
59  ////////////////////////// Se 74 30, 35 MeV
60
61  A = 74;
62  Z = 34;
63  E = 35;
64 
65  nucl.SetParameters( A, Z );
66  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
67  n=p=d=t=h3=h4=g=0;
68
69  for ( j = 0 ; j < rounds ; j++)
70    {
71      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
72      pc = bert.BreakItUp( nucl );
73     
74      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
75        {
76          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
77          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
78          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
79          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
80          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
81          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
82          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
83          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
84         
85          delete pc->GetSecondary( i );
86        } // loop over particle change vector
87     
88      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
89      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
90      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
91           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
92      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
93      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
94      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
95               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
96      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
97      //x  G4cout << "% n p pn 2n 2p p2n * " << en
98      pc->Clear();
99      delete pc;
100     
101    }
102 
103  G4cout << "Se$^{74}$     &  " << E << "   &   pn/2n   &   1.7  & " << (G4double) pn/((G4double)_2n)
104         << " \\\\ % " << pn << "  " << _2n << endl;
105
106  ///////////////////////// Ge 68 32, E=20
107
108  A = 68;
109  Z = 32;
110  E = 20;
111
112  nucl.SetParameters( A, Z );
113  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
114  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
115  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=0;
116
117  for ( j = 0 ; j < rounds ; j++)
118    {
119      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
120      pc = bert.BreakItUp( nucl );
121     
122      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
123        {
124          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
125          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
126          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
127          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
128          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
129          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
130          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
131          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
132         
133          delete pc->GetSecondary( i );
134        } // loop over particle change vector
135     
136      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
137      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
138      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
139           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
140      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
141      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
142      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
143               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
144      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
145      //x  G4cout << "% n p pn 2n 2p p2n * " << en
146      pc->Clear();
147      delete pc;
148     
149    }
150 
151  G4cout << "Ge$^{68}$   &  " << E << "   &   p/n     &   1.76 & " << (G4double) p/((G4double)n)
152         << " \\\\ % " << p << "  " << n << endl;
153
154  ///////////////////////// Ge 68 32, E=35
155
156  A = 68;
157  Z = 32;
158  E = 35;
159
160  nucl.SetParameters( A, Z );
161  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
162  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
163  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=0;
164
165  for ( j = 0 ; j < rounds ; j++)
166    {
167      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
168      pc = bert.BreakItUp( nucl );
169     
170      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
171        {
172          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
173          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
174          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
175          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
176          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
177          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
178          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
179          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
180         
181          delete pc->GetSecondary( i );
182        } // loop over particle change vector
183     
184      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
185      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
186      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
187           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
188      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
189      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
190      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
191               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
192      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
193      //x  G4cout << "% n p pn 2n 2p p2n * " << en
194      pc->Clear();
195      delete pc;
196     
197    }
198 
199  G4cout << "            &  " << E << "   &   pn/2n   &   8.4  & " << (G4double) pn/((G4double)_2n)
200         << " \\\\ % " << pn << "  " << _2n << endl;
201
202  ///////////////////////// Ge 68 32, E=40
203
204  A = 68;
205  Z = 32;
206  E = 40;
207
208  nucl.SetParameters( A, Z );
209  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
210  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
211  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
212
213
214  for ( j = 0 ; j < rounds ; j++)
215    {
216      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
217      pc = bert.BreakItUp( nucl );
218     
219      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
220        {
221          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
222          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
223          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
224          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
225          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
226          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
227          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
228          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
229         
230          delete pc->GetSecondary( i );
231        } // loop over particle change vector
232     
233      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
234      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
235      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
236           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
237      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
238      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
239      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
240               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
241      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
242               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
243      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
244      //x  G4cout << "% n p pn 2n 2p p2n * " << en
245      pc->Clear();
246      delete pc;
247     
248    }
249 
250  G4cout << "            &  " << E << "   &   2pn/p2n &   5.8  & " << (G4double) _2pn/((G4double)p2n)
251         << " \\\\ % " << _2pn << "  " << p2n << endl;
252
253  ///////////////////////// Ge 67 31, E=35
254
255  A = 67;
256  Z = 31;
257  E = 35;
258
259  nucl.SetParameters( A, Z );
260  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
261  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
262  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
263
264
265  for ( j = 0 ; j < rounds ; j++)
266    {
267      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
268      pc = bert.BreakItUp( nucl );
269     
270      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
271        {
272          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
273          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
274          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
275          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
276          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
277          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
278          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
279          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
280         
281          delete pc->GetSecondary( i );
282        } // loop over particle change vector
283     
284      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
285      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
286      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
287           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
288      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
289      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
290      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
291               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
292      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
293               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
294      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
295      //x  G4cout << "% n p pn 2n 2p p2n * " << en
296      pc->Clear();
297      delete pc;
298     
299    }
300 
301  G4cout << "Ga^{67}     &  " << E << "   &   pn/2n   &   3.3  & " << (G4double) pn/((G4double)_2n)
302         << " \\\\ % " << pn << "  " << _2n << endl;
303
304 ///////////////////////// Ni 58 28
305
306  A = 58;
307  Z = 28;
308  E = 25;
309
310  nucl.SetParameters( A, Z );
311  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
312  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
313  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
314
315
316  for ( j = 0 ; j < rounds ; j++)
317    {
318      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
319      pc = bert.BreakItUp( nucl );
320     
321      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
322        {
323          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
324          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
325          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
326          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
327          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
328          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
329          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
330          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
331         
332          delete pc->GetSecondary( i );
333        } // loop over particle change vector
334     
335      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
336      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
337      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
338           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
339      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
340      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
341      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
342               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
343      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
344               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
345      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
346      //x  G4cout << "% n p pn 2n 2p p2n * " << en
347      pc->Clear();
348      delete pc;
349     
350    }
351 
352  G4cout << "Ni^{58}     &  " << E << "   &   p/n     &   3.8  & " << (G4double) p/((G4double)n)
353         << " \\\\ % " << p << "  " << n << endl;
354
355 /////////////////////////
356
357  ///////////////////////// Ni 58 28, 35
358
359  A = 58;
360  Z = 28;
361  E = 35;
362
363  nucl.SetParameters( A, Z );
364  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
365  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
366  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
367
368
369  for ( j = 0 ; j < rounds ; j++)
370    {
371      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
372      pc = bert.BreakItUp( nucl );
373     
374      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
375        {
376          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
377          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
378          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
379          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
380          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
381          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
382          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
383          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
384         
385          delete pc->GetSecondary( i );
386        } // loop over particle change vector
387     
388      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
389      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
390      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
391           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
392      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
393      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
394      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
395               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
396      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
397               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
398      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
399      //x  G4cout << "% n p pn 2n 2p p2n * " << en
400      pc->Clear();
401      delete pc;
402     
403    }
404 
405  G4cout << "Ni^{58}     &  " << E << "   &   pn/2n   &   67   & " << (G4double) pn/((G4double)_2n)
406         << " \\\\ % " << pn << "  " << _2n << endl;
407
408  ///////////////////////// Ni 58 28 40
409
410  A = 58;
411  Z = 28;
412  E = 40;
413
414  nucl.SetParameters( A, Z );
415  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
416  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
417  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
418
419
420  for ( j = 0 ; j < rounds ; j++)
421    {
422      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
423      pc = bert.BreakItUp( nucl );
424     
425      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
426        {
427          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
428          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
429          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
430          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
431          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
432          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
433          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
434          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
435         
436          delete pc->GetSecondary( i );
437        } // loop over particle change vector
438     
439      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
440      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
441      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
442           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
443      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
444      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
445      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
446               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
447      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
448               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
449      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
450      //x  G4cout << "% n p pn 2n 2p p2n * " << en
451      pc->Clear();
452      delete pc;
453     
454    }
455 
456  G4cout << "Ni^{58}     &  " << E << "   &   2pn/p2n &   7.3  & " << (G4double) _2pn/((G4double)p2n)
457         << " \\\\ % " << _2pn << "  " << p2n << endl;
458
459  ///////////////////////// Fe 54 26 35
460
461  A = 54;
462  Z = 26;
463  E = 35;
464
465  nucl.SetParameters( A, Z );
466  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
467  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
468  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
469
470
471  for ( j = 0 ; j < rounds ; j++)
472    {
473      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
474      pc = bert.BreakItUp( nucl );
475     
476      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
477        {
478          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
479          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
480          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
481          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
482          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
483          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
484          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
485          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
486         
487          delete pc->GetSecondary( i );
488        } // loop over particle change vector
489     
490      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
491      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
492      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
493           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
494      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
495      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
496      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
497               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
498      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
499               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
500      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
501      //x  G4cout << "% n p pn 2n 2p p2n * " << en
502      pc->Clear();
503      delete pc;
504     
505    }
506 
507  G4cout << "Fe^{54}     &  " << E << "   &   pn/2n   &   31   & " << (G4double) pn/((G4double)_2n)
508         << " \\\\ % " << pn << "  " << _2n << endl;
509
510  ///
511  ///////////////////////// Mn 54 25 30
512
513  A = 54;
514  Z = 25;
515  E = 30;
516
517  nucl.SetParameters( A, Z );
518  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
519  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
520  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
521
522
523  for ( j = 0 ; j < rounds ; j++)
524    {
525      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
526      pc = bert.BreakItUp( nucl );
527     
528      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
529        {
530          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
531          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
532          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
533          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
534          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
535          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
536          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
537          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
538         
539          delete pc->GetSecondary( i );
540        } // loop over particle change vector
541     
542      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
543      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
544      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
545           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
546      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
547      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
548      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
549               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
550      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
551               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
552      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
553      //x  G4cout << "% n p pn 2n 2p p2n * " << en
554      pc->Clear();
555      delete pc;
556     
557    }
558 
559  G4cout << "Mn^{54}     &  " << E << "   &   pn/2n   &   <1   & " << (G4double) pn/((G4double)_2n)
560         << " \\\\ % " << pn << "  " << _2n << endl;
561
562 ///////////////////////// Cr 50 24
563
564  A = 50;
565  Z = 24;
566  E = 25;
567
568  nucl.SetParameters( A, Z );
569  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
570  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
571  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
572
573
574  for ( j = 0 ; j < rounds ; j++)
575    {
576      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
577      pc = bert.BreakItUp( nucl );
578     
579      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
580        {
581          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
582          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
583          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
584          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
585          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
586          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
587          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
588          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
589         
590          delete pc->GetSecondary( i );
591        } // loop over particle change vector
592     
593      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
594      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
595      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
596           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
597      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
598      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
599      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
600               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
601      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
602               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
603      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
604      //x  G4cout << "% n p pn 2n 2p p2n * " << en
605      pc->Clear();
606      delete pc;
607     
608    }
609 
610  G4cout << "Cr^{50}     &  " << E << "   &   p/n     &   0.51 & " << (G4double) p/((G4double)n)
611         << " \\\\ % " << p << "  " << n << endl;
612
613  ///////////////////////// Cr 50 24 35
614
615  A = 50;
616  Z = 24;
617  E = 35;
618
619  nucl.SetParameters( A, Z );
620  nucl.AddExcitationEnergy( E - nucl.GetEnergyDeposit() );
621  ntemp=0, ptemp=0, dtemp=0, ttemp=0, h3temp=0, h4temp=0, gtemp=0;
622  n=p=d=t=h3=h4=g=pn=_2n=_2p=alfan=p2n=_2pn=0;
623
624
625  for ( j = 0 ; j < rounds ; j++)
626    {
627      ntemp=ptemp=dtemp=ttemp=h3temp=h4temp=gtemp=0;
628      pc = bert.BreakItUp( nucl );
629     
630      for ( G4int i = 0 ; i < pc->GetNumberOfSecondaries() ;  i++ )
631        {
632          const char * name = pc->GetSecondary( i )->GetDefinition()->GetParticleName() ;
633          if ( strcmp ( name , "neutron" ) == 0 ) ntemp++;
634          else if ( strcmp ( name , "proton" ) == 0 ) ptemp++;
635          else if ( strcmp ( name , "deuteron" ) == 0 ) dtemp++;
636          else if ( strcmp ( name , "triton" ) == 0 ) ttemp++;
637          else if ( strcmp ( name , "He3" ) == 0 ) h3temp++;
638          else if ( strcmp ( name , "alpha" ) == 0 ) h4temp++;
639          else if ( strcmp ( name , "gamma" ) == 0 ) gtemp++;
640         
641          delete pc->GetSecondary( i );
642        } // loop over particle change vector
643     
644      if ( ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) n++;
645      if ( ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) p++;
646      if ( ( ntemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ||
647           ( dtemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 1 ) ) pn++;
648      if ( ntemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2n++;
649      if ( ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 2 ) _2p++;
650      if ( ( ntemp == 2 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
651               ( dtemp == 1 && ntemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) ) p2n++;
652      if ( ( ntemp == 1 && ptemp == 2 && pc->GetNumberOfSecondaries() - gtemp == 3 ) ||
653               ( dtemp == 1 && ptemp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) )_2pn++;
654      if ( ntemp == 1 && h4temp == 1 && pc->GetNumberOfSecondaries() - gtemp == 2 ) alfan++;
655      //x  G4cout << "% n p pn 2n 2p p2n * " << en
656      pc->Clear();
657      delete pc;
658     
659    }
660 
661  G4cout << "            &  " << E << "   &   pn/2n   &   46   & " << (G4double) pn/((G4double)_2n)
662         << " \\\\ % " << pn << "  " << _2n << endl;
663
664 
665  return 0;
666}
667
668
669
670
671
672
673
Note: See TracBrowser for help on using the repository browser.