Enable more cuts

* respect CUTGAM, CUTELE, CUTNEU, CUTHAD, CUTMUO

* enable loading cut parameters form JSON for more complex setting

* small fixes

Author: Benedikt Volkel

CMakeLists.txt

@@ -109,4 +109,4 @@ configure_file(
     "${PROJECT_BINARY_DIR}/MCStepLoggerConfig.cmake" @ONLY)
 install(FILES
     "${PROJECT_BINARY_DIR}/MCStepLoggerConfig.cmake"
-    DESTINATION ${INSTALL_CMAKE_DIR}) 
+    DESTINATION ${INSTALL_CMAKE_DIR})

MCReplay/CMakeLists.txt

@@ -33,6 +33,7 @@ set(SRCS
    ${CXX_SOURCE_DIR}/MCReplayGenericApplication.cxx
    ${CXX_SOURCE_DIR}/MCReplayGenericStack.cxx
    ${CXX_SOURCE_DIR}/MCReplayEvGen.cxx
+   ${CXX_SOURCE_DIR}/MCReplayPhysics.cxx
 )
 set(HEADERS
    ${CXX_INCLUDE_DIR}/${MODULE_NAME}/MCReplayEngine.h
@@ -49,9 +50,9 @@ set(EXECUTABLE_SRCS
 mcsl_add_library(
   TARGETNAME Core
   BASENAME ${MODULE_NAME}
-  DEPENDENCIES ROOT::Core ROOT::Hist ROOT::Graf ROOT::Gpad ROOT::Tree ROOT::VMC ROOT::RIO ROOT::Geom ROOT::EG MCStepLoggerCore MCStepLoggerAnalysis
+  DEPENDENCIES ROOT::Core ROOT::Hist ROOT::Graf ROOT::Gpad ROOT::Tree ROOT::VMC ROOT::RIO ROOT::Geom ROOT::EG MCStepLoggerCore
   SOURCES ${SRCS}
-  INCLUDE_DIRECTORIES ${CMAKE_CURRENT_LIST_DIR}/include
+  INCLUDE_DIRECTORIES ${CMAKE_CURRENT_LIST_DIR}/include ${Boost_INCLUDE_DIR}
   LINKDEFDIR ${CMAKE_CURRENT_LIST_DIR}/src
   ROOT_DICTIONARY_HEADERS ${HEADERS})
 

MCReplay/include/MCReplay/MCReplayEngine.h

@@ -1100,6 +1100,13 @@ class MCReplayEngine : public TVirtualMC
     mUserKeepStepMacroPath = path;
   }
 
+  void blockSetProcessesCuts(bool value = true)
+  {
+    mUpdateProcessesCutsBlocked = value;
+  }
+
+  void cutsFromConfig(std::string const& path);
+
  private:
   // init the run, used to guarantee that both ProcessRun and ProcessEvent
   // work just fine
@@ -1138,6 +1145,13 @@ class MCReplayEngine : public TVirtualMC
   bool insertProcessOrCut(std::vector<std::vector<P>*>& insertInto, const std::array<T, N>& allParamsNames, const std::vector<P>& defaultParams, Int_t mediumId, const char* paramName, P parval)
   {
     auto paramIndex = physics::paramToIndex(allParamsNames, paramName);
+    return insertProcessOrCut(insertInto, allParamsNames, defaultParams, mediumId, paramIndex, parval);
+  }
+
+  // add process or cut values based on name and value
+  template <typename P, typename T, std::size_t N>
+  bool insertProcessOrCut(std::vector<std::vector<P>*>& insertInto, const std::array<T, N>& allParamsNames, const std::vector<P>& defaultParams, Int_t mediumId, int paramIndex, P parval)
+  {
     if (paramIndex < 0) {
       return false;
     }
@@ -1158,6 +1172,12 @@ class MCReplayEngine : public TVirtualMC
   bool insertProcessOrCut(std::vector<P>& insertInto, const std::array<T, N>& allParamsNames, const char* paramName, P parval)
   {
     auto paramIndex = physics::paramToIndex(allParamsNames, paramName);
+    return insertProcessOrCut(insertInto, allParamsNames, paramIndex, parval);
+  }
+
+  template <typename P, typename T, std::size_t N>
+  bool insertProcessOrCut(std::vector<P>& insertInto, const std::array<T, N>& allParamsNames, int paramIndex, P parval)
+  {
     if (paramIndex < 0) {
       return false;
     }
@@ -1213,6 +1233,8 @@ class MCReplayEngine : public TVirtualMC
   // increment the current track length
   float mCurrentTrackLength = 0.;
 
+  // block in case framework should be prevented from setting cuts or processes
+  bool mUpdateProcessesCutsBlocked = false;
   // Preliminary process structure
   std::vector<std::vector<Int_t>*> mProcesses;
   // Preliminary cut structure

MCReplay/include/MCReplay/MCReplayPhysics.h

@@ -13,6 +13,8 @@
 #define MC_REPLAY_PHYSICS_H
 
 #include <array>
+#include <algorithm>
+#include <cstring>
 
 // Fix process and cut names
 namespace mcreplay
@@ -23,8 +25,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):
-// 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"};
+// CUTALL: energy cut applying to all tracks independent of PDG
+constexpr std::array<const char*, 12> namesCuts = {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX", "CUTALL"};
 constexpr const char* unknownParam = "UNKNOWN";
 
 template <typename T, std::size_t N>
@@ -46,6 +48,14 @@ const char* indexToParam(const std::array<T, N>& allParams, std::size_t index)
   return allParams[index];
 }
 
+bool isPhoton(int pdg);
+
+bool isElectronPositron(int pdg);
+
+bool isMuonAntiMuon(int pdg);
+
+bool isHadron(int pdg);
+
 } // namespace physics
 } // end namespace mcreplay
 

MCReplay/src/MCReplayEngine.cxx

@@ -12,6 +12,10 @@
 #include <algorithm>
 #include <iostream>
 
+#include <boost/program_options.hpp>
+#include <boost/property_tree/json_parser.hpp>
+#include <boost/property_tree/ptree_fwd.hpp>
+
 // For now include all TGeo headers here
 #include <TROOT.h>
 #include <TInterpreter.h>
@@ -35,8 +39,7 @@
 #include <TGeoEltu.h>
 #include <TGeoHype.h>
 #include <TGeoArb8.h>
-
-#include "TParticle.h"
+#include <TStopwatch.h>
 
 #include "MCReplay/MCReplayEngine.h"
 
@@ -739,11 +742,19 @@ int MCReplayEngine::getMediumId(int volId) const
 
 Bool_t MCReplayEngine::SetProcess(const char* flagName, Int_t flagValue)
 {
+  if (mUpdateProcessesCutsBlocked) {
+    ::Info("MCReplayEngine::Gstpar", "Parameter setting closed, nothing is changed.");
+    return false;
+  }
   return insertProcessOrCut(mProcessesGlobal, physics::namesCuts, flagName, flagValue);
 }
 
 Bool_t MCReplayEngine::SetCut(const char* cutName, Double_t cutValue)
 {
+  if (mUpdateProcessesCutsBlocked) {
+    ::Info("MCReplayEngine::Gstpar", "Parameter setting closed, nothing is changed.");
+    return false;
+  }
   return insertProcessOrCut(mCutsGlobal, physics::namesCuts, cutName, cutValue);
 }
 
@@ -842,8 +853,46 @@ void MCReplayEngine::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
   }
 }
 
+void MCReplayEngine::cutsFromConfig(std::string const& path)
+{
+  if (mUpdateProcessesCutsBlocked) {
+    ::Info("MCReplayEngine::paramsFromConfig", "Parameter setting closed, nothing is changed.");
+    return;
+  }
+
+  ::Info("MCReplayEngine::paramsFromConfig", "Read Params from %s", path.c_str());
+  if (path.empty()) {
+    return;
+  }
+  boost::property_tree::ptree pt;
+  boost::property_tree::json_parser::read_json(path, pt);
+
+  for (auto& it : pt) {
+    auto mediumID = std::stoi(it.first);
+    int i = 0;
+    for (auto& v : it.second) {
+      auto value = v.second.get_value<double>();
+      if (value < 0.) {
+        continue;
+      }
+      // TODO Going back and forth between parameter name and index. Implement so that we use the index directly
+      auto param = physics::indexToParam(physics::namesCuts, i);
+      if (mediumID < 0) {
+        insertProcessOrCut(mCutsGlobal, physics::namesCuts, i, value);
+      } else {
+        insertProcessOrCut(mCuts, physics::namesCuts, mCutsGlobal, mediumID, i, value);
+      }
+      i++;
+    }
+  }
+}
+
 void MCReplayEngine::Gstpar(Int_t itmed, const char* param, Double_t parval)
 {
+  if (mUpdateProcessesCutsBlocked) {
+    ::Info("MCReplayEngine::Gstpar", "Parameter setting closed, nothing is changed.");
+    return;
+  }
   if (insertProcessOrCut(mProcesses, physics::namesProcesses, mProcessesGlobal, itmed, param, (int)parval)) {
     return;
   }
@@ -878,6 +927,24 @@ bool MCReplayEngine::keepDueToCuts(const o2::StepInfo& step) const
     // check global energy cut
     return false;
   }
+  auto pdg = mCurrentLookups->tracktopdg[step.trackID];
+  if ((*mCurrentCuts)[0] > 0. && physics::isPhoton(pdg) && step.E < (*mCurrentCuts)[0]) {
+    return false;
+  }
+  if ((*mCurrentCuts)[1] > 0. && physics::isElectronPositron(pdg) && step.E < (*mCurrentCuts)[1]) {
+    return false;
+  }
+  if (physics::isHadron(pdg)) {
+    if (mCurrentLookups->tracktocharge[step.trackID] && (*mCurrentCuts)[3] > 0. && step.E < (*mCurrentCuts)[3]) {
+      return false;
+    }
+    if ((*mCurrentCuts)[2] > 0. && step.E < (*mCurrentCuts)[2]) {
+      return false;
+    }
+  }
+  if ((*mCurrentCuts)[4] > 0. && physics::isMuonAntiMuon(pdg) && step.E < (*mCurrentCuts)[4]) {
+    return false;
+  }
   return true;
 }
 
@@ -959,6 +1026,8 @@ void MCReplayEngine::ProcessEvent(Int_t eventId)
   unsigned int nUserTracks{0};
   unsigned int nStopTrack{0};
 
+  TStopwatch stopwatch{};
+
   for (auto& step : *steps) {
 
     // TODO This should not happen. (see comment in header file for these flags)
@@ -1068,10 +1137,13 @@ void MCReplayEngine::ProcessEvent(Int_t eventId)
 
   fApplication->FinishEvent();
 
+  stopwatch.Stop();
+
   auto nStepsSkipped = steps->size() - nStepsKept;
   std::cout << "Original number, skipped, kept, skipped fraction and kept fraction of steps: " << steps->size() << " " << nStepsSkipped << " " << nStepsKept << " " << static_cast<float>(nStepsSkipped) / steps->size() << " " << static_cast<float>(nStepsKept) / steps->size() << "\n";
   std::cout << "In addition, " << nUserTracks << " tracks produced during hit creation were ignored\n";
   std::cout << "TVirtualMC::StopTrack was ignored " << nStopTrack << " times\n";
+  std::cout << "Real time: " << stopwatch.RealTime() << ", CPU time: " << stopwatch.CpuTime() << "\n";
 
   delete steps;
   delete mCurrentLookups;

MCReplay/src/MCReplayPhysics.cxx

@@ -0,0 +1,45 @@
+
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include <cmath>
+#include "MCReplay/MCReplayPhysics.h"
+
+bool mcreplay::physics::isPhoton(int pdg)
+{
+  return pdg == 22;
+}
+
+bool mcreplay::physics::isElectronPositron(int pdg)
+{
+  return std::abs(pdg) == 11;
+}
+
+bool mcreplay::physics::isMuonAntiMuon(int pdg)
+{
+  return std::abs(pdg) == 13;
+}
+
+bool mcreplay::physics::isHadron(int pdg)
+{
+  // Taken from Pythia8 ParticleDataEntry::isHadron
+
+  if (pdg <= 100 || (pdg >= 1000000 && pdg <= 9000000) || pdg >= 9900000) {
+    return false;
+  }
+  if (pdg == 130 || pdg == 310) {
+    return true;
+  }
+  if (pdg % 10 == 0 || (pdg / 10) % 10 == 0 || (pdg / 100) % 10 == 0) {
+    return false;
+  }
+  return true;
+}