#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" usingnamespacestd;/// Standard ConstructorAtDSFMTGenSvc::AtDSFMTGenSvc(conststd::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 inputdeclareProperty("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 valuesm_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 DestructorAtDSFMTGenSvc::~AtDSFMTGenSvc(){engineIteri(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!!!StatusCodeAtDSFMTGenSvc::queryInterface(constInterfaceID&riid,void**ppvInterface){if(IAtRndmGenSvc::interfaceID().versionMatch(riid)){*ppvInterface=(IAtRndmGenSvc*)this;}else{// Interface is not directly available: try out a base classreturnAthService::queryInterface(riid,ppvInterface);}addRef();returnStatusCode::SUCCESS;}StatusCodeAtDSFMTGenSvc::initialize(){ATH_MSG_INFO("Initializing "<<name()<<" - package version "<<PACKAGE_VERSION<<"\n INITIALISING RANDOM NUMBER STREAMS. ");/// Incident ServiceServiceHandle<IIncidentSvc>pIncSvc("IncidentSvc",name());// set up the incident service:if(!(pIncSvc.retrieve()).isSuccess()){ATH_MSG_ERROR("Could not locate IncidentSvc ");returnStatusCode::FAILURE;}//start listening to "EndEvent"staticconstintPRIORITY=100;pIncSvc->addListener(this,"EndEvent",PRIORITY);//and to BeginEvent if we are reseedingif(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 fileifstreaminfile(m_file_to_read.c_str());if(!infile){ATH_MSG_ERROR(" Unable to open: "<<m_file_to_read);returnStatusCode::FAILURE;}else{std::stringbuffer;while(std::getline(infile,buffer)){stringstream;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_iteratori=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;returnStatusCode::FAILURE;}}else{msg(MSG::ERROR)<<"bad line\n"<<buffer<<"\n in input file "<<m_file_to_read<<endreq;returnStatusCode::FAILURE;}}}}// Create the various streams according to user's requestfor(VStrings::const_iteratori=m_streams_seeds.begin();i!=m_streams_seeds.end();++i){stringstream;uint32_tseed1,seed2,offset(0);//parse the stream property stringshortll(0);// temp copy so we don't overwrite defaultif(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));returnStatusCode::FAILURE;}// Check if stream already generated (e.g. from reading a file)boolnot_found(true);if(number_of_streams()!=0){engineConstItersf(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 insertedATH_MSG_DEBUG("Applying reseeding offset "<<offset<<" to stream "<<stream);this->setOnDefinedSeeds(seed1,seed2,stream);}}}returnStatusCode::SUCCESS;}voidAtDSFMTGenSvc::handle(constIncident&inc){ATH_MSG_DEBUG(" Handle EndEvent ");if(inc.type()=="EndEvent"){m_engines_copy.clear();engineConstIteriE(begin()),eE(end());while(iE!=eE){HepRandomEngine*engine=GetEngine(iE->first);std::vector<unsignedlong>v=engine->put();std::vector<uint32_t>tseeds(v.size());for(size_ti=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();}elseif(inc.type()=="BeginEvent"){ATH_MSG_DEBUG(" Handle BeginEvent ");EventID*pei((dynamic_cast<constEventIncident&>(inc)).eventInfo().event_ID());//loop over generator streams, combining the stream name to the hashvector<string>::const_iteratori(m_reseedStreamNames.begin());vector<string>::const_iteratore(m_reseedStreamNames.end());//by default (when no streams are specified in streamNames, seed all//streamsif(i==e){if(!(this->setAllOnDefinedSeeds(pei->event_number(),pei->run_number())))throwGaudiException("can not reseed all streams ",name(),StatusCode::FAILURE);}else{while(i!=e){conststring&strName(*i++);if(0==this->setOnDefinedSeeds(pei->event_number(),pei->run_number(),strName)){throwGaudiException(string("can not reseed stream ")+strName,name(),StatusCode::FAILURE);}else{msg()<<MSG::VERBOSE<<"Reseeded stream "<<strName<<" for random service "<<endmsg;}}}}}StatusCodeAtDSFMTGenSvc::finalize(){ATH_MSG_INFO(" FINALISING ");if(m_save_to_file){// Write the status of the Service to filestd::ofstreamoutfile(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_iteratori=m_engines_copy.begin();i!=m_engines_copy.end();++i){outfile<<(*i).first<<" ";for(std::vector<uint32_t>::const_iteratorj=(*i).second.begin();j!=(*i).second.end();++j){outfile<<(*j)<<" ";}outfile<<endl;}ATH_MSG_DEBUG(" wrote seeds to "<<m_file_to_write);}}returnAthService::finalize();}HepRandomEngine*AtDSFMTGenSvc::GetEngine(conststd::string&streamName){engineConstIterciter=m_engines.find(streamName);if(citer==m_engines.end()){m_engines.insert(engineValType(streamName,newCLHEP::dSFMTEngine()));SetStreamSeeds(streamName);}engineIteriter=m_engines.find(streamName);return(HepRandomEngine*)(*iter).second;}voidAtDSFMTGenSvc::CreateStream(uint32_tseed1,uint32_tseed2,conststd::string&streamName){longseeds[3]={seed1,seed2,0};constlong*s=seeds;engineConstIterciter=m_engines.find(streamName);if(citer==m_engines.end())m_engines.insert(engineValType(streamName,newCLHEP::dSFMTEngine(s)));engineIteriter=m_engines.find(streamName);((*iter).second)->setSeeds(s,0);}boolAtDSFMTGenSvc::CreateStream(conststd::vector<uint32_t>&seeds,conststd::string&streamName){engineConstIterciter=m_engines.find(streamName);if(citer==m_engines.end())m_engines.insert(engineValType(streamName,newCLHEP::dSFMTEngine()));engineIteriter=m_engines.find(streamName);std::vector<unsignedlong>longSeeds(seeds.size());for(size_ti=0;i<seeds.size();++i)longSeeds[i]=seeds[i];return(((*iter).second)->getState(longSeeds));}voidAtDSFMTGenSvc::SetStreamSeeds(conststd::string&StreamName){uint32_tseed1;uint32_tseed2;if(StreamName=="PYTHIA"){seed1=m_PYTHIA_default_seed1;seed2=m_PYTHIA_default_seed2;}elseif(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);longseeds[3]={seed1,seed2,0};constlong*s=seeds;engineIteriter=m_engines.find(StreamName);((*iter).second)->setSeeds(s,0);}voidAtDSFMTGenSvc::print(conststd::string&StreamName){engineConstIterciter=m_engines.find(StreamName);if(citer==m_engines.end()){ATH_MSG_WARNING(" Stream = "<<StreamName<<" NOT FOUND");}else{std::vector<unsignedlong>v=((*citer).second)->put();msg(MSG::INFO)<<" Stream = "<<StreamName<<" ";for(std::vector<unsignedlong>::const_iteratori=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;}}voidAtDSFMTGenSvc::print(void){for(engineConstIteri=m_engines.begin();i!=m_engines.end();++i)print((*i).first);}HepRandomEngine*AtDSFMTGenSvc::setOnDefinedSeeds(uint32_teventNumber,uint32_trunNumber,conststd::string&streamName){uint32_ttheHash(eventNumber);map<string,uint32_t>::const_iteratorciter(m_reseedingOffsets.find(streamName));boolhasOffset(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);returnthis->setOnDefinedSeeds(theHash,streamName);}HepRandomEngine*AtDSFMTGenSvc::setOnDefinedSeeds(uint32_ttheSeed,conststd::string&streamName){engineConstIterciter=m_engines.find(streamName);///create engine if not found. FIXME this may not be a good ideaif(citer==m_engines.end())m_engines.insert(engineValType(streamName,newCLHEP::dSFMTEngine()));engineIteriter=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;}boolAtDSFMTGenSvc::setAllOnDefinedSeeds(uint32_teventNumber,uint32_trunNumber){boolallOK(true);engineIteri(m_engines.begin()),e(m_engines.end());while(i!=e&&(allOK=(0!=this->setOnDefinedSeeds(eventNumber,runNumber,(*i++).first)))){/*empty*/}returnallOK;}boolAtDSFMTGenSvc::setAllOnDefinedSeeds(uint32_ttheSeed){boolallOK(true);engineIteri(m_engines.begin()),e(m_engines.end());while(i!=e&&(allOK=(0!=this->setOnDefinedSeeds(theSeed,(*i++).first)))){/*empty*/}returnallOK;}