add TOF diagnostic from ccdb in digitisation (empty crates, TRM errors)

DataFormats/Detectors/TOF/include/DataFormatsTOF/Diagnostic.h

@@ -33,10 +33,10 @@ class Diagnostic
   Diagnostic() = default;
   int fill(ULong64_t pattern);
   int fill(ULong64_t pattern, int frequency);
-  int getFrequency(ULong64_t pattern);                                                    // Get frequency
-  int getFrequencyROW() { return getFrequency(0); }                                       // Readout window frequency
-  int getFrequencyEmptyCrate(int crate) { return getFrequency(getEmptyCrateKey(crate)); } // empty crate frequency
-  int getFrequencyEmptyTOF() { return getFrequency(1); }                                  // empty crate frequency
+  int getFrequency(ULong64_t pattern) const;                                                    // Get frequency
+  int getFrequencyROW() const { return getFrequency(0); }                                       // Readout window frequency
+  int getFrequencyEmptyCrate(int crate) const { return getFrequency(getEmptyCrateKey(crate)); } // empty crate frequency
+  int getFrequencyEmptyTOF() const { return getFrequency(1); }                                  // empty crate frequency
   int fillNoisy(int channel, int frequency = 1) { return fill(getNoisyChannelKey(channel), frequency); }
   int fillROW() { return fill(0); }
   int fillEmptyCrate(int crate, int frequency = 1) { return fill(getEmptyCrateKey(crate), frequency); }
@@ -49,7 +49,8 @@ class Diagnostic
   void fill(const Diagnostic& diag);                       // for calibration
   void fill(const gsl::span<const o2::tof::Diagnostic>){}; // for calibration
   void merge(const Diagnostic* prev);
-  void getNoisyMap(Bool_t* output); // set true in output channel array
+  void getNoisyMap(Bool_t* output, int noisyThr = 1) const; // set true in output channel array
+  void getNoisyLevelMap(Char_t* output) const;              // set true in output channel array
   unsigned long size() const { return mVector.size(); }
   ULong64_t getPattern(int i) const
   {
@@ -64,6 +65,8 @@ class Diagnostic
   int getChannel(ULong64_t pattern) const;
   int getNoisyLevel(ULong64_t pattern) const;
 
+  const std::map<ULong64_t, uint32_t>& getVector() const { return mVector; }
+
  private:
   std::map<ULong64_t, uint32_t> mVector; // diagnostic frequency vector (key/pattern , frequency)
 

DataFormats/Detectors/TOF/src/Diagnostic.cxx

@@ -49,7 +49,7 @@ int Diagnostic::fill(ULong64_t pattern, int frequency)
   return frequency;
 }
 
-int Diagnostic::getFrequency(ULong64_t pattern)
+int Diagnostic::getFrequency(ULong64_t pattern) const
 {
   auto pairC = mVector.find(pattern);
   if (pairC != mVector.end()) {
@@ -61,7 +61,7 @@ int Diagnostic::getFrequency(ULong64_t pattern)
 
 void Diagnostic::print() const
 {
-  LOG(INFO) << "Diagnostic patterns";
+  LOG(info) << "Diagnostic patterns";
   for (const auto& [key, value] : mVector) {
     std::cout << key << " = " << value << "; ";
   }
@@ -120,22 +120,22 @@ int Diagnostic::getNoisyLevel(ULong64_t pattern) const
 
 void Diagnostic::fill(const Diagnostic& diag)
 {
-  LOG(DEBUG) << "Filling diagnostic word";
+  LOG(debug) << "Filling diagnostic word";
   for (auto const& el : diag.mVector) {
-    LOG(DEBUG) << "Filling diagnostic pattern " << el.first << " adding " << el.second << " to " << getFrequency(el.first) << " --> " << el.second + getFrequency(el.first);
+    LOG(debug) << "Filling diagnostic pattern " << el.first << " adding " << el.second << " to " << getFrequency(el.first) << " --> " << el.second + getFrequency(el.first);
     fill(el.first, el.second);
   }
 }
 
 void Diagnostic::merge(const Diagnostic* prev)
 {
-  LOG(DEBUG) << "Merging diagnostic words";
+  LOG(debug) << "Merging diagnostic words";
   for (auto const& el : prev->mVector) {
     fill(el.first, el.second + getFrequency(el.first));
   }
 }
 
-void Diagnostic::getNoisyMap(Bool_t* output)
+void Diagnostic::getNoisyLevelMap(Char_t* output) const
 {
   // set true in output channel array
   for (auto pair : mVector) {
@@ -146,6 +146,23 @@ void Diagnostic::getNoisyMap(Bool_t* output)
       continue;
     }
 
-    output[getChannel(key)] = true;
+    output[getChannel(key)] = getNoisyLevel(key);
+  }
+}
+
+void Diagnostic::getNoisyMap(Bool_t* output, int noisyThr) const
+{
+  // set true in output channel array
+  for (auto pair : mVector) {
+    auto key = pair.first;
+    int slot = getSlot(key);
+
+    if (slot != 14) {
+      continue;
+    }
+
+    if (getNoisyLevel(key) >= noisyThr) {
+      output[getChannel(key)] = true;
+    }
   }
 }

Detectors/TOF/base/include/TOFBase/CalibTOFapi.h

@@ -20,6 +20,8 @@
 #include "CCDB/BasicCCDBManager.h"
 #include "DataFormatsTOF/CalibLHCphaseTOF.h"
 #include "DataFormatsTOF/CalibTimeSlewingParamTOF.h"
+#include "TOFBase/Geo.h"
+#include "DataFormatsTOF/Diagnostic.h"
 
 namespace o2
 {
@@ -35,9 +37,10 @@ class CalibTOFapi
   using CcdbApi = o2::ccdb::CcdbApi;
 
  public:
-  CalibTOFapi() = default;
+  void resetDia();
+  CalibTOFapi();
   CalibTOFapi(const std::string url);
-  CalibTOFapi(long timestamp, o2::dataformats::CalibLHCphaseTOF* phase, o2::dataformats::CalibTimeSlewingParamTOF* slew) : mTimeStamp(timestamp), mLHCphase(phase), mSlewParam(slew) {}
+  CalibTOFapi(long timestamp, o2::dataformats::CalibLHCphaseTOF* phase, o2::dataformats::CalibTimeSlewingParamTOF* slew) : mTimeStamp(timestamp), mLHCphase(phase), mSlewParam(slew) { CalibTOFapi(); }
   ~CalibTOFapi() = default;
   void setTimeStamp(long t)
   {
@@ -50,6 +53,7 @@ class CalibTOFapi
   }
   void readLHCphase();
   void readTimeSlewingParam();
+  void readDiagnosticFrequencies();
   void writeLHCphase(LhcPhase* phase, std::map<std::string, std::string> metadataLHCphase, uint64_t minTimeSTamp, uint64_t maxTimeStamp);
   void writeTimeSlewingParam(SlewParam* param, std::map<std::string, std::string> metadataChannelCalib, uint64_t minTimeSTamp, uint64_t maxTimeStamp = 0);
   float getTimeCalibration(int ich, float tot);
@@ -61,10 +65,27 @@ class CalibTOFapi
   SlewParam& getSlewParamObj() { return *mSlewParam; }
   LhcPhase* getLhcPhase() { return mLHCphase; }
 
+  int getNoisyThreshold() const { return mNoisyThreshold; }
+  void setNoisyThreshold(int val) { mNoisyThreshold = val; }
+  float getEmptyTOFProb() const { return mEmptyTOF; }
+  const float* getEmptyCratesProb() const { return mEmptyCrateProb; }
+  const std::vector<std::pair<int, float>>& getNoisyProb() const { return mNoisy; }
+  const std::vector<std::pair<int, float>>& getTRMerrorProb() const { return mTRMerrorProb; }
+  const std::vector<int>& getTRMmask() const { return mTRMmask; }
+
  private:
   long mTimeStamp;                 ///< timeStamp for queries
   LhcPhase* mLHCphase = nullptr;   ///< object for LHC phase
   SlewParam* mSlewParam = nullptr; ///< object for timeslewing (containing info also for offset and problematic)
+  Diagnostic* mDiaFreq = nullptr;  ///< object for Diagnostic Frequency
+
+  // info from diagnostic
+  int mNoisyThreshold = 1;                          ///< threshold to be noisy
+  float mEmptyTOF = 0;                              ///< probability to have TOF fully empty
+  float mEmptyCrateProb[Geo::kNCrate];              ///< probability to have an empty crate in the current readout window
+  std::vector<std::pair<int, float>> mNoisy;        ///< probTRMerror
+  std::vector<std::pair<int, float>> mTRMerrorProb; ///< probTRMerror
+  std::vector<int> mTRMmask;                        ///< mask error for TRM
 
   ClassDefNV(CalibTOFapi, 1);
 };

Detectors/TOF/base/include/TOFBase/Strip.h

@@ -85,6 +85,7 @@ class Strip
   void fillOutputContainer(std::vector<o2::tof::Digit>& digits);
 
   static int mDigitMerged;
+  std::map<ULong64_t, o2::tof::Digit>& getDigitMap() { return mDigits; }
 
  protected:
   Int_t mStripIndex = -1;                      ///< Strip ID

Detectors/TOF/base/src/CalibTOFapi.cxx

@@ -16,9 +16,24 @@ using namespace o2::tof;
 
 ClassImp(o2::tof::CalibTOFapi);
 
-CalibTOFapi::CalibTOFapi(const std::string url)
+void CalibTOFapi::resetDia()
+{
+  memset(mEmptyCrateProb, 0., Geo::kNCrate * 4);
+  mTRMerrorProb.clear();
+  mTRMmask.clear();
+  mNoisy.clear();
+}
+
+CalibTOFapi::CalibTOFapi()
 {
+  resetDia();
+}
+
+//______________________________________________________________________
 
+CalibTOFapi::CalibTOFapi(const std::string url)
+{
+  CalibTOFapi();
   // setting the URL to the CCDB manager
 
   setURL(url);
@@ -49,6 +64,70 @@ void CalibTOFapi::readTimeSlewingParam()
 
 //______________________________________________________________________
 
+void CalibTOFapi::readDiagnosticFrequencies()
+{
+  static const int NCH_PER_CRATE = Geo::NSTRIPXSECTOR * Geo::NPADS;
+  // getting the Diagnostic Frequency calibration
+  // needed for simulation
+
+  auto& mgr = CcdbManager::instance();
+  mDiaFreq = mgr.getForTimeStamp<Diagnostic>("TOF/Calib/Diagnostic", mTimeStamp);
+
+  resetDia();
+
+  if (!mDiaFreq->getFrequencyROW()) {
+    return;
+  }
+
+  float nrow = (float)mDiaFreq->getFrequencyROW();
+  mEmptyTOF = mDiaFreq->getFrequencyEmptyTOF() / nrow;
+
+  nrow -= mDiaFreq->getFrequencyEmptyTOF();
+
+  if (nrow < 1) {
+    return;
+  }
+
+  // fill empty crates
+  int ncrate[Geo::kNCrate];
+  for (int i = 0; i < Geo::kNCrate; i++) {
+    ncrate[i] = mDiaFreq->getFrequencyEmptyCrate(i) - mDiaFreq->getFrequencyEmptyTOF(); // counts of empty crate for non-empty event
+    mEmptyCrateProb[i] = ncrate[i] / nrow;
+  }
+
+  const auto vectorDia = mDiaFreq->getVector();
+  // fill TRM errors and noisy
+  for (auto pair : vectorDia) {
+    auto key = pair.first;
+    int slot = mDiaFreq->getSlot(key);
+
+    if (slot < 13 && slot > 2) { // is TRM
+      int icrate = mDiaFreq->getCrate(key);
+      int crateslot = icrate * 100 + slot;
+      mTRMerrorProb.push_back(std::make_pair(crateslot, pair.second / (nrow - ncrate[icrate])));
+      mTRMmask.push_back(key - mDiaFreq->getTRMKey(icrate, slot)); // remove crate and slot from the key (28 bit errors remaining)
+      continue;
+    }
+
+    int channel = mDiaFreq->getChannel(key);
+    if (channel > -1) { // noisy
+      int crate = channel / NCH_PER_CRATE;
+      mNoisy.push_back(std::make_pair(channel, pair.second / (nrow - ncrate[crate])));
+      continue;
+    }
+  }
+
+  std::sort(mTRMerrorProb.begin(), mTRMerrorProb.end(), [](const auto& a, const auto& b) {
+    return a.first < b.first;
+  });
+
+  std::sort(mNoisy.begin(), mNoisy.end(), [](const auto& a, const auto& b) {
+    return a.first < b.first;
+  });
+}
+
+//______________________________________________________________________
+
 void CalibTOFapi::writeLHCphase(LhcPhase* phase, std::map<std::string, std::string> metadataLHCphase, uint64_t minTimeStamp, uint64_t maxTimeStamp)
 {
 

Detectors/TOF/base/src/WindowFiller.cxx

@@ -456,27 +456,27 @@ void WindowFiller::fillDiagnosticFrequency()
       int slot = 0;
       if (mReadoutWindowData[j].isEmptyCrate(ic)) {
         mDiagnosticFrequency.fillEmptyCrate(ic);
-      } else {
-        isTOFempty = false;
-      }
-      if (dia) {
-        int fd = mReadoutWindowData[j].firstDia();
-        int lastdia = fd + dia;
-
-        ULong64_t key;
-        for (int dd = fd; dd < lastdia; dd++) {
-          if (mPatterns[dd] >= 28) {
-            slot = mPatterns[dd] - 28;
-            key = mDiagnosticFrequency.getTRMKey(ic, slot);
-            continue;
-          }
-
-          key += (1 << mPatterns[dd]);
-
-          if (dd + 1 == lastdia || mPatterns[dd + 1] >= 28) {
-            mDiagnosticFrequency.fill(key);
+        if (dia) {
+          int fd = mReadoutWindowData[j].firstDia();
+          int lastdia = fd + dia;
+
+          ULong64_t key;
+          for (int dd = fd; dd < lastdia; dd++) {
+            if (mPatterns[dd] >= 28) {
+              slot = mPatterns[dd] - 28;
+              key = mDiagnosticFrequency.getTRMKey(ic, slot);
+              continue;
+            }
+
+            key += (1 << mPatterns[dd]);
+
+            if (dd + 1 == lastdia || mPatterns[dd + 1] >= 28) {
+              mDiagnosticFrequency.fill(key);
+            }
           }
         }
+      } else {
+        isTOFempty = false;
       }
     }
   }