ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ranlux64Engine.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Ranlux64Engine.cc
1 // -*- C++ -*-
2 //
3 // -----------------------------------------------------------------------
4 // HEP Random
5 // --- Ranlux64Engine ---
6 // class implementation file
7 // -----------------------------------------------------------------------
8 // A double-precision implementation of the RanluxEngine generator as
9 // decsribed by the notes of the original ranlux author (Martin Luscher)
10 //
11 // See the note by Martin Luscher, December 1997, entitiled
12 // Double-precision implementation of the random number generator ranlux
13 //
14 // =======================================================================
15 // Ken Smith - Initial draft: 14th Jul 1998
16 // - Removed pow() from flat method 14th Jul 1998
17 // - Added conversion operators: 6th Aug 1998
18 //
19 // Mark Fischler The following were modified mostly to make the routine
20 // exactly match the Luscher algorithm in generating 48-bit
21 // randoms:
22 // 9/9/98 - Substantial changes in what used to be flat() to match
23 // algorithm in Luscher's ranlxd.c
24 // - Added update() method for 12 numbers, making flat() trivial
25 // - Added advance() method to hold the unrolled loop for update
26 // - Distinction between three forms of seeding such that it
27 // is impossible to get same sequence from different forms -
28 // done by discarding some fraction of one macro cycle which
29 // is different for the three cases
30 // - Change the misnomer "seed_table" to the more accurate
31 // "randoms"
32 // - Removed the no longer needed count12, i_lag, j_lag, etc.
33 // - Corrected seed procedure which had been filling bits past
34 // 2^-48. This actually was very bad, invalidating the
35 // number theory behind the proof that ranlxd is good.
36 // - Addition of 2**(-49) to generated number to prevent zero
37 // from being returned; this does not affect the sequence
38 // itself.
39 // - Corrected ecu seeding, which had been supplying only
40 // numbers less than 1/2. This is probably moot.
41 // 9/15/98 - Modified use of the various exponents of 2
42 // to avoid per-instance space overhead. Note that these
43 // are initialized in setSeed, which EVERY constructor
44 // must invoke.
45 // J. Marraffino - Remove dependence on hepString class 13 May 1999
46 // M. Fischler - In restore, checkFile for file not found 03 Dec 2004
47 // M. Fischler - put get Methods for distrib instance save/restore 12/8/04
48 // M. Fischler - split get() into tag validation and
49 // getState() for anonymous restores 12/27/04
50 // M. Fischler - put/get for vectors of ulongs 3/14/05
51 // M. Fischler - State-saving using only ints, for portability 4/12/05
52 //
53 // =======================================================================
54 
55 #include "CLHEP/Random/Random.h"
58 #include "CLHEP/Random/DoubConv.h"
60 
61 #include <string.h> // for strcmp
62 #include <cstdlib> // for std::abs(int)
63 #include <limits> // for numeric_limits
64 
65 namespace CLHEP {
66 
67 namespace {
68  // Number of instances with automatic seed selection
69  CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
70 
71  // Maximum index into the seed table
72  const int maxIndex = 215;
73 }
74 
75 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
76 
77 
78 #ifndef WIN32
79 namespace detail {
80 
81 template< std::size_t n,
82  bool = n < std::size_t(std::numeric_limits<unsigned long>::digits) >
83  struct do_right_shift;
84 template< std::size_t n >
85  struct do_right_shift<n,true>
86 {
87  unsigned long operator()(unsigned long value) { return value >> n; }
88 };
89 template< std::size_t n >
90  struct do_right_shift<n,false>
91 {
92  unsigned long operator()(unsigned long) { return 0ul; }
93 };
94 
95 template< std::size_t nbits >
96  unsigned long rshift( unsigned long value )
97 { return do_right_shift<nbits>()(value); }
98 
99 } // namespace detail
100 #endif
101 
102 std::string Ranlux64Engine::name() const {return "Ranlux64Engine";}
103 
105 : HepRandomEngine()
106 {
107  luxury = 1;
108  int numEngines = numberOfEngines++;
109  int cycle = std::abs(int(numEngines/maxIndex));
110  int curIndex = std::abs(int(numEngines%maxIndex));
111 
112  long mask = ((cycle & 0x007fffff) << 8);
113  long seedlist[2];
114  HepRandom::getTheTableSeeds( seedlist, curIndex );
115  seedlist[0] ^= mask;
116  seedlist[1] = 0;
117 
118  setSeeds(seedlist, luxury);
119  advance ( 8 ); // Discard some iterations and ensure that
120  // this sequence won't match one where seeds
121  // were provided.
122 }
123 
125 : HepRandomEngine()
126 {
127  luxury = lux;
128  long seedlist[2]={seed,0};
129  setSeeds(seedlist, lux);
130  advance ( 2*lux + 1 ); // Discard some iterations to use a different
131  // point in the sequence.
132 }
133 
134 Ranlux64Engine::Ranlux64Engine(int rowIndex, int, int lux)
135 : HepRandomEngine()
136 {
137  luxury = lux;
138  int cycle = std::abs(int(rowIndex/maxIndex));
139  int row = std::abs(int(rowIndex%maxIndex));
140  long mask = (( cycle & 0x000007ff ) << 20 );
141  long seedlist[2];
142  HepRandom::getTheTableSeeds( seedlist, row );
143  seedlist[0] ^= mask;
144  seedlist[1]= 0;
145  setSeeds(seedlist, lux);
146 }
147 
149 : HepRandomEngine()
150 {
151  is >> *this;
152 }
153 
155 
157  // Luscher improves the speed by computing several numbers in a shot,
158  // in a manner similar to that of the Tausworth in DualRand or the Hurd
159  // engines. Thus, the real work is done in update(). Here we merely ensure
160  // that zero, which the algorithm can produce, is never returned by flat().
161 
162  if (index <= 0) update();
163  return randoms[--index] + twoToMinus_49();
164 }
165 
167  // Update the stash of twelve random numbers.
168  // When this routione is entered, index is always 0. The randoms
169  // contains the last 12 numbers in the sequents: s[0] is x[a+11],
170  // s[1] is x[a+10] ... and s[11] is x[a] for some a. Carry contains
171  // the last carry value (c[a+11]).
172  //
173  // The recursion relation (3) in Luscher's note says
174  // delta[n] = x[n-s] = x[n-r] -c[n-1] or for n=a+12,
175  // delta[a+12] = x[a+7] - x[a] -c[a+11] where we use r=12, s=5 per eqn. (7)
176  // This reduces to
177  // s[11] = s[4] - s[11] - carry.
178  // The next number similarly will be given by s[10] = s[3] - s[10] - carry,
179  // and so forth until s[0] is filled.
180  //
181  // However, we need to skip 397, 202 or 109 numbers - these are not divisible
182  // by 12 - to "fare well in the spectral test".
183 
184  advance(pDozens);
185 
186  // Since we wish at the end to have the 12 last numbers in the order of
187  // s[11] first, till s[0] last, we will have to do 1, 10, or 1 iterations
188  // and then re-arrange to place to get the oldest one in s[11].
189  // Generically, this will imply re-arranging the s array at the end,
190  // but we can treat the special case of endIters = 1 separately for superior
191  // efficiency in the cases of levels 0 and 2.
192 
193  double y1;
194 
195  if ( endIters == 1 ) { // Luxury levels 0 and 2 will go here
196  y1 = randoms[ 4] - randoms[11] - carry;
197  if ( y1 < 0.0 ) {
198  y1 += 1.0;
199  carry = twoToMinus_48();
200  } else {
201  carry = 0.0;
202  }
203  randoms[11] = randoms[10];
204  randoms[10] = randoms[ 9];
205  randoms[ 9] = randoms[ 8];
206  randoms[ 8] = randoms[ 7];
207  randoms[ 7] = randoms[ 6];
208  randoms[ 6] = randoms[ 5];
209  randoms[ 5] = randoms[ 4];
210  randoms[ 4] = randoms[ 3];
211  randoms[ 3] = randoms[ 2];
212  randoms[ 2] = randoms[ 1];
213  randoms[ 1] = randoms[ 0];
214  randoms[ 0] = y1;
215 
216  } else {
217 
218  int m, nr, ns;
219  for ( m = 0, nr = 11, ns = 4; m < endIters; ++m, --nr ) {
220  y1 = randoms [ns] - randoms[nr] - carry;
221  if ( y1 < 0.0 ) {
222  y1 += 1.0;
223  carry = twoToMinus_48();
224  } else {
225  carry = 0.0;
226  }
227  randoms[nr] = y1;
228  --ns;
229  if ( ns < 0 ) {
230  ns = 11;
231  }
232  } // loop on m
233 
234  double temp[12];
235  for (m=0; m<12; m++) {
236  temp[m]=randoms[m];
237  }
238 
239  ns = 11 - endIters;
240  for (m=11; m>=0; --m) {
241  randoms[m] = temp[ns];
242  --ns;
243  if ( ns < 0 ) {
244  ns = 11;
245  }
246  }
247 
248  }
249 
250  // Now when we return, there are 12 fresh usable numbers in s[11] ... s[0]
251 
252  index = 11;
253 
254 } // update()
255 
256 void Ranlux64Engine::advance(int dozens) {
257 
258  double y1, y2, y3;
259  double cValue = twoToMinus_48();
260  double zero = 0.0;
261  double one = 1.0;
262 
263  // Technical note: We use Luscher's trick to only do the
264  // carry subtraction when we really have to. Like him, we use
265  // three registers instead of two so that we avoid sequences
266  // like storing y1 then immediately replacing its value:
267  // some architectures lose time when this is done.
268 
269  // Luscher's ranlxd.c fills the stash going
270  // upward. We fill it downward to save a bit of time in the
271  // flat() routine at no cost later. This means that while
272  // Luscher's ir is jr+5, our n-r is (n-s)-5. (Note that
273  // though ranlxd.c initializes ir and jr to 11 and 7, ir as
274  // used is 5 more than jr because update is entered after
275  // incrementing ir.)
276  //
277 
278  // I have CAREFULLY checked that the algorithms do match
279  // in all details.
280 
281  int k;
282  for ( k = dozens; k > 0; --k ) {
283 
284  y1 = randoms[ 4] - randoms[11] - carry;
285 
286  y2 = randoms[ 3] - randoms[10];
287  if ( y1 < zero ) {
288  y1 += one;
289  y2 -= cValue;
290  }
291  randoms[11] = y1;
292 
293  y3 = randoms[ 2] - randoms[ 9];
294  if ( y2 < zero ) {
295  y2 += one;
296  y3 -= cValue;
297  }
298  randoms[10] = y2;
299 
300  y1 = randoms[ 1] - randoms[ 8];
301  if ( y3 < zero ) {
302  y3 += one;
303  y1 -= cValue;
304  }
305  randoms[ 9] = y3;
306 
307  y2 = randoms[ 0] - randoms[ 7];
308  if ( y1 < zero ) {
309  y1 += one;
310  y2 -= cValue;
311  }
312  randoms[ 8] = y1;
313 
314  y3 = randoms[11] - randoms[ 6];
315  if ( y2 < zero ) {
316  y2 += one;
317  y3 -= cValue;
318  }
319  randoms[ 7] = y2;
320 
321  y1 = randoms[10] - randoms[ 5];
322  if ( y3 < zero ) {
323  y3 += one;
324  y1 -= cValue;
325  }
326  randoms[ 6] = y3;
327 
328  y2 = randoms[ 9] - randoms[ 4];
329  if ( y1 < zero ) {
330  y1 += one;
331  y2 -= cValue;
332  }
333  randoms[ 5] = y1;
334 
335  y3 = randoms[ 8] - randoms[ 3];
336  if ( y2 < zero ) {
337  y2 += one;
338  y3 -= cValue;
339  }
340  randoms[ 4] = y2;
341 
342  y1 = randoms[ 7] - randoms[ 2];
343  if ( y3 < zero ) {
344  y3 += one;
345  y1 -= cValue;
346  }
347  randoms[ 3] = y3;
348 
349  y2 = randoms[ 6] - randoms[ 1];
350  if ( y1 < zero ) {
351  y1 += one;
352  y2 -= cValue;
353  }
354  randoms[ 2] = y1;
355 
356  y3 = randoms[ 5] - randoms[ 0];
357  if ( y2 < zero ) {
358  y2 += one;
359  y3 -= cValue;
360  }
361  randoms[ 1] = y2;
362 
363  if ( y3 < zero ) {
364  y3 += one;
365  carry = cValue;
366  }
367  randoms[ 0] = y3;
368 
369  } // End of major k loop doing 12 numbers at each cycle
370 
371 } // advance(dozens)
372 
373 void Ranlux64Engine::flatArray(const int size, double* vect) {
374  for( int i=0; i < size; ++i ) {
375  vect[i] = flat();
376  }
377 }
378 
380 
381 // The initialization is carried out using a Multiplicative
382 // Congruential generator using formula constants of L'Ecuyer
383 // as described in "A review of pseudorandom number generators"
384 // (Fred James) published in Computer Physics Communications 60 (1990)
385 // pages 329-344
386 
387  const int ecuyer_a(53668);
388  const int ecuyer_b(40014);
389  const int ecuyer_c(12211);
390  const int ecuyer_d(2147483563);
391 
392  const int lux_levels[3] = {109, 202, 397};
393  theSeed = seed;
394 
395  if( (lux > 2)||(lux < 0) ){
396  pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
397  }else{
398  pDiscard = lux_levels[luxury];
399  }
400  pDozens = pDiscard / 12;
401  endIters = pDiscard % 12;
402 
403  long init_table[24];
404  long next_seed = seed;
405  long k_multiple;
406  int i;
407  next_seed &= 0xffffffff;
408  while( next_seed >= ecuyer_d ) {
409  next_seed -= ecuyer_d;
410  }
411 
412  for(i = 0;i != 24;i++){
413  k_multiple = next_seed / ecuyer_a;
414  next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
415  - k_multiple * ecuyer_c;
416  if(next_seed < 0) {
417  next_seed += ecuyer_d;
418  }
419  next_seed &= 0xffffffff;
420  init_table[i] = next_seed;
421  }
422  // are we on a 64bit machine?
423  if( sizeof(long) >= 8 ) {
424  int64_t topbits1, topbits2;
425 #ifdef WIN32
426  topbits1 = ( (int64_t) seed >> 32) & 0xffff ;
427  topbits2 = ( (int64_t) seed >> 48) & 0xffff ;
428 #else
429  topbits1 = detail::rshift<32>(seed) & 0xffff ;
430  topbits2 = detail::rshift<48>(seed) & 0xffff ;
431 #endif
432  init_table[0] ^= topbits1;
433  init_table[2] ^= topbits2;
434  //std::cout << " init_table[0] " << init_table[0] << " from " << topbits1 << std::endl;
435  //std::cout << " init_table[2] " << init_table[2] << " from " << topbits2 << std::endl;
436  }
437 
438  for(i = 0;i < 12; i++){
439  randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
440  (init_table[2*i+1] >> 15) * twoToMinus_48();
441  //if( randoms[i] < 0. || randoms[i] > 1. ) {
442  //std::cout << "setSeed: init_table " << init_table[2*i ] << std::endl;
443  //std::cout << "setSeed: init_table " << init_table[2*i+1] << std::endl;
444  //std::cout << "setSeed: random " << i << " is " << randoms[i] << std::endl;
445  //}
446  }
447 
448  carry = 0.0;
449  if ( randoms[11] == 0. ) carry = twoToMinus_48();
450  index = 11;
451 
452 } // setSeed()
453 
454 void Ranlux64Engine::setSeeds(const long * seeds, int lux) {
455 // old code only uses the first long in seeds
456 // setSeed( *seeds ? *seeds : 32767, lux );
457 // theSeeds = seeds;
458 
459 // using code from Ranlux - even those are 32bit seeds,
460 // that is good enough to completely differentiate the sequences
461 
462  const int ecuyer_a = 53668;
463  const int ecuyer_b = 40014;
464  const int ecuyer_c = 12211;
465  const int ecuyer_d = 2147483563;
466 
467  const int lux_levels[3] = {109, 202, 397};
468  const long *seedptr;
469 
470  theSeeds = seeds;
471  seedptr = seeds;
472 
473  if(seeds == 0){
474  setSeed(theSeed,lux);
475  theSeeds = &theSeed;
476  return;
477  }
478 
479  theSeed = *seeds;
480 
481 // number of additional random numbers that need to be 'thrown away'
482 // every 24 numbers is set using luxury level variable.
483 
484  if( (lux > 2)||(lux < 0) ){
485  pDiscard = (lux >= 12) ? (lux-12) : lux_levels[1];
486  }else{
487  pDiscard = lux_levels[luxury];
488  }
489  pDozens = pDiscard / 12;
490  endIters = pDiscard % 12;
491 
492  long init_table[24];
493  long next_seed = *seeds;
494  long k_multiple;
495  int i;
496 
497  for( i = 0;(i != 24)&&(*seedptr != 0);i++){
498  init_table[i] = *seedptr & 0xffffffff;
499  seedptr++;
500  }
501 
502  if(i != 24){
503  next_seed = init_table[i-1];
504  for(;i != 24;i++){
505  k_multiple = next_seed / ecuyer_a;
506  next_seed = ecuyer_b * (next_seed - k_multiple * ecuyer_a)
507  - k_multiple * ecuyer_c;
508  if(next_seed < 0) {
509  next_seed += ecuyer_d;
510  }
511  next_seed &= 0xffffffff;
512  init_table[i] = next_seed;
513  }
514  }
515 
516  for(i = 0;i < 12; i++){
517  randoms[i] = (init_table[2*i ] ) * 2.0 * twoToMinus_32() +
518  (init_table[2*i+1] >> 15) * twoToMinus_48();
519  }
520 
521  carry = 0.0;
522  if ( randoms[11] == 0. ) carry = twoToMinus_48();
523  index = 11;
524 
525 }
526 
527 void Ranlux64Engine::saveStatus( const char filename[] ) const
528 {
529  std::ofstream outFile( filename, std::ios::out ) ;
530  if (!outFile.bad()) {
531  outFile << "Uvec\n";
532  std::vector<unsigned long> v = put();
533  for (unsigned int i=0; i<v.size(); ++i) {
534  outFile << v[i] << "\n";
535  }
536  }
537 }
538 
540 {
541  std::ifstream inFile( filename, std::ios::in);
542  if (!checkFile ( inFile, filename, engineName(), "restoreStatus" )) {
543  std::cerr << " -- Engine state remains unchanged\n";
544  return;
545  }
546  if ( possibleKeywordInput ( inFile, "Uvec", theSeed ) ) {
547  std::vector<unsigned long> v;
548  unsigned long xin;
549  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
550  inFile >> xin;
551  if (!inFile) {
552  inFile.clear(std::ios::badbit | inFile.rdstate());
553  std::cerr << "\nJamesRandom state (vector) description improper."
554  << "\nrestoreStatus has failed."
555  << "\nInput stream is probably mispositioned now." << std::endl;
556  return;
557  }
558  v.push_back(xin);
559  }
560  getState(v);
561  return;
562  }
563 
564  if (!inFile.bad() && !inFile.eof()) {
565 // inFile >> theSeed; removed -- encompased by possibleKeywordInput
566  for (int i=0; i<12; ++i) {
567  inFile >> randoms[i];
568  }
569  inFile >> carry; inFile >> index;
570  inFile >> luxury; inFile >> pDiscard;
571  pDozens = pDiscard / 12;
572  endIters = pDiscard % 12;
573  }
574 }
575 
577 {
578  std::cout << std::endl;
579  std::cout << "--------- Ranlux engine status ---------" << std::endl;
580  std::cout << " Initial seed = " << theSeed << std::endl;
581  std::cout << " randoms[] = ";
582  for (int i=0; i<12; ++i) {
583  std::cout << randoms[i] << std::endl;
584  }
585  std::cout << std::endl;
586  std::cout << " carry = " << carry << ", index = " << index << std::endl;
587  std::cout << " luxury = " << luxury << " pDiscard = "
588  << pDiscard << std::endl;
589  std::cout << "----------------------------------------" << std::endl;
590 }
591 
592 std::ostream & Ranlux64Engine::put( std::ostream& os ) const
593 {
594  char beginMarker[] = "Ranlux64Engine-begin";
595  os << beginMarker << "\nUvec\n";
596  std::vector<unsigned long> v = put();
597  for (unsigned int i=0; i<v.size(); ++i) {
598  os << v[i] << "\n";
599  }
600  return os;
601 }
602 
603 std::vector<unsigned long> Ranlux64Engine::put () const {
604  std::vector<unsigned long> v;
605  v.push_back (engineIDulong<Ranlux64Engine>());
606  std::vector<unsigned long> t;
607  for (int i=0; i<12; ++i) {
609  v.push_back(t[0]); v.push_back(t[1]);
610  }
612  v.push_back(t[0]); v.push_back(t[1]);
613  v.push_back(static_cast<unsigned long>(index));
614  v.push_back(static_cast<unsigned long>(luxury));
615  v.push_back(static_cast<unsigned long>(pDiscard));
616  return v;
617 }
618 
619 std::istream & Ranlux64Engine::get ( std::istream& is )
620 {
621  char beginMarker [MarkerLen];
622  is >> std::ws;
623  is.width(MarkerLen); // causes the next read to the char* to be <=
624  // that many bytes, INCLUDING A TERMINATION \0
625  // (Stroustrup, section 21.3.2)
626  is >> beginMarker;
627  if (strcmp(beginMarker,"Ranlux64Engine-begin")) {
628  is.clear(std::ios::badbit | is.rdstate());
629  std::cerr << "\nInput stream mispositioned or"
630  << "\nRanlux64Engine state description missing or"
631  << "\nwrong engine type found." << std::endl;
632  return is;
633  }
634  return getState(is);
635 }
636 
637 std::string Ranlux64Engine::beginTag ( ) {
638  return "Ranlux64Engine-begin";
639 }
640 
641 std::istream & Ranlux64Engine::getState ( std::istream& is )
642 {
643  if ( possibleKeywordInput ( is, "Uvec", theSeed ) ) {
644  std::vector<unsigned long> v;
645  unsigned long uu;
646  for (unsigned int ivec=0; ivec < VECTOR_STATE_SIZE; ++ivec) {
647  is >> uu;
648  if (!is) {
649  is.clear(std::ios::badbit | is.rdstate());
650  std::cerr << "\nRanlux64Engine state (vector) description improper."
651  << "\ngetState() has failed."
652  << "\nInput stream is probably mispositioned now." << std::endl;
653  return is;
654  }
655  v.push_back(uu);
656  }
657  getState(v);
658  return (is);
659  }
660 
661 // is >> theSeed; Removed, encompassed by possibleKeywordInput()
662 
663  char endMarker [MarkerLen];
664  for (int i=0; i<12; ++i) {
665  is >> randoms[i];
666  }
667  is >> carry; is >> index;
668  is >> luxury; is >> pDiscard;
669  pDozens = pDiscard / 12;
670  endIters = pDiscard % 12;
671  is >> std::ws;
672  is.width(MarkerLen);
673  is >> endMarker;
674  if (strcmp(endMarker,"Ranlux64Engine-end")) {
675  is.clear(std::ios::badbit | is.rdstate());
676  std::cerr << "\nRanlux64Engine state description incomplete."
677  << "\nInput stream is probably mispositioned now." << std::endl;
678  return is;
679  }
680  return is;
681 }
682 
683 bool Ranlux64Engine::get (const std::vector<unsigned long> & v) {
684  if ((v[0] & 0xffffffffUL) != engineIDulong<Ranlux64Engine>()) {
685  std::cerr <<
686  "\nRanlux64Engine get:state vector has wrong ID word - state unchanged\n";
687  return false;
688  }
689  return getState(v);
690 }
691 
692 bool Ranlux64Engine::getState (const std::vector<unsigned long> & v) {
693  if (v.size() != VECTOR_STATE_SIZE ) {
694  std::cerr <<
695  "\nRanlux64Engine get:state vector has wrong length - state unchanged\n";
696  return false;
697  }
698  std::vector<unsigned long> t(2);
699  for (int i=0; i<12; ++i) {
700  t[0] = v[2*i+1]; t[1] = v[2*i+2];
702  }
703  t[0] = v[25]; t[1] = v[26];
705  index = v[27];
706  luxury = v[28];
707  pDiscard = v[29];
708  return true;
709 }
710 
711 } // namespace CLHEP