diff --git a/Common/Utils/include/CommonUtils/BoostHistogramUtils.h b/Common/Utils/include/CommonUtils/BoostHistogramUtils.h index f8d6e2d62fdf7..baa3c934a8bd4 100644 --- a/Common/Utils/include/CommonUtils/BoostHistogramUtils.h +++ b/Common/Utils/include/CommonUtils/BoostHistogramUtils.h @@ -24,6 +24,7 @@ #include #include +#include "Framework/Logger.h" #include "Rtypes.h" #include "TLinearFitter.h" #include "TVectorD.h" @@ -43,6 +44,9 @@ #include #include +using boostHisto2d = boost::histogram::histogram, boost::histogram::axis::regular>, boost::histogram::unlimited_storage>>; +using boostHisto1d = boost::histogram::histogram>, boost::histogram::unlimited_storage>>; + namespace o2 { namespace utils @@ -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 -std::vector fitGaus(Iterator first, Iterator last, BinCenterView axisfirst, const bool ignoreUnderOerflowBin = true) +std::vector 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(); fitter.StoreData(kFALSE); fitter.ClearPoints(); TVectorD par(3); @@ -261,14 +265,16 @@ std::vector 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; } @@ -330,12 +336,12 @@ std::vector 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 @@ -344,8 +350,10 @@ std::vector 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; } @@ -378,6 +386,7 @@ std::vector fitGaus(Iterator first, Iterator last, BinCenterView axisfir result.at(2) = binWidth / TMath::Sqrt(12); result.at(4) = -1; } + return result; } @@ -388,10 +397,10 @@ std::vector fitBoostHistoWithGaus(boost::histogram::histogram& } /// \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 diff --git a/Common/Utils/src/BoostHistogramUtils.cxx b/Common/Utils/src/BoostHistogramUtils.cxx index 231fdb3f0ebd9..6d34371981d69 100644 --- a/Common/Utils/src/BoostHistogramUtils.cxx +++ b/Common/Utils/src/BoostHistogramUtils.cxx @@ -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(); @@ -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."; -} \ No newline at end of file +} +} // namespace utils +} // namespace o2 \ No newline at end of file diff --git a/Detectors/EMCAL/calibration/CMakeLists.txt b/Detectors/EMCAL/calibration/CMakeLists.txt index 6c5b1c60cd2fc..471df7c472efa 100644 --- a/Detectors/EMCAL/calibration/CMakeLists.txt +++ b/Detectors/EMCAL/calibration/CMakeLists.txt @@ -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 @@ -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() \ No newline at end of file diff --git a/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h b/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h index a6f649f2919b1..0d8e1df8c08fa 100644 --- a/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h +++ b/Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h @@ -26,6 +26,10 @@ #include "EMCALBase/Geometry.h" #include +#if (defined(WITH_OPENMP) && !defined(__CLING__)) +#include +#endif + namespace o2 { namespace emcal @@ -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. @@ -88,16 +95,46 @@ class EMCALCalibExtractor /// \brief For now a dummy function to calibrate the time. template - o2::emcal::TimeCalibrationParams calibrateTime(boost::histogram::histogram& hist) + o2::emcal::TimeCalibrationParams calibrateTime(boost::histogram::histogram& 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); + } + 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(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); }; diff --git a/Detectors/EMCAL/calibration/run/runCalibOffline.cxx b/Detectors/EMCAL/calibration/run/runCalibOffline.cxx new file mode 100644 index 0000000000000..0270c20b6f6a3 --- /dev/null +++ b/Detectors/EMCAL/calibration/run/runCalibOffline.cxx @@ -0,0 +1,168 @@ +// 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 runTimeCalibOffline.cxx +/// \author + +#include +#include + +#include "FairLogger.h" + +#include "CommonUtils/BoostHistogramUtils.h" +#include "EMCALCalib/TimeCalibrationParams.h" +#include "EMCALCalib/CalibDB.h" +#include "EMCALCalibration/EMCALCalibExtractor.h" + +#include "CCDB/CcdbApi.h" + +#include "TFile.h" +#include "TH1.h" +#include "TH2.h" + +using CcdbApi = o2::ccdb::CcdbApi; +namespace bpo = boost::program_options; + +int main(int argc, char** argv) +{ + + bpo::variables_map vm; + bpo::options_description opt_general(""); + bpo::options_description opt_hidden(""); + bpo::options_description opt_all; + bpo::positional_options_description opt_pos; + + std::string CalibInputPath; + std::string ccdbServerPath; + bool doBadChannelCalib; + std::string nameCalibInputHist; // hCellIdVsTimeAbove300 for time, hCellIdVsEnergy for bad channel + + unsigned int nthreads; // number of threads used by openMP + + unsigned long rangestart; //30/10/2021, 01:02:32 for run 505566 -> 1635548552000 + unsigned long rangeend; // 30/10/2021, 02:31:10 for run 505566 -> 1635553870000 + + try { + bpo::options_description desc("Allowed options"); + desc.add_options()("help", "Print this help message")("CalibInputPath", bpo::value()->required(), "Set root input histogram")("ccdbServerPath", bpo::value()->default_value(o2::base::NameConf::getCCDBServer()), "Set path to ccdb server")("mode", bpo::value()->required(), "Set if time or bad channel calib")("nameInputHisto", bpo::value()->default_value("hCellIdVsTimeAbove300"), "Set name of input histogram")("nthreads", bpo::value()->default_value(1), "Set number of threads for OpenMP")("timestampStart", bpo::value()->default_value(1635548552000), "Set timestamp from start of run")("timestampEnd", bpo::value()->default_value(1635553870000), "Set timestamp from end of run"); + + bpo::store(bpo::parse_command_line(argc, argv, desc), vm); + + if (vm.count("help")) { + std::cout << "Please specify: CalibInputPath :Path to TFile with input histograms: \n mode: time or badchannel \n nameInputHisto: name of input calibration histogram\n"; + } + + if (vm.count("CalibInputPath")) { + std::cout << "CalibInputPath was set to " + << vm["CalibInputPath"].as() << ".\n"; + CalibInputPath = vm["CalibInputPath"].as(); + } else { + std::cout << "CalibInputPath was not set...\n"; + } + + if (vm.count("ccdbServerPath")) { + std::cout << "ccdbServerPath was set to " + << vm["ccdbServerPath"].as() << ".\n"; + ccdbServerPath = vm["ccdbServerPath"].as(); + } else { + printf("ccdbServerPath was not set.\nWill use standard path %s", ccdbServerPath.c_str()); + } + + if (vm.count("mode")) { + std::cout << "mode was set to " + << vm["mode"].as() << ".\n"; + std::string smode = vm["mode"].as(); + if (smode.find("time") != std::string::npos) { + std::cout << "performing time calibration" << std::endl; + doBadChannelCalib = false; + } else if (smode.find("badchannel") != std::string::npos) { + std::cout << "performing bad channel calibration" << std::endl; + doBadChannelCalib = true; + } else { + std::cout << "mode not set... returning\n"; + return 0; + } + } + + if (vm.count("nameInputHisto")) { + std::cout << "nameInputHisto was set to " + << vm["nameInputHisto"].as() << ".\n"; + nameCalibInputHist = vm["nameInputHisto"].as(); + } + + if (vm.count("nthreads")) { + std::cout << "number of threads was set to " + << vm["nthreads"].as() << ".\n"; + nthreads = vm["nthreads"].as(); + } + + if (vm.count("timestampStart")) { + std::cout << "timestampStart was set to " + << vm["timestampStart"].as() << ".\n"; + rangestart = vm["timestampStart"].as(); + } + + if (vm.count("timestampEnd")) { + std::cout << "timestampEnd was set to " + << vm["timestampEnd"].as() << ".\n"; + rangeend = vm["timestampEnd"].as(); + } + + } catch (bpo::error& e) { + std::cerr << "ERROR: " << e.what() << std::endl + << std::endl; + std::cerr << opt_general << std::endl; + exit(1); + } catch (std::exception& e) { + std::cerr << e.what() << ", application will now exit" << std::endl; + exit(2); + } + + // Set input file and get histogram + TFile* fTimeCalibInput = TFile::Open(CalibInputPath.c_str()); + if (!fTimeCalibInput) { + printf("%s not there... returning\n", CalibInputPath.c_str()); + return 0; + } + + TH2D* hCalibInputHist_ROOT = (TH2D*)fTimeCalibInput->Get(nameCalibInputHist.c_str()); + if (!hCalibInputHist_ROOT) { + printf("%s not there... returning\n", nameCalibInputHist.c_str()); + return 0; + } + + // instance of the calib extractor + o2::emcal::EMCALCalibExtractor CalibExtractor; + CalibExtractor.setNThreads(nthreads); + + // get boost histo from root input histogram + auto hCalibInputHist = o2::utils::boostHistoFromRoot_2D(hCalibInputHist_ROOT); + + // instance of CalibDB + o2::emcal::CalibDB calibdb(ccdbServerPath); + + if (doBadChannelCalib) { + printf("perform bad channel analysis\n"); + o2::emcal::BadChannelMap BCMap; + + // BCMap = CalibExtractor.calibrateBadChannels(hCalibInputHist); + } else { + printf("perform time calibration analysis\n"); + + // calibrate the time + o2::emcal::TimeCalibrationParams TCparams; + TCparams = CalibExtractor.calibrateTime(hCalibInputHist); + + // store parameters in ccdb via emcal calibdb + std::map metadata; + calibdb.storeTimeCalibParam(&TCparams, metadata, rangestart, rangeend); + } +} \ No newline at end of file diff --git a/Detectors/EMCAL/calibration/src/EMCALCalibExtractor.cxx b/Detectors/EMCAL/calibration/src/EMCALCalibExtractor.cxx index 6aa924655cece..d58794487c940 100644 --- a/Detectors/EMCAL/calibration/src/EMCALCalibExtractor.cxx +++ b/Detectors/EMCAL/calibration/src/EMCALCalibExtractor.cxx @@ -35,8 +35,9 @@ boostHisto EMCALCalibExtractor::buildHitAndEnergyMean(double emin, double emax, double meanVal = result.at(1); double sumVal = result.at(3); //..Set the values only for cells that are not yet marked as bad - if (sumVal > 0.) + if (sumVal > 0.) { eSumHisto(cellID, meanVal / (sumVal)); //..average energy per hit + } } return eSumHisto; } @@ -170,13 +171,15 @@ boostHisto EMCALCalibExtractor::buildHitAndEnergyMeanScaled(double emin, double auto energyAxis = boost::histogram::axis::regular<>(100, 0, 100, "t-texp"); Double_t E = energyAxis.value(j); Double_t N = hEnergyScaled.at(j, cell); - if (E < emin || E > emax) + if (E < emin || E > emax) { continue; + } Esum += E * N; Nsum += N; } - if (Nsum > 0.) + if (Nsum > 0.) { eSumHistoScaled(cell, Esum / (Nsum)); //..average energy per hit + } } return eSumHistoScaled;