Source

mana-core-athenaservices / src / AtRndmGenSvc.cxx

The default branch has multiple heads

Full commit
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/ServiceHandle.h"

#include "CLHEP/Random/RanecuEngine.h"
#include "CLHEP/Random/RandGauss.h"

#include "EventInfo/EventIncident.h"
#include "EventInfo/EventInfo.h"
#include "EventInfo/EventID.h"

#include "interpretSeeds.h"
#include "AtRndmGenSvc.h"
#include "crc_combine.h"

#include <cassert>
#include <iostream>

/*FIXME temporarily needed for old seeding scheme*/
#include "StoreGate/tools/hash_functions.h" 

using namespace std;

/// Standard Constructor
AtRndmGenSvc::AtRndmGenSvc(const std::string& name,ISvcLocator* svc)
  : AthService(name,svc), 
    m_read_from_file(false),
    m_file_to_read(name + ".out"),    
    m_save_to_file(true),
    m_file_to_write(name + ".out"),
    m_useOldBrokenSeeding(false),
    m_eventReseed(true),
    m_reseedStreamNames(),
    m_reseedingOffsets(),
    m_engines(), m_engines_copy()
{
    // Get user's input
    declareProperty("Seeds", m_streams_seeds,
		    "seeds for the engines, this is a vector of strings of the form ['SequenceName [OFFSET num] Seed1 Seed2', ...] where OFFSET is an optional integer that allows to change the sequence of randoms for a given run/event no and SequenceName combination. Notice that Seed1/Seed2 are dummy when EventReseeding is used");
    declareProperty("ReadFromFile", m_read_from_file,
		    "set/restore the status of the engine from file");
    declareProperty("FileToRead",   m_file_to_read,
		    "name of a ASCII file, usually produced by AtRndmGenSvc itself at the end of a job, containing the information to fully set/restore the status of Ranecu");
    declareProperty("SaveToFile", m_save_to_file,
		    "save the status of the engine to file");
    declareProperty("FileToWrite",   m_file_to_write,
		    "name of an ASCII file which will be produced on finalize, containing the information to fully set/restore the status of Ranecu");
    declareProperty("UseOldBrokenSeeding",   m_useOldBrokenSeeding,
		    "use old seeding mechanism. This is broken in that the same seeds generate different sequences for 32 and 64 bit architectures");
    declareProperty("EventReseeding", m_eventReseed, "reseed every event using a hash of run and event numbers");
    declareProperty("ReseedStreamNames", m_reseedStreamNames, "the streams we are going to set the seeds of (default: all streams)");
    

    // Set Default values
    m_default_seed1		=	3591;
    m_default_seed2		=	2309736;
    m_PYTHIA_default_seed1	=	93453591;
    m_PYTHIA_default_seed2	=	73436;
    m_HERWIG_default_seed1	=	355391;
    m_HERWIG_default_seed2	=	97336;
}


/// Standard Destructor
AtRndmGenSvc::~AtRndmGenSvc()  
{
  engineIter i(m_engines.begin()), e(m_engines.end());
  while (i != e) delete (i++)->second;
}

// Query the interfaces.
//   Input: riid, Requested interface ID
//          ppvInterface, Pointer to requested interface
//   Return: StatusCode indicating SUCCESS or FAILURE.
// N.B. Don't forget to release the interface after use!!!
StatusCode 
AtRndmGenSvc::queryInterface(const InterfaceID& riid, void** ppvInterface) 
{
    if ( IAtRndmGenSvc::interfaceID().versionMatch(riid) )    {
        *ppvInterface = (IAtRndmGenSvc*)this;
    }
    else  {
	// Interface is not directly available: try out a base class
	return AthService::queryInterface(riid, ppvInterface);
    }
    addRef();
    return StatusCode::SUCCESS;
}

StatusCode 
AtRndmGenSvc::initialize()
{
  ATH_MSG_INFO
    ("Initializing " << name()
     << " - package version " << PACKAGE_VERSION 
     << "\n INITIALISING RANDOM NUMBER STREAMS. ");


  /// Incident Service
  ServiceHandle<IIncidentSvc> pIncSvc("IncidentSvc", name());

  // set up the incident service:
  if (!(pIncSvc.retrieve()).isSuccess()) {
    ATH_MSG_ERROR ("Could not locate IncidentSvc ");
    return StatusCode::FAILURE;
  }
  
  //start listening to "EndEvent"
  static const int PRIORITY = 100;
  pIncSvc->addListener(this, "EndEvent", PRIORITY);

  //and to BeginEvent if we are reseeding
  if (m_eventReseed) {
      ATH_MSG_INFO ("will be reseeded for every event");
      pIncSvc->addListener(this, "BeginEvent");
  }
  pIncSvc.release().ignore();

  if (m_read_from_file)
    {
      // Read from a file
      ifstream	infile( m_file_to_read.c_str() );
      if ( !infile ) {
	ATH_MSG_ERROR (" Unable to open: " << m_file_to_read);
	return StatusCode::FAILURE;
      } else {
	std::string buffer;
	while (std::getline(infile, buffer)) {
	  string stream; 
	  uint32_t seed1, seed2;
	  if (interpretSeeds(buffer, stream, seed1, seed2)) {
	    ATH_MSG_DEBUG 
	      (" INITIALISING " << stream << " stream with seeds "
	       << seed1 << "  " << seed2 
	       << " read from file " << m_file_to_read);
	    CreateStream(seed1, seed2, stream);
	  } else {		
		ATH_MSG_ERROR
		  ("bad line\n" << buffer 
		   << "\n in input file " << m_file_to_read);
		return StatusCode::FAILURE;
	  }		
	}
	
      }
    }

    // Create the various streams according to user's request
    for (VStrings::const_iterator i = m_streams_seeds.begin(); i != m_streams_seeds.end(); ++i) {
      string stream; 
      uint32_t seed1, seed2, offset(0);
      short dummy;
      if (interpretSeeds(*i, stream, seed1, seed2, dummy, offset)) {
	ATH_MSG_VERBOSE
	  ("Seeds property: stream " << stream 
	   << " seeds " << seed1 << ' ' << seed2
	   << ", reseeding offset " << offset);
      } else {
	ATH_MSG_ERROR ("bad Seeds property\n" << (*i));
	return StatusCode::FAILURE;
      }		
      
      // Check if stream already generated (e.g. from reading a file)
      bool				not_found	=	true;
      if ( number_of_streams() != 0 ) {
	AtRndmGenSvc::engineConstIter sf 		= 	begin();
	do {
	  if ((*sf).first == stream) not_found = false;
	  ++sf;
	} while (sf != end() && not_found);
      }
	
      if (not_found) {
	ATH_MSG_DEBUG
	  (" INITIALISING " << stream << " stream with seeds "
	   << seed1 << "  " << seed2);
	CreateStream(seed1, seed2, stream);
	if (m_eventReseed) {
	  m_reseedingOffsets.insert(std::make_pair(stream, offset));
	  // apply the offset we just inserted
	  ATH_MSG_DEBUG("Applying reseeding  offset " << offset << 
			" to stream " << stream);
	  this->setOnDefinedSeeds(seed1, seed2, stream);
	}
      }
	
    }
    return StatusCode::SUCCESS;
}

void 
AtRndmGenSvc::handle(const Incident &inc) {

  if ( inc.type() == "EndEvent") 
  {
    ATH_MSG_DEBUG (" Handle EndEvent ");
    m_engines_copy.clear();
    for (engineConstIter i = begin(); i != end(); ++i)
    {
      HepRandomEngine* engine = GetEngine((*i).first);
      const long*	    s =	engine->getSeeds();
      std::vector<long int> tseeds;
      tseeds.push_back(s[0]);
      tseeds.push_back(s[1]);
      m_engines_copy.insert( 
	std::map<std::string, std::vector<long int> >::value_type( (*i).first,
						                   tseeds ) );
    }
    print();    
  } else if (inc.type() == "BeginEvent") {
    ATH_MSG_DEBUG (" Handle BeginEvent ");
    EventID* pei((dynamic_cast<const EventIncident&>(inc)).eventInfo().event_ID());
    //clear static RandGauss cache (generates two numbers per call to shoot()
    RandGauss::setFlag(false);
    //loop over generator streams, combining the stream name to the hash
    vector<string>::const_iterator i(m_reseedStreamNames.begin());
    vector<string>::const_iterator e(m_reseedStreamNames.end());
    //by default (when no streams are specified in streamNames, seed all
    //streams
    if (i == e) {
      if (!(this->setAllOnDefinedSeeds(pei->event_number(), 
				       pei->run_number())))
	throw GaudiException("can not reseed all streams ", name(), StatusCode::FAILURE);
    } else {
      while (i != e) {
	const string& strName(*i++);
	if (0 == this->setOnDefinedSeeds(pei->event_number(), 
					 pei->run_number(),
					 strName)) { 
	  throw GaudiException(string("can not reseed stream ") + strName,  
			       name(), StatusCode::FAILURE);
	} else {
	  msg() << MSG::VERBOSE << "Reseeded stream " << strName 
		<< " for random service " << endmsg;
	}
      }
    }
  }
}

StatusCode 
AtRndmGenSvc::finalize()
{
  ATH_MSG_INFO (" FINALISING ");

  if (m_save_to_file) {
    // Write the status of the Service into a file
    std::ofstream outfile( m_file_to_write.c_str() );
    if ( !outfile ) {
      ATH_MSG_ERROR ("error: unable to open: " << m_file_to_write);
    } else {
      for  (std::map<std::string, std::vector<long int> >::const_iterator i = m_engines_copy.begin();
	    i != m_engines_copy.end();
	    ++i) { 
	outfile << (*i).first << " " << (*i).second[0] << " " << (*i).second[1] << "\n";
      }
    ATH_MSG_DEBUG (" wrote seeds to " << m_file_to_write );
    }
  }
  return AthService::finalize();
}

HepRandomEngine*
AtRndmGenSvc::GetEngine	( const std::string& streamName )
{
    engineConstIter citer = m_engines.find(streamName);
    if ( citer == m_engines.end() )
    {
	m_engines.insert( engineValType( streamName, new RanecuEngine() ) );
	SetStreamSeeds	( streamName );
    }
    
    engineIter	iter	=	m_engines.find(streamName);
    return	(HepRandomEngine*)(*iter).second;
}

void
AtRndmGenSvc::CreateStream(uint32_t seed1, uint32_t seed2, 
			   const std::string& streamName )
{
    engineConstIter citer = m_engines.find(streamName);
    if (citer == m_engines.end())
      m_engines.insert(engineValType(streamName, new RanecuEngine()));
    engineIter iter = m_engines.find(streamName);
    long seeds[3] = { seed1, seed2, 0 };
    ((*iter).second)->setSeeds( seeds, -1 );
}

void
AtRndmGenSvc::SetStreamSeeds	( const std::string& streamName )
{
    long seed1;
    long seed2;
    if (streamName == "PYTHIA")
    {
	seed1 = m_PYTHIA_default_seed1;
	seed2 = m_PYTHIA_default_seed2;
    }
    else if (streamName == "HERWIG")
    {
	seed1 = m_HERWIG_default_seed1;
	seed2 = m_HERWIG_default_seed2;
    }
    else
    {
	seed1 = m_default_seed1;
	seed2 = m_default_seed2;
    }
    ATH_MSG_WARNING
      (" INITIALISING " << streamName << " stream with DEFAULT seeds "
       << seed1 << "  " << seed2);
    
    long seeds[3] = { seed1, seed2, 0 };
    const long* s = seeds;
    engineIter	iter	=	m_engines.find(streamName);
    ((*iter).second)->setSeeds( s, -1 );
}

void
AtRndmGenSvc::print	( const std::string& streamName )
{
    engineConstIter citer = m_engines.find(streamName);
    if ( citer == m_engines.end() )
    {
      ATH_MSG_WARNING (" Stream =  " << streamName << " NOT FOUND");
    }
    else
    {
	const long* s	=	((*citer).second)->getSeeds();
	ATH_MSG_INFO (" Stream =  " << streamName << ", Seed1 =  "
		      << s[0] << ", Seed2 = " << s[1]);
    }
}

void
AtRndmGenSvc::print	( void )
{
    for (engineConstIter i = m_engines.begin(); i != m_engines.end(); ++i)
	print( (*i).first );
}

HepRandomEngine*
AtRndmGenSvc::setOnDefinedSeeds(uint32_t eventNumber, uint32_t runNumber,
				const std::string& streamName)
{
  if (m_useOldBrokenSeeding) 
    return oldSetOnDefinedSeeds(eventNumber, runNumber, streamName);
  //do it properly
  uint32_t theHash(eventNumber);
  map<string, uint32_t>::const_iterator citer(m_reseedingOffsets.find(streamName));
  bool hasOffset(citer != m_reseedingOffsets.end() && 0 != citer->second);
  if (hasOffset) theHash=crc_combine(theHash, citer->second);

  theHash=crc_combine(theHash, runNumber);
  ATH_MSG_VERBOSE( "Reseeding stream " << streamName 
		   << " with eventNumber " << eventNumber 
		   << " runNumber " << runNumber);
  if (hasOffset) ATH_MSG_VERBOSE("Applied offset " << citer->second);
  return this->setOnDefinedSeeds(theHash, streamName);
}

HepRandomEngine*
AtRndmGenSvc::setOnDefinedSeeds(uint32_t theSeed, 
				const std::string& streamName){
  if (m_useOldBrokenSeeding) 
    return oldSetOnDefinedSeeds(theSeed, streamName);
  //do it properly
  engineConstIter citer = m_engines.find(streamName);
  ///create engine if not found. FIXME this may not be a good idea
  if ( citer == m_engines.end() ) 
    m_engines.insert(engineValType(streamName, 
				   new RanecuEngine() ) );
  engineIter iter = m_engines.find(streamName);
  theSeed=crc_combine(theSeed, streamName);
  ATH_MSG_DEBUG("Reseeding stream " << streamName << " with " << theSeed);
  //Ranecu takes a long as seed and makes a test on the sign of the seed
  //so let's make sure that our seed is presented to Ranecu as a 32 bit
  //signed int
  ((*iter).second)->setSeed( (int32_t)theSeed, 0);
  return (HepRandomEngine*)(*iter).second;
}

HepRandomEngine*
AtRndmGenSvc::oldSetOnDefinedSeeds(uint32_t eventNumber, uint32_t runNumber,
				   const std::string& streamName)
{
  engineConstIter citer = m_engines.find(streamName); 
  if ( citer == m_engines.end() ) 
    m_engines.insert( engineValType(streamName, new RanecuEngine() ) ); 
  engineIter iter = m_engines.find(streamName); 
  int hashedStream(SG::simpleStringHash(streamName)); 
  long seeds[3] = { 1000*runNumber + hashedStream,  
		    eventNumber, 0 }; 
  assert( seeds[0] > 0 ); 
  assert( seeds[1] > 0 ); 
  const long* s = seeds; 
  ((*iter).second)->setSeeds( s, -1 ); 
  return (HepRandomEngine*)(*iter).second;
}

HepRandomEngine*
AtRndmGenSvc::oldSetOnDefinedSeeds(uint32_t theSeed, 
				const std::string& streamName){
  engineConstIter citer = m_engines.find(streamName); 
  if ( citer == m_engines.end() ) 
    m_engines.insert( engineValType(streamName, new RanecuEngine() ) ); 
  engineIter iter = m_engines.find(streamName); 
  int hashedStream(SG::simpleStringHash(streamName)); 
  long seeds[3] = { hashedStream % (theSeed+13), theSeed, 0 }; 
  assert( seeds[0] > 0 ); 
  assert( seeds[1] > 0 ); 
  const long* s = seeds; 
  ((*iter).second)->setSeeds( s, -1 ); 
  return (HepRandomEngine*)(*iter).second; 
}

bool
AtRndmGenSvc::setAllOnDefinedSeeds(uint32_t eventNumber, uint32_t runNumber)
{
  bool allOK(true);
  engineIter i(m_engines.begin()), e(m_engines.end());
  while (i!=e && 
         (allOK=(0 != this->setOnDefinedSeeds(eventNumber,
                                              runNumber,
                                              (*i++).first)))) {
    /*empty*/
  }
  return allOK;
}

bool
AtRndmGenSvc::setAllOnDefinedSeeds(uint32_t theSeed){
  bool allOK(true);
  engineIter i(m_engines.begin()), e(m_engines.end());
  while (i!=e && 
         (allOK=(0 != this->setOnDefinedSeeds(theSeed, (*i++).first)))) {
    /* empty */
    }
  return allOK;
}