[ALICE3] Add method to specify dead phi region in layer

ALICE3/Core/DetLayer.cxx

@@ -16,11 +16,15 @@
 /// \brief  Basic struct to hold information regarding a detector layer to be used in fast simulation
 ///
 
-#include <vector>
-#include <string>
-
 #include "DetLayer.h"
 
+#include <CommonConstants/MathConstants.h>
+
+#include <fairlogger/Logger.h>
+
+#include <string>
+#include <vector>
+
 namespace o2::fastsim
 {
 
@@ -51,6 +55,62 @@
 {
 }
 
+void DetLayer::addDeadPhiRegion(float phiStart, float phiEnd)
+{
+  static constexpr float kDefaultValue = 2.f;
+  static constexpr float kPhiTolerance = 1e-4f;
+  if (mDeadPhiRegions == nullptr) {
+    mDeadPhiRegions = new TGraph();
+    mDeadPhiRegions->SetNameTitle(Form("deadPhiRegions_%s", name.Data()), Form("Dead phi regions for layer %s", name.Data()));
+    mDeadPhiRegions->AddPoint(0, kDefaultValue);
+    mDeadPhiRegions->AddPoint(o2::constants::math::TwoPI, kDefaultValue);
+  }
+  if (phiStart < 0 || phiStart >= o2::constants::math::TwoPI || phiEnd < 0 || phiEnd >= o2::constants::math::TwoPI) {
+    LOG(fatal) << "Cannot add dead phi region with invalid range [" << phiStart << ", " << phiEnd << "] to layer " << name;
+    return;
+  }
+  mDeadPhiRegions->AddPoint(phiStart, kDefaultValue);
+  mDeadPhiRegions->AddPoint(phiEnd, kDefaultValue);
+  mDeadPhiRegions->AddPoint(phiStart + kPhiTolerance, 0.f);
+  mDeadPhiRegions->AddPoint(phiEnd - kPhiTolerance, 0.f);
+  mDeadPhiRegions->Sort();
+}
+
+void DetLayer::setDeadPhiRegions(TGraph* graph)
+{
+  LOG(debug) << "Setting dead phi regions for layer " << name << " with graph " << (graph ? graph->GetName() : "nullptr");
+  if (mDeadPhiRegions != nullptr) {
+    LOG(warning) << "Overriding existing dead phi regions for layer " << name;
+    delete mDeadPhiRegions;
+  }
+  mDeadPhiRegions = graph;
+  if (mDeadPhiRegions->GetN() == 0) {
+    LOG(warning) << "Dead phi regions graph for layer " << name << " is empty, clearing dead regions";
+    mDeadPhiRegions = nullptr;
+    return; // cleared the dead regions
+  }
+  // Check sanity of the graph
+  if (mDeadPhiRegions != nullptr) {
+    for (int i = 0; i < mDeadPhiRegions->GetN(); i++) {
+      const float x = mDeadPhiRegions->GetX()[i];
+      const float y = mDeadPhiRegions->GetY()[i];
+      // First point has to be at 0, last point has to be at 2PI
+      if ((i == 0 && x != 0.f) || (i == mDeadPhiRegions->GetN() - 1 && x != o2::constants::math::TwoPI)) {
+        LOG(fatal) << "Dead phi regions graph for layer " << name << " has invalid x value " << x << " at point " << i << ", first point should be 0 and last point should be 2PI";
+      }
+      LOG(debug) << "Point " << i << ": (" << x << ", " << y << ")";
+      if (x < 0 || x > o2::constants::math::TwoPI) {
+        LOG(fatal) << "Dead phi regions graph for layer " << name << " has invalid x value " << x << " at point " << i;
+      }
+      if (y != 0.f && y != 2.f) {
+        LOG(fatal) << "Dead phi regions graph for layer " << name << " has invalid y value " << y << " at point " << i << ", should be 0 or 2";
+      }
+    }
+  } else {
+    LOG(info) << "Cleared dead phi regions for layer " << name;
+  }
+}
+
 std::string DetLayer::toString() const
 {
   std::string out = "";

ALICE3/Core/DetLayer.h

@@ -19,9 +19,10 @@
 #ifndef ALICE3_CORE_DETLAYER_H_
 #define ALICE3_CORE_DETLAYER_H_
 
-#include <string>
+#include <TGraph.h>
+#include <TString.h>
 
-#include "TString.h"
+#include <string>
 
 namespace o2::fastsim
 {
@@ -32,7 +33,7 @@
   DetLayer() = default;
   // Parametric constructor
   DetLayer(const TString& name_, float r_, float z_, float x0_, float xrho_,
            float resRPhi_ = 0.0f, float resZ_ = 0.0f, float eff_ = 0.0f, int type_ = layerInert);
   // Copy constructor
   DetLayer(const DetLayer& other);
 
@@ -47,6 +48,17 @@
   void setEfficiency(float eff_) { eff = eff_; }
   void setType(int type_) { type = type_; }
 
+  // Dead areas
+
+  /// @brief Add a dead region in phi for this layer
+  /// @param phiStart starting angle in radians of the dead region
+  /// @param phiEnd ending angle in radians of the dead region
+  void addDeadPhiRegion(float phiStart, float phiEnd);
+
+  /// @brief Set the dead regions in phi for this layer with a TGraph containing all regions. The graph should have y=2 for dead regions and y=0 for alive regions.
+  /// @param graph graph of the dead regions. Can be nullptr to clear the dead regions.
+  void setDeadPhiRegions(TGraph* graph);
+
   // Getters
   float getRadius() const { return r; }
   float getZ() const { return z; }
@@ -57,6 +69,7 @@
   float getEfficiency() const { return eff; }
   int getType() const { return type; }
   const TString& getName() const { return name; }
+  const TGraph* getDeadPhiRegions() const { return mDeadPhiRegions; }
 
   // Check layer type
   bool isInert() const { return type == layerInert; }
@@ -70,6 +83,15 @@
     os << layer.toString();
     return os;
   }
+  /// @brief Check if a given phi angle is in a dead region
+  /// @param phi The phi angle to check
+  /// @return True if the phi angle is in a dead region, false otherwise
+  bool isInDeadPhiRegion(float phi) const
+  {
+    if (mDeadPhiRegions == nullptr)
+      return false;
+    return mDeadPhiRegions->Eval(phi) > 1.f;
+  };
 
  private:
   // TString for holding name
@@ -90,11 +112,14 @@
   // efficiency
   float eff; // detection efficiency
 
+  // dead regions in phi (in radians)
+  TGraph* mDeadPhiRegions = nullptr;
+
   // layer type
   int type;                              // 0: undefined/inert, 1: silicon, 2: gas/tpc
   static constexpr int layerInert = 0;   // inert/undefined layer
   static constexpr int layerSilicon = 1; // silicon layer
   static constexpr int layerGas = 2;     // gas/tpc layer
 };
 
 } // namespace o2::fastsim

ALICE3/Core/FastTracker.cxx

@@ -22,7 +22,6 @@
 #include <TObject.h>
 
 #include <fstream>
-#include <map>
 #include <string>
 #include <vector>
 
@@ -33,7 +32,7 @@
 
 // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
 
-void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
+DetLayer* FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
 {
   LOG(debug) << "Adding layer " << name << " r=" << r << " z=" << z << " x0=" << x0 << " xrho=" << xrho << " resRPhi=" << resRPhi << " resZ=" << resZ << " eff=" << eff << " type=" << type;
   DetLayer newLayer(name, r, z, x0, xrho, resRPhi, resZ, eff, type);
@@ -53,6 +52,18 @@
   }
   // Add the new layer to the layers vector
   layers.push_back(newLayer);
+  // Return the last added layer
+  return &layers.back();
+}
+
+void FastTracker::addDeadPhiRegionInLayer(const std::string& layerName, float phiStart, float phiEnd)
+{
+  const int layerIdx = GetLayerIndex(layerName);
+  if (layerIdx < 0) {
+    LOG(fatal) << "Cannot add dead phi region to non-existing layer " << layerName;
+    return;
+  }
+  layers[layerIdx].addDeadPhiRegion(phiStart, phiEnd);
 }
 
 DetLayer FastTracker::GetLayer(int layer, bool ignoreBarrelLayers) const
@@ -185,7 +196,7 @@
   AddLayer("B03", 7., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
   AddLayer("B04", 9., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
   AddLayer("B05", 12., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
-  AddLayer("iTOF", 19, 250, x0iTOF, xrhoiTOF, resRPhiOT, resZOT, 0.0f, 0);
+  AddLayer("iTOF", 19, 250, x0iTOF, xrhoiTOF, resRPhiOT, resZOT, eff, 0);
   AddLayer("B06", 20., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
   AddLayer("B07", 30., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
   AddLayer("B08", 45., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1);
@@ -225,8 +236,8 @@
   for (int i = 0; i < kNPassiveBound; i++) {
     AddLayer(Form("tpc_boundary%d", i), rBoundary[i], zLength, radLBoundary[i], xrhoBoundary[i], 0); // dummy errors
   }
   for (Int_t k = 0; k < tpcRows; k++) {
     Float_t rowRadius = 0;
     if (k < innerRows)
       rowRadius = rowOneRadius + k * tpcInnerRadialPitch;
     else if (k >= innerRows && k < (innerRows + middleRows))
@@ -285,12 +296,15 @@
   for (const auto& layer : layers) {
     LOG(info) << " Reading layer " << layer;
 
-    auto getKey = [&layer, &env](const std::string& name) {
+    auto getKey = [&layer, &env](const std::string& name, const bool required = true) {
       std::string key = layer + "." + name;
       if (!env.Defined(key.c_str())) {
-        LOG(warning) << "Key " << key << " not defined in configuration file, getting the default value";
+        if (required) {
+          LOG(fatal) << "Key " << key << " not defined in configuration file";
+        }
+        LOG(debug) << "Key " << key << " not defined in configuration file, getting the default value";
       }
-      LOG(info) << " Getting key " << key;
+      LOG(debug) << " Getting key " << key << " from configuration file";
       return key;
     };
     const float r = env.GetValue(getKey("r").c_str(), -1.0f);
@@ -302,9 +316,36 @@
     const float resZ = env.GetValue(getKey("resZ").c_str(), 0.0f);
     const float eff = env.GetValue(getKey("eff").c_str(), 0.0f);
     const int type = env.GetValue(getKey("type").c_str(), 0);
+    const char* deadPhiRegions = env.GetValue(getKey("deadPhiRegions", false).c_str(), "");
 
     // void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
-    AddLayer(layer.c_str(), r, z, x0, xrho, resRPhi, resZ, eff, type);
+    LOG(info) << " Adding layer " << layer << " r=" << r << " z=" << z << " x0=" << x0 << " xrho=" << xrho << " resRPhi=" << resRPhi << " resZ=" << resZ << " eff=" << eff << " type=" << type << " deadPhiRegions=" << deadPhiRegions;
+
+    DetLayer* addedLayer = AddLayer(layer.c_str(), r, z, x0, xrho, resRPhi, resZ, eff, type);
+    if (strlen(deadPhiRegions) > 0) { // Taking it as ccdb path or local file
+                                      // Check if it begins with ccdb:
+      if (std::string(deadPhiRegions).rfind("ccdb:", 0) == 0) {
+        std::string ccdbPath = std::string(deadPhiRegions).substr(5); // remove "ccdb:" prefix
+        if (ccdbManager == nullptr) {
+          LOG(fatal) << "CCDB manager is null, cannot retrieve file " << ccdbPath;
+          return;
+        }
+        TGraph* g = ccdbManager->getForTimeStamp<TGraph>(ccdbPath, -1);
+        addedLayer->setDeadPhiRegions(g);
+      } else {
+        // Taking it as local file
+        TFile infile(deadPhiRegions, "READ");
+        if (!infile.IsOpen()) {
+          LOG(fatal) << "Cannot open dead phi regions file " << deadPhiRegions;
+          return;
+        }
+        TGraph* g = (TGraph*)infile.Get(infile.GetListOfKeys()->At(0)->GetName());
+        infile.Close();
+        addedLayer->setDeadPhiRegions(g);
+      }
+    } else {
+      LOG(debug) << " No dead phi regions for layer " << layer;
+    }
   }
 }
 
@@ -493,6 +534,11 @@
       continue; // inert layer, skip
     }
 
+    if (layers[il].isInDeadPhiRegion(inputTrack.getPhi())) {
+      LOGF(debug, "Track is in dead region of layer %d", il);
+      continue; // dead region, skip
+    }
+
     // layer is reached
     if (firstLayerReached < 0) {
       LOGF(debug, "First layer reached: %d", il);
@@ -656,7 +702,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]));
         }
       }
     }
@@ -694,7 +740,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
@@ -707,7 +753,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

@@ -41,7 +41,13 @@ class FastTracker
   virtual ~FastTracker() {}
 
   // Layer and layer configuration
-  void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
+  DetLayer* AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
+
+  /// Add a dead region in phi for a specific layer
+  /// \param layerName Name of the layer to modify
+  /// \param phiStart Start angle of the dead region (in radians)
+  /// \param phiEnd End angle of the dead region (in radians)
+  void addDeadPhiRegionInLayer(const std::string& layerName, float phiStart, float phiEnd);
   DetLayer GetLayer(const int layer, bool ignoreBarrelLayers = true) const;
   std::vector<DetLayer> GetLayers() const { return layers; }
   int GetLayerIndex(const std::string& name) const;

ALICE3/macros/testFastTracker.cxx → ALICE3/macros/testFastTracker.C

@@ -15,12 +15,16 @@
 
 #include "ALICE3/Core/FastTracker.h"
 
-#include "CCDB/BasicCCDBManager.h"
-#include "DataFormatsParameters/GRPLHCIFData.h"
+#include <CCDB/BasicCCDBManager.h>
+#include <DataFormatsParameters/GRPLHCIFData.h>
+
+#include <fairlogger/Logger.h>
 
 void testFastTracker(std::string geometryFile = "a3geo.ini")
 {
 
+  fair::Logger::SetConsoleSeverity(fair::Severity::debug);
+
   // auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
   // ccdb.setURL("http://alice-ccdb.cern.ch");
   o2::fastsim::FastTracker fastTracker;