source: trunk/source/global/HEPRandom/test/G4RandGlobalTest.cc@ 1215

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

nvx fichiers dans CVS
File size: 9.5 KB
RevLine 
[1199]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: G4RandGlobalTest.cc,v 1.6 2006/06/29 19:00:55 gunter Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// ----------------------------------------------------------------------
32#include "Randomize.hh"
33#include "G4ios.hh"
34
35RandEngine theRandEngine;
36DRand48Engine theDRand48Engine;
37RanluxEngine theRanluxEngine(19780503,4);
38RanecuEngine theRanecuEngine;
39
40//======================== CHI-SQUARED TEST ============================
41
42int box[10], chi[8], more=1;
43int prob[8] = { 1, 5, 25, 50, 75, 95, 99, 100 };
44float perc[8] = { 2.09, 3.33, 5.90, 8.34, 11.39, 16.92, 21.67, 100.0 };
45
46void even_test(int seq)
47{
48 int i,j;
49 char Pause;
50
51 for (i=0; i<8; ++i)
52 chi[i]=0;
53
54 G4cout << "\t\t\t --- CHI-SQUARED TEST ---" << G4endl;
55 int fac = int(seq/100);
56
57 G4cout << G4endl;
58
59 for (int k=0; k<seq; ++k) {
60 for (i=0; i<10; ++i)
61 box[i]=0;
62 for (i=0; i<1000; ++i) {
63 int h = int(100*G4UniformRand());
64 for (j=0; j<10; ++j)
65 if ((10*j <= h) && (h < 10*(j+1))) ++box[j];
66 }
67 float X=0;
68 for (i=0; i<10; ++i) {
69 float t=float(box[i]-100);
70 X += t*t;
71 }
72 X /= 100;
73 if (k == (seq%100)) {
74 G4cout << "\tDistribution of a single run (rand_valuex100) ..."
75 << G4endl << G4endl;
76 G4cout << "\t\t 1 10 20 30 40 50 60 70 80 90 100" << G4endl;
77 G4cout << "\t\t";
78 for (j=0; j<10; ++j)
79 G4cout << box[j] << " ";
80 G4cout << G4endl;
81 G4cout << "\t\t\t Chi-squared = " << X << G4endl << G4endl;
82 }
83 for (i=0; i<8; ++i)
84 if (X <= perc[i]) ++chi[i];
85 }
86 G4cout << "\tDistribution of Chi-squared for " << seq << " sequences ..."
87 << G4endl << G4endl;
88 G4cout << "\t\t\t %" << "\tChi-sq " << "\tExpected" << G4endl;
89 for (j=0; j<8; ++j)
90 G4cout << "\t\t\tX <= " << perc[j] << ":\t" << chi[j]/fac
91 << "\t" << prob[j] << G4endl;
92 G4cout << G4endl;
93 G4cout << G4endl;
94 G4cout << " ----- Press <ENTER> to continue -----";
95 if ( (Pause = G4cin.get()) != '\n') exit(0);
96 if (more == 1)
97 if ( (Pause = G4cin.get()) != '\n') exit(0);
98 G4cout << G4endl;
99} // end even_test()
100
101//====================== SERIAL CORRELATION TEST =======================
102
103void corr_test(int seq)
104{
105 int i,j;
106 char Pause;
107 double U;
108 double UU;
109 double UV;
110 double corr = 0.;
111 double stde = 0.;
112 double mv = -1./999;
113 double std = (-mv)*std::sqrt(double(1000*(997)/(1001)));
114 double next, uni, w, mve, temp[10000];
115
116 G4cout << "\t\t\t--- SERIAL CORRELATION TEST ---" << G4endl << G4endl;
117
118 for (j=0; j<seq; ++j) {
119 U = 0.; UU = 0.; UV = 0.;
120 uni = G4UniformRand();
121 w = uni;
122 for (i=1; i<1000; ++i) {
123 next = G4UniformRand();
124 U += next;
125 UU += next*next;
126 UV += uni*next;
127 uni = next;
128 }
129 U += w;
130 UU += w*w;
131 UV += uni*w;
132 temp[j] = (1000*UV-U*U)/(1000*UU-U*U);
133 corr += temp[j];
134 }
135
136 mve = corr/seq;
137 for (j=0; j<seq; ++j) {
138 temp[j] -= mve;
139 stde += temp[j]*temp[j];
140 }
141 stde = std::sqrt((1./(seq-1))*stde);
142
143 G4cout << "Mean value predicted : " << mv << "\t"
144 << "Standard deviation predicted: " << std << G4endl
145 << "Mean value effective : " << mve << "\t"
146 << "Standard deviation effective: " << stde << G4endl << G4endl;
147 G4cout << " ----- Press <ENTER> to continue -----";
148 if ( (Pause = G4cin.get()) != '\n') exit(0);
149 G4cout << G4endl;
150
151} // end corr_test
152
153//=========================== SPATIAL TEST =============================
154
155#define PI 3.1415926
156
157static unsigned long twotoj[16] = { 0x1L,0x2L,0x4L,0x8L,0x10L,0x20L,0x40L,
158 0x80L,0x100L,0x200L,0x400L,0x800L,
159 0x1000L,0x2000L,0x4000L,0x8000L };
160
161float fnc(float x1, float x2, float x3, float x4)
162{
163 return (float) std::sqrt(x1*x1+x2*x2+x3*x3+x4*x4);
164}
165
166void spat_test ()
167{
168 char Pause;
169 long iy[4],jpower,k;
170 int i,j;
171 float x1,x2,x3,x4,yprob[4];
172
173 G4cout << "\t\t\t --- SPATIAL TEST ---" << G4endl
174 << "\tCalculates PI statistically using volume of unit n-sphere "
175 << "n = 2,3,4"
176 << G4endl << G4endl;
177 for (i=1; i<=3; i++) iy[i]=0;
178 G4cout << "\t\t\t# pts\tPI\t(4/3)PI\t(1/2)PI^2" << G4endl << G4endl;
179 for (j=1; j<=15; j++) {
180 for (k=twotoj[j-1]; k>=0; k--) {
181 x1 = (float) G4UniformRand();
182 x2 = (float) G4UniformRand();
183 x3 = (float) G4UniformRand();
184 x4 = (float) G4UniformRand();
185 if (fnc(x1,x2,0.,0.) < 1.) ++iy[1];
186 if (fnc(x1,x2,x3,0.) < 1.) ++iy[2];
187 if (fnc(x1,x2,x3,x4) < 1.) ++iy[3];
188 }
189 jpower = twotoj[j];
190 for (i=1; i<=3; i++)
191 yprob[i] = (float) twotoj[i+1]*iy[i]/jpower;
192 if ((jpower <= 32768) && (jpower != 0))
193 G4cout << "\t\t\t" << jpower << "\t" << yprob[1] << "\t"
194 << yprob[2] << "\t" << yprob[3] << G4endl;
195 }
196 G4cout << G4endl;
197 G4cout << "\t\t\tactual" << "\t" << PI << "\t"
198 << 4.*PI/3. << "\t" << .5*PI*PI << G4endl << G4endl;
199
200 if (more != 99) {
201 G4cout << " ----- Press <ENTER> to continue -----";
202 if ( (Pause = G4cin.get()) != '\n') exit(0);
203 G4cout << G4endl;
204 }
205} // end spat_test()
206
207//============================= GLOBAL =================================
208
209void init()
210{
211 G4cout << G4endl << G4endl;
212 G4cout << "-------------------------- Random distribution test ---------------------------" << G4endl;
213 G4cout << " ------------------------ " << G4endl;
214 G4cout << " >>> Random Engines available <<< " << G4endl << G4endl;
215 G4cout << " > HepJamesRandom (default)" << G4endl;
216 G4cout << " > Rand" << G4endl;
217 G4cout << " > DRand48" << G4endl;
218 G4cout << " > Ranlux" << G4endl;
219 G4cout << " > Ranecu" << G4endl << G4endl;
220 G4cout << " >>> Tests performed <<< " << G4endl << G4endl;
221 G4cout << " > Even distribution test" << G4endl;
222 G4cout << " > Serial correlation test" << G4endl;
223 G4cout << " > Spatial test" << G4endl;
224 G4cout << G4endl << G4endl;
225
226} // end init()
227
228
229void layout(int seq)
230{
231 even_test(seq);
232 corr_test(seq);
233 spat_test();
234} // end layout()
235
236
237void start_test()
238{
239 char sel;
240 int seq;
241
242 G4cout << " Select the number of random values sequences:" << G4endl;
243 G4cout << "\t a - 100 sequences of 1000 random values" << G4endl;
244 G4cout << "\t b - 1000 sequences of 1000 random values" << G4endl;
245 G4cout << "\t c - 10000 sequences of 1000 random values" << G4endl;
246 G4cout << " > ";
247 G4cin >> sel;
248 if ((sel!='a')&&(sel!='b')&&(sel!='c'))
249 exit(0);
250 switch (sel) {
251 case 'a':
252 seq = 100;
253 break;
254 case 'b':
255 seq = 1000;
256 break;
257 case 'c':
258 seq = 10000;
259 break;
260 default:
261 seq = 100;
262 break;
263 }
264
265 G4cout << G4endl << G4endl;
266 G4cout << "------------------------- Test on HepJamesRandom ----------------------------" << G4endl;
267 G4cout << G4endl;
268 layout(seq);
269 G4cout << G4endl << G4endl;
270 more = 0;
271 G4cout << "--------------------------- Test on RandEngine ------------------------------" << G4endl;
272 G4cout << G4endl;
273 HepRandom::setTheEngine(&theRandEngine);
274 layout(seq);
275 G4cout << G4endl << G4endl;
276 G4cout << "------------------------- Test on DRand48Engine -----------------------------" << G4endl;
277 G4cout << G4endl;
278 HepRandom::setTheEngine(&theDRand48Engine);
279 layout(seq);
280 G4cout << G4endl << G4endl;
281 G4cout << "--------------------- Test on RanluxEngine (luxury 4) ------------------------" << G4endl;
282 G4cout << G4endl;
283 HepRandom::setTheEngine(&theRanluxEngine);
284 layout(seq);
285 G4cout << G4endl << G4endl;
286 more = 99;
287 G4cout << "-------------------------- Test on RanecuEngine ------------------------------" << G4endl;
288 G4cout << G4endl;
289 HepRandom::setTheEngine(&theRanecuEngine);
290 layout(seq);
291} // end start_test()
292
293
294int main() {
295
296 init();
297 start_test();
298
299 return 0;
300}
301
Note: See TracBrowser for help on using the repository browser.