Geant4  10.00.p03
RanecuEngine.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RanecuEngine ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 // This file is part of Geant4 (simulation toolkit for HEP).
10 //
11 // RANECU Random Engine - algorithm originally written in FORTRAN77
12 // as part of the MATHLIB HEP library.
13 
14 // =======================================================================
15 // Gabriele Cosmo - Created - 2nd February 1996
16 // - Minor corrections: 31st October 1996
17 // - Added methods for engine status: 19th November 1996
18 // - Added abs for setting seed index: 11th July 1997
19 // - Modified setSeeds() to handle default index: 16th Oct 1997
20 // - setSeed() now resets the engine status to the original
21 // values in the static table of HepRandom: 19th Mar 1998
22 // J.Marraffino - Added stream operators and related constructor.
23 // Added automatic seed selection from seed table and
24 // engine counter: 16th Feb 1998
25 // Ken Smith - Added conversion operators: 6th Aug 1998
26 // J. Marraffino - Remove dependence on hepString class 13 May 1999
27 // M. Fischler - Add endl to the end of saveStatus 10 Apr 2001
28 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
29 // M. Fischler - Methods for distrib. instance save/restore 12/8/04
30 // M. Fischler - split get() into tag validation and
31 // getState() for anonymous restores 12/27/04
32 // M. Fischler - put/get for vectors of ulongs 3/14/05
33 // M. Fischler - State-saving using only ints, for portability 4/12/05
34 // M. Fischler - Modify ctor and setSeed to utilize all info provided
35 // and avoid coincidence of same state from different
36 // seeds 6/22/10
37 //
38 // =======================================================================
39 
40 #include "CLHEP/Random/Random.h"
41 #include "CLHEP/Random/RanecuEngine.h"
42 #include "CLHEP/Random/engineIDulong.h"
43 #include <string.h> // for strcmp
44 #include <cmath>
45 #include <cstdlib>
46 
47 namespace CLHEP {
48 
49 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
50 
51 static const double prec = 4.6566128E-10;
52 
53 std::string RanecuEngine::name() const {return "RanecuEngine";}
54 
55 void RanecuEngine::further_randomize (int seq1, int col, int index, int modulus)
56 {
57  table[seq1][col] -= (index&0x3FFFFFFF);
58  while (table[seq1][col] <= 0) table[seq1][col] += (modulus-1);
59 } // mf 6/22/10
60 
61 // Number of instances with automatic seed selection
62 int RanecuEngine::numEngines = 0;
63 
64 RanecuEngine::RanecuEngine()
65 : HepRandomEngine()
66 {
67  int cycle = std::abs(int(numEngines/maxSeq));
68  seq = std::abs(int(numEngines%maxSeq));
69  numEngines += 1;
70  theSeed = seq;
71  long mask = ((cycle & 0x007fffff) << 8);
72  for (int i=0; i<2; ++i) {
73  for (int j=0; j<maxSeq; ++j) {
74  HepRandom::getTheTableSeeds(table[j],j);
75  table[j][i] ^= mask;
76  }
77  }
78  theSeeds = &table[seq][0];
79 }
80 
81 RanecuEngine::RanecuEngine(int index)
82 : HepRandomEngine()
83 {
84  int cycle = std::abs(int(index/maxSeq));
85  seq = std::abs(int(index%maxSeq));
86  theSeed = seq;
87  long mask = ((cycle & 0x000007ff) << 20);
88  for (int j=0; j<maxSeq; ++j) {
89  HepRandom::getTheTableSeeds(table[j],j);
90  table[j][0] ^= mask;
91  table[j][1] ^= mask;
92  }
93  theSeeds = &table[seq][0];
94  further_randomize (seq, 0, index, shift1); // mf 6/22/10
95 }
96 
97 RanecuEngine::RanecuEngine(std::istream& is)
98 : HepRandomEngine()
99 {
100  is >> *this;
101 }
102 
103 RanecuEngine::~RanecuEngine() {}
104 
105 void RanecuEngine::setSeed(long index, int dum)
106 {
107  seq = std::abs(int(index%maxSeq));
108  theSeed = seq;
109  HepRandom::getTheTableSeeds(table[seq],seq);
110  theSeeds = &table[seq][0];
111  further_randomize (seq, 0, index, shift1); // mf 6/22/10
112  further_randomize (seq, 1, dum, shift2); // mf 6/22/10
113 }
114 
115 void RanecuEngine::setSeeds(const long* seeds, int pos)
116 {
117  if (pos != -1) {
118  seq = std::abs(int(pos%maxSeq));
119  theSeed = seq;
120  }
121  // only positive seeds are allowed
122  table[seq][0] = std::abs(seeds[0])%shift1;
123  table[seq][1] = std::abs(seeds[1])%shift2;
124  theSeeds = &table[seq][0];
125 }
126 
127 void RanecuEngine::setIndex(long index)
128 {
129  seq = std::abs(int(index%maxSeq));
130  theSeed = seq;
131  theSeeds = &table[seq][0];
132 }
133 
134 void RanecuEngine::saveStatus( const char filename[] ) const
135 {
136  std::ofstream outFile( filename, std::ios::out ) ;
137 
138  if (!outFile.bad()) {
139  outFile << "Uvec\n";
140  std::vector<unsigned long> v = put();
141  for (unsigned int i=0; i<v.size(); ++i) {
142  outFile << v[i] << "\n";
143  }
144  }
145 }
146 
147 void RanecuEngine::restoreStatus( const char filename[] )
148 {
149  std::ifstream inFile( filename, std::ios::in);
150  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
151  std::cerr << " -- Engine state remains unchanged\n";
152  return;
153  }
154  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
155  std::vector<unsigned long> v;
156  unsigned long xin;
157  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
158  inFile >> xin;
159  if (!inFile) {
160  inFile.clear(std::ios::badbit | inFile.rdstate());
161  std::cerr << "\nJamesRandom state (vector) description improper."
162  << "\nrestoreStatus has failed."
163  << "\nInput stream is probably mispositioned now." << std::endl;
164  return;
165  }
166  v.push_back(xin);
167  }
168  getState(v);
169  return;
170  }
171 
172  if (!inFile.bad() && !inFile.eof()) {
173 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
174  for (int i=0; i<2; ++i)
175  inFile >> table[theSeed][i];
176  seq = int(theSeed);
177  }
178 }
179 
180 void RanecuEngine::showStatus() const
181 {
182  std::cout << std::endl;
183  std::cout << "--------- Ranecu engine status ---------" << std::endl;
184  std::cout << " Initial seed (index) = " << theSeed << std::endl;
185  std::cout << " Current couple of seeds = "
186  << table[theSeed][0] << ", "
187  << table[theSeed][1] << std::endl;
188  std::cout << "----------------------------------------" << std::endl;
189 }
190 
191 double RanecuEngine::flat()
192 {
193  const int index = seq;
194  long seed1 = table[index][0];
195  long seed2 = table[index][1];
196 
197  int k1 = (int)(seed1/ecuyer_b);
198  int k2 = (int)(seed2/ecuyer_e);
199 
200  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
201  if (seed1 < 0) seed1 += shift1;
202  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
203  if (seed2 < 0) seed2 += shift2;
204 
205  table[index][0] = seed1;
206  table[index][1] = seed2;
207 
208  long diff = seed1-seed2;
209 
210  if (diff <= 0) diff += (shift1-1);
211  return (double)(diff*prec);
212 }
213 
214 void RanecuEngine::flatArray(const int size, double* vect)
215 {
216  const int index = seq;
217  long seed1 = table[index][0];
218  long seed2 = table[index][1];
219  int k1, k2;
220  int i;
221 
222  for (i=0; i<size; ++i)
223  {
224  k1 = (int)(seed1/ecuyer_b);
225  k2 = (int)(seed2/ecuyer_e);
226 
227  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
228  if (seed1 < 0) seed1 += shift1;
229  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
230  if (seed2 < 0) seed2 += shift2;
231 
232  long diff = seed1-seed2;
233  if (diff <= 0) diff += (shift1-1);
234 
235  vect[i] = (double)(diff*prec);
236  }
237  table[index][0] = seed1;
238  table[index][1] = seed2;
239 }
240 
241 RanecuEngine::operator unsigned int() {
242  const int index = seq;
243  long seed1 = table[index][0];
244  long seed2 = table[index][1];
245 
246  int k1 = (int)(seed1/ecuyer_b);
247  int k2 = (int)(seed2/ecuyer_e);
248 
249  seed1 = ecuyer_a*(seed1-k1*ecuyer_b)-k1*ecuyer_c;
250  if (seed1 < 0) seed1 += shift1;
251  seed2 = ecuyer_d*(seed2-k2*ecuyer_e)-k2*ecuyer_f;
252  if (seed2 < 0) seed2 += shift2;
253 
254  table[index][0] = seed1;
255  table[index][1] = seed2;
256  long diff = seed1-seed2;
257  if( diff <= 0 ) diff += (shift1-1);
258 
259  return ((diff << 1) | (seed1&1))& 0xffffffff;
260 }
261 
262 std::ostream & RanecuEngine::put( std::ostream& os ) const
263 {
264  char beginMarker[] = "RanecuEngine-begin";
265  os << beginMarker << "\nUvec\n";
266  std::vector<unsigned long> v = put();
267  for (unsigned int i=0; i<v.size(); ++i) {
268  os << v[i] << "\n";
269  }
270  return os;
271 }
272 
273 std::vector<unsigned long> RanecuEngine::put () const {
274  std::vector<unsigned long> v;
275  v.push_back (engineIDulong<RanecuEngine>());
276  v.push_back(static_cast<unsigned long>(theSeed));
277  v.push_back(static_cast<unsigned long>(table[theSeed][0]));
278  v.push_back(static_cast<unsigned long>(table[theSeed][1]));
279  return v;
280 }
281 
282 std::istream & RanecuEngine::get ( std::istream& is )
283 {
284  char beginMarker [MarkerLen];
285 
286  is >> std::ws;
287  is.width(MarkerLen); // causes the next read to the char* to be <=
288  // that many bytes, INCLUDING A TERMINATION \0
289  // (Stroustrup, section 21.3.2)
290  is >> beginMarker;
291  if (strcmp(beginMarker,"RanecuEngine-begin")) {
292  is.clear(std::ios::badbit | is.rdstate());
293  std::cerr << "\nInput stream mispositioned or"
294  << "\nRanecuEngine state description missing or"
295  << "\nwrong engine type found." << std::endl;
296  return is;
297  }
298  return getState(is);
299 }
300 
301 std::string RanecuEngine::beginTag ( ) {
302  return "RanecuEngine-begin";
303 }
304 
305 std::istream & RanecuEngine::getState ( std::istream& is )
306 {
307  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
308  std::vector<unsigned long> v;
309  unsigned long uu;
310  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
311  is >> uu;
312  if (!is) {
313  is.clear(std::ios::badbit | is.rdstate());
314  std::cerr << "\nRanecuEngine state (vector) description improper."
315  << "\ngetState() has failed."
316  << "\nInput stream is probably mispositioned now." << std::endl;
317  return is;
318  }
319  v.push_back(uu);
320  }
321  getState(v);
322  return (is);
323  }
324 
325 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
326  char endMarker [MarkerLen];
327  for (int i=0; i<2; ++i) {
328  is >> table[theSeed][i];
329  }
330  is >> std::ws;
331  is.width(MarkerLen);
332  is >> endMarker;
333  if (strcmp(endMarker,"RanecuEngine-end")) {
334  is.clear(std::ios::badbit | is.rdstate());
335  std::cerr << "\nRanecuEngine state description incomplete."
336  << "\nInput stream is probably mispositioned now." << std::endl;
337  return is;
338  }
339 
340  seq = int(theSeed);
341  return is;
342 }
343 
344 bool RanecuEngine::get (const std::vector<unsigned long> & v) {
345  if ((v[0] & 0xffffffffUL) != engineIDulong<RanecuEngine>()) {
346  std::cerr <<
347  "\nRanecuEngine get:state vector has wrong ID word - state unchanged\n";
348  return false;
349  }
350  return getState(v);
351 }
352 
353 bool RanecuEngine::getState (const std::vector<unsigned long> & v) {
354  if (v.size() != VECTOR_STATE_SIZE ) {
355  std::cerr <<
356  "\nRanecuEngine get:state vector has wrong length - state unchanged\n";
357  return false;
358  }
359  theSeed = v[1];
360  table[theSeed][0] = v[2];
361  table[theSeed][1] = v[3];
362  seq = int(theSeed);
363  return true;
364 }
365 
366 
367 } // namespace CLHEP
static const int MarkerLen
Definition: DualRand.cc:60
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4String name
Definition: TRTMaterials.hh:40
static const double prec
Definition: RanecuEngine.cc:51
double flat()
Definition: G4AblaRandom.cc:47
static const G4double pos
void setSeeds(const SeedVector &sv)
Set the seeds of the current generator.
Definition: G4INCLRandom.cc:66