From 3f736762aa7fd956a28b4cd56b89e668f7920021 Mon Sep 17 00:00:00 2001 From: Andreas Molander Date: Wed, 14 May 2025 11:26:45 +0300 Subject: [PATCH] [AFIT-109] Add MIP peak monitoring for FT0 --- .../Common/include/FITCommon/HelperCommon.h | 16 +- Modules/FIT/FIT/CMakeLists.txt | 2 + Modules/FIT/FIT/README.md | 46 +++ Modules/FIT/FIT/include/FIT/LinkDef.h | 1 + Modules/FIT/FIT/include/FIT/MIPCheck.h | 101 +++++ Modules/FIT/FIT/src/MIPCheck.cxx | 382 ++++++++++++++++++ Modules/FIT/FT0/etc/ft0-post-processing.json | 148 +++++++ Modules/FIT/FT0/include/FT0/PostProcTask.h | 24 ++ Modules/FIT/FT0/src/PostProcTask.cxx | 40 +- 9 files changed, 756 insertions(+), 4 deletions(-) create mode 100644 Modules/FIT/FIT/README.md create mode 100644 Modules/FIT/FIT/include/FIT/MIPCheck.h create mode 100644 Modules/FIT/FIT/src/MIPCheck.cxx diff --git a/Modules/FIT/Common/include/FITCommon/HelperCommon.h b/Modules/FIT/Common/include/FITCommon/HelperCommon.h index e6794b12b8..ab95d045bd 100644 --- a/Modules/FIT/Common/include/FITCommon/HelperCommon.h +++ b/Modules/FIT/Common/include/FITCommon/HelperCommon.h @@ -26,6 +26,8 @@ #include #include +#include +#include #include #include "QualityControl/QcInfoLogger.h" @@ -50,7 +52,8 @@ inline ValueType getConfigFromPropertyTree(const boost::property_tree::ptree& co template ::value || - std::is_same::value || (std::is_integral::value && !std::is_same::value)>::type> + std::is_same::value || + std::is_integral::value>> inline auto parseParameters(const std::string& param, const std::string& del) { std::regex reg(del); @@ -58,11 +61,18 @@ inline auto parseParameters(const std::string& param, const std::string& del) std::vector vecResult; for (auto it = first; it != last; it++) { if constexpr (std::is_integral::value && !std::is_same::value) { - vecResult.push_back(std::stoi(*it)); + if (!boost::algorithm::trim_copy(*it).empty()) { + vecResult.push_back(std::stoi(*it)); + } } else if constexpr (std::is_floating_point::value) { - vecResult.push_back(std::stod(*it)); + if (!boost::algorithm::trim_copy(*it).empty()) { + vecResult.push_back(std::stof(*it)); + } } else if constexpr (std::is_same::value) { vecResult.push_back(*it); + } else if constexpr (std::is_same::value) { + std::string trimmed = boost::algorithm::trim_copy(*it); + vecResult.push_back(boost::algorithm::to_lower_copy(trimmed) == "true" || trimmed == "1"); } } return vecResult; diff --git a/Modules/FIT/FIT/CMakeLists.txt b/Modules/FIT/FIT/CMakeLists.txt index e172db45cd..866045ac49 100644 --- a/Modules/FIT/FIT/CMakeLists.txt +++ b/Modules/FIT/FIT/CMakeLists.txt @@ -4,6 +4,7 @@ set(MODULE_NAME "O2QcFIT") add_library(${MODULE_NAME}) target_sources(${MODULE_NAME} PRIVATE src/LevelCheck.cxx) +target_sources(${MODULE_NAME} PRIVATE src/MIPCheck.cxx) target_sources(${MODULE_NAME} PRIVATE src/RawDataMetricTask.cxx) target_sources(${MODULE_NAME} PRIVATE src/RecoFITQcTask.cxx) @@ -28,6 +29,7 @@ install(TARGETS ${MODULE_NAME} add_root_dictionary(${MODULE_NAME} HEADERS include/FIT/LevelCheck.h + include/FIT/MIPCheck.h include/FIT/RawDataMetricTask.h include/FIT/RecoFITQcTask.h LINKDEF include/FIT/LinkDef.h) diff --git a/Modules/FIT/FIT/README.md b/Modules/FIT/FIT/README.md new file mode 100644 index 0000000000..f6cc9d9f10 --- /dev/null +++ b/Modules/FIT/FIT/README.md @@ -0,0 +1,46 @@ +# FIT quality control + +## Checks + +### MIPCheck + +The `MIPCheck` checks peaks in channel amplitude spectra. Essentially it is a check on parameters of Gaussian fits of the MIP peaks. + +#### Check parameters + +The parameters of the check are listed in the table below. For all parameters that take a list as argument, the items in the list corresponds to the different peaks; the first item is for the 1 MIP peak and so on. +If a threshold argument is omitted, that threshold will not be checked. This way, it is also possible to disable checks for a threshold for only some of the peaks. To disable a check for the first peak, but keep it for subsequent peaks, a negative value will also disable a check. + +| Name | Type | Default value | Comments | +|------------------------|------------------|---------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| +| `nameObjectToCheck` | `string` | `""` | Name of the MO to check | +| `nPeaksToFit` | `integer` | `2` | Number of MIP peaks to fit. The resulting fit funcion is a sum of `nPeaksToFit` Gaussians. | +| `gausParamsMeans` | list of `float` | `""` | Initial fit paramemters for the gaussian means. If omitted, the 1 MIP peak will default to the amplitude at the ampltude histogram max value. Other peaks will default to multiples of the initial 1 MIP peak mean (2 MIP peak mean = 2 * 1 MIP peak mean, and so on). The values are used as initial guesses and are not fixed. | +| `gausParamsSigma` | list of `float` | `"3.0, 7.0"` | The initial guesses for the MIP peak sigmas. Defaults to 5.0 for omitted peaks. | +| `fitRangeLow` | `float` | `11.0` | Lower limit of the fit | +| `fitRangeHigh` | `float` | `35.0` | Upper limit of the fit | +| `meanWarningsLow` | list of `float` | `""` | Lower warning thresholds for the MIP peak means | +| `meanWarningsHigh` | list of `float` | `""` | Upper warning thresholds for the MIP peak means | +| `meanErrorsLow` | list of `float` | `""` | Lower error thresholds for the MIP peak means | +| `meanErrorsHigh` | list of `float` | `""` | Upper error thresholds for the MIP peak means | +| `sigmaWarnings` | list of `float` | `""` | Sigma warning thresholds | +| `sigmaErrors` | list of `float` | `""` | Sigma error thresholds | +| `drawMeanWarningsLow` | list of `bool` | `"false, false"` | Whether to draw the lower warning thresholds for the MIP peak means | +| `drawMeanWarningsHigh` | list of `bool` | `"false, false"` | Whether to draw the upper warning thresholds for the MIP peak means | +| `drawMeanErrorsLow` | list of `bool` | `"true, true"` | Whether to draw the lower error thresholds for the MIP peak | +| `drawMeanErrorsHigh` | list of `bool` | `"true, true"` | Whether to draw the upper error thresholds for the MIP peak means | +| `drawSigmaWarnings` | list of `bool` | `"false, false"` | Whether to draw the sigma warning thresholds | +| `drawSigmaErrors` | list of `bool` | `"false, false"` | Whether to draw the sigma error thresholds | +| `labelPos` | list of `double` | `"0.15, 0.2, 0.85, 0.45"` | Position of the check label | + +#### Check logic + +The check produces the worst quality it finds among the listed checks below. For each failed check, a quality flag is added to the quality that explains the problem. + +For each peak (pseudo code): + +- if the fit fails -> `quality = BAD` +- `if not (meanWarningLow < mean < meanWarningHigh)` -> `quality = MEDIUM` +- `if not (meanErrorLow < mean < meanErrorHigh)` -> `quality = BAD` +- `if not sigma < sigmaWarning` -> `quality = WARNING` +- `if not sigma < sigmaError` -> `quality = BAD` diff --git a/Modules/FIT/FIT/include/FIT/LinkDef.h b/Modules/FIT/FIT/include/FIT/LinkDef.h index 8a4d8c74ab..fb2dbdc537 100644 --- a/Modules/FIT/FIT/include/FIT/LinkDef.h +++ b/Modules/FIT/FIT/include/FIT/LinkDef.h @@ -4,6 +4,7 @@ #pragma link off all functions; #pragma link C++ class o2::quality_control_modules::fit::LevelCheck + ; +#pragma link C++ class o2::quality_control_modules::fit::MIPCheck + ; #pragma link C++ class o2::quality_control_modules::fit::RawDataMetricTask + ; #pragma link C++ class o2::quality_control_modules::fit::RecoFITQcTask + ; diff --git a/Modules/FIT/FIT/include/FIT/MIPCheck.h b/Modules/FIT/FIT/include/FIT/MIPCheck.h new file mode 100644 index 0000000000..de2e4950d6 --- /dev/null +++ b/Modules/FIT/FIT/include/FIT/MIPCheck.h @@ -0,0 +1,101 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file MIPCheck.h +/// \author andreas.molander@cern.ch +/// + +#ifndef QC_MODULE_FIT_FITMIPCHECK_H +#define QC_MODULE_FIT_FITMIPCHECK_H + +#include "QualityControl/CheckInterface.h" + +namespace o2::quality_control_modules::fit +{ + +/// \brief QC check on MIP peaks in channel amplitude spectra. Essentially it is a check on +/// parameters of Gaussian fits of the MIP peaks. +/// \author andreas.molander@cern.ch +class MIPCheck : public o2::quality_control::checker::CheckInterface +{ + public: + MIPCheck() = default; + ~MIPCheck() override = default; + + void configure() override; + Quality check(std::map>* moMap) override; + void beautify(std::shared_ptr mo, Quality checkResult = Quality::Null) override; + std::string getAcceptedType() override; + void startOfActivity(const Activity& activity) override; + + private: + std::shared_ptr mActivity; + + std::string mNameObjectToCheck = ""; //< Name of the MO to check + + /// Number of MIP peaks to fit. The resulting fit funcion is a sum of `mNPeaksToFit` Gaussians. + int mNPeaksToFit = 2; + + /// Initial fit paramemters for the gaussian means. + /// If omitted, the 1 MIP peak will default to the amplitude at the ampltude histogram max value. + /// Other peaks will default to multiples of the initial 1 MIP peak mean (2 MIP peak mean = 2 * 1 MIP peak mean, and so on). + /// The values are used as initial guesses and are not fixed. + std::vector mGausParamsMeans; + + /// The initial fit parameters for the MIP peak sigmas. The values are used as initial guesses and are not fixed. + std::vector mGausParamsSigmas; + + float mFitRangeLow = 11.0; //< Lower limit of the fit + + float mFitRangeHigh = 35.0; //< Upper limit of the fit + + /// Lower warning thresholds for the MIP peak means. + /// The first element is for the first peak and so on. + std::vector mMeanWarningsLow; + + /// Upper warning thresholds for the MIP peak means. + /// The first element is for the first peak and so on. + std::vector mMeanWarningsHigh; + + /// Lower error thresholds for the MIP peak means. + /// The first element is for the first peak and so on. + std::vector mMeanErrorsLow; + + /// Upper error thresholds for the MIP peak means. + /// The first element is for the first peak and so on. + std::vector mMeanErrorsHigh; + + /// Sigma warning thresholds. + /// The first element is for the first peak and so on. + std::vector mSigmaWarnings; + + /// Sigma error thresholds. + /// The first element is for the first peak and so on. + std::vector mSigmaErrors; + + /// Whether to draw the threhold lines. + /// The first element is for the first peak and so on. + std::vector mDrawMeanWarningsLow; + std::vector mDrawMeanWarningsHigh; + std::vector mDrawMeanErrorsLow; + std::vector mDrawMeanErrorsHigh; + std::vector mDrawSigmaWarnings; + std::vector mDrawSigmaErrors; + + std::vector mVecLabelPos{ 0.15, 0.2, 0.85, 0.45 }; //< Position of the check label + + ClassDefOverride(MIPCheck, 3); +}; + +} // namespace o2::quality_control_modules::fit + +#endif // QC_MODULE_FIT_FITMIPCHECK_H diff --git a/Modules/FIT/FIT/src/MIPCheck.cxx b/Modules/FIT/FIT/src/MIPCheck.cxx new file mode 100644 index 0000000000..240d8563d0 --- /dev/null +++ b/Modules/FIT/FIT/src/MIPCheck.cxx @@ -0,0 +1,382 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file MIPCheck.cxx +/// \author andreas.molander@cern.ch +/// + +#include "FIT/MIPCheck.h" +#include "FITCommon/HelperCommon.h" +#include "QualityControl/MonitorObject.h" +#include "QualityControl/Quality.h" +#include "QualityControl/QcInfoLogger.h" + +// ROOT +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +using namespace std; +using namespace o2::quality_control; + +namespace o2::quality_control_modules::fit +{ + +void MIPCheck::configure() +{ + // Read configuration parameters + + mNameObjectToCheck = mCustomParameters.atOrDefaultValue("nameObjectToCheck", ""); + if (mNameObjectToCheck.empty()) { + ILOG(Error, Support) << "No MO name provided to check. The check will not do anything. Please provide a valid MO name." << ENDM; + } + mNPeaksToFit = std::stoi(mCustomParameters.atOrDefaultValue("nPeaksToFit", "2")); + mFitRangeLow = std::stof(mCustomParameters.atOrDefaultValue("fitRangeLow", "11.0")); + mFitRangeHigh = std::stof(mCustomParameters.atOrDefaultValue("fitRangeHigh", "35.0")); + + std::string parseableParam; + + parseableParam = mCustomParameters.atOrDefaultValue("gausParamsMean", ""); + mGausParamsMeans = helper::parseParameters(parseableParam, ","); + if (mGausParamsMeans.size() != mNPeaksToFit) { + ILOG(Warning, Support) << "Initial fit arguments for gaussian means not provided for all MIP peaks. Trying to estimate them..." << ENDM; + } + + parseableParam = mCustomParameters.atOrDefaultValue("gausParamsSigma", "3.0, 7.0"); + mGausParamsSigmas = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("meanWarningsLow", ""); + mMeanWarningsLow = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("meanWarningsHigh", ""); + mMeanWarningsHigh = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("meanErrorsLow", ""); + mMeanErrorsLow = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("meanErrorsHigh", ""); + mMeanErrorsHigh = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("sigmaWarnings", ""); + mSigmaWarnings = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("sigmaErrors", ""); + mSigmaErrors = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawMeanWarningsLow", "false, false"); + mDrawMeanWarningsLow = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawMeanWarningsHigh", "false, false"); + mDrawMeanWarningsHigh = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawMeanErrorsLow", "true, true"); + mDrawMeanErrorsLow = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawMeanErrorsHigh", "true, true"); + mDrawMeanErrorsHigh = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawSigmaWarnings", "false, false"); + mDrawSigmaWarnings = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("drawSigmaErrors", "false, false"); + mDrawSigmaErrors = helper::parseParameters(parseableParam, ","); + + parseableParam = mCustomParameters.atOrDefaultValue("labelPos", "0.15, 0.2, 0.85, 0.45"); + mVecLabelPos = helper::parseParameters(parseableParam, ","); + if (mVecLabelPos.size() != 4) { + ILOG(Error, Devel) << "Incorrect label coordinates! Setting default." << ENDM; + mVecLabelPos = { 0.15, 0.2, 0.85, 0.45 }; + } +} + +Quality MIPCheck::check(std::map>* moMap) +{ + Quality result = Quality::Null; + + for (auto& [moName, mo] : *moMap) { + if (mo->getName() == mNameObjectToCheck) { + auto* h = dynamic_cast(mo->getObject()); + if (h == nullptr) { + ILOG(Error, Support) << "Could not cast " << mNameObjectToCheck << " to TH1*, skipping" << ENDM; + continue; + } + // unless we find issues, we assume the quality is good + result = Quality::Good; + + // Fit the histogram to find the MIP peaks + + // Gaussian fit function + auto multiGausLambda = [this](Double_t* x, Double_t* par) -> Double_t { + Double_t sum = 0; + for (int i = 0; i < mNPeaksToFit; i++) { + // a = par[i * 3] + // mean = par[i * 3 + 1] + // sigma = par[i * 3 + 2]; + sum += par[i * 3] * TMath::Gaus(x[0], par[i * 3 + 1], par[i * 3 + 2]); + } + return sum; + }; + + // Construct a TF1 using the lambda + TF1 fitFunc("fitFunc", multiGausLambda, mFitRangeLow, mFitRangeHigh, 3 * mNPeaksToFit); + + // Provide initial guesses for the fit parameters + double mean1 = -1; + if (mGausParamsMeans.size() > 0) { + mean1 = mGausParamsMeans.at(0); + } else { + mean1 = h->GetBinCenter(h->GetMaximumBin()); + ILOG(Info, Devel) << Form("No initial guess for the 1 MIP peak mean provided. Estimated: %f", mean1) << ENDM; + } + + for (int i = 0; i < mNPeaksToFit; i++) { + double mean = -1; + if (i == 0) { + mean = mean1; + } else if (mGausParamsMeans.size() > i) { + mean = mGausParamsMeans.at(i); + } else { + mean = mean1 * (i + 1); + ILOG(Info, Devel) << Form("No initial guess for the %i MIP peak mean provided. Estimated: %f", i + 1, mean) << ENDM; + } + + int meanBin = h->FindBin(mean); + double a = h->GetBinContent(meanBin); + double sigma = mGausParamsSigmas.size() > i ? mGausParamsSigmas.at(i) : 5.0; // Use 5.0 for sigmas not provided + + ILOG(Info, Devel) << "Peak " << i << " initial fit parameters: a: " << a << ", mean: " << mean << ", sigma: " << sigma << ENDM; + + fitFunc.SetParameter(i * 3, a); + fitFunc.SetParameter(i * 3 + 1, mean); + fitFunc.SetParameter(i * 3 + 2, sigma); + } + + TFitResultPtr fitResult = h->Fit(&fitFunc, "SR"); + + // Check if the fit was successful + if (fitResult->Status() != 0) { + ILOG(Error, Devel) << "Fit failed for histogram " << moName << ", status: " << fitResult->Status() << ENDM; + result = Quality::Bad; + result.addFlag(FlagTypeFactory::UnknownQuality(), "Fit failed"); + break; // Don't bother checking the thresholds if the fit failed + } + + // Check the threholds, if some are not specified in the config, simply don't check them + for (int i = 0; i < mNPeaksToFit; i++) { + double mean = fitFunc.GetParameter(i * 3 + 1); + double sigma = fitFunc.GetParameter(i * 3 + 2); + + bool meanBad = false; + bool sigmaBad = false; + + // Check the error threholds + + // Mean + if (mMeanErrorsHigh.size() > i && mMeanErrorsHigh.at(i) >= 0 && mean > mMeanErrorsHigh.at(i)) { + result = Quality::Bad; + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak mean (%.2f) above upper error threshold (%.2f)", i + 1, mean, mMeanErrorsHigh.at(i))); + meanBad = true; + } else if (mMeanErrorsLow.size() > i && mMeanErrorsLow.at(i) >= 0 && mean < mMeanErrorsLow.at(i)) { + result = Quality::Bad; + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak mean (%.2f) below lower error threshold (%.2f)", i + 1, mean, mMeanErrorsLow.at(i))); + meanBad = true; + } + + // Sigma + if (mSigmaErrors.size() > i && mSigmaErrors.at(i) >= 0 && sigma > mSigmaErrors.at(i)) { + result = Quality::Bad; + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak sigma (%.2f) larger than error threshold (%.2f)", i + 1, sigma, mSigmaErrors.at(i))); + sigmaBad = true; + } + + // Check the warning thresholds + + // Mean (no need to check in case the mean already crossed an error threshold) + if (!meanBad && mMeanWarningsHigh.size() > i && mMeanWarningsHigh.at(i) >= 0 && mean > mMeanWarningsHigh.at(i)) { + if (result.isBetterThan(Quality::Medium)) { // Don't make the quality better than it was + result.set(Quality::Medium); + } + // Add a flag even if the quality was bad because of another peak + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak mean (%.2f) above upper warning threshold (%.2f)", i + 1, mean, mMeanWarningsHigh.at(i))); + } else if (!meanBad && mMeanWarningsLow.size() > i && mMeanWarningsLow.at(i) >= 0 && mean < mMeanWarningsLow.at(i)) { + if (result.isBetterThan(Quality::Medium)) { // Don't make the quality better than it was + result.set(Quality::Medium); + } + // Add a flag even if the quality was bad because of another peak + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak mean (%.2f) below lower warning threshold (%.2f)", i + 1, mean, mMeanWarningsLow.at(i))); + } + + // Sigma (no need to check in case the sigma already crossed an error threshold) + if (!sigmaBad && mSigmaWarnings.size() > i && mSigmaWarnings.at(i) >= 0 && sigma > mSigmaWarnings.at(i)) { + if (result.isBetterThan(Quality::Medium)) { // Don't make the quality better than it was + result.set(Quality::Medium); + } + result.addFlag(FlagTypeFactory::Unknown(), Form("%i MIP peak sigma (%.2f) larger than warning threshold (%.2f)", i + 1, sigma, mSigmaWarnings.at(i))); + } + } + + beautify(mo, result); + } + } + + return result; +} + +std::string MIPCheck::getAcceptedType() +{ + // This method is a remnant of early interface prototype and will be removed in the scope of ticket QC-373 + return "TH1"; +} + +void MIPCheck::beautify(std::shared_ptr mo, Quality checkResult) +{ + if (mo->getName() == mNameObjectToCheck) { + auto* h = dynamic_cast(mo->getObject()); + if (h == nullptr) { + ILOG(Error, Support) << "Could not cast `example` to TH1*, skipping" << ENDM; + return; + } + + TPaveText* msg = new TPaveText(mVecLabelPos[0], mVecLabelPos[1], mVecLabelPos[2], mVecLabelPos[3], "NDC"); + + std::vector meanValues; + std::vector sigmaValues; + + auto* func = h->GetFunction("fitFunc"); + if (func == nullptr) { + ILOG(Error, Support) << "Could not find fit function in histogram " << mo->getName() << ENDM; + msg->AddText("Gaussian fit not found!"); + } else { + msg->AddText("Gaussian fit"); + func->SetLineColor(kBlack); + + for (int i = 0; i < mNPeaksToFit; i++) { + double mean = func->GetParameter(i * 3 + 1); + double sigma = func->GetParameter(i * 3 + 2); + meanValues.push_back(mean); + sigmaValues.push_back(sigma); + + msg->AddText(Form("%i MIP: #mu = %.2f, #sigma = %.2f", i + 1, mean, sigma)); + } + } + + if (checkResult == Quality::Good) { + msg->AddText(">> Quality::Good <<"); + msg->SetFillColor(kGreen); + } else if (checkResult == Quality::Bad) { + msg->SetFillColor(kRed); + msg->AddText(">> Quality::Bad <<"); + } else if (checkResult == Quality::Medium) { + msg->SetFillColor(kOrange); + msg->AddText(">> Quality::Medium <<"); + } else if (checkResult == Quality::Null) { + msg->AddText(">> Quality::Null <<"); + msg->SetFillColor(kGray); + } + + for (const auto& [flag, comment] : checkResult.getFlags()) { + if (comment.empty()) { + msg->AddText(Form("#rightarrow Flag: %s", flag.getName().c_str())); + } else { + msg->AddText(Form("#rightarrow Flag: %s: %s", flag.getName().c_str(), comment.c_str())); + } + } + + h->GetListOfFunctions()->Add(msg); + + // Add threshold lines + + for (int i = 0; i < mNPeaksToFit; i++) { + // Warning low + if (mMeanWarningsLow.size() > i && mMeanWarningsLow.at(i) >= 0 && mDrawMeanWarningsLow.size() > i && mDrawMeanWarningsLow.at(i)) { + double meanWarningLow = mMeanWarningsLow.at(i); + auto* line = new TLine(meanWarningLow, 0, meanWarningLow, h->GetMaximum()); + line->SetLineColor(kOrange); + line->SetLineStyle(kDashed); + h->GetListOfFunctions()->Add(line); + } + + // Warning high + if (mMeanWarningsHigh.size() > i && mMeanWarningsHigh.at(i) >= 0 && mDrawMeanWarningsHigh.size() > i && mDrawMeanWarningsHigh.at(i)) { + double meanWarningHigh = mMeanWarningsHigh.at(i); + auto* line = new TLine(meanWarningHigh, 0, meanWarningHigh, h->GetMaximum()); + line->SetLineColor(kOrange); + line->SetLineStyle(kDashed); + h->GetListOfFunctions()->Add(line); + } + + // Error low + if (mMeanErrorsLow.size() > i && mMeanErrorsLow.at(i) >= 0 && mDrawMeanErrorsLow.size() > i && mDrawMeanErrorsLow.at(i)) { + double meanErrorLow = mMeanErrorsLow.at(i); + auto* line = new TLine(meanErrorLow, 0, meanErrorLow, h->GetMaximum()); + line->SetLineColor(kRed); + line->SetLineStyle(kDashed); + h->GetListOfFunctions()->Add(line); + } + + // Error high + if (mMeanErrorsHigh.size() > i && mMeanErrorsHigh.at(i) >= 0 && mDrawMeanErrorsHigh.size() > i && mDrawMeanErrorsHigh.at(i)) { + double meanErrorHigh = mMeanErrorsHigh.at(i); + auto* line = new TLine(meanErrorHigh, 0, meanErrorHigh, h->GetMaximum()); + line->SetLineColor(kRed); + line->SetLineStyle(kDashed); + h->GetListOfFunctions()->Add(line); + } + + // Sigma warning + if (mSigmaWarnings.size() > i && mSigmaWarnings.at(i) >= 0 && mDrawSigmaWarnings.size() > i && mDrawSigmaWarnings.at(i) && meanValues.size() > i) { + double sigmaWarningLow = meanValues.at(i) - mSigmaWarnings.at(i); + auto* lineLow = new TLine(sigmaWarningLow, 0, sigmaWarningLow, h->GetMaximum()); + lineLow->SetLineColor(kOrange); + lineLow->SetLineStyle(kDotted); + h->GetListOfFunctions()->Add(lineLow); + double sigmaWarningHigh = meanValues.at(i) + mSigmaWarnings.at(i); + auto* lineHigh = new TLine(sigmaWarningHigh, 0, sigmaWarningHigh, h->GetMaximum()); + lineLow->SetLineColor(kOrange); + lineLow->SetLineStyle(kDotted); + h->GetListOfFunctions()->Add(lineLow); + } + + // Sigma error + if (mSigmaErrors.size() > i && mSigmaErrors.at(i) >= 0 && mDrawSigmaErrors.size() > i && mDrawSigmaErrors.at(i) && meanValues.size() > i) { + double sigmaErrorLow = meanValues.at(i) - mSigmaErrors.at(i); + auto* lineLow = new TLine(sigmaErrorLow, 0, sigmaErrorLow, h->GetMaximum()); + lineLow->SetLineColor(kRed); + lineLow->SetLineStyle(kDotted); + h->GetListOfFunctions()->Add(lineLow); + double sigmaErrorHigh = meanValues.at(i) + mSigmaErrors.at(i); + auto* lineHigh = new TLine(sigmaErrorHigh, 0, sigmaErrorHigh, h->GetMaximum()); + lineHigh->SetLineColor(kRed); + lineHigh->SetLineStyle(kDotted); + h->GetListOfFunctions()->Add(lineHigh); + } + } + } +} + +void MIPCheck::startOfActivity(const Activity& activity) +{ + mActivity = make_shared(activity); +} + +} // namespace o2::quality_control_modules::fit diff --git a/Modules/FIT/FT0/etc/ft0-post-processing.json b/Modules/FIT/FT0/etc/ft0-post-processing.json index dd2e8d21fc..07546c96e2 100644 --- a/Modules/FIT/FT0/etc/ft0-post-processing.json +++ b/Modules/FIT/FT0/etc/ft0-post-processing.json @@ -23,6 +23,154 @@ } }, "checks": { + "ASideInnerMIPCheck": { + "active": "true", + "className": "o2::quality_control_modules::fit::MIPCheck", + "moduleName": "QcFT0", + "policy": "OnAny", + "detectorName": "FT0", + "dataSource": [ + { + "type": "PostProcessing", + "name": "PostProc", + "MOs": [ + "AmpAInner" + ] + } + ], + "checkParameters": { + "nameObjectToCheck": "AmpAInner", + "nPeaksToFit": "2", + "gausParamsMean": "14, 28", + "gausParamsSigma": "3, 7", + "fitRangeLow": "11", + "fitRangeHigh": "35", + "meanWarningsLow": "12, 26", + "meanWarningsHigh": "16, 30", + "meanErrorsLow": "10, 24", + "meanErrorsHigh": "18, 32", + "sigmaWarnings": "4, -1", + "sigmaErrors": "5, -1", + "drawMeanWarningsLow": "false, false", + "drawMeanWarningsHigh": "false, false", + "drawMeanErrorsLow": "true, true", + "drawMeanErrorsHigh": "true, true", + "drawSigmaWarnings": "false, false", + "drawSigmaErrors": "fasle, false", + "labelPos": "0.15, 0.2, 0.85, 0.45" + } + }, + "ASideOuterMIPCheck": { + "active": "true", + "className": "o2::quality_control_modules::fit::MIPCheck", + "moduleName": "QcFT0", + "policy": "OnAny", + "detectorName": "FT0", + "dataSource": [ + { + "type": "PostProcessing", + "name": "PostProc", + "MOs": [ + "AmpAOuter" + ] + } + ], + "checkParameters": { + "nameObjectToCheck": "AmpAOuter", + "nPeaksToFit": "2", + "gausParamsMean": "14, 28", + "gausParamsSigma": "3, 7", + "fitRangeLow": "11", + "fitRangeHigh": "35", + "meanWarningsLow": "12, 26", + "meanWarningsHigh": "16, 30", + "meanErrorsLow": "10, 24", + "meanErrorsHigh": "18, 32", + "sigmaWarnings": "4, -1", + "sigmaErrors": "5, -1", + "drawMeanWarningsLow": "false, false", + "drawMeanWarningsHigh": "false, false", + "drawMeanErrorsLow": "true, true", + "drawMeanErrorsHigh": "true, true", + "drawSigmaWarnings": "false, false", + "drawSigmaErrors": "fasle, false", + "labelPos": "0.15, 0.2, 0.85, 0.45" + } + }, + "CSideMIPCheck": { + "active": "true", + "className": "o2::quality_control_modules::fit::MIPCheck", + "moduleName": "QcFT0", + "policy": "OnAny", + "detectorName": "FT0", + "dataSource": [ + { + "type": "PostProcessing", + "name": "PostProc", + "MOs": [ + "AmpC" + ] + } + ], + "checkParameters": { + "nameObjectToCheck": "AmpC", + "nPeaksToFit": "2", + "gausParamsMean": "14, 28", + "gausParamsSigma": "3, 7", + "fitRangeLow": "11", + "fitRangeHigh": "35", + "meanWarningsLow": "12, 26", + "meanWarningsHigh": "16, 30", + "meanErrorsLow": "10, 24", + "meanErrorsHigh": "18, 32", + "sigmaWarnings": "4, -1", + "sigmaErrors": "5, -1", + "drawMeanWarningsLow": "false, false", + "drawMeanWarningsHigh": "false, false", + "drawMeanErrorsLow": "true, true", + "drawMeanErrorsHigh": "true, true", + "drawSigmaWarnings": "false, false", + "drawSigmaErrors": "fasle, false", + "labelPos": "0.15, 0.2, 0.85, 0.45" + } + }, + "AllChannelsMIPCheck": { + "active": "true", + "className": "o2::quality_control_modules::fit::MIPCheck", + "moduleName": "QcFT0", + "policy": "OnAny", + "detectorName": "FT0", + "dataSource": [ + { + "type": "PostProcessing", + "name": "PostProc", + "MOs": [ + "AmpNormPerChannel" + ] + } + ], + "checkParameters": { + "nameObjectToCheck": "AmpNormPerChannel", + "nPeaksToFit": "2", + "gausParamsMean": "14, 28", + "gausParamsSigma": "3, 7", + "fitRangeLow": "11", + "fitRangeHigh": "35", + "meanWarningsLow": "12, 26", + "meanWarningsHigh": "16, 30", + "meanErrorsLow": "10, 24", + "meanErrorsHigh": "18, 32", + "sigmaWarnings": "4, -1", + "sigmaErrors": "5, -1", + "drawMeanWarningsLow": "false, false", + "drawMeanWarningsHigh": "false, false", + "drawMeanErrorsLow": "true, true", + "drawMeanErrorsHigh": "true, true", + "drawSigmaWarnings": "false, false", + "drawSigmaErrors": "fasle, false", + "labelPos": "0.15, 0.2, 0.85, 0.45" + } + }, "OutOfBunchCollCheck": { "active": "true", "className": "o2::quality_control_modules::ft0::OutOfBunchCollCheck", diff --git a/Modules/FIT/FT0/include/FT0/PostProcTask.h b/Modules/FIT/FT0/include/FT0/PostProcTask.h index 3458c93c20..7a1b8b0fa3 100644 --- a/Modules/FIT/FT0/include/FT0/PostProcTask.h +++ b/Modules/FIT/FT0/include/FT0/PostProcTask.h @@ -35,6 +35,9 @@ #include #include +#include +#include + class TH1F; class TCanvas; class TLegend; @@ -79,6 +82,27 @@ class PostProcTask final : public quality_control::postprocessing::PostProcessin std::unique_ptr mTime; std::unique_ptr mHistStatsSideA; std::unique_ptr mHistStatsSideC; + + /// Sum of count in bins vs amplitude for detector channels 0 to 31 (A-side inner). + /// In other words a projection of range 1-32 in the 2D histogram amplitude vs channel + /// from the DigitQcTask. + std::unique_ptr mHistAmpAInner; + + /// Sum of count in bins vs amplitude for detector channels 32 to 95 (A-side outer). + /// In other words a projection of range 33-96 in the 2D histogram amplitude vs channel + /// from the DigitQcTask. + std::unique_ptr mHistAmpAOuter; + + /// Sum of count in bins vs amplitude for detector channels 96 to 207 (C-side). + /// In other words a projection of range 97-208 in the 2D histogram amplitude vs channel + /// from the DigitQcTask. + std::unique_ptr mHistAmpC; + + /// Sum of normalized count in bins vs amplitude for all detector channels (0-207). + /// In other words the sum of Y projections for each bin in the 2D histogram amplitude vs channel + /// from the DigitQcTask, where each projection is scaled with 1/(counts in that projection) + std::unique_ptr mHistAmpNormPerChannel; + ChannelGeometry mChannelGeometry; //! // Configurations int mLowTimeThreshold{ -192 }; diff --git a/Modules/FIT/FT0/src/PostProcTask.cxx b/Modules/FIT/FT0/src/PostProcTask.cxx index f77268f572..63571de6b6 100644 --- a/Modules/FIT/FT0/src/PostProcTask.cxx +++ b/Modules/FIT/FT0/src/PostProcTask.cxx @@ -32,6 +32,7 @@ #include "FITCommon/HelperCommon.h" #include +#include using namespace o2::quality_control::postprocessing; using namespace o2::quality_control_modules::fit; @@ -83,6 +84,10 @@ void PostProcTask::reset() mHistBcTrgOutOfBunchColl.reset(); mAmpl.reset(); mTime.reset(); + mHistAmpAInner.reset(); + mHistAmpAOuter.reset(); + mHistAmpC.reset(); + mHistAmpNormPerChannel.reset(); mMapHistsToDecompose.clear(); } @@ -104,7 +109,11 @@ void PostProcTask::initialize(Trigger trg, framework::ServiceRegistryRef service mHistTrgValidation = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "TrgValidation", "SW + HW only to validated triggers fraction", mMapTrgBits); mAmpl = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "MeanAmplPerChannel", "mean ampl per channel;Channel;Ampl #mu #pm #sigma", o2::ft0::Constants::sNCHANNELS_PM, 0, o2::ft0::Constants::sNCHANNELS_PM); mTime = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "MeanTimePerChannel", "mean time per channel;Channel;Time #mu #pm #sigma", o2::ft0::Constants::sNCHANNELS_PM, 0, o2::ft0::Constants::sNCHANNELS_PM); - + mHistAmpAInner = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "AmpAInner", "FT0 A-side inner channel amplitudes (ch 0-31);Channel amplitude (ADC ch);Counts", 4200, -100, 4100); + mHistAmpAOuter = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "AmpAOuter", "FT0 A-side outer channel amplitudes (ch 32-95);Channel amplitude (ADC ch);Counts", 4200, -100, 4100); + mHistAmpC = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "AmpC", "FT0 C-side channel ampltitudes (ch 96-207);Channel amplitude (ADC ch);Counts", 4200, -100, 4100); + mHistAmpNormPerChannel = helper::registerHist(getObjectsManager(), quality_control::core::PublicationPolicy::ThroughStop, "", "AmpNormPerChannel", "FT0 channel amplitudes normalized per channel (ch 0-207);Channel amplitude (ADC ch);#sum_{ch}AmpHist_{ch} #times (1/counts_{ch})", 4200, -100, 4100); + mHistAmpNormPerChannel->Sumw2(kFALSE); mChannelGeometry.init(-200., 200., -200., 200., 10.); // values - borders for hist and margin mHistStatsSideA = mChannelGeometry.makeHistSideA("GeoChannelStatA", "Channel occupancy, side-A"); @@ -175,6 +184,35 @@ void PostProcTask::update(Trigger trg, framework::ServiceRegistryRef serviceReg) std::unique_ptr projNom(hAmpPerChannel->ProjectionX("projNom", hAmpPerChannel->GetYaxis()->FindBin(1.0), -1)); std::unique_ptr projDen(hAmpPerChannel->ProjectionX("projDen")); mHistCFDEff->Divide(projNom.get(), projDen.get()); + + // Channel amplitudes for A-side inner, A-side outer, C-side and all channels + mHistAmpAInner->Reset(); + mHistAmpAOuter->Reset(); + mHistAmpC->Reset(); + mHistAmpNormPerChannel->Reset(); + + std::unique_ptr projAInner(hAmpPerChannel->ProjectionY("projAInner", 1, 32)); + std::unique_ptr projAOuter(hAmpPerChannel->ProjectionY("projAOuter", 33, 96)); + std::unique_ptr projC(hAmpPerChannel->ProjectionY("projC", 97, 208)); + + mHistAmpAInner->Add(projAInner.get()); + mHistAmpAOuter->Add(projAOuter.get()); + mHistAmpC->Add(projC.get()); + + // Iterate over channels + int entries = mHistAmpNormPerChannel->GetEntries(); + for (int iBin = 1; iBin < hAmpPerChannel->GetXaxis()->GetNbins() + 1; iBin++) { + std::unique_ptr ampForChannel(hAmpPerChannel->ProjectionY(Form("ampForChannel%i", iBin - 1), iBin, iBin)); + ampForChannel->Sumw2(kFALSE); + + // Here height of the bins of each detector channel must be divided by the sum of all bins (= number of events) in that channel. + // This is needed to equalize the weight of all channels independently of the load (central channels are loaded more). + const double integral = ampForChannel->Integral(); + const double scale = integral > 0. ? 1. / integral : 0.; + mHistAmpNormPerChannel->Add(ampForChannel.get(), scale); + entries += ampForChannel->GetEntries(); + } + mHistAmpNormPerChannel->SetEntries(entries); } // Times auto hTimePerChannel = mPostProcHelper.template getObject("TimePerChannel");