dev: improvements of pb lumi macro (#14822)

* dev: improvements of pb lumi macro

* forgot clang

* more options

* clang

* clang

* multigraph (thanks to Andreas Molander)

---------

Co-authored-by: Roman Lietava <romanlietava@romans-macbook-pro.home>

Author: lietava

Detectors/CTP/macro/PlotPbLumi.C

@@ -9,10 +9,10 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-/// \file TestCTPScalers.C
+/// \file PlotPbLumi.C
 /// \brief create CTP scalers, test it and add to database
 /// \author Roman Lietava
-// root -b -q "GetScalers.C(\"519499\", 1656286373953)"
+// root "PLotPbLumi.C(519499)"
 #if !defined(__CLING__) || defined(__ROOTCLING__)
 
 #include <fairlogger/Logger.h>
@@ -30,33 +30,48 @@
 #include <iostream>
 #endif
 using namespace o2::ctp;
-void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-test.cern.ch:8080")
+//
+// sum = 0: TCE and TSC separatelly otherwise TCE and (TCE+TSC)
+// qc = 0: takes scalers from CCDB (available only for finished runs) otherwise from QCCDB (available for active runs)
+// t0-tlast: window in seconds counted from beginning of run
+//
+void PlotPbLumi(int runNumber = 567905, bool sum = 0, bool qc = 0, Double_t t0 = 0., Double_t tlast = 0.)
 { //
-  // what = 1: znc rate
-  // what = 2: (TCE+TSC)/ZNC
-  // what = 3: TCE/ZNC
-  std::string mCCDBPathCTPScalers = "CTP/Calib/Scalers";
+  // PLots in one canvas
+  // znc rate/28
+  // R = (TCE+TSC)*TVX*B*/ZNC*28
+  // R = TCE*TVX*B/ZNC*28
+  // R = VCH*TVX*B/ZNC*28
+  std::string ccdbHost = "http://alice-ccdb.cern.ch";
+  std::string mCCDBPathCTPScalers = "/CTP/Calib/Scalers";
+  std::string mCCDBPathCTPScalersQC = "qc/CTP/Scalers";
   std::string mCCDBPathCTPConfig = "CTP/Config/Config";
   auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance();
   // 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 lhcifdata = ccdbMgr.getForRun<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", runNumber);
+  // auto lhcifdata = ccdbMgr.getSpecific<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", timeStamp, metadata);
+  if (!lhcifdata) {
+    throw std::runtime_error("No GRPLHCIFData for run " + std::to_string(runNumber));
+  }
   auto bfilling = lhcifdata->getBunchFilling();
   std::vector<int> bcs = bfilling.getFilledBCs();
   int nbc = bcs.size();
   std::cout << "Number of interacting bc:" << nbc << std::endl;
   // Scalers
   std::string srun = std::to_string(runNumber);
-  metadata.clear(); // can be empty
+  std::map<string, string> metadata;
   metadata["runNumber"] = srun;
-  ccdbMgr.setURL("http://ccdb-test.cern.ch:8080");
-  auto scl = ccdbMgr.getSpecific<CTPRunScalers>(mCCDBPathCTPScalers, timeStamp, metadata);
+  CTPRunScalers* scl = nullptr;
+  if (qc) {
+    ccdbMgr.setURL("http://ali-qcdb-gpn.cern.ch:8083");
+    scl = ccdbMgr.getSpecific<CTPRunScalers>(mCCDBPathCTPScalersQC, timeStamp, metadata);
+  } else {
+    scl = ccdbMgr.getSpecific<CTPRunScalers>(mCCDBPathCTPScalers, timeStamp, metadata);
+  }
   if (scl == nullptr) {
     LOG(info) << "CTPRunScalers not in database, timestamp:" << timeStamp;
     return;
@@ -65,6 +80,7 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   std::vector<CTPScalerRecordO2> recs = scl->getScalerRecordO2();
   //
   // CTPConfiguration ctpcfg;
+  ccdbMgr.setURL("http://alice-ccdb.cern.ch");
   auto ctpcfg = ccdbMgr.getSpecific<CTPConfiguration>(mCCDBPathCTPConfig, timeStamp, metadata);
   if (ctpcfg == nullptr) {
     LOG(info) << "CTPRunConfig not in database, timestamp:" << timeStamp;
@@ -85,6 +101,7 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   int tsc = 255;
   int tce = 255;
   int vch = 255;
+  int zncclsi = 255;
   for (auto const& cls : ctpcls) {
     if (cls.name.find("CMTVXTSC-B-NOPF") != std::string::npos && tsc == 255) {
       int itsc = cls.getIndex();
@@ -104,6 +121,12 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
       // vch = scl->getScalerIndexForClass(ivch);
       std::cout << cls.name << ":" << vch << ":" << ivch << std::endl;
     }
+    if (cls.name.find("C1ZNC-B-NOPF-CRU") != std::string::npos) {
+      int iznc = cls.getIndex();
+      zncclsi = clsIndexToScaler[iznc];
+      // vch = scl->getScalerIndexForClass(ivch);
+      std::cout << cls.name << ":" << zncclsi << ":" << iznc << std::endl;
+    }
   }
   if (tsc == 255 || tce == 255 || vch == 255) {
     std::cout << " One of dcalers not available, check config to find alternative)" << std::endl;
@@ -120,11 +143,39 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   double_t orbit0 = recs[0].intRecord.orbit;
   int n = recs.size() - 1;
   std::cout << " Run duration:" << Trun << " Scalers size:" << n + 1 << std::endl;
-  Double_t x[n], znc[n], zncpp[n];
-  Double_t tcetsctoznc[n], tcetoznc[n], vchtoznc[n];
-  for (int i = 0; i < n; i++) {
-    x[i] = (double_t)(recs[i + 1].intRecord.orbit + recs[i].intRecord.orbit) / 2. - orbit0;
-    x[i] *= 88e-6;
+  //
+  int i0 = 0;
+  int ilast = 0;
+  if (t0 != 0. || tlast != 0.) {
+    for (int i = 0; i < n; i++) {
+      double_t ttime = recs[i].epochTime - time0;
+      if (!i0 && t0 < ttime) {
+        i0 = i;
+      }
+      if (!ilast && tlast < ttime) {
+        ilast = i;
+      }
+    }
+  } else {
+    ilast = n;
+  }
+  n = ilast - i0;
+  std::cout << "i0:" << i0 << " ilast:" << ilast << std::endl;
+  // Double_t x[n], znc[n], zncpp[n];
+  std::vector<Double_t> xvec(n), zncvec(n), zncppvec(n), zncclassvec(n);
+  Double_t* x = xvec.data();
+  Double_t* znc = zncvec.data();
+  Double_t* zncpp = zncppvec.data();
+  Double_t* zncclass = zncclassvec.data();
+  // Double_t tcetsctoznc[n], tcetoznc[n], vchtoznc[n];
+  std::vector<Double_t> tcetsctozncvec(n), tcetozncvec(n), vchtozncvec(n);
+  Double_t* tcetsctoznc = tcetsctozncvec.data();
+  Double_t* tcetoznc = tcetozncvec.data();
+  Double_t* vchtoznc = vchtozncvec.data();
+  for (int i = i0; i < ilast; i++) {
+    int iv = i - i0;
+    x[iv] = (double_t)(recs[i + 1].intRecord.orbit + recs[i].intRecord.orbit) / 2. - orbit0;
+    x[iv] *= 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;
@@ -133,33 +184,53 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
     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;
-    zncpp[i] = zncipp / 28.;
-    znc[i] = znci / 28. / tt;
+    zncpp[iv] = zncipp / 28.;
+    znc[iv] = znci / 28. / tt;
+    // znc class
+    znci = recs[i + 1].scalers[zncclsi].l1Before - recs[i].scalers[zncclsi].l1Before;
+    zncclass[iv] = znci / 28. / tt;
+    // std::cout << znc[i]/zncclass[i] << std::endl;
     //
-    auto had = recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
+    double_t had = 0;
+    if (sum) {
+      had += recs[i + 1].scalers[tce].lmBefore - recs[i].scalers[tce].lmBefore;
+    }
+    double_t mutce = -TMath::Log(1. - had / tt / nbc / frev);
     // 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;
+    tcetsctoznc[iv] = (double_t)(had) / zncpp[iv] / 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;
+    tcetoznc[iv] = (double_t)(had) / zncpp[iv] / tt;
     had = recs[i + 1].scalers[vch].lmBefore - recs[i].scalers[vch].lmBefore;
+    double_t muvch = -TMath::Log(1. - had / tt / nbc / frev);
+
     // rat = (double_t)(had)/double_t(recs[i+1].scalersInps[25] - recs[i].scalersInps[25])*28;
-    vchtoznc[i] = (double_t)(had) / zncpp[i] / tt;
+    vchtoznc[iv] = (double_t)(had) / zncpp[iv] / tt;
+    // std::cout << "muzdc:" << mu << " mu tce:" << mutce << " muvch:" << muvch << std::endl;
   }
   //
   gStyle->SetMarkerSize(0.5);
   TGraph* gr1 = new TGraph(n, x, znc);
+  TGraph* gr11 = new TGraph(n, x, zncpp);    // PileuP corrected
+  TGraph* gr12 = new TGraph(n, x, zncclass); // NOT PileuP corrected
   TGraph* gr2 = new TGraph(n, x, tcetsctoznc);
   TGraph* gr3 = new TGraph(n, x, tcetoznc);
   TGraph* gr4 = new TGraph(n, x, vchtoznc);
   gr1->SetMarkerStyle(20);
+  gr11->SetMarkerStyle(20);
+  gr12->SetMarkerStyle(20);
+  gr11->SetMarkerColor(kRed);
+  gr12->SetMarkerColor(kBlue);
   gr2->SetMarkerStyle(21);
   gr3->SetMarkerStyle(23);
   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");
+  if (sum) {
+    gr2->SetTitle("R=(TSC+TCE)*TVTX*B*28/ZNC; time[sec]; R");
+  } else {
+    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");
@@ -168,10 +239,17 @@ void PlotPbLumi(int runNumber, int fillN, std::string ccdbHost = "http://ccdb-te
   gr4->SetTitle("R=(VCH)*TVTX*B*28/ZNC; time[sec]; R");
   // gr4->GetHistogram()->SetMaximum(0.6);
   // gr4->GetHistogram()->SetMinimum(0.4);
+  TMultiGraph* mg1 = new TMultiGraph();
+  mg1->SetTitle("R=ZNC/28 rate [Hz] (red=PilUp Corrected); time[sec]; R");
+  mg1->Add(gr1);
+  mg1->Add(gr11);
+  mg1->Add(gr12);
   TCanvas* c1 = new TCanvas("c1", srun.c_str(), 200, 10, 800, 500);
+  std::string title = "RUN " + std::to_string(runNumber);
+  c1->SetTitle(title.c_str());
   c1->Divide(2, 2);
   c1->cd(1);
-  gr1->Draw("AP");
+  mg1->Draw("AP");
   c1->cd(2);
   gr2->Draw("AP");
   c1->cd(3);