PhotosRandom.cxx
1#include <iostream>
2#include "PhotosRandom.h"
3#include "Photos.h"
4#include "Log.h"
5
6namespace Photospp
7{
8
9bool PhotosRandom::init = false;
10int PhotosRandom::iseed[2]= { 1802, 9373 };
11int PhotosRandom::i97 = 96;
12int PhotosRandom::j97 = 32;
13double PhotosRandom::uran[97]= { 0.0 };
14double PhotosRandom::cran = 362436.0 /16777216.0;
15const double PhotosRandom::cdran = 7654321.0 /16777216.0;
16const double PhotosRandom::cmran = 16777213.0/16777216.0;
17
18void PhotosRandom::setSeed(int s1,int s2)
19{
20 if(s1<0 || s1>31327) Log::Fatal("PhotosRandom::setSeed(): Seed(1) out of range [0,31327]",8);
21 if(s2<0 || s2>30080) Log::Fatal("PhotosRandom::setSeed(): Seed(2) out of range [0,30080]",9);
22 iseed[0]=s1;
23 iseed[1]=s2;
24}
25
26/*******************************************************************************
27 PHORIN: PHOton radiation in decays RANdom number generator init
28
29 Purpose: Initialse PHORAN with the user specified seeds in the
30 array iseed. For details see also: F. James CERN DD-
31 Report November 1988.
32
33 Author(s): B. van Eijk and F. James Created at: 27/09/89
34 Last Update: 22/02/90
35 Rewritten to C++: 18/10/10
36 by T. Przedzinski (tprzedzi@cern.ch)
37*******************************************************************************/
38void PhotosRandom::initialize()
39{
40 long IS1,IS2,IS3,IS4,IS5;
41 double S,T;
42
43// Calculate Marsaglia and Zaman seeds (by F. James)
44 IS1=(iseed[0]/177)%177+2;
45 IS2= iseed[0]%177+2;
46 IS3=(iseed[1]/169)%178+1;
47 IS4= iseed[1]%169;
48 for(int i=0;i<97;i++)
49 {
50 S=0.0;
51 T=0.5;
52 for(int j=0;j<24;j++)
53 {
54 IS5=( ((IS1*IS2)%179)*IS3 )%179;
55 IS1=IS2;
56 IS2=IS3;
57 IS3=IS5;
58 IS4=(53*IS4+1)%169;
59 if( (IS4*IS5)%64>=32) S=S+T;
60 T=0.5*T;
61 }
62 uran[i]=S;
63 }
64 init=true;
65 Log::Debug(0)<<"PhotosRandom::inititalize(): seed: "<<iseed[0]<<", "<<iseed[1]<<std::endl;
66}
67
68/*******************************************************************************
69 PHORAN: PHOton radiation in decays ret number generator based
70 on Marsaglia Algorithm
71
72 Purpose: Generate uniformly distributed random numbers between
73 0 and 1. Super long period: 2**144. See also:
74 G. Marsaglia and A. Zaman, FSU-SCR-87-50, for seed mo-
75 difications to this version see: F. James DD-Report,
76 November 1988. The generator has to be initialized by
77 a call to PHORIN ( C++ version: initialize() ).
78
79 Author(s): B. van Eijk, G. Marsaglia and Created at: 27/09/89
80 A. Zaman Last Update: 27/09/89
81 Rewritten to C++: 18/10/10
82 by T. Przedzinski (tprzedzi@cern.ch)
83*******************************************************************************/
84double PhotosRandom::randomReal()
85{
86 if(!init) Log::Fatal("PhotosRandom::randomReal(): generator not initialized",1);
87 double ret=0.0;
88 while(true)
89 {
90 ret = uran[i97]-uran[j97];
91 if(ret<0.0) ret+=1.;
92 uran[i97]=ret;
93 i97--;
94 if(i97<0) i97=96;
95 j97--;
96 if(j97<0) j97=96;
97 cran-=cdran;
98 if(cran<0.0) cran+=cmran;
99 ret-=cran;
100 if(ret<0.0) ret+=1.0;
101 if(ret>0.0) break;
102 }
103 return ret;
104}
105
106} // namespace Photospp
107
static ostream & Debug(unsigned short int code=0, bool count=true)
Definition: Log.cxx:33
static void Fatal(string text, unsigned short int code=0)