Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 20 additions & 11 deletions Common/Utils/include/CommonUtils/BoostHistogramUtils.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include <TH2.h>
#include <TFile.h>

#include "Framework/Logger.h"
#include "Rtypes.h"
#include "TLinearFitter.h"
#include "TVectorD.h"
Expand All @@ -43,6 +44,9 @@
#include <boost/histogram/make_histogram.hpp>
#include <boost/histogram/accumulators/mean.hpp>

using boostHisto2d = boost::histogram::histogram<std::tuple<boost::histogram::axis::regular<double, boost::use_default, boost::use_default, boost::use_default>, boost::histogram::axis::regular<double, boost::use_default, boost::use_default, boost::use_default>>, boost::histogram::unlimited_storage<std::allocator<char>>>;
using boostHisto1d = boost::histogram::histogram<std::tuple<boost::histogram::axis::regular<double, boost::use_default, boost::use_default, boost::use_default>>, boost::histogram::unlimited_storage<std::allocator<char>>>;

namespace o2
{
namespace utils
Expand Down Expand Up @@ -215,24 +219,24 @@ std::string createErrorMessage(o2::utils::FitGausError_t errorcode);
/// \param first begin iterator of the histogram
/// \param last end iterator of the histogram
/// \param axisFirst axis iterator over the bin centers
/// \param ignoreUnderOerflowBin switch to disable taking under and overflow bin into fit
/// \param ignoreUnderOverflowBin switch to disable taking under and overflow bin into fit
/// \return result
/// result is of the type fitResult, which contains 4 parameters (0-Constant, 1-Mean, 2-Sigma, 3-Sum)
///
/// ** Temp Note: For now we forgo the templated struct in favor of a std::vector in order to
/// have this compile while we are working out the details
template <typename T, typename Iterator, typename BinCenterView>
std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfirst, const bool ignoreUnderOerflowBin = true)
std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfirst, const bool ignoreUnderOverflowBin = true)
{

if (ignoreUnderOerflowBin) {
if (ignoreUnderOverflowBin) {
first++;
last--;
}

TLinearFitter fitter(3, "pol2");
TMatrixD mat(3, 3);
Double_t kTol = mat.GetTol();
static TLinearFitter fitter(3, "pol2");
static TMatrixD mat(3, 3);
static Double_t kTol = mat.GetTol();
Comment on lines 237 to 239
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suspect this is to circumvent the problem with the multiple definition of TFormula when running with OMP threads > 1. This could become tricky in a multi-threaded environment. As far as I understand the variables should be allocated per thread (thread-local storage) based on our build settings. However would be good if we could find a way to avoid making members static.

fitter.StoreData(kFALSE);
fitter.ClearPoints();
TVectorD par(3);
Expand Down Expand Up @@ -261,14 +265,16 @@ std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfir
}
}

// TODO: Check why this is needed
if (*xMax < 4) {
LOG(error) << "Gaus fit failed! xMax < 4";
throw FitGausError_t::FIT_ERROR;
}
if (entries < 12) {
LOG(error) << "Gaus fit failed! entries < 12";
throw FitGausError_t::FIT_ERROR;
}
if (rms < kTol) {
LOG(error) << "Gaus fit failed! rms < kTol";
throw FitGausError_t::FIT_ERROR;
}

Expand Down Expand Up @@ -330,12 +336,12 @@ std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfir
}

if (TMath::Abs(par[1]) < kTol) {
LOG(error) << "Gaus fit failed! TMath::Abs(par[1]) < kTol";
throw FitGausError_t::FIT_ERROR;
;
}
if (TMath::Abs(par[2]) < kTol) {
LOG(error) << "Gaus fit failed! TMath::Abs(par[2]) < kTol";
throw FitGausError_t::FIT_ERROR;
;
}

// calculate parameters for gaus from pol2 fit
Expand All @@ -344,8 +350,10 @@ std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfir
result.at(2) = T(1. / TMath::Sqrt(TMath::Abs(-2. * par[2])));
auto lnparam0 = par[0] - par[1] * par[1] / (4 * par[2]);
if (lnparam0 > 307) {
LOG(error) << "Gaus fit failed! lnparam0 > 307";
throw FitGausError_t::FIT_ERROR;
}

result.at(0) = T(TMath::Exp(lnparam0));
return result;
}
Expand Down Expand Up @@ -378,6 +386,7 @@ std::vector<double> fitGaus(Iterator first, Iterator last, BinCenterView axisfir
result.at(2) = binWidth / TMath::Sqrt(12);
result.at(4) = -1;
}

return result;
}

Expand All @@ -388,10 +397,10 @@ std::vector<double> fitBoostHistoWithGaus(boost::histogram::histogram<axes...>&
}

/// \brief Convert a 1D root histogram to a Boost histogram
decltype(auto) boosthistoFromRoot_1D(TH1D* inHist1D);
boostHisto1d boosthistoFromRoot_1D(TH1D* inHist1D);

// \brief Convert a 2D root histogram to a Boost histogram
decltype(auto) boostHistoFromRoot_2D(TH2D* inHist2D);
boostHisto2d boostHistoFromRoot_2D(TH2D* inHist2D);

/// \brief Function to project 2d boost histogram onto x-axis
/// \param hist2d 2d boost histogram
Expand Down
48 changes: 26 additions & 22 deletions Common/Utils/src/BoostHistogramUtils.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,28 @@

#include "CommonUtils/BoostHistogramUtils.h"

// \brief Convert a 2D root histogram to a Boost histogram
decltype(auto) boostHistoFromRoot_2D(TH2D* inHist2D)
namespace o2
{
namespace utils
{

boostHisto1d boosthistoFromRoot_1D(TH1D* inHist1D)
{
// first setup the proper boost histogram
int nBins = inHist1D->GetNbinsX();
int xMin = inHist1D->GetXaxis()->GetXmin();
int xMax = inHist1D->GetXaxis()->GetXmax();
const char* title = inHist1D->GetXaxis()->GetTitle();
auto mHisto = boost::histogram::make_histogram(boost::histogram::axis::regular<>(nBins, xMin, xMax, title));

// trasfer the acutal values
for (Int_t x = 1; x < nBins + 1; x++) {
mHisto.at(x - 1) = inHist1D->GetBinContent(x);
}
return mHisto;
}

boostHisto2d boostHistoFromRoot_2D(TH2D* inHist2D)
{
// first setup the proper boost histogram
int nBinsX = inHist2D->GetNbinsX();
Expand All @@ -31,32 +51,16 @@ decltype(auto) boostHistoFromRoot_2D(TH2D* inHist2D)
// trasfer the acutal values
for (Int_t x = 1; x < nBinsX + 1; x++) {
for (Int_t y = 1; y < nBinsY + 1; y++) {
mHisto.at(x, y) = inHist2D->GetBinContent(x, y);
mHisto.at(x - 1, y - 1) = inHist2D->GetBinContent(x, y);
}
}
return mHisto;
}

/// \brief Convert a 1D root histogram to a Boost histogram
decltype(auto) boosthistoFromRoot_1D(TH1D* inHist1D)
{
// first setup the proper boost histogram
int nBins = inHist1D->GetNbinsX();
int xMin = inHist1D->GetXaxis()->GetXmin();
int xMax = inHist1D->GetXaxis()->GetXmax();
const char* title = inHist1D->GetXaxis()->GetTitle();
auto mHisto = boost::histogram::make_histogram(boost::histogram::axis::regular<>(nBins, xMin, xMax, title));

// trasfer the acutal values
for (Int_t x = 1; x < nBins + 1; x++) {
mHisto.at(x) = inHist1D->GetBinContent(x);
}
return mHisto;
}

/// \brief Printing an error message when then fit returns an invalid result
/// \param errorcode Error of the type FitGausError_t, thrown when fit result is invalid.
std::string createErrorMessage(o2::utils::FitGausError_t errorcode)
{
return "[Error]: Fit return an invalid result.";
}
}
} // namespace utils
} // namespace o2
25 changes: 17 additions & 8 deletions Detectors/EMCAL/calibration/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,10 @@
# or submit itself to any jurisdiction.

o2_add_library(EMCALCalibration
SOURCES src/EMCALChannelData.cxx
src/EMCALTimeCalibData.cxx
TARGETVARNAME targetName
SOURCES src/EMCALChannelData.cxx
src/EMCALTimeCalibData.cxx
src/EMCALCalibExtractor.cxx
PUBLIC_LINK_LIBRARIES O2::CCDB O2::EMCALBase
O2::EMCALCalib
O2::CommonUtils
Expand All @@ -21,21 +23,28 @@ o2_add_library(EMCALCalibration
O2::Algorithm
Microsoft.GSL::GSL
)




o2_target_root_dictionary(EMCALCalibration
HEADERS include/EMCALCalibration/EMCALChannelCalibrator.h
include/EMCALCalibration/EMCALChannelData.h
include/EMCALCalibration/EMCALTimeCalibData.h
LINKDEF src/EMCALCalibrationLinkDef.h)




o2_add_executable(emcal-channel-calib-workflow
COMPONENT_NAME calibration
SOURCES testWorkflow/emc-channel-calib-workflow.cxx
PUBLIC_LINK_LIBRARIES O2::Framework
O2::EMCALCalibration
O2::DetectorsCalibration)

o2_add_executable(run-calib-offline
COMPONENT_NAME emcal
TARGETVARNAME targetName
SOURCES run/runCalibOffline.cxx
PUBLIC_LINK_LIBRARIES O2::EMCALCalibration)


if (OpenMP_CXX_FOUND)
target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)
target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX)
endif()
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,10 @@
#include "EMCALBase/Geometry.h"
#include <boost/histogram.hpp>

#if (defined(WITH_OPENMP) && !defined(__CLING__))
#include <omp.h>
#endif

namespace o2
{
namespace emcal
Expand All @@ -45,6 +49,9 @@ class EMCALCalibExtractor
int getNsigma() const { return mSigma; }
void setNsigma(int ns) { mSigma = ns; }

void setNThreads(int n) { mNThreads = std::min(n, NCELLS); }
int getNThreads() const { return mNThreads; }

void setUseScaledHistoForBadChannels(bool useScaledHistoForBadChannels) { mUseScaledHistoForBadChannels = useScaledHistoForBadChannels; }

/// \brief Average energy per hit is caluclated for each cell.
Expand Down Expand Up @@ -88,16 +95,46 @@ class EMCALCalibExtractor

/// \brief For now a dummy function to calibrate the time.
template <typename... axes>
o2::emcal::TimeCalibrationParams calibrateTime(boost::histogram::histogram<axes...>& hist)
o2::emcal::TimeCalibrationParams calibrateTime(boost::histogram::histogram<axes...>& hist, double minTime = 0, double maxTime = 1000)
{

auto histReduced = boost::histogram::algorithm::reduce(hist, boost::histogram::algorithm::shrink(minTime, maxTime), boost::histogram::algorithm::shrink(0, NCELLS));

o2::emcal::TimeCalibrationParams TCP;
TCP.addTimeCalibParam(1234, 600, 0);

#if (defined(WITH_OPENMP) && !defined(__CLING__))
if (mNThreads < 1) {
mNThreads = std::min(omp_get_max_threads(), NCELLS);
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

17644 threads is probably out of order for local machines :) Needs to be checked with @shahor02 and @chiarazampolli where to put the maximum here, I guess omp threads will need to be shared among components.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes NCells is unrealistic. But like this it will use all available resources when calling the omp_get_max_threads(), right?

}
LOG(info) << "Number of threads that will be used = " << mNThreads;
#pragma omp parallel for num_threads(mNThreads)
#else
LOG(info) << "OPEN MP will not be used for the time calibration";
mNThreads = 1;
#endif

for (unsigned int i = 0; i < NCELLS; ++i) {
// project boost histogram to 1d just for 1 cell
auto boostHist1d = o2::utils::ProjectBoostHistoX(histReduced, i, i + 1);
// maybe cut histo to max +- 25 ns

// fit with gaussian to extract mean
try {
auto fitValues = o2::utils::fitBoostHistoWithGaus<double>(boostHist1d);
double mean = fitValues.at(1);
// add mean to time calib params
TCP.addTimeCalibParam(i, mean, 0);
} catch (o2::utils::FitGausError_t) {
TCP.addTimeCalibParam(i, 400, 0); // 400 ns shift default value
}
}
return TCP;
}

private:
bool mUseScaledHistoForBadChannels = false; ///< variable to specify whether or not we want to use the scaled histo for the claibration of bad channels.
int mSigma = 4; ///< number of sigma used in the calibration to define outliers
int mNThreads = 1; ///< number of threads used for calibration

ClassDefNV(EMCALCalibExtractor, 1);
};
Expand Down
Loading