Various updates and fixes

* store TGeoManager::GetPath() in MCStepLogger if track enters to
  restore geometry in replay. This is necessary in particular in case of
  steps on or in the vincinity of boundaries. Floating point precision
  might not be enough to find a node based on spatial coordinates. In
  addition, restoring from the path makes navigation and node finding
  obsolete during replay

* extract TGEoVolume ID to VMC volume ID mapping. These are not the same
  for Geant4_VMC hence it is necessary. Mapping is required in replay,
  otheriwse, hit processing can fail.
  Mapping during replay is skipped if it is the identity.

* make replay work for Geant4_VMC. This works thanks to above additions.

* add possibility to load a custom user keepStep function from a ROOT
  macro

* update MCReplay/README.md

* additional small fixes

Author: Benedikt Volkel

MCReplay/README.md

@@ -4,32 +4,36 @@ This is a pseudo-detector-simulation engine based on the [Virtual Monte Carlo (V
 
 The main objective is to be able to provide an engine to study the impact of parameter variations with. Eventually, this can be used to optimise (full-)simulation parameters in view of enhancing their speed and efficiency.
 
-The functionality is compiled into a separate library `libTMCReplay` which also allows to use the `libMCStepLoggerIntercept` which also allows to run the step logging against the `TMCReplay` engine.
+The functionality is compiled into a separate library `libMCReplayCore`. Hence, `libMCStepLoggerIntercept` is completely independent which also allows to run the step logging against the `MCReplayEngine` engine.
+
+The replay has been tested and is verified against [GEANT3_VMC](https://github.com/vmc-project/geant3) and [GEANT4_VMC](https://github.com/vmc-project/geant4_vmc) VMC interfaces.
 
 ## Setting parameters
 
-The list of tunable parameters is foreseen to correspond to the possible settings implemented for [GEANT3_VMC]() and [GEANT4_VMC](https://github.com/vmc-project/geant4_vmc) and can be found [here](include/TMCReplay/Physics.h). The parameters can be set via the usual interfaces
-* `TVirtualMC::SetCut` and `TVirtualMC::SetProcess` to set global cut and process values
-* `TVirtualMC::Gstpar` to set parameters depending on a certain medium
+The list of tunable parameters is foreseen to correspond to the possible settings implemented for GEANT3_VMC and GEANT4_VMC and can be found [here](include/MCReplay/MCReplayPhysics.h). The parameters can be set via the usual interfaces
+* `TVirtualMC::SetCut` and `TVirtualMC::SetProcess` to set global cut and process values,
+* `TVirtualMC::Gstpar` to set parameters depending on a certain medium.
+The cut parameters act like production cuts.
 
 **Please note** that the parameters can be set but at the moment nothing is actually applied. Therefore, the exact same steps would be reproduced regardless of the parameters. Of course, this is the next major development as it is one of the core objectives of this engine.
+**On the other hand**, there is one special parameter called `CUTALLE` which is an energy production cut on **any** PDG which can also be set through the aforementioned interfaces. Note, that this parameter has no meaning to GEANT3_VMC or GEANT4_VMC.
 
 ## Testing different parameter settings
 
-**Note that this section becomes actually meaningful when the aforementioned functionality will have been implemented.** However, in the current state it can already be verified that the `TMCReplay` engine does what it is supposed to do.
+**Note that this section becomes actually meaningful when the aforementioned functionality will have been implemented.** However, in the current state it can already be verified that the `MCReplayEngine` engine does what it is supposed to do.
 
-### Running/testing on step level only (`TMCReplayDummyApplication` and `TMCReplayDummyStack`)
+### Running/testing on step level only (`MCReplayGenericApplication` and `MCReplayGenericStack`)
 
-If only the impact on the steps is desired, all required functionality is already included here. Simply have a look at [this executable](src/replay.cxx). It makes use of the `TMCReplayDummyApplication` and `TMCReplayDummyStack` just to replay the steps. Once implemented, the impact of the parameters on the stepping can be evaluated with this tool alone. To do so, just run the `MCStepLogger` together with the replay, for instance as follows
+If only the replay of steps is desired, all required functionality is already included here, completely independent of any environment the reference run was done. Simply have a look at [this code](src/replay.cxx) serving as the source for the executable `mcreplay`. It makes use of the `MCReplayGenericApplication` and `MCReplayGenericStack` just to replay the steps. Once implemented, the impact of the parameters on the stepping can be evaluated with this tool alone. To do so, just run the `MCStepLogger` followed by the replay, for instance as follows (assuming  everything was setup with `aliBuild`)
 
 ```bash
-MCSTEPLOG_TTREE=1 LD_PRELOAD=$MCSTEPLOGGER_ROOT/lib/libMCStepLoggerIntercept.so tmcreplay
+MCSTEPLOG_TTREE=1 LD_PRELOAD=$MCSTEPLOGGER_ROOT/lib/libMCStepLoggerIntercept.so mcreplay
 ```
 
 That will pick up the step file `MCStepLoggerOutput.root`. If you gave that a different name (or in case also your contained tree has another name), you will find the following help message useful:
 
 ```bash
-> tmcreplay --help
+> mcreplay --help
 Replaying a previously recorded particle transport step-by-step. Mainly meant for checking and performance measurements/optimisation of the transport:
   --help                                show this help message and exit
   --stepfilename arg (=MCStepLoggerOutput.root)
@@ -40,13 +44,18 @@ Replaying a previously recorded particle transport step-by-step. Mainly meant fo
                                         ROOT geometry filename
   --geokeyname arg (=FAIRGeom)          key name inside geo file where to find
                                         geometry tree
-
+  -n [ --nevents ] arg (=-1)            number of events to replay
+  -e [ --energycut ] arg (=-1)          energy cut to be applied [GeV]
 ```
 
-As you can see the geometry is required by this executable as well since it has otherwise no idea what geometry was used for the original simulation.
+As you can see the geometry is required by this executable as well since it has otherwise no idea what the geometry should look like from the `MCStepLoggerOutput.root` file alone.
 
 ### Running in another environment
 
 In that case the geometry construction is assumed to implemented by that environment's implementation of `TVirtualMCApplication`. Also, no stack is constructed automatically and it has to be constructed and passed by the user.
 
 Such a scenario is useful to evaluate the impact of parameter variation on hits or even digits.
+
+## Using the engine as a full simulation engine
+
+The `MCReplayEngine` has been tested with the [ALICE O2 framework](https://github.com/AliceO2Group/AliceO2) and also hits can be reproduced sufficiently which allows for realistic parameter tuning.

MCReplay/include/MCReplay/MCReplayEngine.h

@@ -13,7 +13,6 @@
 #define MC_REPLAY_ENGINE_H
 
 #include <vector>
-#include <unordered_map>
 #include <string>
 
 #include <TString.h>
@@ -36,6 +35,8 @@ class TGeoNode;
 namespace mcreplay
 {
 
+typedef bool (*user_keep_step_type)(const o2::StepInfo&, const o2::StepLookups*);
+
 class MCReplayEngine : public TVirtualMC
 {
 
@@ -1092,13 +1093,20 @@ class MCReplayEngine : public TVirtualMC
     mStepLoggerTreename = treename;
   }
 
+  // aacept a macro path from where to load a custom function whether or not to keep a step
+  void setUserKeepStepMacroPath(const std::string& path)
+  {
+    mUserKeepStepMacroPath = path;
+  }
+
  private:
   // init the run, used to guarantee that both ProcessRun and ProcessEvent
   // work just fine
   bool initRun();
 
-  // TODO Can most likely be removed
-  void adaptToTGeoName(const char* nameIn, char* nameOut) const;
+  // map volume IDs assigned by the TGeoManager (replay engine) to those assigned by the engine used in the reference run
+  Int_t correctGeoVolID(Int_t geoVolId) const;
+  Int_t correctVMCVolID(Int_t vmcVolId) const;
 
   // That is just a helper function to address some of the
   // interactions with the TGeoManager
@@ -1130,7 +1138,6 @@ class MCReplayEngine : public TVirtualMC
   {
     auto paramIndex = physics::paramToIndex(allParamsNames, paramName);
     if (paramIndex < 0) {
-      Warning("Could not set parameter %s, unknown and therefore skipped", paramName);
       return false;
     }
 
@@ -1151,7 +1158,6 @@ class MCReplayEngine : public TVirtualMC
   {
     auto paramIndex = physics::paramToIndex(allParamsNames, paramName);
     if (paramIndex < 0) {
-      Warning("Could not set parameter %s, unknown and therefore skipped", paramName);
       return false;
     }
     insertInto[paramIndex] = parval;
@@ -1174,6 +1180,11 @@ class MCReplayEngine : public TVirtualMC
   std::string mStepLoggerFilename;
   std::string mStepLoggerTreename;
 
+  // mapping volume IDs correctly; a reference run might assign other volume IDs than those which would be assigned only by the TGeoManager. Hence, we need to make sure that - whatever engine was used before - we mimic the exact same volume IDs. Otherwise, the logged IDs and those assigned by this engine based on the TGeoManager are not aligned
+  std::vector<int> mVolIDMap;
+  std::vector<int> mVolIDMapInverse;
+  bool mNeedVolIdMapping = false;
+
   // pointer to opened step file
   TFile* mStepFile = nullptr;
 
@@ -1190,6 +1201,11 @@ class MCReplayEngine : public TVirtualMC
   o2::StepLookups* mCurrentLookups = nullptr;
   o2::StepInfo* mCurrentStep = nullptr;
 
+  // keep track of tracks to be skipped
+  std::vector<bool> mSkipTrack;
+  // keep track of track ID assigned by user stack
+  std::vector<int> mUserTrackId;
+
   // current event ID
   int mCurrentEvent = 0;
 
@@ -1209,6 +1225,10 @@ class MCReplayEngine : public TVirtualMC
   // point to the current map of cuts
   std::vector<Double_t>* mCurrentCuts = nullptr;
 
+  // a user defined cut function, return true if step should be kept, false otherwise
+  std::string mUserKeepStepMacroPath;
+  user_keep_step_type* mUserKeepStep = nullptr;
+
   // local pointer to ROOT's geometry manager
   TGeoManager* mGeoManager = nullptr;
   // current TGeoVolume

MCReplay/include/MCReplay/MCReplayPhysics.h

@@ -23,8 +23,8 @@ namespace physics
 constexpr std::array<const char*, 15> namesProcesses = {"PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "HADR", "MUNU", "DCAY", "LOSS", "MULS", "CKOV", "RAYL", "LABS"};
 // standard GEANT cut parameter names (1-11)
 // additional Replay cut parameter names (12-n):
-// ALLE: energy cut applying to all tracks independent of PDG
-constexpr std::array<const char*, 12> namesCuts = {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "NCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX", "ALLE"};
+// CUTALLE: energy cut applying to all tracks independent of PDG
+constexpr std::array<const char*, 12> namesCuts = {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "NCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX", "CUTALLE"};
 constexpr const char* unknownParam = "UNKNOWN";
 
 template <typename T, std::size_t N>