diff --git a/DataFormats/Detectors/ZDC/CMakeLists.txt b/DataFormats/Detectors/ZDC/CMakeLists.txt index 59bbb6be124cc..deda5de9cd739 100644 --- a/DataFormats/Detectors/ZDC/CMakeLists.txt +++ b/DataFormats/Detectors/ZDC/CMakeLists.txt @@ -12,8 +12,8 @@ o2_add_library(DataFormatsZDC SOURCES src/ChannelData.cxx src/BCData.cxx src/BCRecData.cxx src/RecEvent.cxx src/RecEventAux.cxx src/RawEventData.cxx src/OrbitRawData.cxx src/OrbitRecData.cxx src/OrbitData.cxx src/ZDCTDCData.cxx src/ZDCEnergy.cxx - src/CTF.cxx src/RecEventFlat.cxx - PUBLIC_LINK_LIBRARIES O2::CommonConstants O2::CommonDataFormat + src/CTF.cxx src/RecEventFlat.cxx src/InterCalibData.cxx + PUBLIC_LINK_LIBRARIES O2::CommonConstants O2::CommonDataFormat O2::DetectorsCalibration O2::ZDCBase ROOT::MathCore FairRoot::Base O2::SimulationDataFormat O2::MathUtils Microsoft.GSL::GSL) @@ -25,4 +25,5 @@ o2_target_root_dictionary(DataFormatsZDC include/DataFormatsZDC/RecEvent.h include/DataFormatsZDC/RecEventAux.h include/DataFormatsZDC/RecEventFlat.h include/DataFormatsZDC/OrbitRawData.h include/DataFormatsZDC/ZDCTDCData.h include/DataFormatsZDC/BCRecData.h include/DataFormatsZDC/ZDCEnergy.h - include/DataFormatsZDC/OrbitRecData.h include/DataFormatsZDC/RawEventData.h) + include/DataFormatsZDC/OrbitRecData.h include/DataFormatsZDC/RawEventData.h + include/DataFormatsZDC/InterCalibData.h) diff --git a/DataFormats/Detectors/ZDC/include/DataFormatsZDC/InterCalibData.h b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/InterCalibData.h new file mode 100644 index 0000000000000..82279ad855572 --- /dev/null +++ b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/InterCalibData.h @@ -0,0 +1,45 @@ +// 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 _ZDC_INTERCALIB_DATA_H_ +#define _ZDC_INTERCALIB_DATA_H_ + +#include "ZDCBase/Constants.h" +#include +#include + +/// \file InterCalibData.h +/// \brief Intercalibration intermediate data +/// \author pietro.cortese@cern.ch + +namespace o2 +{ +namespace zdc +{ + +struct InterCalibData { + static constexpr int NPAR = 6; /// Dimension of matrix (1 + 4 coefficients + offset) + static constexpr int NH = 5; /// ZNA, ZPA, ZNC, ZPC, ZEM + double mSum[NH][NPAR][NPAR] = {0}; /// Cumulated sums + uint64_t mCTimeBeg = 0; /// Time of processed time frame + uint64_t mCTimeEnd = 0; /// Time of processed time frame + static constexpr const char* DN[NH] = {"ZNA", "ZPA", "ZNC", "ZPC", "ZEM"}; + InterCalibData& operator+=(const InterCalibData& other); + int getEntries(int ih) const; + void print() const; + void setCreationTime(uint64_t ctime); + ClassDefNV(InterCalibData, 1); +}; + +} // namespace zdc +} // namespace o2 + +#endif diff --git a/DataFormats/Detectors/ZDC/include/DataFormatsZDC/RecEventFlat.h b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/RecEventFlat.h index 927b1e3b2d8dd..11dd493ad722a 100644 --- a/DataFormats/Detectors/ZDC/include/DataFormatsZDC/RecEventFlat.h +++ b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/RecEventFlat.h @@ -20,6 +20,7 @@ #include "ZDCBase/Constants.h" #include "MathUtils/Cartesian.h" #include +#include #include #include #include @@ -37,31 +38,32 @@ using NElem = int; struct RecEventFlat { // NOLINT: false positive in clang-tidy !! o2::InteractionRecord ir; - uint32_t channels = 0; /// pattern of channels acquired - uint32_t triggers = 0; /// pattern of channels with autotrigger bit - std::map ezdc; /// signal in ZDCs - std::vector TDCVal[NTDCChannels]; /// TDC values - std::vector TDCAmp[NTDCChannels]; /// TDC signal amplitudes - std::vector TDCPile[NTDCChannels]; /// TDC pile-up correction flag (TODO) - std::vector* mRecBC; //! Interaction record and references to data - std::vector* mEnergy; //! ZDC energy - std::vector* mTDCData; //! ZDC TDC - std::vector* mInfo; //! Event quality information - std::vector mDecodedInfo; //! Event quality information (decoded) - uint64_t mEntry = 0; //! Current entry - uint64_t mNEntries = 0; //! Number of entries - FirstEntry mFirstE = 0; //! First energy - FirstEntry mFirstT = 0; //! First TDC - FirstEntry mFirstI = 0; //! First info - FirstEntry mStopE = 0; //! Last + 1 energy - FirstEntry mStopT = 0; //! Last + 1 TDC - FirstEntry mStopI = 0; //! Last + 1 info - NElem mNE = 0; //! N energy - NElem mNT = 0; //! N TDC - NElem mNI = 0; //! N info - std::array isBeg{}; //! Beginning of sequence - std::array isEnd{}; //! End of sequence - o2::zdc::BCRecData mCurB; //! Current BC + uint32_t channels = 0; /// pattern of channels acquired + uint32_t ezdcDecoded = 0; /// pattern of decoded energies + uint32_t triggers = 0; /// pattern of channels with autotrigger bit + std::map ezdc; /// signal in ZDCs + std::vector TDCVal[NTDCChannels]; /// TDC values + std::vector TDCAmp[NTDCChannels]; /// TDC signal amplitudes + std::vector TDCPile[NTDCChannels]; /// TDC pile-up correction flag (TODO) + gsl::span mRecBC; //! Interaction record and references to data + gsl::span mEnergy; //! ZDC energy + gsl::span mTDCData; //! ZDC TDC + gsl::span mInfo; //! Event quality information + std::vector mDecodedInfo; //! Event quality information (decoded) + uint64_t mEntry = 0; //! Current entry + uint64_t mNEntries = 0; //! Number of entries + FirstEntry mFirstE = 0; //! First energy + FirstEntry mFirstT = 0; //! First TDC + FirstEntry mFirstI = 0; //! First info + FirstEntry mStopE = 0; //! Last + 1 energy + FirstEntry mStopT = 0; //! Last + 1 TDC + FirstEntry mStopI = 0; //! Last + 1 info + NElem mNE = 0; //! N energy + NElem mNT = 0; //! N TDC + NElem mNI = 0; //! N info + std::array isBeg{}; //! Beginning of sequence + std::array isEnd{}; //! End of sequence + o2::zdc::BCRecData mCurB; //! Current BC // Reconstruction messages std::array genericE{}; /// 0 Generic error @@ -93,7 +95,8 @@ struct RecEventFlat { // NOLINT: false positive in clang-tidy !! uint8_t mVerbosity = DbgZero; //! Verbosity level uint32_t mTriggerMask = 0; //! Trigger mask for printout - void init(std::vector* RecBC, std::vector* Energy, std::vector* TDCData, std::vector* Info); + void init(const std::vector* RecBC, const std::vector* Energy, const std::vector* TDCData, const std::vector* Info); + void init(const gsl::span RecBC, const gsl::span Energy, const gsl::span TDCData, const gsl::span Info); int next(); @@ -102,6 +105,18 @@ struct RecEventFlat { // NOLINT: false positive in clang-tidy !! return mNE; } + inline bool getEnergy(int32_t i, uint8_t& key, float& val) const + { + if (i < mNE) { + auto it = ezdc.begin(); + std::advance(it, i); + key = it->first; + val = it->second; + return true; + } + return false; + } + inline NElem getNTDC() const { return mNT; diff --git a/DataFormats/Detectors/ZDC/include/DataFormatsZDC/ZDCTDCData.h b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/ZDCTDCData.h index 3c85c17b8da2f..b159ccdb1a3a2 100644 --- a/DataFormats/Detectors/ZDC/include/DataFormatsZDC/ZDCTDCData.h +++ b/DataFormats/Detectors/ZDC/include/DataFormatsZDC/ZDCTDCData.h @@ -36,6 +36,7 @@ struct ZDCTDCData { ZDCTDCData() = default; ZDCTDCData(uint8_t ida, int16_t vala, int16_t ampa, bool isbeg = false, bool isend = false) { + // TDC value and amplitude are encoded externally id = ida & 0x0f; id = id | (isbeg ? 0x80 : 0x00); id = id | (isend ? 0x40 : 0x00); @@ -45,6 +46,7 @@ struct ZDCTDCData { ZDCTDCData(uint8_t ida, float vala, float ampa, bool isbeg = false, bool isend = false) { + // TDC value and amplitude are encoded externally id = ida & 0x0f; id = id | (isbeg ? 0x80 : 0x00); id = id | (isend ? 0x40 : 0x00); @@ -75,10 +77,12 @@ struct ZDCTDCData { inline float amplitude() const { + // Return decoded value return FTDCAmp * amp; } inline float value() const { + // Return decoded value (ns) return FTDCVal * val; } inline int ch() const diff --git a/DataFormats/Detectors/ZDC/src/DataFormatsZDCLinkDef.h b/DataFormats/Detectors/ZDC/src/DataFormatsZDCLinkDef.h index 699234e92bf59..abe1d66311913 100644 --- a/DataFormats/Detectors/ZDC/src/DataFormatsZDCLinkDef.h +++ b/DataFormats/Detectors/ZDC/src/DataFormatsZDCLinkDef.h @@ -32,6 +32,7 @@ #pragma link C++ class o2::zdc::RecEventFlat + ; #pragma link C++ class o2::zdc::ZDCEnergy + ; #pragma link C++ class o2::zdc::ZDCTDCData + ; +#pragma link C++ class o2::zdc::InterCalibData + ; #pragma link C++ class std::vector < o2::zdc::ChannelData> + ; #pragma link C++ class std::vector < o2::zdc::BCData> + ; #pragma link C++ class std::vector < o2::zdc::OrbitData> + ; diff --git a/DataFormats/Detectors/ZDC/src/InterCalibData.cxx b/DataFormats/Detectors/ZDC/src/InterCalibData.cxx new file mode 100644 index 0000000000000..9504c64b93239 --- /dev/null +++ b/DataFormats/Detectors/ZDC/src/InterCalibData.cxx @@ -0,0 +1,72 @@ +// 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. + +#include "Framework/Logger.h" +#include "DataFormatsZDC/InterCalibData.h" + +using namespace o2::zdc; + +void InterCalibData::print() const +{ + for (int i = 0; i < NH; i++) { + LOGF(info, "%s", DN[i]); + for (int j = 0; j < NPAR; j++) { + for (int k = 0; k < NPAR; k++) { + if (k == 0) { + printf("%e", mSum[i][j][k]); + } else { + printf(" %e", mSum[i][j][k]); + } + } + printf("\n"); + } + } +} + +InterCalibData& InterCalibData::operator+=(const InterCalibData& other) +{ + for (int32_t ih = 0; ih < NH; ih++) { + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + mSum[ih][i][j] += other.mSum[ih][i][j]; + } + } + } + if (mCTimeBeg == 0 || other.mCTimeBeg < mCTimeBeg) { + mCTimeBeg = other.mCTimeBeg; + } + if (other.mCTimeEnd > mCTimeEnd) { + mCTimeEnd = other.mCTimeEnd; + } + //#ifdef O2_ZDC_DEBUG + LOGF(info, "InterCalibData [%llu : %llu]: %s=%d %s=%d %s=%d %s=%d %s=%d", mCTimeBeg, mCTimeEnd, DN[0], getEntries(0), DN[1], getEntries(1), + DN[2], getEntries(2), DN[3], getEntries(3), DN[4], getEntries(4)); + //#endif + return *this; +} + +void InterCalibData::setCreationTime(uint64_t ctime) +{ + mCTimeBeg = ctime; + mCTimeEnd = ctime; +#ifdef O2_ZDC_DEBUG + LOGF(info, "InterCalibData::setCreationTime %llu", ctime); +#endif +} + +int InterCalibData::getEntries(int ih) const +{ + if (ih < 0 || ih >= NH) { + LOGF(error, "InterCalibData::getEntries ih = %d is out of range", ih); + return 0; + } + return mSum[ih][5][5]; +} diff --git a/DataFormats/Detectors/ZDC/src/RecEventFlat.cxx b/DataFormats/Detectors/ZDC/src/RecEventFlat.cxx index d84999711b824..440edff22f27a 100644 --- a/DataFormats/Detectors/ZDC/src/RecEventFlat.cxx +++ b/DataFormats/Detectors/ZDC/src/RecEventFlat.cxx @@ -15,14 +15,24 @@ using namespace o2::zdc; -void RecEventFlat::init(std::vector* RecBC, std::vector* Energy, std::vector* TDCData, std::vector* Info) +void RecEventFlat::init(const std::vector* RecBC, const std::vector* Energy, const std::vector* TDCData, const std::vector* Info) +{ + mRecBC = *RecBC; + mEnergy = *Energy; + mTDCData = *TDCData; + mInfo = *Info; + mEntry = 0; + mNEntries = mRecBC.size(); +} + +void RecEventFlat::init(const gsl::span RecBC, const gsl::span Energy, const gsl::span TDCData, const gsl::span Info) { mRecBC = RecBC; mEnergy = Energy; mTDCData = TDCData; mInfo = Info; mEntry = 0; - mNEntries = mRecBC->size(); + mNEntries = mRecBC.size(); } void RecEventFlat::clearBitmaps() @@ -55,6 +65,7 @@ void RecEventFlat::clearBitmaps() int RecEventFlat::next() { + ezdcDecoded = 0; if (mEntry >= mNEntries) { return 0; } @@ -68,7 +79,7 @@ int RecEventFlat::next() } // Get References - mCurB = mRecBC->at(mEntry); + mCurB = mRecBC[mEntry]; mCurB.getRef(mFirstE, mNE, mFirstT, mNT, mFirstI, mNI); mStopE = mFirstE + mNE; mStopT = mFirstT + mNT; @@ -84,7 +95,7 @@ int RecEventFlat::next() uint16_t code = 0; uint32_t map = 0; for (int i = mFirstI; i < mStopI; i++) { - uint16_t info = mInfo->at(i); + uint16_t info = mInfo[i]; if (infoState == 0) { if (info & 0x8000) { LOGF(error, "Inconsistent info stream at word %d: 0x%4u", i, info); @@ -122,18 +133,19 @@ int RecEventFlat::next() // Decode energy for (int i = mFirstE; i < mStopE; i++) { - auto myenergy = mEnergy->at(i); + auto myenergy = mEnergy[i]; auto ch = myenergy.ch(); ezdc[ch] = myenergy.energy(); // Assign implicit event info if (adcPedOr[ch] == false && adcPedQC[ch] == false && adcPedMissing[ch] == false) { adcPedEv[ch] = true; } + ezdcDecoded |= (0x1 << ch); } // Decode TDCs for (int i = mFirstT; i < mStopT; i++) { - auto mytdc = mTDCData->at(i); + auto mytdc = mTDCData[i]; auto ch = mytdc.ch(); if (ch < NTDCChannels) { if (mytdc.isBeg()) { diff --git a/Detectors/ZDC/CMakeLists.txt b/Detectors/ZDC/CMakeLists.txt index c803b3917984d..c9f51fc56c903 100644 --- a/Detectors/ZDC/CMakeLists.txt +++ b/Detectors/ZDC/CMakeLists.txt @@ -17,4 +17,5 @@ add_subdirectory(simulation) add_subdirectory(reconstruction) add_subdirectory(macro) add_subdirectory(raw) -add_subdirectory(workflow) \ No newline at end of file +add_subdirectory(workflow) +add_subdirectory(calib) diff --git a/Detectors/ZDC/base/include/ZDCBase/Constants.h b/Detectors/ZDC/base/include/ZDCBase/Constants.h index 077cc4e2212d7..015187f87d927 100644 --- a/Detectors/ZDC/base/include/ZDCBase/Constants.h +++ b/Detectors/ZDC/base/include/ZDCBase/Constants.h @@ -167,6 +167,16 @@ constexpr int IdZPC3 = 23; constexpr int IdZPC4 = 24; constexpr int IdZPCSum = 25; +constexpr uint32_t MaskZNA = 0x0000001f; +constexpr uint32_t MaskAllZNA = 0x0000003f; +constexpr uint32_t MaskZPA = 0x000007c0; +constexpr uint32_t MaskAllZPA = 0x00000fc0; +constexpr uint32_t MaskZEM = 0x00003000; +constexpr uint32_t MaskZNC = 0x000fc000; +constexpr uint32_t MaskAllZNC = 0x0007f000; +constexpr uint32_t MaskZPC = 0x01f00000; +constexpr uint32_t MaskAllZPC = 0x03f00000; + constexpr std::string_view ChannelNames[] = { "ZNAC", "ZNA1", @@ -229,6 +239,7 @@ const std::string CCDBPathTDCCalib = "ZDC/Calib/TDCCalib"; const std::string CCDBPathTDCCorr = "ZDC/Calib/TDCCorr"; const std::string CCDBPathEnergyCalib = "ZDC/Calib/EnergyCalib"; const std::string CCDBPathTowerCalib = "ZDC/Calib/TowerCalib"; +const std::string CCDBPathInterCalibConfig = "ZDC/Calib/InterCalibConfig"; enum Ped { PedND = 0, PedEv = 1, @@ -287,10 +298,11 @@ constexpr std::array ChEnergyCalib{IdZNAC, IdZNASum, IdZPAC, IdZPASum, IdZEM1, IdZEM2, IdZNCC, IdZNCSum, IdZPCC, IdZPCSum}; -constexpr std::array ChTowerCalib{IdZNA1, IdZNA2, IdZNA3, IdZNA4, +constexpr std::array ChTowerCalib{IdZNA1, IdZNA2, IdZNA3, IdZNA4, IdZPA1, IdZPA2, IdZPA3, IdZPA4, IdZNC1, IdZNC2, IdZNC3, IdZNC4, - IdZPC1, IdZPC2, IdZPC3, IdZPC4}; + IdZPC1, IdZPC2, IdZPC3, IdZPC4, + IdZEM2}; constexpr std::array CaloCommonPM{IdZNAC, IdZNAC, IdZNAC, IdZNAC, IdZNAC, IdZNAC, IdZPAC, IdZPAC, IdZPAC, IdZPAC, IdZPAC, IdZPAC, diff --git a/Detectors/ZDC/calib/CMakeLists.txt b/Detectors/ZDC/calib/CMakeLists.txt new file mode 100644 index 0000000000000..d163e440c7af8 --- /dev/null +++ b/Detectors/ZDC/calib/CMakeLists.txt @@ -0,0 +1,51 @@ +# 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. + +o2_add_library(ZDCCalib + SOURCES + src/InterCalib.cxx + src/InterCalibConfig.cxx + src/InterCalibSpec.cxx + src/InterCalibEPN.cxx + src/InterCalibEPNSpec.cxx + PUBLIC_LINK_LIBRARIES + O2::CCDB + O2::DPLUtils + O2::DataFormatsZDC + O2::DetectorsRaw + O2::Headers + O2::SimConfig + O2::SimulationDataFormat + O2::ZDCBase + O2::ZDCSimulation + O2::ZDCReconstruction + O2::DetectorsCalibration + ROOT::Core + ROOT::Hist + ROOT::Minuit) + +o2_target_root_dictionary(ZDCCalib + HEADERS + include/ZDCCalib/InterCalib.h + include/ZDCCalib/InterCalibConfig.h + include/ZDCCalib/InterCalibSpec.h + include/ZDCCalib/InterCalibEPN.h + include/ZDCCalib/InterCalibEPNSpec.h) + +o2_add_executable(intercalib-workflow + COMPONENT_NAME zdc + SOURCES src/zdc-intercalib-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::ZDCWorkflow O2::ZDCCalib O2::DetectorsCalibration) + +o2_add_executable(intercalib-epn-workflow + COMPONENT_NAME zdc + SOURCES src/zdc-intercalib-epn-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::ZDCWorkflow O2::ZDCCalib O2::DetectorsCalibration) diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalib.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalib.h new file mode 100644 index 0000000000000..b32711f9edc6d --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalib.h @@ -0,0 +1,110 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include "ZDCBase/Constants.h" +#include "CommonDataFormat/FlatHisto1D.h" +#include "CommonDataFormat/FlatHisto2D.h" +#include "DataFormatsZDC/RecEvent.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "ZDCCalib/InterCalibConfig.h" +#include "CCDB/CcdbObjectInfo.h" +#ifndef ALICEO2_ZDC_INTERCALIB_H_ +#define ALICEO2_ZDC_INTERCALIB_H_ +namespace o2 +{ +namespace zdc +{ +class InterCalib +{ + using CcdbObjectInfo = o2::ccdb::CcdbObjectInfo; + + public: + InterCalib() = default; + int init(); + static constexpr int HidZNA = 0; + static constexpr int HidZPA = 1; + static constexpr int HidZNC = 2; + static constexpr int HidZPC = 3; + static constexpr int HidZEM = 4; + static constexpr int NH = InterCalibData::NH; + static constexpr int NPAR = InterCalibData::NPAR; + void clear(int ih = -1); + int process(const gsl::span& bcrec, + const gsl::span& energy, + const gsl::span& tdc, + const gsl::span& info); // Calibration of RUN3 data - direct + int process(const InterCalibData& data); // Calibration of RUN3 data - aggregator node + int endOfRun(); // Perform minimization + int process(const char* hname, int ic); // Calibration of RUN2 data + void replay(int ih, THnSparse* hs, int ic); // Test of calibration using RUN2 data + int mini(int ih); + void add(int ih, o2::dataformats::FlatHisto1D& h1); + void add(int ih, o2::dataformats::FlatHisto2D& h2); + int write(const std::string fn = "ZDCInterCalib.root"); + + static double mAdd[NPAR][NPAR]; /// Temporary copy of cumulated sums + static void fcn(int& npar, double* gin, double& chi, double* par, int iflag); + void cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w); + + const ZDCTowerParam& getTowerParamUpd() const { return mTowerParamUpd; }; + CcdbObjectInfo& getCcdbObjectInfo() { return mInfo; } + + void setEnergyParam(const ZDCEnergyParam* param) { mEnergyParam = param; }; + const ZDCEnergyParam* getEnergyParam() const { return mEnergyParam; }; + void setTowerParam(const ZDCTowerParam* param) { mTowerParam = param; }; + const ZDCTowerParam* getTowerParam() const { return mTowerParam; }; + void setInterCalibConfig(const InterCalibConfig* param) { mInterCalibConfig = param; }; + const InterCalibConfig* getInterCalibConfig() const { return mInterCalibConfig; }; + + void setVerbosity(int v) { mVerbosity = v; } + int getVerbosity() const { return mVerbosity; } + + void setSaveDebugHistos() { mSaveDebugHistos = true; } + void setDontSaveDebugHistos() { mSaveDebugHistos = false; } + + static constexpr const char* mHUncN[2 * NH] = {"hZNAS", "hZPAS", "hZNCS", "hZPCS", "hZEM2", "hZNAC", "hZPAC", "hZNCC", "hZPCC", "hZEM1"}; + static constexpr const char* mHUncT[2 * NH] = {"ZNA sum", "ZPA sum", "ZNC sum", "ZPC sum", "ZEM2", "ZNA TC", "ZPA TC", "ZNC TC", "ZPC TC", "ZEM1"}; + static constexpr const char* mCUncN[NH] = {"cZNA", "cZPA", "cZNC", "cZPC", "cZEM"}; + static constexpr const char* mCUncT[NH] = {"ZNA;TC;SUM", "ZPA;TC;SUM", "ZNC;TC;SUM", "ZPC;TC;SUM", "ZEM;ZEM1;ZEM2"}; + + private: + std::array*, 2 * NH> mHUnc{}; + std::array*, NH> mCUnc{}; + std::array, NH> mHCorr{}; + std::array, NH> mCCorr{}; + std::array, NH> mMn{}; + InterCalibData mData; + bool mInitDone = false; + bool mSaveDebugHistos = false; + int32_t mVerbosity = DbgMinimal; + static std::mutex mMtx; /// mutex for critical section + double mPar[NH][NPAR] = {0}; + double mErr[NH][NPAR] = {0}; + const InterCalibConfig* mInterCalibConfig = nullptr; /// Configuration of intercalibration + const ZDCEnergyParam* mEnergyParam = nullptr; /// Energy calibration object + const ZDCTowerParam* mTowerParam = nullptr; /// Tower calibration object + ZDCTowerParam mTowerParamUpd; /// Updated tower calibration object + CcdbObjectInfo mInfo; /// CCDB Info + void assign(int ih, bool ismod); /// Assign updated calibration object +}; +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalibConfig.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibConfig.h new file mode 100644 index 0000000000000..5df736a5a815a --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibConfig.h @@ -0,0 +1,79 @@ +// 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_ZDC_INTERCALIBCONFIG_H +#define O2_ZDC_INTERCALIBCONFIG_H + +#include "ZDCBase/Constants.h" +#include +#include +#include + +/// \file InterCalibConfig.h +/// \brief Configuration of ZDC Tower intercalibration procedure +/// \author P. Cortese + +namespace o2 +{ +namespace zdc +{ +struct InterCalibConfig { + static constexpr int NH = 5; /// ZNA, ZPA, ZNC, ZPC, ZEM + bool enabled[NH] = {true, true, true, true, true}; + double cutLow[NH] = {-std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity()}; + double cutHigh[NH] = {std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + int nb1[NH] = {0}; /// 1D histogram: number of bins + double amin1[NH] = {0}; /// minimum + double amax1[NH] = {0}; /// maximum + int nb2[NH] = {0}; /// 2D histogram: number of bins + double amin2[NH] = {0}; /// minimum + double amax2[NH] = {0}; /// maximum + double l_bnd[NH] = {0.1, 0.1, 0.1, 0.1, 0.1}; + double u_bnd[NH] = {10., 10., 10., 10., 10.}; + double l_bnd_o[NH] = {-20., -20., -20., -20., -20.}; + double u_bnd_o[NH] = {20., 20., 20., 20., 20.}; + double step_o[NH] = {0., 0., 0., 0., 0.}; + double min_e[NH] = {0., 0., 0., 0., 0.}; + std::string desc = ""; + + void print() const; + void resetCuts(); + void resetCutLow(); + void resetCutHigh(); + void resetCutLow(int ih); + void resetCutHigh(int ih); + void setMinEntries(double val); + void setMinEntries(int ih, double val); + void setCutLow(double val); + void setCutHigh(double val); + void setCutLow(int ih, double val); + void setCutHigh(int ih, double val); + void setCuts(double low, double high); + void setCuts(int ih, double low, double high); + void setBinning1D(int nb, double amin, double amax); + void setBinning2D(int nb, double amin, double amax); + void setBinning1D(int ih, int nb, double amin, double amax); + void setBinning2D(int ih, int nb, double amin, double amax); + void setDescription(std::string d) { desc = d; } + void enable(bool c0, bool c1, bool c2, bool c3, bool c4) + { + enabled[0] = c0; + enabled[1] = c1; + enabled[2] = c2; + enabled[3] = c3; + enabled[4] = c4; + } + ClassDefNV(InterCalibConfig, 3); +}; +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPN.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPN.h new file mode 100644 index 0000000000000..01fbfd49239d3 --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPN.h @@ -0,0 +1,72 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include "ZDCBase/Constants.h" +#include "CommonDataFormat/FlatHisto1D.h" +#include "CommonDataFormat/FlatHisto2D.h" +#include "DataFormatsZDC/RecEvent.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "ZDCCalib/InterCalibConfig.h" +#ifndef ALICEO2_ZDC_INTERCALIBEPN_H_ +#define ALICEO2_ZDC_INTERCALIBEPN_H_ +namespace o2 +{ +namespace zdc +{ +class InterCalibEPN +{ + public: + InterCalibEPN() = default; + int init(); + static constexpr int HidZNA = 0; + static constexpr int HidZPA = 1; + static constexpr int HidZNC = 2; + static constexpr int HidZPC = 3; + static constexpr int HidZEM = 4; + static constexpr int NH = InterCalibData::NH; + static constexpr int NPAR = InterCalibData::NPAR; + void clear(int ih = -1); + int process(const gsl::span& bcrec, + const gsl::span& energy, + const gsl::span& tdc, + const gsl::span& info); // Calibration of RUN3 data + int endOfRun(); // Perform minimization + int process(const char* hname, int ic); // Calibration of RUN2 data + int write(const std::string fn = "ZDCInterCalibEPN.root"); + void cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w); + void setInterCalibConfig(const InterCalibConfig* param) { mInterCalibConfig = param; }; + const InterCalibConfig* getInterCalibConfig() const { return mInterCalibConfig; }; + void setSaveDebugHistos() { mSaveDebugHistos = true; } + void setDontSaveDebugHistos() { mSaveDebugHistos = false; } + InterCalibData& getData() { return mData; } + InterCalibData mData; + std::array*, 2 * NH> mH{}; + std::array*, NH> mC{}; + + private: + bool mInitDone = false; + bool mSaveDebugHistos = false; + int32_t mVerbosity = DbgMinimal; + const InterCalibConfig* mInterCalibConfig = nullptr; /// Configuration of intercalibration +}; +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPNSpec.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPNSpec.h new file mode 100644 index 0000000000000..d6e545092350a --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibEPNSpec.h @@ -0,0 +1,56 @@ +// 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 InterCalibEPNSpec.h +/// @brief ZDC intercalibration pre-processing on EPN +/// @author pietro.cortese@cern.ch + +#ifndef O2_ZDC_INTERCALIBEPN_SPEC +#define O2_ZDC_INTERCALIBEPN_SPEC + +#include +#include "Framework/Logger.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "CommonUtils/NameConf.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "ZDCCalib/InterCalibEPN.h" +#include "ZDCCalib/InterCalibConfig.h" + +namespace o2 +{ +namespace zdc +{ + +class InterCalibEPNSpec : public o2::framework::Task +{ + public: + InterCalibEPNSpec(); + InterCalibEPNSpec(const int verbosity); + ~InterCalibEPNSpec() override = default; + void init(o2::framework::InitContext& ic) final; + void updateTimeDependentParams(o2::framework::ProcessingContext& pc); + void run(o2::framework::ProcessingContext& pc) final; + void endOfStream(o2::framework::EndOfStreamContext& ec) final; + + private: + int mVerbosity = 0; // Verbosity level + bool mInitialized = false; // Connect once to CCDB during initialization + InterCalibEPN mInterCalibEPN; // Intercalibration object + TStopwatch mTimer; +}; + +framework::DataProcessorSpec getInterCalibEPNSpec(); + +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalibSpec.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibSpec.h new file mode 100644 index 0000000000000..e6551ea86b726 --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibSpec.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 InterCalibSpec.h +/// @brief ZDC intercalibration +/// @author pietro.cortese@cern.ch + +#ifndef O2_ZDC_INTERCALIB_SPEC +#define O2_ZDC_INTERCALIB_SPEC + +#include +#include "Framework/Logger.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/DataAllocator.h" +#include "Framework/Task.h" +#include "CommonDataFormat/FlatHisto1D.h" +#include "CommonDataFormat/FlatHisto2D.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "CommonUtils/NameConf.h" +#include "ZDCCalib/InterCalib.h" +#include "ZDCCalib/InterCalibConfig.h" +#include "DetectorsCalibration/Utils.h" +#include "CCDB/CcdbObjectInfo.h" + +namespace o2 +{ +namespace zdc +{ + +class InterCalibSpec : public o2::framework::Task +{ + public: + InterCalibSpec(); + InterCalibSpec(const int verbosity); + ~InterCalibSpec() override = default; + void init(o2::framework::InitContext& ic) final; + void updateTimeDependentParams(o2::framework::ProcessingContext& pc); + void run(o2::framework::ProcessingContext& pc) final; + void endOfStream(o2::framework::EndOfStreamContext& ec) final; + void sendOutput(o2::framework::DataAllocator& output); + + private: + int mVerbosity = DbgMinimal; // Verbosity level + InterCalib mInterCalib; // Intercalibration object + TStopwatch mTimer; +}; + +framework::DataProcessorSpec getInterCalibSpec(); + +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/InterCalibWorkflow.h b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibWorkflow.h new file mode 100644 index 0000000000000..6870ba306501d --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/InterCalibWorkflow.h @@ -0,0 +1,26 @@ +// 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_ZDC_INTERCALIBWORKFLOW_H +#define O2_ZDC_INTERCALIBWORKFLOW_H + +/// @file InterCalibWorkflow.h + +#include "Framework/WorkflowSpec.h" + +namespace o2 +{ +namespace zdc +{ +framework::WorkflowSpec getInterCalibWorkflow(const int verbosity); +} // namespace zdc +} // namespace o2 +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/TDCCalib.h b/Detectors/ZDC/calib/include/ZDCCalib/TDCCalib.h new file mode 100644 index 0000000000000..04f2d012c5b66 --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/TDCCalib.h @@ -0,0 +1,63 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include "ZDCBase/Constants.h" +#include "DataFormatsZDC/RecEvent.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCCalib/InterCalibConfig.h" +#ifndef ALICEO2_ZDC_TDCCALIB_H_ +#define ALICEO2_ZDC_TDCCALIB_H_ +namespace o2 +{ +namespace zdc +{ +class InterCalib +{ + public: + TDCCalib() = default; + int init(); + void clear(int ih = -1); + int process(const gsl::span& bcrec, + const gsl::span& energy, + const gsl::span& tdc, + const gsl::span& info); // Calibration of RUN3 data + int endOfRun(); // Perform minimization + int process(const char* hname, int ic); // Calibration of RUN2 data + void replay(int ih, THnSparse* hs, int ic); // Test of calibration using RUN2 data + int mini(int ih); + int write(const std::string fn = "ZDCTDCCalib.root"); + + static void fcn(int& npar, double* gin, double& chi, double* par, int iflag); + + void setTDCParam(const ZDCTDCParam* param) { mTDCParam = param; }; + const ZDCTDCParam* getTDCParam() const { return mTDCParam; }; + void setTDCCalibConfig(const TDCCalibConfig* param) { mTDCCalibConfig = param; }; + const TDCCalibConfig* getTDCCalibConfig() const { return mTDCCalibConfig; }; + + private: + std::array, NTDCChannels> mHTDC{}; + std::array, NTDCChannels> mMn{}; + bool mInitDone = false; + static std::mutex mMtx; /// mutex for critical section + const ZDCTDCParam* mTDCParam = nullptr; /// TDC calibration object + const TDCCalibConfig* mTDCCalibConfig = nullptr; /// Configuration of TDC calibration +}; +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/include/ZDCCalib/TDCCalibConfig.h b/Detectors/ZDC/calib/include/ZDCCalib/TDCCalibConfig.h new file mode 100644 index 0000000000000..42dc863b8e97e --- /dev/null +++ b/Detectors/ZDC/calib/include/ZDCCalib/TDCCalibConfig.h @@ -0,0 +1,65 @@ +// 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_ZDC_INTERCALIBCONFIG_H +#define O2_ZDC_TDCCALIBCONFIG_H + +#include "ZDCBase/Constants.h" +#include +#include + +/// \file TDCCalibConfig.h +/// \brief Configuration of ZDC TDC calibration procedure +/// \author P. Cortese + +namespace o2 +{ +namespace zdc +{ +struct TDCCalibConfig { + double cutLow[NTDCChannels] = {-std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity(), -std::numeric_limits::infinity()}; + double cutHigh[NTDCChannels] = {std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), + std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity()}; + int nb[NH] = {0}; /// Number of bins + double amin[NH] = {0}; /// minimum + double amax[NH] = {0}; /// maximum + double l_bnd[NH] = {0.1, 0.1, 0.1, 0.1, 0.1}; + double u_bnd[NH] = {10., 10., 10., 10., 10.}; + double l_bnd_o[NH] = {-20., -20., -20., -20., -20.}; + double u_bnd_o[NH] = {20., 20., 20., 20., 20.}; + double step_o[NH] = {0., 0., 0., 0., 0.}; + double min_e[NH] = {0., 0., 0., 0., 0.}; + + void print() const; + void resetCuts(); + void resetCutLow(); + void resetCutHigh(); + void resetCutLow(int ih); + void resetCutHigh(int ih); + void setMinEntries(double val); + void setMinEntries(int ih, double val); + void setCutLow(double val); + void setCutHigh(double val); + void setCutLow(int ih, double val); + void setCutHigh(int ih, double val); + void setCuts(double low, double high); + void setCuts(int ih, double low, double high); + void setBinning1D(int nb, double amin, double amax); + void setBinning2D(int nb, double amin, double amax); + void setBinning1D(int ih, int nb, double amin, double amax); + void setBinning2D(int ih, int nb, double amin, double amax); + ClassDefNV(TDCCalibConfig, 1); +}; +} // namespace zdc +} // namespace o2 + +#endif diff --git a/Detectors/ZDC/calib/src/InterCalib.cxx b/Detectors/ZDC/calib/src/InterCalib.cxx new file mode 100644 index 0000000000000..dd00d0fc65df7 --- /dev/null +++ b/Detectors/ZDC/calib/src/InterCalib.cxx @@ -0,0 +1,509 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include +#include "CommonUtils/MemFileHelper.h" +#include "ZDCCalib/InterCalib.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "Framework/Logger.h" +#include "CCDB/CcdbApi.h" + +using namespace o2::zdc; + +double InterCalib::mAdd[InterCalib::NPAR][InterCalib::NPAR] = {0}; +std::mutex InterCalib::mMtx; + +int InterCalib::init() +{ + if (mInterCalibConfig == nullptr) { + LOG(fatal) << "o2::zdc::InterCalib: missing configuration object"; + return -1; + } + clear(); + auto* cfg = mInterCalibConfig; + int ih; + // clang-format off + ih = 0; mHUnc[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 1; mHUnc[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 2; mHUnc[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 3; mHUnc[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 4; mHUnc[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 0; mHUnc[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 1; mHUnc[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 2; mHUnc[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 3; mHUnc[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 4; mHUnc[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 0; mCUnc[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 1; mCUnc[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 2; mCUnc[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 3; mCUnc[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 4; mCUnc[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + // clang-format on + mInitDone = true; + return 0; +} + +int InterCalib::process(const gsl::span& RecBC, + const gsl::span& Energy, + const gsl::span& TDCData, + const gsl::span& Info) +{ + if (!mInitDone) { + init(); + } + LOG(info) << "o2::zdc::InterCalib processing " << RecBC.size() << " b.c."; + o2::zdc::RecEventFlat ev; + ev.init(RecBC, Energy, TDCData, Info); + while (ev.next()) { + if (ev.getNInfo() > 0) { + auto& decodedInfo = ev.getDecodedInfo(); + for (uint16_t info : decodedInfo) { + uint8_t ch = (info >> 10) & 0x1f; + uint16_t code = info & 0x03ff; + // hmsg->Fill(ch, code); + } + if (mVerbosity > DbgMinimal) { + ev.print(); + } + // Need clean data (no messages) + // We are sure there is no pile-up in any channel (too restrictive?) + continue; + } + if (ev.getNEnergy() > 0 && ev.mCurB.triggers == 0) { + LOGF(info, "%9u.%04u Untriggered bunch", ev.mCurB.ir.orbit, ev.mCurB.ir.bc); + // Skip! + continue; + } + if ((ev.ezdcDecoded & MaskZNA) == MaskZNA) { + cumulate(HidZNA, ev.EZDC(IdZNAC), ev.EZDC(IdZNA1), ev.EZDC(IdZNA2), ev.EZDC(IdZNA3), ev.EZDC(IdZNA4), 1.); + } + if ((ev.ezdcDecoded & MaskZPA) == MaskZPA) { + cumulate(HidZPA, ev.EZDC(IdZPAC), ev.EZDC(IdZPA1), ev.EZDC(IdZPA2), ev.EZDC(IdZPA3), ev.EZDC(IdZPA4), 1.); + } + if ((ev.ezdcDecoded & MaskZNC) == MaskZNC) { + cumulate(HidZNC, ev.EZDC(IdZNCC), ev.EZDC(IdZNC1), ev.EZDC(IdZNC2), ev.EZDC(IdZNC3), ev.EZDC(IdZNC4), 1.); + } + if ((ev.ezdcDecoded & MaskZPC) == MaskZPC) { + cumulate(HidZPC, ev.EZDC(IdZPCC), ev.EZDC(IdZPC1), ev.EZDC(IdZPC2), ev.EZDC(IdZPC3), ev.EZDC(IdZPC4), 1.); + } + if ((ev.ezdcDecoded & MaskZEM) == MaskZEM) { + cumulate(HidZEM, ev.EZDC(IdZEM1), ev.EZDC(IdZEM2), 0., 0., 0., 1.); + } + } + return 0; +} + +//______________________________________________________________________________ +// Update calibration coefficients +int InterCalib::endOfRun() +{ + if (mVerbosity > DbgZero) { + LOGF(info, "Computing intercalibration coefficients"); + } + for (int ih = 0; ih < NH; ih++) { + LOGF(info, "%s %g events and cuts (%g:%g)", InterCalibData::DN[ih], mData.mSum[ih][5][5], mInterCalibConfig->cutLow[ih], mInterCalibConfig->cutHigh[ih]); + if (!mInterCalibConfig->enabled[ih]) { + LOGF(info, "DISABLED processing of RUN3 data for ih = %d: %s", ih, InterCalibData::DN[ih]); + assign(ih, false); + } else if (mData.mSum[ih][5][5] >= mInterCalibConfig->min_e[ih]) { + int ierr = mini(ih); + if (ierr) { + LOGF(error, "FAILED processing RUN3 data for ih = %d: %s", ih, InterCalibData::DN[ih]); + assign(ih, false); + } else { + LOGF(info, "Processed RUN3 data for ih = %d: %s", ih, InterCalibData::DN[ih]); + assign(ih, true); + } + } else { + LOGF(info, "FAILED processing RUN3 data for ih = %d: %s: TOO FEW EVENTS: %g", ih, InterCalibData::DN[ih], mData.mSum[ih][5][5]); + assign(ih, false); + } + } + + auto clName = o2::utils::MemFileHelper::getClassName(mTowerParamUpd); + mInfo.setObjectType(clName); + auto flName = o2::ccdb::CcdbApi::generateFileName(clName); + mInfo.setFileName(flName); + mInfo.setPath(CCDBPathTowerCalib); + std::map md; + md["config"] = mInterCalibConfig->desc; + mInfo.setMetaData(md); + uint64_t starting = mData.mCTimeBeg; + if (starting >= 10000) { + starting = starting - 10000; // start 10 seconds before + } + uint64_t stopping = mData.mCTimeEnd + 10000; // stop 10 seconds after + mInfo.setStartValidityTimestamp(starting); + mInfo.setEndValidityTimestamp(stopping); + + if (mSaveDebugHistos) { + write(); + } + return 0; +} + +//______________________________________________________________________________ +// Update calibration object for the five detectors +// ismod=false if it was not possible to update the calibration coefficients +// due to low statistics or minimization error +// ismod=true if the calibration was updated +void InterCalib::assign(int ih, bool ismod) +{ + int id_0[4] = {IdZNA1, IdZNA2, IdZNA3, IdZNA4}; + int id_1[4] = {IdZPA1, IdZPA2, IdZPA3, IdZPA4}; + int id_2[4] = {IdZNC1, IdZNC2, IdZNC3, IdZNC4}; + int id_3[4] = {IdZPC1, IdZPC2, IdZPC3, IdZPC4}; + int id_4[1] = {IdZEM2}; + int nid = 0; + int* id = nullptr; + if (ih == 0) { + nid = 4; + id = id_0; + } else if (ih == 1) { + nid = 4; + id = id_1; + } else if (ih == 2) { + nid = 4; + id = id_2; + } else if (ih == 3) { + nid = 4; + id = id_3; + } else if (ih == 4) { + nid = 1; + id = id_4; + } else { + LOG(fatal) << "InterCalib::assign accessing not existing ih = " << ih; + } + for (int iid = 0; iid < nid; iid++) { + auto ich = id[iid]; + auto oldval = mTowerParam->getTowerCalib(ich); + if (ismod == true) { + auto val = oldval; + if (oldval > 0) { + val = val * mPar[ih][iid + 1]; + } + if (mVerbosity > DbgZero) { + LOGF(info, "%s updated %8.6f -> %8.6f", ChannelNames[ich].data(), oldval, val); + } + mTowerParamUpd.setTowerCalib(ich, val, true); + } else { + if (mVerbosity > DbgZero) { + LOGF(info, "%s NOT CHANGED %8.6f", ChannelNames[ich].data(), oldval); + } + mTowerParamUpd.setTowerCalib(ich, oldval, false); + } + } +} + +int InterCalib::process(const char* hname, int ic) +{ + // Run 2 ZDC calibration is based on multi-dimensional histograms + // with dimensions: + // TC, T1, T2, T3, T4, Trigger + // ic is the number of the selected trigger class + THnSparse* hs = (THnSparse*)gROOT->FindObject(hname); + if (hs == 0) { + LOGF(error, "Not found: %s\n", hname); + return -1; + } + if (!hs->InheritsFrom(THnSparse::Class())) { + LOGF(error, "Not a THnSparse: %s\n", hname); + hs->IsA()->Print(); + return -1; + } + if (!mInitDone) { + init(); + } + TString hn = hname; + int ih = -1; + if (hn.EqualTo("hZNA")) { + ih = HidZNA; + } else if (hn.EqualTo("hZPA")) { + ih = HidZPA; + } else if (hn.EqualTo("hZNC")) { + ih = HidZNC; + } else if (hn.EqualTo("hZPC")) { + ih = HidZPC; + } else if (hn.EqualTo("hZEM")) { + ih = HidZEM; + } else { + LOGF(error, "Not recognized histogram name: %s\n", hname); + return -1; + } + clear(ih); + const int32_t dim = 6; + double x[dim]; + int32_t bins[dim]; + int64_t nb = hs->GetNbins(); + int64_t nn = 0; + LOGF(info, "Histogram %s has %ld bins\n", hname, nb); + double cutl = mInterCalibConfig->cutLow[ih]; + double cuth = mInterCalibConfig->cutHigh[ih]; + double contt = 0; + for (int64_t i = 0; i < nb; i++) { + double cont = hs->GetBinContent(i, bins); + if (cont <= 0) { + continue; + } + for (int32_t d = 0; d < dim; ++d) { + x[d] = hs->GetAxis(d)->GetBinCenter(bins[d]); + } + if (TMath::Nint(x[5] - ic) == 0 && x[0] > cutl && x[0] < cuth) { + nn++; + contt += cont; + cumulate(ih, x[0], x[1], x[2], x[3], x[4], cont); + } + } + int ierr = mini(ih); + if (ierr) { + LOGF(error, "FAILED processing RUN2 data for %s ih = %d\n", hname, ih); + } else { + LOGF(info, "Processed RUN2 data for %s ih = %d\n", hname, ih); + replay(ih, hs, ic); + } + LOGF(info, "Trigger class selection %d and %d bins %g events and cuts (%g:%g): %ld\n", ic, nn, contt, cutl, cuth); + return 0; +} + +void InterCalib::replay(int ih, THnSparse* hs, int ic) +{ + auto* cfg = mInterCalibConfig; + // clang-format off + if(ih == 0)mHCorr[ih] = std::make_unique("hZNASc","ZNA sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 1)mHCorr[ih] = std::make_unique("hZPASc","ZPA sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 2)mHCorr[ih] = std::make_unique("hZNCSc","ZNC sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 3)mHCorr[ih] = std::make_unique("hZPCSc","ZPC sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 4)mHCorr[ih] = std::make_unique("hZEM2c","ZEM2" ,cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 0)mCCorr[ih] = std::make_unique("cZNAc","ZNA;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 1)mCCorr[ih] = std::make_unique("cZPAc","ZPA;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 2)mCCorr[ih] = std::make_unique("cZNCc","ZNC;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 3)mCCorr[ih] = std::make_unique("cZPCc","ZPC;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 4)mCCorr[ih] = std::make_unique("cZEMc","ZEM;ZEM1;ZEM2 corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + // clang-format on + const int32_t dim = 6; + double x[dim]; + int32_t bins[dim]; + int64_t nb = hs->GetNbins(); + double cutl = cfg->cutLow[ih]; + double cuth = cfg->cutHigh[ih]; + double c1 = mPar[ih][1]; + double c2 = mPar[ih][2]; + double c3 = mPar[ih][3]; + double c4 = mPar[ih][4]; + double of = mPar[ih][5]; + for (int64_t i = 0; i < nb; i++) { + double cont = hs->GetBinContent(i, bins); + if (cont <= 0) { + continue; + } + for (int32_t d = 0; d < dim; ++d) { + x[d] = hs->GetAxis(d)->GetBinCenter(bins[d]); + } + if (TMath::Nint(x[5] - ic) == 0 && x[0] > cutl && x[0] < cuth) { + double sumquad = c1 * x[1] + c2 * x[2] + c3 * x[3] + c4 * x[4] + of; + mHCorr[ih]->Fill(sumquad, cont); + mCCorr[ih]->Fill(x[0], sumquad, cont); + } + } +} + +void InterCalib::clear(int ih) +{ + int ihstart = 0; + int ihstop = NH; + if (ih >= 0 && ih < NH) { + ihstart = ih; + ihstop = ih + 1; + } + for (int32_t ii = ihstart; ii < ihstop; ii++) { + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + mData.mSum[ii][i][j] = 0; + } + } + if (mHUnc[ii]) { + mHUnc[ii]->clear(); + } + if (mCUnc[ii]) { + mCUnc[ii]->clear(); + } + } +} + +int InterCalib::process(const InterCalibData& data) +{ + if (!mInitDone) { + init(); + } + mData += data; + return 0; +} + +void InterCalib::add(int ih, o2::dataformats::FlatHisto1D& h1) +{ + if (!mInitDone) { + init(); + } + constexpr int nh = 2 * InterCalibData::NH; + if (ih >= 0 && ih < nh) { + mHUnc[ih]->add(h1); + } else { + LOG(error) << "InterCalib::add: unsupported FlatHisto1D " << ih; + } +} + +void InterCalib::add(int ih, o2::dataformats::FlatHisto2D& h2) +{ + if (!mInitDone) { + init(); + } + if (ih >= 0 && ih < InterCalibData::NH) { + mCUnc[ih]->add(h2); + } else { + LOG(error) << "InterCalib::add: unsupported FlatHisto2D " << ih; + } +} + +void InterCalib::cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w = 1) +{ + if (tc < mInterCalibConfig->cutLow[ih] || tc > mInterCalibConfig->cutHigh[ih]) { + return; + } + double val[NPAR] = {0, 0, 0, 0, 0, 1}; + val[0] = tc; + val[1] = t1; + val[2] = t2; + val[3] = t3; + val[4] = t4; + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = i; j < NPAR; j++) { + mData.mSum[ih][i][j] += val[i] * val[j] * w; + } + } + // mData.mSum[ih][5][5] contains the number of analyzed events + double sumquad = val[1] + val[2] + val[3] + val[4]; + mHUnc[ih]->fill(sumquad, w); + mHUnc[ih + NH]->fill(val[0]); + mCUnc[ih]->fill(val[0], sumquad, w); +} + +//______________________________________________________________________________ +// Compute chi2 for minimization (static function) +void InterCalib::fcn(int& npar, double* gin, double& chi, double* par, int iflag) +{ + chi = 0; + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + chi += (i == 0 ? par[i] : -par[i]) * (j == 0 ? par[j] : -par[j]) * mAdd[i][j]; + } + } +} + +int InterCalib::mini(int ih) +{ + // Copy to static object and symmetrize matrix + // We use a static function and therefore we can do only one minimization at a time + mMtx.lock(); + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + if (j < i) { + mAdd[i][j] = mData.mSum[ih][j][i]; + } else { + mAdd[i][j] = mData.mSum[ih][i][j]; + } + } + } + double arglist[10]; + int ierflg = 0; + double l_bnd = mInterCalibConfig->l_bnd[ih]; + double u_bnd = mInterCalibConfig->u_bnd[ih]; + double start = 1.0; + double step = 0.1; + mMn[ih] = std::make_unique(NPAR); + mMn[ih]->SetFCN(fcn); + mMn[ih]->mnparm(0, "c0", 1., 0., 1., 1., ierflg); + mMn[ih]->mnparm(1, "c1", start, step, l_bnd, u_bnd, ierflg); + if (ih == 4) { + // Only two ZEM calorimeters: equalize response + l_bnd = 0; + u_bnd = 0; + start = 0; + step = 0; + } + mMn[ih]->mnparm(2, "c2", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnparm(3, "c3", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnparm(4, "c4", start, step, l_bnd, u_bnd, ierflg); + l_bnd = mInterCalibConfig->l_bnd_o[ih]; + u_bnd = mInterCalibConfig->u_bnd_o[ih]; + start = 0; + step = mInterCalibConfig->step_o[ih]; + mMn[ih]->mnparm(5, "offset", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnexcm("MIGRAD", arglist, 0, ierflg); + for (Int_t i = 0; i < NPAR; i++) { + mMn[ih]->GetParameter(i, mPar[ih][i], mErr[ih][i]); + } + mMtx.unlock(); + return ierflg; +} + +int InterCalib::write(const std::string fn) +{ + TDirectory* cwd = gDirectory; + TFile* f = new TFile(fn.data(), "recreate"); + if (f->IsZombie()) { + LOG(error) << "Cannot create file: " << fn; + return 1; + } + for (int32_t ih = 0; ih < (2 * NH); ih++) { + if (mHUnc[ih]) { + auto p = mHUnc[ih]->createTH1F(InterCalib::mHUncN[ih]); + p->SetTitle(InterCalib::mHUncT[ih]); + p->Write("", TObject::kOverwrite); + } + } + for (int32_t ih = 0; ih < NH; ih++) { + if (mCUnc[ih]) { + auto p = mCUnc[ih]->createTH2F(InterCalib::mCUncN[ih]); + p->SetTitle(InterCalib::mCUncT[ih]); + p->Write("", TObject::kOverwrite); + } + } + // Only after replay of RUN2 data + for (int32_t ih = 0; ih < NH; ih++) { + if (mHCorr[ih]) { + mHCorr[ih]->Write("", TObject::kOverwrite); + } + } + for (int32_t ih = 0; ih < NH; ih++) { + if (mCCorr[ih]) { + mCCorr[ih]->Write("", TObject::kOverwrite); + } + } + // Minimization output + const char* mntit[NH] = {"mZNA", "mZPA", "mZNC", "mZPC", "mZEM"}; + for (int32_t ih = 0; ih < NH; ih++) { + if (mMn[ih]) { + mMn[ih]->Write(mntit[ih], TObject::kOverwrite); + } + } + f->Close(); + cwd->cd(); + return 0; +} diff --git a/Detectors/ZDC/calib/src/InterCalibConfig.cxx b/Detectors/ZDC/calib/src/InterCalibConfig.cxx new file mode 100644 index 0000000000000..b112764fd188a --- /dev/null +++ b/Detectors/ZDC/calib/src/InterCalibConfig.cxx @@ -0,0 +1,143 @@ +// 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. + +#include "Framework/Logger.h" +#include "ZDCCalib/InterCalibConfig.h" + +using namespace o2::zdc; + +void InterCalibConfig::print() const +{ + const char* hn[NH] = {"ZNA", "ZPA", "ZNC", "ZPC", "ZEM"}; + for (Int_t ih = 0; ih < NH; ih++) { + LOG(info) << hn[ih] << " limits = (" << cutLow[ih] << " : " << cutHigh[ih] << ")"; + } + for (Int_t ih = 0; ih < NH; ih++) { + LOG(info) << hn[ih] << " booking 1D = (" << nb1[ih] << ", " << amin1[ih] << ", " << amax1[ih] << ")"; + } + for (Int_t ih = 0; ih < NH; ih++) { + LOG(info) << hn[ih] << " booking 2D = (" << nb2[ih] << ", " << amin2[ih] << ", " << amax2[ih] << ")"; + } +} + +void InterCalibConfig::setMinEntries(double val) +{ + for (int32_t ih = 0; ih < NH; ih++) { + min_e[ih] = val; + } +} + +void InterCalibConfig::setMinEntries(int ih, double val) +{ + min_e[ih] = val; +} + +void InterCalibConfig::resetCuts() +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutLow[ih] = -std::numeric_limits::infinity(); + cutHigh[ih] = std::numeric_limits::infinity(); + } +} + +void InterCalibConfig::resetCutLow() +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutLow[ih] = -std::numeric_limits::infinity(); + } +} + +void InterCalibConfig::resetCutHigh() +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutHigh[ih] = std::numeric_limits::infinity(); + } +} + +void InterCalibConfig::resetCutLow(int ih) +{ + cutLow[ih] = -std::numeric_limits::infinity(); +} + +void InterCalibConfig::resetCutHigh(int ih) +{ + cutHigh[ih] = std::numeric_limits::infinity(); +} + +void InterCalibConfig::setCutLow(double val) +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutLow[ih] = val; + } +} + +void InterCalibConfig::setCutHigh(double val) +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutHigh[ih] = val; + } +} + +void InterCalibConfig::setCutLow(int ih, double val) +{ + cutLow[ih] = val; +} + +void InterCalibConfig::setCutHigh(int ih, double val) +{ + cutHigh[ih] = val; +} + +void InterCalibConfig::setCuts(double low, double high) +{ + for (int32_t ih = 0; ih < NH; ih++) { + cutLow[ih] = low; + cutHigh[ih] = high; + } +} + +void InterCalibConfig::setCuts(int ih, double low, double high) +{ + cutHigh[ih] = low; + cutLow[ih] = high; +} + +void InterCalibConfig::setBinning1D(int nb, double amin, double amax) +{ + for (int32_t ih = 0; ih < NH; ih++) { + nb1[ih] = nb; + amin1[ih] = amin; + amax1[ih] = amax; + } +} + +void InterCalibConfig::setBinning2D(int nb, double amin, double amax) +{ + for (int32_t ih = 0; ih < NH; ih++) { + nb2[ih] = nb; + amin2[ih] = amin; + amax2[ih] = amax; + } +} + +void InterCalibConfig::setBinning1D(int ih, int nb, double amin, double amax) +{ + nb1[ih] = nb; + amin1[ih] = amin; + amax1[ih] = amax; +} + +void InterCalibConfig::setBinning2D(int ih, int nb, double amin, double amax) +{ + nb2[ih] = nb; + amin2[ih] = amin; + amax2[ih] = amax; +} diff --git a/Detectors/ZDC/calib/src/InterCalibEPN.cxx b/Detectors/ZDC/calib/src/InterCalibEPN.cxx new file mode 100644 index 0000000000000..d3f36374c5e05 --- /dev/null +++ b/Detectors/ZDC/calib/src/InterCalibEPN.cxx @@ -0,0 +1,254 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include +#include "ZDCCalib/InterCalibEPN.h" +#include "ZDCCalib/InterCalib.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "DataFormatsZDC/InterCalibData.h" +#include "Framework/Logger.h" + +using namespace o2::zdc; + +int InterCalibEPN::init() +{ + if (mInterCalibConfig == nullptr) { + LOG(fatal) << "o2::zdc::InterCalibEPN: missing configuration object"; + return -1; + } + clear(); + auto* cfg = mInterCalibConfig; + int ih; + // clang-format off + ih = 0; mH[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 1; mH[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 2; mH[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 3; mH[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 4; mH[ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 0; mH[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 1; mH[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 2; mH[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 3; mH[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 4; mH[NH+ih] = new o2::dataformats::FlatHisto1D(cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + ih = 0; mC[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 1; mC[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 2; mC[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 3; mC[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + ih = 4; mC[ih] = new o2::dataformats::FlatHisto2D(cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + // clang-format on + mInitDone = true; + return 0; +} + +int InterCalibEPN::process(const gsl::span& RecBC, + const gsl::span& Energy, + const gsl::span& TDCData, + const gsl::span& Info) +{ + if (!mInitDone) { + init(); + } + LOG(info) << "o2::zdc::InterCalibEPN processing " << RecBC.size() << " b.c. @ TS " << mData.mCTimeBeg << " : " << mData.mCTimeEnd; + o2::zdc::RecEventFlat ev; + ev.init(RecBC, Energy, TDCData, Info); + while (ev.next()) { + if (ev.getNInfo() > 0) { + auto& decodedInfo = ev.getDecodedInfo(); + for (uint16_t info : decodedInfo) { + uint8_t ch = (info >> 10) & 0x1f; + uint16_t code = info & 0x03ff; + // hmsg->Fill(ch, code); + } + if (mVerbosity > DbgMinimal) { + ev.print(); + } + // Need clean data (no messages) + // We are sure there is no pile-up in any channel (too restrictive?) + continue; + } + if (ev.getNEnergy() > 0 && ev.mCurB.triggers == 0) { + LOGF(info, "%9u.%04u Untriggered bunch", ev.mCurB.ir.orbit, ev.mCurB.ir.bc); + // Skip! + continue; + } + if ((ev.ezdcDecoded & MaskZNA) == MaskZNA) { + cumulate(HidZNA, ev.EZDC(IdZNAC), ev.EZDC(IdZNA1), ev.EZDC(IdZNA2), ev.EZDC(IdZNA3), ev.EZDC(IdZNA4), 1.); + } + if ((ev.ezdcDecoded & MaskZPA) == MaskZPA) { + cumulate(HidZPA, ev.EZDC(IdZPAC), ev.EZDC(IdZPA1), ev.EZDC(IdZPA2), ev.EZDC(IdZPA3), ev.EZDC(IdZPA4), 1.); + } + if ((ev.ezdcDecoded & MaskZNC) == MaskZNC) { + cumulate(HidZNC, ev.EZDC(IdZNCC), ev.EZDC(IdZNC1), ev.EZDC(IdZNC2), ev.EZDC(IdZNC3), ev.EZDC(IdZNC4), 1.); + } + if ((ev.ezdcDecoded & MaskZPC) == MaskZPC) { + cumulate(HidZPC, ev.EZDC(IdZPCC), ev.EZDC(IdZPC1), ev.EZDC(IdZPC2), ev.EZDC(IdZPC3), ev.EZDC(IdZPC4), 1.); + } + if ((ev.ezdcDecoded & MaskZEM) == MaskZEM) { + cumulate(HidZEM, ev.EZDC(IdZEM1), ev.EZDC(IdZEM2), 0., 0., 0., 1.); + } + } + return 0; +} + +int InterCalibEPN::endOfRun() +{ + if (mVerbosity > DbgZero) { + LOGF(info, "InterCalibEPN::endOfRun ts (%llu:%llu)", mData.mCTimeBeg, mData.mCTimeEnd); + for (int ih = 0; ih < NH; ih++) { + LOGF(info, "%s %g events and cuts (%g:%g)", InterCalibData::DN[ih], mData.mSum[ih][5][5], mInterCalibConfig->cutLow[ih], mInterCalibConfig->cutHigh[ih]); + } + } + if (mSaveDebugHistos) { + write(); + } + return 0; +} + +int InterCalibEPN::process(const char* hname, int ic) +{ + // Run 2 ZDC calibration is based on multi-dimensional histograms + // with dimensions: + // TC, T1, T2, T3, T4, Trigger + // ic is the number of the selected trigger class + THnSparse* hs = (THnSparse*)gROOT->FindObject(hname); + if (hs == 0) { + LOGF(error, "Not found: %s\n", hname); + return -1; + } + if (!hs->InheritsFrom(THnSparse::Class())) { + LOGF(error, "Not a THnSparse: %s\n", hname); + hs->IsA()->Print(); + return -1; + } + TString hn = hname; + int ih = -1; + if (hn.EqualTo("hZNA")) { + ih = HidZNA; + } else if (hn.EqualTo("hZPA")) { + ih = HidZPA; + } else if (hn.EqualTo("hZNC")) { + ih = HidZNC; + } else if (hn.EqualTo("hZPC")) { + ih = HidZPC; + } else if (hn.EqualTo("hZEM")) { + ih = HidZEM; + } else { + LOGF(error, "Not recognized histogram name: %s\n", hname); + return -1; + } + clear(ih); + const int32_t dim = 6; + double x[dim]; + int32_t bins[dim]; + int64_t nb = hs->GetNbins(); + int64_t nn = 0; + LOGF(info, "Histogram %s has %ld bins\n", hname, nb); + double cutl = mInterCalibConfig->cutLow[ih]; + double cuth = mInterCalibConfig->cutHigh[ih]; + double contt = 0; + for (int64_t i = 0; i < nb; i++) { + double cont = hs->GetBinContent(i, bins); + if (cont <= 0) { + continue; + } + for (int32_t d = 0; d < dim; ++d) { + x[d] = hs->GetAxis(d)->GetBinCenter(bins[d]); + } + if (TMath::Nint(x[5] - ic) == 0 && x[0] > cutl && x[0] < cuth) { + nn++; + contt += cont; + cumulate(ih, x[0], x[1], x[2], x[3], x[4], cont); + } + } + LOGF(info, "Trigger class selection %d and %d bins %g events and cuts (%g:%g): %ld", ic, nn, contt, cutl, cuth); + return 0; +} + +void InterCalibEPN::clear(int ih) +{ + int ihstart = 0; + int ihstop = NH; + if (ih >= 0 && ih < NH) { + ihstart = ih; + ihstop = ih + 1; + } + for (int32_t ii = ihstart; ii < ihstop; ii++) { + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + mData.mSum[ii][i][j] = 0; + } + } + if (mH[ii]) { + mH[ii]->clear(); + } + if (mC[ii]) { + mC[ii]->clear(); + } + } +} + +void InterCalibEPN::cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w = 1) +{ + if (tc < mInterCalibConfig->cutLow[ih] || tc > mInterCalibConfig->cutHigh[ih]) { + return; + } + double val[NPAR] = {0, 0, 0, 0, 0, 1}; + val[0] = tc; + val[1] = t1; + val[2] = t2; + val[3] = t3; + val[4] = t4; + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = i; j < NPAR; j++) { + mData.mSum[ih][i][j] += val[i] * val[j] * w; + } + } + // mData.mSum[ih][5][5] contains the number of analyzed events + double sumquad = val[1] + val[2] + val[3] + val[4]; + mH[ih]->fill(sumquad, w); + mH[ih + NH]->fill(val[0]); + mC[ih]->fill(val[0], sumquad, w); +} + +int InterCalibEPN::write(const std::string fn) +{ + TDirectory* cwd = gDirectory; + TFile* f = new TFile(fn.data(), "recreate"); + if (f->IsZombie()) { + LOG(error) << "Cannot create file: " << fn; + return 1; + } + for (int32_t ih = 0; ih < (2 * NH); ih++) { + if (mH[ih]) { + auto p = mH[ih]->createTH1F(InterCalib::mHUncN[ih]); + p->SetTitle(InterCalib::mHUncT[ih]); + p->Write("", TObject::kOverwrite); + } + } + for (int32_t ih = 0; ih < NH; ih++) { + if (mC[ih]) { + auto p = mC[ih]->createTH2F(InterCalib::mCUncN[ih]); + p->SetTitle(InterCalib::mCUncT[ih]); + p->Write("", TObject::kOverwrite); + } + } + f->Close(); + cwd->cd(); + return 0; +} diff --git a/Detectors/ZDC/calib/src/InterCalibEPNSpec.cxx b/Detectors/ZDC/calib/src/InterCalibEPNSpec.cxx new file mode 100644 index 0000000000000..99d305d2c7031 --- /dev/null +++ b/Detectors/ZDC/calib/src/InterCalibEPNSpec.cxx @@ -0,0 +1,160 @@ +// 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 InterCalibEPNSpec.cxx +/// @brief ZDC reconstruction +/// @author pietro.cortese@cern.ch + +#include +#include +#include +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "Framework/Logger.h" +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/DataRefUtils.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DataFormatsZDC/BCData.h" +#include "DataFormatsZDC/ChannelData.h" +#include "DataFormatsZDC/OrbitData.h" +#include "DataFormatsZDC/RecEvent.h" +#include "ZDCBase/ModuleConfig.h" +#include "CommonUtils/NameConf.h" +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "ZDCReconstruction/RecoConfigZDC.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "ZDCCalib/InterCalibEPNSpec.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace zdc +{ + +InterCalibEPNSpec::InterCalibEPNSpec() +{ + mTimer.Stop(); + mTimer.Reset(); +} + +InterCalibEPNSpec::InterCalibEPNSpec(const int verbosity) : mVerbosity(verbosity) +{ + mTimer.Stop(); + mTimer.Reset(); +} + +void InterCalibEPNSpec::init(o2::framework::InitContext& ic) +{ + mVerbosity = ic.options().get("verbosity-level"); +} + +void InterCalibEPNSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + // we call these methods just to trigger finaliseCCDB callback + pc.inputs().get("intercalibconfig"); +} + +void InterCalibEPNSpec::run(ProcessingContext& pc) +{ + updateTimeDependentParams(pc); + if (!mInitialized) { + mInitialized = true; + std::string loadedConfFiles = "Loaded ZDC configuration files:"; + // InterCalib configuration + auto interConfig = pc.inputs().get("intercalibconfig"); + if (!interConfig) { + LOG(fatal) << "Missing InterCalibConfig calibration InterCalibConfig"; + return; + } else { + loadedConfFiles += " InterCalibConfig"; + if (mVerbosity > DbgZero) { + LOG(info) << "Loaded InterCalib configuration object"; + interConfig->print(); + } + } + + mInterCalibEPN.setInterCalibConfig(interConfig.get()); + + LOG(info) << loadedConfFiles; + mTimer.CpuTime(); + mTimer.Start(false); + } + + const auto ref = pc.inputs().getFirstValid(true); + auto creationTime = DataRefUtils::getHeader(ref)->creation; // approximate time in ms + mInterCalibEPN.getData().setCreationTime(creationTime); + + auto bcrec = pc.inputs().get>("bcrec"); + auto energy = pc.inputs().get>("energy"); + auto tdc = pc.inputs().get>("tdc"); + auto info = pc.inputs().get>("info"); + + // Process reconstructed data + mInterCalibEPN.process(bcrec, energy, tdc, info); + + // Send intermediate calibration data and debug histograms + o2::framework::Output output("ZDC", "INTERCALIBDATA", 0, Lifetime::Timeframe); + pc.outputs().snapshot(output, mInterCalibEPN.mData); + char outputd[o2::header::gSizeDataDescriptionString]; + for (int ih = 0; ih < (2 * InterCalibData::NH); ih++) { + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_1DH%d", ih); + o2::framework::Output output("ZDC", outputd, 0, Lifetime::Timeframe); + pc.outputs().snapshot(output, mInterCalibEPN.mH[ih]->getBase()); + } + for (int ih = 0; ih < InterCalibData::NH; ih++) { + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_2DH%d", ih); + o2::framework::Output output("ZDC", outputd, 0, Lifetime::Timeframe); + pc.outputs().snapshot(output, mInterCalibEPN.mC[ih]->getBase()); + } +} + +void InterCalibEPNSpec::endOfStream(EndOfStreamContext& ec) +{ + mInterCalibEPN.endOfRun(); + mTimer.Stop(); + LOGF(info, "ZDC EPN Intercalibration total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +framework::DataProcessorSpec getInterCalibEPNSpec() +{ + std::vector inputs; + inputs.emplace_back("bcrec", "ZDC", "BCREC", 0, Lifetime::Timeframe); + inputs.emplace_back("energy", "ZDC", "ENERGY", 0, Lifetime::Timeframe); + inputs.emplace_back("tdc", "ZDC", "TDCDATA", 0, Lifetime::Timeframe); + inputs.emplace_back("info", "ZDC", "INFO", 0, Lifetime::Timeframe); + inputs.emplace_back("intercalibconfig", "ZDC", "INTERCALIBCONFIG", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathInterCalibConfig.data()))); + + std::vector outputs; + outputs.emplace_back("ZDC", "INTERCALIBDATA", 0, Lifetime::Timeframe); + char outputd[o2::header::gSizeDataDescriptionString]; + for (int ih = 0; ih < (2 * InterCalibData::NH); ih++) { + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_1DH%d", ih); + outputs.emplace_back("ZDC", outputd, 0, Lifetime::Timeframe); + } + for (int ih = 0; ih < InterCalibData::NH; ih++) { + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_2DH%d", ih); + outputs.emplace_back("ZDC", outputd, 0, Lifetime::Timeframe); + } + return DataProcessorSpec{ + "zdc-intercalib-epn", + inputs, + outputs, + AlgorithmSpec{adaptFromTask()}, + o2::framework::Options{{"verbosity-level", o2::framework::VariantType::Int, 0, {"Verbosity level"}}}}; +} + +} // namespace zdc +} // namespace o2 diff --git a/Detectors/ZDC/calib/src/InterCalibSpec.cxx b/Detectors/ZDC/calib/src/InterCalibSpec.cxx new file mode 100644 index 0000000000000..3e0f02303f70f --- /dev/null +++ b/Detectors/ZDC/calib/src/InterCalibSpec.cxx @@ -0,0 +1,206 @@ +// 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 InterCalibSpec.cxx +/// @brief ZDC reconstruction +/// @author pietro.cortese@cern.ch + +#include +#include +#include +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "CCDB/CcdbApi.h" +#include "Framework/Logger.h" +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/CCDBParamSpec.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DataFormatsZDC/BCData.h" +#include "DataFormatsZDC/ChannelData.h" +#include "DataFormatsZDC/OrbitData.h" +#include "DataFormatsZDC/RecEvent.h" +#include "ZDCBase/ModuleConfig.h" +#include "CommonUtils/NameConf.h" +#include "CommonUtils/MemFileHelper.h" +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "ZDCReconstruction/RecoConfigZDC.h" +#include "ZDCReconstruction/ZDCEnergyParam.h" +#include "ZDCReconstruction/ZDCTowerParam.h" +#include "ZDCCalib/InterCalibSpec.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace zdc +{ + +InterCalibSpec::InterCalibSpec() +{ + mTimer.Stop(); + mTimer.Reset(); +} + +InterCalibSpec::InterCalibSpec(const int verbosity) : mVerbosity(verbosity) +{ + mTimer.Stop(); + mTimer.Reset(); +} + +void InterCalibSpec::init(o2::framework::InitContext& ic) +{ + // int minEnt = std::max(300, ic.options().get("min-entries")); + // int nb = std::max(500, ic.options().get("nbins")); + // int slotL = ic.options().get("tf-per-slot"); + // int delay = ic.options().get("max-delay"); + mVerbosity = ic.options().get("verbosity-level"); + mInterCalib.setVerbosity(mVerbosity); + mTimer.CpuTime(); + mTimer.Start(false); +} + +void InterCalibSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + // we call these methods just to trigger finaliseCCDB callback + std::string loadedConfFiles = "Loaded ZDC configuration files:"; + // Energy calibration + auto energyParam = pc.inputs().get("energycalib"); + if (!energyParam) { + LOG(fatal) << "Missing ZDCEnergyParam calibration object"; + return; + } else { + loadedConfFiles += " ZDCEnergyParam"; + if (mVerbosity > DbgMinimal) { + LOG(info) << "Loaded Energy calibration ZDCEnergyParam"; + energyParam->print(); + } + } + + // Tower calibration + auto towerParam = pc.inputs().get("towercalib"); + if (!towerParam) { + LOG(fatal) << "Missing ZDCTowerParam calibration object"; + return; + } else { + loadedConfFiles += " ZDCTowerParam"; + if (mVerbosity > DbgMinimal) { + LOG(info) << "Loaded Tower calibration ZDCTowerParam"; + towerParam->print(); + } + } + + // InterCalib configuration + auto interConfig = pc.inputs().get("intercalibconfig"); + if (!interConfig) { + LOG(fatal) << "Missing InterCalibConfig calibration InterCalibConfig"; + return; + } else { + loadedConfFiles += " InterCalibConfig"; + if (mVerbosity > DbgMinimal) { + LOG(info) << "Loaded InterCalib configuration object"; + interConfig->print(); + } + } + + LOG(info) << loadedConfFiles; + + mInterCalib.setEnergyParam(energyParam.get()); + mInterCalib.setTowerParam(towerParam.get()); + mInterCalib.setInterCalibConfig(interConfig.get()); +} + +void InterCalibSpec::run(ProcessingContext& pc) +{ + updateTimeDependentParams(pc); + auto data = pc.inputs().get("intercalibdata"); + mInterCalib.process(data); + for (int ih = 0; ih < (2 * InterCalibData::NH); ih++) { + o2::dataformats::FlatHisto1D histoView(pc.inputs().get>(fmt::format("inter_1dh{}", ih).data())); + mInterCalib.add(ih, histoView); + } + for (int ih = 0; ih < InterCalibData::NH; ih++) { + o2::dataformats::FlatHisto2D histoView(pc.inputs().get>(fmt::format("inter_2dh{}", ih).data())); + mInterCalib.add(ih, histoView); + } +} + +void InterCalibSpec::endOfStream(EndOfStreamContext& ec) +{ + mInterCalib.endOfRun(); + mTimer.Stop(); + sendOutput(ec.outputs()); + LOGF(info, "ZDC Intercalibration total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +//________________________________________________________________ +void InterCalibSpec::sendOutput(o2::framework::DataAllocator& output) +{ + // extract CCDB infos and calibration objects, convert it to TMemFile and send them to the output + // TODO in principle, this routine is generic, can be moved to Utils.h + using clbUtils = o2::calibration::Utils; + const auto& payload = mInterCalib.getTowerParamUpd(); + auto& info = mInterCalib.getCcdbObjectInfo(); + auto image = o2::ccdb::CcdbApi::createObjectImage(&payload, &info); + if (mVerbosity > DbgMinimal) { + payload.print(); + } + LOG(info) << "Sending object " << info.getPath() << "/" << info.getFileName() << " of size " << image->size() + << " bytes, valid for " << info.getStartValidityTimestamp() << " : " << info.getEndValidityTimestamp(); + output.snapshot(Output{o2::calibration::Utils::gDataOriginCDBPayload, "ZDC_Intercalib", 0}, *image.get()); // vector + output.snapshot(Output{o2::calibration::Utils::gDataOriginCDBWrapper, "ZDC_Intercalib", 0}, info); // root-serialized + // TODO: reset the outputs once they are already sent (is it necessary?) + // mInterCalib.init(); +} + +framework::DataProcessorSpec getInterCalibSpec() +{ + using device = o2::zdc::InterCalibSpec; + using clbUtils = o2::calibration::Utils; + + std::vector inputs; + inputs.emplace_back("intercalibconfig", "ZDC", "INTERCALIBCONFIG", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathInterCalibConfig.data()))); + inputs.emplace_back("energycalib", "ZDC", "ENERGYCALIB", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathEnergyCalib.data()))); + inputs.emplace_back("towercalib", "ZDC", "TOWERCALIB", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathTowerCalib.data()))); + inputs.emplace_back("intercalibdata", "ZDC", "INTERCALIBDATA", 0, Lifetime::Timeframe); + + char outputa[o2::header::gSizeDataDescriptionString]; + char outputd[o2::header::gSizeDataDescriptionString]; + for (int ih = 0; ih < (2 * InterCalibData::NH); ih++) { + snprintf(outputa, o2::header::gSizeDataDescriptionString, "inter_1dh%d", ih); + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_1DH%d", ih); + inputs.emplace_back(outputa, "ZDC", outputd, 0, Lifetime::Timeframe); + } + for (int ih = 0; ih < InterCalibData::NH; ih++) { + snprintf(outputa, o2::header::gSizeDataDescriptionString, "inter_2dh%d", ih); + snprintf(outputd, o2::header::gSizeDataDescriptionString, "INTER_2DH%d", ih); + inputs.emplace_back(outputa, "ZDC", outputd, 0, Lifetime::Timeframe); + } + + std::vector outputs; + outputs.emplace_back(ConcreteDataTypeMatcher{o2::calibration::Utils::gDataOriginCDBPayload, "ZDC_Intercalib"}, Lifetime::Sporadic); + outputs.emplace_back(ConcreteDataTypeMatcher{o2::calibration::Utils::gDataOriginCDBWrapper, "ZDC_Intercalib"}, Lifetime::Sporadic); + + return DataProcessorSpec{ + "zdc-calib-towers", + inputs, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{ + {"tf-per-slot", VariantType::Int, 5, {"number of TFs per calibration time slot"}}, + {"max-delay", VariantType::Int, 3, {"number of slots in past to consider"}}, + {"min-entries", VariantType::Int, 500, {"minimum number of entries to fit single time slot"}}, + {"verbosity-level", o2::framework::VariantType::Int, 1, {"Verbosity level"}}}}; +} + +} // namespace zdc +} // namespace o2 diff --git a/Detectors/ZDC/calib/src/TDCCalib.cxx b/Detectors/ZDC/calib/src/TDCCalib.cxx new file mode 100644 index 0000000000000..337f6a92c3537 --- /dev/null +++ b/Detectors/ZDC/calib/src/TDCCalib.cxx @@ -0,0 +1,353 @@ +// 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. + +#include +#include +#include +#include +#include +#include +#include +#include +#include "ZDCCalib/TDCCalib.h" +#include "ZDCReconstruction/ZDCTDCParam.h" +#include "Framework/Logger.h" + +using namespace o2::zdc; + +std::mutex TDCCalib::mMtx; + +int TDCCalib::init() +{ + if (mTDCCalibConfig == nullptr) { + LOG(fatal) << "o2::zdc::TDCCalib: missing configuration object"; + return -1; + } + clear(); + for (int ih = 0; ih < NTDCChannels; ih++) { + std::string hn = fmt::format("h{}", ChannelNames[TDCSignal[ih]]); + std::string ht = fmt::format("TDC {};(ch)", ChannelNames[TDCSignal[ih]]); + mHTDC[ih] = std::make_unique(hn.data(), ht.data(), mTDCCalibConfig->nb[ih], mTDCCalibConfig->amin[ih], mTDCCalibConfig->amax[ih]); + } + mInitDone = true; + return 0; +} + +int TDCCalib::process(const gsl::span& RecBC, + const gsl::span& Energy, + const gsl::span& TDCData, + const gsl::span& Info) +{ + if (!mInitDone) { + init(); + } + LOG(info) << "o2::zdc::TDCCalib processing " << RecBC.size() << " b.c."; + o2::zdc::RecEventFlat ev; + ev.init(RecBC, Energy, TDCData, Info); + while (ev.next()) { + if (ev.getNInfo() > 0) { + auto& decodedInfo = ev.getDecodedInfo(); + // for (uint16_t info : decodedInfo) { + // uint8_t ch = (info >> 10) & 0x1f; + // uint16_t code = info & 0x03ff; + // } + ev.print(); + // Need clean data (no messages) + // We are sure there is no pile-up in any channel (too restrictive?) + continue; + } + // TDC + for (int32_t itdc = 0; itdc < o2::zdc::NTDCChannels; itdc++) { + int nhit = ev.NtdcV(itdc); + if (nhit == 1) { + // TDC without in-bunch pile-up + int ich = o2::zdc::TDCSignal[itdc]; + if (ev.genericE[ich] || ev.tdcPileEvC[ich] || ev.tdcPileEvE[ich] || + ev.tdcPileM1C[ich] || ev.tdcPileM1E[ich] || ev.tdcPileM2C[ich] || ev.tdcPileM2E[ich] || ev.tdcPileM3C[ich] || ev.tdcPileM3E[ich] || + ev.tdcSigE[ich]) { + continue; + } + if (ev.NtdcA(itdc) != nhit) { + LOGF(error, "Mismatch in TDC %d data Val=%d Amp=%d", itdc, ev.NtdcV(itdc), ev.NtdcA(ich)); + continue; + } + auto amp = ev.TDCAmp[itdc][0]; + if (amp < mTDCCalibConfig.cutLow[itdc] || mTDCCalibConfig > cutHigh[itdc]) { + continue; + } + mHTDC[itdc]->Fill(ev.TDCVal[itdc][0]); + } + } + } + return 0; +} + +int TDCCalib::endOfRun() +{ + for (int ih = 0; ih < NH; ih++) { + if (mSum[ih][5][5] >= mTDCCalibConfig->min_e[ih]) { + int ierr = mini(ih); + if (ierr) { + LOGF(error, "FAILED processing RUN3 data for ih = %d - ", ih); + } else { + LOGF(info, "Processed RUN3 data for ih = %d: ", ih); + } + } else { + LOGF(info, "FAILED processing RUN3 data for ih = %d: TOO FEW EVENTS: ", ih); + } + LOGF(info, "%g events and cuts (%g:%g)\n", mSum[ih][5][5], mTDCCalibConfig->cutLow[ih], mTDCCalibConfig->cutHigh[ih]); + } + write(); + return 0; +} + +int TDCCalib::process(const char* hname, int ic) +{ + // Run 2 ZDC calibration is based on multi-dimensional histograms + // with dimensions: + // TC, T1, T2, T3, T4, Trigger + // ic is the number of the selected trigger class + THnSparse* hs = (THnSparse*)gROOT->FindObject(hname); + if (hs == 0) { + LOGF(error, "Not found: %s\n", hname); + return -1; + } + if (!hs->InheritsFrom(THnSparse::Class())) { + LOGF(error, "Not a THnSparse: %s\n", hname); + hs->IsA()->Print(); + return -1; + } + TString hn = hname; + int ih = -1; + if (hn.EqualTo("hZNA")) { + ih = HidZNA; + } else if (hn.EqualTo("hZPA")) { + ih = HidZPA; + } else if (hn.EqualTo("hZNC")) { + ih = HidZNC; + } else if (hn.EqualTo("hZPC")) { + ih = HidZPC; + } else if (hn.EqualTo("hZEM")) { + ih = HidZEM; + } else { + LOGF(error, "Not recognized histogram name: %s\n", hname); + return -1; + } + clear(ih); + const int32_t dim = 6; + double x[dim]; + int32_t bins[dim]; + int64_t nb = hs->GetNbins(); + int64_t nn = 0; + LOGF(info, "Histogram %s has %ld bins\n", hname, nb); + double cutl = mTDCCalibConfig->cutLow[ih]; + double cuth = mTDCCalibConfig->cutHigh[ih]; + double contt = 0; + for (int64_t i = 0; i < nb; i++) { + double cont = hs->GetBinContent(i, bins); + if (cont <= 0) { + continue; + } + for (int32_t d = 0; d < dim; ++d) { + x[d] = hs->GetAxis(d)->GetBinCenter(bins[d]); + } + if (TMath::Nint(x[5] - ic) == 0 && x[0] > cutl && x[0] < cuth) { + nn++; + contt += cont; + cumulate(ih, x[0], x[1], x[2], x[3], x[4], cont); + } + } + int ierr = mini(ih); + if (ierr) { + LOGF(error, "FAILED processing RUN2 data for %s ih = %d\n", hname, ih); + } else { + LOGF(info, "Processed RUN2 data for %s ih = %d\n", hname, ih); + replay(ih, hs, ic); + } + LOGF(info, "Trigger class selection %d and %d bins %g events and cuts (%g:%g): %ld\n", ic, nn, contt, cutl, cuth); + return 0; +} + +void TDCCalib::replay(int ih, THnSparse* hs, int ic) +{ + auto* cfg = mTDCCalibConfig; + // clang-format off + if(ih == 0)mHCorr[ih] = std::make_unique("hZNASc","ZNA sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 1)mHCorr[ih] = std::make_unique("hZPASc","ZPA sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 2)mHCorr[ih] = std::make_unique("hZNCSc","ZNC sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 3)mHCorr[ih] = std::make_unique("hZPCSc","ZPC sum corr",cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 4)mHCorr[ih] = std::make_unique("hZEM2c","ZEM2" ,cfg->nb1[ih],cfg->amin1[ih],cfg->amax1[ih]); + if(ih == 0)mCCorr[ih] = std::make_unique("cZNAc","ZNA;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 1)mCCorr[ih] = std::make_unique("cZPAc","ZPA;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 2)mCCorr[ih] = std::make_unique("cZNCc","ZNC;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 3)mCCorr[ih] = std::make_unique("cZPCc","ZPC;TC;SUM corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + if(ih == 4)mCCorr[ih] = std::make_unique("cZEMc","ZEM;ZEM1;ZEM2 corr",cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih],cfg->nb2[ih],cfg->amin2[ih],cfg->amax2[ih]); + // clang-format on + const int32_t dim = 6; + double x[dim]; + int32_t bins[dim]; + int64_t nb = hs->GetNbins(); + double cutl = cfg->cutLow[ih]; + double cuth = cfg->cutHigh[ih]; + double c1 = mPar[ih][1]; + double c2 = mPar[ih][2]; + double c3 = mPar[ih][3]; + double c4 = mPar[ih][4]; + double of = mPar[ih][5]; + for (int64_t i = 0; i < nb; i++) { + double cont = hs->GetBinContent(i, bins); + if (cont <= 0) { + continue; + } + for (int32_t d = 0; d < dim; ++d) { + x[d] = hs->GetAxis(d)->GetBinCenter(bins[d]); + } + if (TMath::Nint(x[5] - ic) == 0 && x[0] > cutl && x[0] < cuth) { + double sumquad = c1 * x[1] + c2 * x[2] + c3 * x[3] + c4 * x[4] + of; + mHCorr[ih]->Fill(sumquad, cont); + mCCorr[ih]->Fill(x[0], sumquad, cont); + } + } +} + +void TDCCalib::clear(int ih) +{ + int ihstart = 0; + int ihstop = NH; + if (ih >= 0 && ih < NH) { + ihstart = ih; + ihstop = ih + 1; + } + for (int32_t ii = ihstart; ii < ihstop; ii++) { + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + mSum[ii][i][j] = 0; + } + } + if (mHUnc[ii]) { + mHUnc[ii]->Reset(); + } + if (mCUnc[ii]) { + mCUnc[ii]->Reset(); + } + } +} + +void TDCCalib::cumulate(int ih, double tc, double t1, double t2, double t3, double t4, double w = 1) +{ + // TODO: add histogram + // TODO: store data to redraw histograms + if (tc < mTDCCalibConfig->cutLow[ih] || tc > mTDCCalibConfig->cutHigh[ih]) { + return; + } + double val[NPAR] = {0, 0, 0, 0, 0, 1}; + val[0] = tc; + val[1] = t1; + val[2] = t2; + val[3] = t3; + val[4] = t4; + for (int32_t i = 0; i < 6; i++) { + for (int32_t j = i; j < 6; j++) { + mSum[ih][i][j] += val[i] * val[j] * w; + } + } + // mSum[ih][5][5] contains the number of analyzed events + double sumquad = val[1] + val[2] + val[3] + val[4]; + mHUnc[ih]->Fill(sumquad, w); + mHUnc[ih + NH]->Fill(val[0]); + mCUnc[ih]->Fill(val[0], sumquad, w); +} + +void TDCCalib::fcn(int& npar, double* gin, double& chi, double* par, int iflag) +{ + chi = 0; + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + chi += (i == 0 ? par[i] : -par[i]) * (j == 0 ? par[j] : -par[j]) * mAdd[i][j]; + } + } +} + +int TDCCalib::mini(int ih) +{ + // Copy to static object and symmetrize matrix + // We use a static function and therefore we can do only one minimization at a time + mMtx.lock(); + for (int32_t i = 0; i < NPAR; i++) { + for (int32_t j = 0; j < NPAR; j++) { + if (j < i) { + mAdd[i][j] = mSum[ih][j][i]; + } else { + mAdd[i][j] = mSum[ih][i][j]; + } + } + } + double arglist[10]; + int ierflg = 0; + double l_bnd = mTDCCalibConfig->l_bnd[ih]; + double u_bnd = mTDCCalibConfig->u_bnd[ih]; + double start = 1.0; + double step = 0.1; + mMn[ih] = std::make_unique(NPAR); + mMn[ih]->SetFCN(fcn); + mMn[ih]->mnparm(0, "c0", 1., 0., 1., 1., ierflg); + mMn[ih]->mnparm(1, "c1", start, step, l_bnd, u_bnd, ierflg); + if (ih == 4) { + // Only two ZEM calorimeters: equalize response + l_bnd = 0; + u_bnd = 0; + start = 0; + step = 0; + } + mMn[ih]->mnparm(2, "c2", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnparm(3, "c3", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnparm(4, "c4", start, step, l_bnd, u_bnd, ierflg); + l_bnd = mTDCCalibConfig->l_bnd_o[ih]; + u_bnd = mTDCCalibConfig->u_bnd_o[ih]; + start = 0; + step = mTDCCalibConfig->step_o[ih]; + mMn[ih]->mnparm(5, "offset", start, step, l_bnd, u_bnd, ierflg); + mMn[ih]->mnexcm("MIGRAD", arglist, 0, ierflg); + for (Int_t i = 0; i < NPAR; i++) { + mMn[ih]->GetParameter(i, mPar[ih][i], mErr[ih][i]); + } + mMtx.unlock(); + return ierflg; +} + +int TDCCalib::write(const std::string fn) +{ + TDirectory* cwd = gDirectory; + TFile* f = new TFile(fn.data(), "recreate"); + if (f->IsZombie()) { + LOG(error) << "Cannot create file: " << fn; + return 1; + } + for (int32_t ih = 0; ih < (2 * NH); ih++) { + if (mHUnc[ih]) { + mHUnc[ih]->Write("", TObject::kOverwrite); + } + } + for (int32_t ih = 0; ih < NH; ih++) { + if (mCUnc[ih]) { + mCUnc[ih]->Write("", TObject::kOverwrite); + } + } + const char* mntit[NH] = {"mZNA", "mZPA", "mZNC", "mZPC", "mZEM"}; + for (int32_t ih = 0; ih < NH; ih++) { + if (mMn[ih]) { + mMn[ih]->Write(mntit[ih], TObject::kOverwrite); + } + } + f->Close(); + cwd->cd(); + return 0; +} diff --git a/Detectors/ZDC/calib/src/ZDCCalibLinkDef.h b/Detectors/ZDC/calib/src/ZDCCalibLinkDef.h new file mode 100644 index 0000000000000..800fef2f11d83 --- /dev/null +++ b/Detectors/ZDC/calib/src/ZDCCalibLinkDef.h @@ -0,0 +1,20 @@ +// 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::zdc::InterCalibConfig + ; + +#endif diff --git a/Detectors/ZDC/calib/src/zdc-intercalib-epn-workflow.cxx b/Detectors/ZDC/calib/src/zdc-intercalib-epn-workflow.cxx new file mode 100644 index 0000000000000..059d50f345a27 --- /dev/null +++ b/Detectors/ZDC/calib/src/zdc-intercalib-epn-workflow.cxx @@ -0,0 +1,42 @@ +// 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. + +#include "Framework/DataProcessorSpec.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "ZDCCalib/InterCalibEPNSpec.h" + +using namespace o2::framework; + +// ------------------------------------------------------------------ +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + WorkflowSpec specs; + specs.emplace_back(o2::zdc::getInterCalibEPNSpec()); + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + return specs; +} diff --git a/Detectors/ZDC/calib/src/zdc-intercalib-workflow.cxx b/Detectors/ZDC/calib/src/zdc-intercalib-workflow.cxx new file mode 100644 index 0000000000000..4e644d77c54a2 --- /dev/null +++ b/Detectors/ZDC/calib/src/zdc-intercalib-workflow.cxx @@ -0,0 +1,36 @@ +// 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. + +#include "Framework/DataProcessorSpec.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "ZDCCalib/InterCalibSpec.h" + +using namespace o2::framework; + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + WorkflowSpec specs; + specs.emplace_back(o2::zdc::getInterCalibSpec()); + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + return specs; +} diff --git a/Detectors/ZDC/macro/CMakeLists.txt b/Detectors/ZDC/macro/CMakeLists.txt index 54c81d703f929..f0c572e17620e 100644 --- a/Detectors/ZDC/macro/CMakeLists.txt +++ b/Detectors/ZDC/macro/CMakeLists.txt @@ -52,3 +52,8 @@ o2_add_test_root_macro(CreateTDCCorr.C PUBLIC_LINK_LIBRARIES O2::ZDCBase O2::ZDCReconstruction O2::SimulationDataFormat O2::CCDB LABELS zdc) + +o2_add_test_root_macro(CreateInterCalibConfig.C + PUBLIC_LINK_LIBRARIES O2::ZDCBase O2::ZDCReconstruction + O2::SimulationDataFormat O2::CCDB + LABELS zdc) diff --git a/Detectors/ZDC/macro/CreateInterCalibConfig.C b/Detectors/ZDC/macro/CreateInterCalibConfig.C new file mode 100644 index 0000000000000..a615a7c093e36 --- /dev/null +++ b/Detectors/ZDC/macro/CreateInterCalibConfig.C @@ -0,0 +1,70 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) + +#include "Framework/Logger.h" +#include "CCDB/CcdbApi.h" +#include "ZDCCalib/InterCalib.h" +#include "ZDCCalib/InterCalibConfig.h" +#include "ZDCBase/Constants.h" +#include +#include +#include + +#endif + +using namespace o2::zdc; +using namespace std; + +void CreateInterCalibConfig(long tmin = 0, long tmax = -1, std::string ccdbHost = "") +{ + + // This object allows for the configuration of the intercalibration of the 4 towers of each calorimeter + // and the calibration of ZEM2 relative to ZEM1 + + InterCalibConfig conf; + + // Enable intercalibration for all calorimeters + // If intercalibration is disabled the intercalibration coefficients + // are copied from previous valid object and flagged as not modified + // ZNA ZPA ZNC ZPC ZEM2 + conf.enable(true, true, true, true, true); + + // The version for this macro considers NO energy calibration, i.e. all coefficients = 1 + // It is necessary to set the binning + conf.setBinning1D(1200, 0, 12000); + conf.setBinning2D(300, 0, 12000); + + conf.setDescription("Simulated, no energy scaling"); + + conf.setMinEntries(100); + + conf.print(); + + o2::ccdb::CcdbApi api; + map metadata; // can be empty + if (ccdbHost.size() == 0 || ccdbHost == "external") { + ccdbHost = "http://alice-ccdb.cern.ch:8080"; + } else if (ccdbHost == "internal") { + ccdbHost = "http://o2-ccdb.internal/"; + } else if (ccdbHost == "test") { + ccdbHost = "http://ccdb-test.cern.ch:8080"; + } else if (ccdbHost == "local") { + ccdbHost = "http://localhost:8080"; + } + api.init(ccdbHost.c_str()); + LOG(info) << "CCDB server: " << api.getURL(); + // store abitrary user object in strongly typed manner + api.storeAsTFileAny(&conf, CCDBPathInterCalibConfig, metadata, tmin, tmax); + + // return conf; +} diff --git a/Detectors/ZDC/macro/CreateSimCondition.C b/Detectors/ZDC/macro/CreateSimCondition.C index fbcf8f26f28ab..48b85852b7ef0 100644 --- a/Detectors/ZDC/macro/CreateSimCondition.C +++ b/Detectors/ZDC/macro/CreateSimCondition.C @@ -32,7 +32,7 @@ void CreateSimCondition(long tmin = 0, long tmax = -1, std::string ccdbHost = "" o2::zdc::SimCondition conf; const float Gains[5] = {15.e-3, 30.e-3, 100.e-3, 15.e-3, 30.e-3}; // gain (response per photoelectron) - const float fudgeFactor = 4.0; // ad hoc factor to tune the gain in the MC + const float fudgeFactor = 5.0; // ad hoc factor to tune the gain in the MC // Source of line shapes, pedestal and noise for each channel // Missing histos for: towers 1-4 of all calorimeters, zem1, all towers of zpc @@ -44,6 +44,38 @@ void CreateSimCondition(long tmin = 0, long tmax = -1, std::string ccdbHost = "" "zpatc", "zpatc", "zpatc", "zpatc", "zpatc", "zpatc" // ZPCC, ZPC1, ZPC2, ZPC3, ZPC4, ZPCS (shape not used) }; + // clang-format off + // Compensation of signal arrival time (ns) + float pos[o2::zdc::NChannels+1]={ +-1, // pos_ZNAC +-1, // pos_ZNA1 +-1, // pos_ZNA2 +-1, // pos_ZNA3 +-1, // pos_ZNA4 +-1, // pos_ZNAS +-1, // pos_ZPAC +-1, // pos_ZPA1 +-1, // pos_ZPA2 +-1, // pos_ZPA3 +-1, // pos_ZPA4 +-1, // pos_ZPAS +-1, // pos_ZEM1 +-1, // pos_ZEM2 +-1, // pos_ZNCC +-1, // pos_ZNC1 +-1, // pos_ZNC2 +-1, // pos_ZNC3 +-1, // pos_ZNC4 +-1, // pos_ZNCS +-1, // pos_ZPCC +-1, // pos_ZPC1 +-1, // pos_ZPC2 +-1, // pos_ZPC3 +-1, // pos_ZPC4 +-1, // pos_ZPCS +0}; // pos_END + // clang-format on + for (int ic = 0; ic < o2::zdc::NChannels; ic++) { auto& channel = conf.channels[ic]; diff --git a/Detectors/ZDC/macro/CreateTowerCalib.C b/Detectors/ZDC/macro/CreateTowerCalib.C index 079d951231f5b..42cb9aec471a0 100644 --- a/Detectors/ZDC/macro/CreateTowerCalib.C +++ b/Detectors/ZDC/macro/CreateTowerCalib.C @@ -32,6 +32,8 @@ void CreateTowerCalib(long tmin = 0, long tmax = -1, std::string ccdbHost = "") // This object allows for the calibration of the 4 towers of each calorimeter // The relative calibration coefficients of towers w.r.t. the common PM // need to be provided + // I.e. energy calibration is the product of Common PM calibration (or ZEM1) + // and tower intercalibration coefficient (or ZEM2) conf.setTowerCalib(IdZNA1, 1.); conf.setTowerCalib(IdZNA2, 1.); @@ -53,6 +55,11 @@ void CreateTowerCalib(long tmin = 0, long tmax = -1, std::string ccdbHost = "") conf.setTowerCalib(IdZPC3, 1.); conf.setTowerCalib(IdZPC4, 1.); + // ZEM2 has special calibration: can be calibrated + // as a common PM and as a tower (equalized to ZEM1) + // The coefficient applied is the product of the two + conf.setTowerCalib(IdZEM2, 1.); + conf.print(); o2::ccdb::CcdbApi api; diff --git a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/DigiReco.h b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/DigiReco.h index 34a304c374601..447cf8ca76740 100644 --- a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/DigiReco.h +++ b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/DigiReco.h @@ -148,7 +148,7 @@ class DigiReco bool mCorrBackground = true; /// Enable TDC pile-up correction bool mCorrBackgroundSet = false; /// TDC pile-up correction set via function call - int correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FTDCVal, float& FTDCAmp, bool isbeg, bool isend); /// Correct TDC single signal + int correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& fTDCVal, float& fTDCAmp, bool isbeg, bool isend); /// Correct TDC single signal int correctTDCBackground(int ibc, int itdc, std::deque& tdc); /// TDC amplitude and time corrections due to pile-up from previous bunches O2_ZDC_DIGIRECO_FLT getPoint(int itdc, int ibeg, int iend, int i); /// Interpolation for current TDC diff --git a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCEnergyParam.h b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCEnergyParam.h index 2a420814a553e..60331d80b896c 100644 --- a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCEnergyParam.h +++ b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCEnergyParam.h @@ -28,7 +28,7 @@ struct ZDCEnergyParam { float energy_calib[NChannels] = {0}; // Energy calibration coefficients void setEnergyCalib(uint32_t ich, float val); float getEnergyCalib(uint32_t ich) const; - void print(); + void print() const; ClassDefNV(ZDCEnergyParam, 1); }; } // namespace zdc diff --git a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTDCParam.h b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTDCParam.h index ae609863e7b82..a5c53b624fa25 100644 --- a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTDCParam.h +++ b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTDCParam.h @@ -32,7 +32,7 @@ struct ZDCTDCParam { float getShift(uint32_t ich) const; void setFactor(uint32_t ich, float val); float getFactor(uint32_t ich) const; - void print(); + void print() const; ClassDefNV(ZDCTDCParam, 2); }; } // namespace zdc diff --git a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTowerParam.h b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTowerParam.h index 7e2eb82b11671..544c26753839b 100644 --- a/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTowerParam.h +++ b/Detectors/ZDC/reconstruction/include/ZDCReconstruction/ZDCTowerParam.h @@ -26,10 +26,16 @@ namespace zdc { struct ZDCTowerParam { float tower_calib[NChannels] = {0}; // Tower calibration coefficients - void setTowerCalib(uint32_t ich, float val); + std::array modified{}; + ZDCTowerParam() + { + modified.fill(false); + } + void clearFlags(); + void setTowerCalib(uint32_t ich, float val, bool ismodified = true); float getTowerCalib(uint32_t ich) const; - void print(); - ClassDefNV(ZDCTowerParam, 1); + void print() const; + ClassDefNV(ZDCTowerParam, 2); }; } // namespace zdc } // namespace o2 diff --git a/Detectors/ZDC/reconstruction/src/DigiReco.cxx b/Detectors/ZDC/reconstruction/src/DigiReco.cxx index fc14f14ca598f..b98b18e508e89 100644 --- a/Detectors/ZDC/reconstruction/src/DigiReco.cxx +++ b/Detectors/ZDC/reconstruction/src/DigiReco.cxx @@ -127,7 +127,8 @@ void DigiReco::init() if (!mTDCParam) { LOG(fatal) << "TDC " << itdc << " missing configuration object and no manual override"; } else { - fval = mTDCParam->getShift(itdc) / FTDCVal; + // Encode TDC shift + fval = mTDCParam->getShift(itdc) / o2::zdc::FTDCVal; } } auto val = std::nearbyint(fval); @@ -139,7 +140,7 @@ void DigiReco::init() } tdc_shift[itdc] = val; if (mVerbosity > DbgZero) { - LOG(info) << itdc << " " << ChannelNames[TDCSignal[itdc]] << " shift= " << tdc_shift[itdc] << " i.s. = " << val * FTDCVal << " ns"; + LOG(info) << itdc << " " << ChannelNames[TDCSignal[itdc]] << " shift= " << tdc_shift[itdc] << " i.s. = " << val * o2::zdc::FTDCVal << " ns"; } } // Amplitude calibration @@ -174,7 +175,7 @@ void DigiReco::init() } } if (mVerbosity > DbgZero) { - LOG(info) << itdc << " " << ChannelNames[TDCSignal[itdc]] << " search= " << ropt.tdc_search[itdc] << " i.s. = " << ropt.tdc_search[itdc] * FTDCVal << " ns"; + LOG(info) << itdc << " " << ChannelNames[TDCSignal[itdc]] << " search= " << ropt.tdc_search[itdc] << " i.s. = " << ropt.tdc_search[itdc] * o2::zdc::FTDCVal << " ns"; } } @@ -1413,7 +1414,7 @@ void DigiReco::assignTDC(int ibun, int ibeg, int iend, int itdc, int tdc, float } // Encode amplitude and assign - auto myamp = TDCAmp / FTDCAmp; + auto myamp = TDCAmp / o2::zdc::FTDCAmp; rec.TDCVal[itdc].push_back(TDCVal); rec.TDCAmp[itdc].push_back(myamp); int& ihit = mReco[ibun].ntdc[itdc]; @@ -1423,7 +1424,7 @@ void DigiReco::assignTDC(int ibun, int ibeg, int iend, int itdc, int tdc, float rec.tdcAmp[itdc][ihit] = myamp; } else { LOG(error) << rec.ir.orbit << "." << rec.ir.bc << " " - << "ibun=" << ibun << " itdc=" << itdc << " tdc=" << tdc << " TDCVal=" << TDCVal * FTDCVal << " TDCAmp=" << TDCAmp * FTDCAmp << " OVERFLOW"; + << "ibun=" << ibun << " itdc=" << itdc << " tdc=" << tdc << " TDCVal=" << TDCVal * o2::zdc::FTDCVal << " TDCAmp=" << TDCAmp * o2::zdc::FTDCAmp << " OVERFLOW"; } #endif // Assign info about pedestal subtration @@ -1439,7 +1440,7 @@ void DigiReco::assignTDC(int ibun, int ibeg, int iend, int itdc, int tdc, float } #ifdef O2_ZDC_DEBUG LOG(info) << __func__ << " itdc=" << itdc << " " << ChannelNames[isig] << " @ ibun=" << ibun << " " << mReco[ibun].ir.orbit << "." << mReco[ibun].ir.bc << " " - << " tdc=" << tdc << " -> " << TDCValCorr << " shift=" << tdc_shift[itdc] << " -> TDCVal=" << TDCVal << "=" << TDCVal * FTDCVal + << " tdc=" << tdc << " -> " << TDCValCorr << " shift=" << tdc_shift[itdc] << " -> TDCVal=" << TDCVal << "=" << TDCVal * o2::zdc::FTDCVal << " mSource[" << isig << "] = " << unsigned(mSource[isig]) << " = " << mOffset[isig] << " amp=" << amp << " -> " << TDCAmpCorr << " calib=" << tdc_calib[itdc] << " -> TDCAmp=" << TDCAmp << "=" << myamp << (ibun == ibeg ? " B" : "") << (ibun == iend ? " E" : ""); @@ -1644,7 +1645,7 @@ void DigiReco::correctTDCPile() //#define O2_ZDC_DEBUG_CORR -int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FTDCVal, float& FTDCAmp, bool isbeg, bool isend) +int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& fTDCVal, float& fTDCAmp, bool isbeg, bool isend) { // Correction of single TDC signals // This function takes into account the position of the signal in the sequence @@ -1653,8 +1654,8 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT constexpr int TDCMax = TDCRange - TSNH - 1; // Fallback is no correction appliead - FTDCVal = TDCVal; - FTDCAmp = TDCAmp; + fTDCVal = TDCVal; + fTDCAmp = TDCAmp; if (mTDCCorr == nullptr) { #ifdef O2_ZDC_DEBUG @@ -1665,8 +1666,8 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT if (isbeg == false && isend == false) { // Mid bunch - FTDCAmp = TDCAmp / mTDCCorr->mAFMidC[itdc][0]; - FTDCVal = TDCVal - mTDCCorr->mTSMidC[itdc][0]; + fTDCAmp = TDCAmp / mTDCCorr->mAFMidC[itdc][0]; + fTDCVal = TDCVal - mTDCCorr->mTSMidC[itdc][0]; } else if (isbeg == true) { { auto p0 = mTDCCorr->mTSBegC[itdc][0]; @@ -1675,9 +1676,9 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT auto diff = TDCVal - p0; auto p2 = mTDCCorr->mTSBegC[itdc][2]; auto p3 = mTDCCorr->mTSBegC[itdc][3]; - FTDCVal = TDCVal - (p1 + p2 * diff + p3 * diff * diff); + fTDCVal = TDCVal - (p1 + p2 * diff + p3 * diff * diff); } else { - FTDCVal = TDCVal - p1; + fTDCVal = TDCVal - p1; } } { @@ -1687,9 +1688,9 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT auto diff = TDCVal - p0; auto p2 = mTDCCorr->mAFBegC[itdc][2]; auto p3 = mTDCCorr->mAFBegC[itdc][3]; - FTDCAmp = TDCAmp / (p1 + p2 * diff + p3 * diff * diff); + fTDCAmp = TDCAmp / (p1 + p2 * diff + p3 * diff * diff); } else { - FTDCAmp = TDCAmp / p1; + fTDCAmp = TDCAmp / p1; } } } else if (isend == true) { @@ -1700,9 +1701,9 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT auto diff = TDCVal - p0; auto p2 = mTDCCorr->mTSEndC[itdc][2]; auto p3 = mTDCCorr->mTSEndC[itdc][3]; - FTDCVal = TDCVal - (p1 + p2 * diff + p3 * diff * diff); + fTDCVal = TDCVal - (p1 + p2 * diff + p3 * diff * diff); } else { - FTDCVal = TDCVal - p1; + fTDCVal = TDCVal - p1; } } { @@ -1712,9 +1713,9 @@ int DigiReco::correctTDCSignal(int itdc, int16_t TDCVal, float TDCAmp, float& FT auto diff = TDCVal - p0; auto p2 = mTDCCorr->mAFEndC[itdc][2]; auto p3 = mTDCCorr->mAFEndC[itdc][3]; - FTDCAmp = TDCAmp / (p1 + p2 * diff + p3 * diff * diff); + fTDCAmp = TDCAmp / (p1 + p2 * diff + p3 * diff * diff); } else { - FTDCAmp = TDCAmp / p1; + fTDCAmp = TDCAmp / p1; } } } else { diff --git a/Detectors/ZDC/reconstruction/src/ZDCEnergyParam.cxx b/Detectors/ZDC/reconstruction/src/ZDCEnergyParam.cxx index ef67712c0f4d7..eecf797b2df24 100644 --- a/Detectors/ZDC/reconstruction/src/ZDCEnergyParam.cxx +++ b/Detectors/ZDC/reconstruction/src/ZDCEnergyParam.cxx @@ -43,7 +43,7 @@ float ZDCEnergyParam::getEnergyCalib(uint32_t ich) const } } -void ZDCEnergyParam::print() +void ZDCEnergyParam::print() const { for (Int_t ich = 0; ich < NChannels; ich++) { if (energy_calib[ich] > 0) { diff --git a/Detectors/ZDC/reconstruction/src/ZDCTDCParam.cxx b/Detectors/ZDC/reconstruction/src/ZDCTDCParam.cxx index 8081d6601e649..653898b088b32 100644 --- a/Detectors/ZDC/reconstruction/src/ZDCTDCParam.cxx +++ b/Detectors/ZDC/reconstruction/src/ZDCTDCParam.cxx @@ -52,7 +52,7 @@ float ZDCTDCParam::getFactor(uint32_t ich) const } } -void ZDCTDCParam::print() +void ZDCTDCParam::print() const { for (int itdc = 0; itdc < o2::zdc::NTDCChannels; itdc++) { LOG(info) << ChannelNames[TDCSignal[itdc]] << " shift = " << tdc_shift[itdc] << " ns factor = " << tdc_calib[itdc]; diff --git a/Detectors/ZDC/reconstruction/src/ZDCTowerParam.cxx b/Detectors/ZDC/reconstruction/src/ZDCTowerParam.cxx index d2c17e4ff64dd..3a95ef79a61e3 100644 --- a/Detectors/ZDC/reconstruction/src/ZDCTowerParam.cxx +++ b/Detectors/ZDC/reconstruction/src/ZDCTowerParam.cxx @@ -14,7 +14,12 @@ using namespace o2::zdc; -void ZDCTowerParam::setTowerCalib(uint32_t ich, float val) +void ZDCTowerParam::clearFlags() +{ + modified.fill(false); +} + +void ZDCTowerParam::setTowerCalib(uint32_t ich, float val, bool ismodified) { bool in_list = false; for (int il = 0; il < ChTowerCalib.size(); il++) { @@ -25,6 +30,7 @@ void ZDCTowerParam::setTowerCalib(uint32_t ich, float val) } if (in_list) { tower_calib[ich] = val; + modified[ich] = ismodified; } else { LOG(fatal) << __func__ << " channel " << ich << " not in allowed range"; for (int il = 0; il < ChTowerCalib.size(); il++) { @@ -43,11 +49,11 @@ float ZDCTowerParam::getTowerCalib(uint32_t ich) const } } -void ZDCTowerParam::print() +void ZDCTowerParam::print() const { for (Int_t ich = 0; ich < NChannels; ich++) { if (tower_calib[ich] > 0) { - LOG(info) << ChannelNames[ich] << " calibration factor = " << tower_calib[ich]; + LOG(info) << ChannelNames[ich] << (modified[ich] ? " NEW" : " OLD") << " calibration factor = " << tower_calib[ich]; } } } diff --git a/Detectors/ZDC/workflow/CMakeLists.txt b/Detectors/ZDC/workflow/CMakeLists.txt index ad4a5f4abd63e..727f9b57996be 100644 --- a/Detectors/ZDC/workflow/CMakeLists.txt +++ b/Detectors/ZDC/workflow/CMakeLists.txt @@ -16,13 +16,15 @@ o2_add_library(ZDCWorkflow src/ZDCDataReaderDPLSpec.cxx src/ZDCDigitWriterDPLSpec.cxx src/ZDCRecoWriterDPLSpec.cxx - src/RecEventReaderSpec.cxx + src/RecEventReaderSpec.cxx src/RecoWorkflow.cxx src/DigitRecoSpec.cxx src/DigitReaderSpec.cxx + src/RecoReaderSpec.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DataFormatsZDC O2::ZDCRaw + O2::ZDCCalib O2::SimulationDataFormat O2::DPLUtils O2::ZDCReconstruction @@ -43,6 +45,11 @@ o2_add_executable(digits-read SOURCES src/digits-reader-workflow.cxx PUBLIC_LINK_LIBRARIES O2::ZDCWorkflow) +o2_add_executable(reco-read + COMPONENT_NAME zdc + SOURCES src/reco-reader-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::ZDCWorkflow) + o2_add_executable(digits-reco COMPONENT_NAME zdc SOURCES src/zdc-reco-workflow.cxx diff --git a/Detectors/ZDC/workflow/include/ZDCWorkflow/RecoReaderSpec.h b/Detectors/ZDC/workflow/include/ZDCWorkflow/RecoReaderSpec.h new file mode 100644 index 0000000000000..95b17cf00e392 --- /dev/null +++ b/Detectors/ZDC/workflow/include/ZDCWorkflow/RecoReaderSpec.h @@ -0,0 +1,45 @@ +// 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_ZDC_RECOREADERSPEC_H +#define O2_ZDC_RECOREADERSPEC_H + +#include "TFile.h" +#include "TTree.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +namespace o2 +{ +namespace zdc +{ + +class RecoReader : public framework::Task +{ + public: + RecoReader() {} + ~RecoReader() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + + private: + std::unique_ptr mTree; + std::unique_ptr mFile; +}; + +/// create a processor spec +/// read reconstructed ZDC data from a root file +framework::DataProcessorSpec getRecoReaderSpec(); + +} // end namespace zdc +} // end namespace o2 + +#endif // O2_ZDC_RECOREADERSPEC_H diff --git a/Detectors/ZDC/workflow/src/RecoReaderSpec.cxx b/Detectors/ZDC/workflow/src/RecoReaderSpec.cxx new file mode 100644 index 0000000000000..396e9fb38bda0 --- /dev/null +++ b/Detectors/ZDC/workflow/src/RecoReaderSpec.cxx @@ -0,0 +1,96 @@ +// 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 RecoReaderSpec.cxx + +#include + +#include "TTree.h" + +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/Logger.h" +#include "ZDCWorkflow/RecoReaderSpec.h" +#include "DataFormatsZDC/OrbitData.h" +#include "DataFormatsZDC/BCData.h" +#include "DataFormatsZDC/ChannelData.h" +#include "DataFormatsZDC/RecEvent.h" +#include "CommonUtils/NameConf.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace zdc +{ + +void RecoReader::init(InitContext& ic) +{ + auto filename = o2::utils::Str::concat_string(o2::utils::Str::rectifyDirectory(ic.options().get("input-dir")), + ic.options().get("zdc-reco-infile")); + mFile = std::make_unique(filename.c_str()); + if (!mFile->IsOpen()) { + LOG(error) << "Cannot open the " << filename.c_str() << " file !"; + throw std::runtime_error("cannot open input digits file"); + } + mTree.reset((TTree*)mFile->Get("o2rec")); + if (!mTree) { + LOG(error) << "Did not find o2sim tree in " << filename.c_str(); + throw std::runtime_error("Did not find o2rec tree in ZDC reco file"); + } +} + +void RecoReader::run(ProcessingContext& pc) +{ + + std::vector RecBC, *RecBCPtr = &RecBC; + std::vector Energy, *EnergyPtr = &Energy; + std::vector TDCData, *TDCDataPtr = &TDCData; + std::vector Info, *InfoPtr = &Info; + + mTree->SetBranchAddress("ZDCRecBC", &RecBCPtr); + mTree->SetBranchAddress("ZDCRecE", &EnergyPtr); + mTree->SetBranchAddress("ZDCRecTDC", &TDCDataPtr); + mTree->SetBranchAddress("ZDCRecInfo", &InfoPtr); + + auto ent = mTree->GetReadEntry() + 1; + assert(ent < mTree->GetEntries()); // this should not happen + mTree->GetEntry(ent); + LOG(info) << "ZDCRecoReader pushed " << RecBC.size() << " b.c. " << Energy.size() << " Energies " << TDCData.size() << " TDCs " << Info.size() << " Infos"; + pc.outputs().snapshot(Output{"ZDC", "BCREC", 0, Lifetime::Timeframe}, RecBC); + pc.outputs().snapshot(Output{"ZDC", "ENERGY", 0, Lifetime::Timeframe}, Energy); + pc.outputs().snapshot(Output{"ZDC", "TDCDATA", 0, Lifetime::Timeframe}, TDCData); + pc.outputs().snapshot(Output{"ZDC", "INFO", 0, Lifetime::Timeframe}, Info); + if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } +} + +DataProcessorSpec getRecoReaderSpec() +{ + std::vector outputs; + outputs.emplace_back("ZDC", "BCREC", 0, Lifetime::Timeframe); + outputs.emplace_back("ZDC", "ENERGY", 0, Lifetime::Timeframe); + outputs.emplace_back("ZDC", "TDCDATA", 0, Lifetime::Timeframe); + outputs.emplace_back("ZDC", "INFO", 0, Lifetime::Timeframe); + return DataProcessorSpec{ + "zdc-reco-reader", + Inputs{}, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{ + {"zdc-reco-infile", VariantType::String, "zdcreco.root", {"Name of the input file"}}, + {"input-dir", VariantType::String, "none", {"Input directory"}}}}; +} + +} // namespace zdc +} // namespace o2 diff --git a/Detectors/ZDC/workflow/src/reco-reader-workflow.cxx b/Detectors/ZDC/workflow/src/reco-reader-workflow.cxx new file mode 100644 index 0000000000000..94ee7f0e8df87 --- /dev/null +++ b/Detectors/ZDC/workflow/src/reco-reader-workflow.cxx @@ -0,0 +1,55 @@ +// 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 reco-reader-workflow.cxx +/// \brief Implementation of ZDC reco reader +/// +/// \author ruben,shahoian@cern.ch, pietro.cortese@cern.ch + +#include "Framework/CallbackService.h" +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/Task.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "ZDCWorkflow/RecoReaderSpec.h" +#include "CommonUtils/ConfigurableParam.h" + +using namespace o2::framework; + +// ------------------------------------------------------------------ +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + std::string keyvaluehelp("Semicolon separated key=value strings ..."); + workflowOptions.push_back(ConfigParamSpec{"configKeyValues", VariantType::String, "", {keyvaluehelp}}); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + WorkflowSpec specs; + + DataProcessorSpec producer = o2::zdc::getRecoReaderSpec(); + specs.push_back(producer); + + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + + return specs; +}