Skip to content
Snippets Groups Projects
Commit 3068b333 authored by Volker Friese's avatar Volker Friese Committed by Felix Weiglhofer
Browse files

Simple QA of DigiEvents implemented in algo, with offline integration into cbmreco_fairrun.

parent b25a81d5
Branches
Tags
1 merge request!1199Simple QA of DigiEvents implemented in algo, with offline integration into cbmreco_fairrun.
Pipeline #23215 failed
......@@ -38,6 +38,8 @@ set(SRCS
detectors/rich/RichReadoutConfig.cxx
detectors/rich/UnpackRich.cxx
global/Reco.cxx
qa/DigiEventQa.cxx
qa/Histo1D.cxx
)
set(BUILD_INFO_CXX ${CMAKE_CURRENT_BINARY_DIR}/base/BuildInfo.cxx)
......@@ -68,6 +70,7 @@ target_include_directories(Algo
${CMAKE_CURRENT_SOURCE_DIR}/evselector
${CMAKE_CURRENT_SOURCE_DIR}/unpack
${CMAKE_CURRENT_SOURCE_DIR}/detectors
${CMAKE_CURRENT_SOURCE_DIR}/qa
)
target_link_libraries(Algo
......
/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Volker Friese [committer] */
#include "DigiEventQa.h"
using std::string;
using std::vector;
namespace cbm::algo::qa
{
// --- Execution --------------------------------------------------------
DigiEventQaData DigiEventQa::operator()(const vector<CbmDigiEvent>& events, const DigiEventQaConfig& config) const
{
// --- Instantiate return object
DigiEventQaData result;
for (const auto& entry : config.fData) {
ECbmModuleId subsystem = entry.first;
string hist_name = "digi_time_" + ToString(subsystem);
auto& detConfig = config.fData.at(subsystem);
result.fDigiTimeHistos.try_emplace(subsystem, detConfig.fNumBins, detConfig.fMinValue, detConfig.fMaxValue,
hist_name);
}
// --- Event loop. Fill histograms.
for (const auto& event : events)
for (auto& subsystem : config.fData)
QaDigiTimeInEvent(event, subsystem.first, result.fDigiTimeHistos.at(subsystem.first));
result.fNumEvents = events.size();
return result;
}
// --------------------------------------------------------------------------
// --- QA: digi time within event ----------------------------------------
void DigiEventQa::QaDigiTimeInEvent(const CbmDigiEvent& event, ECbmModuleId system, Histo1D& histo) const
{
switch (system) {
case ECbmModuleId::kT0: FillDeltaT<CbmBmonDigi>(event.fData.fT0.fDigis, event.fTime, histo); break;
case ECbmModuleId::kSts: FillDeltaT<CbmStsDigi>(event.fData.fSts.fDigis, event.fTime, histo); break;
case ECbmModuleId::kMuch: FillDeltaT<CbmMuchDigi>(event.fData.fMuch.fDigis, event.fTime, histo); break;
case ECbmModuleId::kRich: FillDeltaT<CbmRichDigi>(event.fData.fRich.fDigis, event.fTime, histo); break;
case ECbmModuleId::kTrd: FillDeltaT<CbmTrdDigi>(event.fData.fTrd.fDigis, event.fTime, histo); break;
case ECbmModuleId::kTrd2d: FillDeltaT<CbmTrdDigi>(event.fData.fTrd2d.fDigis, event.fTime, histo); break;
case ECbmModuleId::kTof: FillDeltaT<CbmTofDigi>(event.fData.fTof.fDigis, event.fTime, histo); break;
case ECbmModuleId::kPsd: FillDeltaT<CbmPsdDigi>(event.fData.fPsd.fDigis, event.fTime, histo); break;
default: throw std::runtime_error("DigiEventQa: Invalid system Id " + ToString(system));
}
}
// --------------------------------------------------------------------------
} /* namespace cbm::algo::qa */
/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Volker Friese [committer] */
#ifndef ALGO_QA_DIGIEVENTQA_H
#define ALGO_QA_DIGIEVENTQA_H 1
#include "CbmDefs.h"
#include "CbmDigiEvent.h"
#include <map>
#include <vector>
#include "Histo1D.h"
namespace cbm::algo::qa
{
/** @struct DigiEventQaData
** @brief QA results for CbmDigiEvent objects
** @author Volker Friese <v.friese@gsi.de>
** @since 16 June 2023
**/
struct DigiEventQaData {
std::map<ECbmModuleId, Histo1D> fDigiTimeHistos = {};
size_t fNumEvents = 0;
};
/** @struct DigiEventQaDetConfig
** @brief Configuration data for the QA of CbmDigiEvents for a given detector
** @author Volker Friese <v.friese@gsi.de>
** @since 16 June 2023
**/
struct DigiEventQaDetConfig {
uint32_t fNumBins;
double fMinValue;
double fMaxValue;
std::string ToString() const
{
std::stringstream ss;
ss << "nbins " << fNumBins << " min " << fMinValue << " max " << fMaxValue << "\n";
return ss.str();
}
};
/** @struct DigiEventQaConfig
** @brief Configuration data for the QA of CbmDigiEvents
** @author Volker Friese <v.friese@gsi.de>
** @since 16 June 2023
**/
struct DigiEventQaConfig {
std::map<ECbmModuleId, DigiEventQaDetConfig> fData;
std::string ToString() const
{
std::stringstream ss;
for (const auto& entry : fData)
ss << "\n Subsystem " << ::ToString(entry.first) << " " << entry.second.ToString();
return ss.str();
}
};
/** @class DigiEventQa
** @brief QA for CbmDigiEvent objects
** @author Volker Friese <v.friese@gsi.de>
** @since 16 June 2023
**/
class DigiEventQa {
public:
/** @brief Constructor **/
DigiEventQa() = default;
/** @brief Destructor **/
virtual ~DigiEventQa() = default;
/** @brief Execution
** @param events Vector of DigiEvents to analyse
** @return QA data object
**/
DigiEventQaData operator()(const std::vector<CbmDigiEvent>& events, const DigiEventQaConfig& config) const;
private:
/** @brief Fill histogram with digi time within event
** @param digis Vector with digi objects
** @param eventTime Time of event
** @param histo Histogram to be filled
**
** The templated class is required to implement the method double GetTime().
**/
template<class Digi>
void FillDeltaT(const std::vector<Digi>& digis, double eventTime, Histo1D& histo) const
{
for (const Digi& digi : digis)
histo.Add(digi.GetTime() - eventTime);
}
/** @brief Fill histogram with digi time within event
** @param digis Vector with digi objects
** @param eventTime Time of event
** @param histo Histogram to be filled
**/
void QaDigiTimeInEvent(const CbmDigiEvent& event, ECbmModuleId system, Histo1D& histo) const;
};
} /* namespace cbm::algo::qa */
#endif /* ALGO_QA_DIGIEVENTQA_H */
/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Volker Friese [committer] */
#include "Histo1D.h"
#include <cassert>
#include <cstdint>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
#include <cmath>
namespace cbm::algo
{
// ----- Standard constructor
Histo1D::Histo1D(uint32_t numBins, double minValue, double maxValue, const std::string& name)
: fName(name)
, fNumBins(numBins)
, fMinValue(minValue)
, fMaxValue(maxValue)
{
if (numBins == 0 || maxValue <= minValue)
throw std::runtime_error("Histo1D: Invalid specifications " + std::to_string(numBins) + " "
+ std::to_string(minValue) + " " + std::to_string(maxValue));
fContent.resize(numBins);
}
// ----- Add entry
void Histo1D::Add(double value, double weight)
{
if (value < fMinValue) fUnderflow += weight;
else if (!(value < fMaxValue))
fOverflow += weight;
else {
uint32_t bin = uint32_t(double(fNumBins) * (value - fMinValue) / (fMaxValue - fMinValue));
assert(bin < fNumBins);
fContent[bin] += weight;
fNumEntries++;
}
}
// ----- Clear content
void Histo1D::Clear()
{
fContent.assign(fNumBins, 0.);
fUnderflow = 0;
fOverflow = 0;
fNumEntries = 0;
}
// ----- Content access
double Histo1D::Content(uint32_t bin) const
{
if (bin >= fNumBins) return 0.;
else
return fContent[bin];
}
// ----- First moment of distribution
double Histo1D::Mean() const
{
double sum1 = 0.;
double sum2 = 0.;
double binsize = (fMaxValue - fMinValue) / double(fNumBins);
for (uint32_t bin = 0; bin < fNumBins; bin++) {
double x = fMinValue + (bin + 0.5) * binsize;
sum1 += x * fContent[bin];
sum2 += fContent[bin];
}
return sum1 / sum2;
}
// ----- Operator +=
// TODO: Comparison of floating point numbers; probably fishy.
Histo1D& Histo1D::operator+=(const Histo1D& other)
{
if (other.fNumBins == fNumBins && other.fMinValue == fMinValue && other.fMaxValue == fMaxValue) {
fUnderflow += other.fUnderflow;
fOverflow += other.fOverflow;
for (uint32_t bin = 0; bin < fNumBins; bin++)
fContent[bin] += other.fContent[bin];
fNumEntries += other.fNumEntries;
}
else
throw std::runtime_error("Histo1D: Trying to add incompatible histograms");
return *this;
}
// ----- Second moment of distribution
double Histo1D::Stddev() const
{
double sum1 = 0.;
double sum2 = 0.;
double binsize = (fMaxValue - fMinValue) / double(fNumBins);
for (uint32_t bin = 0; bin < fNumBins; bin++) {
double x = fMinValue + (bin + 0.5) * binsize;
sum1 += x * x * fContent[bin];
sum2 += fContent[bin];
}
double mean = Mean();
return sqrt(sum1 / sum2 - mean * mean);
}
// ----- Properties to string
std::string Histo1D::ToString() const
{
std::stringstream ss;
ss << fName << ": bins " << fNumBins << " range " << fMinValue << " to " << fMaxValue << " entries " << fNumEntries
<< " mean " << Mean() << " stddev " << Stddev() << " out of range " << fUnderflow << " , " << fOverflow;
return ss.str();
}
} /* namespace cbm::algo */
/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Volker Friese [committer] */
#ifndef ALGO_QA_HISTO1D_H
#define ALGO_QA_HISTO1D_H 1
#include <cstdint>
#include <string>
#include <vector>
namespace cbm::algo
{
/** @class Histo1D
** @author Volker Friese <v.friese@gsi.de>
** @since 15 June 2023
**
** Lightweight one-dimensional histogram class
**/
class Histo1D {
public:
/** @brief No default constructor **/
Histo1D() = delete;
/** @brief Standard constructor
** @param numBins Number of bins
** @param minValue Lower edge of histogram range
** @param maxValue Upper edge of histogram range
** @param name Histogram name
**/
Histo1D(uint32_t numBins, double minValue, double maxValue, const std::string& name = "");
/** @brief Destructor **/
virtual ~Histo1D() = default;
/** @brief Add an entry to the histogram
** @param value Value to be histogrammed
** @param weight Weight for the entry
**/
void Add(double value, double weight = 1.);
/** @brief Clear histogram contents **/
void Clear();
/** @brief Histogram content in a bin
** @param bin Bin index
** @return Histogram content. Zero if bin index is out of scope.
**/
double Content(uint32_t bin) const;
/** @brief Number of entries **/
double NumEntries() const { return fNumEntries; }
/** @brief Upper edge **/
double MaxValue() const { return fMaxValue; }
/** @brief First moment of distribution **/
double Mean() const;
/** @brief Lower edge **/
double MinValue() const { return fMinValue; }
/** @brief Histogram name **/
const std::string& Name() const { return fName; }
/** @brief Number of bins **/
uint32_t NumBins() const { return fNumBins; }
/** @brief Add another histogram to an existing one
** @param other Histogram object to be added
**/
Histo1D& operator+=(const Histo1D& other);
/** @brief Overflow **/
double Overflow() const { return fOverflow; }
/** @brief Second moment of distribution **/
double Stddev() const;
/** @brief Properties to string **/
std::string ToString() const;
/** @brief Underflow **/
double UnderFlow() const { return fUnderflow; }
private:
/** Properties **/
const std::string fName;
const uint32_t fNumBins;
const double fMinValue;
const double fMaxValue;
/** Content **/
std::vector<double> fContent = {};
double fUnderflow = 0;
double fOverflow = 0;
size_t fNumEntries = 0;
};
/** @brief Adding two histograms
** @param h1,h2 Histograms to be added
** @return Sum histogram
**/
inline Histo1D operator+(Histo1D h1, const Histo1D& h2)
{
h1 += h2;
return h1;
}
} /* namespace cbm::algo */
#endif /* ALGO_QA_HISTO1D_H */
......@@ -3,6 +3,7 @@
set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}
${CBMROOT_SOURCE_DIR}
)
set(SRCS
......
......@@ -4,23 +4,28 @@
#include "CbmTaskDigiEventQa.h"
#include "CbmDigiBranchBase.h"
#include "CbmDigiManager.h"
#include "CbmDigiTimeslice.h"
#include "CbmModuleList.h"
#include "CbmReco.h" // for CbmRecoConfig
#include <FairRunOnline.h>
#include <Logger.h>
#include <TH1F.h>
#include <TH1D.h>
#include <THttpServer.h>
#include <TStopwatch.h>
#include <cassert>
#include <iomanip>
#include "algo/qa/Histo1D.h"
using namespace std;
using cbm::algo::qa::DigiEventQaConfig;
using cbm::algo::qa::DigiEventQaData;
using cbm::algo::qa::DigiEventQaDetConfig;
#define NUM_BINS 100
#define BORDER 10.
// ----- Constructor -----------------------------------------------------
CbmTaskDigiEventQa::CbmTaskDigiEventQa() : FairTask("DigiEventQa") {}
......@@ -35,11 +40,24 @@ CbmTaskDigiEventQa::~CbmTaskDigiEventQa() {}
// ----- Configuration ---------------------------------------------------
void CbmTaskDigiEventQa::Config(const CbmRecoConfig& config)
{
auto limitIt = config.fEvtbuildWindows.find(ECbmModuleId::kSts);
if (limitIt != config.fEvtbuildWindows.end()) {
double lower = limitIt->second.first - 10.;
double upper = limitIt->second.second + 10.;
fHistDigiTimeSts = new TH1F("hDigiTimeSts", "STS digi time in event", 100, lower, upper);
// The histogram ranges are defined by the event building windows. The number of bins
// is hard-coded. To be changed on request.
for (const auto& entry : config.fEvtbuildWindows) {
ECbmModuleId system = entry.first;
// --- Create histogram
string name = "hDigiTime" + ToString(entry.first);
string title = ToString(entry.first) + " digi time in event";
double lower = entry.second.first - BORDER; // Lower edge of histogram
double upper = entry.second.second + BORDER; // Upper edge of histogram
assert(fDigiTimeHistos.count(system) == 0);
fDigiTimeHistos[system] = new TH1D(name.c_str(), title.c_str(), NUM_BINS, lower, upper);
// --- Set algo configuration
assert(fConfig.fData.count(system) == 0);
fConfig.fData[system] = {NUM_BINS, lower, upper};
}
}
// ---------------------------------------------------------------------------
......@@ -52,55 +70,32 @@ void CbmTaskDigiEventQa::Exec(Option_t*)
// --- Timer and counters
TStopwatch timer;
timer.Start();
double sumT = 0.;
double sumTsq = 0.;
double sumNumDigis = 0.;
double sumNumDigisSq = 0.;
size_t numEvents = 0;
size_t numDigis = 0;
// --- Event loop
for (const auto& event : *fEvents) {
numEvents++;
double numDigisEvt = double(event.fData.fSts.fDigis.size());
sumNumDigis += numDigisEvt;
sumNumDigisSq += numDigisEvt * numDigisEvt;
for (const auto& digi : event.fData.fSts.fDigis) {
numDigis++;
const double tDigi = digi.GetTime() - event.fTime;
if (fHistDigiTimeSts) fHistDigiTimeSts->Fill(tDigi);
sumT += tDigi;
sumTsq += tDigi * tDigi;
}
}
// --- First and second moments of digi times and digi numbers
double tMean = sumT / double(numDigis);
double tSqMean = sumTsq / double(numDigis);
double tRms = sqrt(tSqMean - tMean * tMean);
double numDigisMean = sumNumDigis / double(numEvents);
double numDigisSqMean = sumNumDigisSq / double(numEvents);
double numDigisRms = sqrt(numDigisSqMean - numDigisMean * numDigisMean);
// --- Algo execution
DigiEventQaConfig config;
DigiEventQaData result = fAlgo(*fEvents, fConfig);
// --- Copy QA results (Histo1D) into ROOT output histograms
// TODO: Probably not the most efficient implementation. Creates first a ROOT histogram from a CBM one (data copy),
// which is then added to the member histogram (another data copy). Should implement a method for direct addition TH1D + Histo1D.
for (const auto& entry : result.fDigiTimeHistos) {
ECbmModuleId subsystem = entry.first;
fDigiTimeHistos[subsystem]->Add(ToTH1D(entry.second));
}
// --- Timeslice log
timer.Stop();
fExecTime += timer.RealTime();
size_t numEvents = result.fNumEvents;
stringstream logOut;
logOut << setw(15) << left << GetName() << " [";
logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] ";
logOut << "TS " << fNumTs << ", events " << numEvents << ", digis " << numDigis << ", digis per event ("
<< numDigisMean << " +- " << numDigisRms << ")"
<< ", digi time (" << tMean << " +- " << tRms << ") ns";
logOut << "TS " << fNumTs << ", events " << numEvents;
LOG(info) << logOut.str();
// --- Run statistics
fNumTs++;
fNumEvents += numEvents;
fNumDigis += numDigis;
// See: macro/run/run_unpack_online.C
// See : recp/detectors/sts/unpack/CbmStsUnpackMonitor
}
// ----------------------------------------------------------------------------
......@@ -112,13 +107,15 @@ void CbmTaskDigiEventQa::Finish()
LOG(info) << GetName() << ": Run summary";
LOG(info) << "Timeslices : " << fNumTs;
LOG(info) << "Events : " << fNumEvents;
LOG(info) << "Digis : " << fNumDigis;
LOG(info) << "Exec time : " << fixed << setprecision(2) << 1000. * fExecTime / double(fNumTs) << " ms / TS";
if (fHistDigiTimeSts)
LOG(info) << "STS digi times : entries " << fHistDigiTimeSts->GetEntries() << ", mean "
<< fHistDigiTimeSts->GetMean() << ", stddev " << fHistDigiTimeSts->GetStdDev();
for (const auto& entry : fDigiTimeHistos) {
ECbmModuleId subsystem = entry.first;
TH1D* histo = entry.second;
LOG(info) << ToString(subsystem) << " digi times: entries " << histo->GetEntries() << ", mean " << histo->GetMean()
<< ", stddev " << histo->GetStdDev();
histo->Write();
}
LOG(info) << "=====================================";
if (fHistDigiTimeSts) fHistDigiTimeSts->Write();
}
// ----------------------------------------------------------------------------
......@@ -145,11 +142,14 @@ InitStatus CbmTaskDigiEventQa::Init()
THttpServer* server = FairRunOnline::Instance()->GetHttpServer();
if (server) {
LOG(info) << "--- Http server present; registering histograms";
if (fHistDigiTimeSts) server->Register("DigiEvent", fHistDigiTimeSts);
for (const auto& entry : fDigiTimeHistos)
server->Register("DigiEvent", entry.second);
}
else
LOG(info) << "--- No Http server present";
// --- Log configuration
LOG(info) << "--- Algo configuration for " << fConfig.fData.size() << " subsystem(s): " << fConfig.ToString();
LOG(info) << "==================================================";
return kSUCCESS;
......@@ -157,4 +157,18 @@ InitStatus CbmTaskDigiEventQa::Init()
// ----------------------------------------------------------------------------
// ----- Convert CBM histogram to ROOT histogram --------------------------
TH1D* CbmTaskDigiEventQa::ToTH1D(const cbm::algo::Histo1D& source)
{
TH1D* result =
new TH1D(source.Name().c_str(), source.Name().c_str(), source.NumBins(), source.MinValue(), source.MaxValue());
for (uint32_t bin = 0; bin < source.NumBins(); bin++) {
result->SetBinContent(bin, source.Content(bin));
}
result->SetEntries(source.NumEntries());
return result;
}
// ----------------------------------------------------------------------------
ClassImp(CbmTaskDigiEventQa)
......@@ -12,19 +12,20 @@
#include <vector>
//#include "EventBuilder.h"
#include "algo/qa/DigiEventQa.h"
#include "algo/qa/Histo1D.h"
class CbmDigiManager;
class CbmRecoConfig;
class TH1F;
class TH1D;
/** @class CbmTaskDigiEventQa
** @brief QA task class for digi events produced by the event builder
** @author Volker Friese <v.friese@gsi.de>
** @since 15.03.2022
**
** Currently implemented functionality: histogram the STS digi time within each event.
** To be expanded for other detector systems and probably more QA figures.
** Currently implemented functionality: histogram the digi time within each event.
** To be expanded for more QA figures.
** The histograms are published to the THttpServer.
**/
class CbmTaskDigiEventQa : public FairTask {
......@@ -67,13 +68,25 @@ private: // methods
virtual InitStatus Init();
/** @brief Create a ROOT TH1D from a Histo1D object
** @param Source histogram
** @param ROOT histogram
*/
TH1D* ToTH1D(const cbm::algo::Histo1D& source);
private: // members
const std::vector<CbmDigiEvent>* fEvents = nullptr; //! Input data (events)
size_t fNumTs = 0; ///< Number of processed timeslices
size_t fNumEvents = 0; ///< Number of analysed events
size_t fNumDigis = 0; ///< Number of analysed digis
double fExecTime = 0.; ///< Execution time [s]
TH1F* fHistDigiTimeSts = nullptr; //!
cbm::algo::qa::DigiEventQa fAlgo = {};
cbm::algo::qa::DigiEventQaConfig fConfig = {};
// ---- Histograms with digi times
std::map<ECbmModuleId, TH1D*> fDigiTimeHistos = {};
ClassDef(CbmTaskDigiEventQa, 1);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment