ctpdev: lumi per run added to fetcher (AliceO2Group#13945)

* dev: lumi calcalation in fetcher

* clang

Author: lietava

DataFormats/Detectors/CTP/include/DataFormatsCTP/CTPRateFetcher.h

@@ -33,15 +33,18 @@ class CTPRateFetcher
   void setupRun(int runNumber, o2::ccdb::BasicCCDBManager* ccdb, uint64_t timeStamp, bool initScalers);
   void updateScalers(ctp::CTPRunScalers& scalers);
   int getRates(std::array<double, 3>& rates, o2::ccdb::BasicCCDBManager* ccdb, int runNumber, const std::string sourceName); // rates at start,stop and middle of the run
-  void setOrbit(bool orb) { mOrbit = orb; }
-  void setOutsideLimits(bool qc) { mOutsideLimits = qc; }
+  double getLumi(o2::ccdb::BasicCCDBManager* ccdb, int runNumber, const std::string sourceName, int puCorr = 0);             // total lumi for a run
+  double getLumiNoPuCorr(const std::string& classname, int type = 1);
+  double getLumiWPuCorr(const std::string& classname, int type = 1);
+  void setOrbit(bool orb) { mOrbit = orb; }               // use orbit instead of time
+  void setOutsideLimits(bool qc) { mOutsideLimits = qc; } // return first/last rate of time outside of run
 
  private:
   double fetchCTPratesInputs(uint64_t timeStamp, int input);
   double fetchCTPratesClasses(uint64_t timeStamp, const std::string& className, int inputType = 1);
   double fetchCTPratesInputsNoPuCorr(uint64_t timeStamp, int input);
   double fetchCTPratesClassesNoPuCorr(uint64_t timeStamp, const std::string& className, int inputType = 1);
-
+  double getLumi(const std::string& classname, int type = 1, int puCorr = 0);
   double pileUpCorrection(double rate);
   int mRunNumber = -1;
   bool mOutsideLimits = 0;

DataFormats/Detectors/CTP/include/DataFormatsCTP/Scalers.h

@@ -143,6 +143,10 @@ class CTPRunScalers
       mScalerRecordO2[mScalerRecordO2.size() - 1].scalers[i].l1After - mScalerRecordO2[0].scalers[i].l1After,
     };
   }
+  /// retrieves integral - same interface as getRate, no pileup correction
+  uint64_t getLumiNoPuCorr(int classindex, int type) const;
+  /// retrieves vector of counters - same interface as getRate, needed for
+  std::vector<std::pair<double_t, double_t>> getRatesForIndex(int classindex, int type) const;
   /// retrieves time boundaries of this scaler object from O2 scalers
   std::pair<unsigned long, unsigned long> getTimeLimit() const
   {

DataFormats/Detectors/CTP/src/CTPRateFetcher.cxx

@@ -41,7 +41,7 @@ double CTPRateFetcher::fetchNoPuCorr(o2::ccdb::BasicCCDBManager* ccdb, uint64_t
     if (runNumber < 534202) {
       return fetchCTPratesClassesNoPuCorr(timeStamp, "minbias_TVX_L0", 3); // 2022
     } else {
-      double_t ret = fetchCTPratesClassesNoPuCorr(timeStamp, "CMTVX-B-NOPF");
+      double ret = fetchCTPratesClassesNoPuCorr(timeStamp, "CMTVX-B-NOPF");
       if (ret == -2.) {
         LOG(info) << "Trying different class";
         ret = fetchCTPratesClassesNoPuCorr(timeStamp, "CMTVX-NONE");
@@ -77,6 +77,94 @@ int CTPRateFetcher::getRates(std::array<double, 3>& rates, o2::ccdb::BasicCCDBMa
   rates[2] = rateM;
   return 0;
 }
+double CTPRateFetcher::getLumiNoPuCorr(const std::string& classname, int type)
+{
+  if (classname == "zncinp") {
+    return mScalers.getLumiNoPuCorr(26, 7);
+  }
+  std::vector<ctp::CTPClass>& ctpcls = mConfig.getCTPClasses();
+  std::vector<int> clslist = mConfig.getTriggerClassList();
+  int classIndex = -1;
+  for (size_t i = 0; i < clslist.size(); i++) {
+    if (ctpcls[i].name.find(classname) != std::string::npos) {
+      classIndex = i;
+      break;
+    }
+  }
+  if (classIndex == -1) {
+    LOG(warn) << "Trigger class " << classname << " not found in CTPConfiguration";
+    return -1;
+  }
+  return mScalers.getLumiNoPuCorr(classIndex, type);
+}
+double CTPRateFetcher::getLumiWPuCorr(const std::string& classname, int type)
+{
+  std::vector<std::pair<double, double>> scals;
+  if (classname == "zncinp") {
+    scals = mScalers.getRatesForIndex(26, 7);
+  } else {
+    std::vector<ctp::CTPClass>& ctpcls = mConfig.getCTPClasses();
+    std::vector<int> clslist = mConfig.getTriggerClassList();
+    int classIndex = -1;
+    for (size_t i = 0; i < clslist.size(); i++) {
+      if (ctpcls[i].name.find(classname) != std::string::npos) {
+        classIndex = i;
+        break;
+      }
+    }
+    if (classIndex == -1) {
+      LOG(warn) << "Trigger class " << classname << " not found in CTPConfiguration";
+      return -1;
+    }
+    scals = mScalers.getRatesForIndex(classIndex, type);
+  }
+  double lumi = 0;
+  for (auto const& ss : scals) {
+    // std::cout << ss.first << " " << ss.second << " " << pileUpCorrection(ss.first/ss.second) << std::endl;
+    lumi += pileUpCorrection(ss.first / ss.second) * ss.second;
+  }
+  return lumi;
+}
+double CTPRateFetcher::getLumi(const std::string& classname, int type, int puCorr)
+{
+  if (puCorr) {
+    return getLumiWPuCorr(classname, type);
+  } else {
+    return getLumiNoPuCorr(classname, type);
+  }
+}
+
+double CTPRateFetcher::getLumi(o2::ccdb::BasicCCDBManager* ccdb, int runNumber, const std::string sourceName, int puCorr)
+{
+  // setupRun(runNumber, ccdb, timeStamp, 1);
+  if (sourceName.find("ZNC") != std::string::npos) {
+    if (runNumber < 544448) {
+      return getLumi("zncinp", 1, puCorr) / (sourceName.find("hadronic") != std::string::npos ? 28. : 1.);
+    } else {
+      return getLumi("C1ZNC-B-NOPF-CRU", 6, puCorr) / (sourceName.find("hadronic") != std::string::npos ? 28. : 1.);
+    }
+  } else if (sourceName == "T0CE") {
+    return getLumi("CMTVXTCE-B-NOPF", 1, puCorr);
+  } else if (sourceName == "T0SC") {
+    return getLumi("CMTVXTSC-B-NOPF", 1, puCorr);
+  } else if (sourceName == "T0VTX") {
+    if (runNumber < 534202) {
+      return getLumi("minbias_TVX_L0", 3, puCorr); // 2022
+    } else {
+      double ret = getLumi("CMTVX-B-NOPF", 1, puCorr);
+      if (ret == -1.) {
+        LOG(info) << "Trying different class";
+        ret = getLumi("CMTVX-NONE", 1, puCorr);
+        if (ret < 0) {
+          LOG(fatal) << "None of the classes used for lumi found";
+        }
+      }
+      return ret;
+    }
+  }
+  LOG(error) << "CTP Lumi for " << sourceName << " not available";
+  return 0;
+}
 //
 double CTPRateFetcher::fetchCTPratesClasses(uint64_t timeStamp, const std::string& className, int inputType)
 {

DataFormats/Detectors/CTP/src/Scalers.cxx

@@ -657,7 +657,77 @@ void CTPRunScalers::printLMBRateVsT() const
     }
   }
 }
-
+//
+uint64_t CTPRunScalers::getLumiNoPuCorr(int classindex, int type) const
+{
+  if (type < 7) {
+    const auto s0 = mScalerRecordO2[0].scalers[classindex];
+    const auto s1 = mScalerRecordO2[mScalerRecordO2.size() - 1].scalers[classindex];
+    switch (type) {
+      case 1:
+        return (s1.lmBefore - s0.lmBefore);
+      case 2:
+        return (s1.lmAfter - s0.lmAfter);
+      case 3:
+        return (s1.l0Before - s0.l0Before);
+      case 4:
+        return (s1.l0After - s0.l0After);
+      case 5:
+        return (s1.l1Before - s0.l1Before);
+      case 6:
+        return (s1.l1After - s0.l1After);
+      default:
+        LOG(error) << "Wrong type:" << type;
+        return -1; // wrong type
+    }
+  } else if (type == 7) {
+    auto s0 = mScalerRecordO2[0].scalersInps[classindex]; // type CTPScalerO2*
+    auto s1 = mScalerRecordO2[mScalerRecordO2.size() - 1].scalersInps[classindex];
+    return (s1 - s0);
+  } else {
+    LOG(error) << "Wrong type:" << type;
+    return -1; // wrong type
+  }
+};
+//
+std::vector<std::pair<double_t, double_t>> CTPRunScalers::getRatesForIndex(int classindex, int type) const
+{
+  std::vector<std::pair<double_t, double_t>> scals;
+  for (int i = 0; i < mScalerRecordO2.size() - 1; i++) {
+    double_t diff = 0;
+    // double_t timeDiff = mScalerRecordO2[i + 1].epochTime -  mScalerRecordO2[i].epochTime;
+    double_t timeDiff = (mScalerRecordO2[i + 1].intRecord.orbit - mScalerRecordO2[i].intRecord.orbit) * o2::constants::lhc::LHCOrbitMUS / 1.e6;
+    if (type < 7) {
+      const auto s0 = mScalerRecordO2[i].scalers[classindex];
+      const auto s1 = mScalerRecordO2[i + 1].scalers[classindex];
+      if (type == 1) {
+        diff = s1.lmBefore - s0.lmBefore;
+      } else if (type == 2) {
+        diff = s1.lmAfter - s0.lmAfter;
+      } else if (type == 3) {
+        diff = s1.l0Before - s0.l0Before;
+      } else if (type == 4) {
+        diff = s1.l0After - s0.l0After;
+      } else if (type == 5) {
+        diff = s1.l1Before - s0.l1Before;
+      } else if (type == 6) {
+        diff = s1.l1After - s0.l1After;
+      } else {
+        LOG(error) << "Wrong type:" << type;
+        return scals; // wrong type
+      }
+    } else if (type == 7) {
+      auto s0 = mScalerRecordO2[i].scalersInps[classindex]; // type CTPScalerO2*
+      auto s1 = mScalerRecordO2[i + 1].scalersInps[classindex];
+      diff = s1 - s0;
+    } else {
+      LOG(error) << "Wrong type:" << type;
+      return scals; // wrong type
+    }
+    scals.emplace_back(std::pair<double_t, double_t>{diff, timeDiff});
+  }
+  return scals;
+};
 // returns the pair of global (levelled) interaction rate, as well as instantaneous interpolated
 // rate in Hz at a certain orbit number within the run
 // type - 7 : inputs

Detectors/CTP/macro/GetRates.C

@@ -12,17 +12,44 @@
 #if !defined(__CLING__) || defined(__ROOTCLING__)
 #include <CCDB/BasicCCDBManager.h>
 #include <DataFormatsCTP/Configuration.h>
-#include <DataFormatsCTP/CTPRateFetcher.h>
+#include "CTPWorkflowScalers/ctpCCDBManager.h"
+#include "Framework/Logger.h"
 #endif
 using namespace o2::ctp;
 
-void TestFetcher(int runNumber = 535087)
+void GetRates(int run = 559617)
 {
-  auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
-  std::pair<int64_t, int64_t> pp = ccdb.getRunDuration(runNumber);
-  long ts = pp.first + 60;
-  std::cout << "Run duration:" << pp.first << " " << pp.second << std::endl;
-  // Opening run
-  CTPRateFetcher fetcher;
-  fetcher.setupRun(runNumber, &ccdb, ts, 1);
+  uint64_t inputmaskCum = 0, classmackCum = 0;
+  int ntrigSel = 0;
+
+  auto& cmb = o2::ccdb::BasicCCDBManager::instance();
+  auto ctpcfg = cmb.getSpecificForRun<o2::ctp::CTPConfiguration>("CTP/Config/Config", run);
+  if (!ctpcfg) {
+    LOGP(error, "Can not get config for run {}", run);
+    return;
+  }
+  CTPConfiguration ctpconfig;
+  ctpconfig.loadConfigurationRun3(ctpcfg->getConfigString());
+  ctpconfig.printStream(std::cout);
+  auto& triggerclasses = ctpconfig.getCTPClasses();
+  LOGP(info, "Found {} trigger classes", triggerclasses.size());
+  int indexInList = 0;
+  for (const auto& trgclass : triggerclasses) {
+    uint64_t inputmask = 0;
+    if (trgclass.descriptor != nullptr) {
+      inputmask = trgclass.descriptor->getInputsMask();
+      // LOGP(info, "inputmask: {:#x}", inputmask);
+    }
+    trgclass.printStream(std::cout);
+    //    std::cout << indexInList << ": " << trgclass.name << ", input mask 0x" << std::hex << inputmask << ", class mask 0x" << trgclass.classMask << std::dec << std::endl;
+    indexInList++;
+    if (trgclass.cluster->getClusterDetNames().find("TRD") != std::string::npos || trgclass.cluster->getClusterDetNames().find("trd") != std::string::npos) {
+      LOGP(info, "Found TRD trigger cluster, class mask: {:#x}, input mask: {:#x}", trgclass.classMask, inputmask);
+      inputmaskCum |= inputmask;
+      classmackCum |= trgclass.classMask;
+      ntrigSel++;
+    }
+  }
+
+  LOGP(info, "Found {} triggers with TRD: classMasks: {:#x}  inputMasks: {:#x}", ntrigSel, classmackCum, inputmaskCum);
 }

Detectors/CTP/macro/PlotPbLumi.C

@@ -20,6 +20,11 @@
 #include "CCDB/BasicCCDBManager.h"
 #include "DataFormatsCTP/Scalers.h"
 #include "DataFormatsCTP/Configuration.h"
+#include "DataFormatsParameters/GRPLHCIFData.h"
+#include "TGraph.h"
+#include "TMath.h"
+#include "TCanvas.h"
+#include "TStyle.h"
 #include <string>
 #include <map>
 #include <iostream>

Detectors/CTP/macro/TestGetRates.C

@@ -16,17 +16,33 @@
 #endif
 using namespace o2::ctp;
 
-void TestGetRates(int runNumber = 557251)
+void TestGetRates(int runN = 0)
 {
+  std::vector<int> runs;
+  std::vector<std::string> codes = {"T0VTX", "T0VTX", "ZNChadronic", "ZNChadronic", "T0VTX"};
+  if (runN == 0) {
+    runs = {529066, 539218, 544013, 544518, 557251};
+  } else {
+    runs.push_back(runN);
+  }
   auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
-  // Opening run
-  std::pair<int64_t, int64_t> pp = ccdb.getRunDuration(runNumber);
-  long ts = pp.first + 60;
-  std::cout << "Run duration:" << pp.first << " " << pp.second << std::endl;
-  CTPRateFetcher fetcher;
-  fetcher.setupRun(runNumber, &ccdb, ts, 1);
-  fetcher.setOrbit(1);
-  std::array<double, 3> rates;
-  fetcher.getRates(rates, &ccdb, runNumber, "T0VTX");
-  std::cout << "Start:" << rates[0] << " End:" << rates[1] << " Middle:" << rates[2] << std::endl;
+  int i = 0;
+  for (auto const& runNumber : runs) {
+    // Opening run
+    std::pair<int64_t, int64_t> pp = ccdb.getRunDuration(runNumber);
+    long ts = pp.first + 60;
+    // std::cout << "Run duration:" << pp.first << " " << pp.second << std::endl;
+    std::cout << "===> RUN:" << runNumber << " duration:" << (pp.second - pp.first) / 1000. << std::endl;
+
+    CTPRateFetcher fetcher;
+    fetcher.setupRun(runNumber, &ccdb, ts, 1);
+    fetcher.setOrbit(1);
+    std::array<double, 3> rates;
+    fetcher.getRates(rates, &ccdb, runNumber, codes[i]);
+    std::cout << "Start:" << rates[0] << " End:" << rates[1] << " Middle:" << rates[2] << " code:" << codes[i] << std::endl;
+    double lumi1 = fetcher.getLumi(&ccdb, runNumber, codes[i], 0);
+    double lumi2 = fetcher.getLumi(&ccdb, runNumber, codes[i], 1);
+    std::cout << " Lumi NO pile up corr:" << lumi1 << " Lumi with pile upcorr:" << lumi2 << " code:" << codes[i] << std::endl;
+    i++;
+  }
 }