TRD update noisy MCM calibration (#8701)

* TRD fix noisy MCM identifyer

* TRD TrackBasedCalib uses DPL CCDB access

* Do not crash if noise map missing but raise alarm

Author: martenole

DataFormats/Detectors/TRD/include/DataFormatsTRD/NoiseCalibration.h

@@ -44,11 +44,11 @@ class NoiseStatusMCM
   static constexpr int getMcmIdxGlb(int hcid, int rob, int mcm) { return hcid * constants::NMCMHCMAX + (rob / 2) * constants::NMCMROB + mcm; }
 
   // setters
-  void setIsNoisy(int hcid, int rob, int mcm) { mNoiseFlag.set(hcid * constants::NMCMHCMAX + rob * constants::NMCMROB + mcm); }
+  void setIsNoisy(int hcid, int rob, int mcm) { mNoiseFlag.set(getMcmIdxGlb(hcid, rob, mcm)); }
   void setIsNoisy(int mcmIdxGlb) { mNoiseFlag.set(mcmIdxGlb); }
 
   // getters
-  bool getIsNoisy(int hcid, int rob, int mcm) const { return mNoiseFlag.test(hcid * constants::NMCMHCMAX + (rob / 2) * constants::NMCMROB + mcm); }
+  bool getIsNoisy(int hcid, int rob, int mcm) const { return mNoiseFlag.test(getMcmIdxGlb(hcid, rob, mcm)); }
   bool getIsNoisy(int mcmIdxGlb) const { return mNoiseFlag.test(mcmIdxGlb); }
   auto getNumberOfNoisyMCMs() const { return mNoiseFlag.count(); }
   bool isTrackletFromNoisyMCM(const Tracklet64& trklt) const { return getIsNoisy(trklt.getHCID(), trklt.getROB(), trklt.getMCM()); }

Detectors/TRD/calibration/include/TRDCalibration/TrackBasedCalib.h

@@ -54,6 +54,9 @@ class TrackBasedCalib
   /// Initialize the input arrays
   void setInput(const o2::globaltracking::RecoContainer& input);
 
+  /// Set the MCM noise map
+  void setNoiseMapMCM(const NoiseStatusMCM* map) { mNoiseCalib = map; };
+
   /// Reset the output
   void reset();
 
@@ -78,7 +81,7 @@ class TrackBasedCalib
   MatCorrType mMatCorr{MatCorrType::USEMatCorrNONE}; ///< if material correction should be done
   RecoParam mRecoParam;                              ///< parameters required for TRD reconstruction
   AngularResidHistos mAngResHistos;                  ///< aggregated data for the track based calibration
-  NoiseStatusMCM* mNoiseCalib{nullptr};              ///< CCDB object with information of noisy MCMs
+  const NoiseStatusMCM* mNoiseCalib{nullptr};        ///< CCDB object with information of noisy MCMs
   // input arrays which should not be modified since they are provided externally
   gsl::span<const TrackTRD> mTracksInITSTPCTRD;        ///< TRD tracks reconstructed from TPC or ITS-TPC seeds
   gsl::span<const TrackTRD> mTracksInTPCTRD;           ///< TRD tracks reconstructed from TPC or TPC seeds

Detectors/TRD/calibration/macros/checkNoisyMCMs.C

@@ -12,6 +12,21 @@
 /// \file checkNoisyMCMs.C
 /// \brief macro identify noisy MCMs based on the number of found tracklets
 
+/*
+  Provided an input file with tracklets and trigger records this macro produces a bitset
+  where noisy MCMs are flagged.
+  In case useMean is requested the mean number of tracklets sent from each MCM is calculated.
+  MCMs which don't send any tracklets are not taken into account. The 10% of MCMs which send
+  the highest number of tracklets are also not considered for the mean calculation. Then, MCMs
+  which sent more than 100 times as many tracklets as the mean are flagged as noisy. If useMean
+  is not used simply the MCMs are flagged which sent on average at least one tracklet every 10th
+  trigger.
+
+  Important: Compare the number of tracklets sent from the first and the last masked MCM. In case
+  there is a large discrepancy between the two numbers the noise threshold might need to be adjusted
+  or not enough statistics is used.
+*/
+
 #if !defined(__CLING__) || defined(__ROOTCLING__)
 #include "CCDB/CcdbApi.h"
 #include "DataFormatsTRD/NoiseCalibration.h"
@@ -26,10 +41,11 @@
 #include <utility>
 #include <map>
 #include <algorithm>
+#include <locale.h>
 
 #endif
 
-void checkNoisyMCMs(std::string inpFile = "trdtracklets.root", bool enableWriteCCDB = false)
+void checkNoisyMCMs(std::string inpFile = "trdtracklets.root", bool useMean = true, bool detailedOutput = true, bool enableWriteCCDB = false)
 {
   using namespace o2::trd;
 
@@ -51,13 +67,15 @@ void checkNoisyMCMs(std::string inpFile = "trdtracklets.root", bool enableWriteC
   tree->SetBranchAddress("TrackTrg", &trigInPtr);
   tree->SetBranchAddress("Tracklet", &trackletsInPtr);
 
-  std::array<std::pair<unsigned int, unsigned short>, constants::MAXHALFCHAMBER * constants::NMCMHCMAX> trackletCounter{};
+  std::array<std::pair<unsigned int, unsigned int>, constants::MAXHALFCHAMBER * constants::NMCMHCMAX> trackletCounter{};
   std::for_each(trackletCounter.begin(), trackletCounter.end(), [](auto& counter) { static int idx = 0; counter.second = idx++; });
 
   // count number of tracklets per MCM
   size_t totalTrackletCounter = 0;
+  size_t totalTriggerCounter = 0;
   for (int iEntry = 0; iEntry < tree->GetEntries(); ++iEntry) {
     tree->GetEntry(iEntry);
+    totalTriggerCounter += trigIn.size();
     for (const auto& tracklet : tracklets) {
       totalTrackletCounter++;
       int mcmGlb = NoiseStatusMCM::getMcmIdxGlb(tracklet.getHCID(), tracklet.getROB(), tracklet.getMCM());
@@ -67,29 +85,53 @@ void checkNoisyMCMs(std::string inpFile = "trdtracklets.root", bool enableWriteC
 
   // estimate noise threshold
   std::sort(trackletCounter.begin(), trackletCounter.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); // sort by number of tracklets per MCM
-  float mean = 0;
-  int nActiveMcms = 0;
-  for (int idx = 0; idx < static_cast<int>(0.9 * (constants::MAXHALFCHAMBER * constants::NMCMHCMAX)); ++idx) {
-    // get average number of tracklets discarding the 10% of MCMs with most entries
+  float mean90 = 0;                                                                                                          // mean excluding 10% of MCMs with most entries
+  int nActiveMcms90 = 0;
+  std::vector<int> nTrkltsPerActiveMcm; // for median calculation
+  for (int idx = 0; idx < constants::MAXHALFCHAMBER * constants::NMCMHCMAX; ++idx) {
     if (trackletCounter[idx].first > 0) {
-      nActiveMcms++;
-      mean += trackletCounter[idx].first;
+      nTrkltsPerActiveMcm.push_back(trackletCounter[idx].first);
+      if (idx < static_cast<int>(0.9 * (constants::MAXHALFCHAMBER * constants::NMCMHCMAX))) {
+        // for the average number of tracklets we exclude the 10% of MCMs with most entries
+        nActiveMcms90++;
+        mean90 += trackletCounter[idx].first;
+      }
     }
   }
-  if (!nActiveMcms) {
+  if (!nActiveMcms90) {
     printf("ERROR: did not find any MCMs which sent tracklets. Aborting\n");
     return;
   }
-  mean /= nActiveMcms;
-  float noiseThreshold = 10 * mean;
 
-  for (const auto& counter : trackletCounter) {
-    if (counter.first > noiseThreshold) {
-      noiseMap.setIsNoisy(counter.second);
+  auto n = nTrkltsPerActiveMcm.size() / 2;
+  auto median = nTrkltsPerActiveMcm[n]; // the distinction between odd/even number of entries in the vector is irrelevant
+  mean90 /= nActiveMcms90;
+  auto noiseThreshold = (useMean) ? 100.f * mean90 : (float)totalTriggerCounter / 10.f;
+
+  setlocale(LC_NUMERIC, "en_US.utf-8");
+  printf("Info: Found in total %'lu MCMs which sent tracklets with a median of %i tracklets per MCM\n", nTrkltsPerActiveMcm.size(), median);
+  printf("Info: Excluding the 10%% of MCMs which sent the highest number of tracklets %'i MCMs remain with on average %.2f tracklets per MCM\n", nActiveMcms90, mean90);
+  printf("Important: Masking MCMs which sent more than %.2f tracklets for given period\n", noiseThreshold);
+
+  std::vector<int> nTrackletsFromNoisyMcm;
+  bool hasPrintedFirstNoisy = false;
+  for (int idx = 0; idx < constants::MAXHALFCHAMBER * constants::NMCMHCMAX; ++idx) {
+    auto mcmIdx = trackletCounter[idx].second;
+    if (trackletCounter[idx].first > noiseThreshold) {
+      if (!hasPrintedFirstNoisy) {
+        printf("Info: The first masked MCM idx(%i) with glb idx %i sent %i trackelts\n", idx, mcmIdx, trackletCounter[idx].first);
+        hasPrintedFirstNoisy = true;
+      }
+      noiseMap.setIsNoisy(mcmIdx);
+      nTrackletsFromNoisyMcm.push_back(trackletCounter[idx].first);
     }
   }
+  if (noiseMap.getNumberOfNoisyMCMs() > 0) {
+    printf("Info: Last masked MCM sent %i tracklets\n", nTrackletsFromNoisyMcm.back());
+  }
 
   if (enableWriteCCDB) {
+    printf("Info: Uploading to CCDB (by default only to ccdb-test)\n");
     o2::ccdb::CcdbApi ccdb;
     ccdb.init("http://ccdb-test.cern.ch:8080");
     std::map<std::string, std::string> metadata;
@@ -98,7 +140,31 @@ void checkNoisyMCMs(std::string inpFile = "trdtracklets.root", bool enableWriteC
     auto timeStampEnd = timeStampStart;
     timeStampEnd += 1e3 * 60 * 60 * 24 * 60; // 60 days
     ccdb.storeAsTFileAny(&noiseMap, "TRD/Calib/NoiseMapMCM", metadata, timeStampStart, timeStampEnd);
+  } else {
+    // write to local file
+    printf("Info: Writing to local file mcmNoiseMap.root\n");
+    auto fOut = new TFile("mcmNoiseMap.root", "recreate");
+    fOut->WriteObjectAny(&noiseMap, "o2::trd::NoiseStatusMCM", "map");
+    fOut->Close();
   }
 
-  printf("Found in total %lu noisy MCMs for %lu tracklets from %lu trigger records\n", noiseMap.getNumberOfNoisyMCMs(), totalTrackletCounter, trigIn.size());
+  printf("Info: Found in total %lu noisy MCMs for %'lu tracklets from %'lu triggers\n", noiseMap.getNumberOfNoisyMCMs(), totalTrackletCounter, totalTriggerCounter);
+
+  if (detailedOutput) {
+    size_t countTrackletsFromNoisyMcms = 0;
+    printf("Info: Number of tracklets sent per masked MCM: \n");
+    printf("----->\n");
+    for (int idx = 0; idx < constants::MAXHALFCHAMBER * constants::NMCMHCMAX; ++idx) {
+      auto mcmIdx = trackletCounter[idx].second;
+      if (noiseMap.getIsNoisy(mcmIdx)) {
+        countTrackletsFromNoisyMcms += trackletCounter[idx].first;
+        int hcid, rob, mcm;
+        NoiseStatusMCM::convertMcmIdxGlb(mcmIdx, hcid, rob, mcm);
+        printf("Masked MCM idx (%i), glbIdx(%i). HCID(%i), ROB(%i), MCM(%i). nTracklets: %i\n", idx, mcmIdx, hcid, rob, mcm, trackletCounter[idx].first);
+      }
+    }
+    printf("\n<-----\n");
+    printf("Info: Number of tracklets sent from masked MCMs: %'lu (%.2f%%)\n", countTrackletsFromNoisyMcms, (float)countTrackletsFromNoisyMcms / totalTrackletCounter * 100);
+  }
+  printf("Info: Done\n");
 }

Detectors/TRD/calibration/src/TrackBasedCalib.cxx

@@ -19,7 +19,6 @@
 #include "DetectorsBase/GeometryManager.h"
 #include "TRDBase/Geometry.h"
 #include "TRDBase/PadPlane.h"
-#include "CCDB/BasicCCDBManager.h"
 #include "CommonUtils/NameConf.h"
 #include <fairlogger/Logger.h>
 
@@ -34,12 +33,6 @@ void TrackBasedCalib::reset()
 void TrackBasedCalib::init()
 {
   mRecoParam.setBfield(o2::base::Propagator::Instance()->getNominalBz());
-  auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
-  ccdbmgr.setURL(o2::base::NameConf::getCCDBServer());
-  mNoiseCalib = ccdbmgr.get<o2::trd::NoiseStatusMCM>("TRD/Calib/NoiseMapMCM");
-  if (!mNoiseCalib) {
-    LOG(error) << "Map of noisy MCMs not available";
-  }
 }
 
 void TrackBasedCalib::setInput(const o2::globaltracking::RecoContainer& input)
@@ -69,6 +62,10 @@ void TrackBasedCalib::calculateAngResHistos()
     return;
   }
 
+  if (!mNoiseCalib) {
+    LOG(alarm) << "No MCM noise map available. Please upload valid object to CCDB.";
+  }
+
   LOGF(info, "As input tracks are available: %lu ITS-TPC-TRD tracks and %lu TPC-TRD tracks", mTracksInITSTPCTRD.size(), mTracksInTPCTRD.size());
 
   int nTracksSuccessITSTPCTRD = doTrdOnlyTrackFits(mTracksInITSTPCTRD);
@@ -99,7 +96,7 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       if (trkWork.getTrackletIndex(iLayer) == -1) {
         continue;
       }
-      if (mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
+      if (mNoiseCalib && mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
         // ignore tracklets which originate from noisy MCMs
         continue;
       }
@@ -118,7 +115,7 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       if (trkWork.getTrackletIndex(iLayer) == -1) {
         continue;
       }
-      if (mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
+      if (mNoiseCalib && mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
         // ignore tracklets which originate from noisy MCMs
         continue;
       }
@@ -137,7 +134,7 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span<const TrackTRD>& tracks)
       if (trkWork.getTrackletIndex(iLayer) == -1) {
         continue;
       }
-      if (mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
+      if (mNoiseCalib && mNoiseCalib->isTrackletFromNoisyMCM(mTrackletsRaw[trkWork.getTrackletIndex(iLayer)])) {
         // ignore tracklets which originate from noisy MCMs
         continue;
       }

Detectors/TRD/workflow/src/TrackBasedCalibSpec.cxx

@@ -17,6 +17,7 @@
 #include "TRDCalibration/TrackBasedCalib.h"
 #include "Framework/Task.h"
 #include "Framework/ConfigParamRegistry.h"
+#include "Framework/CCDBParamSpec.h"
 #include "DetectorsBase/GeometryManager.h"
 #include "DetectorsBase/Propagator.h"
 #include "CommonUtils/NameConf.h"
@@ -40,9 +41,12 @@ class TRDTrackBasedCalibDevice : public Task
   ~TRDTrackBasedCalibDevice() override = default;
   void init(InitContext& ic) final;
   void run(ProcessingContext& pc) final;
+  void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) final;
   void endOfStream(framework::EndOfStreamContext& ec) final;
 
  private:
+  void updateTimeDependentParams(framework::ProcessingContext& pc);
+
   std::shared_ptr<DataRequest> mDataRequest;
   TrackBasedCalib mCalibrator; // gather input data for calibration of vD, ExB and gain
   std::unique_ptr<Output> mOutput;
@@ -61,6 +65,7 @@ void TRDTrackBasedCalibDevice::init(InitContext& ic)
 
 void TRDTrackBasedCalibDevice::run(ProcessingContext& pc)
 {
+  updateTimeDependentParams(pc);
   if (!mDataHeaderSet) {
     mOutput = std::make_unique<Output>(o2::header::gDataOriginTRD, "ANGRESHISTS", 0, Lifetime::Timeframe);
     mDataHeaderSet = true;
@@ -78,6 +83,19 @@ void TRDTrackBasedCalibDevice::run(ProcessingContext& pc)
   }
 }
 
+void TRDTrackBasedCalibDevice::updateTimeDependentParams(ProcessingContext& pc)
+{
+  pc.inputs().get<o2::trd::NoiseStatusMCM*>("mcmnoisemap"); // just to trigger the finaliseCCDB
+}
+
+void TRDTrackBasedCalibDevice::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
+{
+  if (matcher == ConcreteDataMatcher("TRD", "MCMNOISEMAP", 0)) {
+    LOG(info) << "NoiseStatusMCM object has been updated";
+    mCalibrator.setNoiseMapMCM((const o2::trd::NoiseStatusMCM*)obj);
+  }
+}
+
 void TRDTrackBasedCalibDevice::endOfStream(EndOfStreamContext& ec)
 {
   if (mNumberOfProcessedTFs > 0) {
@@ -105,11 +123,14 @@ DataProcessorSpec getTRDTrackBasedCalibSpec(o2::dataformats::GlobalTrackID::mask
   dataRequest->requestTracks(srcTrk, false);
   dataRequest->requestClusters(srcClu, false);
 
+  auto& inputs = dataRequest->inputs;
+  inputs.emplace_back("mcmnoisemap", "TRD", "MCMNOISEMAP", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/NoiseMapMCM"));
+
   outputs.emplace_back(o2::header::gDataOriginTRD, "ANGRESHISTS", 0, Lifetime::Timeframe);
 
   return DataProcessorSpec{
     "trd-trackbased-calib",
-    dataRequest->inputs,
+    inputs,
     outputs,
     AlgorithmSpec{adaptFromTask<TRDTrackBasedCalibDevice>(dataRequest)},
     Options{}};