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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 9.5 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: G4RandGlobalTest.cc,v 1.6 2006/06/29 19:00:55 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
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.