From 40783ece1038e3551b1bd04ea4b0e4548a7c7b31 Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Wed, 22 Nov 2023 23:07:17 +0100 Subject: [PATCH 1/2] [MCH] new ROF clustering method --- .../MUON/MCH/TimeClustering/CMakeLists.txt | 4 +- .../ROFTimeClusterFinderV2.h | 116 ++++ .../TimeClusterFinderSpecV2.h | 37 ++ .../TimeClusterizerParamV2.h | 46 ++ .../src/MCHTimeClusteringLinkDef.h | 2 + .../src/ROFTimeClusterFinder.cxx | 4 +- .../src/ROFTimeClusterFinderV2.cxx | 596 ++++++++++++++++++ .../src/TimeClusterFinderSpecV2.cxx | 346 ++++++++++ .../src/TimeClusterizerParamV2.cxx | 15 + .../src/digits-to-timeclusters-workflow.cxx | 11 + .../MUON/MCH/Workflow/src/reco-workflow.cxx | 12 + 11 files changed, 1185 insertions(+), 4 deletions(-) create mode 100644 Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h create mode 100644 Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterFinderSpecV2.h create mode 100644 Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterizerParamV2.h create mode 100644 Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx create mode 100644 Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx create mode 100644 Detectors/MUON/MCH/TimeClustering/src/TimeClusterizerParamV2.cxx diff --git a/Detectors/MUON/MCH/TimeClustering/CMakeLists.txt b/Detectors/MUON/MCH/TimeClustering/CMakeLists.txt index 28fada0e94907..ff94aa3b54bed 100644 --- a/Detectors/MUON/MCH/TimeClustering/CMakeLists.txt +++ b/Detectors/MUON/MCH/TimeClustering/CMakeLists.txt @@ -10,10 +10,10 @@ # or submit itself to any jurisdiction. o2_add_library(MCHTimeClustering - SOURCES src/ROFTimeClusterFinder.cxx src/TimeClusterizerParam.cxx src/TimeClusterFinderSpec.cxx + SOURCES src/ROFTimeClusterFinder.cxx src/TimeClusterizerParam.cxx src/TimeClusterFinderSpec.cxx src/ROFTimeClusterFinderV2.cxx src/TimeClusterizerParamV2.cxx src/TimeClusterFinderSpecV2.cxx PUBLIC_LINK_LIBRARIES O2::MCHBase O2::Framework O2::MCHDigitFiltering O2::MCHROFFiltering) -o2_target_root_dictionary(MCHTimeClustering HEADERS include/MCHTimeClustering/TimeClusterizerParam.h) +o2_target_root_dictionary(MCHTimeClustering HEADERS include/MCHTimeClustering/TimeClusterizerParam.h include/MCHTimeClustering/TimeClusterizerParamV2.h) o2_add_executable( digits-to-timeclusters-workflow diff --git a/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h new file mode 100644 index 0000000000000..9cf05e6b32cc9 --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h @@ -0,0 +1,116 @@ +// 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 ROFTimeClusterFinder.h +/// \brief Class to group the fired pads according to their time stamp +/// +/// \author Andrea Ferrero, CEA + +#ifndef O2_MCH_ROFTIMECLUSTERFINDERV2_H_ +#define O2_MCH_ROFTIMECLUSTERFINDERV2_H_ + +#include +#include +#include +#include + +#include "DataFormatsMCH/Digit.h" +#include "DataFormatsMCH/ROFRecord.h" +#include "MCHDigitFiltering/DigitFilter.h" + +namespace o2 +{ +namespace mch +{ + +class ROFTimeClusterFinderV2 +{ + public: + using ROFVector = std::vector; + + ROFTimeClusterFinderV2(gsl::span rofs, + gsl::span digits, + uint32_t timeClusterSize, + uint32_t peakSearchWindow, + uint32_t nBins, + uint32_t nDigitsMin, + bool improvePeakSearch, + bool mergeROFs, + bool debug); + ~ROFTimeClusterFinderV2() = default; + + /// process the digit ROFs and create the time clusters + void process(); + + /// return the vector of time-cluster ROFs + const ROFVector& getROFRecords() const { return mOutputROFs; } + + /// stores the output ROFs into a flat memory buffer + char* saveROFRsToBuffer(size_t& bufSize); + + /// debugging output + void dumpInputROFs(); + void dumpOutputROFs(); + + private: + /// structure representing one bin in the time histogram used for the peak search algorithm + struct TimeBin { + int32_t mFirstIdx{-1}; ///< index of the first digit ROF in the bin + int32_t mLastIdx{-1}; ///< index of the last digit ROF in the bin + uint32_t mNDigitsPS{0}; ///< number of digits in the bin for the peak search step + + TimeBin() = default; + + bool empty() { return mNDigitsPS == 0; } + + bool operator<(const TimeBin& other) { return (mNDigitsPS < other.mNDigitsPS); } + bool operator>(const TimeBin& other) { return (mNDigitsPS > other.mNDigitsPS); } + bool operator<=(const TimeBin& other) { return (mNDigitsPS <= other.mNDigitsPS); } + bool operator>=(const TimeBin& other) { return (mNDigitsPS >= other.mNDigitsPS); } + }; + + // peak search parameters + uint32_t mTimeClusterSize; ///< maximum size of one time cluster, in bunch crossings + uint32_t mPeakSearchWindow; ///< width of the peak search window, in BC units + uint32_t mPeakSearchNbins; ///< number of time bins for the peak search algorithm (must be an odd number >= 3) + uint32_t mPeakSearchNDigitsMin; ///< minimum number of digits for peak candidates + uint32_t mBinWidth; ///< width of one time bin in the peak search algorithm, in bunch crossings + uint32_t mNbinsInOneTF; ///< maximum number of peak search bins in one time frame + DigitFilter mIsGoodDigit; ///< function to select only digits that are likely signal + bool mImprovePeakSearch; ///< whether to only use signal-like digits in the peak search + bool mMergeROFs; ///< whether to merge consecutive ROFs + + /// fill the time histogram for the peak search algorithm + void initTimeBins(); + /// search for the next peak in the time histogram + int32_t getNextPeak(); + /// search for the end of a given peak + int32_t getPeakEnd(int32_t peak, int32_t max); + /// create an output ROF containing all the digits in the [firstBin,lastBin] range of the time histogram + void storeROF(int32_t firstBin, int32_t lastBin, int32_t peak); + + std::vector mTimeBins; ///< time histogram for the peak search algorithm + int32_t mLastSavedTimeBin; ///< index of the last bin that has been stored in the output ROFs + int32_t mLastSavedInputRof; ///< index of the last bin that has been stored in the output ROFs + int32_t mLastPeak; ///< index of the last found peak + int32_t mLastPeakEnd; ///< index of the end bin of last found peak, used in the merging step + + gsl::span mInputROFs; ///< input digit ROFs + gsl::span mDigits; ///< input digits + ROFVector mOutputROFs{}; ///< output time cluster ROFs + std::vector mRofHasPeak; + bool mDebug{false}; ///< be more verbose +}; + +} // namespace mch +} // namespace o2 + +#endif // O2_MCH_ROFTIMECLUSTERFINDERV2_H_ diff --git a/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterFinderSpecV2.h b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterFinderSpecV2.h new file mode 100644 index 0000000000000..0a29e99e98f94 --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterFinderSpecV2.h @@ -0,0 +1,37 @@ +// 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 TimeClusterFinderSpecV2.h +/// \brief Definition of a data processor to run the time clusterizer +/// +/// \author Andrea Ferrero, CEA + +#ifndef O2_MCH_TIMECLUSTERFINDERSPECV2_H_ +#define O2_MCH_TIMECLUSTERFINDERSPECV2_H_ + +#include "Framework/DataProcessorSpec.h" + +namespace o2 +{ +namespace mch +{ + +o2::framework::DataProcessorSpec + getTimeClusterFinderSpecV2(const char* specName = "mch-time-cluster-finder-new", + std::string_view inputDigitDataDescription = "F-DIGITS", + std::string_view inputDigitRofDataDescription = "F-DIGITROFS", + std::string_view outputDigitRofDataDescription = "TC-F-DIGITROFS", + std::string_view inputIRFrameDataDescription = "ITS/IRFRAMES"); + +} // end namespace mch +} // end namespace o2 + +#endif // O2_MCH_TIMECLUSTERFINDERSPECV2_H_ diff --git a/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterizerParamV2.h b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterizerParamV2.h new file mode 100644 index 0000000000000..995221b32200c --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/TimeClusterizerParamV2.h @@ -0,0 +1,46 @@ +// 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_MCH_TIMECLUSTERING_TIMECLUSTERIZER_PARAM_V2_H_ +#define O2_MCH_TIMECLUSTERING_TIMECLUSTERIZER_PARAM_V2_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" + +namespace o2::mch +{ + +/** + * @class TimeClusterizerParamV2 + * @brief Configurable parameters for the time clustering + */ +struct TimeClusterizerParamV2 : public o2::conf::ConfigurableParamHelper { + + bool onlyTrackable = true; ///< only output ROFs that match the trackable condition @see MCHROFFiltering/TrackableFilter + + int maxClusterWidth = 1000 / 25; ///< initial width of time clusters, in BC units (default 1us) + int peakSearchWindow = 20; ///< width of the peak search window, in BC units + int peakSearchNbins = 5; ///< number of time bins for the peak search algorithm (must be an odd number >= 3) + int peakSearchNDigitsMin = 10; ///< minimum number of digits for peak candidates + bool peakSearchSignalOnly = true; ///< only use signal-like hits in peak search + bool mergeROFs = true; ///< whether to merge consecutive ROFs + int mergeGapMax = 11; ///< maximum allowed gap between ROFs, above which no merging is done + int minDigitsPerROF = 0; ///< minimum number of digits per ROF (below that threshold ROF is discarded) + bool irFramesOnly = false; ///< only output ROFs that overlap one of the IRFrames (provided externally, e.g. by ITS) @see MCHROFFiltering/IRFrameFilter + + float rofRejectionFraction = 0; ///< fraction of output (i.e. time-clusterized) ROFs to discard. If 0 (default) keep them all. WARNING: use a non zero value only at Pt2 for sync reco, if needed. + + O2ParamDef(TimeClusterizerParamV2, "MCHTimeClusterizerV2"); +}; + +} // namespace o2::mch + +#endif diff --git a/Detectors/MUON/MCH/TimeClustering/src/MCHTimeClusteringLinkDef.h b/Detectors/MUON/MCH/TimeClustering/src/MCHTimeClusteringLinkDef.h index d8969288004d6..c5221b9e68bbb 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/MCHTimeClusteringLinkDef.h +++ b/Detectors/MUON/MCH/TimeClustering/src/MCHTimeClusteringLinkDef.h @@ -17,5 +17,7 @@ #pragma link C++ class o2::mch::TimeClusterizerParam + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::mch::TimeClusterizerParam> + ; +#pragma link C++ class o2::mch::TimeClusterizerParamV2 + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::mch::TimeClusterizerParamV2> + ; #endif diff --git a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx index b8099852c8215..7187aaea19243 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx @@ -235,10 +235,10 @@ void ROFTimeClusterFinder::storeROF(int32_t firstBin, int32_t lastBin) void ROFTimeClusterFinder::process() { - if (mDebug) { + //if (mDebug) { std::cout << "\n\n==================\n[ROFTimeClusterFinder] processing new TF\n" << std::endl; - } + //} initTimeBins(); mOutputROFs.clear(); diff --git a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx new file mode 100644 index 0000000000000..0dc3e99e050e3 --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx @@ -0,0 +1,596 @@ +// 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 "MCHTimeClustering/ROFTimeClusterFinderV2.h" + +#include +#include +#include "Framework/Logger.h" + +namespace o2 +{ +namespace mch +{ + +using namespace std; + +//_________________________________________________________________________________________________ + +ROFTimeClusterFinderV2::ROFTimeClusterFinderV2(gsl::span rofs, + gsl::span digits, + uint32_t timeClusterSize, + uint32_t peakSearchWindow, + uint32_t nBins, + uint32_t nDigitsMin, + bool improvePeakSearch, + bool mergeROFs, + bool debug) + : mTimeClusterSize(timeClusterSize), + mPeakSearchWindow(peakSearchWindow), + mPeakSearchNbins(nBins), + mBinWidth(peakSearchWindow / nBins), + mPeakSearchNDigitsMin(nDigitsMin), + mNbinsInOneTF(0), + mIsGoodDigit(createDigitFilter(20, true, true)), + mImprovePeakSearch(improvePeakSearch), + mMergeROFs(mergeROFs), + mTimeBins{}, + mLastSavedTimeBin(-1), + mLastPeak(-1), + mInputROFs(rofs), + mDigits(digits), + mOutputROFs{}, + mDebug(debug) +{ +} + +//_________________________________________________________________________________________________ + +void ROFTimeClusterFinderV2::initTimeBins() +{ + static constexpr uint32_t maxNTimeBins = 1e7; + + //std::cout << "[TOTO] mTimeClusterSize " << mTimeClusterSize << " mPeakSearchNbins " << mPeakSearchNbins << " mBinWidth " << mBinWidth << std::endl; + + mTimeBins.clear(); + mNbinsInOneTF = 0; + + if (mInputROFs.empty()) { + return; + } + + // initialize the time bins vector + o2::InteractionRecord mFirstIR = mInputROFs.front().getBCData(); + auto tfSize = mInputROFs.back().getBCData().differenceInBC(mFirstIR); + mNbinsInOneTF = tfSize / mBinWidth + 1; + if (mNbinsInOneTF > maxNTimeBins) { + LOGP(alarm, "Number of time bins exceeding the limit"); + mNbinsInOneTF = maxNTimeBins; + } + mTimeBins.resize(mNbinsInOneTF); + + // store the number of digits in each bin + int64_t previousROFbc = -1; + for (size_t iRof = 0; iRof < mInputROFs.size(); iRof++) { + const auto& rof = mInputROFs[iRof]; + const auto& ir = rof.getBCData(); + auto rofBc = ir.differenceInBC(mFirstIR); + auto binIdx = rofBc / mBinWidth; + + // sanity checks: ROFs must be ordered in time and not be empty + if (rofBc <= previousROFbc || rofBc > tfSize) { + LOGP(alarm, "Wrong ROF ordering"); + break; + } + previousROFbc = rofBc; + if (rof.getNEntries() < 1) { + LOGP(alarm, "Empty ROF"); + break; + } + + // stop here if the number of bins exceeds the limit + if (binIdx >= mNbinsInOneTF) { + break; + } + + auto& timeBin = mTimeBins[binIdx]; + + if (timeBin.mFirstIdx < 0) { + timeBin.mFirstIdx = iRof; + } + timeBin.mLastIdx = iRof; + + int nDigitsPS = 0; + if (mImprovePeakSearch) { + auto rofDigits = mDigits.subspan(rof.getFirstIdx(), rof.getNEntries()); + for (auto& digit : rofDigits) { + if (mIsGoodDigit(digit)) { + nDigitsPS += 1; + } + } + } else { + nDigitsPS = rof.getNEntries(); + } + timeBin.mNDigitsPS += nDigitsPS; + } + + if (mDebug) { + std::cout << "Peak search histogram:" << std::endl; + for (int32_t i = 0; i < mNbinsInOneTF; i++) { + if (mTimeBins[i].mFirstIdx >= 0) { + auto nDigits = mInputROFs[mTimeBins[i].mLastIdx].getLastIdx() - mInputROFs[mTimeBins[i].mFirstIdx].getFirstIdx() + 1; + std::cout << fmt::format("bin {}: {}/{}", i, mTimeBins[i].mNDigitsPS, nDigits) << std::endl; + } + } + } +} + +//_________________________________________________________________________________________________ + +static int targetOrbit = 376784078; + +int32_t ROFTimeClusterFinderV2::getNextPeak() +{ + // number of time bins before and after the peak seeds that are considered in the peak search + int32_t sPadding = mPeakSearchNbins / 2; + + if (mDebug) { + std::cout << "Searching peak from " << mLastSavedTimeBin + 1 << std::endl; + } + + // loop over the bins and search for local maxima + // a local maxima is defined as a bin that is higher than all the surrounding bins, in a window given by mPeakSearchNbins + for (int32_t i = mLastSavedTimeBin + 1; i < mNbinsInOneTF; i++) { + auto& peak = mTimeBins[i]; + // only time bins with a number of digits higher or equal to the required minimum are considered as peak seeds + if (peak.empty() || peak.mNDigitsPS < mPeakSearchNDigitsMin) { + continue; + } + + int peakOrbit = mInputROFs[peak.mFirstIdx].getBCData().orbit; + if (peakOrbit == targetOrbit) { + int peakBc = mInputROFs[peak.mFirstIdx].getBCData().bc; + std::cout << "[TOTO] peak seed " << i << " " << peak.mNDigitsPS << " bc " << peakBc << " padding " << sPadding << std::endl; + } + + + bool found{true}; + // the peak must be strictly higher than previous bins + for (int j = i - sPadding; j < i; j++) { + if (j > mLastSavedTimeBin && peakOrbit == targetOrbit) { + int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; + int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; + std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << std::endl; + } + if (j > mLastSavedTimeBin && peak <= mTimeBins[j]) { + found = false; + break; + } + } + // the peak must be higher than or equal to next bins + for (int j = i + 1; j <= i + sPadding; j++) { + if (j < mNbinsInOneTF && peakOrbit == targetOrbit) { + int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; + int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; + std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << std::endl; + } + if (j < mNbinsInOneTF && peak < mTimeBins[j]) { + found = false; + break; + } + } + if (peakOrbit == targetOrbit) { + std::cout << "[TOTO] peak found " << found << std::endl; + } + + if (!found) { + continue; + } + + if (mDebug) { + auto nDigits = mInputROFs[peak.mLastIdx].getLastIdx() - mInputROFs[peak.mFirstIdx].getFirstIdx() + 1; + std::cout << fmt::format("new peak found at bin {}, entries = {}/{}", i, peak.mNDigitsPS, nDigits) << std::endl; + } + + return i; + } + return -1; +} + +//_________________________________________________________________________________________________ + + +int32_t ROFTimeClusterFinderV2::getPeakEnd(int32_t peak, int32_t max) +{ + // number of time bins before and after the peak seeds that are considered in the peak search + int32_t sPadding = mPeakSearchNbins / 2; + + if (peak < 0) { + return -1; + } + + if (mDebug) { + std::cout << "Searching peak end from " << peak + 1 << std::endl; + } + + int peakOrbit = mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().orbit; + if (peakOrbit == targetOrbit) { + //int peakBc = mInputROFs[peak.mFirstIdx].getBCData().bc; + //std::cout << "[TOTO] peak tail " << i << " " << peak.mNDigitsPS << " bc " << peakBc << " padding " << sPadding << std::endl; + } + + // loop over the bins and search for local maxima + // a local maxima is defined as a bin that is higher than all the surrounding bins, in a window given by mPeakSearchNbins + for (int32_t i = peak + 1; i < max; i++) { + if (i >= mNbinsInOneTF) { + return i - 1; + } + + auto& tail = mTimeBins[i]; + // only time bins with a number of digits higher or equal to the required minimum are considered as peak seeds + if (tail.empty()) { + if (peakOrbit == targetOrbit) { + int bc = (mTimeBins[peak].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().bc; + std::cout << "[TOTO] found peak tail " << i << " orbit " << peakOrbit << " bc " << bc << " nDigits " << 0 << std::endl; + } + return i; + } + + // the tail must be higher than or equal to the next bins + for (int j = i + 1; j < i + sPadding; j++) { + //if (peakOrbit == targetOrbit) { + // int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; + // int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; + // std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << std::endl; + //} + if (j < mNbinsInOneTF && tail < mTimeBins[j]) { + if (peakOrbit == targetOrbit) { + int orbit = (mTimeBins[i].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[i].mFirstIdx].getBCData().orbit; + int bc = (mTimeBins[i].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[i].mFirstIdx].getBCData().bc; + std::cout << "[TOTO] found peak tail " << i << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[i].mNDigitsPS << std::endl; + } + return i; + } + } + } + return max; +} + +//_________________________________________________________________________________________________ + +void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t peak) +{ + if (mDebug) { + std::cout << fmt::format("Storing ROF from bin range [{},{}] with peak {}", firstBin, lastBin, peak) << std::endl; + } + + bool hasPeak = peak >= 0; + + if (firstBin < 0) { + firstBin = 0; + } + if (lastBin >= mNbinsInOneTF) { + lastBin = mNbinsInOneTF - 1; + } + + // check if part of this range can be appended to the last saved ROF + int firstBinNew = firstBin; + if (!hasPeak && !mOutputROFs.empty() && (mOutputROFs.back().getBCWidth() < mTimeClusterSize)) { + // this ROF is a dummy one without peak, so first we prolong the last saved one if it is smaller than the maximum allowed + auto& prevRof = mOutputROFs.back(); + int32_t bcWidth = prevRof.getBCWidth(); + int32_t nDigits = prevRof.getNEntries(); + // loop over bins in the range to find those that can be attached to the last saved ROF + for (int32_t j = firstBin; j <= lastBin; j++) { + auto& timeBin = mTimeBins[j]; + if (timeBin.mFirstIdx < 0) { + continue; + } + + const auto& rof = mInputROFs[timeBin.mFirstIdx]; + auto newWidth = rof.getBCData().differenceInBC(prevRof.getBCData()) + rof.getBCWidth(); + if (prevRof.getBCData().orbit == targetOrbit) { + std::cout << "[TOTO] trying to extend " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc + << " with rof " << rof.getBCData().orbit << " / " << rof.getBCData().bc << " newWidth " << newWidth << std::endl; + } + // we stop if the total widht after merging would exceed the maximum allowd + if (newWidth > mTimeClusterSize) { + break; + } + + // update ROF parameters + bcWidth = newWidth; + nDigits += rof.getNEntries(); + + // remove the current bin from the range and update the last saved bin + firstBinNew = j + 1; + mLastSavedTimeBin = j; + } + if (firstBinNew > firstBin) { + prevRof = ROFRecord(prevRof.getBCData(), prevRof.getFirstIdx(), nDigits, bcWidth); + if (prevRof.getBCData().orbit == targetOrbit) { + std::cout << "[TOTO] ROF " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc + << " extended to " << prevRof.getBCWidth() << std::endl; + } + } + } + + firstBin = firstBinNew; + + int32_t rofFirstIdx{-1}; + int32_t rofLastIdx{-1}; + for (int32_t j = firstBin; j <= lastBin; j++) { + auto& timeBin = mTimeBins[j]; + if (timeBin.mFirstIdx < 0) { + continue; + } + + if (rofFirstIdx < 0) { + rofFirstIdx = timeBin.mFirstIdx; + }; + rofLastIdx = timeBin.mLastIdx; + + if (mDebug) { + auto size = timeBin.mLastIdx - timeBin.mFirstIdx + 1; + auto nDigits = mInputROFs[timeBin.mLastIdx].getLastIdx() - mInputROFs[timeBin.mFirstIdx].getFirstIdx() + 1; + std::cout << fmt::format(" bin {}: firstIdx={} size={} ndigits={}/{}", j, timeBin.mFirstIdx, size, timeBin.mNDigitsPS, nDigits) << std::endl; + } + } + + if (rofFirstIdx < 0) { + // the range is empty, no ROF to store + mLastSavedTimeBin = lastBin; + return; + } + + // get the indexes of the first and last ROFs in this time cluster, and the corresponding interaction records + auto& firstRofInCluster = mInputROFs[rofFirstIdx]; + auto irFirst = firstRofInCluster.getBCData(); + auto& lastRofInCluster = mInputROFs[rofLastIdx]; + auto& irLast = lastRofInCluster.getBCData(); + + // a new time ROF is stored only if there are some digit ROFs in the interval + if (rofFirstIdx >= 0) { + // get the indexes of the first and last ROFs in this time cluster, and the corresponding interaction records + auto& firstRofInCluster = mInputROFs[rofFirstIdx]; + auto irFirst = firstRofInCluster.getBCData(); + auto& lastRofInCluster = mInputROFs[rofLastIdx]; + auto& irLast = lastRofInCluster.getBCData(); + + // get the index of the first digit and the number of digits in this time cluster + auto firstDigitIdx = firstRofInCluster.getFirstIdx(); + auto nDigits = lastRofInCluster.getLastIdx() - firstDigitIdx + 1; + + // compute the width in BC units of the time-cluster ROF + auto bcDiff = irLast.differenceInBC(irFirst); + auto bcWidth = bcDiff + lastRofInCluster.getBCWidth(); + + if (irFirst.orbit == targetOrbit) { + int peakBC = (peak < 0) ? -1 : mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().bc; + std::cout << "[TOTO] storing ROF from " << irFirst.orbit << " / " << irFirst.bc + << " to " << irLast.orbit << " / " << irLast.bc << std::endl; + std::cout << "[TOTO] mLastSavedInputRof " << mLastSavedInputRof << " mRofHasPeak.back() " << mRofHasPeak.back() + << " peak " << peak << " peakBC " << peakBC << std::endl; + } + + bool merged = false; + // check the gap with respect to the previously saved ROF + // if smaller than two ADC clock cycles, the current range is attached to the previous ROF + bool doMerge = false; + if (irFirst.orbit == targetOrbit) { + auto& prevRof = mOutputROFs.back(); + std::cout << "[TOTO] checking " << irFirst.orbit << " / " << irFirst.bc << " and " + << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc + << " peak " << peak << " firstBin " << firstBin << " mLastPeakEnd " << mLastPeakEnd + << std::endl; + } + if (mMergeROFs && hasPeak && peak >= firstBin && mLastPeakEnd >= 0) { + int32_t peakGap = (peak - mLastPeakEnd) * mBinWidth; + + auto& prevRof = mOutputROFs.back(); + if (irFirst.orbit == targetOrbit) { + std::cout << "[TOTO] checking " << irFirst.orbit << " / " << irFirst.bc << " and " + << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc + << " peak " << peak << " mLastPeakEnd " << mLastPeakEnd << " peakGap " << peakGap + << std::endl; + } + + // merge only if the gap is smaller than threshold + if (peakGap <= 12) { + doMerge = true; + } + } + if (doMerge) { + auto& lastRofPrevCluster = mInputROFs[mLastSavedInputRof]; + auto& irPrev = lastRofPrevCluster.getBCData(); + auto bcGap = irFirst.differenceInBC(irPrev); + + if (true) { + auto& prevRof = mOutputROFs.back(); + if (irFirst.orbit == targetOrbit) { + std::cout << "[TOTO] merging " << irFirst.orbit << " / " << irFirst.bc << " into " + << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc + << " bcWidth " << bcWidth << " " << prevRof.getBCWidth() + << " nDigits " << nDigits << " " << prevRof.getNEntries() + << std::endl; + } + bcDiff = irLast.differenceInBC(prevRof.getBCData()); + bcWidth = bcDiff + lastRofInCluster.getBCWidth(); + nDigits += prevRof.getNEntries(); + prevRof = ROFRecord(prevRof.getBCData(), prevRof.getFirstIdx(), nDigits, bcWidth); + merged = true; + if (irFirst.orbit == targetOrbit) { + std::cout << "[TOTO] " << irFirst.orbit << " / " << irFirst.bc << " merged into " + << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << std::endl; + } + } + } + + if (!merged) { + // create a ROF that includes all the digits in this time cluster + mOutputROFs.emplace_back(irFirst, firstDigitIdx, nDigits, bcWidth); + mRofHasPeak.emplace_back(hasPeak); + if (irFirst.orbit == targetOrbit) { + std::cout << "[TOTO] added ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc + << " bcWidth " << mOutputROFs.back().getBCWidth() + << " nDigits " << mOutputROFs.back().getNEntries() + << std::endl; + } + } + + if (mDebug) { + std::cout << fmt::format("new ROF stored: firstDigitIdx={} size={} bcWidth={}", firstDigitIdx, nDigits, bcWidth) << std::endl; + } + mLastSavedInputRof = rofLastIdx; + } + + mLastSavedTimeBin = lastBin; +} + +//_________________________________________________________________________________________________ + +void ROFTimeClusterFinderV2::process() +{ + if (mDebug) { + std::cout << "\n\n==================\n[ROFTimeClusterFinderV2] processing new TF\n" + << std::endl; + } + + int clusterSizeBins = mTimeClusterSize / mBinWidth; + + initTimeBins(); + mOutputROFs.clear(); + mRofHasPeak.clear(); + + mLastSavedTimeBin = -1; + mLastSavedInputRof = -1; + mLastPeak = -1; + mLastPeakEnd = -1; + int32_t peak{-1}; + while ((peak = getNextPeak()) >= 0) { + int32_t peakStart = peak - clusterSizeBins / 2; + if (peakStart <= mLastSavedTimeBin) { + peakStart = mLastSavedTimeBin + 1; + } + int32_t peakEnd = getPeakEnd(peak, peakStart + clusterSizeBins - 1); + if (peakEnd < peakStart) { + continue; + } + int peakOrbit = mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().orbit; + + // peak found, we add the corresponding rof(s) + // first we fill the gap between the last peak and the current one, if needed + if (mDebug) { + std::cout << fmt::format("peakStart={} mLastSavedTimeBin={}", peakStart, mLastSavedTimeBin) << std::endl; + } + while ((peakStart - mLastSavedTimeBin) > 1) { + int32_t firstBin = mLastSavedTimeBin + 1; + int32_t lastBin = firstBin + clusterSizeBins - 1; + if (lastBin >= peakStart) { + lastBin = peakStart - 1; + } + storeROF(firstBin, lastBin, -1); + if (mOutputROFs.back().getBCData().orbit == targetOrbit) { + std::cout << "[TOTO] last stored ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc + << " bcWidth " << mOutputROFs.back().getBCWidth() + << " nDigits " << mOutputROFs.back().getNEntries() + << std::endl; + } + } + bool print = false; + if (peakOrbit == targetOrbit) { + print = true; + int peakBc = mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().bc; + //std::cout << "[TOTO] peak bc " << peakBc << std::endl; + } + storeROF(peakStart, peakEnd, peak); + mLastPeak = peak; + mLastPeakEnd = peakEnd; + if (mOutputROFs.back().getBCData().orbit == targetOrbit) { + std::cout << "[TOTO] last stored ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc + << " bcWidth " << mOutputROFs.back().getBCWidth() + << " nDigits " << mOutputROFs.back().getNEntries() + << std::endl; + } + } + + // fill the gap between the last peak and the end of the TF + int32_t firstBin = mLastSavedTimeBin + 1; + while (firstBin < mNbinsInOneTF) { + int32_t lastBin = firstBin + clusterSizeBins - 1; + storeROF(firstBin, lastBin, -1); + + firstBin = lastBin + 1; + } + + for (auto& rof : mOutputROFs) { + if (rof.getBCData().orbit == targetOrbit) { + std::cout << "[TOTO] ROF " << rof.getBCData().orbit << " / " << rof.getBCData().bc + << " bcWidth " << rof.getBCWidth() + << " nDigits " << rof.getNEntries() + << std::endl; + } + + } +} + +//_________________________________________________________________________________________________ + +char* ROFTimeClusterFinderV2::saveROFRsToBuffer(size_t& bufSize) +{ + static constexpr size_t sizeOfROFRecord = sizeof(o2::mch::ROFRecord); + + if (mDebug) { + dumpOutputROFs(); + } + + bufSize = mOutputROFs.size() * sizeOfROFRecord; + o2::mch::ROFRecord* buf = reinterpret_cast(malloc(bufSize)); + if (!buf) { + bufSize = 0; + return nullptr; + } + + o2::mch::ROFRecord* p = buf; + for (size_t i = 0; i < mOutputROFs.size(); i++) { + auto& rof = mOutputROFs[i]; + memcpy(p, &(rof), sizeOfROFRecord); + p += 1; + } + + return reinterpret_cast(buf); +} + +//_________________________________________________________________________________________________ + +void ROFTimeClusterFinderV2::dumpInputROFs() +{ + for (size_t i = 0; i < mInputROFs.size(); i++) { + auto& rof = mInputROFs[i]; + const auto ir = rof.getBCData(); + std::cout << fmt::format(" ROF {} RANGE {}-{} SIZE {} IR {}-{},{}\n", i, rof.getFirstIdx(), rof.getLastIdx(), + rof.getNEntries(), ir.orbit, ir.bc, ir.toLong()); + } +} + +//_________________________________________________________________________________________________ + +void ROFTimeClusterFinderV2::dumpOutputROFs() +{ + for (size_t i = 0; i < mOutputROFs.size(); i++) { + auto& rof = mOutputROFs[i]; + std::cout << fmt::format(" ROF {} RANGE {}-{} SIZE {} IR {}-{},{}\n", i, rof.getFirstIdx(), rof.getLastIdx(), + rof.getNEntries(), rof.getBCData().orbit, rof.getBCData().bc, rof.getBCData().toLong()); + } +} + +} // namespace mch +} // namespace o2 diff --git a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx new file mode 100644 index 0000000000000..ef1ae201a3012 --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx @@ -0,0 +1,346 @@ +// 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 TimeClusterFinderSpecV2.cxx +/// \brief Implementation of a data processor to run the time clusterizer +/// +/// \author Andrea Ferrero, CEA + +#include "MCHTimeClustering/TimeClusterFinderSpecV2.h" + +#include +#include +#include +#include +#include +#include + +#include + +#include "CommonDataFormat/IRFrame.h" + +#include "Framework/CallbackService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Logger.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" + +#include "MCHBase/TrackerParam.h" +#include "MCHDigitFiltering/DigitFilter.h" +#include "MCHROFFiltering/IRFrameFilter.h" +#include "MCHROFFiltering/MultiplicityFilter.h" +#include "MCHROFFiltering/TrackableFilter.h" +#include "MCHTimeClustering/ROFTimeClusterFinderV2.h" +#include "MCHTimeClustering/TimeClusterizerParamV2.h" + +namespace o2 +{ +namespace mch +{ + +using namespace std; +using namespace o2::framework; + +class TimeClusterFinderTaskV2 +{ + public: + //_________________________________________________________________________________________________ + void init(framework::InitContext& ic) + { + const auto& param = TimeClusterizerParamV2::Instance(); + mTimeClusterWidth = param.maxClusterWidth; + mPeakSearchWindow = param.peakSearchWindow; + mPeakSearchNbins = param.peakSearchNbins; + mPeakSearchNDigitsMin = param.peakSearchNDigitsMin; + mMinDigitPerROF = param.minDigitsPerROF; + mOnlyTrackable = param.onlyTrackable; + mPeakSearchSignalOnly = param.peakSearchSignalOnly; + mMergeROFs = param.mergeROFs; + mIRFramesOnly = param.irFramesOnly; + mDebug = ic.options().get("mch-debug"); + mROFRejectionFraction = param.rofRejectionFraction; + + if (mDebug) { + fair::Logger::SetConsoleColor(true); + } + // number of bins must be >= 3 + if (mPeakSearchNbins < 3) { + mPeakSearchNbins = 3; + } + // number of bins must be odd + if ((mPeakSearchNbins % 2) == 0) { + mPeakSearchNbins += 1; + } + if (mROFRejectionFraction > 0) { + std::random_device rd; + mGenerator = std::mt19937(rd()); + } + + LOGP(info, "TimeClusterWidth : {}", mTimeClusterWidth); + LOGP(info, "PeakSearchNbins : {}", mPeakSearchNbins); + LOGP(info, "MinDigitPerROF : {}", mMinDigitPerROF); + LOGP(info, "OnlyTrackable : {}", mOnlyTrackable); + LOGP(info, "PeakSearchSignalOnly : {}", mPeakSearchSignalOnly); + LOGP(info, "IRFramesOnly : {}", mIRFramesOnly); + LOGP(info, "ROFRejectionFraction : {}", mROFRejectionFraction); + + auto stop = [this]() { + if (mTFcount) { + LOGP(info, "duration = {} us / TF", mTimeProcess.count() * 1000 / mTFcount); + } + }; + ic.services().get().set(stop); + } + + ROFFilter createRandomRejectionFilter(float rejectionFraction) + { + return [this, rejectionFraction](const ROFRecord& /*rof*/) { + double rnd = mDistribution(mGenerator); + return rnd > rejectionFraction; + }; + } + + //_________________________________________________________________________________________________ + void run(framework::ProcessingContext& pc) + { + auto rofs = pc.inputs().get>("rofs"); + auto digits = pc.inputs().get>("digits"); + + o2::mch::ROFTimeClusterFinderV2 rofProcessor(rofs, + digits, + mTimeClusterWidth, + mPeakSearchWindow, + mPeakSearchNbins, + mPeakSearchNDigitsMin, + mPeakSearchSignalOnly, + mMergeROFs, + mDebug); + + if (mDebug) { + LOGP(warning, "{:=>60} ", fmt::format("{:6d} Input ROFS", rofs.size())); + } + + auto tStart = std::chrono::high_resolution_clock::now(); + rofProcessor.process(); + auto tEnd = std::chrono::high_resolution_clock::now(); + mTimeProcess += tEnd - tStart; + + if (mDebug) { + LOGP(warning, "{:=>60} ", fmt::format("{:6d} Output ROFS", rofProcessor.getROFRecords().size())); + } + + auto& outRofs = pc.outputs().make>(OutputRef{"rofs"}); + const auto& pRofs = rofProcessor.getROFRecords(); + + // prepare the list of filters we want to apply to ROFs + std::vector filters; + + if (mOnlyTrackable) { + // selects only ROFs that are trackable + const auto& trackerParam = TrackerParam::Instance(); + std::array requestStation{ + trackerParam.requestStation[0], + trackerParam.requestStation[1], + trackerParam.requestStation[2], + trackerParam.requestStation[3], + trackerParam.requestStation[4]}; + filters.emplace_back(createTrackableFilter(digits, + requestStation, + trackerParam.moreCandidates)); + } + if (mMinDigitPerROF > 0) { + // selects only those ROFs have that minimum number of digits + filters.emplace_back(createMultiplicityFilter(mMinDigitPerROF)); + } + if (mIRFramesOnly) { + // selects only those ROFs that overlop some IRFrame + auto irFrames = pc.inputs().get>("irframes"); + filters.emplace_back(createIRFrameFilter(irFrames)); + } + + std::string extraMsg = ""; + + if (mROFRejectionFraction > 0) { + filters.emplace_back(createRandomRejectionFilter(mROFRejectionFraction)); + extraMsg = fmt::format(" (CAUTION : hard-rejected {:3.0f}% of the output ROFs)", mROFRejectionFraction * 100); + } + + // a single filter which is the AND combination of the elements of the filters vector + auto filter = createROFFilter(filters); + + std::copy_if(begin(pRofs), + end(pRofs), + std::back_inserter(outRofs), + filter); + + const float p1 = rofs.size() > 0 ? 100. * pRofs.size() / rofs.size() : 0; + const float p2 = rofs.size() > 0 ? 100. * outRofs.size() / rofs.size() : 0; + + LOGP(info, + "NEW TF {} Processed {} input ROFs, " + "time-clusterized them into {} ROFs ({:3.0f}%) " + "and output {} ({:3.0f}%) of them{}", + mTFcount, rofs.size(), + pRofs.size(), p1, + outRofs.size(), p2, extraMsg); + mTFcount += 1; + + std::vector> orbits = { + //{376285608, 1455}, + //{376285618, 1075}, + //{376293180, 2152}, + //{376353595, 1807}, + //{376376230, 1441}, + //{376285604, 828}, + //{376429114, 2550}, + //{376670762, 2548}, + //{376678313, 453}, + //{376867115, 2335}, + //{378166083, 1029} + //{376270496, 2752}, + //{376300693, 3249}, + //{376376238, 352}, + //{376376238, 1843}, + //{376421564, 1702} + //{376361141, 2781}, + //{376784078, 3156} + }; + //int orbit = 376293180; //376285618; //376285608; + for (auto o : orbits) { + int orbit = o.first; + int bc = o.second; + std::ofstream rofsOutput; + int id = 0; + for (auto& rof : pRofs) { + if (rof.getBCData().orbit != orbit) { + continue; + } + if (rof.getBCData().bc < (bc - 500)) { + continue; + } + if (rof.getBCData().bc > (bc + 500)) { + continue; + } + if (!rofsOutput.is_open()) { + rofsOutput.open(fmt::format("rofs-{}.txt", orbit)); + } + for (auto& irof : rofs) { + if (irof.getBCData().orbit != orbit) { + continue; + } + if (irof.getBCData().bc < rof.getBCData().bc) { + continue; + } + if (irof.getBCData().bc >= (rof.getBCData().bc + rof.getBCWidth())) { + break; + } + rofsOutput << fmt::format("{:>3}", id) << " " << fmt::format("{:>3}", rof.getBCWidth()) << (filter(rof) ? " + " : " - ") + << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4} ", irof.getNEntries()); + for (int i = 0; i < irof.getNEntries(); i++) { + rofsOutput << "*"; + if (i > 100) break; + } + rofsOutput << std::endl; + } + rofsOutput << "------" << std::endl; + id += 1; + } + rofsOutput.close(); + } + for (auto& rof : pRofs) { + std::ofstream rofsOutput; + if (rof.getBCWidth() < 200) { + continue; + } + int orbit = rof.getBCData().orbit; + if (!rofsOutput.is_open()) { + rofsOutput.open(fmt::format("rofs-{}-{}.txt", orbit, rof.getBCData().bc)); + } + for (auto& irof : rofs) { + if (irof.getBCData().orbit != orbit) { + continue; + } + if (irof.getBCData().bc < rof.getBCData().bc) { + continue; + } + if (irof.getBCData().bc >= (rof.getBCData().bc + rof.getBCWidth())) { + break; + } + rofsOutput << fmt::format("{:>3}", rof.getBCWidth()) << (filter(rof) ? " + " : " - ") + << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4} ", irof.getNEntries()); + for (int i = 0; i < irof.getNEntries(); i++) { + rofsOutput << "*"; + if (i > 100) break; + } + rofsOutput << std::endl; + } + rofsOutput.close(); + } + } + + private: + std::chrono::duration + mTimeProcess{}; ///< timer + + uint32_t mTimeClusterWidth; ///< maximum size of one time cluster, in bunch crossings + uint32_t mPeakSearchWindow; ///< width of the peak search window, in BC units + uint32_t mPeakSearchNbins; ///< number of time bins for the peak search algorithm (must be an odd number >= 3) + uint32_t mPeakSearchNDigitsMin; ///< minimum number of digits for peak candidates + + int mTFcount{0}; ///< number of processed time frames + int mDebug{0}; ///< verbosity flag + int mMinDigitPerROF; ///< minimum digit per ROF threshold + bool mPeakSearchSignalOnly; ///< only use signal-like hits in peak search + bool mMergeROFs; ///< whether to merge consecutive ROFs + bool mOnlyTrackable; ///< only keep ROFs that are trackable + bool mIRFramesOnly; ///< only keep ROFs that overlap some IRFrame + float mROFRejectionFraction; ///< fraction of output ROFs to reject (to save time in sync reco). MUST BE 0 for anything but Pt2 reco ! + std::uniform_real_distribution mDistribution{0.0, 1.0}; + std::mt19937 mGenerator; +}; + +//_________________________________________________________________________________________________ +o2::framework::DataProcessorSpec + getTimeClusterFinderSpecV2(const char* specName, + std::string_view inputDigitDataDescription, + std::string_view inputDigitRofDataDescription, + std::string_view outputDigitRofDataDescription, + std::string_view inputIRFrameDataDescription) +{ + std::string input = fmt::format("rofs:MCH/{}/0;digits:MCH/{}/0", + inputDigitRofDataDescription.data(), + inputDigitDataDescription.data()); + if (TimeClusterizerParamV2::Instance().irFramesOnly && inputIRFrameDataDescription.size()) { + LOGP(info, "will select IRFrames from {}", inputIRFrameDataDescription); + input += ";irframes:"; + input += inputIRFrameDataDescription; + } + std::string output = fmt::format("rofs:MCH/{}/0", outputDigitRofDataDescription.data()); + + std::vector outputs; + auto matchers = select(output.c_str()); + for (auto& matcher : matchers) { + outputs.emplace_back(DataSpecUtils::asOutputSpec(matcher)); + } + + return DataProcessorSpec{ + specName, + Inputs{select(input.c_str())}, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{{"mch-debug", VariantType::Bool, false, {"enable verbose output"}}}}; +} +} // end namespace mch +} // end namespace o2 diff --git a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterizerParamV2.cxx b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterizerParamV2.cxx new file mode 100644 index 0000000000000..4b0df0f952edf --- /dev/null +++ b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterizerParamV2.cxx @@ -0,0 +1,15 @@ +// 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 "MCHTimeClustering/TimeClusterizerParamV2.h" +#include "CommonUtils/ConfigurableParam.h" + +O2ParamImpl(o2::mch::TimeClusterizerParamV2) diff --git a/Detectors/MUON/MCH/TimeClustering/src/digits-to-timeclusters-workflow.cxx b/Detectors/MUON/MCH/TimeClustering/src/digits-to-timeclusters-workflow.cxx index 9c049a30c69e8..0665b5f6b4044 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/digits-to-timeclusters-workflow.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/digits-to-timeclusters-workflow.cxx @@ -21,6 +21,7 @@ #include "Framework/ControlService.h" #include "Framework/Task.h" #include "MCHTimeClustering/TimeClusterFinderSpec.h" +#include "MCHTimeClustering/TimeClusterFinderSpecV2.h" using namespace o2; using namespace o2::framework; @@ -30,6 +31,7 @@ void customize(std::vector& workflowOptions) workflowOptions.push_back(ConfigParamSpec{"input-digits-data-description", VariantType::String, "F-DIGITS", {"description string for the input digit data"}}); workflowOptions.push_back(ConfigParamSpec{"input-digitrofs-data-description", VariantType::String, "F-DIGITROFS", {"description string for the input digit rofs data"}}); workflowOptions.push_back(ConfigParamSpec{"output-digitrofs-data-description", VariantType::String, "TC-F-DIGITROFS", {"description string for the output digit rofs data"}}); + workflowOptions.push_back(ConfigParamSpec{"new-rof-clustering", VariantType::Bool, false, {"enable new version of rof clustering"}}); workflowOptions.push_back(ConfigParamSpec{ "configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}); } @@ -40,6 +42,15 @@ WorkflowSpec defineDataProcessing(const ConfigContext& cc) { o2::conf::ConfigurableParam::updateFromString(cc.options().get("configKeyValues")); + if (cc.options().get("new-rof-clustering")) { + return { + o2::mch::getTimeClusterFinderSpecV2( + "mch-time-clustering", + cc.options().get("input-digits-data-description"), + cc.options().get("input-digitrofs-data-description"), + cc.options().get("output-digitrofs-data-description"))}; + } + return { o2::mch::getTimeClusterFinderSpec( "mch-time-clustering", diff --git a/Detectors/MUON/MCH/Workflow/src/reco-workflow.cxx b/Detectors/MUON/MCH/Workflow/src/reco-workflow.cxx index 29e34f254bf28..e5cb71f79d047 100644 --- a/Detectors/MUON/MCH/Workflow/src/reco-workflow.cxx +++ b/Detectors/MUON/MCH/Workflow/src/reco-workflow.cxx @@ -30,6 +30,7 @@ #include "MCHStatus/StatusMapCreatorParam.h" #include "MCHStatus/StatusMapCreatorSpec.h" #include "MCHTimeClustering/TimeClusterFinderSpec.h" +#include "MCHTimeClustering/TimeClusterFinderSpecV2.h" #include "MCHTracking/TrackFinderSpec.h" #include "MCHWorkflow/ClusterFinderOriginalSpec.h" #include "MCHWorkflow/ErrorWriterSpec.h" @@ -59,6 +60,7 @@ void customize(std::vector& workflowOptions) {"disable-tracking", o2::framework::VariantType::Bool, false, {"disable tracking step (for debug)"}}, {"digits", VariantType::Bool, false, {"Write digits associated to tracks"}}, {"triggered", VariantType::Bool, false, {"use MID to trigger the MCH reconstruction"}}, + {"new-rof-clustering", VariantType::Bool, false, {"enable new version of rof clustering"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); std::swap(workflowOptions, options); @@ -74,6 +76,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) auto disableRootInput = configcontext.options().get("disable-root-input"); auto digits = configcontext.options().get("digits"); auto triggered = configcontext.options().get("triggered"); + auto newRofClustering = configcontext.options().get("new-rof-clustering"); auto useMC = !configcontext.options().get("disable-mc"); auto disableClustering = configcontext.options().get("disable-clustering"); auto disableTracking = disableClustering || configcontext.options().get("disable-tracking"); @@ -100,10 +103,19 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) "F-DIGITROFS", "E-F-DIGITROFS", "F-DIGITLABELS", "E-F-DIGITLABELS")); } else { + if (newRofClustering) { + std::cout << "[TOTO] calling getTimeClusterFinderSpecV2()" << std::endl; + specs.emplace_back(o2::mch::getTimeClusterFinderSpecV2("mch-time-cluster-finder-new", + "F-DIGITS", + "F-DIGITROFS", + "TC-F-DIGITROFS")); + } else { + std::cout << "[TOTO] calling getTimeClusterFinderSpec()" << std::endl; specs.emplace_back(o2::mch::getTimeClusterFinderSpec("mch-time-cluster-finder", "F-DIGITS", "F-DIGITROFS", "TC-F-DIGITROFS")); + } } specs.emplace_back(o2::mch::getPreClusterFinderSpec("mch-precluster-finder", From 6f5b8a9ae74e0032900493388d24cf102fdb994d Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Fri, 29 Dec 2023 17:16:58 +0100 Subject: [PATCH 2/2] Temporary commit --- .../DigitFiltering/src/DigitFilteringSpec.cxx | 4 +- .../src/PreClusterFinderSpec.cxx | 26 ++- .../ROFTimeClusterFinderV2.h | 2 + .../src/ROFTimeClusterFinder.cxx | 8 +- .../src/ROFTimeClusterFinderV2.cxx | 141 ++++++++++++-- .../src/TimeClusterFinderSpec.cxx | 4 + .../src/TimeClusterFinderSpecV2.cxx | 176 ++++++++++++++++-- .../src/ClusterFinderOriginalSpec.cxx | 4 +- .../MUON/MCH/Workflow/src/DataDecoderSpec.cxx | 1 + 9 files changed, 319 insertions(+), 47 deletions(-) diff --git a/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx b/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx index ea173bd1fc7cc..ccac657e82ea0 100644 --- a/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx +++ b/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx @@ -143,10 +143,10 @@ class DigitFilteringTask auto labelMsg = mUseMC ? fmt::format("| {} labels (out of {})", oLabels->getNElements(), iLabels->getNElements()) : ""; - LOGP(info, "Kept after filtering : {} rofs (out of {}) | {} digits (out of {}) {}", + /*LOGP(info, "Kept after filtering : {} rofs (out of {}) | {} digits (out of {}) {}", oRofs.size(), iRofs.size(), oDigits.size(), iDigits.size(), - labelMsg); + labelMsg);*/ if (mTimeCalib != 0) { shiftDigitsTime(oRofs, oDigits); diff --git a/Detectors/MUON/MCH/PreClustering/src/PreClusterFinderSpec.cxx b/Detectors/MUON/MCH/PreClustering/src/PreClusterFinderSpec.cxx index 2989401b55628..9c05fe5e06d82 100644 --- a/Detectors/MUON/MCH/PreClustering/src/PreClusterFinderSpec.cxx +++ b/Detectors/MUON/MCH/PreClustering/src/PreClusterFinderSpec.cxx @@ -40,6 +40,7 @@ #include "MCHBase/PreCluster.h" #include "MCHBase/SanityCheck.h" #include "MCHPreClustering/PreClusterFinder.h" +#include "MCHMappingInterface/Segmentation.h" #include #include @@ -210,11 +211,34 @@ class PreClusterFinderTask }); mErrorMap.add(errorMap); - LOGP(info, "Processed {} digit rofs with {} digits and output {} precluster rofs with {} preclusters and {} digits", + /*LOGP(info, "Processed {} digit rofs with {} digits and output {} precluster rofs with {} preclusters and {} digits", digitROFs.size(), nDigitsInRofs, preClusterROFs.size(), mPreClusters.size(), mUsedDigits.size()); + for (auto& preCluster: mPreClusters) { + double totADC = 0; + for (int i = preCluster.firstDigit; i <= preCluster.lastDigit(); i++) { + auto& d = mUsedDigits[i]; + totADC += d.getADC(); + } + if (totADC < 7000) continue; + for (int i = preCluster.firstDigit; i <= preCluster.lastDigit(); i++) { + auto& d = mUsedDigits[i]; + auto detID = d.getDetID(); + auto padID = d.getPadID(); + if (padID < 0) { + continue; + } + const o2::mch::mapping::Segmentation& segment = o2::mch::mapping::segmentation(detID); + bool bend = segment.isBendingPad(padID); + float X = segment.padPositionX(padID); + float Y = segment.padPositionY(padID); + uint32_t time = d.getTime(); + LOGP(info, " [TOTO] DE {:4d} PAD {:5d} ADC {:6d} TIME {:4d}\tC {} PAD_XY {:+2.2f} , {:+2.2f}", + detID, padID, d.getADC(), time, (bend ? (int)0 : (int)1), X, Y); + } + }*/ } private: diff --git a/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h index 9cf05e6b32cc9..313c337fc53bc 100644 --- a/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h +++ b/Detectors/MUON/MCH/TimeClustering/include/MCHTimeClustering/ROFTimeClusterFinderV2.h @@ -66,6 +66,8 @@ class ROFTimeClusterFinderV2 int32_t mFirstIdx{-1}; ///< index of the first digit ROF in the bin int32_t mLastIdx{-1}; ///< index of the last digit ROF in the bin uint32_t mNDigitsPS{0}; ///< number of digits in the bin for the peak search step + std::array mHasChamber; + std::array mHasStation; TimeBin() = default; diff --git a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx index 7187aaea19243..fae54fb0eab19 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinder.cxx @@ -128,12 +128,12 @@ int32_t ROFTimeClusterFinder::getNextPeak() int32_t sPadding = mNbinsInOneWindow / 2; if (mDebug) { - std::cout << "Searching peak from " << mLastSavedTimeBin + 1 << std::endl; + std::cout << "Searching peak from " << mLastSavedTimeBin + 1 << " mNbinsInOneTF " << mNbinsInOneTF << std::endl; } // loop over the bins and search for local maxima // a local maxima is defined as a bin tht is higher than all the surrounding 4 bins (2 below and 2 above) - for (int32_t i = mLastSavedTimeBin + sPadding + 1; i < mNbinsInOneTF; i++) { + for (int32_t i = mLastSavedTimeBin + 1; i < mNbinsInOneTF; i++) { auto& peak = mTimeBins[i]; if (peak.empty()) { continue; @@ -235,10 +235,10 @@ void ROFTimeClusterFinder::storeROF(int32_t firstBin, int32_t lastBin) void ROFTimeClusterFinder::process() { - //if (mDebug) { + if (mDebug) { std::cout << "\n\n==================\n[ROFTimeClusterFinder] processing new TF\n" << std::endl; - //} + } initTimeBins(); mOutputROFs.clear(); diff --git a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx index 0dc3e99e050e3..2c759283e09a4 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/ROFTimeClusterFinderV2.cxx @@ -102,6 +102,13 @@ void ROFTimeClusterFinderV2::initTimeBins() } auto& timeBin = mTimeBins[binIdx]; + timeBin.mHasChamber = { + false, false, false, false, false, + false, false, false, false, false + }; + timeBin.mHasStation = { + false, false, false, false, false + }; if (timeBin.mFirstIdx < 0) { timeBin.mFirstIdx = iRof; @@ -112,9 +119,20 @@ void ROFTimeClusterFinderV2::initTimeBins() if (mImprovePeakSearch) { auto rofDigits = mDigits.subspan(rof.getFirstIdx(), rof.getNEntries()); for (auto& digit : rofDigits) { - if (mIsGoodDigit(digit)) { - nDigitsPS += 1; + if (!mIsGoodDigit(digit)) { + continue; } + + nDigitsPS += 1; + int deId = digit.getDetID(); + int chId = (deId / 100) - 1; + if (chId < 0) continue; + if (chId >= 10) continue; + timeBin.mHasChamber[chId] = true; + int stId = chId / 2; + if (stId < 0) continue; + if (stId >= 5) continue; + timeBin.mHasStation[stId] = true; } } else { nDigitsPS = rof.getNEntries(); @@ -135,7 +153,7 @@ void ROFTimeClusterFinderV2::initTimeBins() //_________________________________________________________________________________________________ -static int targetOrbit = 376784078; +static int targetOrbit = 74410449; int32_t ROFTimeClusterFinderV2::getNextPeak() { @@ -155,45 +173,122 @@ int32_t ROFTimeClusterFinderV2::getNextPeak() continue; } + // check number of stations associated to the digits in the search window + std::array hasChamber = mTimeBins[i].mHasChamber; + std::array hasStation = mTimeBins[i].mHasStation; + int peakOrbit = mInputROFs[peak.mFirstIdx].getBCData().orbit; if (peakOrbit == targetOrbit) { int peakBc = mInputROFs[peak.mFirstIdx].getBCData().bc; - std::cout << "[TOTO] peak seed " << i << " " << peak.mNDigitsPS << " bc " << peakBc << " padding " << sPadding << std::endl; + std::cout << "[TOTO] ==== " << std::endl; + std::cout << "[TOTO] checking peak seed " << i << " " << peak.mNDigitsPS << " bc " << peakBc << " stations [ "; + for (auto st : mTimeBins[i].mHasStation) { + if (st) std::cout<<"1 "; + else std::cout<<"0 "; + } + std::cout << "]" << std::endl; } bool found{true}; // the peak must be strictly higher than previous bins for (int j = i - sPadding; j < i; j++) { + if (j < 0) { continue; } + //if (j <= mLastSavedTimeBin) { continue; } + if (j > mLastSavedTimeBin && peakOrbit == targetOrbit) { int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; - std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << std::endl; + std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << " stations [ "; + for (auto st : mTimeBins[j].mHasStation) { + if (st) std::cout<<"1 "; + else std::cout<<"0 "; + } + std::cout << "]" << std::endl; } if (j > mLastSavedTimeBin && peak <= mTimeBins[j]) { found = false; break; } + + for (size_t k = 0; k < hasChamber.size(); k++) { + if (mTimeBins[j].mHasChamber[k]) { + hasChamber[k] = true; + } + } + + for (size_t k = 0; k < hasStation.size(); k++) { + if (mTimeBins[j].mHasStation[k]) { + hasStation[k] = true; + if (peakOrbit == targetOrbit) { + int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; + int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; + //std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " adding station " << k << std::endl; + } + } + } } // the peak must be higher than or equal to next bins for (int j = i + 1; j <= i + sPadding; j++) { + if (j >= mNbinsInOneTF) { break; } + if (j < mNbinsInOneTF && peakOrbit == targetOrbit) { int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; - std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << std::endl; + std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " nDigits " << mTimeBins[j].mNDigitsPS << " stations [ "; + for (auto st : mTimeBins[j].mHasStation) { + if (st) std::cout<<"1 "; + else std::cout<<"0 "; + } + std::cout << "]" << std::endl; } if (j < mNbinsInOneTF && peak < mTimeBins[j]) { found = false; break; } + + for (size_t k = 0; k < hasChamber.size(); k++) { + if (mTimeBins[j].mHasChamber[k]) { + hasChamber[k] = true; + } + } + + for (size_t k = 0; k < hasStation.size(); k++) { + if (mTimeBins[j].mHasStation[k]) { + hasStation[k] = true; + if (peakOrbit == targetOrbit) { + int orbit = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().orbit; + int bc = (mTimeBins[j].mFirstIdx < 0) ? -1 : mInputROFs[mTimeBins[j].mFirstIdx].getBCData().bc; + //std::cout << "[TOTO] " << j << " orbit " << orbit << " bc " << bc << " adding station " << k << std::endl; + } + } + } + } + + int nStations = 0; + for (auto st : hasStation) { + if (st) { + nStations += 1; + } } if (peakOrbit == targetOrbit) { - std::cout << "[TOTO] peak found " << found << std::endl; + int peakOrbit = mInputROFs[peak.mFirstIdx].getBCData().orbit; + int peakBc = mInputROFs[peak.mFirstIdx].getBCData().bc; + std::cout << "[TOTO] peak seed " << i << " orbit " << peakOrbit << " bc " << peakBc << " nDigits " << mTimeBins[i].mNDigitsPS << " found=" << found << " nStations " << nStations << " [ "; + for (auto st : hasStation) { + if (st) std::cout<<"1 "; + else std::cout<<"0 "; + } + std::cout << "]" << std::endl; } - if (!found) { + if (!found || (nStations < 4)) { continue; } + if (peakOrbit == targetOrbit) { + auto nDigits = mInputROFs[peak.mLastIdx].getLastIdx() - mInputROFs[peak.mFirstIdx].getFirstIdx() + 1; + std::cout << fmt::format("[TOTO] new peak found at bin {}, entries = {}/{}", i, peak.mNDigitsPS, nDigits) << std::endl; + } if (mDebug) { auto nDigits = mInputROFs[peak.mLastIdx].getLastIdx() - mInputROFs[peak.mFirstIdx].getFirstIdx() + 1; @@ -297,7 +392,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t const auto& rof = mInputROFs[timeBin.mFirstIdx]; auto newWidth = rof.getBCData().differenceInBC(prevRof.getBCData()) + rof.getBCWidth(); - if (prevRof.getBCData().orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || prevRof.getBCData().orbit == targetOrbit) { std::cout << "[TOTO] trying to extend " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << " with rof " << rof.getBCData().orbit << " / " << rof.getBCData().bc << " newWidth " << newWidth << std::endl; } @@ -316,7 +411,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t } if (firstBinNew > firstBin) { prevRof = ROFRecord(prevRof.getBCData(), prevRof.getFirstIdx(), nDigits, bcWidth); - if (prevRof.getBCData().orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || prevRof.getBCData().orbit == targetOrbit) { std::cout << "[TOTO] ROF " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << " extended to " << prevRof.getBCWidth() << std::endl; } @@ -348,6 +443,9 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t if (rofFirstIdx < 0) { // the range is empty, no ROF to store mLastSavedTimeBin = lastBin; + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit)) { + std::cout << "[TOTO] empty ROF " << std::endl; + } return; } @@ -373,7 +471,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t auto bcDiff = irLast.differenceInBC(irFirst); auto bcWidth = bcDiff + lastRofInCluster.getBCWidth(); - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { int peakBC = (peak < 0) ? -1 : mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().bc; std::cout << "[TOTO] storing ROF from " << irFirst.orbit << " / " << irFirst.bc << " to " << irLast.orbit << " / " << irLast.bc << std::endl; @@ -385,7 +483,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t // check the gap with respect to the previously saved ROF // if smaller than two ADC clock cycles, the current range is attached to the previous ROF bool doMerge = false; - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { auto& prevRof = mOutputROFs.back(); std::cout << "[TOTO] checking " << irFirst.orbit << " / " << irFirst.bc << " and " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc @@ -396,7 +494,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t int32_t peakGap = (peak - mLastPeakEnd) * mBinWidth; auto& prevRof = mOutputROFs.back(); - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { std::cout << "[TOTO] checking " << irFirst.orbit << " / " << irFirst.bc << " and " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << " peak " << peak << " mLastPeakEnd " << mLastPeakEnd << " peakGap " << peakGap @@ -415,7 +513,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t if (true) { auto& prevRof = mOutputROFs.back(); - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { std::cout << "[TOTO] merging " << irFirst.orbit << " / " << irFirst.bc << " into " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << " bcWidth " << bcWidth << " " << prevRof.getBCWidth() @@ -427,7 +525,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t nDigits += prevRof.getNEntries(); prevRof = ROFRecord(prevRof.getBCData(), prevRof.getFirstIdx(), nDigits, bcWidth); merged = true; - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { std::cout << "[TOTO] " << irFirst.orbit << " / " << irFirst.bc << " merged into " << prevRof.getBCData().orbit << " / " << prevRof.getBCData().bc << std::endl; } @@ -438,7 +536,7 @@ void ROFTimeClusterFinderV2::storeROF(int32_t firstBin, int32_t lastBin, int32_t // create a ROF that includes all the digits in this time cluster mOutputROFs.emplace_back(irFirst, firstDigitIdx, nDigits, bcWidth); mRofHasPeak.emplace_back(hasPeak); - if (irFirst.orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || irFirst.orbit == targetOrbit) { std::cout << "[TOTO] added ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc << " bcWidth " << mOutputROFs.back().getBCWidth() << " nDigits " << mOutputROFs.back().getNEntries() @@ -497,8 +595,11 @@ void ROFTimeClusterFinderV2::process() if (lastBin >= peakStart) { lastBin = peakStart - 1; } + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || peakOrbit == targetOrbit) { + std::cout << "[TOTO] storing dummy ROF, start " << firstBin << " end " << lastBin << std::endl; + } storeROF(firstBin, lastBin, -1); - if (mOutputROFs.back().getBCData().orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || peakOrbit == targetOrbit) { std::cout << "[TOTO] last stored ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc << " bcWidth " << mOutputROFs.back().getBCWidth() << " nDigits " << mOutputROFs.back().getNEntries() @@ -506,15 +607,15 @@ void ROFTimeClusterFinderV2::process() } } bool print = false; - if (peakOrbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || peakOrbit == targetOrbit) { print = true; int peakBc = mInputROFs[mTimeBins[peak].mFirstIdx].getBCData().bc; - //std::cout << "[TOTO] peak bc " << peakBc << std::endl; + std::cout << "[TOTO] storing collision ROF, peak " << peak << " start " << peakStart << " end " << peakEnd << " bc " << peakBc << std::endl; } storeROF(peakStart, peakEnd, peak); mLastPeak = peak; mLastPeakEnd = peakEnd; - if (mOutputROFs.back().getBCData().orbit == targetOrbit) { + if ((mOutputROFs.size() > 0 && mOutputROFs.back().getBCData().orbit == targetOrbit) || peakOrbit == targetOrbit) { std::cout << "[TOTO] last stored ROF " << mOutputROFs.back().getBCData().orbit << " / " << mOutputROFs.back().getBCData().bc << " bcWidth " << mOutputROFs.back().getBCWidth() << " nDigits " << mOutputROFs.back().getNEntries() diff --git a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpec.cxx b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpec.cxx index 456098a4cff41..c2d0d9cdc6d1e 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpec.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpec.cxx @@ -115,6 +115,8 @@ class TimeClusterFinderTask auto rofs = pc.inputs().get>("rofs"); auto digits = pc.inputs().get>("digits"); + const auto& tinfo = pc.services().get(); + auto firstTForbit = tinfo.firstTForbit; o2::mch::ROFTimeClusterFinder rofProcessor(rofs, digits, mTimeClusterWidth, mNbinsInOneWindow, mPeakSearchSignalOnly, mDebug); if (mDebug) { @@ -177,6 +179,7 @@ class TimeClusterFinderTask const float p1 = rofs.size() > 0 ? 100. * pRofs.size() / rofs.size() : 0; const float p2 = rofs.size() > 0 ? 100. * outRofs.size() / rofs.size() : 0; + if (mDebug) { LOGP(info, "TF {} Processed {} input ROFs, " "time-clusterized them into {} ROFs ({:3.0f}%) " @@ -184,6 +187,7 @@ class TimeClusterFinderTask mTFcount, rofs.size(), pRofs.size(), p1, outRofs.size(), p2, extraMsg); + } mTFcount += 1; } diff --git a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx index ef1ae201a3012..4f1096b5881bf 100644 --- a/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx +++ b/Detectors/MUON/MCH/TimeClustering/src/TimeClusterFinderSpecV2.cxx @@ -45,6 +45,8 @@ #include "MCHTimeClustering/ROFTimeClusterFinderV2.h" #include "MCHTimeClustering/TimeClusterizerParamV2.h" +#include "MCHDigitFiltering/DigitFilter.h" + namespace o2 { namespace mch @@ -128,9 +130,17 @@ class TimeClusterFinderTaskV2 mMergeROFs, mDebug); + const auto& tinfo = pc.services().get(); + auto firstTForbit = tinfo.firstTForbit; if (mDebug) { LOGP(warning, "{:=>60} ", fmt::format("{:6d} Input ROFS", rofs.size())); } + if (!rofs.empty()) { + uint32_t firstOrbit = rofs.front().getBCData().orbit; + uint32_t lastOrbit = rofs.back().getBCData().orbit; + uint32_t delta = lastOrbit - firstOrbit; + std::cout << fmt::format("[TOTO] orbits {:10d} {:10d} {:10d} delta {:10d} ", firstTForbit, firstOrbit, lastOrbit, delta) << std::endl; + } auto tStart = std::chrono::high_resolution_clock::now(); rofProcessor.process(); @@ -216,8 +226,25 @@ class TimeClusterFinderTaskV2 //{376421564, 1702} //{376361141, 2781}, //{376784078, 3156} + + //{74406400, 475}, + //{74442508, 1541}, + //{74442523, 1608}, + //{74442565, 1802}, + //{74442586, 712}, + //{74442635, 2260}, + //{74478663, 36}, + //{74479236, 935}, + //{74514817, 1670}, + //{74514951, 1609} + //{74442860, 1610}, + //{74478601, 1169} + {74479380, 1122}, + {74480393, 1644}, + {74552579, 1807} }; //int orbit = 376293180; //376285618; //376285608; + auto isGoodDigit = createDigitFilter(20, true, true); for (auto o : orbits) { int orbit = o.first; int bc = o.second; @@ -227,16 +254,23 @@ class TimeClusterFinderTaskV2 if (rof.getBCData().orbit != orbit) { continue; } - if (rof.getBCData().bc < (bc - 500)) { - continue; - } - if (rof.getBCData().bc > (bc + 500)) { + if (rof.getBCData().orbit < firstTForbit) { continue; } + std::cout << "[TOTO] " << fmt::format("{:>3}", id) << " " + << fmt::format("{:>10} {:>10} {:>4} {:>3}", firstTForbit, rof.getBCData().orbit, rof.getBCData().bc, rof.getBCWidth()) + << (filter(rof) ? " + " : " - ") << std::endl; + //if (rof.getBCData().bc < (bc - 500)) { + // continue; + //} + //if (rof.getBCData().bc > (bc + 500)) { + // continue; + //} if (!rofsOutput.is_open()) { rofsOutput.open(fmt::format("rofs-{}.txt", orbit)); } - for (auto& irof : rofs) { + for (int i = 0; i < rofs.size(); i++) { + auto& irof = rofs[i]; if (irof.getBCData().orbit != orbit) { continue; } @@ -246,12 +280,61 @@ class TimeClusterFinderTaskV2 if (irof.getBCData().bc >= (rof.getBCData().bc + rof.getBCWidth())) { break; } + std::cout << "[TOTO] " << fmt::format("{:>10} {:>4} {:>4}", irof.getBCData().orbit, irof.getBCData().bc, irof.getNEntries()) << std::endl; + std::array hasChamber = { + false, false, false, false, false, + false, false, false, false, false + }; + std::array hasStation = { + false, false, false, false, false + }; + for (int j = i - 2; j <= i + 2; j++) { + if (j < 0) continue; + if (j >= rofs.size()) continue; + auto& irof2 = rofs[i]; + for (int d = irof2.getFirstIdx(); d <= irof2.getLastIdx(); d++) { + const auto& digit = digits[d]; + if (!isGoodDigit(digit)) continue; + int deId = digit.getDetID(); + int chId = (deId / 100) - 1; + if (chId < 0) continue; + if (chId >= 10) continue; + hasChamber[chId] = true; + int stId = chId / 2; + if (stId < 0) continue; + if (stId >= 5) continue; + hasStation[stId] = true; + } + } + int nCh = 0; + for (auto ch : hasChamber) { + if (ch) nCh += 1; + } + int nSt = 0; + for (auto st : hasStation) { + if (st) nSt += 1; + } + int nGood = 0; + for (int d = irof.getFirstIdx(); d <= irof .getLastIdx(); d++) { + const auto& digit = digits[d]; + if (!isGoodDigit(digit)) continue; + nGood += 1; + } rofsOutput << fmt::format("{:>3}", id) << " " << fmt::format("{:>3}", rof.getBCWidth()) << (filter(rof) ? " + " : " - ") - << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4} ", irof.getNEntries()); - for (int i = 0; i < irof.getNEntries(); i++) { + << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4}/{:>4} ", irof.getNEntries(), nGood); + for (int i = 0; i < nGood; i++) { rofsOutput << "*"; - if (i > 100) break; + if (i > 100) { + break; + } + } + for (int i = nGood; i < irof.getNEntries(); i++) { + rofsOutput << "-"; + if (i > 100) { + break; + } } + rofsOutput << "| " << nCh << " / " << nSt; rofsOutput << std::endl; } rofsOutput << "------" << std::endl; @@ -261,29 +344,86 @@ class TimeClusterFinderTaskV2 } for (auto& rof : pRofs) { std::ofstream rofsOutput; - if (rof.getBCWidth() < 200) { + if (rof.getBCWidth() < 100) { continue; } int orbit = rof.getBCData().orbit; if (!rofsOutput.is_open()) { rofsOutput.open(fmt::format("rofs-{}-{}.txt", orbit, rof.getBCData().bc)); + rofsOutput << firstTForbit << std::endl; } - for (auto& irof : rofs) { - if (irof.getBCData().orbit != orbit) { - continue; - } - if (irof.getBCData().bc < rof.getBCData().bc) { + for (int i = 0; i < rofs.size(); i++) { + auto& irof = rofs[i]; + auto dBC = irof.getBCData().differenceInBC(rof.getBCData()); + if (dBC < 0) { continue; } - if (irof.getBCData().bc >= (rof.getBCData().bc + rof.getBCWidth())) { + if (dBC >= rof.getBCWidth()) { break; } + std::array hasChamber = { + false, false, false, false, false, + false, false, false, false, false + }; + std::array hasStation = { + false, false, false, false, false + }; + for (int j = i - 2; j <= i + 2; j++) { + if (j < 0) continue; + if (j >= rofs.size()) continue; + auto& irof2 = rofs[j]; + for (int d = irof2.getFirstIdx(); d <= irof2.getLastIdx(); d++) { + const auto& digit = digits[d]; + int deId = digit.getDetID(); + int chId = (deId / 100) - 1; + if (chId < 0) continue; + if (chId >= 10) continue; + hasChamber[chId] = true; + int stId = chId / 2; + if (stId < 0) continue; + if (stId >= 5) continue; + hasStation[stId] = true; + } + } + int nCh = 0; + for (auto ch : hasChamber) { + if (ch) nCh += 1; + } + int nSt = 0; + for (auto st : hasStation) { + if (st) nSt += 1; + } + + int nGood = 0; + std::array hasStation2 = { + false, false, false, false, false + }; + for (int d = irof.getFirstIdx(); d <= irof .getLastIdx(); d++) { + const auto& digit = digits[d]; + int deId = digit.getDetID(); + int chId = (deId / 100) - 1; + if (chId < 0) continue; + if (chId >= 10) continue; + int stId = chId / 2; + hasStation2[stId] = true; + if (!isGoodDigit(digit)) continue; + nGood += 1; + } rofsOutput << fmt::format("{:>3}", rof.getBCWidth()) << (filter(rof) ? " + " : " - ") - << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4} ", irof.getNEntries()); - for (int i = 0; i < irof.getNEntries(); i++) { + << fmt::format("{:>10}", irof.getBCData().orbit) << " " << fmt::format("{:>4}", irof.getBCData().bc) << " " << fmt::format("{:>4}/{:>4} ", nGood, irof.getNEntries()); + for (int i = 0; i < nGood; i++) { rofsOutput << "*"; - if (i > 100) break; + if (i > 100) { + break; + } + } + for (int i = nGood; i < irof.getNEntries(); i++) { + rofsOutput << "-"; + if (i > 100) { + break; + } } + rofsOutput << "| " << nCh << " / " << nSt; rofsOutput << std::endl; } rofsOutput.close(); diff --git a/Detectors/MUON/MCH/Workflow/src/ClusterFinderOriginalSpec.cxx b/Detectors/MUON/MCH/Workflow/src/ClusterFinderOriginalSpec.cxx index 2007f0634960f..6917e8cd84474 100644 --- a/Detectors/MUON/MCH/Workflow/src/ClusterFinderOriginalSpec.cxx +++ b/Detectors/MUON/MCH/Workflow/src/ClusterFinderOriginalSpec.cxx @@ -137,8 +137,8 @@ class ClusterFinderOriginalTask }); mErrorMap.add(errorMap); - LOGP(info, "Found {:4d} clusters from {:4d} preclusters in {:2d} ROFs", - clusters.size(), preClusters.size(), preClusterROFs.size()); + //LOGP(info, "Found {:4d} clusters from {:4d} preclusters in {:2d} ROFs", + // clusters.size(), preClusters.size(), preClusterROFs.size()); } private: diff --git a/Detectors/MUON/MCH/Workflow/src/DataDecoderSpec.cxx b/Detectors/MUON/MCH/Workflow/src/DataDecoderSpec.cxx index 0f5745af1f8f1..08712886ce71b 100644 --- a/Detectors/MUON/MCH/Workflow/src/DataDecoderSpec.cxx +++ b/Detectors/MUON/MCH/Workflow/src/DataDecoderSpec.cxx @@ -254,6 +254,7 @@ class DataDecoderTask auto tStart = std::chrono::high_resolution_clock::now(); mDecoder->reset(); + mDecoder->setOrbitsInTF(256); for (auto&& input : pc.inputs()) { if (input.spec->binding == "readout") { decodeReadout(input);