source: trunk/source/processes/hadronic/models/de_excitation/evaporation/test/G4EvaporationProbabilityTest.cc @ 1199

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

nvx fichiers dans CVS
File size: 18.9 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// -------------------------------------------------------------------
28//      GEANT 4 class file --- Copyright CERN 1998
29//      CERN Geneva Switzerland
30//
31//
32//      File name:     G4EvaporationProbabilityTest.cc
33//
34//      Author:        Maria Grazia Pia (pia@genova.infn.it),
35//
36//      Creation date: 27 October 1998
37//
38//      Modifications:
39//
40// -------------------------------------------------------------------
41
42#include "globals.hh"
43
44#include "G4ios.hh"
45#include <fstream>
46#include <iomanip>
47#include <iostream>
48
49#include "CLHEP/Hist/TupleManager.h"
50#include "CLHEP/Hist/HBookFile.h"
51#include "CLHEP/Hist/Histogram.h"
52#include "CLHEP/Hist/Tuple.h"
53#include "CLHEP/String/Strings.h"
54#include "G4VEvaporationChannel.hh"
55#include "G4CompetitiveFission.hh"
56
57#include "G4VEvaporationChannel.hh"
58#include "G4EvaporationChannel.hh"
59#include "G4PhotonEvaporation.hh"
60#include "G4DataVector.hh"
61#include "G4LorentzVector.hh"
62#include "G4ThreeVector.hh"
63#include "G4NucleiProperties.hh"
64#include "G4Fragment.hh"
65#include "Randomize.hh"
66#include "G4FragmentVector.hh"
67
68#include "g4rw/tvvector.h"
69
70int main()
71{
72  // MGP ---- HBOOK initialization
73  HepTupleManager* hbookManager;
74  hbookManager = new HBookFile("prob.hbook", 58);
75  assert (hbookManager != 0);
76
77  // MGP ---- Book histograms
78
79  HepHistogram* hGammaE;
80  hGammaE = hbookManager->histogram("Gamma energy", 100,0.,10.);
81  assert (hGammaE != 0); 
82
83  HepHistogram* hNProducts;
84  hNProducts = hbookManager->histogram("Number of products", 20,0.,20.);
85  assert (hNProducts != 0); 
86
87  HepHistogram* hProb;
88  hProb = hbookManager->histogram("Probability * 1.e25", 100,0.,10.);
89  assert (hProb != 0); 
90
91  // MGP ---- Book a ntuple
92  HepTuple* ntuple;
93  ntuple = hbookManager->ntuple("G4EvaporationProbability");
94  assert (ntuple != 0);
95
96  G4int Z;
97  G4int A;
98
99  G4cout << "Enter Z and A" << G4endl;
100  G4cin >> Z >> A;
101
102  assert (Z > 0);
103  assert (A > 0);
104  assert (A > Z);
105 
106  G4int iter;
107  G4cout << "Enter number of iterations " << G4endl;
108  G4cin >> iter;
109  if (iter <1) iter = 1;
110
111  G4int nProbs;
112  G4cout << "Enter max number of channels (0 - 32) " << G4endl;
113  G4cin >> nProbs;
114  if (nProbs < 0 || nProbs > 32) nProbs = 32;
115
116  G4double excMin;
117  G4double excMax;   
118  G4cout << "Enter initial min an max excitation energy" << G4endl;
119  G4cin >> excMin >> excMax;
120  assert (excMin >= 0.);
121  assert (excMax > 0.);
122  assert (excMax >= excMin);
123
124
125  // G4EvaporationProbability for this (Z,A) material
126
127  // Excitation energy levels for each channel
128
129  G4int nExcitedStates = 10;
130
131  G4double zero = 0.;
132
133  G4RWTValVector<G4double> ExcitEnergyChann00(nExcitedStates,zero);
134  G4RWTValVector<G4double> ExcitEnergyChann01(nExcitedStates,zero);
135  G4RWTValVector<G4double> ExcitEnergyChann02(nExcitedStates,zero);
136  G4RWTValVector<G4double> ExcitEnergyChann03(nExcitedStates,zero);
137  G4RWTValVector<G4double> ExcitEnergyChann04(nExcitedStates,zero);
138  G4RWTValVector<G4double> ExcitEnergyChann05(nExcitedStates,zero);
139  G4RWTValVector<G4double> ExcitEnergyChann06(nExcitedStates,zero);
140  G4RWTValVector<G4double> ExcitEnergyChann07(nExcitedStates,zero);
141  G4RWTValVector<G4double> ExcitEnergyChann08(nExcitedStates,zero);
142  G4RWTValVector<G4double> ExcitEnergyChann09(nExcitedStates,zero);
143  G4RWTValVector<G4double> ExcitEnergyChann10(nExcitedStates,zero);
144  G4RWTValVector<G4double> ExcitEnergyChann11(nExcitedStates,zero);
145  G4RWTValVector<G4double> ExcitEnergyChann12(nExcitedStates,zero);
146  G4RWTValVector<G4double> ExcitEnergyChann13(nExcitedStates,zero);
147  G4RWTValVector<G4double> ExcitEnergyChann14(nExcitedStates,zero);
148  G4RWTValVector<G4double> ExcitEnergyChann15(nExcitedStates,zero);
149  G4RWTValVector<G4double> ExcitEnergyChann16(nExcitedStates,zero);
150  G4RWTValVector<G4double> ExcitEnergyChann17(nExcitedStates,zero);
151  G4RWTValVector<G4double> ExcitEnergyChann18(nExcitedStates,zero);
152  G4RWTValVector<G4double> ExcitEnergyChann19(nExcitedStates,zero);
153  G4RWTValVector<G4double> ExcitEnergyChann20(nExcitedStates,zero);
154  G4RWTValVector<G4double> ExcitEnergyChann21(nExcitedStates,zero);
155  G4RWTValVector<G4double> ExcitEnergyChann22(nExcitedStates,zero);
156  G4RWTValVector<G4double> ExcitEnergyChann23(nExcitedStates,zero);
157  G4RWTValVector<G4double> ExcitEnergyChann24(nExcitedStates,zero);
158  G4RWTValVector<G4double> ExcitEnergyChann25(nExcitedStates,zero);
159  G4RWTValVector<G4double> ExcitEnergyChann26(nExcitedStates,zero);
160  G4RWTValVector<G4double> ExcitEnergyChann27(nExcitedStates,zero);
161  G4RWTValVector<G4double> ExcitEnergyChann28(nExcitedStates,zero);
162  G4RWTValVector<G4double> ExcitEnergyChann29(nExcitedStates,zero);
163  G4RWTValVector<G4double> ExcitEnergyChann30(nExcitedStates,zero);
164  G4RWTValVector<G4double> ExcitEnergyChann31(nExcitedStates,zero);
165
166
167  // Spin of excitation energy levels for each channel
168
169  G4int izero = 0;
170
171  G4RWTValVector<G4int> ExcitSpinChann00(nExcitedStates,izero);
172  G4RWTValVector<G4int> ExcitSpinChann01(nExcitedStates,izero);
173  G4RWTValVector<G4int> ExcitSpinChann02(nExcitedStates,izero);
174  G4RWTValVector<G4int> ExcitSpinChann03(nExcitedStates,izero);
175  G4RWTValVector<G4int> ExcitSpinChann04(nExcitedStates,izero);
176  G4RWTValVector<G4int> ExcitSpinChann05(nExcitedStates,izero);
177  G4RWTValVector<G4int> ExcitSpinChann06(nExcitedStates,izero);
178  G4RWTValVector<G4int> ExcitSpinChann07(nExcitedStates,izero);
179  G4RWTValVector<G4int> ExcitSpinChann08(nExcitedStates,izero);
180  G4RWTValVector<G4int> ExcitSpinChann09(nExcitedStates,izero);
181  G4RWTValVector<G4int> ExcitSpinChann10(nExcitedStates,izero);
182  G4RWTValVector<G4int> ExcitSpinChann11(nExcitedStates,izero);
183  G4RWTValVector<G4int> ExcitSpinChann12(nExcitedStates,izero);
184  G4RWTValVector<G4int> ExcitSpinChann13(nExcitedStates,izero);
185  G4RWTValVector<G4int> ExcitSpinChann14(nExcitedStates,izero);
186  G4RWTValVector<G4int> ExcitSpinChann15(nExcitedStates,izero);
187  G4RWTValVector<G4int> ExcitSpinChann16(nExcitedStates,izero);
188  G4RWTValVector<G4int> ExcitSpinChann17(nExcitedStates,izero);
189  G4RWTValVector<G4int> ExcitSpinChann18(nExcitedStates,izero);
190  G4RWTValVector<G4int> ExcitSpinChann19(nExcitedStates,izero);
191  G4RWTValVector<G4int> ExcitSpinChann20(nExcitedStates,izero);
192  G4RWTValVector<G4int> ExcitSpinChann21(nExcitedStates,izero);
193  G4RWTValVector<G4int> ExcitSpinChann22(nExcitedStates,izero);
194  G4RWTValVector<G4int> ExcitSpinChann23(nExcitedStates,izero);
195  G4RWTValVector<G4int> ExcitSpinChann24(nExcitedStates,izero);
196  G4RWTValVector<G4int> ExcitSpinChann25(nExcitedStates,izero);
197  G4RWTValVector<G4int> ExcitSpinChann26(nExcitedStates,izero);
198  G4RWTValVector<G4int> ExcitSpinChann27(nExcitedStates,izero);
199  G4RWTValVector<G4int> ExcitSpinChann28(nExcitedStates,izero);
200  G4RWTValVector<G4int> ExcitSpinChann29(nExcitedStates,izero);
201  G4RWTValVector<G4int> ExcitSpinChann30(nExcitedStates,izero);
202  G4RWTValVector<G4int> ExcitSpinChann31(nExcitedStates,izero);
203
204
205  // (in MeV)
206
207  ExcitEnergyChann09(0) = 3.56;
208
209  ExcitEnergyChann10(0) = 0.48;
210
211  ExcitEnergyChann11(0) = 0.98;
212
213  ExcitEnergyChann12(0) = 0.43;
214 
215  ExcitEnergyChann15(0) = 3.37;
216  ExcitEnergyChann15(1) = 5.96;
217  ExcitEnergyChann15(2) = 6.18;
218  ExcitEnergyChann15(3) = 6.26;
219
220  ExcitEnergyChann17(0) = 0.72;
221  ExcitEnergyChann17(1) = 1.74;
222  ExcitEnergyChann17(2) = 2.15;
223  ExcitEnergyChann17(3) = 3.59;
224
225  ExcitEnergyChann18(0) = 2.13;
226  ExcitEnergyChann18(1) = 4.44;
227  ExcitEnergyChann18(2) = 5.02;
228  ExcitEnergyChann18(3) = 6.76;
229  ExcitEnergyChann18(4) = 7.29;
230  ExcitEnergyChann18(5) = 7.98;
231  ExcitEnergyChann18(6) = 8.56;
232
233  ExcitEnergyChann19(0) = 0.95;
234  ExcitEnergyChann19(1) = 1.67;
235  ExcitEnergyChann19(2) = 6.25;
236
237  ExcitEnergyChann20(0) = 2.00;
238  ExcitEnergyChann20(1) = 4.32;
239  ExcitEnergyChann20(2) = 4.80;
240  ExcitEnergyChann20(3) = 6.34;
241  ExcitEnergyChann20(4) = 6.48;
242  ExcitEnergyChann20(5) = 6.90;
243  ExcitEnergyChann20(6) = 7.50;
244  ExcitEnergyChann20(7) = 8.10;
245  ExcitEnergyChann20(8) = 8.42;
246  ExcitEnergyChann20(9) = 8.66;
247
248  ExcitEnergyChann21(0) = 4.44;
249
250  ExcitEnergyChann22(0) = 3.09;
251  ExcitEnergyChann22(1) = 3.68;
252  ExcitEnergyChann22(2) = 3.85;
253
254  ExcitEnergyChann23(0) = 6.09;
255  ExcitEnergyChann23(1) = 6.69;
256  ExcitEnergyChann23(2) = 6.96;
257  ExcitEnergyChann23(3) = 7.34;
258
259  ExcitEnergyChann25(0) = 2.31;
260  ExcitEnergyChann25(1) = 3.95;
261  ExcitEnergyChann25(2) = 4.92;
262  ExcitEnergyChann25(3) = 5.11;
263  ExcitEnergyChann25(4) = 5.69;
264  ExcitEnergyChann25(5) = 5.83;
265  ExcitEnergyChann25(6) = 6.20;
266  ExcitEnergyChann25(7) = 6.44;
267  ExcitEnergyChann25(8) = 7.03;
268
269  ExcitEnergyChann26(0) = 5.28;
270  ExcitEnergyChann26(1) = 6.32;
271  ExcitEnergyChann26(2) = 7.22;
272  ExcitEnergyChann26(3) = 7.57;
273  ExcitEnergyChann26(4) = 8.31;
274  ExcitEnergyChann26(5) = 8.57;
275  ExcitEnergyChann26(6) = 9.15;
276  ExcitEnergyChann26(7) = 9.79;
277  ExcitEnergyChann26(8) = 10.0;
278
279  ExcitEnergyChann27(0) = 0.12;
280  ExcitEnergyChann27(1) = 0.30;
281  ExcitEnergyChann27(2) = 0.40;
282
283  ExcitEnergyChann28(0) = 5.22;
284  ExcitEnergyChann28(1) = 6.18;
285  ExcitEnergyChann28(2) = 6.83;
286  ExcitEnergyChann28(3) = 7.28;
287
288  ExcitEnergyChann29(0) = 6.10;
289  ExcitEnergyChann29(1) = 6.92;
290  ExcitEnergyChann29(2) = 7.12;
291
292  ExcitEnergyChann30(0) = 0.87;
293  ExcitEnergyChann30(1) = 3.06;
294  ExcitEnergyChann30(2) = 3.84;
295
296  ExcitEnergyChann31(0) = 1.98;
297  ExcitEnergyChann31(1) = 3.57;
298  ExcitEnergyChann31(2) = 3.92;
299  ExcitEnergyChann31(3) = 4.46;
300  ExcitEnergyChann31(4) = 5.10;
301  ExcitEnergyChann31(5) = 5.33;
302  ExcitEnergyChann31(6) = 5.53;
303  ExcitEnergyChann31(7) = 6.20;
304  ExcitEnergyChann31(8) = 6.38;
305  ExcitEnergyChann31(9) = 6.88;
306
307
308  ExcitSpinChann09(0) = 1;
309
310  ExcitSpinChann10(0) = 2;
311
312  ExcitSpinChann11(0) = 3;
313
314  ExcitSpinChann12(0) = 2;
315
316  ExcitSpinChann15(0) = 5;
317  ExcitSpinChann15(1) = 8;
318  ExcitSpinChann15(2) = 1;
319  ExcitSpinChann15(3) = 5;
320
321  ExcitSpinChann17(0) = 3;
322  ExcitSpinChann17(1) = 1;
323  ExcitSpinChann17(2) = 3;
324  ExcitSpinChann17(3) = 5;
325
326  ExcitSpinChann18(0) = 2;
327  ExcitSpinChann18(1) = 6;
328  ExcitSpinChann18(2) = 4;
329  ExcitSpinChann18(3) = 10;
330  ExcitSpinChann18(4) = 6;
331  ExcitSpinChann18(5) = 4;
332  ExcitSpinChann18(6) = 6;
333
334  ExcitSpinChann19(0) = 5;
335  ExcitSpinChann19(1) = 5;
336  ExcitSpinChann19(2) = 4;
337
338  ExcitSpinChann20(0) = 2;
339  ExcitSpinChann20(1) = 6;
340  ExcitSpinChann20(2) = 4;
341  ExcitSpinChann20(3) = 2;
342  ExcitSpinChann20(4) = 8;
343  ExcitSpinChann20(5) = 6;
344  ExcitSpinChann20(6) = 4;
345  ExcitSpinChann20(7) = 4;
346  ExcitSpinChann20(8) = 6;
347  ExcitSpinChann20(9) = 8;
348
349  ExcitSpinChann21(0) = 5;
350
351  ExcitSpinChann22(0) = 2;
352  ExcitSpinChann22(1) = 4;
353  ExcitSpinChann22(2) = 6;
354
355  ExcitSpinChann23(0) = 3;
356  ExcitSpinChann23(1) = 8;
357  ExcitSpinChann23(2) = 6;
358  ExcitSpinChann23(3) = 5;
359
360  ExcitSpinChann25(0) = 1;
361  ExcitSpinChann25(1) = 3;
362  ExcitSpinChann25(2) = 1;
363  ExcitSpinChann25(3) = 5;
364  ExcitSpinChann25(4) = 3;
365  ExcitSpinChann25(5) = 7;
366  ExcitSpinChann25(6) = 3;
367  ExcitSpinChann25(7) = 7;
368  ExcitSpinChann25(8) = 5;
369
370  ExcitSpinChann26(0) = 8;
371  ExcitSpinChann26(1) = 4;
372  ExcitSpinChann26(2) = 10;
373  ExcitSpinChann26(3) = 8;
374  ExcitSpinChann26(4) = 2;
375  ExcitSpinChann26(5) = 4;
376  ExcitSpinChann26(6) = 14;
377  ExcitSpinChann26(7) = 14;
378  ExcitSpinChann26(8) = 8;
379
380  ExcitSpinChann27(0) = 1;
381  ExcitSpinChann27(1) = 7;
382  ExcitSpinChann27(2) = 3;
383
384  ExcitSpinChann28(0) = 8;
385  ExcitSpinChann28(1) = 4;
386  ExcitSpinChann28(2) = 10;
387  ExcitSpinChann28(3) = 8;
388
389  ExcitSpinChann29(0) = 8;
390  ExcitSpinChann29(1) = 5;
391  ExcitSpinChann29(2) = 3;
392
393  ExcitSpinChann30(0) = 2;
394  ExcitSpinChann30(1) = 2;
395  ExcitSpinChann30(2) = 6;
396
397  ExcitSpinChann31(0) = 5;
398  ExcitSpinChann31(1) = 10;
399  ExcitSpinChann31(2) = 5;
400  ExcitSpinChann31(3) = 3;
401  ExcitSpinChann31(4) = 7;
402  ExcitSpinChann31(5) = 13;
403  ExcitSpinChann31(6) = 5;
404  ExcitSpinChann31(7) = 3;
405  ExcitSpinChann31(8) = 12;
406  ExcitSpinChann31(9) = 1;
407
408  G4RWTPtrOrderedVector<G4VEvaporationChannel>* theChannels = new G4RWTPtrOrderedVector<G4VEvaporationChannel>;
409
410  //                                        |Gamma|A| Z|
411  //                                        +-----+-+--+
412
413  theChannels->insert(new G4EvaporationChannel(  2,  1, 0, &ExcitEnergyChann00, &ExcitSpinChann00)); // n
414  theChannels->insert(new G4EvaporationChannel(  2,  1, 1, &ExcitEnergyChann01, &ExcitSpinChann01)); // p
415  theChannels->insert(new G4EvaporationChannel(  6,  2, 1, &ExcitEnergyChann02, &ExcitSpinChann02)); // H2
416  theChannels->insert(new G4EvaporationChannel(  6,  3, 1, &ExcitEnergyChann03, &ExcitSpinChann03)); // H3
417  theChannels->insert(new G4EvaporationChannel(  6,  3, 2, &ExcitEnergyChann04, &ExcitSpinChann04)); // He3
418  theChannels->insert(new G4EvaporationChannel(  4,  4, 2, &ExcitEnergyChann05, &ExcitSpinChann05)); // He4
419  theChannels->insert(new G4EvaporationChannel( 20,  5, 2, &ExcitEnergyChann06, &ExcitSpinChann06)); // He5
420  theChannels->insert(new G4EvaporationChannel( 30,  6, 2, &ExcitEnergyChann07, &ExcitSpinChann07)); // He6
421  theChannels->insert(new G4EvaporationChannel( 20,  5, 3, &ExcitEnergyChann08, &ExcitSpinChann08)); // Li5
422  theChannels->insert(new G4EvaporationChannel( 54,  6, 3, &ExcitEnergyChann09, &ExcitSpinChann09)); // Li6
423  theChannels->insert(new G4EvaporationChannel( 73,  7, 3, &ExcitEnergyChann10, &ExcitSpinChann10)); // Li7
424  theChannels->insert(new G4EvaporationChannel(101,  8, 3, &ExcitEnergyChann11, &ExcitSpinChann11)); // Li8
425  theChannels->insert(new G4EvaporationChannel( 73,  7, 4, &ExcitEnergyChann12, &ExcitSpinChann12)); // Be7
426  theChannels->insert(new G4EvaporationChannel(  8,  8, 4, &ExcitEnergyChann13, &ExcitSpinChann13)); // Be8
427  theChannels->insert(new G4EvaporationChannel(146,  9, 4, &ExcitEnergyChann14, &ExcitSpinChann14)); // Be9
428  theChannels->insert(new G4EvaporationChannel(100, 10, 4, &ExcitEnergyChann15, &ExcitSpinChann15)); // Be10
429  theChannels->insert(new G4EvaporationChannel(100,  9, 5, &ExcitEnergyChann16, &ExcitSpinChann16)); // B9
430  theChannels->insert(new G4EvaporationChannel(343, 10, 5, &ExcitEnergyChann17, &ExcitSpinChann17)); // B10
431  theChannels->insert(new G4EvaporationChannel(174, 11, 5, &ExcitEnergyChann18, &ExcitSpinChann18)); // B11
432  theChannels->insert(new G4EvaporationChannel(393, 12, 5, &ExcitEnergyChann19, &ExcitSpinChann19)); // B12
433  theChannels->insert(new G4EvaporationChannel(186, 11, 6, &ExcitEnergyChann20, &ExcitSpinChann20)); // C11
434  theChannels->insert(new G4EvaporationChannel( 61, 12, 6, &ExcitEnergyChann21, &ExcitSpinChann21)); // C12
435  theChannels->insert(new G4EvaporationChannel(202, 13, 6, &ExcitEnergyChann22, &ExcitSpinChann22)); // C13
436  theChannels->insert(new G4EvaporationChannel(113, 14, 6, &ExcitEnergyChann23, &ExcitSpinChann23)); // C14
437  theChannels->insert(new G4EvaporationChannel(213, 13, 7, &ExcitEnergyChann24, &ExcitSpinChann24)); // N13
438  theChannels->insert(new G4EvaporationChannel(233, 14, 7, &ExcitEnergyChann25, &ExcitSpinChann25)); // N14
439  theChannels->insert(new G4EvaporationChannel(180, 15, 7, &ExcitEnergyChann26, &ExcitSpinChann26)); // N15
440  theChannels->insert(new G4EvaporationChannel(696, 16, 7, &ExcitEnergyChann27, &ExcitSpinChann27)); // N16
441  theChannels->insert(new G4EvaporationChannel(194, 15, 8, &ExcitEnergyChann28, &ExcitSpinChann28)); // O15
442  theChannels->insert(new G4EvaporationChannel(120, 16, 8, &ExcitEnergyChann29, &ExcitSpinChann29)); // O16
443  theChannels->insert(new G4EvaporationChannel(458, 17, 8, &ExcitEnergyChann30, &ExcitSpinChann30)); // O17
444  theChannels->insert(new G4EvaporationChannel(590, 18, 8, &ExcitEnergyChann31, &ExcitSpinChann31)); // O18
445  theChannels->insert(new G4CompetitiveFission()); // Fission Channel
446  theChannels->insert(new G4PhotonEvaporation()); // Photon Channel
447
448  G4int nChannels = theChannels->entries();
449
450  G4PhotonEvaporation* photonEvaporation = new G4PhotonEvaporation;
451
452  G4int i;
453  for (i=0; i<iter; i++)
454    {
455      G4double excitation = excMin + G4UniformRand() * (excMax - excMin);
456
457      G4cout << G4endl << "TEST >>>>>>>>> Iteration " << i
458                 << " <<<<<<<<< Initial excitation " << excitation << G4endl;
459         
460      G4LorentzVector p4(0.,0.,0.,G4NucleiProperties::GetAtomicMass(A,Z)+excitation);
461      G4Fragment nucleus(A,Z,p4);
462
463      ntuple->column("exc",excitation);
464
465      // Photon evaporation probability
466      photonEvaporation->Initialize(nucleus);
467      G4FragmentVector* products = photonEvaporation->BreakUp(nucleus);
468      G4double probPhoton = photonEvaporation->GetEmissionProbability();
469      G4cout << "TEST: Photon probability = " << probPhoton << G4endl;
470      ntuple->column("probg",probPhoton);
471
472      // Probability for each channel
473
474      G4int n;
475      for (n=0; n<=nProbs; n++)
476        {
477          G4Fragment initialNucleus(A,Z,p4);
478          theChannels->at(n)->Initialize(initialNucleus);
479          G4double prob = theChannels->at(n)->GetEmissionProbability();
480          G4cout << "TEST: probability(" << n << ") = " << prob << G4endl;
481          ntuple->column("prob"+HepString(n),prob);
482        }
483
484
485      // Fill histograms
486
487      G4int nProducts = 0;
488      if (products !=0) nProducts = products->entries();
489
490      G4cout << "TEST: " << nProducts << " products generated: gammaE ";
491     
492      G4int ig = 0;
493      for (ig=0; ig<nProducts; ig++)
494        {
495          G4double productE = products->at(ig)->GetMomentum().e();
496          G4ThreeVector pProd(products->at(ig)->GetMomentum());
497          if (products->at(ig)->GetA() < 1)
498            {
499              G4cout << productE << " ";
500              hGammaE->accumulate(productE);
501              ntuple->column("egamma",productE);
502            }
503          else
504            { G4cout << "(" << products->at(ig)->GetZ() << "," << products->at(ig)->GetA() << ") "; }
505        }
506      G4cout << G4endl;
507     
508      ntuple->dumpData();
509     
510      delete products;
511
512    }
513
514  hbookManager->write();
515
516  theChannels->clearAndDestroy();
517  delete theChannels;
518  delete photonEvaporation;
519 
520  return EXIT_SUCCESS;
521}
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
Note: See TracBrowser for help on using the repository browser.