From 77aa39a8f52a29a04f8fddad2a240d1f6ed76569 Mon Sep 17 00:00:00 2001 From: Luisa Bergmann Date: Tue, 25 Jul 2023 14:14:46 +0200 Subject: [PATCH] Adding TRD t0 fitting procedure --- DataFormats/Detectors/TRD/CMakeLists.txt | 2 + .../TRD/include/DataFormatsTRD/CalT0.h | 23 ++ .../TRD/include/DataFormatsTRD/Constants.h | 1 + .../TRD/include/DataFormatsTRD/PHData.h | 1 + .../TRD/include/DataFormatsTRD/T0FitHistos.h | 61 ++++ .../Detectors/TRD/src/DataFormatsTRDLinkDef.h | 2 + DataFormats/Detectors/TRD/src/T0FitHistos.cxx | 75 +++++ Detectors/TRD/calibration/CMakeLists.txt | 2 + .../TRDCalibration/CalibrationParams.h | 4 + .../include/TRDCalibration/T0Fit.h | 111 +++++++ Detectors/TRD/calibration/src/T0Fit.cxx | 270 ++++++++++++++++++ .../calibration/src/TRDCalibrationLinkDef.h | 3 + Detectors/TRD/workflow/CMakeLists.txt | 1 + .../workflow/include/TRDWorkflow/T0FitSpec.h | 177 ++++++++++++ Detectors/TRD/workflow/io/CMakeLists.txt | 1 + .../include/TRDWorkflowIO/TRDPHReaderSpec.h | 53 ++++ .../TRD/workflow/io/src/TRDPHReaderSpec.cxx | 83 ++++++ .../TRD/workflow/src/trd-calib-workflow.cxx | 10 + 18 files changed, 880 insertions(+) create mode 100644 DataFormats/Detectors/TRD/include/DataFormatsTRD/T0FitHistos.h create mode 100644 DataFormats/Detectors/TRD/src/T0FitHistos.cxx create mode 100644 Detectors/TRD/calibration/include/TRDCalibration/T0Fit.h create mode 100644 Detectors/TRD/calibration/src/T0Fit.cxx create mode 100644 Detectors/TRD/workflow/include/TRDWorkflow/T0FitSpec.h create mode 100644 Detectors/TRD/workflow/io/include/TRDWorkflowIO/TRDPHReaderSpec.h create mode 100644 Detectors/TRD/workflow/io/src/TRDPHReaderSpec.cxx diff --git a/DataFormats/Detectors/TRD/CMakeLists.txt b/DataFormats/Detectors/TRD/CMakeLists.txt index 7697d12fec1b2..4689fa6a40fc9 100644 --- a/DataFormats/Detectors/TRD/CMakeLists.txt +++ b/DataFormats/Detectors/TRD/CMakeLists.txt @@ -14,6 +14,7 @@ o2_add_library(DataFormatsTRD src/LinkRecord.cxx src/AngularResidHistos.cxx src/GainCalibHistos.cxx + src/T0FitHistos.cxx src/DcsCcdbObjects.cxx src/Tracklet64.cxx src/RawData.cxx @@ -37,6 +38,7 @@ o2_target_root_dictionary(DataFormatsTRD include/DataFormatsTRD/DcsCcdbObjects.h include/DataFormatsTRD/AngularResidHistos.h include/DataFormatsTRD/GainCalibHistos.h + include/DataFormatsTRD/T0FitHistos.h include/DataFormatsTRD/HelperMethods.h include/DataFormatsTRD/Hit.h include/DataFormatsTRD/Digit.h diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalT0.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalT0.h index 30474888a7424..91f251b14b4d5 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalT0.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalT0.h @@ -32,11 +32,34 @@ class CalT0 ~CalT0() = default; void setT0(int iDet, float t0) { mT0[iDet] = t0; } + void setT0av(float t0) { mT0av = t0; } float getT0(int iDet) const { return mT0[iDet]; } + float getT0av() const { return mT0av; } + float calcT0av() const + { + if (mT0.size() == 0) { + return -1; + } + float sum = 0; + int counts = 0; + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + if (mT0[iDet] > -5) { + sum += mT0[iDet]; + ++counts; + } + } + + if (counts > 0) { + return sum / counts; + } else { + return -2; + } + } private: std::array mT0{}; ///< calibrated T0 per TRD chamber + float mT0av{-1}; ///< calibrated average T0 ClassDefNV(CalT0, 1); }; diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h index 04f05c402e667..103b830715fff 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h @@ -77,6 +77,7 @@ constexpr double VDRIFTDEFAULT = 1.546; ///< default value for vDrift constexpr double EXBDEFAULT = 0.0; ///< default value for LorentzAngle constexpr int NBINSGAINCALIB = 320; ///< number of bins in the charge (Q0+Q1+Q2) histogram for gain calibration constexpr float MPVDEDXDEFAULT = 42.; ///< default Most Probable Value of TRD dEdx +constexpr float T0DEFAULT = 1.2; ///< default value for t0 // array size to store incoming half cru payload. constexpr int HBFBUFFERMAX = 1048576; ///< max buffer size for data read from a half cru, (all events) diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h index de939ef0a0091..dbbc3161f9b02 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h @@ -13,6 +13,7 @@ #define ALICEO2_TRD_PHDATA_H_ #include +#include "Rtypes.h" namespace o2::trd { diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/T0FitHistos.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/T0FitHistos.h new file mode 100644 index 0000000000000..eff39bbbac438 --- /dev/null +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/T0FitHistos.h @@ -0,0 +1,61 @@ +// 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 T0FitHistos.h +/// \brief Class to store the TRD PH values for each TRD chamber + +#ifndef ALICEO2_T0FITHISTOS_H +#define ALICEO2_T0FITHISTOS_H + +#include "DataFormatsTRD/Constants.h" +#include "DataFormatsTRD/PHData.h" +#include "Framework/InputRecord.h" +#include "Rtypes.h" +#include +#include +#include + +namespace o2 +{ +namespace trd +{ + +class T0FitHistos +{ + public: + T0FitHistos() = default; + T0FitHistos(const T0FitHistos&) = default; + ~T0FitHistos() = default; + void init(); + void reset(); + auto getDetector(int index) const { return mDet[index]; } + auto getTimeBin(int index) const { return mTB[index]; } + auto getADC(int index) const { return mADC[index]; } + auto getNEntries() const { return mNEntriesTot; } + + void fill(const std::vector data); + void merge(const T0FitHistos* prev); + void print(); + + private: + std::vector mDet{}; + std::vector mTB{}; + std::vector mADC{}; + size_t mNEntriesTot{0}; + bool mInitialized{false}; + + ClassDefNV(T0FitHistos, 1); +}; + +} // namespace trd +} // namespace o2 + +#endif // ALICEO2_T0FITHISTOS_H diff --git a/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h b/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h index 8f258914d5a48..3c3b93e58ce40 100644 --- a/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h +++ b/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h @@ -33,6 +33,7 @@ #pragma link C++ class o2::trd::NoiseStatusMCM + ; #pragma link C++ class o2::trd::AngularResidHistos + ; #pragma link C++ class o2::trd::GainCalibHistos + ; +#pragma link C++ class o2::trd::T0FitHistos + ; #pragma link C++ class o2::trd::CalVdriftExB + ; #pragma link C++ class o2::trd::CalGain + ; #pragma link C++ class o2::trd::CalT0 + ; @@ -51,6 +52,7 @@ #pragma link C++ class std::vector < o2::trd::Digit> + ; #pragma link C++ class std::vector < o2::trd::AngularResidHistos> + ; #pragma link C++ class std::vector < o2::trd::GainCalibHistos> + ; +#pragma link C++ class std::vector < o2::trd::T0FitHistos> + ; #pragma link C++ class std::vector < o2::trd::PHData> + ; #pragma link C++ class std::vector < o2::trd::KrCluster> + ; #pragma link C++ class std::vector < o2::trd::KrClusterTriggerRecord> + ; diff --git a/DataFormats/Detectors/TRD/src/T0FitHistos.cxx b/DataFormats/Detectors/TRD/src/T0FitHistos.cxx new file mode 100644 index 0000000000000..3b62525393fe8 --- /dev/null +++ b/DataFormats/Detectors/TRD/src/T0FitHistos.cxx @@ -0,0 +1,75 @@ +// 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 T0FitHistos.cxx +/// \brief Class to store the TRD PH values for each chamber + +#include "DataFormatsTRD/T0FitHistos.h" +#include +#include + +using namespace o2::trd; +using namespace o2::trd::constants; + +void T0FitHistos::init() +{ + mDet.resize(0); + mTB.resize(0); + mADC.resize(0); + mInitialized = true; +} + +void T0FitHistos::reset() +{ + mDet.resize(0); + mTB.resize(0); + mADC.resize(0); + mNEntriesTot = 0; +} + +void T0FitHistos::fill(const std::vector data) +{ + + if (!mInitialized) { + init(); + } + + for (auto ph : data) { + int det = ph.getDetector(); + int tb = ph.getTimebin(); + int adc = ph.getADC(); + + if (ph.getNneighbours() != 2) { + continue; + } + + mDet.push_back(det); + mTB.push_back(tb); + mADC.push_back(adc); + ++mNEntriesTot; + } +} + +void T0FitHistos::merge(const T0FitHistos* prev) +{ + auto sizePrev = (int)prev->getNEntries(); + + for (int i = 0; i < sizePrev; ++i) { + mDet.push_back(prev->getDetector(i)); + mTB.push_back(prev->getTimeBin(i)); + mADC.push_back(prev->getADC(i)); + } +} + +void T0FitHistos::print() +{ + LOG(info) << "There are " << mNEntriesTot << " entries in the container"; +} diff --git a/Detectors/TRD/calibration/CMakeLists.txt b/Detectors/TRD/calibration/CMakeLists.txt index 4a9af0739a328..52444d2855b1f 100644 --- a/Detectors/TRD/calibration/CMakeLists.txt +++ b/Detectors/TRD/calibration/CMakeLists.txt @@ -16,6 +16,7 @@ o2_add_library(TRDCalibration src/CalibratorVdExB.cxx src/CalibratorGain.cxx src/CalibratorNoise.cxx + src/T0Fit.cxx src/PadCalibCCDBBuilder.cxx src/KrClusterFinder.cxx src/DCSProcessor.cxx @@ -35,6 +36,7 @@ o2_add_library(TRDCalibration include/TRDCalibration/CalibratorVdExB.h include/TRDCalibration/CalibratorGain.h include/TRDCalibration/CalibratorNoise.h + include/TRDCalibration/T0Fit.h include/TRDCalibration/CalibrationParams.h include/TRDCalibration/PadCalibCCDBBuilder.h include/TRDCalibration/KrClusterFinder.h diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h index 6ae7c6c727363..677673a7f85f3 100644 --- a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h +++ b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h @@ -39,6 +39,10 @@ struct TRDCalibParams : public o2::conf::ConfigurableParamHelper float dEdxTPCMin = 30.; float dEdxTPCMax = 70.; + // For t0 fits + size_t minEntriesTotalT0Fit = 1'000'000; ///< minimum number of entries in inclusive distribution for (meaningful) t0 fit + size_t minEntriesChamberT0Fit = 30'000; ///< minimum number of entries in one chamber for (meaningful) t0 fit + // Creation of PH plots float pileupCut = 0.7; diff --git a/Detectors/TRD/calibration/include/TRDCalibration/T0Fit.h b/Detectors/TRD/calibration/include/TRDCalibration/T0Fit.h new file mode 100644 index 0000000000000..ef08b3b9f80f6 --- /dev/null +++ b/Detectors/TRD/calibration/include/TRDCalibration/T0Fit.h @@ -0,0 +1,111 @@ +// 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 T0Fit.h +/// \brief Fits the TRD PH distributions to extract the t0 value +/// \author Luisa Bergmann + +#ifndef O2_TRD_T0FIT_H +#define O2_TRD_T0FIT_H + +#include "DataFormatsTRD/T0FitHistos.h" +#include "DetectorsCalibration/TimeSlotCalibration.h" +#include "DetectorsCalibration/TimeSlot.h" +#include "DataFormatsTRD/Constants.h" +#include "CCDB/CcdbObjectInfo.h" +#include "DataFormatsTRD/CalT0.h" +#include "TRDCalibration/CalibrationParams.h" + +#include "Rtypes.h" +#include "TH2.h" +#include "TH1.h" +#include "TProfile.h" +#include "TF1.h" +#include "Fit/Fitter.h" +#include "TFile.h" +#include "TTree.h" + +#include +#include +#include + +namespace o2 +{ +namespace trd +{ +//______________________________________________________________________________________________ +struct ErfLandauChi2Functor { + double operator()(const double* par) const; + std::vector x; ///< x-value (time-bin) of adc profile + std::vector y; ///< y-value (av. adc) for corresp. time-bin + float lowerBoundFit; ///< lower bound for fit + float upperBoundFit; ///< upper bound for fit +}; + +//______________________________________________________________________________________________ +class T0Fit final : public o2::calibration::TimeSlotCalibration +{ + using Slot = o2::calibration::TimeSlot; + + public: + T0Fit() = default; + ~T0Fit() final = default; + + bool hasEnoughData(const Slot& slot) const final { return slot.getContainer()->getNEntries() >= mMinEntriesTotal; } + void initOutput() final; + void finalizeSlot(Slot& slot) final; + Slot& emplaceNewSlot(bool front, TFType tStart, TFType tEnd) final; + + /// (Re-)Creates a file "trd_t0fit.root". This lets continually fill + /// a tree with the fit results. + void createOutputFile(); + + /// Close the output file. E.g. First write the tree to the file and let the + /// smart pointers take care of closing the file. + void closeOutputFile(); + + const std::vector& getCcdbObjectVector() const { return mObjectVector; } + std::vector& getCcdbObjectInfoVector() { return mInfoVector; } + + void initProcessing(); + + /// Initialize the fit values once with the previous valid ones if they are + /// available. + void retrievePrev(o2::framework::ProcessingContext& pc); + + private: + bool mInitDone{false}; ///< flag to avoid creating output etc multiple times + const TRDCalibParams& mParams{TRDCalibParams::Instance()}; ///< reference to calibration parameters + size_t mMinEntriesTotal{mParams.minEntriesTotalT0Fit}; ///< minimum total number of angular deviations + size_t mMinEntriesChamber{mParams.minEntriesChamberT0Fit}; ///< minimum number of angular deviations per chamber for accepting refitted value + bool mEnableOutput{false}; ///< enable output in a root file instead of the ccdb + std::unique_ptr mOutFile{nullptr}; ///< output file + std::unique_ptr mOutTree{nullptr}; ///< output tree + ErfLandauChi2Functor mFitFunctor; ///< used for minimization process, provides chi2 estimate + ROOT::Fit::Fitter mFitter; ///< instance of the ROOT fitter + std::array mParamsStart; ///< Starting parameters for fit + std::unique_ptr mFuncErfLandau; ///< helper function to calculate the t0 value after the fitting procedure + std::array t0_chambers; ///< t0 values of the individual chambers + float t0_average{-5}; ///< average t0 value across all chambers + + std::vector mInfoVector; ///< vector of CCDB infos; each element is filled with CCDB description of accompanying CCDB calibration object + std::vector mObjectVector; ///< vector of CCDB calibration objects; the extracted t0 per chamber and average for given slot + + std::unique_ptr adcProfIncl; ///< profile that holds inclusive PH spectrum + std::array, o2::trd::constants::MAXCHAMBER> adcProfDet; ///< array of profiles for PH spectrum of each chamber + + ClassDefNV(T0Fit, 1); +}; + +} // namespace trd +} // namespace o2 + +#endif // O2_TRD_T0FIT_H diff --git a/Detectors/TRD/calibration/src/T0Fit.cxx b/Detectors/TRD/calibration/src/T0Fit.cxx new file mode 100644 index 0000000000000..17030842bd981 --- /dev/null +++ b/Detectors/TRD/calibration/src/T0Fit.cxx @@ -0,0 +1,270 @@ +// 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 T0Fit.cxx +/// \brief Fits the TRD PH distributions to extract the t0 value +/// \author Luisa Bergmann + +#include "TRDCalibration/T0Fit.h" +#include "Framework/ProcessingContext.h" +#include "Framework/TimingInfo.h" +#include "Framework/InputRecord.h" +#include "DataFormatsTRD/Constants.h" +#include "TRDBase/GeometryBase.h" +#include "TStopwatch.h" +#include "CCDB/CcdbApi.h" +#include "CCDB/BasicCCDBManager.h" +#include "CommonUtils/NameConf.h" +#include "CommonUtils/MemFileHelper.h" +#include +#include + +#include "TMath.h" +#include "TCanvas.h" + +#include +#include +#include +#include + +using namespace o2::trd::constants; + +namespace o2::trd +{ +//______________________________________________________________________________________________ +double ErfLandauChi2Functor::operator()(const double* par) const +{ + // provides chi2 estimate of PH profile comparing the y-value of profiles to + // par[0]*ROOT::Math::erf(x) + par[1]*TMath::Landau(x,par[2],par[3]) + // + // par[0] : offset + // par[1] : amplitude + // par[2] : mean + // par[3] : sigma + + double sum2 = 0; + + for (int i = lowerBoundFit; i <= upperBoundFit; ++i) { + + if (TMath::Abs(y[i]) < 1e-7) { + continue; + } + + double funcVal = par[0] * TMath::Erf(x[i]) + par[1] * TMath::Landau(x[i], par[2], par[3]); + + sum2 += TMath::Power(y[i] - funcVal, 2) / y[i]; + // sum2 += TMath::Power((y[i] - funcVal)/yErr[i],2); + } + + return sum2; +} + +//______________________________________________________________________________________________ +using Slot = o2::calibration::TimeSlot; + +void T0Fit::initOutput() +{ + // reset the CCDB output vectors + mInfoVector.clear(); + mObjectVector.clear(); +} + +void T0Fit::initProcessing() +{ + if (mInitDone) { + return; + } + + adcProfIncl = std::make_unique("adcProfIncl", "adcProfIncl", 30, -0.5, 29.5); + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + adcProfDet[iDet] = std::make_unique(Form("adcProfDet_%i", iDet), Form("adcProfDet_%i", iDet), 30, -0.5, 29.5); + } + + mFitFunctor.lowerBoundFit = 0; + mFitFunctor.upperBoundFit = 5; + + double mParamsStart[4]; + mParamsStart[0] = 100; + mParamsStart[1] = 500; + mParamsStart[2] = 0.5; + mParamsStart[3] = 0.5; + + mFitter.SetFCN(4, mFitFunctor, mParamsStart); + + ROOT::Math::MinimizerOptions opt; + opt.SetMinimizerType("Minuit2"); + opt.SetMinimizerAlgorithm("Migrad"); + opt.SetPrintLevel(0); + opt.SetMaxFunctionCalls(1000); + opt.SetTolerance(.001); + mFitter.Config().SetMinimizerOptions(opt); + + mFuncErfLandau = std::make_unique( + "fErfLandau", [&](double* x, double* par) { return par[0] * TMath::Erf(x[0]) + par[1] * TMath::Landau(x[0], par[2], par[3]); }, mFitFunctor.lowerBoundFit, mFitFunctor.upperBoundFit, 4); + + // set tree addresses + if (mEnableOutput) { + mOutTree->Branch("t0_chambers", &t0_chambers); + mOutTree->Branch("t0_average", &t0_average); + } + + mInitDone = true; +} + +void T0Fit::retrievePrev(o2::framework::ProcessingContext& pc) +{ + // We either get a pointer to a valid object from the last ~hour or to the default object + // which is always present. The first has precedence over the latter. + auto dataCalT0 = pc.inputs().get("calt0"); + std::string msg = "Default Object"; + // We check if the object we got is the default one by comparing it to the defaults. + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + if (dataCalT0->getT0(iDet) != constants::T0DEFAULT) { + msg = "Previous Object"; + break; + } + } + LOG(info) << "FitInstance: From CCDB retrieved " << msg; + + // Here we set each entry regardless if it is the default or not. + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + t0_chambers[iDet] = dataCalT0->getT0(iDet); + } + t0_average = dataCalT0->getT0av(); +} + +void T0Fit::finalizeSlot(Slot& slot) +{ + // do actual fits for the data provided in the given slot + + TStopwatch timer; + timer.Start(); + initProcessing(); + + // get data and fill profiles + auto dataPH = slot.getContainer(); + + for (int iEntry = 0; iEntry < dataPH->getNEntries(); ++iEntry) { + int iDet = dataPH->getDetector(iEntry); + adcProfIncl->Fill(dataPH->getTimeBin(iEntry), dataPH->getADC(iEntry)); + adcProfDet[iDet]->Fill(dataPH->getTimeBin(iEntry), dataPH->getADC(iEntry)); + } + + // do fits + // inclusive distribution + mFitFunctor.x.clear(); + mFitFunctor.y.clear(); + for (int i = 1; i <= 30; ++i) { + mFitFunctor.x.push_back(i - 1); + mFitFunctor.y.push_back(adcProfIncl->GetBinContent(i)); + } + + mFitter.FitFCN(); + auto fitResult = mFitter.Result(); + + if (fitResult.IsValid()) { + mFuncErfLandau->SetParameters(fitResult.GetParams()[0], fitResult.GetParams()[1], fitResult.GetParams()[2], fitResult.GetParams()[3]); + t0_average = mFuncErfLandau->GetMaximumX(); + } else { + LOG(info) << "t0 fit for inclusive distribtion is not valid, set to -5"; + t0_average = -5; + } + + // single chambers + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + if (adcProfDet[iDet]->GetEntries() < mMinEntriesChamber) { + LOG(info) << "not enough entries in chamber " << iDet << " for t0 fit, set to -5"; + t0_chambers[iDet] = -5; + continue; + } + + mFitFunctor.x.clear(); + mFitFunctor.y.clear(); + for (int i = 1; i <= 30; ++i) { + mFitFunctor.x.push_back(i - 1); + mFitFunctor.y.push_back(adcProfDet[iDet]->GetBinContent(i)); + } + + mFitter.FitFCN(); + fitResult = mFitter.Result(); + + if (fitResult.IsValid()) { + mFuncErfLandau->SetParameters(fitResult.GetParams()[0], fitResult.GetParams()[1], fitResult.GetParams()[2], fitResult.GetParams()[3]); + t0_chambers[iDet] = mFuncErfLandau->GetMaximumX(); + } else { + LOG(info) << "t0 fit for chamber " << iDet << " is not valid, set to -5"; + t0_chambers[iDet] = -5; + } + } + + // fill tree + if (mEnableOutput) { + mOutTree->Fill(); + + LOGF(debug, "Fit result for inclusive distribution: t0 = ", t0_average); + for (int iDet = 0; iDet < MAXCHAMBER; ++iDet) { + LOGF(debug, "Fit result for chamber %i: t0 = ", iDet, t0_chambers[iDet]); + } + } + + // assemble CCDB object + CalT0 t0Object; + t0Object.setT0av(t0_average); + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + t0Object.setT0(iDet, t0_chambers[iDet]); + } + auto clName = o2::utils::MemFileHelper::getClassName(t0Object); + auto flName = o2::ccdb::CcdbApi::generateFileName(clName); + std::map metadata; // TODO: do we want to store any meta data? + long startValidity = slot.getStartTimeMS(); + mInfoVector.emplace_back("TRD/Calib/CalT0", clName, flName, metadata, startValidity, startValidity + o2::ccdb::CcdbObjectInfo::HOUR); + mObjectVector.push_back(t0Object); +} + +Slot& T0Fit::emplaceNewSlot(bool front, TFType tStart, TFType tEnd) +{ + auto& container = getSlots(); + auto& slot = front ? container.emplace_front(tStart, tEnd) : container.emplace_back(tStart, tEnd); + slot.setContainer(std::make_unique()); + return slot; +} + +void T0Fit::createOutputFile() +{ + mEnableOutput = true; + mOutFile = std::make_unique("trd_calt0.root", "RECREATE"); + if (mOutFile->IsZombie()) { + LOG(error) << "Failed to create output file!"; + mEnableOutput = false; + return; + } + mOutTree = std::make_unique("calib", "t0 values"); + LOG(info) << "Created output file trd_calt0.root"; +} + +void T0Fit::closeOutputFile() +{ + if (!mEnableOutput) { + return; + } + + try { + mOutFile->cd(); + mOutTree->Write(); + mOutTree.reset(); + mOutFile->Close(); + mOutFile.reset(); + } catch (std::exception const& e) { + LOG(error) << "Failed to write t0-value data file, reason: " << e.what(); + } + mEnableOutput = false; +} +} // namespace o2::trd diff --git a/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h b/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h index 52556ab90cebf..26a808fd15860 100644 --- a/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h +++ b/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h @@ -20,8 +20,11 @@ #pragma link C++ class o2::calibration::TimeSlotCalibration < o2::trd::AngularResidHistos> + ; #pragma link C++ class o2::calibration::TimeSlot < o2::trd::GainCalibHistos> + ; #pragma link C++ class o2::calibration::TimeSlotCalibration < o2::trd::GainCalibHistos> + ; +#pragma link C++ class o2::calibration::TimeSlot < o2::trd::T0FitHistos> + ; +#pragma link C++ class o2::calibration::TimeSlotCalibration < o2::trd::T0FitHistos> + ; #pragma link C++ class o2::trd::CalibratorNoise + ; #pragma link C++ class o2::trd::CalibratorGain + ; +#pragma link C++ class o2::trd::T0Fit + ; #pragma link C++ class o2::trd::ChannelInfoDetailed + ; #pragma link C++ class o2::trd::TrackBasedCalib + ; #pragma link C++ class o2::trd::PadCalibCCDBBuilder + ; diff --git a/Detectors/TRD/workflow/CMakeLists.txt b/Detectors/TRD/workflow/CMakeLists.txt index f71995a8fed00..069e5f11dc00e 100644 --- a/Detectors/TRD/workflow/CMakeLists.txt +++ b/Detectors/TRD/workflow/CMakeLists.txt @@ -27,6 +27,7 @@ o2_add_library(TRDWorkflow include/TRDWorkflow/TRDPulseHeightSpec.h include/TRDWorkflow/TRDGlobalTrackingQCSpec.h include/TRDWorkflow/NoiseCalibSpec.h + include/TRDWorkflow/T0FitSpec.h PUBLIC_LINK_LIBRARIES O2::Framework O2::DPLUtils O2::Steer O2::Algorithm diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/T0FitSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/T0FitSpec.h new file mode 100644 index 0000000000000..decd9b8220f28 --- /dev/null +++ b/Detectors/TRD/workflow/include/TRDWorkflow/T0FitSpec.h @@ -0,0 +1,177 @@ +// 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. + +#ifndef O2_TRD_T0FITSPEC_H +#define O2_TRD_T0FITSPEC_H + +/// \file T0FitSpec.h +/// \brief DPL device for steering the TRD t0 fits +/// \author Luisa Bergmann + +#include "TRDCalibration/T0Fit.h" +#include "DetectorsCalibration/Utils.h" +#include "DataFormatsTRD/PHData.h" +#include "DataFormatsTRD/CalT0.h" +#include "CommonUtils/MemFileHelper.h" +#include "Framework/Task.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/WorkflowSpec.h" +#include "CCDB/CcdbApi.h" +#include "CCDB/CcdbObjectInfo.h" +#include "Framework/CCDBParamSpec.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include + +#include "TH2.h" +#include "TH1.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace calibration +{ + +class T0FitDevice : public o2::framework::Task +{ + public: + T0FitDevice(std::shared_ptr req) : mCCDBRequest(req) {} + void init(o2::framework::InitContext& ic) final + { + o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest); + auto slotL = ic.options().get("sec-per-slot"); + auto delay = ic.options().get("max-delay"); + + mFitInstance = std::make_unique(); + mFitInstance->setSlotLengthInSeconds(slotL); + mFitInstance->setMaxSlotsDelay(delay); + if (ic.options().get("enable-root-output")) { + mFitInstance->createOutputFile(); + } + } + + void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final + { + o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj); + } + + void run(o2::framework::ProcessingContext& pc) final + { + const auto& tinfo = pc.services().get(); + if (tinfo.globalRunNumberChanged) { // new run is starting + mRunStopRequested = false; + mFitInstance->retrievePrev(pc); // SOR initialization is performed here + } + if (mRunStopRequested) { + return; + } + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + + auto dataT0Fit = pc.inputs().get>("input"); + o2::base::TFIDInfoHelper::fillTFIDInfo(pc, mFitInstance->getCurrentTFInfo()); + LOG(detail) << "Processing TF " << mFitInstance->getCurrentTFInfo().tfCounter << " with " << dataT0Fit.size() << " PHData entries"; + mFitInstance->process(dataT0Fit); + + if (pc.transitionState() == TransitionHandlingState::Requested) { + LOG(info) << "Run stop requested, finalizing"; + mRunStopRequested = true; + mFitInstance->checkSlotsToFinalize(o2::calibration::INFINITE_TF); + mFitInstance->closeOutputFile(); + } + sendOutput(pc.outputs()); + } + + void endOfStream(o2::framework::EndOfStreamContext& ec) final + { + if (mRunStopRequested) { + return; + } + mFitInstance->checkSlotsToFinalize(o2::calibration::INFINITE_TF); + mFitInstance->closeOutputFile(); + sendOutput(ec.outputs()); + } + + void stop() final + { + mFitInstance->closeOutputFile(); + } + + private: + std::unique_ptr mFitInstance; + std::shared_ptr mCCDBRequest; + bool mRunStopRequested = false; // flag that run was stopped (and the last output is sent) + //________________________________________________________________ + void sendOutput(DataAllocator& output) + { + // extract CCDB infos and calibration objects, convert it to TMemFile and send them to the output + + using clbUtils = o2::calibration::Utils; + const auto& payloadVec = mFitInstance->getCcdbObjectVector(); + auto& infoVec = mFitInstance->getCcdbObjectInfoVector(); // use non-const version as we update it + + assert(payloadVec.size() == infoVec.size()); + + for (uint32_t i = 0; i < payloadVec.size(); i++) { + auto& w = infoVec[i]; + auto image = o2::ccdb::CcdbApi::createObjectImage(&payloadVec[i], &w); + LOG(info) << "Sending object " << w.getPath() << "/" << w.getFileName() << " of size " << image->size() + << " bytes, valid for " << w.getStartValidityTimestamp() << " : " << w.getEndValidityTimestamp(); + + output.snapshot(Output{clbUtils::gDataOriginCDBPayload, "CALT0", i}, *image.get()); // vector + output.snapshot(Output{clbUtils::gDataOriginCDBWrapper, "CALT0", i}, w); // root-serialized + } + if (payloadVec.size()) { + mFitInstance->initOutput(); // reset the outputs once they are already sent + } + } +}; + +} // namespace calibration + +namespace framework +{ + +DataProcessorSpec getTRDT0FitSpec() +{ + using device = o2::calibration::T0FitDevice; + using clbUtils = o2::calibration::Utils; + + std::vector outputs; + outputs.emplace_back(ConcreteDataTypeMatcher{o2::calibration::Utils::gDataOriginCDBPayload, "CALT0"}, Lifetime::Sporadic); + outputs.emplace_back(ConcreteDataTypeMatcher{o2::calibration::Utils::gDataOriginCDBWrapper, "CALT0"}, Lifetime::Sporadic); + std::vector inputs; + inputs.emplace_back("input", "TRD", "PULSEHEIGHT"); + inputs.emplace_back("calt0", "TRD", "CALT0", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/CalT0")); + auto ccdbRequest = std::make_shared(true, // orbitResetTime + true, // GRPECS=true + false, // GRPLHCIF + false, // GRPMagField + false, // askMatLUT + o2::base::GRPGeomRequest::None, // geometry + inputs); + + return DataProcessorSpec{ + "calib-t0-fit", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(ccdbRequest)}, + Options{ + {"sec-per-slot", VariantType::UInt32, 900u, {"number of seconds per calibration time slot"}}, + {"max-delay", VariantType::UInt32, 2u, {"number of slots in past to consider"}}, + {"enable-root-output", VariantType::Bool, false, {"output t0 values to root file"}}, + }}; +} + +} // namespace framework +} // namespace o2 + +#endif // O2_TRD_T0FITSPEC_H diff --git a/Detectors/TRD/workflow/io/CMakeLists.txt b/Detectors/TRD/workflow/io/CMakeLists.txt index dfdae1edcc706..e27dba98b477b 100644 --- a/Detectors/TRD/workflow/io/CMakeLists.txt +++ b/Detectors/TRD/workflow/io/CMakeLists.txt @@ -20,6 +20,7 @@ o2_add_library(TRDWorkflowIO src/TRDTrackWriterSpec.cxx src/TRDTrackReaderSpec.cxx src/TRDCalibWriterSpec.cxx + src/TRDPHReaderSpec.cxx src/KrClusterWriterSpec.cxx PUBLIC_LINK_LIBRARIES O2::DataFormatsTRD O2::SimulationDataFormat O2::DPLUtils O2::GPUDataTypeHeaders O2::DataFormatsTPC) diff --git a/Detectors/TRD/workflow/io/include/TRDWorkflowIO/TRDPHReaderSpec.h b/Detectors/TRD/workflow/io/include/TRDWorkflowIO/TRDPHReaderSpec.h new file mode 100644 index 0000000000000..5e71f42f9fb19 --- /dev/null +++ b/Detectors/TRD/workflow/io/include/TRDWorkflowIO/TRDPHReaderSpec.h @@ -0,0 +1,53 @@ +// 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 TRDPHReaderSpec.h + +#ifndef O2_TRD_PHREADER +#define O2_TRD_PHREADER + +#include "TFile.h" +#include "TTree.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "DataFormatsTRD/PHData.h" + +namespace o2 +{ +namespace trd +{ + +class TRDPHReader : public o2::framework::Task +{ + public: + TRDPHReader() = default; + ~TRDPHReader() override = default; + void init(o2::framework::InitContext& ic) final; + void run(o2::framework::ProcessingContext& pc) final; + + private: + void connectTree(); + std::unique_ptr mFile; + std::unique_ptr mTree; + std::string mInFileName{"trd_PH.root"}; + std::string mInTreeName{"ph"}; + std::vector mPHValues, *mPHValuesPtr = &mPHValues; ///< to be used for branch address +}; + +/// create a processor spec +/// read TRD calibration data from a root file +framework::DataProcessorSpec getTRDPHReaderSpec(); + +} // namespace trd +} // namespace o2 + +#endif /* O2_TRD_PHREADER */ diff --git a/Detectors/TRD/workflow/io/src/TRDPHReaderSpec.cxx b/Detectors/TRD/workflow/io/src/TRDPHReaderSpec.cxx new file mode 100644 index 0000000000000..235c58c2ad302 --- /dev/null +++ b/Detectors/TRD/workflow/io/src/TRDPHReaderSpec.cxx @@ -0,0 +1,83 @@ +// 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 TRDPHReaderSpec.cxx + +#include "TRDWorkflowIO/TRDPHReaderSpec.h" + +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "fairlogger/Logger.h" +#include "CommonUtils/NameConf.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace trd +{ + +void TRDPHReader::init(InitContext& ic) +{ + // get the option from the init context + LOG(info) << "Init TRD PH reader!"; + mInFileName = o2::utils::Str::concat_string(o2::utils::Str::rectifyDirectory(ic.options().get("input-dir")), + ic.options().get("trd-ph-infile")); + mInTreeName = o2::utils::Str::concat_string(o2::utils::Str::rectifyDirectory(ic.options().get("input-dir")), + ic.options().get("treename")); + connectTree(); +} + +void TRDPHReader::connectTree() +{ + mTree.reset(nullptr); // in case it was already loaded + mFile.reset(TFile::Open(mInFileName.c_str())); + assert(mFile && !mFile->IsZombie()); + mTree.reset((TTree*)mFile->Get(mInTreeName.c_str())); + assert(mTree); + mTree->SetBranchAddress("values", &mPHValuesPtr); + LOG(info) << "Loaded tree from " << mInFileName << " with " << mTree->GetEntries() << " entries"; +} + +void TRDPHReader::run(ProcessingContext& pc) +{ + auto currEntry = mTree->GetReadEntry() + 1; + assert(currEntry < mTree->GetEntries()); + mTree->GetEntry(currEntry); + + LOG(info) << "Pushing vector of PH values filled with " << mPHValues.size() << " entries at tree entry " << currEntry; + pc.outputs().snapshot(Output{o2::header::gDataOriginTRD, "PULSEHEIGHT", 0, Lifetime::Timeframe}, mPHValues); + + if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } +} + +DataProcessorSpec getTRDPHReaderSpec() +{ + std::vector outputs; + outputs.emplace_back(o2::header::gDataOriginTRD, "PULSEHEIGHT", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "TRDPHReader", + Inputs{}, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{ + {"trd-ph-infile", VariantType::String, "trd_PH.root", {"Name of the input file"}}, + {"input-dir", VariantType::String, "none", {"Input directory"}}, + {"treename", VariantType::String, "ph", {"Name of top-level TTree"}}, + }}; +} + +} // namespace trd +} // namespace o2 diff --git a/Detectors/TRD/workflow/src/trd-calib-workflow.cxx b/Detectors/TRD/workflow/src/trd-calib-workflow.cxx index f3eb520873dbd..6a119d59272d7 100644 --- a/Detectors/TRD/workflow/src/trd-calib-workflow.cxx +++ b/Detectors/TRD/workflow/src/trd-calib-workflow.cxx @@ -12,9 +12,11 @@ #include "Framework/DataProcessorSpec.h" #include "TRDWorkflowIO/TRDCalibReaderSpec.h" #include "TRDWorkflowIO/TRDDigitReaderSpec.h" +#include "TRDWorkflowIO/TRDPHReaderSpec.h" #include "TRDWorkflow/VdAndExBCalibSpec.h" #include "TRDWorkflow/GainCalibSpec.h" #include "TRDWorkflow/NoiseCalibSpec.h" +#include "TRDWorkflow/T0FitSpec.h" #include "CommonUtils/ConfigurableParam.h" using namespace o2::framework; @@ -28,6 +30,7 @@ void customize(std::vector& workflowOptions) {"vDriftAndExB", o2::framework::VariantType::Bool, false, {"enable vDrift and ExB calibration"}}, {"noise", o2::framework::VariantType::Bool, false, {"enable noise and pad status calibration"}}, {"gain", o2::framework::VariantType::Bool, false, {"enable gain calibration"}}, + {"t0", o2::framework::VariantType::Bool, false, {"enable t0 fit"}}, {"calib-dds-collection-index", VariantType::Int, -1, {"allow only single collection to produce calibration objects (use -1 for no limit)"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}}; @@ -77,5 +80,12 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(getTRDGainCalibSpec()); } + if (configcontext.options().get("t0")) { + if (enableRootInp) { + specs.emplace_back(o2::trd::getTRDPHReaderSpec()); + } + specs.emplace_back(getTRDT0FitSpec()); + } + return specs; }