From f025ce83fc2aa691877b767024546cfaed29e098 Mon Sep 17 00:00:00 2001
From: Pierre-Alain Loizeau <p.-a.loizeau@gsi.de>
Date: Thu, 23 Nov 2023 12:49:15 +0000
Subject: [PATCH] DC Nov23 backports: cbmreco: Enable storage of STS clusters
and digis and extend archive plotter.
---
algo/base/Definitions.h | 4 +-
algo/data/sts/Cluster.h | 16 +++
algo/detectors/sts/HitfinderChain.cxx | 46 ++++++-
algo/detectors/sts/HitfinderChain.h | 9 +-
algo/global/Reco.cxx | 16 ++-
algo/global/RecoResults.h | 11 +-
algo/global/StorableRecoResults.h | 11 ++
reco/app/cbmreco/main.cxx | 3 +
services/archive_explorer/CMakeLists.txt | 2 +
services/archive_explorer/app/Application.cxx | 23 +++-
services/archive_explorer/app/Histograms.cxx | 90 +++++++++++++
services/archive_explorer/app/Histograms.h | 46 +++++++
services/archive_explorer/app/Options.cxx | 8 ++
services/archive_explorer/app/Options.h | 30 +++--
.../lib/HistogramCollection.cxx | 54 ++++++++
.../lib/HistogramCollection.h | 78 ++++++++++++
services/archive_explorer/lib/Server.cxx | 119 +++++++-----------
services/archive_explorer/lib/Server.h | 40 +++---
18 files changed, 482 insertions(+), 124 deletions(-)
create mode 100644 services/archive_explorer/app/Histograms.cxx
create mode 100644 services/archive_explorer/app/Histograms.h
create mode 100644 services/archive_explorer/lib/HistogramCollection.cxx
create mode 100644 services/archive_explorer/lib/HistogramCollection.h
diff --git a/algo/base/Definitions.h b/algo/base/Definitions.h
index 490edc62a4..cf2742fcfe 100644
--- a/algo/base/Definitions.h
+++ b/algo/base/Definitions.h
@@ -41,7 +41,7 @@ namespace cbm::algo
enum class RecoData
{
- Digi, //< Raw output from unpackers
+ RawDigi, //< Raw output from unpackers
DigiEvent, //< Digis after event building
Cluster,
Hit,
@@ -81,7 +81,7 @@ CBM_ENUM_DICT(cbm::algo::Step,
);
CBM_ENUM_DICT(cbm::algo::RecoData,
- {"Digi", RecoData::Digi},
+ {"RawDigi", RecoData::RawDigi},
{"DigiEvent", RecoData::DigiEvent},
{"Cluster", RecoData::Cluster},
{"Hit", RecoData::Hit},
diff --git a/algo/data/sts/Cluster.h b/algo/data/sts/Cluster.h
index 1e102b32e4..56b89092f0 100644
--- a/algo/data/sts/Cluster.h
+++ b/algo/data/sts/Cluster.h
@@ -4,6 +4,8 @@
#ifndef CBM_ALGO_DATA_STS_CLUSTER_H
#define CBM_ALGO_DATA_STS_CLUSTER_H
+#include <boost/serialization/access.hpp>
+
#include "Definitions.h"
namespace cbm::algo::sts
@@ -16,6 +18,20 @@ namespace cbm::algo::sts
real fPositionError; ///< Cluster centre error (r.m.s.) in channel number units
u32 fTime; ///< cluster time [ns]
real fTimeError; ///< Error of cluster time [ns]
+
+ private: // serialization
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, unsigned int /*version*/)
+ {
+ ar& fCharge;
+ ar& fSize;
+ ar& fPosition;
+ ar& fPositionError;
+ ar& fTime;
+ ar& fTimeError;
+ }
};
} // namespace cbm::algo::sts
diff --git a/algo/detectors/sts/HitfinderChain.cxx b/algo/detectors/sts/HitfinderChain.cxx
index d0a2d5aa5c..65ade190e7 100644
--- a/algo/detectors/sts/HitfinderChain.cxx
+++ b/algo/detectors/sts/HitfinderChain.cxx
@@ -44,7 +44,7 @@ void sts::HitfinderChain::Finalize()
xpu::set<TheHitfinder>(fHitfinder);
}
-sts::HitfinderChain::Result sts::HitfinderChain::operator()(gsl::span<const CbmStsDigi> digis)
+sts::HitfinderChain::Result sts::HitfinderChain::operator()(gsl::span<const CbmStsDigi> digis, bool storeClusters)
{
EnsureParameters();
@@ -169,7 +169,7 @@ sts::HitfinderChain::Result sts::HitfinderChain::operator()(gsl::span<const CbmS
// xpu::run_kernel<FindClustersBasic>(xpu::grid::n_blocks(hfc.nModules * 2));
// xpu::run_kernel<CalculateClustersBasic>(xpu::grid::n_blocks(hfc.nModules * 2));
L_(debug) << "STS Hitfinder Chain: Sort Clusters...";
- queue.launch<SortClusters>(xpu::n_blocks(nModuleSides)); // FIXME n_blocks(nModuleSides) once debugging is done
+ queue.launch<SortClusters>(xpu::n_blocks(nModuleSides));
L_(debug) << "STS Hitfinder Chain: Find Hits...";
queue.copy(hfc.nClustersPerModule, xpu::d2h);
@@ -231,6 +231,9 @@ sts::HitfinderChain::Result sts::HitfinderChain::operator()(gsl::span<const CbmS
}
}
+ PartitionedVector<sts::Cluster> clusters;
+ if (storeClusters) clusters = FlattenClusters(queue);
+
size_t nHitsTotal = std::accumulate(nHits.begin(), nHits.end(), 0);
L_(info) << "Timeslice contains " << nHitsTotal << " STS hits and " << nClustersTotal << " STS clusters";
@@ -241,6 +244,7 @@ sts::HitfinderChain::Result sts::HitfinderChain::operator()(gsl::span<const CbmS
Result result;
result.hits = std::move(hits);
+ result.clusters = std::move(clusters);
result.monitor = monitor[0];
return result;
}
@@ -578,6 +582,44 @@ PartitionedSpan<sts::Hit> sts::HitfinderChain::FlattenHits(xpu::queue queue)
return partitionedHits;
}
+PartitionedVector<sts::Cluster> sts::HitfinderChain::FlattenClusters(xpu::queue queue)
+{
+ auto& hfc = fHitfinder;
+
+ xpu::h_view nClusters(hfc.nClustersPerModule);
+ int maxClustersPerModule = hfc.maxClustersPerModule;
+ sts::Cluster* clusters = hfc.clusterDataPerModule.get();
+
+ int nModuleSides = nClusters.size();
+
+ L_(debug) << "Copy clusters from " << nModuleSides << " modules sides";
+
+ std::vector<u32> addresses;
+ for (size_t i = 0; i < 2; i++) { // Repeat twice for front and back
+ for (auto& m : fPars->setup.modules)
+ addresses.push_back(m.address);
+ }
+
+ size_t nClustersTotal = std::accumulate(nClusters.begin(), nClusters.end(), 0);
+
+ // TODO: copy could be made much more efficient by using pinned memory
+ // But clusters are stored for debugging only, so no effort is made to optimize this
+ std::vector<sts::Cluster> clustersFlat;
+ clustersFlat.resize(nClustersTotal);
+ size_t offset = 0;
+ for (int m = 0; m < nModuleSides; m++) {
+ size_t nClustersInModule = nClusters[m];
+ queue.copy(clusters + maxClustersPerModule * m, clustersFlat.data() + offset,
+ nClustersInModule * sizeof(sts::Cluster));
+ offset += nClustersInModule;
+ }
+ queue.wait();
+
+ std::vector<size_t> sizes(nClusters.begin(), nClusters.end());
+ PartitionedVector<sts::Cluster> out(std::move(clustersFlat), gsl::make_span(sizes), gsl::make_span(addresses));
+ return out;
+}
+
size_t sts::HitfinderChain::GetNHits(xpu::h_view<int> nHitsPerModule, int module)
{
return std::min<size_t>(nHitsPerModule[module], fHitfinder.maxHitsPerModule);
diff --git a/algo/detectors/sts/HitfinderChain.h b/algo/detectors/sts/HitfinderChain.h
index ddcf955a65..891a3d6187 100644
--- a/algo/detectors/sts/HitfinderChain.h
+++ b/algo/detectors/sts/HitfinderChain.h
@@ -15,6 +15,7 @@
#include <vector>
#include "PartitionedSpan.h"
+#include "PartitionedVector.h"
#include "SubChain.h"
#include "gpu/xpu_legacy.h"
#include "sts/Hitfinder.h"
@@ -40,6 +41,7 @@ namespace cbm::algo::sts
public:
struct Result {
PartitionedSpan<sts::Hit> hits;
+ PartitionedVector<sts::Cluster> clusters;
sts::HitfinderMonitor monitor;
};
@@ -51,7 +53,7 @@ namespace cbm::algo::sts
*/
void Finalize();
- Result operator()(gsl::span<const CbmStsDigi>);
+ Result operator()(gsl::span<const CbmStsDigi>, bool storeClusters = false);
private:
static constexpr u16 InvalidModule = 0xFFFF;
@@ -104,6 +106,11 @@ namespace cbm::algo::sts
*/
PartitionedSpan<sts::Hit> FlattenHits(xpu::queue);
+ /**
+ * Copy Clusters back to host into a flat array.
+ */
+ PartitionedVector<sts::Cluster> FlattenClusters(xpu::queue);
+
/**
* @brief Returns the number of hits of a module.
*
diff --git a/algo/global/Reco.cxx b/algo/global/Reco.cxx
index 56be1fbbc7..41965fe098 100644
--- a/algo/global/Reco.cxx
+++ b/algo/global/Reco.cxx
@@ -27,9 +27,6 @@ void Reco::Validate(const Options& opts)
{
if (!fs::exists(opts.ParamsDir())) throw FatalError("ParamsDir does not exist: {}", opts.ParamsDir().string());
- if (opts.HasOutput(RecoData::Digi)) throw FatalError("Storing raw digis via 'Digi' option not supported");
- if (opts.HasOutput(RecoData::Cluster)) throw FatalError("Storing clusters via 'Cluster' option not supported");
-
bool hasOutputFile = !opts.OutputFile().empty();
bool hasOutputType = !opts.OutputTypes().empty();
@@ -154,9 +151,12 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
sts::HitfinderMonitor stsHitfinderMonitor;
PartitionedSpan<sts::Hit> stsHits;
+ PartitionedVector<sts::Cluster> stsClusters;
if (Opts().HasStep(Step::LocalReco) && Opts().HasDetector(fles::Subsystem::STS)) {
- auto stsResults = fStsHitFinder(unpackResult.first.fSts);
+ bool storeClusters = Opts().HasOutput(RecoData::Cluster);
+ auto stsResults = fStsHitFinder(unpackResult.first.fSts, storeClusters);
stsHits = stsResults.hits;
+ stsClusters = std::move(stsResults.clusters);
stsHitfinderMonitor = stsResults.monitor;
QueueStsRecoMetrics(stsHitfinderMonitor);
}
@@ -185,9 +185,13 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
QueueTrackingMetrics(trackingOutput.monitorData);
}
+ if (Opts().HasOutput(RecoData::RawDigi)) results.stsDigis = std::move(unpackResult.first.fSts);
if (Opts().HasOutput(RecoData::DigiEvent)) results.events = std::move(events);
- if (Opts().HasOutput(RecoData::Hit)) results.stsHits = std::move(stsHits);
- if (Opts().HasOutput(RecoData::Hit)) results.tofHits = std::move(tofHits);
+ if (Opts().HasOutput(RecoData::Cluster)) results.stsClusters = std::move(stsClusters);
+ if (Opts().HasOutput(RecoData::Hit)) {
+ results.stsHits = std::move(stsHits);
+ results.tofHits = std::move(tofHits);
+ }
}
PrintTimings(tsTimes);
diff --git a/algo/global/RecoResults.h b/algo/global/RecoResults.h
index 42f1aa04a0..ecc34dd934 100644
--- a/algo/global/RecoResults.h
+++ b/algo/global/RecoResults.h
@@ -1,6 +1,6 @@
/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergei Zharko [committer] */
+ Authors: Sergei Zharko [committer], Felix Weiglhofer */
/// \file RecoResults.h
/// \date 22.10.2023
@@ -15,14 +15,17 @@
#include "PartitionedSpan.h"
#include "PartitionedVector.h"
#include "ca/core/utils/CaVector.h"
-#include "data/sts/Hit.h"
+#include "sts/Cluster.h"
+#include "sts/Hit.h"
namespace cbm::algo
{
- /// \struct RecoResults
- /// \brief Structure to keep reconstructed results: digi-events, hits and tracks
+ /// @name RecoResults
+ /// @brief Structure to keep reconstructed results: digi-events, hits and tracks
struct RecoResults {
std::vector<DigiEvent> events;
+ PODVector<CbmStsDigi> stsDigis;
+ PartitionedVector<sts::Cluster> stsClusters;
PartitionedSpan<sts::Hit> stsHits;
PartitionedVector<tof::Hit> tofHits;
ca::Vector<ca::Track> tracks;
diff --git a/algo/global/StorableRecoResults.h b/algo/global/StorableRecoResults.h
index 87d99fc8b3..8a11c4522b 100644
--- a/algo/global/StorableRecoResults.h
+++ b/algo/global/StorableRecoResults.h
@@ -17,6 +17,7 @@
#include "PartitionedVector.h"
#include "ca/core/data/CaTrack.h"
#include "ca/core/utils/CaVector.h"
+#include "sts/Cluster.h"
#include "sts/Hit.h"
namespace cbm::algo
@@ -47,6 +48,12 @@ namespace cbm::algo
std::vector<CbmDigiEvent>& DigiEvents() { return fDigiEvents; }
const std::vector<CbmDigiEvent>& DigiEvents() const { return fDigiEvents; }
+ std::vector<CbmStsDigi>& StsDigis() { return fStsDigis; }
+ const std::vector<CbmStsDigi>& StsDigis() const { return fStsDigis; }
+
+ PartitionedVector<sts::Cluster>& StsClusters() { return fStsClusters; }
+ const PartitionedVector<sts::Cluster>& StsClusters() const { return fStsClusters; }
+
PartitionedVector<sts::Hit>& StsHits() { return fStsHits; }
const PartitionedVector<sts::Hit>& StsHits() const { return fStsHits; }
@@ -67,6 +74,8 @@ namespace cbm::algo
uint64_t fTsStartTime = UINT64_MAX;
std::vector<CbmDigiEvent> fDigiEvents;
+ std::vector<CbmStsDigi> fStsDigis;
+ PartitionedVector<sts::Cluster> fStsClusters;
PartitionedVector<sts::Hit> fStsHits;
PartitionedVector<tof::Hit> fTofHits;
@@ -90,6 +99,8 @@ namespace cbm::algo
ar& fTsStartTime;
ar& fDigiEvents;
+ ar& fStsDigis;
+ ar& fStsClusters;
ar& fStsHits;
ar& fTofHits;
ar& fTracks;
diff --git a/reco/app/cbmreco/main.cxx b/reco/app/cbmreco/main.cxx
index f5e2f84d8d..296671a491 100644
--- a/reco/app/cbmreco/main.cxx
+++ b/reco/app/cbmreco/main.cxx
@@ -29,6 +29,9 @@ std::shared_ptr<StorableRecoResults> makeStorableRecoResults(const fles::Timesli
storable->DigiEvents().emplace_back(digiEvent.ToStorable());
}
+ // TODO: some of these copies can be avoided / made into moves
+ storable->StsDigis() = ToStdVector(results.stsDigis);
+ storable->StsClusters() = results.stsClusters;
storable->StsHits() = results.stsHits;
storable->TofHits() = results.tofHits;
storable->Tracks() = results.tracks;
diff --git a/services/archive_explorer/CMakeLists.txt b/services/archive_explorer/CMakeLists.txt
index 7631ccca32..9e48002750 100644
--- a/services/archive_explorer/CMakeLists.txt
+++ b/services/archive_explorer/CMakeLists.txt
@@ -5,6 +5,7 @@ set(APP cbm-archive-explorer)
set(APP_SRCS
app/Application.cxx
+ app/Histograms.cxx
app/Options.cxx
app/main.cxx
)
@@ -30,6 +31,7 @@ set(INCLUDE_DIRECTORIES
)
set(SRCS
+ lib/HistogramCollection.cxx
lib/Server.cxx
)
diff --git a/services/archive_explorer/app/Application.cxx b/services/archive_explorer/app/Application.cxx
index c3cc8e17e2..31d86f4479 100644
--- a/services/archive_explorer/app/Application.cxx
+++ b/services/archive_explorer/app/Application.cxx
@@ -5,6 +5,7 @@
#include <log.hpp>
+#include "Histograms.h"
#include "RecoResultsInputArchive.h"
#include "Server.h"
@@ -18,8 +19,26 @@ Application::Application(const Options& options) : fOpts(options)
L_(info) << "Opening archive " << fOpts.Input();
auto archive = std::make_shared<RecoResultsInputArchive>(fOpts.Input().string());
-
- fServer = std::make_unique<Server>(fOpts.Port(), archive);
+ auto histograms = std::make_shared<Histograms>();
+
+ std::shared_ptr<RecoResultsInputArchive> archive2;
+ std::shared_ptr<Histograms> histograms2;
+ if (not fOpts.Input2().empty()) {
+ L_(info) << "Opening second archive " << fOpts.Input2();
+ archive2 = std::make_shared<RecoResultsInputArchive>(fOpts.Input2().string());
+ histograms2 = std::make_shared<Histograms>();
+ }
+
+ Server::Settings settings {
+ .port = fOpts.Port(),
+ .sensor = fOpts.Sensor(),
+ .archive = archive,
+ .histograms = histograms,
+ .archive2 = archive2,
+ .histograms2 = histograms2,
+ };
+
+ fServer = std::make_unique<Server>(settings);
}
Application::~Application() {}
diff --git a/services/archive_explorer/app/Histograms.cxx b/services/archive_explorer/app/Histograms.cxx
new file mode 100644
index 0000000000..df937e0518
--- /dev/null
+++ b/services/archive_explorer/app/Histograms.cxx
@@ -0,0 +1,90 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#include "Histograms.h"
+
+#include <TH1.h>
+
+#include <log.hpp>
+#include <unordered_map>
+
+#include "StorableRecoResults.h"
+
+using namespace cbm::explore;
+
+Histograms::Histograms()
+{
+ CreateFolder("/sts", "STS");
+ CreateFolder("/sts/digis", "Digis");
+ CreateHisto(fHStsDigisTime, "/sts/digis", "hStsDigisTime", "Sts Digis Time", 1000, 0, 128000000);
+
+ CreateFolder("/sts/clusters", "Clusters");
+ CreateHisto(fHStsClustersTime, "/sts/clusters", "hStsClustersTime", "Sts Clusters Time", 1000, 0, 128000000);
+
+ CreateFolder("/sts/hits", "Hits");
+ CreateHisto(fHStsHitsX, "/sts/hits", "hStsHitsX", "Sts Hits X", 100, -10, 10);
+ CreateHisto(fHStsHitsY, "/sts/hits", "hStsHitsY", "Sts Hits Y", 100, -10, 10);
+ CreateHisto(fHStsHitsZ, "/sts/hits", "hStsHitsZ", "Sts Hits Z", 100, 0, 50);
+ CreateHisto(fHStsHitsTime, "/sts/hits", "hStsHitsTime", "Sts Hits Time", 1000, 0, 128000000);
+ CreateHisto(fHStsHitsDx, "/sts/hits", "hStsHitsDx", "Sts Hits Dx", 100, 0, 0.1);
+ CreateHisto(fHStsHitsDy, "/sts/hits", "hStsHitsDy", "Sts Hits Dy", 100, 0, 0.1);
+ CreateHisto(fHStsHitsDz, "/sts/hits", "hStsHitsDz", "Sts Hits Dz", 100, 0, 0.1);
+ CreateHisto(fHStsHitsDxy, "/sts/hits", "hStsHitsDxy", "Sts Hits Dxy", 100, 0, 0.1);
+ CreateHisto(fHStsHitsTimeError, "/sts/hits", "hStsHitsTimeError", "Sts Hits Time Error", 100, 0, 0.1);
+ CreateHisto(fHStsHitsDu, "/sts/hits", "hStsHitsDu", "Sts Hits Du", 100, 0, 0.1);
+ CreateHisto(fHStsHitsDv, "/sts/hits", "hStsHitsDv", "Sts Hits Dv", 100, 0, 0.1);
+
+ CreateHisto(fHNHitsFromCluster, "/sts/hits", "hNHitsFromCluster", "Number of hits from cluster", 100, 0, 100);
+}
+
+void Histograms::FillHistos(HistoData fill)
+{
+ auto& stsDigis = fill.data->StsDigis();
+ L_(info) << "Filling histos with " << stsDigis.size() << " STS digis";
+ for (auto& digi : stsDigis) {
+ if (fill.Skip(digi.GetAddress())) continue;
+ fHStsDigisTime->Fill(digi.GetTimeU32());
+ }
+
+ auto& stsClusters = fill.data->StsClusters();
+ for (size_t p = 0; p < stsClusters.NPartitions(); p++) {
+ auto [sensor, address] = stsClusters.Partition(p);
+ if (fill.Skip(address)) continue;
+ L_(info) << "Filling STS sensor " << address << " with " << sensor.size() << " clusters";
+ for (auto& cluster : sensor) {
+ fHStsClustersTime->Fill(cluster.fTime);
+ }
+ }
+
+ auto& stsHits = fill.data->StsHits();
+ for (size_t p = 0; p < stsHits.NPartitions(); p++) {
+ auto [sensor, address] = stsHits.Partition(p);
+ if (fill.Skip(address)) continue;
+ L_(info) << "Filling STS sensor " << address << " with " << sensor.size() << " hits";
+
+ std::unordered_map<int, int> nHitsFromClusterF, nHitsFromClusterB;
+
+ for (auto& hit : sensor) {
+ fHStsHitsX->Fill(hit.X());
+ fHStsHitsY->Fill(hit.Y());
+ fHStsHitsZ->Fill(hit.Z());
+ fHStsHitsTime->Fill(hit.Time());
+ fHStsHitsDx->Fill(hit.Dx());
+ fHStsHitsDy->Fill(hit.Dy());
+ fHStsHitsDz->Fill(hit.fDz);
+ fHStsHitsDxy->Fill(hit.fDxy);
+ fHStsHitsTimeError->Fill(hit.TimeError());
+ fHStsHitsDu->Fill(hit.fDu);
+ fHStsHitsDv->Fill(hit.fDv);
+
+ nHitsFromClusterF[hit.fFrontClusterId]++;
+ nHitsFromClusterB[hit.fBackClusterId]++;
+ }
+
+ for (auto& [clusterId, nHits] : nHitsFromClusterF)
+ fHNHitsFromCluster->Fill(nHits);
+
+ for (auto& [clusterId, nHits] : nHitsFromClusterB)
+ fHNHitsFromCluster->Fill(nHits);
+ }
+}
diff --git a/services/archive_explorer/app/Histograms.h b/services/archive_explorer/app/Histograms.h
new file mode 100644
index 0000000000..1a4561fbd8
--- /dev/null
+++ b/services/archive_explorer/app/Histograms.h
@@ -0,0 +1,46 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#pragma once
+
+#include "HistogramCollection.h"
+
+class TH1F;
+class TH1I;
+
+namespace cbm::explore
+{
+
+ class Histograms : public HistogramCollection {
+
+ public:
+ Histograms();
+
+ virtual ~Histograms() {}
+
+ void FillHistos(HistoData fill) override;
+
+ private:
+ // Sts Digis
+ TH1I* fHStsDigisTime = nullptr;
+
+ // Sts Clusters
+ TH1I* fHStsClustersTime = nullptr;
+
+ // Sts Hits
+ TH1F* fHStsHitsX = nullptr;
+ TH1F* fHStsHitsY = nullptr;
+ TH1F* fHStsHitsZ = nullptr;
+ TH1I* fHStsHitsTime = nullptr;
+ TH1F* fHStsHitsDx = nullptr;
+ TH1F* fHStsHitsDy = nullptr;
+ TH1F* fHStsHitsDz = nullptr;
+ TH1F* fHStsHitsDxy = nullptr;
+ TH1F* fHStsHitsTimeError = nullptr;
+ TH1F* fHStsHitsDu = nullptr;
+ TH1F* fHStsHitsDv = nullptr;
+
+ TH1I* fHNHitsFromCluster = nullptr;
+ };
+
+} // namespace cbm::explore
diff --git a/services/archive_explorer/app/Options.cxx b/services/archive_explorer/app/Options.cxx
index 9520cc81bd..50ad358a2c 100644
--- a/services/archive_explorer/app/Options.cxx
+++ b/services/archive_explorer/app/Options.cxx
@@ -9,6 +9,8 @@
namespace po = boost::program_options;
+using namespace cbm::explore;
+
Options::Options(int argc, char** argv)
{
po::options_description options("Options");
@@ -17,8 +19,12 @@ Options::Options(int argc, char** argv)
options.add_options()
("input,i", po::value(&fInput)->value_name("<input>")->required(),
"input archive")
+ ("input2,j", po::value(&fInput2)->value_name("<input2>"),
+ "second input archive, used for comparison (optional)")
("port,p", po::value(&fPort)->value_name("<port>")->default_value(8080),
"port to listen on")
+ ("sensor", po::value(&fSensor)->value_name("<sensor>"),
+ "sensor to filter on")
("help,h",
"produce help message")
;
@@ -42,6 +48,8 @@ Options::Options(int argc, char** argv)
std::exit(EXIT_SUCCESS);
}
+ fFilterSensor = vm.count("sensor") > 0;
+
try {
po::notify(vm);
}
diff --git a/services/archive_explorer/app/Options.h b/services/archive_explorer/app/Options.h
index 9ec407b889..d8a5f0734c 100644
--- a/services/archive_explorer/app/Options.h
+++ b/services/archive_explorer/app/Options.h
@@ -3,18 +3,28 @@
Authors: Felix Weiglhofer [committer] */
#pragma once
-#include <compat/Filesystem.h>
+#include <optional>
-class Options {
-public:
- Options(int argc, char** argv);
+#include "compat/Filesystem.h"
- int Port() const { return fPort; }
+namespace cbm::explore
+{
+ class Options {
+ public:
+ Options(int argc, char** argv);
- cbm::algo::fs::path Input() const { return fInput; }
+ int Port() const { return fPort; }
+ algo::fs::path Input() const { return fInput; }
+ algo::fs::path Input2() const { return fInput2; }
-private:
- int fPort;
- std::string fInput;
-};
+ std::optional<uint32_t> Sensor() const { return fFilterSensor ? std::make_optional(fSensor) : std::nullopt; }
+
+ private:
+ int fPort;
+ uint32_t fSensor;
+ bool fFilterSensor = false;
+ std::string fInput;
+ std::string fInput2;
+ };
+} // namespace cbm::explore
diff --git a/services/archive_explorer/lib/HistogramCollection.cxx b/services/archive_explorer/lib/HistogramCollection.cxx
new file mode 100644
index 0000000000..527f51de88
--- /dev/null
+++ b/services/archive_explorer/lib/HistogramCollection.cxx
@@ -0,0 +1,54 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#include "HistogramCollection.h"
+
+#include <TH1.h>
+
+#include <fmt/format.h>
+
+using namespace cbm::explore;
+
+HistogramCollection::~HistogramCollection() {}
+
+void HistogramCollection::Reset()
+{
+ for (auto& histo : fHistos) {
+ histo.histo->Reset();
+ }
+}
+
+TH1* HistogramCollection::GetHisto(const std::string& name) const
+{
+ for (auto& histo : fHistos) {
+ if (histo.histo->GetName() == name) { return histo.histo; }
+ }
+ return nullptr;
+}
+
+std::string HistogramCollection::GetHistoPath(const std::string& name) const
+{
+ for (auto& histo : fHistos) {
+ if (histo.histo->GetName() == name) { return histo.path; }
+ }
+ return "";
+}
+
+void HistogramCollection::Div(const HistogramCollection& other)
+{
+ for (auto& histo : fHistos) {
+ auto otherHisto = other.GetHisto(histo.histo->GetName());
+ if (not otherHisto) {
+ throw std::runtime_error(
+ fmt::format("Cannot divide histograms, because they do not match: No histogram with name "
+ "{} found in other collection.",
+ histo.histo->GetName()));
+ }
+ histo.histo->Divide(otherHisto);
+ }
+}
+
+void HistogramCollection::CreateFolder(const char* path, const char* name)
+{
+ fFolders.push_back({.path = path, .name = name});
+}
diff --git a/services/archive_explorer/lib/HistogramCollection.h b/services/archive_explorer/lib/HistogramCollection.h
new file mode 100644
index 0000000000..8a5760c74d
--- /dev/null
+++ b/services/archive_explorer/lib/HistogramCollection.h
@@ -0,0 +1,78 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#pragma once
+
+#include <memory>
+#include <optional>
+#include <string>
+#include <vector>
+
+class TH1;
+
+namespace cbm::algo
+{
+ class StorableRecoResults;
+}
+
+namespace cbm::explore
+{
+
+ struct Folder {
+ std::string path;
+ std::string name;
+ };
+
+ struct Histo {
+ std::string path;
+ TH1* histo;
+ };
+
+ struct HistoData {
+ algo::StorableRecoResults* data = nullptr;
+ std::optional<uint32_t> sensor;
+
+ bool Skip(uint32_t address) const
+ {
+ if (not sensor) return false;
+ return address != *sensor;
+ }
+ };
+
+ class HistogramCollection {
+
+ public:
+ HistogramCollection() = default;
+ virtual ~HistogramCollection();
+
+ /// @brief Clear all histograms.
+ void Reset();
+
+ TH1* GetHisto(const std::string& name) const;
+ std::string GetHistoPath(const std::string& name) const;
+ std::vector<Histo> GetHistos() const { return fHistos; }
+
+ std::vector<Folder> GetFolders() const { return fFolders; }
+
+ void Div(const HistogramCollection& other);
+
+ virtual void FillHistos(HistoData fill) = 0;
+
+ protected:
+ template<typename Histo_t>
+ void CreateHisto(Histo_t*& histo, const char* folder, const char* name, const char* title, int nbins, double xlow,
+ double xmax)
+ {
+ histo = new Histo_t(name, title, nbins, xlow, xmax);
+ histo->SetDirectory(nullptr);
+ fHistos.push_back({.path = folder, .histo = histo});
+ }
+
+ void CreateFolder(const char* path, const char* name);
+
+ private:
+ std::vector<Folder> fFolders;
+ std::vector<Histo> fHistos;
+ };
+
+} // namespace cbm::explore
diff --git a/services/archive_explorer/lib/Server.cxx b/services/archive_explorer/lib/Server.cxx
index 331e4894db..cacdad97f3 100644
--- a/services/archive_explorer/lib/Server.cxx
+++ b/services/archive_explorer/lib/Server.cxx
@@ -15,23 +15,21 @@ using namespace cbm::explore;
ClassImp(Server);
-Server::Server(int port, std::shared_ptr<algo::RecoResultsInputArchive> archive)
- : TNamed("server", "server")
- , fArchive(archive)
+Server::Server(Settings settings) : TNamed("server", "server"), fS(settings)
{
- L_(info) << "Starting Server with port " << port;
+ L_(info) << "Starting Server with port " << fS.port;
- fServer = new THttpServer(fmt::format("http:{}", port).c_str());
+ if (fS.sensor) L_(info) << "Filtering on sensor " << *fS.sensor;
+ else {
+ L_(info) << "Not filtering on sensor";
+ }
+
+ fServer = new THttpServer(fmt::format("http:{}", fS.port).c_str());
fServer->SetTimer(0, true);
// FIXME: Why doesn't root use the local jsroot by default???
fServer->SetJSROOT("https://jsroot.gsi.de/7.1.0/");
-}
-
-Server::~Server() {}
-int Server::Run()
-{
// register this object in http server
fServer->Register("/", this);
fServer->Hide("/server");
@@ -46,7 +44,30 @@ int Server::Run()
fServer->RegisterCommand("/NextTS", "/server/->RequestNextTS()");
fServer->SetItemField("/NextTS", "_title", "Next Timeslice");
- CreateHistos();
+ if (fS.archive == nullptr) throw std::runtime_error("Archive must not be null");
+ if (fS.histograms == nullptr) throw std::runtime_error("Histograms must not be null");
+
+ if (fS.archive2 != nullptr || fS.histograms2 != nullptr) {
+ if (fS.archive2 == nullptr) throw std::runtime_error("Archive2 must not be null");
+ if (fS.histograms2 == nullptr) throw std::runtime_error("Histograms2 must not be null");
+ }
+
+ // Setup folders
+ for (auto& folder : fS.histograms->GetFolders()) {
+ fServer->CreateItem(folder.path.c_str(), folder.name.c_str());
+ fServer->SetItemField(folder.path.c_str(), "_kind", "Folder");
+ }
+
+ // Setup histograms
+ for (auto& [path, histo] : fS.histograms->GetHistos()) {
+ fServer->Register(path.c_str(), histo);
+ }
+}
+
+Server::~Server() {}
+
+int Server::Run()
+{
NextTS();
L_(info) << "Starting event loop";
@@ -77,73 +98,27 @@ void Server::NextTS()
{
L_(info) << "Fetching next Timeslice...";
- auto recoResults = fArchive->get();
-
- if (fArchive->eos()) {
+ auto recoResults = fS.archive->get();
+ if (fS.archive->eos()) {
L_(info) << "End of archive reached";
return;
}
- ResetHistos();
-
- auto& stsHits = recoResults->StsHits();
- for (size_t p = 0; p < stsHits.NPartitions(); p++) {
- auto [sensor, address] = stsHits.Partition(p);
- L_(info) << "Filling sensor " << address << " with " << sensor.size() << " hits";
- for (auto& hit : sensor) {
- fHStsHitsX->Fill(hit.X());
- fHStsHitsY->Fill(hit.Y());
- fHStsHitsZ->Fill(hit.Z());
- fHStsHitsTime->Fill(hit.Time());
- fHStsHitsDx->Fill(hit.Dx());
- fHStsHitsDy->Fill(hit.Dy());
- fHStsHitsDz->Fill(hit.fDz);
- fHStsHitsDxy->Fill(hit.fDxy);
- fHStsHitsTimeError->Fill(hit.TimeError());
- fHStsHitsDu->Fill(hit.fDu);
- fHStsHitsDv->Fill(hit.fDv);
- }
- }
+ HistoData fill {.data = recoResults.get(), .sensor = fS.sensor};
- L_(info) << "Finished Timeslice";
-}
+ fS.histograms->Reset();
+ fS.histograms->FillHistos(fill);
-void Server::ResetHistos()
-{
- for (auto* histo : fHistos) {
- histo->Reset();
+ bool divide = fS.archive2 != nullptr;
+ if (divide) {
+ recoResults = fS.archive2->get();
+ if (fS.archive2->eos()) throw std::runtime_error("End of archive2 reached too early");
+ fS.histograms2->Reset();
+ fill = {.data = recoResults.get(), .sensor = fS.sensor};
+ fS.histograms2->FillHistos(fill);
+ L_(info) << "Dividing histograms";
+ fS.histograms->Div(*fS.histograms2);
}
-}
-void Server::CreateHistos()
-{
- CreateFolder("/sts", "STS");
- CreateFolder("/sts/hits", "Hits");
- CreateHisto(fHStsHitsX, "/sts/hits", "hStsHitsX", "Sts Hits X", 100, -10, 10);
- CreateHisto(fHStsHitsY, "/sts/hits", "hStsHitsY", "Sts Hits Y", 100, -10, 10);
- CreateHisto(fHStsHitsZ, "/sts/hits", "hStsHitsZ", "Sts Hits Z", 100, 0, 50);
- CreateHisto(fHStsHitsTime, "/sts/hits", "hStsHitsTime", "Sts Hits Time", 1000, 0, 128000000);
- CreateHisto(fHStsHitsDx, "/sts/hits", "hStsHitsDx", "Sts Hits Dx", 100, 0, 0.1);
- CreateHisto(fHStsHitsDy, "/sts/hits", "hStsHitsDy", "Sts Hits Dy", 100, 0, 0.1);
- CreateHisto(fHStsHitsDz, "/sts/hits", "hStsHitsDz", "Sts Hits Dz", 100, 0, 0.1);
- CreateHisto(fHStsHitsDxy, "/sts/hits", "hStsHitsDxy", "Sts Hits Dxy", 100, 0, 0.1);
- CreateHisto(fHStsHitsTimeError, "/sts/hits", "hStsHitsTimeError", "Sts Hits Time Error", 100, 0, 0.1);
- CreateHisto(fHStsHitsDu, "/sts/hits", "hStsHitsDu", "Sts Hits Du", 100, 0, 0.1);
- CreateHisto(fHStsHitsDv, "/sts/hits", "hStsHitsDv", "Sts Hits Dv", 100, 0, 0.1);
-}
-
-template<typename Histo_t>
-void Server::CreateHisto(Histo_t*& histo, const char* folder, const char* name, const char* title, int nbins,
- double xlow, double xmax)
-{
- histo = new Histo_t(name, title, nbins, xlow, xmax);
- histo->SetDirectory(nullptr);
- fServer->Register(folder, histo);
- fHistos.push_back(histo);
-}
-
-void Server::CreateFolder(const char* path, const char* name)
-{
- fServer->CreateItem(path, name);
- fServer->SetItemField(path, "_kind", "Folder");
+ L_(info) << "Finished Timeslice";
}
diff --git a/services/archive_explorer/lib/Server.h b/services/archive_explorer/lib/Server.h
index d20c7e3e8a..03f0cf2f5c 100644
--- a/services/archive_explorer/lib/Server.h
+++ b/services/archive_explorer/lib/Server.h
@@ -8,6 +8,7 @@
#include <memory>
+#include "HistogramCollection.h"
#include "RecoResultsInputArchive.h"
class THttpServer;
@@ -21,7 +22,19 @@ namespace cbm::explore
class Server : public TNamed {
public:
- Server(int port, std::shared_ptr<algo::RecoResultsInputArchive> archive);
+ struct Settings {
+ int port;
+ std::optional<uint32_t> sensor; // optional sensor to filter on
+ std::shared_ptr<algo::RecoResultsInputArchive> archive;
+ std::shared_ptr<HistogramCollection> histograms;
+
+ // Optional second archive and histograms for comparison
+ // histograms will be divided by this one, if provided
+ std::shared_ptr<algo::RecoResultsInputArchive> archive2;
+ std::shared_ptr<HistogramCollection> histograms2;
+ };
+
+ Server(Settings settings);
virtual ~Server();
@@ -35,35 +48,12 @@ namespace cbm::explore
// Internal state
THttpServer* fServer = nullptr; //!
- std::shared_ptr<algo::RecoResultsInputArchive> fArchive; //!
- std::vector<TH1*> fHistos; //!
+ Settings fS; //!
// Server commands
bool fRequestNextTS = false;
- // Histograms
- TH1F* fHStsHitsX = nullptr;
- TH1F* fHStsHitsY = nullptr;
- TH1F* fHStsHitsZ = nullptr;
- TH1I* fHStsHitsTime = nullptr;
- TH1F* fHStsHitsDx = nullptr;
- TH1F* fHStsHitsDy = nullptr;
- TH1F* fHStsHitsDz = nullptr;
- TH1F* fHStsHitsDxy = nullptr;
- TH1F* fHStsHitsTimeError = nullptr;
- TH1F* fHStsHitsDu = nullptr;
- TH1F* fHStsHitsDv = nullptr;
-
void NextTS();
- void ResetHistos();
-
- void CreateHistos();
-
- template<typename Histo_t>
- void CreateHisto(Histo_t*& histo, const char* folder, const char* name, const char* title, int nbins, double xlow,
- double xmax);
-
- void CreateFolder(const char* path, const char* name);
ClassDef(Server, 1);
};
--
GitLab