From d92e3a5e55fa430c0b2ada9ae10baa48b936721d Mon Sep 17 00:00:00 2001
From: Andre Sailer <andre.philippe.sailer@cern.ch>
Date: Fri, 1 Jul 2016 14:52:48 +0000
Subject: [PATCH] Add plugin to skip events and still stay reproducible
Use same algorithm as used in Marlin EventSeeder
---
DDG4/include/DDG4/Geant4EventSeed.h | 205 ++++++++++++++++++++++++++++
DDG4/plugins/Geant4EventSeed.cpp | 89 ++++++++++++
doc/doxygen/DD4hepGroups.h | 3 +
3 files changed, 297 insertions(+)
create mode 100644 DDG4/include/DDG4/Geant4EventSeed.h
create mode 100644 DDG4/plugins/Geant4EventSeed.cpp
diff --git a/DDG4/include/DDG4/Geant4EventSeed.h b/DDG4/include/DDG4/Geant4EventSeed.h
new file mode 100644
index 000000000..4c4c70475
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4EventSeed.h
@@ -0,0 +1,205 @@
+// $Id: $
+//==========================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : A.Sailer
+//
+//==========================================================================
+#ifndef DD4HEP_DDG4_GEANT4EVENTSEED_H
+#define DD4HEP_DDG4_GEANT4EVENTSEED_H
+
+// Framework include files
+#include "DDG4/Geant4RunAction.h"
+
+/// Namespace for the AIDA detector description toolkit
+namespace DD4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace Simulation {
+
+ /**
+ * \addtogroup Geant4RunActions
+ * @{
+ * \package EventSeed
+ * \brief Set the event seed for each event
+ *
+ * This plugins allows one to skip events and still be reproducible.
+ * We are using jenkins_hash from an intial seed, runID and eventID
+ * see http://burtleburtle.net/bob/hash/evahash.html
+ *
+ * \author A.Sailer
+ * \version 1.0
+ * @}
+ */
+ class Geant4EventSeed: public Geant4RunAction {
+
+ protected:
+ unsigned int m_initialSeed;
+ unsigned int m_runID;
+ std::string m_type;
+ bool m_initialised;
+ public:
+ /// Standard constructor with initializing arguments
+ Geant4EventSeed(Geant4Context*, const std::string& );
+ /// Default destructor
+ virtual ~Geant4EventSeed();
+ /// begin-of-run callback
+ void begin(const G4Run*);
+ /// begin-of-event callback
+ void beginEvent(const G4Event*);
+ };
+
+ /*
+
+ Hashing for random seed as used in Marlin EventSeeder processor
+
+ Original source by Bob Jenkins
+
+ http://www.burtleburtle.net/bob/hash/doobs.html
+
+ Hash a variable-length key into a 32-bit value
+
+ */
+
+#define hashsize(n) ( 1U << (n) )
+#define hashmask(n) ( hashsize ( n ) - 1 )
+
+
+ /*
+ --------------------------------------------------------------------
+ mix -- mix 3 32-bit values reversibly.
+ For every delta with one or two bits set, and the deltas of all three
+ high bits or all three low bits, whether the original value of a,b,c
+ is almost all zero or is uniformly distributed,
+ * If mix() is run forward or backward, at least 32 bits in a,b,c
+ have at least 1/4 probability of changing.
+ * If mix() is run forward, every bit of c will change between 1/3 and
+ 2/3 of the time. (Well, 22/100 and 78/100 for some 2-bit deltas.)
+ mix() was built out of 36 single-cycle latency instructions in a
+ structure that could supported 2x parallelism, like so:
+ a -= b;
+ a -= c; x = (c>>13);
+ b -= c; a ^= x;
+ b -= a; x = (a<<8);
+ c -= a; b ^= x;
+ c -= b; x = (b>>13);
+ ...
+ Unfortunately, superscalar Pentiums and Sparcs can't take advantage
+ of that parallelism. They've also turned some of those single-cycle
+ latency instructions into multi-cycle latency instructions. Still,
+ this is the fastest good hash I could find. There were about 2^^68
+ to choose from. I only looked at a billion or so.
+ --------------------------------------------------------------------
+ */
+#define mix(a,b,c) \
+ { \
+ a -= b; a -= c; a ^= (c>>13); \
+ b -= c; b -= a; b ^= (a<<8); \
+ c -= a; c -= b; c ^= (b>>13); \
+ a -= b; a -= c; a ^= (c>>12); \
+ b -= c; b -= a; b ^= (a<<16); \
+ c -= a; c -= b; c ^= (b>>5); \
+ a -= b; a -= c; a ^= (c>>3); \
+ b -= c; b -= a; b ^= (a<<10); \
+ c -= a; c -= b; c ^= (b>>15); \
+ }
+
+ /*
+ --------------------------------------------------------------------
+ jenkins_hash() -- hash a variable-length key into a 32-bit value
+ k : the key (the unaligned variable-length array of bytes)
+ len : the length of the key, counting by bytes
+ initval : can be any 4-byte value
+ Returns a 32-bit value. Every bit of the key affects every bit of
+ the return value. Every 1-bit and 2-bit delta achieves avalanche.
+ About 6*len+35 instructions.
+
+ The best hash table sizes are powers of 2. There is no need to do
+ mod a prime (mod is sooo slow!). If you need less than 32 bits,
+ use a bitmask. For example, if you need only 10 bits, do
+ h = (h & hashmask(10));
+ In which case, the hash table should have hashsize(10) elements.
+
+ If you are hashing n strings (ub1 **)k, do it like this:
+ for (i=0, h=0; i<n; ++i) h = hash( k[i], len[i], h);
+
+ By Bob Jenkins, 1996. bob_jenkins@burtleburtle.net. You may use this
+ code any way you wish, private, educational, or commercial. It's free.
+
+ See http://burtleburtle.net/bob/hash/evahash.html
+ Use for hash table lookup, or anything where one collision in 2^^32 is
+ acceptable. Do NOT use for cryptographic purposes.
+ --------------------------------------------------------------------
+ */
+ unsigned jenkins_hash ( unsigned char *k, unsigned length, unsigned initval )
+ {
+ unsigned a, b;
+ unsigned c = initval;
+ unsigned len = length;
+
+ a = b = 0x9e3779b9;
+
+ while ( len >= 12 ) {
+ a += ( k[0] + ( (unsigned)k[1] << 8 )
+ + ( (unsigned)k[2] << 16 )
+ + ( (unsigned)k[3] << 24 ) );
+ b += ( k[4] + ( (unsigned)k[5] << 8 )
+ + ( (unsigned)k[6] << 16 )
+ + ( (unsigned)k[7] << 24 ) );
+ c += ( k[8] + ( (unsigned)k[9] << 8 )
+ + ( (unsigned)k[10] << 16 )
+ + ( (unsigned)k[11] << 24 ) );
+
+ mix ( a, b, c );
+
+ k += 12;
+ len -= 12;
+ }
+
+ c += length;
+
+ switch ( len ) {
+ case 11: c += ( (unsigned)k[10] << 24 );
+ case 10: c += ( (unsigned)k[9] << 16 );
+ case 9 : c += ( (unsigned)k[8] << 8 );
+ /* First byte of c reserved for length */
+ case 8 : b += ( (unsigned)k[7] << 24 );
+ case 7 : b += ( (unsigned)k[6] << 16 );
+ case 6 : b += ( (unsigned)k[5] << 8 );
+ case 5 : b += k[4];
+ case 4 : a += ( (unsigned)k[3] << 24 );
+ case 3 : a += ( (unsigned)k[2] << 16 );
+ case 2 : a += ( (unsigned)k[1] << 8 );
+ case 1 : a += k[0];
+ }
+
+ mix ( a, b, c );
+
+ return c;
+ }
+
+ /// calculate hash from initialSeed, eventID and runID
+ unsigned int hash( unsigned int initialSeed, unsigned int eventNumber, unsigned int runNumber ){
+ unsigned int seed = 0;
+ unsigned char * c = (unsigned char *) &eventNumber ;
+ seed = jenkins_hash( c, sizeof eventNumber, seed) ;
+
+ c = (unsigned char *) &runNumber ;
+ seed = jenkins_hash( c, sizeof runNumber, seed) ;
+
+ c = (unsigned char *) &initialSeed ;
+ seed = jenkins_hash( c, sizeof initialSeed, seed) ;
+
+ return seed;
+ }
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4EVENTSEED_H
diff --git a/DDG4/plugins/Geant4EventSeed.cpp b/DDG4/plugins/Geant4EventSeed.cpp
new file mode 100644
index 000000000..ccc9cd5a9
--- /dev/null
+++ b/DDG4/plugins/Geant4EventSeed.cpp
@@ -0,0 +1,89 @@
+// $Id: $
+//==========================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author : A.Sailer
+//
+//=========================================================================
+
+//Class include file
+#include "DDG4/Geant4EventSeed.h"
+
+// Framework include files
+#include "DD4hep/InstanceCount.h"
+#include "DD4hep/Printout.h"
+
+#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4Random.h"
+#include "DDG4/Factories.h"
+
+#include "CLHEP/Random/EngineFactory.h"
+
+//Geant includes
+#include <G4Run.hh>
+#include <G4Event.hh>
+
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4EventSeed::Geant4EventSeed(Geant4Context* c, const std::string& typ) : Geant4RunAction(c, typ),
+ m_initialSeed(0),
+ m_runID(0),
+ m_type(typ),
+ m_initialised(false)
+{
+ Geant4Action::runAction().callAtBegin(this,&Geant4EventSeed::begin);
+ Geant4Action::eventAction().callAtBegin(this,&Geant4EventSeed::beginEvent);
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4EventSeed::~Geant4EventSeed() {
+ InstanceCount::decrement(this);
+}
+
+/// begin-of-run callback
+void Geant4EventSeed::begin(const G4Run* run) {
+
+ if(not m_initialised){
+ m_initialised = true;
+ m_initialSeed = Geant4Random::instance()->engine()->getSeed();
+ }
+
+ m_runID = run->GetRunID();
+
+ DD4hep::printout( PrintLevel::INFO, m_type,
+ "Get RunID: runID=%u\n",
+ m_runID );
+
+}
+
+/// begin-of-event callback
+void Geant4EventSeed::beginEvent(const G4Event* evt) {
+
+ Geant4Random *rndm = Geant4Random::instance();
+
+ unsigned int eventID = evt->GetEventID();
+ unsigned int newSeed = hash( m_initialSeed, eventID, m_runID );
+
+ DD4hep::printout( PrintLevel::INFO, m_type,
+ "At beginEvent: eventID=%u, runID=%u" \
+ " initialSeed=%u, newSeed=%u" ,
+ evt->GetEventID(), m_runID,
+ m_initialSeed, newSeed );
+
+ rndm->setSeed( newSeed );
+
+ if ( DD4hep::printLevel() <= PrintLevel::DEBUG ) {
+ rndm->showStatus();
+ }
+
+}
+
+DECLARE_GEANT4ACTION(Geant4EventSeed)
diff --git a/doc/doxygen/DD4hepGroups.h b/doc/doxygen/DD4hepGroups.h
index 939a97110..a95bb5272 100644
--- a/doc/doxygen/DD4hepGroups.h
+++ b/doc/doxygen/DD4hepGroups.h
@@ -93,6 +93,9 @@ namespace UTIL {}
* \defgroup Geant4ActionPlugin
\brief Plugins that are a Geant4Actions
+ * \defgroup Geant4RunActions
+ \brief Plugins that are a Geant4RunActions
+
* \defgroup SurfacePlugin
\brief Plugins to manipulate surfaces automatically
--
GitLab