Source

mana-core-athenaservices / src / AtDSFMTGenSvc.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
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/ServiceHandle.h"

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

#include "interpretSeeds.h"
#include "AtDSFMTGenSvc.h"
#include "crc_combine.h"

#include "AtlasCLHEP_RandomGenerators/dSFMTEngine.h"

#include <cassert>
#include <iostream>

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

using namespace std;

/// Standard Constructor
AtDSFMTGenSvc::AtDSFMTGenSvc(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_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 AtDSFMTGenSvc itself at the end of a job, containing the information to fully set/restore the status");
    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");
    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
AtDSFMTGenSvc::~AtDSFMTGenSvc()  
{
  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 AtDSFMTGenSvc::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 AtDSFMTGenSvc::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", PRIORITY);
  }
  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; 
        std::vector<uint32_t> seeds;
        //split the space-separated string in 3 words:
        if (interpretSeeds(buffer, stream, seeds)) {
          msg (MSG::DEBUG) 
            << " INITIALISING " << stream << " stream with seeds ";
          for (std::vector<uint32_t>::const_iterator i = seeds.begin(); i != seeds.end(); ++i){
            // The state returned is a set of 32-bit numbers.
            // On 64-bit platforms, though, the vector holds 64-bit ints.
            // For the first word, one gets garbage in the high 32 bits.
            // (This is because the crc32 routine used in clhep
            // to hash the engine names doesn't mask down to 32 bits.)
            // So mask off the garbage so that we get consistent results
            // across platforms.
            msg() << ((*i) & 0xffffffffu) << " ";
          }
          msg() << " read from file " << m_file_to_read << endreq;
          if (CreateStream(seeds, stream)) {
            msg(MSG::DEBUG)
              << stream << " stream initialized succesfully" <<endreq;
          } else {
            msg(MSG::ERROR)
              << stream << " stream FAILED to initialize" <<endreq;
            return StatusCode::FAILURE;
          }
        } else {                
          msg(MSG::ERROR)
            << "bad line\n" << buffer 
            << "\n in input file " << m_file_to_read << endreq;
          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);
    //parse the stream property string
    short ll(0); // temp copy so we don't overwrite default
    if (interpretSeeds(*i, stream, seed1, seed2, ll, offset)) {
      ATH_MSG_VERBOSE("Seeds property: stream " << stream 
                      << " seeds " << seed1 << ' ' << seed2 
                      << ", luxury level " << ll
                      << ", 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 ) {
      engineConstIter sf(begin());
      while (sf != end() && (not_found=((*sf).first != stream))) ++sf;
    }
        
    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 AtDSFMTGenSvc::handle(const Incident &inc) {
  ATH_MSG_DEBUG (" Handle EndEvent ");

  if ( inc.type() == "EndEvent") 
  {
    m_engines_copy.clear();
    engineConstIter iE(begin()), eE(end());
    while(iE != eE) {
      HepRandomEngine* engine = GetEngine(iE->first);
      std::vector<unsigned long> v = engine->put();
      std::vector<uint32_t> tseeds(v.size());
      for (size_t i=0; i<v.size(); ++i) {
        // The state returned is a set of 32-bit numbers.
        // On 64-bit platforms, though, the vector holds 64-bit ints.
        // For the first word, one gets garbage in the high 32 bits.
        // (This is because the crc32 routine used in clhep
        // to hash the engine names doesn't mask down to 32 bits.)
        // Mask off the garbage to get consistent results
        // across platforms.
        tseeds[i] = (v[i] & 0xffffffffu);
      }
      m_engines_copy.insert(std::make_pair(iE->first, tseeds));
      ++iE;
    }
    
    print();    
  } else if (inc.type() == "BeginEvent") {
    ATH_MSG_DEBUG (" Handle BeginEvent ");
    EventID* pei((dynamic_cast<const EventIncident&>(inc)).eventInfo().event_ID());
    //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 AtDSFMTGenSvc::finalize()
{
  ATH_MSG_INFO (" FINALISING ");

  if (m_save_to_file) {
    // Write the status of the Service to 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<uint32_t> >::const_iterator i = m_engines_copy.begin();
           i != m_engines_copy.end();
           ++i) {
        outfile << (*i).first << " ";
        for (std::vector<uint32_t>::const_iterator j = (*i).second.begin(); j !=(*i).second.end(); ++j){
          outfile << (*j) << " ";
        }
        outfile << endl;
      }
      ATH_MSG_DEBUG (" wrote seeds to " << m_file_to_write );
    }
  }
  return AthService::finalize();
}

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

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

bool AtDSFMTGenSvc::CreateStream(const std::vector<uint32_t>& seeds, const std::string& streamName)
{
    engineConstIter citer = m_engines.find(streamName);
    if ( citer == m_engines.end() ) m_engines.insert(engineValType(streamName, new CLHEP::dSFMTEngine() ) );
    engineIter iter = m_engines.find(streamName);
    std::vector<unsigned long> longSeeds(seeds.size());
    for (size_t i=0; i<seeds.size(); ++i) longSeeds[i]=seeds[i];
    return (((*iter).second)->getState( longSeeds ));
}

void AtDSFMTGenSvc::SetStreamSeeds( const std::string& StreamName )
{
    uint32_t seed1;
    uint32_t 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,0);
}

void AtDSFMTGenSvc::print(const std::string& StreamName )
{
    engineConstIter citer = m_engines.find(StreamName);
    if ( citer == m_engines.end() ) {
      ATH_MSG_WARNING (" Stream =  " << StreamName << " NOT FOUND");
    } else {
      std::vector<unsigned long> v = ((*citer).second)->put();
      msg(MSG::INFO) << " Stream =  " << StreamName << " ";
      for (std::vector<unsigned long>::const_iterator i = v.begin(); i != v.end(); ++i){
        // The state returned is a set of 32-bit numbers.
        // On 64-bit platforms, though, the vector holds 64-bit ints.
        // For the first word, one gets garbage in the high 32 bits.
        // (This is because the crc32 routine used in clhep
        // to hash the engine names doesn't mask down to 32 bits.)
        // So mask off the garbage so that we get consistent results
        // across platforms.
        msg() << ((*i) & 0xffffffffu) << " ";
      }
      msg() << endreq;
    }
}

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

HepRandomEngine* AtDSFMTGenSvc::setOnDefinedSeeds(uint32_t eventNumber, uint32_t runNumber,
                                                  const std::string& streamName)
{
  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* AtDSFMTGenSvc::setOnDefinedSeeds(uint32_t theSeed, 
                                                  const std::string& streamName){
  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 CLHEP::dSFMTEngine() ) );

  engineIter iter = m_engines.find(streamName);
  theSeed=crc_combine(theSeed, streamName);
  ATH_MSG_DEBUG("Reseeding stream " << streamName << " with " << theSeed);
  HepRandomEngine* eng = (*iter).second;
  eng->setSeed( theSeed, 0 );
  return (HepRandomEngine*)eng;
}

bool
AtDSFMTGenSvc::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
AtDSFMTGenSvc::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;
}