Add interface for sim parameters to JSON (#8017)

benedikt-voelkel · web-flow · commit ee545bdaea27 · 2022-03-03T11:39:40.000+01:00
first step to give cut and process setting from the outside more
structure

introduce JSON parsing into MaterialManager for that matter

leave the previous interface for now where cuts and processes are set from
bare text files

introduce configurable params

MaterialManagerParam.paramFile to pass and overwrite parameters
from command line
MaterialManagerParam.outputFile to write parameters used in the sim to
JSON for further usage and inspection
introduced policy:
parameters set later will overwrite previous settings

usage to overwrite from cmd

```
o2-sim &lt;args&gt; --configKeyValues "MaterialManagerParam.inputFile=&lt;path/to/file.json&gt;"
to write to med_params.json

o2-sim &lt;args&gt; --configKeyValues "MaterialManagerParam.outputFile=med_params.json"
Of course, both configKeyValues can be passed simultaneously.

The output in med_params.json looks like

{
  "module_1": [
    {
      "local_id": 1,
      "global_id": 5,
      "medium_name": "med_1_1",
      "material_name": "mat_1_1",
      "cuts": {
        "CUTGAM": 0.1,
        "CUTELE": 0.1,
        "CUTNEU": 0.1,
        "CUTHAD": 0.1,
        "CUTMUO": 0.1,
        "BCUTE": 0.1,
        "BCUTM": 0.1,
        "...": 0.1},
      "processes": {
        "PAIR": 0,
        "COMP": 1,
        "PHOT": -1,
        "PFIS": 1,
        "DRAY": 0,
        "ANNI": 0,
        "BREM": 1,
        "...": -1}
    },
    {
      "local_id": 2,
      "global_id": 3,
      "medium_name": "med_1_2",
      "material_name": "mat_1_2",
      "cuts": {"...": "..."},
      "processes": {"...": "..."}
    }
  ],
  "module_N": [
    {
      "local_id": 10,
      "global_id": 50,
      "medium_name": "med_N_10",
      "material_name": "mat_N_10",
      "cuts": {"...": "..."},
      "processes": {"...": "..."}
    },
    {
      "local_id": 20,
      "global_id": 30,
      "medium_name": "med_N_20",
      "material_name": "mat_N_20",
      "cuts": {"...": "..."},
      "processes": {"...": "..."}
    }
  ],
  "default":
    {
      "cuts": {"...": "..."},
      "processes": {"...": "..."}
    }
    "enableSpecialCuts": true,
    "enableSpecialProcesses": false
}
```
This can be edited/extended and used again as input as well. The global_id is the ID assigned by ROOT's TGeo and it is indeed useful to have that available when internal information of module-local medium IDs is not available.

In addition, the latter is extremely useful to be used in parameter optimisation together with the MCReplay engine
diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt
@@ -13,6 +13,7 @@ o2_add_library(DetectorsBase
                SOURCES src/Detector.cxx
                        src/GeometryManager.cxx
                        src/MaterialManager.cxx
+                       src/MaterialManagerParam.cxx
                        src/Propagator.cxx
                        src/MatLayerCyl.cxx
                        src/MatLayerCylSet.cxx
@@ -41,6 +42,7 @@ o2_target_root_dictionary(DetectorsBase
                           HEADERS include/DetectorsBase/Detector.h
                                   include/DetectorsBase/GeometryManager.h
                                   include/DetectorsBase/MaterialManager.h
+                                  include/DetectorsBase/MaterialManagerParam.h
                                   include/DetectorsBase/Propagator.h
                                   include/DetectorsBase/Ray.h
                                   include/DetectorsBase/MatCell.h
diff --git a/Detectors/Base/include/DetectorsBase/MaterialManager.h b/Detectors/Base/include/DetectorsBase/MaterialManager.h
@@ -133,6 +133,8 @@ class MaterialManager
 
   /// load cuts and process flags from a data file (like AliRoot did)
   void loadCutsAndProcessesFromFile(const char* modname, const char* filename);
+  void loadCutsAndProcessesFromJSON(ESpecial special = ESpecial::kFALSE, std::string const& filename = "");
+  void writeCutsAndProcessesToJSON(std::string const& filename = "");
 
   /// Set flags whether to use special cuts and process settings
   void enableSpecialProcesses(bool val = true) { mApplySpecialProcesses = val; }
diff --git a/Detectors/Base/include/DetectorsBase/MaterialManagerParam.h b/Detectors/Base/include/DetectorsBase/MaterialManagerParam.h
@@ -0,0 +1,30 @@
+// 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.
+
+#ifndef DETECTORS_BASE_INCLUDE_MATERIALMANAGERPARAM_H_
+#define DETECTORS_BASE_INCLUDE_MATERIALMANAGERPARAM_H_
+
+#include "CommonUtils/ConfigurableParam.h"
+#include "CommonUtils/ConfigurableParamHelper.h"
+
+namespace o2
+{
+
+struct MaterialManagerParam : public o2::conf::ConfigurableParamHelper<MaterialManagerParam> {
+  std::string inputFile = "";
+  std::string outputFile = "";
+  // boilerplate stuff + make principal key "MaterialManagerParam"
+  O2ParamDef(MaterialManagerParam, "MaterialManagerParam");
+};
+
+} // namespace o2
+
+#endif /* DETECTORS_BASE_INCLUDE_MATERIALMANAGERPARAM_H_ */
diff --git a/Detectors/Base/src/DetectorsBaseLinkDef.h b/Detectors/Base/src/DetectorsBaseLinkDef.h
@@ -24,6 +24,7 @@
 #pragma link C++ class o2::base::GeometryManager + ;
 #pragma link C++ class o2::base::GeometryManager::MatBudgetExt + ;
 #pragma link C++ class o2::base::MaterialManager + ;
+#pragma link C++ class o2::MaterialManagerParam + ;
 
 #pragma link C++ class o2::base::Ray + ;
 #pragma link C++ class o2::base::MatCell + ;
diff --git a/Detectors/Base/src/MaterialManager.cxx b/Detectors/Base/src/MaterialManager.cxx
@@ -13,6 +13,7 @@
 /// \brief Implementation of the MaterialManager class
 
 #include "DetectorsBase/MaterialManager.h"
+#include "DetectorsBase/MaterialManagerParam.h"
 #include "TVirtualMC.h"
 #include "TString.h" // for TString
 #include <TGeoMedium.h>
@@ -26,8 +27,39 @@
 #endif
 #include <cassert>
 #include <set>
+#include "rapidjson/document.h"
+#include "rapidjson/istreamwrapper.h"
+#include "rapidjson/ostreamwrapper.h"
+#include "rapidjson/prettywriter.h"
 
 using namespace o2::base;
+namespace rj = rapidjson;
+
+namespace
+{
+/// helper to read/write cuts and processes from/to JSON
+template <typename K, typename V>
+void writeSingleJSONParamBatch(std::unordered_map<K, const char*> const& idToName, std::map<K, V> const& valMap, V defaultValue, rapidjson::Value& parent, rapidjson::Document::AllocatorType& a)
+{
+  for (auto& itName : idToName) {
+    auto itVal = valMap.find(itName.first);
+    if (itVal != valMap.end()) {
+      parent.AddMember(rj::Value(itName.second, std::strlen(itName.second), a), rj::Value(itVal->second), a);
+      continue;
+    }
+    parent.AddMember(rj::Value(itName.second, std::strlen(itName.second), a), rj::Value(defaultValue), a);
+  }
+}
+
+/// specific names of keys wo expect and write in cut and process JSON files
+static constexpr const char* jsonKeyID = "local_id";
+static constexpr const char* jsonKeyIDGlobal = "global_id";
+static constexpr const char* jsonKeyDefault = "default";
+static constexpr const char* jsonKeyCuts = "cuts";
+static constexpr const char* jsonKeyProcesses = "processes";
+static constexpr const char* jsonKeyEnableSpecialCuts = "enableSpecialCuts";
+static constexpr const char* jsonKeyEnableSpecialProcesses = "enableSpecialProcesses";
+} // namespace
 
 const std::unordered_map<EProc, const char*> MaterialManager::mProcessIDToName = {
   {EProc::kPAIR, "PAIR"},
@@ -178,12 +210,16 @@ void MaterialManager::Cut(ESpecial special, int globalindex, ECut cut, Float_t v
     return;
   }
   if (special == ESpecial::kFALSE) {
-    mDefaultCutMap[cut] = val;
+    auto ins = mDefaultCutMap.insert({cut, val});
+    if (ins.second) {
+      TVirtualMC::GetMC()->SetCut(it->second, val);
+    }
     /// Explicit template definition to cover this which differs from global cut setting
-    TVirtualMC::GetMC()->SetCut(it->second, val);
   } else if (mApplySpecialCuts) {
-    mMediumCutMap[globalindex][cut] = val;
-    TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
+    auto ins = mMediumCutMap[globalindex].insert({cut, val});
+    if (ins.second) {
+      TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
+    }
   }
 }
 
@@ -198,12 +234,16 @@ void MaterialManager::Process(ESpecial special, int globalindex, EProc process,
     return;
   }
   if (special == ESpecial::kFALSE) {
-    mDefaultProcessMap[process] = val;
+    auto ins = mDefaultProcessMap.insert({process, val});
+    if (ins.second) {
+      TVirtualMC::GetMC()->SetProcess(it->second, val);
+    }
     /// Explicit template definition to cover this which differs from global process setting
-    TVirtualMC::GetMC()->SetProcess(it->second, val);
   } else if (mApplySpecialProcesses) {
-    mMediumProcessMap[globalindex][process] = val;
-    TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
+    auto ins = mMediumProcessMap[globalindex].insert({process, val});
+    if (ins.second) {
+      TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
+    }
   }
 }
 
@@ -326,6 +366,169 @@ TGeoMedium* MaterialManager::getTGeoMedium(const char* mediumname)
   return med;
 }
 
+void MaterialManager::loadCutsAndProcessesFromJSON(ESpecial special, std::string const& filename)
+{
+  const std::string filenameIn = filename.empty() ? o2::MaterialManagerParam::Instance().inputFile : filename;
+  if (filenameIn.empty()) {
+    return;
+  }
+  std::ifstream is(filenameIn);
+  if (!is.is_open()) {
+    LOG(error) << "Cannot open file " << filenameIn;
+    return;
+  }
+  auto digestCutsFromJSON = [this](int globalindex, rj::Value& cuts) {
+    auto special = globalindex < 0 ? ESpecial::kFALSE : ESpecial::kTRUE;
+    for (auto& cut : cuts.GetObject()) {
+      auto name = cut.name.GetString();
+      bool found = false;
+      for (auto& cn : mCutIDToName) {
+        if (std::strcmp(name, cn.second) == 0) {
+          Cut(special, globalindex, cn.first, cut.value.GetFloat());
+          found = true;
+        }
+      }
+      if (!found) {
+        LOG(warn) << "Unknown cut parameter " << name;
+      }
+    }
+  };
+  auto digestProcessesFromJSON = [this](int globalindex, rj::Value& processes) {
+    auto special = globalindex < 0 ? ESpecial::kFALSE : ESpecial::kTRUE;
+    for (auto& proc : processes.GetObject()) {
+      auto name = proc.name.GetString();
+      for (auto& pn : mProcessIDToName) {
+        if (std::strcmp(name, pn.second) == 0) {
+          Process(special, globalindex, pn.first, proc.value.GetInt());
+        }
+      }
+    }
+  };
+
+  rj::IStreamWrapper isw(is);
+  rj::Document d;
+  d.ParseStream(isw);
+
+  if (special == ESpecial::kFALSE && d.HasMember(jsonKeyDefault)) {
+    // defaults
+    auto& defaultParams = d[jsonKeyDefault];
+    if (defaultParams.HasMember(jsonKeyCuts)) {
+      digestCutsFromJSON(-1, defaultParams[jsonKeyCuts]);
+    }
+    if (defaultParams.HasMember(jsonKeyProcesses)) {
+      digestProcessesFromJSON(-1, defaultParams[jsonKeyProcesses]);
+    }
+  } else if (special == ESpecial::kTRUE) {
+    // read whether to apply special cuts and processes at all
+    if (d.HasMember(jsonKeyEnableSpecialCuts)) {
+      enableSpecialCuts(d[jsonKeyEnableSpecialCuts].GetBool());
+    }
+    if (d.HasMember(jsonKeyEnableSpecialProcesses)) {
+      enableSpecialProcesses(d[jsonKeyEnableSpecialProcesses].GetBool());
+    }
+    // special
+    for (auto& m : d.GetObject()) {
+      if (m.name.GetString()[0] == '\0' || !m.value.IsArray()) {
+        // do not parse anything with empty key, these at the most meant to be comments
+        continue;
+      }
+      for (auto& batch : m.value.GetArray()) {
+        if (std::strcmp(m.name.GetString(), jsonKeyDefault) == 0) {
+          // don't do defaults here
+          continue;
+        }
+        // set via their global indices
+        auto index = getMediumID(m.name.GetString(), batch[jsonKeyID].GetInt());
+        if (index < 0) {
+          continue;
+        }
+        if (batch.HasMember(jsonKeyCuts)) {
+          digestCutsFromJSON(index, batch[jsonKeyCuts]);
+        }
+        if (batch.HasMember(jsonKeyProcesses)) {
+          digestProcessesFromJSON(index, batch[jsonKeyProcesses]);
+        }
+      }
+    }
+  }
+}
+
+void MaterialManager::writeCutsAndProcessesToJSON(std::string const& filename)
+{
+  const std::string filenameOut = filename.empty() ? o2::MaterialManagerParam::Instance().outputFile : filename;
+  if (filenameOut.empty()) {
+    return;
+  }
+
+  // write parameters as global AND module specific
+  std::ofstream os(filenameOut);
+  if (!os.is_open()) {
+    LOG(error) << "Cannot create file " << filenameOut;
+    return;
+  }
+
+  rj::Document d;
+  rj::Document::AllocatorType& a = d.GetAllocator();
+  d.SetObject();
+
+  // add each local medium with params per module
+  for (auto& itMed : mMediumMap) {
+    // prepare array for module
+    rj::Value toAdd(rj::kArrayType);
+    // extract each medium's local and global index
+    for (auto& locToGlob : itMed.second) {
+      auto globalindex = locToGlob.second;
+      auto itCut = mMediumCutMap.find(globalindex);
+      auto itProc = mMediumProcessMap.find(globalindex);
+      // prepare a batch summarising localID, globaldID, cuts and processes
+      rj::Value oLoc(rj::kObjectType);
+      // IDs
+      oLoc.AddMember(rj::Value(jsonKeyID, std::strlen(jsonKeyID), a), rj::Value(locToGlob.first), a);
+      oLoc.AddMember(rj::Value(jsonKeyIDGlobal, std::strlen(jsonKeyIDGlobal)), rj::Value(locToGlob.second), a);
+      // add medium and material name
+      auto mediumIt = mTGeoMediumMap.find({itMed.first, locToGlob.first});
+      const char* medName = mediumIt->second->GetName();
+      const char* matName = mediumIt->second->GetMaterial()->GetName();
+      // not using variables for key names cause they are only written for info but not read
+      oLoc.AddMember(rj::Value("medium_name", 11, a), rj::Value(medName, std::strlen(medName), a), a);
+      oLoc.AddMember(rj::Value("material_name", 13, a), rj::Value(matName, std::strlen(matName), a), a);
+      // prepare for cuts
+      if (itCut != mMediumCutMap.end()) {
+        rj::Value cutMap(rj::kObjectType);
+        writeSingleJSONParamBatch(mCutIDToName, itCut->second, -1.f, cutMap, a);
+        oLoc.AddMember(rj::Value(jsonKeyCuts, std::strlen(jsonKeyCuts), a), cutMap, a);
+      }
+      // prepare for processes
+      if (itProc != mMediumProcessMap.end()) {
+        rj::Value procMap(rj::kObjectType);
+        writeSingleJSONParamBatch(mProcessIDToName, itProc->second, -1, procMap, a);
+        oLoc.AddMember(rj::Value(jsonKeyProcesses, std::strlen(jsonKeyProcesses), a), procMap, a);
+      }
+      // append this medium to module array
+      toAdd.PushBack(oLoc, a);
+    }
+    // append the entire module array
+    d.AddMember(rj::Value(itMed.first.c_str(), itMed.first.size(), a), toAdd, a);
+  }
+  // also add default parameters
+  rj::Value cutMapDef(rj::kObjectType);
+  rj::Value procMapDef(rj::kObjectType);
+  writeSingleJSONParamBatch(mCutIDToName, mDefaultCutMap, -1.f, cutMapDef, a);
+  writeSingleJSONParamBatch(mProcessIDToName, mDefaultProcessMap, -1, procMapDef, a);
+  rj::Value defaultParams(rj::kObjectType);
+  defaultParams.AddMember(rj::Value(jsonKeyCuts, std::strlen(jsonKeyCuts), a), cutMapDef, a);
+  defaultParams.AddMember(rj::Value(jsonKeyProcesses, std::strlen(jsonKeyProcesses), a), procMapDef, a);
+  d.AddMember(rj::Value(jsonKeyDefault, std::strlen(jsonKeyDefault), a), defaultParams, a);
+
+  d.AddMember(rj::Value(jsonKeyEnableSpecialCuts, std::strlen(jsonKeyEnableSpecialCuts), a), rj::Value(mApplySpecialCuts), a);
+  d.AddMember(rj::Value(jsonKeyEnableSpecialProcesses, std::strlen(jsonKeyEnableSpecialProcesses), a), rj::Value(mApplySpecialProcesses), a);
+  // now write to file
+  rj::OStreamWrapper osw(os);
+  rj::PrettyWriter<rj::OStreamWrapper> writer(osw);
+  writer.SetIndent(' ', 2);
+  d.Accept(writer);
+}
+
 void MaterialManager::loadCutsAndProcessesFromFile(const char* modname, const char* filename)
 {
   // Implementation of a method to set cuts and processes as done in AliRoot.
@@ -446,6 +649,8 @@ void MaterialManager::SpecialProcess(const char* modname, int localindex, EProc
   int globalindex = getMediumID(modname, localindex);
   if (globalindex != -1) {
     Process(ESpecial::kTRUE, globalindex, parID, val);
+  } else {
+    LOG(warn) << "SpecialProcess: NO GLOBALINDEX FOUND FOR " << modname << " " << localindex;
   }
 }
 
diff --git a/Detectors/Base/src/MaterialManagerParam.cxx b/Detectors/Base/src/MaterialManagerParam.cxx
@@ -0,0 +1,13 @@
+// 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 "DetectorsBase/MaterialManagerParam.h"
+O2ParamImpl(o2::MaterialManagerParam);
diff --git a/Detectors/gconfig/src/SetCuts.cxx b/Detectors/gconfig/src/SetCuts.cxx
@@ -37,6 +37,9 @@ void SetCuts()
   // \note All following settings could also be set in Cave since it is always loaded.
   // Use MaterialManager to set processes and cuts
   auto& mgr = o2::base::MaterialManager::Instance();
+  // This loads default cuts and processes if they are defined in the MaterialManagerParam.inputFile
+  // The cuts and processes below will only be set if they were not defined in the JSON
+  mgr.loadCutsAndProcessesFromJSON();
   auto& params = o2::GlobalProcessCutSimParam::Instance();
 
   LOG(info) << "Set default settings for processes and cuts.";
diff --git a/Steer/CMakeLists.txt b/Steer/CMakeLists.txt
@@ -14,6 +14,7 @@ o2_add_library(Steer
                        src/HitProcessingManager.cxx src/MCKinematicsReader.cxx
 		       PUBLIC_LINK_LIBRARIES O2::CommonDataFormat
 		                     O2::CommonConstants
+                             O2::DetectorsBase
                                      O2::SimulationDataFormat
                                      O2::DetectorsCommonDataFormats)
 
diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx
@@ -20,6 +20,7 @@
 #include <SimConfig/SimConfig.h>
 #include <DetectorsBase/Detector.h>
 #include "DetectorsBase/Aligner.h"
+#include "DetectorsBase/MaterialManager.h"
 #include <CommonUtils/ShmManager.h>
 #include <cassert>
 #include <SimulationDataFormat/MCEventHeader.h>
@@ -122,7 +123,12 @@ void O2MCApplicationBase::ConstructGeometry()
 
 void O2MCApplicationBase::InitGeometry()
 {
+  // load special cuts which might be given from the outside first.
+  auto& matMgr = o2::base::MaterialManager::Instance();
+  matMgr.loadCutsAndProcessesFromJSON(o2::base::MaterialManager::ESpecial::kTRUE);
+  // During the following, FairModule::SetSpecialPhysicsCuts will be called for each module
   FairMCApplication::InitGeometry();
+  matMgr.writeCutsAndProcessesToJSON();
   // now the sensitive volumes are set up in fVolMap and we can query them
   for (auto e : fVolMap) {
     // since fVolMap contains multiple entries (if multiple copies), this may
