[ALICE3] Add A3 geometry provider task

ALICE3/Core/DelphesO2TrackSmearer.cxx

@@ -66,6 +66,7 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
     std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix
     const std::string outPath = "/tmp/LUTs/";
     filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str());
+    LOG(info) << " --- Local LUT filename will be: " << filename;
     std::ifstream checkFile(filename); // Check if file already exists
     if (!checkFile.is_open()) {        // File does not exist, retrieve from CCDB
       LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path;

ALICE3/Core/FastTracker.cxx

@@ -11,15 +11,18 @@
 
 #include "FastTracker.h"
 
-#include "ReconstructionDataFormats/TrackParametrization.h"
+#include "Common/Core/TableHelper.h"
+
+#include <ReconstructionDataFormats/TrackParametrization.h>
 
-#include "TMath.h"
-#include "TMatrixD.h"
-#include "TMatrixDSymEigen.h"
-#include "TRandom.h"
 #include <TEnv.h>
 #include <THashList.h>
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TMatrixDSymEigen.h>
 #include <TObject.h>
+#include <TRandom.h>
+#include <TSystem.h>
 
 #include <fstream>
 #include <map>
@@ -31,6 +34,81 @@
 namespace fastsim
 {
 
+std::map<std::string, std::map<std::string, std::string>> GeometryContainer::parseTEnvConfiguration(std::string filename, std::vector<std::string>& layers)
+{
+  std::map<std::string, std::map<std::string, std::string>> configMap;
+  filename = gSystem->ExpandPathName(filename.c_str());
+  TEnv env(filename.c_str());
+  THashList* table = env.GetTable();
+  layers.clear();
+  for (int i = 0; i < table->GetEntries(); ++i) {
+    const std::string key = table->At(i)->GetName();
+    // key should contain exactly one dot
+    if (key.find('.') == std::string::npos || key.find('.') != key.rfind('.')) {
+      LOG(fatal) << "Key " << key << " does not contain exactly one dot";
+      continue;
+    }
+    const std::string firstPart = key.substr(0, key.find('.'));
+    if (std::find(layers.begin(), layers.end(), firstPart) == layers.end()) {
+      layers.push_back(firstPart);
+    }
+  }
+  env.Print();
+  // Layers
+  for (const auto& layer : layers) {
+    LOG(info) << " Reading layer " << layer;
+    for (int i = 0; i < table->GetEntries(); ++i) {
+      const std::string key = table->At(i)->GetName();
+      if (key.find(layer + ".") == 0) {
+        const std::string paramName = key.substr(key.find('.') + 1);
+        const std::string value = env.GetValue(key.c_str(), "");
+        configMap[layer][paramName] = value;
+      }
+    }
+  }
+  return configMap;
+}
+
+void GeometryContainer::init(o2::framework::InitContext& initContext)
+{
+  std::vector<std::string> detectorConfiguration;
+  const bool found = common::core::getTaskOptionValue(initContext, "on-the-fly-detector-geometry-provider", "detectorConfiguration", detectorConfiguration, false);
+  if (!found) {
+    LOG(fatal) << "Could not retrieve detector configuration from OnTheFlyDetectorGeometryProvider task.";
+    return;
+  }
+  LOG(info) << "Size of detector configuration: " << detectorConfiguration.size();
+  for (const auto& configFile : detectorConfiguration) {
+    LOG(info) << "Detector geometry configuration file used: " << configFile;
+    addEntry(configFile);
+  }
+}
+
+std::map<std::string, std::string> GeometryContainer::GeometryEntry::getConfiguration(const std::string& layerName) const
+{
+  auto it = mConfigurations.find(layerName);
+  if (it != mConfigurations.end()) {
+    return it->second;
+  } else {
+    LOG(fatal) << "Layer " << layerName << " not found in geometry configurations.";
+    return {};
+  }
+}
+
+std::string GeometryContainer::GeometryEntry::getValue(const std::string& layerName, const std::string& key, bool require) const
+{
+  auto layer = getConfiguration(layerName);
+  auto entry = layer.find(key);
+  if (entry != layer.end()) {
+    return layer.at(key);
+  } else if (require) {
+    LOG(fatal) << "Key " << key << " not found in layer " << layerName << " configurations.";
+    return "";
+  } else {
+    return "";
+  }
+}
+
 // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
 
 DetLayer* FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
@@ -360,7 +438,7 @@
       index = 1;
     z0 = -4 * sigmaD + i * dz0;
     dist += index * (dz0 / 3.) * (1 / o2::math_utils::sqrt(o2::constants::math::TwoPI) / sigmaD) * std::exp(-z0 * z0 / 2. / sigmaD / sigmaD) * (1 / o2::math_utils::sqrt((z - z0) * (z - z0) + r * r));
     if (index != 4)
       index = 4;
     else
       index = 2;
@@ -490,7 +568,7 @@
     // check if layer is reached
     float targetX = 1e+3;
     inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField);
     if (targetX > 999.f) {
       LOGF(debug, "Failed to find intercept for layer %d at radius %.2f cm", il, layers[il].getRadius());
       break; // failed to find intercept
     }
@@ -572,7 +650,7 @@
 
     float targetX = 1e+3;
     inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField);
     if (targetX > 999)
       continue; // failed to find intercept
 
     if (!inputTrack.propagateTo(targetX, magneticField)) {
@@ -644,7 +722,7 @@
   // backpropagate to original radius
   float finalX = 1e+3;
   bool inPropStatus = inwardTrack.getXatLabR(initialRadius, finalX, magneticField);
   if (finalX > 999) {
     LOG(debug) << "Failed to find intercept for initial radius " << initialRadius << " cm, x = " << finalX << " and status " << inPropStatus << " and sn = " << inwardTrack.getSnp() << " r = " << inwardTrack.getY() * inwardTrack.getY();
     return -3; // failed to find intercept
   }
@@ -654,7 +732,7 @@
   }
 
   // only attempt to continue if intercepts are at least four
   if (nIntercepts < 4)
     return nIntercepts;
 
   // generate efficiency
@@ -681,7 +759,7 @@
   TMatrixDSym m(5);
   double fcovm[5][5]; // double precision is needed for regularisation
 
   for (int ii = 0, k = 0; ii < 5; ++ii) {
     for (int j = 0; j < ii + 1; ++j, ++k) {
       fcovm[ii][j] = covMat[k];
       fcovm[j][ii] = covMat[k];
@@ -689,7 +767,7 @@
   }
 
   // evaluate ruben's conditional, regularise
   const bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix
   bool rubenConditional = false;
   for (int ii = 0; ii < 5; ii++) {
     for (int jj = 0; jj < 5; jj++) {
@@ -698,7 +776,7 @@
       if (fcovm[ii][jj] * fcovm[ii][jj] > std::abs(fcovm[ii][ii] * fcovm[jj][jj])) {
         rubenConditional = true;
         if (makePositiveDefinite) {
           fcovm[ii][jj] = TMath::Sign(1, fcovm[ii][jj]) * covMatFactor * sqrt(std::abs(fcovm[ii][ii] * fcovm[jj][jj]));
         }
       }
     }
@@ -736,7 +814,7 @@
     for (int j = 0; j < 5; ++j)
       val += eigVec[j][ii] * outputTrack.getParam(j);
     // smear parameters according to eigenvalues
     params_[ii] = gRandom->Gaus(val, sqrt(eigVal[ii]));
   }
 
   // invert eigenvector matrix
@@ -749,7 +827,7 @@
     outputTrack.setParam(val, ii);
   }
   // should make a sanity check that par[2] sin(phi) is in [-1, 1]
   if (fabs(outputTrack.getParam(2)) > 1.) {
     LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << outputTrack.getParam(2);
     return -2;
   }

ALICE3/Core/FastTracker.h

@@ -15,10 +15,10 @@
 #include "DetLayer.h"
 
 #include <CCDB/BasicCCDBManager.h>
+#include <Framework/InitContext.h>
+#include <Framework/Logger.h>
 #include <ReconstructionDataFormats/Track.h>
 
-#include <fairlogger/Logger.h>
-
 #include <map>
 #include <string>
 #include <vector>
@@ -28,6 +28,63 @@ namespace o2
 namespace fastsim
 {
 
+class GeometryContainer
+{
+ public:
+  GeometryContainer() = default;
+  virtual ~GeometryContainer() = default;
+
+  void init(o2::framework::InitContext& initContext);
+
+  /**
+   * @brief Parses a TEnv configuration file and returns the key-value pairs split per entry
+   * @param filename Path to the TEnv configuration file
+   * @param layers Vector to store the order of the layers as they appear in the file
+   * @return A map where each key is a layer name and the value is another map of key-value pairs for that layer
+   */
+  static std::map<std::string, std::map<std::string, std::string>> parseTEnvConfiguration(std::string filename, std::vector<std::string>& layers);
+
+  // A container for the geometry info
+  struct GeometryEntry {
+    // Default constructor
+    GeometryEntry() = default;
+    explicit GeometryEntry(std::string filename) : name(filename)
+    {
+      mConfigurations = GeometryContainer::parseTEnvConfiguration(filename, layerNames);
+    }
+    std::map<std::string, std::map<std::string, std::string>> getConfigurations() const { return mConfigurations; }
+    std::map<std::string, std::string> getConfiguration(const std::string& layerName) const;
+    std::vector<std::string> getLayerNames() const { return layerNames; }
+    std::string getValue(const std::string& layerName, const std::string& key, bool require = true) const;
+    float getFloatValue(const std::string& layerName, const std::string& key) const { return std::stof(getValue(layerName, key)); }
+    int getIntValue(const std::string& layerName, const std::string& key) const { return std::stoi(getValue(layerName, key)); }
+
+   private:
+    std::string name; // Filename of the geometry
+    std::map<std::string, std::map<std::string, std::string>> mConfigurations;
+    std::vector<std::string> layerNames; // Ordered names of the layers
+  };
+
+  // Add a geometry entry from a configuration file
+  void addEntry(const std::string& filename) { entries.emplace_back(filename); }
+
+  // Getters
+  int getNumberOfConfigurations() const { return entries.size(); }
+  const std::vector<GeometryEntry>& getEntries() const { return entries; }
+  const GeometryEntry& getEntry(const int id) const { return entries.at(id); }
+  GeometryEntry getGeometryEntry(const int id) const { return entries.at(id); }
+
+  // Get configuration maps
+  std::map<std::string, std::map<std::string, std::string>> getConfigurations(const int id) const { return entries.at(id).getConfigurations(); }
+  std::map<std::string, std::string> getConfiguration(const int id, const std::string& layerName) const { return entries.at(id).getConfiguration(layerName); }
+
+  // Get specific values
+  float getFloatValue(const int id, const std::string& layerName, const std::string& key) const { return entries.at(id).getFloatValue(layerName, key); }
+
+ private:
+  std::vector<GeometryEntry> entries;
+};
+
 // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
 
 // this class implements a synthetic smearer that allows

ALICE3/Core/FastTrackerLinkDef.h

@@ -16,6 +16,7 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
+#pragma link C++ class o2::fastsim::GeometryContainer + ;
 #pragma link C++ class o2::fastsim::FastTracker + ;
 #pragma link C++ class o2::fastsim::DelphesO2LutWriter + ;
 

ALICE3/TableProducer/OTF/CMakeLists.txt

@@ -28,3 +28,8 @@ o2physics_add_dpl_workflow(on-the-fly-tracker-pid
                     SOURCES onTheFlyTrackerPid.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(on-the-fly-detector-geometry-provider
+                    SOURCES onTheFlyDetectorGeometryProvider.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::ALICE3Core O2Physics::FastTracker
+                    COMPONENT_NAME Analysis)

ALICE3/TableProducer/OTF/onTheFlyDetectorGeometryProvider.cxx

@@ -0,0 +1,100 @@
+// 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.
+
+/// \file onTheFlyDetectorGeometryProvider.cxx
+///
+/// \brief Reader and and provider of the detector geometry for on-the-fly simulation
+///
+/// \author Nicolò Jacazio <nicolo.jacazio@cern.ch>, Universita del Piemonte Orientale (IT)
+///
+
+#include "ALICE3/Core/FastTracker.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/runDataProcessing.h>
+
+#include <map>
+#include <string>
+#include <vector>
+
+struct OnTheFlyDetectorGeometryProvider {
+
+  o2::framework::Configurable<std::vector<std::string>> detectorConfiguration{"detectorConfiguration",
+                                                                              std::vector<std::string>{"$O2PHYSICS_ROOT/share/alice3/a3geometry_v3.ini"},
+                                                                              "Paths of the detector geometry configuration files"};
+  o2::framework::Service<o2::ccdb::BasicCCDBManager> ccdb;
+  void init(o2::framework::InitContext&)
+  {
+    ccdb->setURL("http://alice-ccdb.cern.ch");
+    ccdb->setTimestamp(-1);
+  }
+  void run(o2::framework::ProcessingContext& pc)
+  {
+    o2::fastsim::GeometryContainer geometryContainer; // Checking that the geometry files can be accessed and loaded
+    LOG(info) << "On-the-fly detector geometry provider running.";
+    if (detectorConfiguration.value.empty()) {
+      LOG(fatal) << "No detector configuration files provided.";
+      return;
+    }
+    int idx = 0;
+    for (auto& configFile : detectorConfiguration.value) {
+      LOG(info) << "Loading detector geometry from configuration file: " << configFile;
+      // If the filename starts with ccdb: then take the file from the ccdb
+      if (configFile.rfind("ccdb:", 0) == 0) {
+        std::string ccdbPath = configFile.substr(5); // remove "ccdb:" prefix
+        const std::string outPath = "/tmp/DetGeo/";
+        configFile = Form("%s/%s/snapshot.root", outPath.c_str(), ccdbPath.c_str());
+        std::ifstream checkFile(configFile); // Check if file already exists
+        if (!checkFile.is_open()) {          // File does not exist, retrieve from CCDB
+          LOG(info) << " --- CCDB source detected for detector geometry " << configFile;
+          std::map<std::string, std::string> metadata;
+          ccdb->getCCDBAccessor().retrieveBlob(ccdbPath, outPath, metadata, 1);
+          LOG(info) << " --- Now retrieving geometry configuration from CCDB to: " << configFile;
+        } else { // File exists, proceed to load
+          LOG(info) << " --- Geometry configuration file already exists: " << configFile << ". Skipping download.";
+          checkFile.close();
+        }
+        detectorConfiguration.value[idx] = configFile; // Update the filename to the local file
+      }
+      geometryContainer.addEntry(configFile);
+      idx++;
+    }
+    pc.services().get<o2::framework::ControlService>().endOfStream();
+    pc.services().get<o2::framework::ControlService>().readyToQuit(o2::framework::QuitRequest::Me);
+  }
+};
+
+// #define VERIFY_ALICE3
+#ifdef VERIFY_ALICE3
+struct OnTheFlyDetectorGeometryUser {
+  void init(o2::framework::InitContext& initContext)
+  {
+    o2::fastsim::GeometryContainer geometryContainer; // Checking that the configuration can be inherited
+    geometryContainer.init(initContext);
+  }
+  void run(o2::framework::ProcessingContext& pc)
+  {
+    pc.services().get<o2::framework::ControlService>().endOfStream();
+    pc.services().get<o2::framework::ControlService>().readyToQuit(o2::framework::QuitRequest::Me);
+  }
+};
+#endif
+
+o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)
+{
+  o2::framework::WorkflowSpec spec;
+  spec.push_back(adaptAnalysisTask<OnTheFlyDetectorGeometryProvider>(cfgc));
+#ifdef VERIFY_ALICE3
+  spec.push_back(adaptAnalysisTask<OnTheFlyDetectorGeometryUser>(cfgc));
+#endif
+  return spec;
+}