Reads TF files and prints rate per BC

Author: Ella Taylor

Detectors/CTP/workflow/CMakeLists.txt

@@ -19,7 +19,8 @@ o2_add_library(CTPWorkflow
                                      O2::DetectorsRaw
                                      O2::Algorithm
                                      O2::CTPReconstruction
-                                     O2::CTPWorkflowIO)
+                                     O2::CTPWorkflowIO
+                                     O2::DataFormatsParameters)
 o2_add_executable(reco-workflow
                   COMPONENT_NAME ctp
                   SOURCES src/ctp-raw-decoder.cxx

Detectors/CTP/workflowLumi/include/CTPWorkflowLumi/RawDecoderSpec.h

@@ -52,7 +52,13 @@ class RawDecoderSpec : public framework::Task
   /// Output HW errors: {"CTP", "RAWHWERRORS", 0, Lifetime::Timeframe} -later
   void run(framework::ProcessingContext& ctx) final;
   void updateTimeDependentParams(framework::ProcessingContext& pc);
-
+  /// \brief Compute per BC luminosity from the interaction counts from CTP digits
+  /// \param ctpdigits Vector of CTP digits to be processed
+  /// \return Array of luminosity values for each BC
+  std::pair<std::array<double, o2::constants::lhc::LHCMaxBunches>, std::array<double, o2::constants::lhc::LHCMaxBunches>> computeLumiPerBC(const o2::pmr::vector<CTPDigit>& ctpdigits);
+  /// \brief Integrate luminosity per BC over multiple time frames
+  /// \param perTF Array of luminosity values for each BC from a single time frame
+  void integrateLumi(const std::array<double, o2::constants::lhc::LHCMaxBunches>& perTFInp1, const std::array<double, o2::constants::lhc::LHCMaxBunches>& perTFInp2);
  protected:
  private:
   // for digits
@@ -85,6 +91,13 @@ class RawDecoderSpec : public framework::Task
   std::array<uint64_t, o2::ctp::CTP_NCLASSES> mClsA{};
   std::array<uint64_t, o2::ctp::CTP_NCLASSES> mClsB{}; // from inputs
   bool mCheckConsistency = false;
+  std::array<double, o2::constants::lhc::LHCMaxBunches> mCountsPerBC1{};
+  std::array<double, o2::constants::lhc::LHCMaxBunches> mCountsPerBC2{};
+  double totalTime = 0.0;
+  const double orbitsPerTF = 32;
+  const double tfTime = orbitsPerTF * o2::constants::lhc::LHCOrbitMUS * 1e-6; // total time in seconds for one timeframe
+  std::bitset<3564> mLHCBCs;
+  const double timeInterval = o2::constants::lhc::LHCOrbitMUS * 1e-6; // one HBF
 };
 
 /// \brief Creating DataProcessorSpec for the CTP

Detectors/CTP/workflowLumi/src/RawDecoderSpec.cxx

@@ -21,6 +21,8 @@
 #include "DataFormatsCTP/TriggerOffsetsParam.h"
 #include "Framework/CCDBParamSpec.h"
 #include "DataFormatsCTP/Configuration.h"
+#include "CommonConstants/LHCConstants.h"
+#include <DataFormatsParameters/GRPLHCIFData.h>
 
 using namespace o2::ctp::reco_workflow;
 
@@ -181,6 +183,9 @@ void RawDecoderSpec::run(framework::ProcessingContext& ctx)
     mErrorTCR += mDecoder.getErrorTCR();
     mIRRejected += mDecoder.getIRRejected();
     mTCRRejected += mDecoder.getTCRRejected();
+    // Luminosity per bunch crossing
+    const auto [countsPerBC1, countsPerBC2] = computeLumiPerBC(mOutputDigits);
+    integrateLumi(countsPerBC1, countsPerBC2);
   }
   if (mDoLumi) {
     uint32_t tfCountsT = 0;
@@ -224,6 +229,89 @@ void RawDecoderSpec::run(framework::ProcessingContext& ctx)
     ctx.outputs().snapshot(o2::framework::Output{"CTP", "LUMI", 0}, mOutputLumiInfo);
   }
 }
+// Function to compute luminosity per BC from the interaction counts from CTP digits
+std::pair<std::array<double, o2::constants::lhc::LHCMaxBunches>, std::array<double, o2::constants::lhc::LHCMaxBunches>> RawDecoderSpec::computeLumiPerBC(const o2::pmr::vector<CTPDigit>& ctpdigits)
+{
+  int inp1 = mOutputLumiInfo.inp1;
+  int inp2 = mOutputLumiInfo.inp2;
+
+  uint64_t inputMask1 = 1ull << (inp1 - 1); // TVX
+  uint64_t inputMask2 = 1ull << (inp2 - 1); // VBA
+
+  double integratedRate = 0.0;
+  std::cout << "Lumi called" << std::endl;
+
+  std::array<double, o2::constants::lhc::LHCMaxBunches> countsPerBC1{};
+  std::array<double, o2::constants::lhc::LHCMaxBunches> countsPerBC2{};
+
+  for (const auto& digit : ctpdigits) {
+    uint64_t mask = digit.CTPInputMask.to_ullong();
+    uint16_t bc = digit.intRecord.bc;
+
+    if (bc < o2::constants::lhc::LHCMaxBunches) {
+      if (mask & inputMask1) {
+        countsPerBC1[bc] += 1.0;
+        integratedRate += 1.0;
+        //  std::cout << "Orbit: " << std::dec << digit.intRecord.orbit << " Orbit: 0x" << std::hex << digit.intRecord.orbit << std::endl;
+      }
+
+      if (mask & inputMask2) {
+        countsPerBC2[bc] += 1.0;
+      }
+    }
+  }
+  totalTime += tfTime; // Accumulate total time for all processed time frames
+  for (size_t bc = 0; bc < countsPerBC1.size(); ++bc) {
+    if (countsPerBC1[bc] > 0) {
+      // LOG(info) << " BC " << bc << ": " << lumiPerBC[bc]/totalTime;
+    }
+  }
+  // std::cout << "Integrated luminosity over all BCs: " << integratedRate/totalTime << std::endl;
+  return {countsPerBC1, countsPerBC2};
+}
+// Accumulate luminosity per BC over multiple time frames
+void RawDecoderSpec::integrateLumi(const std::array<double, o2::constants::lhc::LHCMaxBunches>& perTFInp1,
+                                   const std::array<double, o2::constants::lhc::LHCMaxBunches>& perTFInp2)
+{
+
+  for (size_t bc = 0; bc < mCountsPerBC1.size(); ++bc) {
+    mCountsPerBC1[bc] += perTFInp1[bc];
+  }
+
+  for (size_t bc = 0; bc < mCountsPerBC2.size(); ++bc) {
+    mCountsPerBC2[bc] += perTFInp2[bc];
+  }
+
+  // Count number of filled BCs
+  size_t filledBCs = mLHCBCs.count();
+
+  for (size_t bc = 0; bc < mCountsPerBC1.size(); ++bc) {
+    if (mLHCBCs.test(bc)) { // Only print filled BCs
+      LOG(info) << " Filled BC " << bc
+                << ": Input1 Lumi: " << mCountsPerBC1[bc] / (totalTime * filledBCs)
+                << ", Input2 Lumi: " << mCountsPerBC2[bc] / (totalTime * filledBCs)
+                << "; Accumulated Counts Input1: " << mCountsPerBC1[bc]
+                << ", Input2: " << mCountsPerBC2[bc];
+    } else if (mCountsPerBC1[bc] > 0) { // Only print non-zero luminosity
+      LOG(info) << " BC " << bc
+                << ": Input1 Lumi: " << mCountsPerBC1[bc] / (totalTime * filledBCs)
+                << ", Input2 Lumi: " << mCountsPerBC2[bc] / (totalTime * filledBCs)
+                << "; Accumulated Counts Input1: " << mCountsPerBC1[bc]
+                << ", Input2: " << mCountsPerBC2[bc];
+    }
+  }
+  // Calculate and print the total integrated luminosity
+  int totalCountsInp1 = 0;
+  int totalCountsInp2 = 0;
+  for (const auto& count : mCountsPerBC1) {
+    totalCountsInp1 += count;
+  }
+  LOG(info) << "Total Integrated Luminosity Input 1: " << totalCountsInp1 / totalTime;
+  for (const auto& count : mCountsPerBC2) {
+    totalCountsInp2 += count;
+  }
+  LOG(info) << "Total Integrated Luminosity Input 2: " << totalCountsInp2 / totalTime;
+}
 o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool askDISTSTF, bool digits, bool lumi)
 {
   if (!digits && !lumi) {
@@ -237,6 +325,7 @@ o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool
 
   std::vector<o2::framework::OutputSpec> outputs;
   inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("CTP/Config/Config", 1));
+  inputs.emplace_back("grplhcif", "GLO", "GRPLHCIF", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("GLO/Config/GRPLHCIF"));
   inputs.emplace_back("trigoffset", "CTP", "Trig_Offset", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("CTP/Config/TriggerOffsets"));
   if (digits) {
     outputs.emplace_back("CTP", "DIGITS", 0, o2::framework::Lifetime::Timeframe);
@@ -245,7 +334,7 @@ o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool
     outputs.emplace_back("CTP", "LUMI", 0, o2::framework::Lifetime::Timeframe);
   }
   return o2::framework::DataProcessorSpec{
-    "ctp-raw-decoder",
+    "ctp-raw-decoder-lumi",
     inputs,
     outputs,
     o2::framework::AlgorithmSpec{o2::framework::adaptFromTask<o2::ctp::reco_workflow::RawDecoderSpec>(digits, lumi)},
@@ -271,5 +360,17 @@ void RawDecoderSpec::updateTimeDependentParams(framework::ProcessingContext& pc)
       mDecoder.setCTPConfig(*ctpcfg);
       LOG(info) << "ctpconfig for run done:" << mDecoder.getCTPConfig().getRunNumber();
     }
+    const auto grplhcif = pc.inputs().get<o2::parameters::GRPLHCIFData*>("grplhcif");
+    if (grplhcif != nullptr) {
+      LOG(info) << "GRPLHCIF injection scheme: " << grplhcif->getInjectionScheme();
+
+      // Get filled bunches
+      auto bfilling = grplhcif->getBunchFilling();
+      std::vector<int> bcs = bfilling.getFilledBCs();
+      mLHCBCs.reset();
+      for (auto const& bc : bcs) {
+        mLHCBCs.set(bc, 1);
+      }
+    }
   }
 }