dev: decoder: checking only trigger class bits which belongs to run

Author: lietava

Detectors/CTP/macro/PlotPbLumi.C

@@ -20,37 +20,56 @@
 #include "CCDB/BasicCCDBManager.h"
 #include "DataFormatsCTP/Scalers.h"
 #include "DataFormatsCTP/Configuration.h"
+#include <DataFormatsParameters/GRPLHCIFData.h>
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TStyle.h"
+#include "TF1.h"
 #include <string>
 #include <map>
 #include <iostream>
 #endif
 using namespace o2::ctp;
-void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-test.cern.ch:8080")
-{ //
-  // what = 1: znc rate
-  // what = 2: (TCE+TSC)/ZNC
-  // what = 3: TCE/ZNC
+using namespace std::chrono;
+//
+// if fillN = 0: pileup correction not done
+// QCDB =1 : use for ongoing run
+//
+void PlotPbLumi(int runNumber, int fillN = 0, int QCDB = 0,  std::string ccdbHost = "http://alice-ccdb.cern.ch" )
+{ // "http://ccdb-test.cern.ch:8080"
   std::string mCCDBPathCTPScalers = "CTP/Calib/Scalers";
+  std::string mQCDBPathCTPScalers = "qc/CTP/Scalers";
   std::string mCCDBPathCTPConfig = "CTP/Config/Config";
   auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance();
+  ccdbMgr.setURL(ccdbHost);
   // Timestamp
   auto soreor = ccdbMgr.getRunDuration(runNumber);
   uint64_t timeStamp = (soreor.second - soreor.first) / 2 + soreor.first;
   std::cout << "Timestamp:" << timeStamp << std::endl;
   // Filling
-  std::string sfill = std::to_string(fillN);
   std::map<string, string> metadata;
-  metadata["fillNumber"] = sfill;
-  auto lhcifdata = ccdbMgr.getSpecific<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", timeStamp, metadata);
-  auto bfilling = lhcifdata->getBunchFilling();
-  std::vector<int> bcs = bfilling.getFilledBCs();
-  int nbc = bcs.size();
-  std::cout << "Number of interacting bc:" << nbc << std::endl;
+  int nbc = 0;
+  if(fillN) {
+    std::string sfill = std::to_string(fillN);
+    metadata["fillNumber"] = sfill;
+    auto lhcifdata = ccdbMgr.getSpecific<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", timeStamp, metadata);
+    auto bfilling = lhcifdata->getBunchFilling();
+    std::vector<int> bcs = bfilling.getFilledBCs();
+    nbc = bcs.size();
+    std::cout << "Number of interacting bc:" << nbc << std::endl;
+  }
+   if(QCDB) {  // use this option for ongoing run
+    mCCDBPathCTPScalers = mQCDBPathCTPScalers;
+    ccdbMgr.setURL("http://ccdb-test.cern.ch:8080");
+    timeStamp = duration_cast< milliseconds >(system_clock::now().time_since_epoch()).count();
+    std::cout << "For scalers using Current time:" << timeStamp << std::endl;
+  }
   // Scalers
   std::string srun = std::to_string(runNumber);
   metadata.clear(); // can be empty
   metadata["runNumber"] = srun;
-  ccdbMgr.setURL("http://ccdb-test.cern.ch:8080");
   auto scl = ccdbMgr.getSpecific<CTPRunScalers>(mCCDBPathCTPScalers, timeStamp, metadata);
   if (scl == nullptr) {
     LOG(info) << "CTPRunScalers not in database, timestamp:" << timeStamp;
@@ -105,48 +124,72 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
     return;
   }
   //
-  // Anal
-  //
   // Times
   double_t frev = 11245;
   double_t time0 = recs[0].epochTime;
   double_t timeL = recs[recs.size() - 1].epochTime;
   double_t Trun = timeL - time0;
   double_t orbit0 = recs[0].intRecord.orbit;
   int n = recs.size() - 1;
+  if(runNumber == 559143) {
+    n = 400;
+  }
+  if(runNumber == 559561 ) {
+    n = n - 3;  // rate drops at the end
+  }
+  if(runNumber == 559575){
+    n = n - 6;
+  }
+  if(runNumber == 559617) {
+    n = n - 5;
+  }
+  if(runNumber == 559632) {
+   n = n - 6;
+  }
   std::cout << " Run duration:" << Trun << " Scalers size:" << n + 1 << std::endl;
-  Double_t x[n], znc[n], zncpp[n];
+  Double_t x[n], znc[n], zncpp[n], ex[n], eznc[n];
   Double_t tcetsctoznc[n], tcetoznc[n], vchtoznc[n];
+  Double_t etcetsctoznc[n], etcetoznc[n], evchtoznc[n];
+
   for (int i = 0; i < n; i++) {
+    ex[i] = 0;
     x[i] = (double_t)(recs[i + 1].intRecord.orbit + recs[i].intRecord.orbit) / 2. - orbit0;
     x[i] *= 88e-6;
     // x[i] = (double_t)(recs[i+1].epochTime + recs[i].epochTime)/2.;
     double_t tt = (double_t)(recs[i + 1].intRecord.orbit - recs[i].intRecord.orbit);
     tt = tt * 88e-6;
     //
-    // std::cout << recs[i+1].scalersInps[25] << std::endl;
     double_t znci = (double_t)(recs[i + 1].scalersInps[25] - recs[i].scalersInps[25]);
-    double_t mu = -TMath::Log(1. - znci / tt / nbc / frev);
-    double_t zncipp = mu * nbc * frev;
+    double_t zncipp = znci;
+    double_t mu = 0;
+    if(fillN) {
+       mu = -TMath::Log(1. - znci / tt / nbc / frev);
+      zncipp = mu * nbc * frev * tt;
+    }
     zncpp[i] = zncipp / 28.;
-    znc[i] = znci / 28. / tt;
-    //
-    auto had = recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
-    // std::cout << recs[i+1].scalers[tce].lmBefore << std::endl;
-    had += recs[i + 1].scalers[tsc].lmBefore - recs[i].scalers[tsc].lmBefore;
-    // rat = (double_t)(had)/double_t(recs[i+1].scalersInps[25] - recs[i].scalersInps[25])*28;
-    tcetsctoznc[i] = (double_t)(had) / zncpp[i] / tt;
-    had = recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
-    // rat = (double_t)(had)/double_t(recs[i+1].scalersInps[25] - recs[i].scalersInps[25])*28;
-    tcetoznc[i] = (double_t)(had) / zncpp[i] / tt;
-    had = recs[i + 1].scalers[vch].lmBefore - recs[i].scalers[vch].lmBefore;
-    // rat = (double_t)(had)/double_t(recs[i+1].scalersInps[25] - recs[i].scalersInps[25])*28;
-    vchtoznc[i] = (double_t)(had) / zncpp[i] / tt;
+    znc[i] = zncipp / 28. / tt;
+    eznc[i] = TMath::Sqrt(zncipp) / 28. / tt;
+    if(1){
+      //
+      auto had = recs[i + 1].scalers[tsc].lmBefore - recs[i].scalers[tsc].lmBefore;
+      had += recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
+      tcetsctoznc[i] = (double_t)(had) / zncpp[i];
+      etcetsctoznc[i] = TMath::Sqrt(tcetsctoznc[i]*(1-tcetsctoznc[i])/ zncpp[i]);
+      had = recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
+      double_t tcec = had;
+      tcetoznc[i] = (double_t)(had) / zncpp[i];
+      had = recs[i + 1].scalers[vch].lmBefore - recs[i].scalers[vch].lmBefore;
+      vchtoznc[i] = (double_t)(had) / zncpp[i] ;
+      //std::cout << "mu:" << mu << " zncpp corr:" <<  zncipp << "  zncraw:" << znci << " tce:" << tcec << " tce/had" << tcec/zncpp[i] << std::endl;  
+  } 
   }
   //
+  TFile myfile("file.root","RECREATE");
   gStyle->SetMarkerSize(0.5);
-  TGraph* gr1 = new TGraph(n, x, znc);
+  //TGraph* gr1 = new TGraph(n, x, znc);
+  TGraphErrors* gr1 = new TGraphErrors(n, x, znc, ex, eznc);
   TGraph* gr2 = new TGraph(n, x, tcetsctoznc);
+  //TGraphErrors* gr2 = new TGraphErrors(n, x, tcetsctoznc, ex, etcetsctoznc); // nom and denom are strongly correlated
   TGraph* gr3 = new TGraph(n, x, tcetoznc);
   TGraph* gr4 = new TGraph(n, x, vchtoznc);
   gr1->SetMarkerStyle(20);
@@ -155,6 +198,7 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   gr4->SetMarkerStyle(23);
   gr1->SetTitle("R=ZNC/28 rate [Hz]; time[sec]; R");
   gr2->SetTitle("R=(TSC+TCE)*TVTX*B*28/ZNC; time[sec]; R");
+  //gr2->SetTitle("R=TSC*TVTX*B*28/ZNC; time[sec]; R");
   // gr2->GetHistogram()->SetMaximum(1.1);
   // gr2->GetHistogram()->SetMinimum(0.9);
   gr3->SetTitle("R=(TCE)*TVTX*B*28/ZNC; time[sec]; R");
@@ -164,19 +208,29 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   // gr4->GetHistogram()->SetMaximum(0.6);
   // gr4->GetHistogram()->SetMinimum(0.4);
   TCanvas* c1 = new TCanvas("c1", srun.c_str(), 200, 10, 800, 500);
+  TF1 *fun = new TF1("poly0","[0]+x*[1]");
   c1->Divide(2, 2);
   c1->cd(1);
   gr1->Draw("AP");
   c1->cd(2);
+  fun->SetParameter(0,1);
+  gr2->Fit("poly0","FM");
   gr2->Draw("AP");
   c1->cd(3);
+  fun->SetParameter(0,0.5);
+  gr3->Fit("poly0");
   gr3->Draw("AP");
   c1->cd(4);
+  gr4->Fit("poly0");
   gr4->Draw("AP");
   // getRate test:
   double tt = timeStamp / 1000.;
   std::pair<double, double> r1 = scl->getRateGivenT(tt, 25, 7);
   std::cout << "ZDC input getRateGivetT:" << r1.first / 28. << " " << r1.second / 28. << std::endl;
   std::pair<double, double> r2 = scl->getRateGivenT(tt, tce, 1);
   std::cout << "LM before TCE class getRateGivetT:" << r2.first << " " << r2.second << std::endl;
+  gr1->Write();
+  gr2->Write();
+  gr3->Write();
+  gr4->Write();
 }

Detectors/CTP/macro/TestConfig.C

@@ -16,13 +16,13 @@
 #endif
 using namespace o2::ctp;
 
-void TestConfig(bool test = 0)
+void TestConfig(bool test = 1)
 {
   if (test == 0) {
     return;
   }
-  uint64_t timestamp = 1660196771632;
-  std::string run = "523148";
+  uint64_t timestamp = 1660276134898+10000;
+  std::string run = "523186";
   o2::ctp::ctpCCDBManager::setCCDBHost("https://alice-ccdb.cern.ch");
   bool ok;
   auto ctpcfg = o2::ctp::ctpCCDBManager::getConfigFromCCDB(timestamp, run, ok);
@@ -31,13 +31,17 @@ void TestConfig(bool test = 0)
   }
   CTPConfiguration ctpconfig;
   ctpconfig.loadConfigurationRun3(ctpcfg.getConfigString());
-  // ctpconfig.printStream(std::cout);
-  // return;
-  // ctpconfig.assignDescriptors();
-  // ctpconfig.createInputsInDecriptorsFromNames();
   ctpconfig.printStream(std::cout);
   auto& triggerclasses = ctpconfig.getCTPClasses();
   std::cout << "Found " << triggerclasses.size() << " trigger classes" << std::endl;
+  for (const auto& trg : triggerclasses) {
+      if (trg.cluster->maskCluster[o2::detectors::DetID::EMC]) {
+        // Class triggering EMCAL cluster
+        LOG(info) << "Found trigger class for EMCAL cluster: " << trg.name << " with input mask " << std::bitset<64>(trg.descriptor->getInputsMask());
+        trg.descriptor->getInputsMask();
+      }
+    }
+  return;
   int indexInList = 0;
   for (const auto& trgclass : triggerclasses) {
     uint64_t inputmask = 0;

Detectors/CTP/macro/TestFetcher.C

@@ -16,7 +16,7 @@
 #endif
 using namespace o2::ctp;
 
-void TestFetcher(int runNumber = 557251)
+void TestFetcher(int runNumber = 559456)
 {
   auto& ccdb = o2::ccdb::BasicCCDBManager::instance();
   std::pair<int64_t, int64_t> pp = ccdb.getRunDuration(runNumber);
@@ -25,7 +25,9 @@ void TestFetcher(int runNumber = 557251)
   // Opening run
   CTPRateFetcher fetcher;
   fetcher.setupRun(runNumber, &ccdb, ts, 0);
-  ccdb.setURL("http://ali-qcdb-gpn.cern.ch:8083/");
+  //ccdb.setURL("http://ali-qcdb-gpn.cern.ch:8083/");
+  ccdb.setURL("http://ccdb-test.cern.ch:8080/");
+
   std::string QCDBPathCTPScalers = "qc/CTP/Scalers";
   map<string, string> metadata; // can be empty
   std::string run = std::to_string(runNumber);

Detectors/CTP/reconstruction/include/CTPReconstruction/RawDataDecoder.h

@@ -22,6 +22,7 @@
 #include "Framework/InputRecord.h"
 #include "DataFormatsCTP/Digits.h"
 #include "DataFormatsCTP/LumiInfo.h"
+#include "DataFormatsCTP/Configuration.h"
 
 namespace o2
 {
@@ -43,14 +44,15 @@ class RawDataDecoder
   void setVerbose(bool v) { mVerbose = v; }
   void setMAXErrors(int m) { mErrorMax = m; }
   int setLumiInp(int lumiinp, std::string inp);
+  void setCTPConfig( CTPConfiguration cfg) { mCTPConfig = std::move(cfg); };
   uint32_t getIRRejected() const { return mIRRejected; }
   uint32_t getTCRRejected() const { return mTCRRejected; }
   std::vector<uint32_t>& getTFOrbits() { return mTFOrbits; }
   int getErrorIR() { return mErrorIR; }
   int getErrorTCR() { return mErrorTCR; }
   int init();
   static int shiftNew(const o2::InteractionRecord& irin, uint32_t TFOrbit, std::bitset<48>& inpmask, int64_t shift, int level, std::map<o2::InteractionRecord, CTPDigit>& digmap);
-  static int shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digitsMap, o2::pmr::vector<CTPDigit>& digits, uint32_t TFOrbit);
+  static int shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digitsMap, o2::pmr::vector<CTPDigit>& digits, uint32_t TFOrbit, uint64_t trgclassmask = 0xffffffffffffffff);
 
  private:
   static constexpr uint32_t TF_TRIGGERTYPE_MASK = 0x800;
@@ -79,6 +81,7 @@ class RawDataDecoder
   int mErrorTCR = 0;
   int mErrorMax = 3;
   bool mStickyError = false;
+  CTPConfiguration mCTPConfig;
 };
 } // namespace ctp
 } // namespace o2

Detectors/CTP/reconstruction/src/RawDataDecoder.cxx

@@ -293,7 +293,9 @@ int RawDataDecoder::decodeRaw(o2::framework::InputRecord& inputs, std::vector<o2
     // std::cout << "last lumi:" << nhb  << std::endl;
   }
   if (mDoDigits & mDecodeInps) {
-    shiftInputs(digitsMap, digits, mTFOrbit);
+    uint64_t trgclassmask = mCTPConfig.getTriggerClassMask();
+    std::cout << "trgclassmask:" << std::hex << trgclassmask << std::dec << std::endl;
+    shiftInputs(digitsMap, digits, mTFOrbit, trgclassmask);
   }
   if (mDoDigits && !mDecodeInps) {
     for (auto const& dig : digitsMap) {
@@ -519,7 +521,7 @@ int RawDataDecoder::shiftNew(const o2::InteractionRecord& irin, uint32_t TFOrbit
 }
 //
 
-int RawDataDecoder::shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digitsMap, o2::pmr::vector<CTPDigit>& digits, uint32_t TFOrbit)
+int RawDataDecoder::shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digitsMap, o2::pmr::vector<CTPDigit>& digits, uint32_t TFOrbit, uint64_t trgclassmask)
 {
   // int nClasswoInp = 0; // counting classes without input which should never happen
   int nLM = 0;
@@ -594,7 +596,7 @@ int RawDataDecoder::shiftInputs(std::map<o2::InteractionRecord, CTPDigit>& digit
     if ((d.CTPInputMask & L1MASKInputs).count()) {
       nL1++;
     }
-    if (d.CTPClassMask.count()) {
+    if ((d.CTPClassMask).to_ulong() & trgclassmask) {
       if (d.CTPInputMask.count()) {
         nTwI++;
       } else {

Detectors/CTP/workflow/include/CTPWorkflow/RawDecoderSpec.h

@@ -68,6 +68,7 @@ class RawDecoderSpec : public framework::Task
   uint32_t mNTFToIntegrate = 1;
   uint32_t mNHBIntegratedT = 0;
   uint32_t mNHBIntegratedV = 0;
+  bool mDecodeinputs = 0;
   std::deque<size_t> mHistoryT;
   std::deque<size_t> mHistoryV;
   RawDataDecoder mDecoder;

Detectors/CTP/workflow/src/RawDecoderSpec.cxx

@@ -20,13 +20,15 @@
 #include "CommonUtils/VerbosityConfig.h"
 #include "Framework/InputRecord.h"
 #include "DataFormatsCTP/TriggerOffsetsParam.h"
+#include "Framework/CCDBParamSpec.h"
+#include "DataFormatsCTP/Configuration.h"
 
 using namespace o2::ctp::reco_workflow;
 
 void RawDecoderSpec::init(framework::InitContext& ctx)
 {
-  bool decodeinps = ctx.options().get<bool>("ctpinputs-decoding");
-  mDecoder.setDecodeInps(decodeinps);
+  mDecodeinputs = ctx.options().get<bool>("ctpinputs-decoding");
+  mDecoder.setDecodeInps(mDecodeinputs);
   mNTFToIntegrate = ctx.options().get<int>("ntf-to-average");
   mVerbose = ctx.options().get<bool>("use-verbose-mode");
   int maxerrors = ctx.options().get<int>("print-errors-num");
@@ -42,7 +44,7 @@ void RawDecoderSpec::init(framework::InitContext& ctx)
   mOutputLumiInfo.inp2 = inp2;
   mMaxInputSize = ctx.options().get<int>("max-input-size");
   mMaxInputSizeFatal = ctx.options().get<bool>("max-input-size-fatal");
-  LOG(info) << "CTP reco init done. Inputs decoding here:" << decodeinps << " DoLumi:" << mDoLumi << " DoDigits:" << mDoDigits << " NTF:" << mNTFToIntegrate << " Lumi inputs:" << lumiinp1 << ":" << inp1 << " " << lumiinp2 << ":" << inp2 << " Max errors:" << maxerrors << " Max input size:" << mMaxInputSize << " MaxInputSizeFatal:" << mMaxInputSizeFatal;
+  LOG(info) << "CTP reco init done. Inputs decoding here:" << mDecodeinputs << " DoLumi:" << mDoLumi << " DoDigits:" << mDoDigits << " NTF:" << mNTFToIntegrate << " Lumi inputs:" << lumiinp1 << ":" << inp1 << " " << lumiinp2 << ":" << inp2 << " Max errors:" << maxerrors << " Max input size:" << mMaxInputSize << " MaxInputSizeFatal:" << mMaxInputSizeFatal;
   // mOutputLumiInfo.printInputs();
 }
 void RawDecoderSpec::endOfStream(framework::EndOfStreamContext& ec)
@@ -135,6 +137,11 @@ void RawDecoderSpec::run(framework::ProcessingContext& ctx)
   if (fatal_flag) {
     ret = mDecoder.decodeRawFatal(inputs, filter);
   } else {
+    if(mDecodeinputs) {
+     const auto ctpcfg = inputs.get<o2::ctp::CTPConfiguration*>("ctpconfig");
+     //ctpcfg->printStream(std::cout);
+     mDecoder.setCTPConfig(*ctpcfg);
+    }
     ret = mDecoder.decodeRaw(inputs, filter, mOutputDigits, lumiPointsHBF1);
   }
   if (ret == 1) {
@@ -176,6 +183,8 @@ void RawDecoderSpec::run(framework::ProcessingContext& ctx)
       mOutputLumiInfo.orbit = lumiPointsHBF1[0].orbit;
     }
     mOutputLumiInfo.counts = mCountsT;
+
+
     mOutputLumiInfo.countsFV0 = mCountsV;
     mOutputLumiInfo.nHBFCounted = mNHBIntegratedT;
     mOutputLumiInfo.nHBFCountedFV0 = mNHBIntegratedV;
@@ -199,6 +208,7 @@ o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool
 
   std::vector<o2::framework::OutputSpec> outputs;
   if (digits) {
+    inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("CTP/Config/Config", 1));
     outputs.emplace_back("CTP", "DIGITS", 0, o2::framework::Lifetime::Timeframe);
   }
   if (lumi) {