diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9f450c880316d42cc974484c257d85cecd07c3c1..ee2ec3a8e6ad07fdd7ec95666893d5e2854b50c4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -184,6 +184,11 @@ if(fmt_FOUND)
list(APPEND packages fmt)
endif()
+#Searching for ONNXRuntime
+find_package(ONNXRuntime CONFIG)
+if (ONNXRuntime_FOUND)
+ list(APPEND packages ONNXRuntime)
+endif()
#Searching for FairMQ
find_package(FairMQ CONFIG REQUIRED)
diff --git a/core/base/CbmHistManager.cxx b/core/base/CbmHistManager.cxx
index c16b33d9df40b9498097edda57db1b81c0dc1293..1cc8bc50aa1e21b3e431495d4466fb94d3c4bf37 100755
--- a/core/base/CbmHistManager.cxx
+++ b/core/base/CbmHistManager.cxx
@@ -127,6 +127,13 @@ void CbmHistManager::WriteToFile()
}
}
+void CbmHistManager::WriteCanvasToFile()
+{
+ for (auto& canvas : fCanvases) {
+ canvas->Write();
+ }
+}
+
void CbmHistManager::ReadFromFile(TFile* file)
{
assert(file != nullptr);
diff --git a/core/base/CbmHistManager.h b/core/base/CbmHistManager.h
index fae5f81a52d80b7cd59a0a0e2cd12784f151584b..f5247fe88e8346efa5c01d773fa35ca3a08d5a9a 100644
--- a/core/base/CbmHistManager.h
+++ b/core/base/CbmHistManager.h
@@ -429,6 +429,11 @@ public:
*/
void WriteToFile();
+ /**
+ * \brief Write all canvas to current opened file.
+ */
+ void WriteCanvasToFile();
+
/**
* \brief Read histograms from file.
* \param[in] file Pointer to file with histograms.
diff --git a/core/data/rich/CbmRichHit.cxx b/core/data/rich/CbmRichHit.cxx
index 93de285ad09bf2253a8466ad7b53317a47311cc3..5a63c6827456fe0a3a3c8a018861919d53e9753c 100644
--- a/core/data/rich/CbmRichHit.cxx
+++ b/core/data/rich/CbmRichHit.cxx
@@ -22,6 +22,7 @@ CbmRichHit::CbmRichHit()
// fNPhotons(0),
// fAmplitude(0.),
fToT(0.)
+ , fIsNoiseNN(false)
{
SetType(kRICHHIT);
SetTime(0.);
@@ -34,6 +35,7 @@ CbmRichHit::CbmRichHit(double x, double y)
// fNPhotons(0),
// fAmplitude(0.),
fToT(0.)
+ , fIsNoiseNN(false)
{
SetType(kRICHHIT);
SetX(x);
@@ -48,6 +50,7 @@ CbmRichHit::CbmRichHit(double x, double y, double ts, double tot)
// fNPhotons(0),
// fAmplitude(0.),
fToT(tot)
+ , fIsNoiseNN(false)
{
SetType(kRICHHIT);
SetX(x);
@@ -61,8 +64,8 @@ std::string CbmRichHit::ToString() const
{
stringstream ss;
ss << "CbmRichHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=("
- << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() << " refId="
- << GetRefId()
+ << GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() << " refId=" << GetRefId() << " isNoiseNN="
+ << GetIsNoiseNN()
// << " pmtId=" << GetPmtId() << " nofPhotons=" << GetNPhotons()
// << " amplitude=" << GetAmplitude()
<< endl;
diff --git a/core/data/rich/CbmRichHit.h b/core/data/rich/CbmRichHit.h
index b31e48c37d0090e1a24b0ec93f3a72f1383adc16..9e0dc922f8dc795c0335deaa48e82b48c654987e 100644
--- a/core/data/rich/CbmRichHit.h
+++ b/core/data/rich/CbmRichHit.h
@@ -58,12 +58,14 @@ public:
//virtual void SetNPhotons (int32_t n) { fNPhotons = n; }
//virtual void SetAmplitude(double amp) { fAmplitude = amp; }
void SetToT(double tot) { fToT = tot; }
+ void SetIsNoiseNN(bool isNoiseNN) { fIsNoiseNN = isNoiseNN; }
/** Accessors **/
virtual int32_t GetPmtId() const { return fPmtId; }
//virtual int32_t GetNPhotons() const { return fNPhotons; }
//virtual double GetAmplitude() const { return fAmplitude; }
double GetToT() const { return fToT; }
+ bool GetIsNoiseNN() const { return fIsNoiseNN; }
/** Outdated. Use CbmHit::GetTime() and SetTime() instead. **/
// double GetTimestamp() const { return GetTime(); }
@@ -75,8 +77,11 @@ private:
//double fAmplitude; // hit amplitude
double fToT; // hit time-over-threshold
+ // Flag for mRICH noise hit classification
+ // if true -> hit is classified as noise and excluded from ringfinding
+ bool fIsNoiseNN;
- ClassDef(CbmRichHit, 3)
+ ClassDef(CbmRichHit, 4)
};
#endif // CBMRICHHIT_H_
diff --git a/reco/detectors/rich/CMakeLists.txt b/reco/detectors/rich/CMakeLists.txt
index becc3e386c136430fc70b1efdd8ecc87c5abab77..7b3c900478f889094e6fd994bcbd63a22ea473a4 100644
--- a/reco/detectors/rich/CMakeLists.txt
+++ b/reco/detectors/rich/CMakeLists.txt
@@ -25,6 +25,8 @@ set(SRCS
mcbm/CbmRichMCbmAerogelAna.cxx
mcbm/CbmRichMCbmToTShifter.cxx
mcbm/CbmRichMCbmSEDisplay.cxx
+ mcbm/CbmRichMCbmDenoiseCnn.cxx
+ mcbm/CbmRichMCbmDenoiseQa.cxx
qa/CbmRichRingFitterQa.cxx
qa/CbmRichGeoOpt.cxx
@@ -118,6 +120,14 @@ set(PRIVATE_DEPENDENCIES
ROOT::Minuit2
)
+if(ONNXRuntime_FOUND)
+ ADD_DEFINITIONS(-DHAVE_ONNXRUNTIME)
+ Message(STATUS "Rich will be compiled with ONNXRuntime support")
+ list(APPEND PRIVATE_DEPENDENCIES ONNXRuntime::ONNXRuntime)
+else()
+ Message(STATUS "Rich will NOT be compiled with ONNXRuntime support")
+endif()
+
set(INTERFACE_DEPENDENCIES
FairLogger::FairLogger
external::fles_ipc
diff --git a/reco/detectors/rich/CbmRichRecoLinkDef.h b/reco/detectors/rich/CbmRichRecoLinkDef.h
index a832c0dd25658ce645712bbbd8b6586ee53ffd7c..8423c69103f74b8e3b82ab0ce906741b836e0a95 100644
--- a/reco/detectors/rich/CbmRichRecoLinkDef.h
+++ b/reco/detectors/rich/CbmRichRecoLinkDef.h
@@ -22,6 +22,10 @@
#pragma link C++ class CbmRichMCbmAerogelAna;
#pragma link C++ class CbmRichMCbmToTShifter;
#pragma link C++ class CbmRichMCbmSEDisplay;
+#ifdef HAVE_ONNXRUNTIME
+#pragma link C++ class CbmRichMCbmDenoiseCnn +;
+#endif
+#pragma link C++ class CbmRichMCbmDenoiseQa +;
//qa
#pragma link C++ class CbmRichGeoTest + ;
diff --git a/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx b/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
index d0474d45b00eb286b44b6baf595f162d8f88db91..c887cd71243b57ab12bb48c449e2018494604b05 100644
--- a/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
+++ b/reco/detectors/rich/finder/CbmRichRingFinderHough.cxx
@@ -82,7 +82,7 @@ Int_t CbmRichRingFinderHough::DoFind(CbmEvent* event, TClonesArray* rHitArray, T
for (Int_t iH0 = 0; iH0 < nofRichHits; iH0++) {
Int_t iH = event ? event->GetIndex(ECbmDataType::kRichHit, iH0) : iH0;
CbmRichHit* hit = static_cast<CbmRichHit*>(rHitArray->At(iH));
- if (hit != nullptr) {
+ if (hit != nullptr && !hit->GetIsNoiseNN()) {
CbmRichHoughHit tempPoint;
tempPoint.fHit.fX = hit->GetX();
tempPoint.fHit.fY = hit->GetY();
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5ecebd00b7c392764e53dde29903fef5a023c119
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.cxx
@@ -0,0 +1,137 @@
+/* Copyright (C) 2024 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+#if HAVE_ONNXRUNTIME
+
+#include "CbmRichMCbmDenoiseCnn.h"
+
+#include "CbmDigiManager.h"
+#include "CbmEvent.h"
+#include "CbmRichDetectorData.h"
+#include "CbmRichDigiMapManager.h"
+#include "CbmRichHit.h"
+
+#include <Logger.h>
+
+#include <TClonesArray.h>
+#include <TStopwatch.h>
+
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+#include <onnxruntime/core/session/onnxruntime_cxx_api.h>
+
+void CbmRichMCbmDenoiseCnn::Init()
+{
+ fOrtEnv = std::make_unique<Ort::Env>(ORT_LOGGING_LEVEL_WARNING, GetName());
+
+ fOrtSessionOptions = std::make_unique<Ort::SessionOptions>();
+ fOrtSession = std::make_unique<Ort::Session>(*fOrtEnv, fOnnxFilePath.c_str(), *fOrtSessionOptions);
+
+ fOrtRunOptions = std::make_unique<Ort::RunOptions>(nullptr);
+ fOrtAllocatorInfo = std::make_unique<Ort::MemoryInfo>(Ort::MemoryInfo::CreateCpu(OrtDeviceAllocator, OrtMemTypeCPU));
+ fOrtInput = std::make_unique<Ort::Value>(Ort::Value::CreateTensor<float>(
+ *fOrtAllocatorInfo, fInput.data(), fInput.size(), fInputShape.data(), fInputShape.size()));
+
+ // Thread numbers need to be explicitly set. With the current Ort version it can't detect
+ // the numbers on the cluster, should be fixed in a later version.
+ fOrtSessionOptions->SetIntraOpNumThreads(1);
+ fOrtSessionOptions->SetInterOpNumThreads(1);
+ fOrtSessionOptions->SetExecutionMode(ORT_SEQUENTIAL);
+ fOrtSessionOptions->SetGraphOptimizationLevel(ORT_ENABLE_ALL);
+
+ fCbmRichDigiMapManager = &CbmRichDigiMapManager::GetInstance();
+}
+
+void CbmRichMCbmDenoiseCnn::Process(CbmEvent* event, const TClonesArray* richHits)
+{
+ int nHits = event ? static_cast<int>(event->GetNofData(ECbmDataType::kRichHit)) : richHits->GetEntriesFast();
+ std::vector<int> timeWindowHitIndices;
+ for (int i = 0; i < nHits; i++) { // Sliding time window loop
+ timeWindowHitIndices.clear();
+ int seedIdx = event ? static_cast<int>(event->GetIndex(ECbmDataType::kRichHit, static_cast<uint32_t>(i))) : i;
+ CbmRichHit* seedHit = static_cast<CbmRichHit*>(richHits->At(seedIdx));
+ if (!seedHit) continue;
+ timeWindowHitIndices.push_back(seedIdx);
+ int hitsInTimeWindow = 1;
+ for (int j = i + 1; j < nHits; j++) { // Search for hits in time window
+ int hitIdx = event ? static_cast<int>(event->GetIndex(ECbmDataType::kRichHit, j)) : j;
+ CbmRichHit* hit = static_cast<CbmRichHit*>(richHits->At(hitIdx));
+ if (!hit) continue;
+ double dt = hit->GetTime() - seedHit->GetTime();
+ if (dt < fTimeWindowLength) {
+ hitsInTimeWindow++;
+ timeWindowHitIndices.push_back(hitIdx);
+ }
+ else {
+ break;
+ }
+ }
+
+ if (hitsInTimeWindow >= fMinHitsInTimeWindow) {
+ ProcessTimeWindow(timeWindowHitIndices, seedHit->GetTime(), richHits);
+ i += hitsInTimeWindow - 1; // Move to last hit inside time window
+ }
+ else {
+ seedHit->SetIsNoiseNN(true);
+ }
+ }
+}
+
+void CbmRichMCbmDenoiseCnn::ProcessTimeWindow(const std::vector<int>& timeWindowHitIndices, const double& seedHitTime,
+ const TClonesArray* richHits)
+{
+ fInput = {}; // Reset all input values to 0.0
+ std::vector<int> gridIndices;
+ gridIndices.reserve(50);
+ for (const auto& hitIdx : timeWindowHitIndices) {
+ CbmRichHit* hit = static_cast<CbmRichHit*>(richHits->At(hitIdx));
+ const auto gridIdx = AddressToGridIndex(hit->GetAddress());
+ if (gridIdx < 0 || gridIdx >= static_cast<int>(fInput.size())) {
+ LOG(error) << GetName() << "::ProcessTimeWindow: Invalid grid index: " << gridIdx;
+ continue;
+ }
+ gridIndices.push_back(gridIdx);
+ // Shift time by 1ns to distinguish empty pixels from seed hit time
+ fInput[gridIdx] = static_cast<float>((hit->GetTime() - seedHitTime + 1.0) / (fTimeWindowLength + 1.0));
+ }
+
+ const auto output = Inference(gridIndices);
+
+ for (std::size_t i = 0; i < timeWindowHitIndices.size(); i++) {
+ CbmRichHit* hit = static_cast<CbmRichHit*>(richHits->At(timeWindowHitIndices[i]));
+ bool isNoise = output[i] < fClassificationThreshold;
+ hit->SetIsNoiseNN(isNoise);
+ }
+}
+
+std::vector<float> CbmRichMCbmDenoiseCnn::Inference(const std::vector<int>& gridIndices)
+{
+ auto output = fOrtSession->Run(*fOrtRunOptions.get(), fInputNames.data(), fOrtInput.get(), 1, fOutputNames.data(), 1);
+ float* intarr = output.front().GetTensorMutableData<float>();
+
+ std::vector<float> out(gridIndices.size());
+ std::transform(gridIndices.begin(), gridIndices.end(), out.begin(),
+ [intarr](int gridIdx) { return intarr[gridIdx]; });
+ return out;
+}
+
+int CbmRichMCbmDenoiseCnn::AddressToGridIndex(int address)
+{
+ CbmRichPixelData* pixel_data = fCbmRichDigiMapManager->GetPixelDataByAddress(address);
+
+ // Calculate local X [0,7],Y [0,7] indices of one MAPMT
+ int pmtUID = pixel_data->fPmtId;
+ int pmtIndX = (pmtUID >> 4) & 0xF; // index ascending from -x to x
+ int pmtIndY = pmtUID & 0xF; // index ascending from y to -y
+
+ // Calculate global X [0,31],Y [0,72] indices
+ int globalIndX = 8 * pmtIndX + pixel_data->fPixelId % 8;
+ int globalIndY = 8 * pmtIndY + pixel_data->fPixelId / 8;
+
+ return 32 * globalIndY + globalIndX;
+}
+
+#endif // HAVE_ONNXRUNTIME
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.h b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.h
new file mode 100644
index 0000000000000000000000000000000000000000..1e6a44be9b23e5d5a0d5465e042a6a2b0b9c387f
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseCnn.h
@@ -0,0 +1,121 @@
+/* Copyright (C) 2024 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+/**
+* \file CbmRichMCbmDenoiseCnn.h
+* \author M.Beyer
+* \date 2024
+**/
+
+#if HAVE_ONNXRUNTIME
+
+#pragma once
+
+#include <TObject.h>
+#include <TSystem.h>
+
+#include <array>
+#include <string>
+#include <vector>
+
+#include <onnxruntime/core/session/onnxruntime_cxx_api.h>
+
+class CbmEvent;
+class CbmRichDigiMapManager;
+class TClonesArray;
+
+/**
+* \class CbmRichMCbmDenoiseCnn
+*
+* \brief MCbm mRICH noise removal using a CNN
+*
+* \author M.Beyer
+* \date 2024
+**/
+class CbmRichMCbmDenoiseCnn : public TObject {
+ public:
+ /** Default constructor */
+ CbmRichMCbmDenoiseCnn() = default;
+
+ /** Destructor */
+ ~CbmRichMCbmDenoiseCnn() = default;
+
+ /** Copy constructor (disabled) */
+ CbmRichMCbmDenoiseCnn(const CbmRichMCbmDenoiseCnn&) = delete;
+
+ /** Assignment operator (disabled) */
+ CbmRichMCbmDenoiseCnn operator=(const CbmRichMCbmDenoiseCnn&) = delete;
+
+ /** Initialize ONNX Runtime structures and .onnx model */
+ void Init();
+
+ /**
+ * \brief Process time sorted RICH hits and created non overlapping sliding time windows
+ * \param event if CbmEvent* is nullptr -> process on full Ts
+ * \param richHits TClonesArray of RICH hits
+ */
+ void Process(CbmEvent* event, const TClonesArray* richHits);
+
+ /**
+ * \brief Preparing hits in a non overlapping sliding time window and doing inference.
+ * Fills the classification result as a flag into CbmRichHit.fIsNoiseNN
+ * \param timeWindowHitIndices vector of hit indices from richHits TClonesArray
+ * \param seedHitTime time of seed hit
+ * \param richHits TClonesArray of RICH hits
+ */
+ void ProcessTimeWindow(const std::vector<int>& timeWindowHitIndices, const double& seedHitTime,
+ const TClonesArray* richHits);
+
+ /**
+ * \brief Mapping from pixel address to row major grid index
+ * \param int Digi/Hit address
+ * \return Grid index [0, 2303]
+ */
+ int AddressToGridIndex(const int address);
+
+ /**
+ * \brief ONNX runtime inference on input array fInput filled in ProcessTimeWindow(...)
+ * \param gridIndices vector of grid indices [0, 2303]
+ * \return vector of classification outputs for the given grid indices
+ */
+ std::vector<float> Inference(const std::vector<int>& gridIndices);
+
+ /**
+ * \brief Set the classification threshold
+ * \param threshold value between (0.0, 1.0). Larger values will remove hits more aggressively
+ */
+ void SetClassificationThreshold(float threshold) { fClassificationThreshold = threshold; }
+
+ private:
+ const double fTimeWindowLength{25.};
+ const int fMinHitsInTimeWindow{7};
+
+ float fClassificationThreshold{0.5};
+
+ const std::string fOnnxFilePath =
+ std::string(gSystem->Getenv("VMCWORKDIR")) + "/parameters/rich/mRich/rich_v21c_mcbm_UNet.onnx";
+
+ // TODO: read names/shapes from onnx file
+ const std::vector<const char*> fInputNames{"input"};
+ const std::vector<const char*> fOutputNames{"output"};
+
+ const std::array<int64_t, 4> fInputShape{1, 1, 72, 32};
+ const std::array<int64_t, 4> fOutputShape{1, 1, 72, 32};
+
+ std::array<float, 2304> fInput{};
+
+ std::unique_ptr<Ort::Env> fOrtEnv{nullptr};
+ std::unique_ptr<Ort::SessionOptions> fOrtSessionOptions{nullptr};
+ std::unique_ptr<Ort::Session> fOrtSession{nullptr};
+ std::unique_ptr<Ort::MemoryInfo> fOrtAllocatorInfo{nullptr};
+ std::unique_ptr<Ort::RunOptions> fOrtRunOptions{nullptr};
+
+ std::unique_ptr<Ort::Value> fOrtInput{nullptr};
+
+ CbmRichDigiMapManager* fCbmRichDigiMapManager{nullptr};
+
+ ClassDef(CbmRichMCbmDenoiseCnn, 1);
+};
+
+#endif
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fe83df3ec16be42145c81dd50687f093e5f699af
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.cxx
@@ -0,0 +1,198 @@
+/* Copyright (C) 2024 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+#include "CbmRichMCbmDenoiseQa.h"
+
+#include "CbmEvent.h"
+#include "CbmHistManager.h"
+#include "CbmRichHit.h"
+#include "CbmRichRing.h"
+
+#include <FairRootManager.h>
+#include <Logger.h>
+
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TDirectory.h>
+#include <TEllipse.h>
+#include <TFile.h>
+#include <TH2D.h>
+
+InitStatus CbmRichMCbmDenoiseQa::Init()
+{
+ FairRootManager* manager = FairRootManager::Instance();
+ if (!manager) LOG(fatal) << GetName() << "::Init: No FairRootManager";
+
+ fCbmEvents = static_cast<TClonesArray*>(manager->GetObject("CbmEvent"));
+ LOG(info) << GetName() << "::Init: CbmEvent array " << (fCbmEvents ? "found" : "not found");
+
+ fRichHits = static_cast<TClonesArray*>(manager->GetObject("RichHit"));
+ if (!fRichHits) LOG(fatal) << GetName() << "::Init: No RichHit array!";
+
+ fRichRings = static_cast<TClonesArray*>(manager->GetObject("RichRing"));
+ if (!fRichRings) LOG(fatal) << GetName() << "::Init: No RichRing array!";
+
+ InitHistograms();
+
+ return kSUCCESS;
+}
+
+void CbmRichMCbmDenoiseQa::InitHistograms()
+{
+ fHM = std::make_unique<CbmHistManager>();
+ fHMSed = std::make_unique<CbmHistManager>();
+
+ // 1D Histograms
+ fHM->Create1<TH1D>("fhRingsPerEvent", "number of ring per event;# rings per event;Entries", 6, -0.5, 5.5);
+
+ fHM->Create1<TH1D>("fhRingRadius", "ring radius;ring radius [cm];Entries", 100, 0., 7.);
+
+ fHM->Create1<TH1D>("fhRingRadiusUp", "ring radius upper half;ring radius [cm]; count", 100, 0., 7.);
+
+ fHM->Create1<TH1D>("fhRingRadiusDown", "ring radius lower half;ring radius [cm]; count", 100, 0., 7.);
+
+ fHM->Create1<TH1D>("fhRingHits", "number of ring hits;# ring hits;Entries", 50, -0.5, 49.5);
+
+ fHM->Create1<TH1D>("fhRingHitsUp", "number of ring hits upper half;# ring hits;Entries", 50, -0.5, 49.5);
+
+ fHM->Create1<TH1D>("fhRingHitsDown", "number of ring hits lower half;# ring hits;Entries", 50, -0.5, 49.5);
+
+ fHM->Create1<TH1D>("fhRingHitsTimeDiff", "ringhit to ring time difference;hittime-ringtime [ns]; count", 100, -10.,
+ 10.);
+
+ // 2D Histograms
+ fHM->Create2<TH2D>("fhRichHitsXY", "position of rich hits;X [cm];Y [cm];Entries", 37, -11.2, 11.2, 82, -24.115,
+ 23.885);
+
+ fHM->Create2<TH2D>("fhRingHitsXY", "position of ring hits;X [cm];Y [cm];Entries", 37, -11.2, 11.2, 82, -24.115,
+ 23.885);
+
+ fHM->Create2<TH2D>("fhRingCenterXY", "position of ring centers;X [cm];Y [cm];Entries", 37, -11.2, 11.2, 82, -24.115,
+ 23.885);
+}
+
+void CbmRichMCbmDenoiseQa::Exec(Option_t* /*option*/)
+{
+ LOG(debug) << GetName() << " TS " << fTsNum;
+ fTsNum++;
+
+ if (fCbmEvents) {
+ Int_t nEvents = fCbmEvents->GetEntriesFast();
+ for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ CbmEvent* event = static_cast<CbmEvent*>(fCbmEvents->At(iEvent));
+ if (!event) continue;
+ if (event->GetNofData(ECbmDataType::kRichRing) > 0) {
+ fHM->H1("fhRingsPerEvent")->Fill(static_cast<double>(event->GetNofData(ECbmDataType::kRichRing)));
+ if (fnSEDs < fMaxSEDs) {
+ DrawSED(event);
+ fnSEDs++;
+ }
+ }
+ Process(event);
+ fEventNum++;
+ }
+ }
+ else {
+ Process(nullptr);
+ }
+}
+
+void CbmRichMCbmDenoiseQa::Process(CbmEvent* event)
+{
+ // RICH Hits
+ Int_t nRichHits = event ? static_cast<Int_t>(event->GetNofData(ECbmDataType::kRichHit)) : fRichHits->GetEntriesFast();
+ for (int i = 0; i < nRichHits; i++) {
+ Int_t hitIndex = event ? static_cast<Int_t>(event->GetIndex(ECbmDataType::kRichHit, static_cast<uint32_t>(i))) : i;
+ CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(hitIndex));
+ if (!hit) continue;
+ fHM->H2("fhRichHitsXY")->Fill(hit->GetX(), hit->GetY());
+ }
+
+ // RICH Rings
+ Int_t nRichRings =
+ event ? static_cast<Int_t>(event->GetNofData(ECbmDataType::kRichRing)) : fRichRings->GetEntriesFast();
+ for (int i = 0; i < nRichRings; i++) {
+ Int_t ringIndex =
+ event ? static_cast<Int_t>(event->GetIndex(ECbmDataType::kRichRing, static_cast<uint32_t>(i))) : i;
+ CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(ringIndex));
+ if (!ring) continue;
+ fHM->H1("fhRingRadius")->Fill(ring->GetRadius());
+ fHM->H1(std::string("fhRingRadius") + std::string((ring->GetCenterY() > 0) ? "Up" : "Down"))
+ ->Fill(ring->GetRadius());
+ fHM->H1("fhRingHits")->Fill(ring->GetNofHits());
+ fHM->H1(std::string("fhRingHits") + std::string((ring->GetCenterY() > 0) ? "Up" : "Down"))
+ ->Fill(ring->GetNofHits());
+ fHM->H2("fhRingCenterXY")->Fill(ring->GetCenterX(), ring->GetCenterY());
+ // Ring hits
+ Int_t nRingHits = ring->GetNofHits();
+ for (int j = 0; j < nRingHits; j++) {
+ CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(static_cast<Int_t>(ring->GetHit(j))));
+ if (!hit) continue;
+ fHM->H1("fhRingHitsTimeDiff")->Fill(hit->GetTime() - ring->GetTime());
+ fHM->H2("fhRingHitsXY")->Fill(hit->GetX(), hit->GetY());
+ }
+ }
+}
+
+void CbmRichMCbmDenoiseQa::DrawSED(CbmEvent* event)
+{
+ std::string hName = "fhSED" + std::to_string(fEventNum);
+ auto c = fHMSed->CreateCanvas(hName, hName, 800, 1505);
+ c->SetMargin(0.15, 0.15, 0.1, 0.1);
+ c->cd();
+ TH2D* h = new TH2D(hName.c_str(), hName.c_str(), 37, -11.2, 11.2, 82, -24.115, 23.885);
+ h->SetXTitle("X [cm]");
+ h->SetYTitle("Y [cm]");
+ h->SetZTitle("NN classification 0: noise, 1: signal");
+ h->GetZaxis()->SetTitleOffset(1.3);
+ fHMSed->Add(hName, h);
+
+ Int_t nHits = static_cast<Int_t>(event->GetNofData(ECbmDataType::kRichHit));
+ for (int i = 0; i < nHits; i++) {
+ Int_t hitIndex = static_cast<Int_t>(event->GetIndex(ECbmDataType::kRichHit, static_cast<uint32_t>(i)));
+ CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(hitIndex));
+ if (!hit) continue;
+ if (h->GetBinContent(h->FindBin(hit->GetX(), hit->GetY())) > 0) {
+ LOG(info) << GetName() << "::DrawSED channel already filled for this event. Skipping.";
+ continue;
+ }
+ if (!hit->GetIsNoiseNN()) {
+ h->Fill(hit->GetX(), hit->GetY(), 1.0);
+ }
+ else {
+ h->Fill(hit->GetX(), hit->GetY(), 0.001);
+ }
+ }
+ h->Draw("colz");
+
+ Int_t nRings = static_cast<Int_t>(event->GetNofData(ECbmDataType::kRichRing));
+ for (int i = 0; i < nRings; i++) {
+ Int_t ringIndex = static_cast<Int_t>(event->GetIndex(ECbmDataType::kRichRing, static_cast<uint32_t>(i)));
+ CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(ringIndex));
+ if (!ring) continue;
+ TEllipse* circle = new TEllipse(ring->GetCenterX(), ring->GetCenterY(), ring->GetAaxis(), ring->GetBaxis(), 0.,
+ 360., ring->GetPhi() * 180. / TMath::Pi());
+ circle->SetFillStyle(0);
+ circle->SetLineWidth(2);
+ circle->Draw("same");
+ }
+}
+
+void CbmRichMCbmDenoiseQa::Finish()
+{
+ TDirectory* oldir = gDirectory;
+ TFile* outFile = FairRootManager::Instance()->GetOutFile();
+ if (outFile) {
+ outFile->mkdir(GetName());
+ outFile->cd(GetName());
+ fHM->WriteToFile();
+ outFile->cd();
+ outFile->mkdir((std::string(GetName()) + std::string("SEDs")).c_str());
+ outFile->cd((std::string(GetName()) + std::string("SEDs")).c_str());
+ fHMSed->WriteCanvasToFile();
+ fHM->Clear();
+ fHMSed->Clear();
+ }
+ gDirectory->cd(oldir->GetPath());
+}
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.h b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.h
new file mode 100644
index 0000000000000000000000000000000000000000..098a2c820048e45853593c7041f1f8515ed9f2f8
--- /dev/null
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmDenoiseQa.h
@@ -0,0 +1,82 @@
+/* Copyright (C) 2024 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+/**
+* \file CbmRichMCbmDenoiseQa.h
+* \author M.Beyer
+* \date 2024
+**/
+
+#pragma once
+
+#include "FairTask.h"
+
+class CbmDigiManager;
+class CbmRichDigi;
+class CbmHistManager;
+class CbmEvent;
+class TClonesArray;
+
+/**
+* \class CbmRichMCbmDenoiseQa
+*
+* \brief QA for MCbm mRICH noise removal
+*
+* \author M.Beyer
+* \date 2024
+**/
+class CbmRichMCbmDenoiseQa : public FairTask {
+ public:
+ /** Default constructor */
+ CbmRichMCbmDenoiseQa() : FairTask("CbmRichMCbmDenoiseQa"){};
+
+ /** Destructor */
+ ~CbmRichMCbmDenoiseQa() = default;
+
+ /** Copy constructor (disabled) */
+ CbmRichMCbmDenoiseQa(const CbmRichMCbmDenoiseQa&) = delete;
+
+ /** Assignment operator (disabled) */
+ CbmRichMCbmDenoiseQa operator=(const CbmRichMCbmDenoiseQa&) = delete;
+
+ /** Inherited from FairTask */
+ InitStatus Init();
+
+ /** Inherited from FairTask */
+ void Exec(Option_t* option);
+
+ /** Inherited from FairTask */
+ void Finish();
+
+ /** Initialize histogram manager */
+ void InitHistograms();
+
+ /**
+ * \brief Process data and fill histograms
+ * \param event if CbmEvent* is nullptr -> process full Ts
+ */
+ void Process(CbmEvent* event);
+
+ /** Draw a Single-event-display */
+ void DrawSED(CbmEvent* event);
+
+ /** Set maximum number of SEDs drawn */
+ void SetMaxSEDs(Int_t maxSEDs) { fMaxSEDs = maxSEDs; }
+
+ private:
+ Int_t fTsNum{};
+ Int_t fEventNum{};
+ Int_t fnSEDs{};
+
+ Int_t fMaxSEDs{20};
+
+ std::unique_ptr<CbmHistManager> fHM{nullptr};
+ std::unique_ptr<CbmHistManager> fHMSed{nullptr}; // For single event displays
+
+ TClonesArray* fCbmEvents{nullptr};
+ TClonesArray* fRichHits{nullptr};
+ TClonesArray* fRichRings{nullptr};
+
+ ClassDef(CbmRichMCbmDenoiseQa, 1)
+};
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
index 1a86ebd6222ca1a38d6b83bc5f1785bb56473e96..7f69b82b4aac27a6055da6608cb012697104e401 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.cxx
@@ -12,21 +12,31 @@
#include "CbmRichGeoManager.h"
#include "CbmRichHit.h"
#include "CbmRichPoint.h"
+#if HAVE_ONNXRUNTIME
+#include "CbmRichMCbmDenoiseCnn.h"
+#endif
#include <FairRootManager.h>
#include <Logger.h>
-#include "TClonesArray.h"
+#include <TClonesArray.h>
+#include <TStopwatch.h>
+#include <iomanip>
#include <iostream>
+using std::fixed;
+using std::right;
+using std::setprecision;
+using std::setw;
+
using namespace std;
CbmRichMCbmHitProducer::CbmRichMCbmHitProducer()
: FairTask("CbmRichMCbmHitProducer")
, fRichHits(NULL)
- , fEventNum(0)
+ , fNofTs(0)
, fHitError(0.6 / sqrt(12.))
,
//fMappingFile("mRICH_Mapping_vert_20190318.geo")
@@ -66,6 +76,15 @@ InitStatus CbmRichMCbmHitProducer::Init()
a = 0.;
if (fDoICD) read_ICD(fICD_offset_read, 0);
+#if HAVE_ONNXRUNTIME
+ if (fUseDenoiseNN) {
+ fDenoiseCnn = std::make_unique<CbmRichMCbmDenoiseCnn>();
+ fDenoiseCnn->Init();
+ fDenoiseCnn->SetClassificationThreshold(fDenoiseCnnThreshold);
+ LOG(info) << "CbmRichMCbmHitProducer: Initialized CbmRichMCbmDenoiseCnn";
+ }
+#endif
+
return kSUCCESS;
}
@@ -109,23 +128,39 @@ void CbmRichMCbmHitProducer::InitMapping()
void CbmRichMCbmHitProducer::Exec(Option_t* /*option*/)
{
- fEventNum++;
- LOG(info) << "CbmRichMCbmHitProducer Event (or TimeSlice) " << fEventNum;
-
+ fNofDigis = 0;
+ TStopwatch timer;
+ timer.Start();
fRichHits->Delete();
// if CbmEvent does not exist then process standard event.
// if CbmEvent exists then proceed all events in time slice.
- Int_t nUnits = (fCbmEvents != nullptr) ? fCbmEvents->GetEntriesFast() : 1;
-
+ Int_t nUnits = fCbmEvents ? fCbmEvents->GetEntriesFast() : 1;
+ if (fCbmEvents) fNofEvents += nUnits;
for (Int_t iUnit = 0; iUnit < nUnits; iUnit++) {
- CbmEvent* event = (fCbmEvents != nullptr) ? static_cast<CbmEvent*>(fCbmEvents->At(iUnit)) : nullptr;
+ CbmEvent* event = fCbmEvents ? static_cast<CbmEvent*>(fCbmEvents->At(iUnit)) : nullptr;
ProcessData(event);
}
+ timer.Stop();
+ fTotalTime += timer.RealTime();
+
+ std::stringstream logOut;
+ logOut << setw(20) << left << GetName() << "[";
+ logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << fNofTs;
+ if (fCbmEvents) logOut << ", events " << nUnits;
+ logOut << ", digis " << fNofDigis;
+ logOut << ", hits " << fRichHits->GetEntriesFast();
+ LOG(info) << logOut.str();
+ fTotalNofDigis += fNofDigis;
+ fTotalNofHits += fRichHits->GetEntriesFast();
+ fNofTs++;
}
void CbmRichMCbmHitProducer::ProcessData(CbmEvent* event)
{
+ TStopwatch timer;
+ timer.Start();
if (event != NULL) {
LOG(debug) << "CbmRichMCbmHitProducer CbmEvent mode. CbmEvent # " << event->GetNumber();
Int_t nofDigis = event->GetNofData(ECbmDataType::kRichDigi);
@@ -145,6 +180,18 @@ void CbmRichMCbmHitProducer::ProcessData(CbmEvent* event)
ProcessDigi(event, iDigi);
}
}
+ timer.Stop();
+ fHitProducerTime += timer.RealTime();
+
+// Process all hits with DenoiseNN, setting value RichHit.fIsNoiseNN
+#if HAVE_ONNXRUNTIME
+ if (fUseDenoiseNN) {
+ timer.Start();
+ fDenoiseCnn->Process(event, fRichHits);
+ timer.Stop();
+ fDenoiseCnnTime += timer.RealTime();
+ }
+#endif
}
void CbmRichMCbmHitProducer::ProcessDigi(CbmEvent* event, Int_t digiIndex)
@@ -154,6 +201,7 @@ void CbmRichMCbmHitProducer::ProcessDigi(CbmEvent* event, Int_t digiIndex)
if (digi->GetAddress() < 0) return;
Int_t DiRICH_Add = (digi->GetAddress() >> 16) & 0xFFFF;
if (DiRICH_Add == 0x7901 || DiRICH_Add == 0x7902) return; // TRBaddresses 0x7901 and 0x7902 are for FSD/NCAL
+ fNofDigis++;
if (isInToT(digi->GetToT())) {
TVector3 posPoint;
CbmRichPixelData* data = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(digi->GetAddress());
@@ -212,7 +260,43 @@ void CbmRichMCbmHitProducer::AddHit(CbmEvent* event, TVector3& posHit, const Cbm
}
-void CbmRichMCbmHitProducer::Finish() { fRichHits->Clear(); }
+void CbmRichMCbmHitProducer::Finish()
+{
+ fRichHits->Clear();
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Time slices : " << fNofTs;
+ LOG(info) << "Digis / TS : " << fixed << setprecision(2) << fTotalNofDigis / Double_t(fNofTs);
+ LOG(info) << "Hits / TS : " << fixed << setprecision(2) << fTotalNofHits / Double_t(fNofTs);
+ LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fTotalTime / Double_t(fNofTs) << " ms";
+ if (fCbmEvents) {
+ LOG(info) << "Events : " << fNofEvents;
+ LOG(info) << "Events / TS : " << fixed << setprecision(2) << fNofEvents / Double_t(fNofTs);
+ if (fNofEvents > 0) {
+ LOG(info) << "Digis / ev : " << fixed << setprecision(2) << fTotalNofDigis / Double_t(fNofEvents);
+ LOG(info) << "Hits / ev : " << fixed << setprecision(2) << fTotalNofHits / Double_t(fNofEvents);
+ LOG(info) << "Time / ev : " << fixed << setprecision(2) << 1000. * fTotalTime / Double_t(fNofEvents)
+ << " ms";
+ }
+ }
+
+ TString eventOrTsStr = fCbmEvents && fNofEvents > 0 ? "ev: " : "TS: ";
+ Double_t eventOrTsValue = Double_t(fCbmEvents && fNofEvents > 0 ? fNofEvents : fNofTs);
+ if (fTotalTime != 0.) {
+ LOG(info) << "===== Time by task (real time) ======";
+ LOG(info) << "HitProducer / " << eventOrTsStr << fixed << setprecision(2) << setw(9) << right
+ << 1000. * fHitProducerTime / eventOrTsValue << " ms [" << setw(5) << right
+ << 100 * fHitProducerTime / fTotalTime << " %]";
+#if HAVE_ONNXRUNTIME
+ if (fUseDenoiseNN)
+ LOG(info) << "Denoise NN / " << eventOrTsStr << fixed << setprecision(2) << setw(9) << right
+ << 1000. * fDenoiseCnnTime / eventOrTsValue << " ms [" << setw(5) << right
+ << 100 * fDenoiseCnnTime / fTotalTime << " %]";
+#endif
+ }
+ LOG(info) << "=====================================\n";
+}
bool CbmRichMCbmHitProducer::isInToT(const double ToT)
diff --git a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
index c51a7a98ac69ef14c38152d8fe1e2c4ceed13a8f..2f843d08e2715167348da1543dcb73ef13f6d1bb 100644
--- a/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
+++ b/reco/detectors/rich/mcbm/CbmRichMCbmHitProducer.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2019-2020 UGiessen/JINR-LIT, Giessen/Dubna
+/* Copyright (C) 2019-2024 UGiessen/JINR-LIT, Giessen/Dubna
SPDX-License-Identifier: GPL-3.0-only
- Authors: Semen Lebedev [committer], Adrian Amatus Weber */
+ Authors: Semen Lebedev [committer], Adrian Amatus Weber, Martin Beyer */
#ifndef CBM_RICH_MCBM_HIT_PRODUCER
#define CBM_RICH_MCBM_HIT_PRODUCER
@@ -14,6 +14,7 @@ class TClonesArray;
class TVector3;
class CbmEvent;
class CbmRichDigi;
+class CbmRichMCbmDenoiseCnn;
class CbmRichMCbmMappingData {
public:
@@ -126,8 +127,19 @@ public:
*/
void DoRestrictToFullAcc(bool val = true) { fRestrictToFullAcc = val; }
+#if HAVE_ONNXRUNTIME
+ void applyDenoiseNN(bool val = true) { fUseDenoiseNN = val; }
+ void SetClassifierThreshold(float val) { fDenoiseCnnThreshold = val; }
+#endif
+
+ private:
+#if HAVE_ONNXRUNTIME
+ bool fUseDenoiseNN = true;
+ float fDenoiseCnnThreshold = 0.5;
+ std::unique_ptr<CbmRichMCbmDenoiseCnn> fDenoiseCnn = nullptr;
+ double fDenoiseCnnTime = 0.;
+#endif
-private:
CbmDigiManager* fDigiMan = nullptr;
TClonesArray* fRichHits; // RICH hits
TClonesArray* fCbmEvents = nullptr; // CbmEvent for time-based simulations
@@ -141,10 +153,18 @@ private:
std::map<Int_t, CbmRichMCbmMappingData> fRichMapping;
- Int_t fEventNum; // event number
+ double fTotalTime = 0.;
+ double fHitProducerTime = 0.;
+ Int_t fNofTs;
+
+ Int_t fNofDigis = 0;
Int_t fNofHits = 0;
+ Int_t fNofEvents = 0;
+ Int_t fTotalNofDigis = 0;
+ Int_t fTotalNofHits = 0;
+
Double_t fHitError;
std::string fMappingFile;