Reproducible simulation on the track level

sawenzel · sawenzel · commit ce995a19485d · 2019-02-22T11:56:45.000+01:00
This commit introduces the possibly to have binary reproducible simulation
on the individual track level, in the sense that:

"If we change simulation of track A, it will never influence the simulation of track B
unless B is a child of A."

This is very useful for systematic optimization studies (for instance studying the effect
of geometry tracking limits on the physics). Before this commit,
a change in tracking limit would have influenced the evolution of the overall simulation
leading to large statistical fluctuations in the hits.
Such fluctuations make it difficult to judge the usefulness of a parameter change or require
a large production study.

We use here the idea of seeding the RNG for each track -- based on a unique hash derived
from its properties and independent of the history of the track. (The hash function is needs review).

By default the option is disabled. It can be enabled via

`o2sim --configKeyValues = "SimCutParams.trackSeed=1`

For the moment, this only works with Geant3.

The commit also provides an updated macro to statistically analyse hits.
This macro can be used to detect if a simulation parameter change has an impact on physics.
diff --git a/Common/SimConfig/include/SimConfig/SimCutParams.h b/Common/SimConfig/include/SimConfig/SimCutParams.h
@@ -22,9 +22,12 @@ namespace conf
 // (mostly used in O2MCApplication stepping)
 struct SimCutParams : public o2::conf::ConfigurableParamHelper<SimCutParams> {
   bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication
+  bool trackSeed = false;    // per track seeding for track-reproducible mode
 
   double maxRTracking = 1E20;    // max R tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
   double maxAbsZTracking = 1E20; // max |Z| tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
+  double ZmaxA = 1E20;           // max Z tracking cut on A side in cm -- applied in the stepping function
+  double ZmaxC = 1E20;           // max Z tracking cut on C side in cm -- applied in the stepping function
 
   O2ParamDef(SimCutParams, "SimCutParams");
 };
diff --git a/DataFormats/simulation/src/Stack.cxx b/DataFormats/simulation/src/Stack.cxx
@@ -16,6 +16,7 @@
 #include "DetectorsBase/Detector.h"
 #include "DetectorsCommonDataFormats/DetID.h"
 #include "SimulationDataFormat/MCTrack.h"
+#include "SimConfig/SimCutParams.h"
 
 #include "FairDetector.h" // for FairDetector
 #include "FairLogger.h"   // for FairLogger
@@ -245,6 +246,28 @@ void Stack::notifyFinishPrimary()
   }
 }
 
+// calculates a hash based on particle properties
+// hash may serve as seed for this track
+ULong_t getHash(TParticle const& p)
+{
+  auto asLong = [](double x) {
+    return (ULong_t) * (reinterpret_cast<ULong_t*>(&x));
+  };
+
+  ULong_t hash;
+  o2::MCTrackT<double> track(p);
+
+  hash = asLong(track.GetStartVertexCoordinatesX());
+  hash ^= asLong(track.GetStartVertexCoordinatesY());
+  hash ^= asLong(track.GetStartVertexCoordinatesZ());
+  hash ^= asLong(track.GetStartVertexCoordinatesT());
+  hash ^= asLong(track.GetStartVertexMomentumX());
+  hash ^= asLong(track.GetStartVertexMomentumY());
+  hash ^= asLong(track.GetStartVertexMomentumZ());
+  hash += (ULong_t)track.GetPdgCode();
+  return hash;
+}
+
 TParticle* Stack::PopNextTrack(Int_t& iTrack)
 {
   // This functions is mainly used by Geant3?
@@ -276,6 +299,16 @@ TParticle* Stack::PopNextTrack(Int_t& iTrack)
   mIndexOfCurrentTrack = mCurrentParticle.GetStatusCode();
   iTrack = mIndexOfCurrentTrack;
 
+  if (o2::conf::SimCutParams::Instance().trackSeed) {
+    auto hash = getHash(mCurrentParticle);
+    // LOG(INFO) << "SEEDING NEW TRACK USING HASH" << hash;
+    // init seed per track
+    gRandom->SetSeed(hash);
+
+    // NOTE: THE BETTER PLACE WOULD BE IN PRETRACK HOOK BUT THIS DOES NOT SEEM TO WORK
+    // WORKS ONLY WITH G3 SINCE G4 DOES NOT CALL THIS FUNCTION
+  }
+
   // LOG(INFO) << "transporting ID " << mIndexOfCurrentTrack << "\n";
   return &mCurrentParticle;
 }
diff --git a/Detectors/Passive/src/Cave.cxx b/Detectors/Passive/src/Cave.cxx
@@ -25,6 +25,7 @@
 #include "DetectorsBase/Detector.h"
 #include "DetectorsPassive/Cave.h"
 #include "SimConfig/SimConfig.h"
+#include "SimConfig/SimCutParams.h"
 #include <TRandom.h>
 #include "FairLogger.h"
 #include "TGeoManager.h"
@@ -105,14 +106,17 @@ void Cave::BeginPrimary()
 {
   static int primcounter = 0;
 
-  auto& conf = o2::conf::SimConfig::Instance();
-  auto chunks = conf.getInternalChunkSize();
-  if (chunks != -1) {
-    if (primcounter % chunks == 0) {
-      static int counter = 1;
-      auto seed = counter + 10;
-      gRandom->SetSeed(seed);
-      counter++;
+  // only do it, if not in pure trackSeeding mode
+  if (!o2::conf::SimCutParams::Instance().trackSeed) {
+    auto& conf = o2::conf::SimConfig::Instance();
+    auto chunks = conf.getInternalChunkSize();
+    if (chunks != -1) {
+      if (primcounter % chunks == 0) {
+        static int counter = 1;
+        auto seed = counter + 10;
+        gRandom->SetSeed(seed);
+        counter++;
+      }
     }
   }
   primcounter++;
diff --git a/Steer/include/Steer/O2MCApplication.h b/Steer/include/Steer/O2MCApplication.h
@@ -61,6 +61,7 @@ class O2MCApplication : public O2MCApplicationBase
       det->EndOfEvent();
     }
     fStack->Reset();
+    LOG(INFO) << "This event/chunk did " << mStepCounter << " steps";
   }
 
   /** Define actions at the end of run */
diff --git a/Steer/include/Steer/O2MCApplicationBase.h b/Steer/include/Steer/O2MCApplicationBase.h
@@ -42,13 +42,17 @@ class O2MCApplicationBase : public FairMCApplication
   ~O2MCApplicationBase() override = default;
 
   void Stepping() override;
+  void PreTrack() override;
+  void BeginEvent() override;
+  void FinishEvent() override;
 
   // specific implementation of our hard geometry limits
   double TrackingRmax() const override { return mCutParams.maxRTracking; }
   double TrackingZmax() const override { return mCutParams.maxAbsZTracking; }
 
  protected:
   o2::conf::SimCutParams const& mCutParams; // reference to parameter system
+  unsigned long long mStepCounter{ 0 };
 
   ClassDefOverride(O2MCApplicationBase, 1)
 };
diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx
@@ -45,10 +45,46 @@ void TypedVectorAttach(const char* name, FairMQChannel& channel, FairMQParts& pa
 
 void O2MCApplicationBase::Stepping()
 {
+  mStepCounter++;
+  if (mCutParams.stepFiltering) {
+    // we can kill tracks here based on our
+    // custom detector specificities
+
+    float x, y, z;
+    fMC->TrackPosition(x, y, z);
+
+    if (z > mCutParams.ZmaxA) {
+      fMC->StopTrack();
+      return;
+    }
+    if (-z > mCutParams.ZmaxC) {
+      fMC->StopTrack();
+      return;
+    }
+  }
+
   // dispatch first to stepping function in FairRoot
   FairMCApplication::Stepping();
+}
 
-  // now the right place to put our specific step/filtering criteria:
+void O2MCApplicationBase::PreTrack()
+{
+  // dispatch first to function in FairRoot
+  FairMCApplication::PreTrack();
+}
+
+void O2MCApplicationBase::FinishEvent()
+{
+  // dispatch first to function in FairRoot
+  FairMCApplication::FinishEvent();
+  LOG(INFO) << "This event/chunk did " << mStepCounter << " steps";
+}
+
+void O2MCApplicationBase::BeginEvent()
+{
+  // dispatch first to function in FairRoot
+  FairMCApplication::BeginEvent();
+  mStepCounter = 0;
 }
 
 void O2MCApplication::initLate()
diff --git a/macro/analyzeHits.C b/macro/analyzeHits.C
diff --git a/macro/compareHits.C b/macro/compareHits.C

Original file line number	Diff line number	Diff line change
`@@ -61,6 +61,7 @@ class O2MCApplication : public O2MCApplicationBase`
`61`	`61`	`det->EndOfEvent();`
`62`	`62`	`}`
`63`	`63`	`fStack->Reset();`
	`64`	`+ LOG(INFO) << "This event/chunk did " << mStepCounter << " steps";`
`64`	`65`	`}`
`65`	`66`
`66`	`67`	`/** Define actions at the end of run */`