[PWGCF] Add lite colshape table (#17410)

Author: ariedel-cern

PWGCF/Femto/Core/cascadeBuilder.h

@@ -793,15 +793,15 @@ struct CascadeBuilderDerivedToDerivedProducts : o2::framework::ProducesGroup {
 class CascadeBuilderDerivedToDerived
 {
  public:
-  CascadeBuilderDerivedToDerived() = default;
-  ~CascadeBuilderDerivedToDerived() = default;
-
   template <typename T>
   void init(T& config)
   {
     mLimitXi = config.limitXi.value;
     mLimitOmega = config.limitOmega.value;
 
+    if (mLimitXi < 0 || mLimitOmega < 0) {
+      LOG(fatal) << "Cascade limits must be non-negative (got " << mLimitXi << " and " << mLimitOmega << "). Breaking...";
+    }
     if (mLimitXi == 0 && mLimitOmega == 0) {
       LOG(fatal) << "Both xi limit and omega limit are 0. Breaking...";
     }
@@ -810,34 +810,40 @@ class CascadeBuilderDerivedToDerived
   template <typename T1, typename T2, typename T3, typename T4>
   bool collisionHasTooFewXis(T1 const& col, T2 const& /*xiTable*/, T3& partitionXi, T4& cache)
   {
+    if (mLimitXi == 0) { // xis disabled, cannot reject on them
+      return false;
+    }
     auto xiSlice = partitionXi->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
-    return xiSlice.size() < mLimitXi;
+    return xiSlice.size() < static_cast<int64_t>(mLimitXi);
   }
 
   template <typename T1, typename T2, typename T3, typename T4>
   bool collisionHasTooFewOmegas(T1 const& col, T2 const& /*omegaTable*/, T3& partitionOmega, T4& cache)
   {
+    if (mLimitOmega == 0) { // omegas disabled, cannot reject on them
+      return false;
+    }
     auto omegaSlice = partitionOmega->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
-    return omegaSlice.size() < mLimitOmega;
+    return omegaSlice.size() < static_cast<int64_t>(mLimitOmega);
   }
 
   template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
   void processXis(T1 const& col, T2 const& /*xiTable*/, T3 const& oldTrackTable, T4& partitionXi, T5& trackBuilder, T6& cache, T7& newXiTable, T8& newTrackTable, T9& newCollisionTable)
   {
+    if (mLimitXi == 0) { // xis disabled
+      return;
+    }
+
     auto xiSlice = partitionXi->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
 
     for (auto const& xi : xiSlice) {
+      auto bachelor = oldTrackTable.rawIteratorAt(utils::daughterRow(xi.bachelorId(), oldTrackTable));
+      auto posDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(xi.posDauId(), oldTrackTable));
+      auto negDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(xi.negDauId(), oldTrackTable));
 
-      // auto bachelor = xi.template bachelor_as<T3>();
-      // auto posDaughter = xi.template posDau_as<T3>();
-      // auto negDaughter = xi.template negDau_as<T3>();
-      auto bachelor = oldTrackTable.rawIteratorAt(xi.bachelorId() - oldTrackTable.offset());
-      auto posDaughter = oldTrackTable.rawIteratorAt(xi.posDauId() - oldTrackTable.offset());
-      auto negDaughter = oldTrackTable.rawIteratorAt(xi.negDauId() - oldTrackTable.offset());
-
-      int bachelorIndex = trackBuilder.getDaughterIndex(bachelor, newTrackTable, newCollisionTable);
-      int posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable);
-      int negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable);
+      int64_t bachelorIndex = trackBuilder.getDaughterIndex(bachelor, newTrackTable, newCollisionTable);
+      int64_t posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable);
+      int64_t negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable);
 
       newXiTable.producedXis(newCollisionTable.producedCollision.lastIndex(),
                              xi.signedPt(),
@@ -854,20 +860,20 @@ class CascadeBuilderDerivedToDerived
   template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
   void processOmegas(T1 const& col, T2 const& /*omegaTable*/, T3 const& oldTrackTable, T4& partitionOmega, T5& trackBuilder, T6& cache, T7& newOmegaTable, T8& newTrackTable, T9& newCollisionTable)
   {
+    if (mLimitOmega == 0) { // omegas disabled
+      return;
+    }
+
     auto omegaSlice = partitionOmega->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
 
     for (auto const& omega : omegaSlice) {
+      auto bachelor = oldTrackTable.rawIteratorAt(utils::daughterRow(omega.bachelorId(), oldTrackTable));
+      auto posDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(omega.posDauId(), oldTrackTable));
+      auto negDaughter = oldTrackTable.rawIteratorAt(utils::daughterRow(omega.negDauId(), oldTrackTable));
 
-      // auto bachelor = omega.template bachelor_as<T3>();
-      // auto posDaughter = omega.template posDau_as<T3>();
-      // auto negDaughter = omega.template negDau_as<T3>();
-      auto bachelor = oldTrackTable.rawIteratorAt(omega.bachelorId() - oldTrackTable.offset());
-      auto posDaughter = oldTrackTable.rawIteratorAt(omega.posDauId() - oldTrackTable.offset());
-      auto negDaughter = oldTrackTable.rawIteratorAt(omega.negDauId() - oldTrackTable.offset());
-
-      int bachelorIndex = trackBuilder.getDaughterIndex(bachelor, newTrackTable, newCollisionTable);
-      int posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable);
-      int negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable);
+      int64_t bachelorIndex = trackBuilder.getDaughterIndex(bachelor, newTrackTable, newCollisionTable);
+      int64_t posDaughterIndex = trackBuilder.getDaughterIndex(posDaughter, newTrackTable, newCollisionTable);
+      int64_t negDaughterIndex = trackBuilder.getDaughterIndex(negDaughter, newTrackTable, newCollisionTable);
 
       newOmegaTable.producedOmegas(newCollisionTable.producedCollision.lastIndex(),
                                    omega.signedPt(),

PWGCF/Femto/Core/collisionBuilder.h

@@ -291,7 +291,7 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
         {collisionFilterNames.at(kFilterSphericityMax), mSphericityMax},
         {collisionFilterNames.at(kFilterRctFlags), mUseRctFlags ? 1.f : 0.f},
       });
-  };
+  }
 
   /// \brief Initialize the Zorro trigger machinery for a new run. No-op if no triggers configured.
   template <typename T1, typename T2>
@@ -320,7 +320,7 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
     mMagField = MagField;
   }
 
-  float getMagneticField()
+  [[nodiscard]] int getMagneticField()
   {
     return mMagField;
   }
@@ -491,14 +491,14 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
     }
 
     this->assembleBitmask<SelectionHistName>();
-  };
+  }
 
  protected:
   template <typename T>
   float computeSphericity(T const& tracks)
   {
-    int minNumberTracks = 2;
-    double maxSphericity = 2.f;
+    const int64_t minNumberTracks = 2;
+    const double maxSphericity = 2.f;
     if (tracks.size() <= minNumberTracks) {
       return maxSphericity;
     }
@@ -529,7 +529,7 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
 
   // filter cuts
   float mVtxZMin = -12.f;
-  float mVtxZMax = -12.f;
+  float mVtxZMax = 12.f;
   float mSphericityMin = 0.f;
   float mSphericityMax = 1.f;
   float mMagFieldMin = -5.f;
@@ -539,7 +539,7 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
   float mCentMin = 0.f;
   float mCentMax = 100.f;
 
-  int mMagField = 0.f;
+  int mMagField = 0;
   float mSphericity = 0.f;
   float mCentrality = 0.f;
   float mMultiplicity = 0.f;
@@ -568,6 +568,7 @@ struct CollisionBuilderProducts : o2::framework::ProducesGroup {
   o2::framework::Produces<o2::aod::FColMults> producedMultiplicityEstimators;
   o2::framework::Produces<o2::aod::FColCents> producedCentralityEstimators;
   o2::framework::Produces<o2::aod::FColShapes> producedShapes;
+  o2::framework::Produces<o2::aod::FLiteColShapes> producedLiteShapes;
 };
 
 struct ConfCollisionTables : o2::framework::ConfigurableGroup {
@@ -580,6 +581,7 @@ struct ConfCollisionTables : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<int> produceMults{"produceMults", -1, "Produce Multiplicities (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> produceCents{"produceCents", -1, "Produce Centralities (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> produceShapes{"produceShapes", -1, "Produce Event shape variables (-1: auto; 0 off; 1 on)"};
+  o2::framework::Configurable<int> produceLiteShapes{"produceLiteShapes", -1, "Produce Lite Event shape variables (-1: auto; 0 off; 1 on)"};
 };
 
 template <auto& SelectionHistName, auto& FilterHistName>
@@ -606,16 +608,22 @@ class CollisionBuilder
     mProducedMultiplicities = utils::enableTable("FColMults_001", confTable.produceMults.value, initContext);
     mProducedCentralities = utils::enableTable("FColCents_001", confTable.produceCents.value, initContext);
     mProducedShapes = utils::enableTable("FColShapes_001", confTable.produceShapes.value, initContext);
+    mProducedLiteShapes = utils::enableTable("FLiteColShapes_001", confTable.produceLiteShapes.value, initContext);
 
     if (mProducedCollisions && mProducedLiteCollisions) {
       LOG(fatal) << "FCols and FLiteCols are mutually exclusive -- enable only one. "
                  << "FLiteCols is meant to only replace FCols at the producer stage (for better compression in derived data); "
                  << "use the dedicated converter task to reconstruct FCols from FLiteCols downstream.";
     }
+    if (mProducedShapes && mProducedLiteShapes) {
+      LOG(fatal) << "FColShapes and FLiteColShapes are mutually exclusive -- enable only one. "
+                 << "FLiteColShapes is meant to only replace FColShapes at the producer stage (for better compression in derived data); "
+                 << "use the dedicated converter task to reconstruct FColShapes from FLiteColShapes downstream.";
+    }
 
     if (mProducedCollisions || mProducedLiteCollisions || mProducedCollisionMasks ||
         mProducedPositions || mProducedSphericities || mProducedMultiplicities ||
-        mProducedCentralities) {
+        mProducedCentralities || mProducedShapes || mProducedLiteShapes) {
       mFillAnyTable = true;
     } else {
       LOG(info) << "No tables configured, Selection object will not be configured...";
@@ -634,15 +642,14 @@ class CollisionBuilder
     if (mRunNumber != bc.runNumber()) {
       mRunNumber = bc.runNumber();
       if (mMagFieldForced == 0) {
-        static o2::parameters::GRPMagField* grpo = nullptr;
-        LOG(info) << "Get magentic field with Path: " << mGrpPath << "; Run number: " << mRunNumber;
-        grpo = ccdb->template getForRun<o2::parameters::GRPMagField>(mGrpPath, mRunNumber);
+        o2::parameters::GRPMagField* grpo = ccdb->template getForRun<o2::parameters::GRPMagField>(mGrpPath, mRunNumber);
+        LOG(info) << "Get magnetic field with Path: " << mGrpPath << "; Run number: " << mRunNumber;
         if (grpo == nullptr) {
           LOG(fatal) << "GRP object not found for Run " << mRunNumber;
         }
         mMagField = static_cast<int>(grpo->getNominalL3Field()); // get magnetic field in kG
       } else {
-        LOG(info) << "Force magentic field to " << mMagFieldForced << "kG";
+        LOG(info) << "Force magnetic field to " << mMagFieldForced << "kG";
         mMagField = mMagFieldForced;
       }
 
@@ -752,6 +759,12 @@ class CollisionBuilder
         mCollisionSelection.getEventPlane());
     }
 
+    if (mProducedLiteShapes) {
+      collisionProducts.producedLiteShapes(
+        o2::aod::femtocollisions::lite::binQvec(mCollisionSelection.getQvector()),
+        o2::aod::femtocollisions::lite::binEventPlaneAngle(mCollisionSelection.getEventPlane()));
+    }
+
     mCollisionAlreadyFilled = true;
   }
 
@@ -795,6 +808,7 @@ class CollisionBuilder
   bool mProducedMultiplicities = false;
   bool mProducedCentralities = false;
   bool mProducedShapes = false;
+  bool mProducedLiteShapes = false;
 };
 
 struct CollisionBuilderDerivedToDerivedProducts : o2::framework::ProducesGroup {

PWGCF/Femto/Core/femtoUtils.h

@@ -47,6 +47,20 @@ inline std::optional<T2> getIndex(const T1& index, const std::unordered_map<T1,
   return std::nullopt;
 }
 
+/// Translate a global daughter index into a row of the current track table frame.
+/// Aborts if the index does not fall inside the frame, which would otherwise
+/// silently produce an out-of-range iterator.
+template <typename T>
+int64_t daughterRow(int64_t daughterId, T const& trackTable)
+{
+  const int64_t row = daughterId - trackTable.offset();
+  if (daughterId < 0 || row < 0 || row >= static_cast<int64_t>(trackTable.size())) {
+    LOG(fatal) << "Daughter index " << daughterId << " out of range for track table (offset "
+               << trackTable.offset() << ", size " << trackTable.size() << "). Breaking...";
+  }
+  return row;
+}
+
 template <typename T>
 float itsSignal(T const& track)
 {
@@ -243,7 +257,7 @@ inline int signum(T x)
 }
 
 template <typename T>
-inline T binLinear(float value, float lo, float hi, float step)
+T binLinear(float value, float lo, float hi, float step)
 {
   float v = std::clamp(value, lo, hi);
   auto idx = static_cast<int64_t>(std::round((v - lo) / step));
@@ -253,46 +267,46 @@ inline T binLinear(float value, float lo, float hi, float step)
 }
 
 template <typename T>
-inline float unBinLinear(T binned, float lo, float step)
+float unBinLinear(T binned, float lo, float step)
 {
   auto idx = static_cast<int64_t>(binned) - static_cast<int64_t>(std::numeric_limits<T>::min());
   return lo + static_cast<float>(idx) * step;
 }
 
 template <typename T>
-inline T binLogSigned(float signedValue, float magMin, float magMax)
+T binLogSigned(float signedValue, float magMin, float magMax)
 {
   static_assert(std::is_unsigned_v<T>, "binLogSigned requires an unsigned storage type");
-  constexpr uint32_t TotalBits = sizeof(T) * 8;
-  constexpr uint32_t HalfLevels = 1u << (TotalBits - 1);
-  uint32_t sign = (signedValue < 0.f) ? 1u : 0u;
-  float mag = std::clamp(std::fabs(signedValue), magMin, magMax);
-  float logLo = std::log(magMin);
-  float logHi = std::log(magMax);
-  float step = (logHi - logLo) / static_cast<float>(HalfLevels - 1);
-  auto idx = static_cast<uint32_t>(std::round((std::log(mag) - logLo) / step));
-  idx = std::clamp(idx, 0u, HalfLevels - 1);
-  return static_cast<T>((sign << (TotalBits - 1)) | idx);
+  constexpr uint64_t TotalBits = sizeof(T) * 8;
+  constexpr uint64_t HalfLevels = uint64_t{1} << (TotalBits - 1);
+  const uint64_t sign = (signedValue < 0.f) ? uint64_t{1} : uint64_t{0};
+  const float mag = std::clamp(std::fabs(signedValue), magMin, magMax);
+  const float logLo = std::log(magMin);
+  const float logHi = std::log(magMax);
+  const float step = (logHi - logLo) / static_cast<float>(HalfLevels - 1);
+  auto idx = static_cast<int64_t>(std::round((std::log(mag) - logLo) / step));
+  idx = std::clamp(idx, int64_t{0}, static_cast<int64_t>(HalfLevels - 1));
+  return static_cast<T>((sign << (TotalBits - 1)) | static_cast<uint64_t>(idx));
 }
 
 template <typename T>
-inline float unBinLogSigned(T binned, float magMin, float magMax)
+float unBinLogSigned(T binned, float magMin, float magMax)
 {
-  constexpr uint32_t TotalBits = sizeof(T) * 8;
-  constexpr uint32_t HalfLevels = 1u << (TotalBits - 1);
-  constexpr T SignMask = static_cast<T>(1u << (TotalBits - 1));
+  static_assert(std::is_unsigned_v<T>, "unBinLogSigned requires an unsigned storage type");
+  constexpr uint64_t TotalBits = sizeof(T) * 8;
+  constexpr uint64_t HalfLevels = uint64_t{1} << (TotalBits - 1);
+  constexpr T SignMask = static_cast<T>(uint64_t{1} << (TotalBits - 1));
   constexpr T MagMask = static_cast<T>(SignMask - 1);
-  float sign = (binned & SignMask) ? -1.f : 1.f;
-  uint32_t idx = binned & MagMask;
-  float logLo = std::log(magMin);
-  float logHi = std::log(magMax);
-  float step = (logHi - logLo) / static_cast<float>(HalfLevels - 1);
-  float mag = std::exp(logLo + static_cast<float>(idx) * step);
-  return sign * mag;
+  const float sign = (binned & SignMask) ? -1.f : 1.f;
+  const auto idx = static_cast<uint64_t>(binned & MagMask);
+  const float logLo = std::log(magMin);
+  const float logHi = std::log(magMax);
+  const float step = (logHi - logLo) / static_cast<float>(HalfLevels - 1);
+  return sign * std::exp(logLo + static_cast<float>(idx) * step);
 }
 
 template <typename T>
-inline int unBinSign(T binned)
+int unBinSign(T binned)
 {
   static_assert(std::is_unsigned_v<T>, "unBinSign requires an unsigned storage type");
   constexpr uint64_t TotalBits = sizeof(T) * 8;
@@ -301,7 +315,7 @@ inline int unBinSign(T binned)
 }
 
 template <typename T>
-inline T binLogUnsigned(float value, float magMin, float magMax)
+T binLogUnsigned(float value, float magMin, float magMax)
 {
   static_assert(std::is_unsigned_v<T>, "binLogUnsigned requires an unsigned storage type");
   constexpr uint64_t TotalBits = sizeof(T) * 8;
@@ -316,7 +330,7 @@ inline T binLogUnsigned(float value, float magMin, float magMax)
 }
 
 template <typename T>
-inline float unBinLogUnsigned(T binned, float magMin, float magMax)
+float unBinLogUnsigned(T binned, float magMin, float magMax)
 {
   constexpr uint64_t TotalBits = sizeof(T) * 8;
   constexpr uint64_t Levels = uint64_t{1} << TotalBits;