From 30513e5b6a13daa87acaba23132f93e5ad5d0a7c Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Mon, 6 Nov 2017 09:15:20 +0100 Subject: [PATCH 01/12] add Cluster MCLabels fix gsl include directory WIP: Cluster MCLabels --- .../include/TPCReconstruction/BoxClusterer.h | 9 +- .../include/TPCReconstruction/Clusterer.h | 31 +++-- .../include/TPCReconstruction/ClustererTask.h | 21 +-- .../include/TPCReconstruction/HwClusterer.h | 64 ++++----- .../TPC/reconstruction/src/BoxClusterer.cxx | 4 +- .../TPC/reconstruction/src/ClustererTask.cxx | 20 ++- .../TPC/reconstruction/src/HwClusterer.cxx | 122 ++++++------------ cmake/O2Dependencies.cmake | 2 + macro/run_clus_tpc.C | 4 +- 9 files changed, 136 insertions(+), 141 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h index d7d6d8fd440b8..d28b9dde07d71 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h @@ -21,6 +21,9 @@ #include "TPCReconstruction/BoxCluster.h" #include "TPCBase/CalDet.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" + namespace o2{ namespace TPC { @@ -40,9 +43,11 @@ namespace o2{ /// Steer conversion of points to digits /// @param digits Container with TPC digits + /// @param mcDigitTruth MC Digit Truth container + /// @param mcClusterTruth MC Cluster Truth container /// @return Container with clusters - void Process(std::vector const & digits) override; - void Process(std::vector>& digits) override; + void Process(std::vector const & digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; + void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; /// Set a pedestal object void setPedestals(CalPad* pedestals) { mPedestals = pedestals; } diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h index b56771755a203..0e602e74b1907 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h @@ -19,14 +19,20 @@ #include "TPCReconstruction/ClusterContainer.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" + namespace o2{ namespace TPC { - + class Digit; /// \class Clusterer /// \brief Base Class for TPC clusterer class Clusterer { + protected: + using MCLabel = o2::dataformats::MCTruthContainer; + public: /// Default Constructor @@ -39,20 +45,22 @@ class Clusterer { /// \param minQMax Minimum peak charge for cluster /// \param requirePositiveCharge Positive charge is required /// \param requireNeighbouringPad Requires at least 2 adjecent pads with charge above threshold - Clusterer(int rowsMax, int padsMax, int timeBinsMax, int minQMax, + Clusterer(int rowsMax, int padsMax, int timeBinsMax, int minQMax, bool requirePositiveCharge, bool requireNeighbouringPad); - + /// Destructor virtual ~Clusterer() = default; - + /// Initialization function for clusterer virtual void Init() = 0; - + /// Processing all digits /// \param digits Container with TPC digits + /// @param mcDigitTruth MC Digit Truth container + /// @param mcClusterTruth MC Cluster Truth container /// \return Container with clusters - virtual void Process(std::vector const &digits) = 0; - virtual void Process(std::vector>& digits) = 0; + virtual void Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) = 0; + virtual void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) = 0; void setRowsMax(int val) { mRowsMax = val; }; void setPadsMax(int val) { mPadsMax = val; }; @@ -67,20 +75,19 @@ class Clusterer { float getMinQMax() const { return mMinQMax; }; bool hasRequirePositiveCharge() const { return mRequirePositiveCharge; }; bool hasRequireNeighbouringPad() const { return mRequireNeighbouringPad; }; - + protected: - - + int mRowsMax; ///< Maximum row number int mPadsMax; ///< Maximum pad number int mTimeBinsMax; ///< Maximum time bin float mMinQMax; ///< Minimun Qmax for cluster bool mRequirePositiveCharge; ///< If true, require charge > 0 bool mRequireNeighbouringPad; ///< If true, require 2+ pads minimum - + }; } } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index 54a08b04a1343..330a7b090f9a8 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -24,16 +24,18 @@ #include "TPCReconstruction/Clusterer.h" // for Clusterer #include "TPCReconstruction/BoxClusterer.h" // for Clusterer #include "TPCReconstruction/HwClusterer.h" // for Clusterer +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" #include namespace o2 { namespace TPC{ - + class ClustererTask : public FairTask{ public: ClustererTask(); ~ClustererTask() override; - + InitStatus Init() override; void Exec(Option_t *option) override; @@ -45,24 +47,24 @@ class ClustererTask : public FairTask{ }; }; - bool isClustererEnable(ClustererType type) const { + bool isClustererEnable(ClustererType type) const { switch (type) { case ClustererType::HW: return mHwClustererEnable; case ClustererType::Box: return mBoxClustererEnable; }; }; - Clusterer* getClusterer(ClustererType type) { + Clusterer* getClusterer(ClustererType type) { switch (type) { case ClustererType::HW: return mHwClusterer; case ClustererType::Box: return mBoxClusterer; }; }; - + BoxClusterer* getBoxClusterer() const { return mBoxClusterer; }; HwClusterer* getHwClusterer() const { return mHwClusterer; }; // Clusterer *GetClusterer() const { return fClusterer; } - + /// Switch for triggered / continuous readout /// \param isContinuous - false for triggered readout, true for continuous readout void setContinuousReadout(bool isContinuous); @@ -74,12 +76,15 @@ class ClustererTask : public FairTask{ BoxClusterer *mBoxClusterer; HwClusterer *mHwClusterer; - + std::vector const *mDigitsArray; + o2::dataformats::MCTruthContainer const *mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray // produced data containers std::vector *mClustersArray; std::vector *mHwClustersArray; - + o2::dataformats::MCTruthContainer mClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mClustersArrays + o2::dataformats::MCTruthContainer mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays + ClassDefOverride(ClustererTask, 1) }; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index 37a0fd491c57c..c5e64f362a3cd 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -16,13 +16,17 @@ #include "TPCReconstruction/Clusterer.h" #include "TPCReconstruction/HwCluster.h" -#include "TPCBase/CalDet.h" +#include "TPCBase/CalDet.h" + +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" #include +#include namespace o2{ namespace TPC { - + class ClusterContainer; class ClustererTask; class HwClusterFinder; @@ -32,6 +36,7 @@ class Digit; /// \class HwClusterer /// \brief Class for TPC HW cluster finding class HwClusterer : public Clusterer { + using MCLabel = o2::dataformats::MCTruthContainer; public: enum class Processing : int { Sequential, Parallel}; @@ -40,7 +45,7 @@ class HwClusterer : public Clusterer { /// \param processingType parallel or sequential /// \param globalTime value of first timebin /// \param cru Number of CRUs to process - /// \param minQDiff Min charge differece + /// \param minQDiff Min charge differece /// \param assignChargeUnique Avoid using same charge for multiple nearby clusters /// \param enableNoiseSim Enables the Noise simulation for empty pads (noise object has to be set) /// \param enablePedestalSubtraction Enables the Pedestal subtraction (pedestal object has to be set) @@ -50,31 +55,33 @@ class HwClusterer : public Clusterer { HwClusterer( std::vector *output, Processing processingType = Processing::Parallel, - int globalTime = 0, + int globalTime = 0, int cruMin = 0, - int cruMax = 360, + int cruMax = 360, float minQDiff = 0, bool assignChargeUnique = false, - bool enableNoiseSim = true, - bool enablePedestalSubtraction = true, - int padsPerCF = 8, - int timebinsPerCF = 8, + bool enableNoiseSim = true, + bool enablePedestalSubtraction = true, + int padsPerCF = 8, + int timebinsPerCF = 8, int cfPerRow = 0); - + /// Destructor ~HwClusterer() override; - + // Should this really be a public member? // Maybe better to just call by process void Init() override; - + /// Steer conversion of points to digits /// @param digits Container with TPC digits + /// @param mcDigitTruth MC Digit Truth container + /// @param mcClusterTruth MC Cluster Truth container /// @return Container with clusters - void Process(std::vector const &digits) override; - void Process(std::vector>& digits) override; + void Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) override; + void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) override; - void setProcessingType(Processing processing) { mProcessingType = processing; }; + void setProcessingType(Processing processing) { mProcessingType = processing; }; void setNoiseObject(CalDet* noiseObject) { mNoiseObject = noiseObject; }; void setPedestalObject(CalDet* pedestalObject) { mPedestalObject = pedestalObject; }; @@ -85,7 +92,7 @@ class HwClusterer : public Clusterer { void setCRUMin(int cru) { mCRUMin = cru; }; void setCRUMax(int cru) { mCRUMax = cru; }; - + private: // To be done /* BoxClusterer(const BoxClusterer &); */ @@ -103,26 +110,21 @@ class HwClusterer : public Clusterer { CalDet* iNoiseObject; CalDet* iPedestalObject; }; - + static void processDigits( - const std::vector>& digits, - const std::vector>& clusterFinder, - std::vector& cluster, + const std::vector>>>& digits, + const std::vector>& clusterFinder, + std::vector& cluster, CfConfig config); -// int iCRU, -// int maxRows, -// int maxPads, -// unsigned minTimeBin, -// unsigned maxTimeBin); - + void ProcessTimeBins(int iTimeBinMin, int iTimeBinMax); std::vector>> mClusterFinder; - std::vector>> mDigitContainer; + std::vector>>>> mDigitContainer; std::vector> mClusterStorage; - - Processing mProcessingType; + + Processing mProcessingType; int mGlobalTime; int mCRUMin; @@ -138,7 +140,7 @@ class HwClusterer : public Clusterer { int mLastTimebin; std::vector* mClusterArray; ///< Internal cluster storage - + CalDet* mNoiseObject; CalDet* mPedestalObject; }; @@ -146,4 +148,4 @@ class HwClusterer : public Clusterer { } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx index e56477347ece2..a9eed507b318d 100644 --- a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx @@ -145,7 +145,7 @@ void BoxClusterer::Init() } //________________________________________________________________________ -void BoxClusterer::Process(std::vector const &digits) +void BoxClusterer::Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) { mClusterArray->clear(); // check this @@ -181,7 +181,7 @@ void BoxClusterer::Process(std::vector const &digits) } //________________________________________________________________________ -void BoxClusterer::Process(std::vector>& digits) +void BoxClusterer::Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) { mClusterArray->clear(); diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 184d7ddc151b2..0829e7b9a651b 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -35,8 +35,11 @@ ClustererTask::ClustererTask() , mBoxClusterer(nullptr) , mHwClusterer(nullptr) , mDigitsArray(nullptr) + , mDigitMCTruthArray(nullptr) , mClustersArray(nullptr) , mHwClustersArray(nullptr) + , mClustersMCTruthArray() + , mHwClustersMCTruthArray() { } @@ -74,12 +77,20 @@ InitStatus ClustererTask::Init() return kERROR; } + mDigitMCTruthArray = mgr->InitObjectAs("TPCDigitMCTruth"); + if( !mDigitMCTruthArray ) { + LOG(ERROR) << "TPC MC truth not registered in the FairRootManager. Exiting ..." << FairLogger::endl; + return kERROR; + } + if (mBoxClustererEnable) { - + // Create and register output container mClustersArray = new std::vector; mgr->RegisterAny("TPCCluster", mClustersArray, kTRUE); + mgr->Register("TPCClusterMCTruth", "TPC", &mClustersMCTruthArray, kTRUE); + // create clusterer and pass output pointer mBoxClusterer = new BoxClusterer(mClustersArray); mBoxClusterer->Init(); @@ -90,6 +101,8 @@ InitStatus ClustererTask::Init() mHwClustersArray = new std::vector; mgr->RegisterAny("TPCClusterHW", mHwClustersArray, kTRUE); + mgr->Register("TPCClusterHWMCTruth", "TPC", &mHwClustersMCTruthArray, kTRUE); + // create clusterer and pass output pointer mHwClusterer = new HwClusterer(mHwClustersArray); mHwClusterer->setContinuousReadout(mIsContinuousReadout); @@ -109,12 +122,13 @@ void ClustererTask::Exec(Option_t *option) if (mBoxClustererEnable) { mClustersArray->clear(); - mBoxClusterer->Process(*mDigitsArray); + mBoxClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mClustersMCTruthArray); } if (mHwClustererEnable) { mHwClustersArray->clear(); - mHwClusterer->Process(*mDigitsArray); + mHwClustersMCTruthArray.clear(); + mHwClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mHwClustersMCTruthArray); LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl; } } diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 0bc3491eedeb2..bcd3118fc0c0e 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -19,7 +19,7 @@ #include "TPCBase/PadPos.h" #include "TPCBase/CRU.h" #include "TPCBase/PadSecPos.h" -#include "TPCBase/CalArray.h" +#include "TPCBase/CalArray.h" #include "FairLogger.h" #include "TMath.h" @@ -72,7 +72,7 @@ HwClusterer::~HwClusterer() //________________________________________________________________________ void HwClusterer::Init() { - LOG(DEBUG) << "Enter Initializer of HwClusterer" << FairLogger::endl; + LOG(DEBUG) << "Enter Initializer of HwClusterer" << FairLogger::endl; /* * initalize all cluster finder @@ -110,9 +110,9 @@ void HwClusterer::Init() mClusterStorage.resize(mCRUMax); - /* + /* * vector of digit vectors, one vector for each CRU (possible thread) to - * store there only those digits which are relevant for this particular + * store there only those digits which are relevant for this particular * CRU (thread) */ mDigitContainer.resize(mCRUMax); @@ -123,15 +123,10 @@ void HwClusterer::Init() //________________________________________________________________________ void HwClusterer::processDigits( - const std::vector>& digits, - const std::vector>& clusterFinder, + const std::vector>>>& digits, + const std::vector>& clusterFinder, std::vector& cluster, const CfConfig config) -// int iCRU, -// int maxRows, -// int maxPads, -// unsigned minTime, -// unsigned maxTime) { int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; if (timeDiff < 0) return; @@ -140,13 +135,6 @@ void HwClusterer::processDigits( // std::cout << "thread " << this_id << " started.\n"; // g_display_mutex.unlock(); -// float iAllBins[maxTime][maxPads]; -// for (float (&timebin)[maxPads] : iAllBins) { -// for (float &pad : timebin) { -// pad = 0.0; -// } -// } - float iAllBins[timeDiff][config.iMaxPads]; for (int iRow = 0; iRow < config.iMaxRows; iRow++){ @@ -170,10 +158,10 @@ void HwClusterer::processDigits( * fill in digits */ for (auto& digit : digits[iRow]){ - const Int_t iTime = digit->getTimeStamp(); - const Int_t iPad = digit->getPad() + 2; // offset to have 2 empty pads on the "left side" - const Float_t charge = digit->getChargeFloat(); - + const Int_t iTime = digit.first->getTimeStamp(); + const Int_t iPad = digit.first->getPad() + 2; // offset to have 2 empty pads on the "left side" + const Float_t charge = digit.first->getChargeFloat(); + // std::cout << iCRU << " " << iRow << " " << iPad << " " << iTime << " (" << iTime-minTime << "," << timeDiff << ") " << charge << std::endl; iAllBins[iTime-config.iMinTimeBin][iPad] += charge; if (config.iEnablePedestalSubtraction && config.iPedestalObject != nullptr) { @@ -182,7 +170,7 @@ void HwClusterer::processDigits( iAllBins[iTime-config.iMinTimeBin][iPad] -= pedestal; } } - + /* * copy data to cluster finders */ @@ -190,17 +178,17 @@ void HwClusterer::processDigits( const Short_t iTimebinsPerCF = clusterFinder[iRow][0]->getNtimebins(); std::vector::const_reverse_iterator> cfWithCluster; unsigned time,pad; - for (time = 0; time < timeDiff; ++time){ + for (time = 0; time < timeDiff; ++time){ // ordering important!! for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { const Short_t cf = pad / (iPadsPerCF-2-2); clusterFinder[iRow][cf]->AddTimebin(&iAllBins[time][pad],time+config.iMinTimeBin,(config.iMaxPads-pad)>=iPadsPerCF?iPadsPerCF:(config.iMaxPads-pad)); } - + /* * search for clusters and store reference to CF if one was found */ if (clusterFinder[iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { - /* + /* * ordering is important: from right to left, so that the CFs could inform each other if cluster was found */ for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { @@ -225,7 +213,7 @@ void HwClusterer::processDigits( * search for clusters and store reference to CF if one was found */ if (clusterFinder[iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { - /* + /* * ordering is important: from right to left, so that the CFs could inform each other if cluster was found */ for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { @@ -239,8 +227,8 @@ void HwClusterer::processDigits( (*rit)->setTimebinsAfterLastProcessing(0); } } - - /* + + /* * collect found cluster */ for (std::vector::const_reverse_iterator &cf_it : cfWithCluster) { @@ -251,15 +239,6 @@ void HwClusterer::processDigits( (*cf_it)->clearClusterContainer(); } -// /* -// * remove digits again from storage -// */ -// for (std::vector::const_iterator it = digits[iRow].begin(); it != digits[iRow].end(); ++it){ -// const Int_t iTime = (*it)->getTimeStamp(); -// const Int_t iPad = (*it)->getPad() + 2; // offset to have 2 empty pads on the "left side" -// -// iAllBins[iTime-minTime][iPad] = 0.0; -// } } // g_display_mutex.lock(); @@ -268,11 +247,12 @@ void HwClusterer::processDigits( } //________________________________________________________________________ -void HwClusterer::Process(std::vector const &digits) +void HwClusterer::Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth + ) { mClusterArray->clear(); - /* + /* * clear old storages */ for (auto& cs : mClusterStorage) cs.clear(); @@ -280,68 +260,68 @@ void HwClusterer::Process(std::vector const &digits) for (auto& dcc : dc) dcc.clear(); } -// int iCRU; -// int iRow; -// int iPad; -// float charge; int iTimeBin; int iTimeBinMin = (mIsContinuousReadout)?mLastTimebin + 1 : 0; //int iTimeBinMin = mLastTimebin + 1; int iTimeBinMax = mLastTimebin; - /* + /* * Loop over digits */ + int digitIndex = 0; for (const auto& digit : digits) { + /* * add current digit to storage */ -// iCRU = digit->getCRU(); -// iRow = digit->getRow(); -// iPad = digit->getPad(); -// charge = digit->getChargeFloat(); -// if ((cru == 179)) {// && iRow == 5)){ -// printf("hw: digi: %d, %d, %d, %d, %.2f\n", cru, iRow, iPad, iTimeBin, charge); -// } iTimeBin = digit.getTimeStamp(); if (iTimeBin < iTimeBinMin) continue; iTimeBinMax = std::max(iTimeBinMax,iTimeBin); - mDigitContainer[digit.getCRU()][digit.getRow()].push_back(&digit); + if (mcDigitTruth == nullptr) + mDigitContainer[digit.getCRU()][digit.getRow()].push_back(std::make_pair(&digit,nullptr)); + else + mDigitContainer[digit.getCRU()][digit.getRow()].push_back(std::make_pair(&digit,mcDigitTruth->getLabels(digitIndex))); + ++digitIndex; } ProcessTimeBins(iTimeBinMin, iTimeBinMax); } //________________________________________________________________________ -void HwClusterer::Process(std::vector>& digits) +void HwClusterer::Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) { mClusterArray->clear(); - /* + /* * clear old storages */ for (auto& cs : mClusterStorage) cs.clear(); for (auto& dc : mDigitContainer ) { - for (auto& dcc : dc) dcc.clear(); + for (auto& dcc : dc) dcc.clear(); } int iTimeBin; int iTimeBinMin = (mIsContinuousReadout)?mLastTimebin + 1 : 0; int iTimeBinMax = mLastTimebin; - /* + /* * Loop over digits */ + int digitIndex = 0; for (auto& digit_ptr : digits) { Digit* digit = digit_ptr.get(); + /* * add current digit to storage */ - iTimeBin = digit->getTimeStamp(); if (iTimeBin < iTimeBinMin) continue; iTimeBinMax = std::max(iTimeBinMax,iTimeBin); - mDigitContainer[digit->getCRU()][digit->getRow()].push_back(digit); + if (mcDigitTruth == nullptr) + mDigitContainer[digit->getCRU()][digit->getRow()].push_back(std::make_pair(digit,nullptr)); + else + mDigitContainer[digit->getCRU()][digit->getRow()].push_back(std::make_pair(digit,mcDigitTruth->getLabels(digitIndex))); + ++digitIndex; } ProcessTimeBins(iTimeBinMin, iTimeBinMax); @@ -386,12 +366,6 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) std::ref(mClusterFinder[iCRU]), // cluster finder for individual CRUs std::ref(mClusterStorage[iCRU]), // container to store found clusters cfConfig -// iCRU, // current CRU for deb. purposes -// mRowsMax, // max. numbers of rows per CRU -// mPadsMax+2+2, // max. numbers of pads in each row (+2 empty ones on each side) -// iTimeBinMin, // Min timebin of digit -// iTimeBinMax // Max timebin of digits - ); else { processDigits( @@ -399,11 +373,6 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) std::ref(mClusterFinder[iCRU]), std::ref(mClusterStorage[iCRU]), cfConfig); -// iCRU, -// mRowsMax, -// mPadsMax+2+2, -// iTimeBinMin, -// iTimeBinMax); } } @@ -419,23 +388,14 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) * collect clusters from individual cluster finder */ for (auto& cc : mClusterStorage) { - if (cc.size() != 0) { + if (cc.size() != 0) { for (auto& c : cc){ ClusterContainer::AddCluster(mClusterArray, c.getCRU(),c.getRow(),c.getQ(),c.getQmax(), c.getPadMean(),c.getTimeMean(),c.getPadSigma(),c.getTimeSigma()); -//// if (c.getPadSigma() > 0.45 && c.getPadSigma() < 0.5) { -// if ((c.getCRU() == 179)){// && c.getRow() == 5)){// && (int)c.getPadMean() == 103 && (int)c.getTimeMean() == 170) || -//// (iCRU == 256 && iRow == 10 && (int)c.getPadMean() == 27 && (int)c.getTimeMean() == 181) ) { -// std::cout << "HwCluster - "; -//// c.Print(std::cout); -// c.PrintDetails(std::cout); -// std::cout << std::endl; -// } -// c.Print(); } } } mLastTimebin = iTimeBinMax; } - + diff --git a/cmake/O2Dependencies.cmake b/cmake/O2Dependencies.cmake index 672a11657bbc3..ed6cf87d51816 100644 --- a/cmake/O2Dependencies.cmake +++ b/cmake/O2Dependencies.cmake @@ -740,6 +740,8 @@ o2_define_bucket( ${CMAKE_SOURCE_DIR}/Detectors/Passive/include ${CMAKE_SOURCE_DIR}/Detectors/TPC/base/include ${CMAKE_SOURCE_DIR}/Detectors/TPC/simulation/include + ${CMAKE_SOURCE_DIR}/DataFormats/simulation/include + ${MS_GSL_INCLUDE_DIR} ) o2_define_bucket( diff --git a/macro/run_clus_tpc.C b/macro/run_clus_tpc.C index dbca39bc0724a..c04a4872fa870 100644 --- a/macro/run_clus_tpc.C +++ b/macro/run_clus_tpc.C @@ -56,8 +56,8 @@ void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isConti // Initialize everything run->Init(); -// clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Parallel); - clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Sequential); + clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Parallel); +// clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Sequential); // Start simulation timer.Start(); From c66033e32a240de8ffad61b2318b20c831eae777 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Tue, 7 Nov 2017 11:30:13 +0100 Subject: [PATCH 02/12] replaced BoxCluster with Cluster, additional information were never used --- .../include/TPCReconstruction/BoxClusterer.h | 40 +++++++++---------- .../include/TPCReconstruction/ClustererTask.h | 2 +- .../reconstruction/macro/RawClusterFinder.C | 2 +- .../TPC/reconstruction/src/BoxClusterer.cxx | 6 +-- .../TPC/reconstruction/src/ClustererTask.cxx | 2 +- 5 files changed, 25 insertions(+), 27 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h index d28b9dde07d71..57dfac5c13d2a 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h @@ -18,29 +18,29 @@ #include "Rtypes.h" #include "TPCReconstruction/Clusterer.h" -#include "TPCReconstruction/BoxCluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCBase/CalDet.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "SimulationDataFormat/MCCompLabel.h" namespace o2{ - + namespace TPC { - + class ClusterContainer; - + class BoxClusterer : public Clusterer { public: - BoxClusterer(std::vector *output); - + BoxClusterer(std::vector *output); + /// Destructor ~BoxClusterer() override; - + // Should this really be a public member? // Maybe better to just call by process void Init() override; - + /// Steer conversion of points to digits /// @param digits Container with TPC digits /// @param mcDigitTruth MC Digit Truth container @@ -48,7 +48,7 @@ namespace o2{ /// @return Container with clusters void Process(std::vector const & digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; - + /// Set a pedestal object void setPedestals(CalPad* pedestals) { mPedestals = pedestals; } @@ -56,21 +56,21 @@ namespace o2{ // To be done /* BoxClusterer(const BoxClusterer &); */ /* BoxClusterer &operator=(const BoxClusterer &); */ - + void FindLocalMaxima(const Int_t iCRU); void CleanArrays(); void GetPadAndTimeBin(Int_t bin, Short_t& iPad, Short_t& iTimeBin); - Int_t Update(const Int_t iCRU, const Int_t iRow, const Int_t iPad, + Int_t Update(const Int_t iCRU, const Int_t iRow, const Int_t iPad, const Int_t iTimeBin, Float_t signal); Float_t GetQ(const Float_t* adcArray, const Short_t pad, - const Short_t time, Short_t& timeMin, Short_t& timeMax, + const Short_t time, Short_t& timeMin, Short_t& timeMax, Short_t& padMin, Short_t& padMax) const; - Bool_t UpdateCluster(Float_t charge, Int_t deltaPad, Int_t deltaTime, - Float_t& qTotal, Double_t& meanPad, - Double_t& sigmaPad, Double_t& meanTime, + Bool_t UpdateCluster(Float_t charge, Int_t deltaPad, Int_t deltaTime, + Float_t& qTotal, Double_t& meanPad, + Double_t& sigmaPad, Double_t& meanTime, Double_t& sigmaTime); - - + + // // Expand buffer // @@ -79,12 +79,12 @@ namespace o2{ Int_t* mAllNSigBins; //!* mClusterArray; ///< Internal cluster storage - + std::vector* mClusterArray; ///< Internal cluster storage + ClassDefNV(BoxClusterer, 1); }; } } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index 330a7b090f9a8..b70ca077e0db6 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -80,7 +80,7 @@ class ClustererTask : public FairTask{ std::vector const *mDigitsArray; o2::dataformats::MCTruthContainer const *mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray // produced data containers - std::vector *mClustersArray; + std::vector *mClustersArray; std::vector *mHwClustersArray; o2::dataformats::MCTruthContainer mClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mClustersArrays o2::dataformats::MCTruthContainer mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays diff --git a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C index 9d7e335504728..ca14a1397894e 100644 --- a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C +++ b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C @@ -144,7 +144,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt // ===| output file and container |=========================================== std::vector arrCluster; - std::vector *arrClusterBox = nullptr; + std::vector *arrClusterBox = nullptr; float cherenkovValue = 0.; int runNumber = 0; diff --git a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx index a9eed507b318d..dc5e1b1b58cc2 100644 --- a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx @@ -93,7 +93,6 @@ #include "TPCReconstruction/BoxClusterer.h" #include "TPCBase/Digit.h" #include "TPCReconstruction/ClusterContainer.h" -#include "TPCReconstruction/BoxCluster.h" #include "FairLogger.h" #include "TMath.h" @@ -104,7 +103,7 @@ ClassImp(o2::TPC::BoxClusterer) using namespace o2::TPC; //________________________________________________________________________ -BoxClusterer::BoxClusterer(std::vector *output): +BoxClusterer::BoxClusterer(std::vector *output): Clusterer(), mClusterArray(output), mAllBins(nullptr), @@ -352,8 +351,7 @@ void BoxClusterer::FindLocalMaxima(const Int_t iCRU) Short_t nPad = maxP-minP+1; Short_t nTimeBins = maxT-minT+1; Short_t size = 10*nPad+nTimeBins; - BoxCluster* cluster = ClusterContainer::AddCluster(mClusterArray, iCRU, iRow, qTot, qMax, meanP, meanT,sigmaP, sigmaT); - cluster->setBoxParameters(pad, timebin, size); + Cluster* cluster = ClusterContainer::AddCluster(mClusterArray, iCRU, iRow, qTot, qMax, meanP, meanT,sigmaP, sigmaT); // if ((iCRU == 179)) {// && iRow == 5)){// && (int)meanP == 103 && (int)meanT == 170) || //// (iCRU == 256 && iRow == 10 && (int)meanP == 27 && (int)meanT == 181) ) { diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 0829e7b9a651b..745a2b1d6b9eb 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -86,7 +86,7 @@ InitStatus ClustererTask::Init() if (mBoxClustererEnable) { // Create and register output container - mClustersArray = new std::vector; + mClustersArray = new std::vector; mgr->RegisterAny("TPCCluster", mClustersArray, kTRUE); mgr->Register("TPCClusterMCTruth", "TPC", &mClustersMCTruthArray, kTRUE); From ed7fd886de2b1d1850eccdaca2cf11969fd24753 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Tue, 7 Nov 2017 13:09:00 +0100 Subject: [PATCH 03/12] replaced HwCluster with Cluster, additional information were never used --- Detectors/TPC/reconstruction/CMakeLists.txt | 4 - .../include/TPCReconstruction/ClustererTask.h | 2 +- .../include/TPCReconstruction/HwCluster.h | 8 +- .../TPCReconstruction/HwClusterFinder.h | 10 +- .../include/TPCReconstruction/HwClusterer.h | 14 +-- .../include/TPCReconstruction/TPCCATracking.h | 6 +- .../reconstruction/macro/RawClusterFinder.C | 2 +- .../TPC/reconstruction/macro/convertTracks.C | 2 +- .../TPC/reconstruction/macro/runCATracking.C | 3 +- .../reconstruction/macro/testClustererData.C | 2 +- .../TPC/reconstruction/macro/testTracks.C | 4 +- .../TPC/reconstruction/src/ClustererTask.cxx | 2 +- .../reconstruction/src/HwClusterFinder.cxx | 119 +++++++++++++----- .../TPC/reconstruction/src/HwClusterer.cxx | 15 +-- .../TPC/reconstruction/src/TPCCATracking.cxx | 3 +- .../src/TPCReconstructionLinkDef.h | 2 - 16 files changed, 121 insertions(+), 77 deletions(-) diff --git a/Detectors/TPC/reconstruction/CMakeLists.txt b/Detectors/TPC/reconstruction/CMakeLists.txt index eef0425f832ac..1d0986175d38b 100644 --- a/Detectors/TPC/reconstruction/CMakeLists.txt +++ b/Detectors/TPC/reconstruction/CMakeLists.txt @@ -7,7 +7,6 @@ link_directories( ${LINK_DIRECTORIES}) set(SRCS src/AdcClockMonitor.cxx - src/BoxCluster.cxx src/BoxClusterer.cxx src/Cluster.cxx src/Clusterer.cxx @@ -16,7 +15,6 @@ set(SRCS src/GBTFrame.cxx src/GBTFrameContainer.cxx src/HalfSAMPAData.cxx - src/HwCluster.cxx src/HwClusterer.cxx src/HwClusterFinder.cxx src/HwFixedPoint.cxx @@ -31,7 +29,6 @@ set(SRCS set(HEADERS include/${MODULE_NAME}/AdcClockMonitor.h - include/${MODULE_NAME}/BoxCluster.h include/${MODULE_NAME}/BoxClusterer.h include/${MODULE_NAME}/Cluster.h include/${MODULE_NAME}/Clusterer.h @@ -40,7 +37,6 @@ set(HEADERS include/${MODULE_NAME}/GBTFrame.h include/${MODULE_NAME}/GBTFrameContainer.h include/${MODULE_NAME}/HalfSAMPAData.h - include/${MODULE_NAME}/HwCluster.h include/${MODULE_NAME}/HwClusterer.h include/${MODULE_NAME}/HwClusterFinder.h include/${MODULE_NAME}/HwFixedPoint.h diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index b70ca077e0db6..14e4d8ee052f6 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -81,7 +81,7 @@ class ClustererTask : public FairTask{ o2::dataformats::MCTruthContainer const *mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray // produced data containers std::vector *mClustersArray; - std::vector *mHwClustersArray; + std::vector *mHwClustersArray; o2::dataformats::MCTruthContainer mClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mClustersArrays o2::dataformats::MCTruthContainer mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h index 7100df03c1540..1c42ab1054abe 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h @@ -53,10 +53,10 @@ class HwCluster : public Cluster { /// \param other HwCluster to be copied HwCluster(const HwCluster& other); - short getPad() const { return mPad; } - short getTime() const { return mTime; } - short getSizeP() const { return mSizeP; } - short getSizeT() const { return mSizeT; } +// short getPad() const { return mPad; } +// short getTime() const { return mTime; } +// short getSizeP() const { return mSizeP; } +// short getSizeT() const { return mSizeT; } /// Set all cluster data /// \param cru CRU diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index 6ab8aac820712..d606c55c8456f 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -20,7 +20,7 @@ namespace o2{ namespace TPC { -class HwCluster; +class Cluster; /// \class HwClusterFinder /// \brief Class for TPC HW cluster finder @@ -88,10 +88,10 @@ class HwClusterFinder { bool getRequirePositiveCharge() const { return mRequirePositiveCharge; } bool getRequireNeighbouringPad() const { return mRequireNeighbouringPad; } bool getRequireNeighbouringTimebin() const { return mRequireNeighbouringTimebin; } - bool getAutoProcessing() const { return mAutoProcessing; } + bool getAutoProcessing() const { return mAutoProcessing; } bool getmAssignChargeUnique() const { return mAssignChargeUnique; } - HwClusterFinder* getNextCF() const { return mNextCF; } - std::vector* getClusterContainer() { return &clusterContainer; } + HwClusterFinder* getNextCF() const { return mNextCF; } + std::vector* getClusterContainer() { return &clusterContainer; } // Setter functions void setTimebinsAfterLastProcessing(int val) { mTimebinsAfterLastProcessing = val; }; @@ -131,7 +131,7 @@ class HwClusterFinder { void printCluster(short time, short pad); // local variables - std::vector clusterContainer; + std::vector clusterContainer; int mTimebinsAfterLastProcessing; float** mData; float** tmpCluster; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index c5e64f362a3cd..efdef3f219c6e 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -15,7 +15,7 @@ #define ALICEO2_TPC_HWClusterer_H_ #include "TPCReconstruction/Clusterer.h" -#include "TPCReconstruction/HwCluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCBase/CalDet.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -27,10 +27,8 @@ namespace o2{ namespace TPC { -class ClusterContainer; class ClustererTask; class HwClusterFinder; -class HwCluster; class Digit; /// \class HwClusterer @@ -53,8 +51,8 @@ class HwClusterer : public Clusterer { /// \param timebinsPerCF Timebins per cluster finder /// \param cfPerRow Number of cluster finder in each row HwClusterer( - std::vector *output, - Processing processingType = Processing::Parallel, + std::vector *output, + Processing processingType = Processing::Parallel, int globalTime = 0, int cruMin = 0, int cruMax = 360, @@ -114,7 +112,7 @@ class HwClusterer : public Clusterer { static void processDigits( const std::vector>>>& digits, const std::vector>& clusterFinder, - std::vector& cluster, + std::vector& cluster, CfConfig config); void ProcessTimeBins(int iTimeBinMin, int iTimeBinMax); @@ -122,7 +120,7 @@ class HwClusterer : public Clusterer { std::vector>> mClusterFinder; std::vector>>>> mDigitContainer; - std::vector> mClusterStorage; + std::vector> mClusterStorage; Processing mProcessingType; @@ -139,7 +137,7 @@ class HwClusterer : public Clusterer { int mCfPerRow; int mLastTimebin; - std::vector* mClusterArray; ///< Internal cluster storage + std::vector* mClusterArray; ///< Internal cluster storage CalDet* mNoiseObject; CalDet* mPedestalObject; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/TPCCATracking.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/TPCCATracking.h index 258d8bcefc4e7..b0b20a08cec8b 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/TPCCATracking.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/TPCCATracking.h @@ -16,7 +16,7 @@ #include #include -#include "TPCReconstruction/HwCluster.h" +#include "TPCReconstruction/Cluster.h" class TChain; class AliHLTTPCCAO2Interface; class AliHLTTPCCAClusterData; @@ -37,11 +37,11 @@ class TPCCATracking int initialize(const char* options = nullptr); void deinitialize(); - int runTracking(const std::vector* inputClusters, std::vector* outputTracks) {return runTracking(nullptr, inputClusters, outputTracks);} + int runTracking(const std::vector* inputClusters, std::vector* outputTracks) {return runTracking(nullptr, inputClusters, outputTracks);} int runTracking(TChain* inputClusters, std::vector* outputTracks) {return runTracking(inputClusters, nullptr, outputTracks);} private: - int runTracking(TChain* inputClustersChain, const std::vector* inputClustersArray, std::vector* outputTracks); + int runTracking(TChain* inputClustersChain, const std::vector* inputClustersArray, std::vector* outputTracks); std::unique_ptr mTrackingCAO2Interface; //Pointer to Interface class in HLT O2 CA Tracking library. //The tracking code itself is not included in the O2 package, but contained in the CA library. diff --git a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C index ca14a1397894e..f4bb0214c6d45 100644 --- a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C +++ b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C @@ -143,7 +143,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt } // ===| output file and container |=========================================== - std::vector arrCluster; + std::vector arrCluster; std::vector *arrClusterBox = nullptr; float cherenkovValue = 0.; int runNumber = 0; diff --git a/Detectors/TPC/reconstruction/macro/convertTracks.C b/Detectors/TPC/reconstruction/macro/convertTracks.C index 95c0badb5635b..b28294f3ebc0a 100644 --- a/Detectors/TPC/reconstruction/macro/convertTracks.C +++ b/Detectors/TPC/reconstruction/macro/convertTracks.C @@ -58,7 +58,7 @@ void convertTracks(TString inputBinaryFile, TString inputClusters, TString chere float cherenkovValue = 0.; int runNumber = 0; - std::vector *clusters=nullptr; + std::vector *clusters=nullptr; c.SetBranchAddress("TPCClusterHW", &clusters); //c.SetBranchAddress("TPC_Cluster", &clusters); c.SetBranchAddress("cherenkovValue", &cherenkovValue); diff --git a/Detectors/TPC/reconstruction/macro/runCATracking.C b/Detectors/TPC/reconstruction/macro/runCATracking.C index 1f4274d100e2c..96357141511ab 100644 --- a/Detectors/TPC/reconstruction/macro/runCATracking.C +++ b/Detectors/TPC/reconstruction/macro/runCATracking.C @@ -21,7 +21,6 @@ #include "TClonesArray.h" #include "TPCReconstruction/Cluster.h" -#include "TPCReconstruction/HwCluster.h" #include "TPCReconstruction/TPCCATracking.h" #include "TPCReconstruction/TrackTPC.h" #include "DetectorsBase/Track.h" @@ -60,7 +59,7 @@ void runCATracking(TString filename, TString outputFile, TString options="", boo } tout.Fill(); } else { - std::vector* clusters=nullptr; + std::vector* clusters=nullptr; c.SetBranchAddress("TPCClusterHW", &clusters); // ===| event ranges |======================================================== diff --git a/Detectors/TPC/reconstruction/macro/testClustererData.C b/Detectors/TPC/reconstruction/macro/testClustererData.C index 6660d6bf691ad..ec7305140ca76 100644 --- a/Detectors/TPC/reconstruction/macro/testClustererData.C +++ b/Detectors/TPC/reconstruction/macro/testClustererData.C @@ -42,7 +42,7 @@ void testClustererData(Int_t maxEvents=50, TString fileInfo="GBTx0_Run005:0:0;GB int mTimeBinsPerCall=500; // ===| output file and container |=========================================== - std::vector arrCluster; + std::vector arrCluster; TFile fout(outputFileName,"recreate"); TTree t("clusters","clusters"); t.Branch("cl", &arrCluster); diff --git a/Detectors/TPC/reconstruction/macro/testTracks.C b/Detectors/TPC/reconstruction/macro/testTracks.C index 33e88aa9cf185..ca7969d0f4c0d 100644 --- a/Detectors/TPC/reconstruction/macro/testTracks.C +++ b/Detectors/TPC/reconstruction/macro/testTracks.C @@ -22,7 +22,7 @@ #include "TPCReconstruction/TrackTPC.h" #include "DetectorsBase/Track.h" #include "TPCSimulation/Cluster.h" -#include "TPCSimulation/HwCluster.h" +#include "TPCSimulation/Cluster.h" #include "TPCBase/Mapper.h" #endif @@ -44,7 +44,7 @@ void testTracks(int checkEvent = 0, TFile *clusFile = TFile::Open(clusterFile.data()); TTree *clusterTree = (TTree *)gDirectory->Get("cbmsim"); - std::vector *clusters=nullptr; + std::vector *clusters=nullptr; clusterTree->SetBranchAddress("TPCClusterHW",&clusters); TGraph *grClusters = new TGraph(); diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 745a2b1d6b9eb..07c1d78a575b0 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -98,7 +98,7 @@ InitStatus ClustererTask::Init() if (mHwClustererEnable) { // Register output container - mHwClustersArray = new std::vector; + mHwClustersArray = new std::vector; mgr->RegisterAny("TPCClusterHW", mHwClustersArray, kTRUE); mgr->Register("TPCClusterHWMCTruth", "TPC", &mHwClustersMCTruthArray, kTRUE); diff --git a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx index 608a24d8e0d95..f89725844671b 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx @@ -14,7 +14,7 @@ #include "TPCReconstruction/HwClusterFinder.h" #include "TPCReconstruction/ClusterContainer.h" -#include "TPCReconstruction/HwCluster.h" +#include "TPCReconstruction/Cluster.h" #include "FairLogger.h" @@ -25,7 +25,7 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterFinder::HwClusterFinder( short cru, short row, short id, - short padOffset, short pads, short timebins, + short padOffset, short pads, short timebins, float diffThreshold, float chargeThreshold, bool requirePositiveCharge) : mGlobalTimeOfLast(0) @@ -50,7 +50,7 @@ HwClusterFinder::HwClusterFinder( , mZeroTimebin(nullptr) , mNextCF(nullptr) { - if (mPads < mClusterSizePads) { + if (mPads < mClusterSizePads) { LOG(ERROR) << "Given width in pad direction is smaller than cluster size in pad direction." << " width in pad direction was increased to cluster size." << FairLogger::endl; mPads = mClusterSizePads; @@ -68,7 +68,7 @@ HwClusterFinder::HwClusterFinder( short t,p; mData = new float*[mTimebins]; for (t = 0; t < mTimebins; ++t){ - mData[t] = new float[mPads]; + mData[t] = new float[mPads]; for (p = 0; p < mPads; ++p){ mData[t][p] = 0; } @@ -105,7 +105,7 @@ HwClusterFinder::HwClusterFinder(const HwClusterFinder& other) short t,p; mData = new float*[mTimebins]; for (t = 0; t < mTimebins; ++t){ - mData[t] = new float[mPads]; + mData[t] = new float[mPads]; for (p = 0; p < mPads; ++p){ mData[t][p] = other.mData[t][p]; } @@ -189,9 +189,9 @@ bool HwClusterFinder::findCluster() // o o|o_o_o_o|o o <- tMin // o o o o o o o o // o o o o o o o o <- 0 - // // - + // + // // In time direction // @@ -208,6 +208,15 @@ bool HwClusterFinder::findCluster() short delPm = -2; short delPp = 2; + double qMax; + double qTot; + double charge; + double meanP, meanT; + double sigmaP, sigmaT; + short minP, minT; + short maxP, maxT; + short deltaP, deltaT; + short clusterSize; // // peak finding // @@ -245,11 +254,11 @@ bool HwClusterFinder::findCluster() // printf("## cluster found at t=%d, p=%d (in CF %d in row %d of CRU %d)\n",t,p,mId,mRow,mCRU); // printf("##\n"); // printCluster(t,p); - + // // cluster was found!! // - + // prepare temp storage for (tt=0; tt older tmpCluster[4][2] = chargeForCluster(&mData[t+2][p ],&mData[t+1][p ]); // t+X -> newer tmpCluster[2][0] = chargeForCluster(&mData[t ][p-2],&mData[t ][p-1]); tmpCluster[2][4] = chargeForCluster(&mData[t ][p+2],&mData[t ][p+1]); - - + + // The cells of the corners have 3 neighbours. // o o o o // o i i o - // C + // C // o i i o // o o o o - + // bottom left tmpCluster[3][0] = chargeForCluster(&mData[t+1][p-2],&mData[t+1][p-1]); tmpCluster[4][0] = chargeForCluster(&mData[t+2][p-2],&mData[t+1][p-1]); @@ -318,21 +326,70 @@ bool HwClusterFinder::findCluster() tmpCluster[0][0] = chargeForCluster(&mData[t-2][p-2],&mData[t-1][p-1]); tmpCluster[1][0] = chargeForCluster(&mData[t-1][p-2],&mData[t-1][p-1]); -// if ((mCRU == 179 && mRow == 1 && p+mPadOffset == 103 && mGlobalTimeOfLast-(mTimebins-1)+t == 170)/* || -// (mCRU == 256 && mRow == 10 && p+mPadOffset == 27 && mGlobalTimeOfLast-(mTimebins-1)+t == 181)*/ ) { -// PrintLocalStorage(); -// } + // + // calculate cluster Properties + // + + qMax = tmpCluster[2][2]; + qTot = 0; + meanP = 0; + meanT = 0; + sigmaP = 0; + sigmaT = 0; + minP = mClusterSizePads; + minT = mClusterSizeTime; + maxP = 0; + maxT = 0; + clusterSize = 0; + for (tt = 0; tt < mClusterSizeTime; ++tt) { + deltaT = tt - mClusterSizeTime/2; + for (pp = 0; pp < mClusterSizePads; ++pp) { + deltaP = pp - mClusterSizePads/2; + + charge = tmpCluster[tt][pp]; + + qTot += charge; + + meanP += charge * deltaP; + meanT += charge * deltaT; + + sigmaP += charge * deltaP*deltaP; + sigmaT += charge * deltaT*deltaT; + + if (charge > 0) { + minP = std::min(minP,pp); maxP = std::max(maxP,pp); + minT = std::min(minT,tt); maxT = std::max(maxT,tt); + } + + } + } + + clusterSize = (maxP-minP+1)*10 + (maxT-minT+1); + + if (qTot > 0) { + meanP /= qTot; + meanT /= qTot; + sigmaP /= qTot; + sigmaT /= qTot; + + sigmaP = std::sqrt(sigmaP - (meanP*meanP)); + sigmaT = std::sqrt(sigmaT - (meanT*meanT)); + + meanP += p+mPadOffset; + meanT += mGlobalTimeOfLast-(mTimebins-1)+t; + } ++foundNclusters; - clusterContainer.emplace_back(mCRU, mRow, mClusterSizePads, mClusterSizeTime, tmpCluster,p+mPadOffset,mGlobalTimeOfLast-(mTimebins-1)+t); +// clusterContainer.emplace_back(mCRU, mRow, mClusterSizePads, mClusterSizeTime, tmpCluster,p+mPadOffset,mGlobalTimeOfLast-(mTimebins-1)+t); + clusterContainer.emplace_back(mCRU, mRow, qTot, qMax, meanP, sigmaP, meanT, sigmaT); if (mAssignChargeUnique) { - if (p < (pMin+4)) { + if (p < (pMin+4)) { // If the cluster peak is in one of the 6 leftmost pads, the Cluster Finder // on the left has to know about it to ignore the already used pads. if (mNextCF != nullptr) mNextCF->clusterAlreadyUsed(t,p+mPadOffset,tmpCluster); } - + // // subtract found cluster from storage @@ -395,7 +452,7 @@ void HwClusterFinder::clusterAlreadyUsed(short time, short pad, float** cluster) if (t < 0 || t >= mTimebins) continue; for (p=localPad-2; p<=localPad+2; ++p){ if (p < 0 || p >= mPads) continue; - + mData[t][p] -= cluster[t-time+2][p-localPad+2]; } } diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index bcd3118fc0c0e..2956a5ae3fc7f 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -11,10 +11,8 @@ /// \file AliTPCUpgradeHwClusterer.cxx /// \brief Hwclusterer for the TPC -#include "TPCReconstruction/HwCluster.h" #include "TPCReconstruction/HwClusterer.h" #include "TPCReconstruction/HwClusterFinder.h" -#include "TPCReconstruction/ClusterContainer.h" #include "TPCBase/Digit.h" #include "TPCBase/PadPos.h" #include "TPCBase/CRU.h" @@ -33,7 +31,7 @@ std::mutex g_display_mutex; using namespace o2::TPC; //________________________________________________________________________ -HwClusterer::HwClusterer(std::vector *output, Processing processingType, int globalTime, int cruMin, int cruMax, float minQDiff, +HwClusterer::HwClusterer(std::vector *output, Processing processingType, int globalTime, int cruMin, int cruMax, float minQDiff, bool assignChargeUnique, bool enableNoiseSim, bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF, int cfPerRow) : Clusterer() , mClusterArray(output) @@ -125,7 +123,7 @@ void HwClusterer::Init() void HwClusterer::processDigits( const std::vector>>>& digits, const std::vector>& clusterFinder, - std::vector& cluster, + std::vector& cluster, const CfConfig config) { int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; @@ -232,8 +230,8 @@ void HwClusterer::processDigits( * collect found cluster */ for (std::vector::const_reverse_iterator &cf_it : cfWithCluster) { - std::vector* cc = (*cf_it)->getClusterContainer(); - for (HwCluster& c : *cc){ + std::vector* cc = (*cf_it)->getClusterContainer(); + for (Cluster& c : *cc){ cluster.push_back(c); } (*cf_it)->clearClusterContainer(); @@ -389,9 +387,8 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) */ for (auto& cc : mClusterStorage) { if (cc.size() != 0) { - for (auto& c : cc){ - ClusterContainer::AddCluster(mClusterArray, c.getCRU(),c.getRow(),c.getQ(),c.getQmax(), - c.getPadMean(),c.getTimeMean(),c.getPadSigma(),c.getTimeSigma()); + for (auto& c : cc) { + mClusterArray->emplace_back(c); } } } diff --git a/Detectors/TPC/reconstruction/src/TPCCATracking.cxx b/Detectors/TPC/reconstruction/src/TPCCATracking.cxx index 917fba3ee7125..04f88f72868cc 100644 --- a/Detectors/TPC/reconstruction/src/TPCCATracking.cxx +++ b/Detectors/TPC/reconstruction/src/TPCCATracking.cxx @@ -22,7 +22,6 @@ #include "TPCBase/Sector.h" #include "TPCReconstruction/TrackTPC.h" #include "TPCReconstruction/Cluster.h" -#include "TPCReconstruction/HwCluster.h" #include "TPCBase/ParameterDetector.h" #include "TPCBase/ParameterGas.h" #include "TPCBase/ParameterElectronics.h" @@ -64,7 +63,7 @@ void TPCCATracking::deinitialize() { mClusterData = nullptr; } -int TPCCATracking::runTracking(TChain* inputClustersChain, const std::vector* inputClustersArray, std::vector* outputTracks) { +int TPCCATracking::runTracking(TChain* inputClustersChain, const std::vector* inputClustersArray, std::vector* outputTracks) { if (mTrackingCAO2Interface == nullptr) return (1); if (inputClustersChain && inputClustersArray) { LOG(FATAL) << "Internal error, must not pass in both TChain and TClonesArray of clusters\n"; diff --git a/Detectors/TPC/reconstruction/src/TPCReconstructionLinkDef.h b/Detectors/TPC/reconstruction/src/TPCReconstructionLinkDef.h index b1a200f3ebe5d..36c4b0b30a6f5 100644 --- a/Detectors/TPC/reconstruction/src/TPCReconstructionLinkDef.h +++ b/Detectors/TPC/reconstruction/src/TPCReconstructionLinkDef.h @@ -25,14 +25,12 @@ #pragma link C++ class o2::TPC::TPCCATracking; #pragma link C++ class o2::TPC::HardwareClusterDecoder; -#pragma link C++ class o2::TPC::BoxCluster+; #pragma link C++ class o2::TPC::BoxClusterer+; #pragma link C++ class o2::TPC::ClusterTimeStamp+; #pragma link C++ class o2::TPC::Cluster+; #pragma link C++ class o2::TPC::Clusterer+; #pragma link C++ class o2::TPC::ClusterContainer+; #pragma link C++ class o2::TPC::ClustererTask+; -#pragma link C++ class o2::TPC::HwCluster+; #pragma link C++ class o2::TPC::HwClusterer+; #pragma link C++ class o2::TPC::HwClusterFinder+; #pragma link C++ class o2::TPC::HwFixedPoint+; From 16632d828d2fc7a89a4ce8ae7f1311d578a1d005 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Fri, 10 Nov 2017 13:09:26 +0100 Subject: [PATCH 04/12] removed Box-/HwCluster files --- .../include/TPCReconstruction/BoxCluster.h | 93 ---------- .../include/TPCReconstruction/HwCluster.h | 96 ---------- .../TPC/reconstruction/src/BoxCluster.cxx | 58 ------ .../TPC/reconstruction/src/HwCluster.cxx | 170 ------------------ 4 files changed, 417 deletions(-) delete mode 100644 Detectors/TPC/reconstruction/include/TPCReconstruction/BoxCluster.h delete mode 100644 Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h delete mode 100644 Detectors/TPC/reconstruction/src/BoxCluster.cxx delete mode 100644 Detectors/TPC/reconstruction/src/HwCluster.cxx diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxCluster.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxCluster.h deleted file mode 100644 index ee3f3ba959225..0000000000000 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxCluster.h +++ /dev/null @@ -1,93 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 BoxCluster.h -/// \brief Class to have some more info about the BoxClusterer clusters -#ifndef ALICEO2_TPC_BOXCLUSTER_H -#define ALICEO2_TPC_BOXCLUSTER_H - -#ifndef __CINT__ -#include // for base_object -#endif - -#include "TPCReconstruction/Cluster.h" -#include "FairTimeStamp.h" // for FairTimeStamp -#include "Rtypes.h" // for Double_t, ULong_t, etc - -namespace boost { namespace serialization { class access; } } - -namespace o2{ - namespace TPC{ - - /// \class BoxCluster - /// \brief Class to have some more info about the BoxClusterer clusters - /// - class BoxCluster : public Cluster { - public: - - /// Default constructor - BoxCluster(); - - /// Constructor, initializing all values - /// @param cru CRU (sector) - /// @param row Row - /// @param pad Pad with the maximum charge - /// @param time Time bin with the maximum charge - /// @param q Total charge of cluster - /// @param qmax Maximum charge in a single cell (pad, time) - /// @param padmean Mean position of cluster in pad direction - /// @param padsigma Sigma of cluster in pad direction - /// @param timemean Mean position of cluster in time direction - /// @param timesigma Sigma of cluster in time direction - /// @param size Number of pads * 10 + Number of time bins - BoxCluster(Short_t cru, Short_t row, Short_t pad, Short_t time, - Float_t q, Float_t qmax, Float_t padmean, - Float_t padsigma, Float_t timemean, Float_t timesigma, - Short_t size); - - /// Destructor - ~BoxCluster() = default; - - /// Setter for special Box cluster parameters - /// @param pad Pad with the maximum charge - /// @param time Time bin with the maximum charge - /// @param size Number of pads * 10 + Number of time bins - void setBoxParameters(Short_t pad, Short_t time, Short_t size); - - Short_t getPad() const { return mPad; } - Short_t getTime() const { return mTime; } - Short_t getSize() const { return mSize; } - // get number of pads covered by the cluster - Short_t getNpads() const { return Short_t(mSize/10); } - // get number of time bins covered by the cluster - Short_t getNtimebins() const { return mSize%10; } - - - /// Print function: Print basic information to the output stream - /// @param output stream - /// @return The output stream - friend std::ostream& operator<< (std::ostream& out, const BoxCluster &c) { return c.print(out); } - std::ostream& print(std::ostream &output) const; - - private: -#ifndef __CINT__ - friend class boost::serialization::access; -#endif - - Short_t mPad; - Short_t mTime; - Short_t mSize; - - ClassDefNV(BoxCluster, 1); - }; - } -} - -#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h deleted file mode 100644 index 1c42ab1054abe..0000000000000 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h +++ /dev/null @@ -1,96 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 HwCluster.h -/// \brief Class to have some more info about the HwClusterer clusters -/// \author Sebastian Klewin -#ifndef ALICEO2_TPC_HWCLUSTER_H -#define ALICEO2_TPC_HWCLUSTER_H - -#include "TPCReconstruction/Cluster.h" -#include - -namespace boost { namespace serialization { class access; } } - -namespace o2 { -namespace TPC{ - -/// \class HwCluster -/// \brief Class to store HW clusters with cluster data -class HwCluster : public Cluster { - public: - - /// Default Constructors - HwCluster(); - - /// Constructor - /// \param sizeP Cluster size in pad direction - /// \param sizeT Cluster size in time direction - HwCluster(short sizeP, short sizeT); - - /// Constructor - /// \param cru CRU - /// \param row Row - /// \param sizeP Cluster size in pad direction - /// \param sizeT Cluster size in time direction - /// \param clusterData 2D array of size sizeT x sizeP with cluster data - /// \param maxPad Pad with max charge value - /// \param maxTime Timebin with max charge value - HwCluster(short cru, short row, short sizeP, short sizeT, - float** clusterData, short maxPad, short maxTime); - - /// Destructor - ~HwCluster() = default; - - /// Copy Constructor - /// \param other HwCluster to be copied - HwCluster(const HwCluster& other); - -// short getPad() const { return mPad; } -// short getTime() const { return mTime; } -// short getSizeP() const { return mSizeP; } -// short getSizeT() const { return mSizeT; } - - /// Set all cluster data - /// \param cru CRU - /// \param row Row - /// \param sizeP Cluster size in pad direction - /// \param sizeT Cluster size in time direction - /// \param clusterData 2D array of size sizeT x sizeP with cluster data - /// \param maxPad Pad with max charge value - /// \param maxTime Timebin with max charge value - void setClusterData(short cru, short row, short sizeP, short sizeT, - float** clusterData, short maxPad, short maxTime); - - /// Print function - /// \param output Stream to put the HwCluster on - /// \return The output stream - friend std::ostream& operator<< (std::ostream& out, const HwCluster &c) { return c.print(out); } - std::ostream& print(std::ostream &output) const; - std::ostream& PrintDetails(std::ostream &output) const; - - private: - /// Calculates the cluster properties according to locally stored data - void calculateClusterProperties(); - - short mPad; ///< Pad with max charge - short mTime; ///< Timebin with max charge - short mSizeP; ///< Cluster size in pad direction - short mSizeT; ///< Cluster size in time direction - short mSize; ///< Actual size of cluster - - std::vector> mClusterData; ///< CLuster data - - ClassDefNV(HwCluster, 1); - }; -} -} - -#endif diff --git a/Detectors/TPC/reconstruction/src/BoxCluster.cxx b/Detectors/TPC/reconstruction/src/BoxCluster.cxx deleted file mode 100644 index ccec3b96019d9..0000000000000 --- a/Detectors/TPC/reconstruction/src/BoxCluster.cxx +++ /dev/null @@ -1,58 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 BoxCluster.cxx -/// \brief Class to have some more info about the BoxClusterer clusters - -#include "TPCReconstruction/BoxCluster.h" -#include "TPCReconstruction/Cluster.h" - -ClassImp(o2::TPC::BoxCluster) - -using namespace o2::TPC; - -//________________________________________________________________________ -BoxCluster::BoxCluster(): - Cluster(), - mPad(-1), - mTime(-1), - mSize(-1) -{ -} - -//________________________________________________________________________ -BoxCluster::BoxCluster(Short_t cru, Short_t row, Short_t pad, Short_t time, - Float_t q, Float_t qmax, Float_t padmean, - Float_t padsigma, Float_t timemean, - Float_t timesigma, Short_t size): - Cluster(cru, row, q, qmax, padmean, padsigma, timemean, timesigma), - mPad(pad), - mTime(time), - mSize(size) -{ -} - -//________________________________________________________________________ -void BoxCluster::setBoxParameters(Short_t pad, Short_t time, Short_t size) -{ - mPad = pad; - mTime = time; - mSize = size; -} - -//________________________________________________________________________ -std::ostream& BoxCluster::print(std::ostream &output) const -{ - Cluster::print(output); - output << " centered at (pad, time) = " << mPad << ", " << mTime - << " covering " << Int_t(mSize/10) << " pads and " << mSize%10 - << " time bins"; - return output; -} diff --git a/Detectors/TPC/reconstruction/src/HwCluster.cxx b/Detectors/TPC/reconstruction/src/HwCluster.cxx deleted file mode 100644 index 98581d913245e..0000000000000 --- a/Detectors/TPC/reconstruction/src/HwCluster.cxx +++ /dev/null @@ -1,170 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 HwCluster.cxx -/// \brief Class to have some more info about the HwClusterer clusters - -#include "TPCReconstruction/HwCluster.h" -#include "FairLogger.h" - -#include - -ClassImp(o2::TPC::HwCluster) - -using namespace o2::TPC; - -//________________________________________________________________________ -HwCluster::HwCluster() - : HwCluster(5,5) -{} - -//________________________________________________________________________ -HwCluster::HwCluster(short sizeP, short sizeT) - : Cluster() - , mPad(-1) - , mTime(-1) - , mSizeP(sizeP) - , mSizeT(sizeT) - , mSize(0) - , mClusterData(mSizeT, std::vector (mSizeP, 0)) -{ -} - -//________________________________________________________________________ -HwCluster::HwCluster(short cru, short row, short sizeP, short sizeT, - float** clusterData, short maxPad, short maxTime) - : Cluster(cru,row,-1,-1,-1,-1,-1,-1) - , mPad(maxPad) - , mTime(maxTime) - , mSizeP(sizeP) - , mSizeT(sizeT) - , mSize(0) - , mClusterData(mSizeT, std::vector (mSizeP)) -{ - short t,p; - for (t=0; t 0) { - minP = std::min(minP,p); maxP = std::max(maxP,p); - minT = std::min(minT,t); maxT = std::max(maxT,t); - } - - } - } - - mSize = (maxP-minP+1)*10 + (maxT-minT+1); - - if (qTot > 0) { - meanP /= qTot; - meanT /= qTot; - sigmaP /= qTot; - sigmaT /= qTot; - - sigmaP = std::sqrt(sigmaP - (meanP*meanP)); - sigmaT = std::sqrt(sigmaT - (meanT*meanT)); - - meanP += mPad; - meanT += mTime; - } - - setParameters(getCRU(),getRow(),(double)qTot,(double)qMax,(double)meanP,sigmaP,(double)meanT,sigmaT); -} - -//________________________________________________________________________ -std::ostream& HwCluster::print(std::ostream &output) const -{ - Cluster::print(output); - output << " centered at (pad, time) = " << mPad << ", " << mTime - << " covering " << (short)(mSize/10) << " pads and " << mSize%10 - << " time bins"; - return output; -} - -//________________________________________________________________________ -std::ostream& HwCluster::PrintDetails(std::ostream &output) const -{ - Cluster::print(output); - output << " centered at (pad, time) = " << mPad << ", " << mTime - << " covering " << (short)(mSize/10) << " pads and " << mSize%10 - << " time bins" << "\n"; - short t,p; - for (t=0; t Date: Thu, 9 Nov 2017 13:05:31 +0100 Subject: [PATCH 05/12] continued add Cluster MCLabels WIP: Cluster MCLabels, cleaning up WIP: Cluster MCLabels, fix buffering of MC labels of previous events --- Detectors/TPC/reconstruction/CMakeLists.txt | 1 - .../include/TPCReconstruction/BoxClusterer.h | 10 +- .../include/TPCReconstruction/Clusterer.h | 31 +-- .../include/TPCReconstruction/ClustererTask.h | 66 ++++-- .../TPCReconstruction/HwClusterFinder.h | 56 +++-- .../include/TPCReconstruction/HwClusterer.h | 153 +++++++------ .../reconstruction/macro/RawClusterFinder.C | 2 +- .../reconstruction/macro/readClusterMCtruth.C | 43 ++++ .../TPC/reconstruction/src/BoxClusterer.cxx | 33 ++- .../TPC/reconstruction/src/Clusterer.cxx | 22 +- .../TPC/reconstruction/src/ClustererTask.cxx | 38 ++-- .../reconstruction/src/HwClusterFinder.cxx | 100 +++++---- .../TPC/reconstruction/src/HwClusterer.cxx | 202 +++++++++++------- 13 files changed, 453 insertions(+), 304 deletions(-) create mode 100644 Detectors/TPC/reconstruction/macro/readClusterMCtruth.C diff --git a/Detectors/TPC/reconstruction/CMakeLists.txt b/Detectors/TPC/reconstruction/CMakeLists.txt index 1d0986175d38b..27c2808652fec 100644 --- a/Detectors/TPC/reconstruction/CMakeLists.txt +++ b/Detectors/TPC/reconstruction/CMakeLists.txt @@ -6,7 +6,6 @@ link_directories( ${LINK_DIRECTORIES}) set(SRCS src/AdcClockMonitor.cxx - src/BoxClusterer.cxx src/Cluster.cxx src/Clusterer.cxx diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h index 57dfac5c13d2a..c5b0e4f487ae8 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h @@ -37,17 +37,13 @@ namespace o2{ /// Destructor ~BoxClusterer() override; - // Should this really be a public member? - // Maybe better to just call by process - void Init() override; - /// Steer conversion of points to digits /// @param digits Container with TPC digits /// @param mcDigitTruth MC Digit Truth container - /// @param mcClusterTruth MC Cluster Truth container + /// @param eventCount event counter /// @return Container with clusters - void Process(std::vector const & digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; - void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) override; + void Process(std::vector const & digits, MCLabelContainer const* mcDigitTruth, int eventCount) override; + void Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) override; /// Set a pedestal object void setPedestals(CalPad* pedestals) { mPedestals = pedestals; } diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h index 0e602e74b1907..8b23f4dae72d6 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h @@ -17,8 +17,6 @@ #include #include -#include "TPCReconstruction/ClusterContainer.h" - #include "SimulationDataFormat/MCTruthContainer.h" #include "SimulationDataFormat/MCCompLabel.h" @@ -31,12 +29,12 @@ class Digit; /// \brief Base Class for TPC clusterer class Clusterer { protected: - using MCLabel = o2::dataformats::MCTruthContainer; + using MCLabelContainer = o2::dataformats::MCTruthContainer; public: /// Default Constructor - Clusterer(); + Clusterer() : Clusterer(18, 138, 1024, 5, true, true) {}; /// Constructor /// \param rowsMax Max number of rows to process @@ -51,16 +49,13 @@ class Clusterer { /// Destructor virtual ~Clusterer() = default; - /// Initialization function for clusterer - virtual void Init() = 0; - /// Processing all digits /// \param digits Container with TPC digits - /// @param mcDigitTruth MC Digit Truth container - /// @param mcClusterTruth MC Cluster Truth container + /// \param mcDigitTruth MC Digit Truth container + /// \param eventCount event counter /// \return Container with clusters - virtual void Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) = 0; - virtual void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) = 0; + virtual void Process(std::vector const &digits, MCLabelContainer const* mcDigitTruth, int eventCount) = 0; + virtual void Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) = 0; void setRowsMax(int val) { mRowsMax = val; }; void setPadsMax(int val) { mPadsMax = val; }; @@ -85,7 +80,19 @@ class Clusterer { bool mRequirePositiveCharge; ///< If true, require charge > 0 bool mRequireNeighbouringPad; ///< If true, require 2+ pads minimum - }; +}; + +//________________________________________________________________________ +inline Clusterer::Clusterer(int rowsMax, int padsMax, int timeBinsMax, int minQMax, + bool requirePositiveCharge, bool requireNeighbouringPad) + : mRowsMax(rowsMax) + , mPadsMax(padsMax) + , mTimeBinsMax(timeBinsMax) + , mMinQMax(minQMax) + , mRequirePositiveCharge(requirePositiveCharge) + , mRequireNeighbouringPad(requireNeighbouringPad) +{} + } } diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index 14e4d8ee052f6..6206bf1df0c4c 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -18,20 +18,23 @@ #ifndef __ALICEO2__ClustererTask__ #define __ALICEO2__ClustererTask__ -#include #include "FairTask.h" // for FairTask, InitStatus #include "Rtypes.h" // for ClustererTask::Class, ClassDef, etc -#include "TPCReconstruction/Clusterer.h" // for Clusterer -#include "TPCReconstruction/BoxClusterer.h" // for Clusterer -#include "TPCReconstruction/HwClusterer.h" // for Clusterer +#include "TPCReconstruction/Clusterer.h" // for Clusterer +#include "TPCReconstruction/BoxClusterer.h" // for Clusterer +#include "TPCReconstruction/HwClusterer.h" // for Clusterer #include "SimulationDataFormat/MCTruthContainer.h" #include "SimulationDataFormat/MCCompLabel.h" #include +#include namespace o2 { namespace TPC{ class ClustererTask : public FairTask{ + + using MCLabelContainer = o2::dataformats::MCTruthContainer; + public: ClustererTask(); ~ClustererTask() override; @@ -40,6 +43,10 @@ class ClustererTask : public FairTask{ void Exec(Option_t *option) override; enum class ClustererType : int { HW, Box}; + + /// Switch to enable individual clusterer + /// \param type - Clusterer type, HW or Box + /// \param val - Enable set to true or false void setClustererEnable(ClustererType type, bool val) { switch (type) { case ClustererType::HW: mHwClustererEnable = val; break; @@ -47,6 +54,9 @@ class ClustererTask : public FairTask{ }; }; + /// Returns status of Cluster enable + /// \param type - Clusterer type, HW or Box + /// \return Enable status bool isClustererEnable(ClustererType type) const { switch (type) { case ClustererType::HW: return mHwClustererEnable; @@ -54,36 +64,46 @@ class ClustererTask : public FairTask{ }; }; + /// Returns pointer to requested Clusterer type + /// \param type - Clusterer type, HW or Box + /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() Clusterer* getClusterer(ClustererType type) { switch (type) { - case ClustererType::HW: return mHwClusterer; - case ClustererType::Box: return mBoxClusterer; + case ClustererType::HW: return mHwClusterer.get(); + case ClustererType::Box: return mBoxClusterer.get(); }; }; - BoxClusterer* getBoxClusterer() const { return mBoxClusterer; }; - HwClusterer* getHwClusterer() const { return mHwClusterer; }; - // Clusterer *GetClusterer() const { return fClusterer; } + /// Returns pointer to Box Clusterer + /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() + BoxClusterer* getBoxClusterer() const { return mBoxClusterer.get(); }; + + /// Returns pointer to Box Clusterer + /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() + HwClusterer* getHwClusterer() const { return mHwClusterer.get(); }; /// Switch for triggered / continuous readout /// \param isContinuous - false for triggered readout, true for continuous readout void setContinuousReadout(bool isContinuous); private: - bool mBoxClustererEnable; - bool mHwClustererEnable; - bool mIsContinuousReadout; ///< Switch for continuous readout - - BoxClusterer *mBoxClusterer; - HwClusterer *mHwClusterer; - - std::vector const *mDigitsArray; - o2::dataformats::MCTruthContainer const *mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray - // produced data containers - std::vector *mClustersArray; - std::vector *mHwClustersArray; - o2::dataformats::MCTruthContainer mClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mClustersArrays - o2::dataformats::MCTruthContainer mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays + bool mBoxClustererEnable; ///< Switch to enable Box Clusterfinder + bool mHwClustererEnable; ///< Switch to enable Hw Clusterfinder + bool mIsContinuousReadout; ///< Switch for continuous readout + int mEventCount; ///< Event counter + + std::unique_ptr mBoxClusterer; ///< Box Clusterfinder instance + std::unique_ptr mHwClusterer; ///< Hw Clusterfinder instance + + // Digit arrays + std::vector const *mDigitsArray; ///< Array of TPC digits + MCLabelContainer const *mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray + + // Cluster arrays + std::vector *mClustersArray; ///< Array of clusters found by Box Clusterfinder + std::vector *mHwClustersArray; ///< Array of clusters found by Hw Clusterfinder + std::unique_ptr mClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mClustersArrays + std::unique_ptr mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays ClassDefOverride(ClustererTask, 1) }; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index d606c55c8456f..2d04d5ded74b9 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -14,18 +14,32 @@ #ifndef ALICEO2_TPC_HWClusterFinder_H_ #define ALICEO2_TPC_HWClusterFinder_H_ +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" + #include +#include #include namespace o2{ namespace TPC { class Cluster; - + /// \class HwClusterFinder /// \brief Class for TPC HW cluster finder class HwClusterFinder { public: + struct MiniDigit { + float charge; + int event; + int index; + + MiniDigit() : charge(0), event(-1), index(-1) {}; + MiniDigit(const MiniDigit& other) : charge(other.charge), event(other.event), index(other.index) {}; + void clear() { charge = 0; event = -1; index = -1; }; + }; + /// Default Constructor /// \param cru CRU of this cluster finder /// \param row Row of this cluster finder @@ -36,7 +50,7 @@ class HwClusterFinder { /// \param diffThreshold Minimum charge difference at neighboring pads /// \param chargeThreshold Minimum charge of cluster peak /// \param requirePositiveCharge Charge >0 required - HwClusterFinder(short cru, short row, short id, + HwClusterFinder(short cru, short row, short id, short padOffset, short pads=8, short timebins=8, float diffThreshold=0, float chargeThreshold=5, bool requirePositiveCharge=true); @@ -50,14 +64,14 @@ class HwClusterFinder { /// \param timebin Array of size "length" with new charges /// \param globalTime Global time of this timebin /// \param length Size of array "timebin" - bool AddTimebin(float* timebin, unsigned globalTime, int length = 8); + bool AddTimebin(MiniDigit* timebin, unsigned globalTime, int length = 8); - /// Add multiple timebins at once - /// \param nBins Number of timebins - /// \param timebins 2D array with new charges - /// \param globalTime Global time of this timebin - /// \param length Size of array "timebin" - bool AddTimebins(int nBins, float** timebins, unsigned globalTimeOfLast, int length = 8); +// /// Add multiple timebins at once +// /// \param nBins Number of timebins +// /// \param timebins 2D array with new charges +// /// \param globalTime Global time of this timebin +// /// \param length Size of array "timebin" +// bool AddTimebins(int nBins, float** timebins, unsigned globalTimeOfLast, int length = 8); /// Add a timebin with charges of 0 /// \param globalTime Global time of this timebin @@ -71,7 +85,7 @@ class HwClusterFinder { /// \param globalTimeAfterReset Global time of the first timebin after reset void reset(unsigned globalTimeAfterReset); - + // Getter functions int getGlobalTimeOfLast() const { return mGlobalTimeOfLast; } int getTimebinsAfterLastProcessing() const { return mTimebinsAfterLastProcessing; }; @@ -88,10 +102,11 @@ class HwClusterFinder { bool getRequirePositiveCharge() const { return mRequirePositiveCharge; } bool getRequireNeighbouringPad() const { return mRequireNeighbouringPad; } bool getRequireNeighbouringTimebin() const { return mRequireNeighbouringTimebin; } - bool getAutoProcessing() const { return mAutoProcessing; } + bool getAutoProcessing() const { return mAutoProcessing; } bool getmAssignChargeUnique() const { return mAssignChargeUnique; } HwClusterFinder* getNextCF() const { return mNextCF; } std::vector* getClusterContainer() { return &clusterContainer; } + std::vector>>* getClusterDigitIndices() { return &clusterDigitIndices; } // Setter functions void setTimebinsAfterLastProcessing(int val) { mTimebinsAfterLastProcessing = val; }; @@ -114,7 +129,7 @@ class HwClusterFinder { /// Clears the local cluster storage - void clearClusterContainer() { clusterContainer.clear(); } + void clearClusterContainer() { clusterContainer.clear(); clusterDigitIndices.clear(); } /// Process the cluster finding bool findCluster(); @@ -123,19 +138,20 @@ class HwClusterFinder { /// \param time Time bin of found cluster peak /// \param pad Pad of found cluster peak /// \param cluster Cluster charges - void clusterAlreadyUsed(short time, short pad, float** cluster); + void clusterAlreadyUsed(short time, short pad, MiniDigit** cluster); private: - float chargeForCluster(float* charge, float* toCompare); + MiniDigit chargeForCluster(MiniDigit* charge, MiniDigit* toCompare); void printCluster(short time, short pad); // local variables std::vector clusterContainer; + std::vector>> clusterDigitIndices; int mTimebinsAfterLastProcessing; - float** mData; - float** tmpCluster; - float* mZeroTimebin; + MiniDigit** mData; + MiniDigit** tmpCluster; + MiniDigit* mZeroTimebin; // configuration @@ -161,7 +177,7 @@ class HwClusterFinder { }; //________________________________________________________________________ -inline bool HwClusterFinder::AddTimebin(float* timebin, unsigned globalTime, int length) +inline bool HwClusterFinder::AddTimebin(MiniDigit* timebin, unsigned globalTime, int length) { mGlobalTimeOfLast = globalTime; ++mTimebinsAfterLastProcessing; @@ -169,7 +185,7 @@ inline bool HwClusterFinder::AddTimebin(float* timebin, unsigned globalTime, int // // reordering of the local arrays // - float* data0 = mData[0]; + MiniDigit* data0 = mData[0]; std::memmove(mData,mData+1,(mTimebins-1)*sizeof(mData[0])); mData[mTimebins-1] = data0; if (length < mPads) { @@ -187,4 +203,4 @@ inline bool HwClusterFinder::AddTimebin(float* timebin, unsigned globalTime, int } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index efdef3f219c6e..4face4f5e87a3 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -22,7 +22,10 @@ #include "SimulationDataFormat/MCCompLabel.h" #include +#include #include +#include +#include namespace o2{ namespace TPC { @@ -34,113 +37,139 @@ class Digit; /// \class HwClusterer /// \brief Class for TPC HW cluster finding class HwClusterer : public Clusterer { - using MCLabel = o2::dataformats::MCTruthContainer; + + using MCLabelContainer = o2::dataformats::MCTruthContainer; + public: enum class Processing : int { Sequential, Parallel}; /// Constructor - /// \param output is pointer to vector to be filled with clusters + /// \param clusterOutput is pointer to vector to be filled with clusters + /// \param labelOutput is reference to storage to be filled with MC labels /// \param processingType parallel or sequential - /// \param globalTime value of first timebin /// \param cru Number of CRUs to process - /// \param minQDiff Min charge differece + /// \param minQDiff Min charge difference /// \param assignChargeUnique Avoid using same charge for multiple nearby clusters /// \param enableNoiseSim Enables the Noise simulation for empty pads (noise object has to be set) /// \param enablePedestalSubtraction Enables the Pedestal subtraction (pedestal object has to be set) /// \param padsPerCF Pads per cluster finder - /// \param timebinsPerCF Timebins per cluster finder - /// \param cfPerRow Number of cluster finder in each row + /// \param timebinsPerCF Time bins per cluster finder HwClusterer( - std::vector *output, + std::vector *clusterOutput, + std::unique_ptr &labelOutput, Processing processingType = Processing::Parallel, - int globalTime = 0, int cruMin = 0, - int cruMax = 360, + int cruMax = 359, float minQDiff = 0, bool assignChargeUnique = false, bool enableNoiseSim = true, bool enablePedestalSubtraction = true, int padsPerCF = 8, - int timebinsPerCF = 8, - int cfPerRow = 0); + int timebinsPerCF = 8); /// Destructor ~HwClusterer() override; - // Should this really be a public member? - // Maybe better to just call by process - void Init() override; - /// Steer conversion of points to digits /// @param digits Container with TPC digits /// @param mcDigitTruth MC Digit Truth container - /// @param mcClusterTruth MC Cluster Truth container + /// @param eventCount event counter /// @return Container with clusters - void Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) override; - void Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) override; + void Process(std::vector const &digits, + MCLabelContainer const* mcDigitTruth, int eventCount) override; + void Process(std::vector>& digits, + MCLabelContainer const* mcDigitTruth, int eventCount) override; - void setProcessingType(Processing processing) { mProcessingType = processing; }; + /// Switch processing type between parallel and sequential + /// \param type Type to be used + void setProcessingType(Processing type) { mProcessingType = type; }; - void setNoiseObject(CalDet* noiseObject) { mNoiseObject = noiseObject; }; - void setPedestalObject(CalDet* pedestalObject) { mPedestalObject = pedestalObject; }; + /// Setter for noise object, noise will be added before cluster finding + /// \param noiseObject CalDet object, containing noise simulation + void setNoiseObject(std::shared_ptr> noiseObject) { mNoiseObject = noiseObject; }; + + /// Setter for pedestal object, pedestal value will be subtracted before cluster finding + /// \param pedestalObject CalDet object, containing pedestals for each pad + void setPedestalObject(std::shared_ptr> pedestalObject) { mPedestalObject = pedestalObject; }; /// Switch for triggered / continuous readout /// \param isContinuous - false for triggered readout, true for continuous readout void setContinuousReadout(bool isContinuous) { mIsContinuousReadout = isContinuous; }; + /// Setters for CRU range to be processed + /// param cru ID of min/max CRU to be processed void setCRUMin(int cru) { mCRUMin = cru; }; void setCRUMax(int cru) { mCRUMax = cru; }; private: - // To be done - /* BoxClusterer(const BoxClusterer &); */ - /* BoxClusterer &operator=(const BoxClusterer &); */ + /* + * Helper functions + */ + + /// Configuration struct for the processDigits function struct CfConfig { - int iCRU; - int iMaxRows; - int iMaxPads; - int iMinTimeBin; - int iMaxTimeBin; - bool iEnableNoiseSim; - bool iEnablePedestalSubtraction; - bool iIsContinuousReadout; - CalDet* iNoiseObject; - CalDet* iPedestalObject; + int iCRU; ///< CRU ID + int iMaxRows; ///< Maximum row number to be processed + int iMaxPads; ///< Maximum number of pads per row + int iMinTimeBin; ///< Minumum digit time bin + int iMaxTimeBin; ///< Maximum digit time bin + bool iEnableNoiseSim; ///< Noise simulation enable switch + bool iEnablePedestalSubtraction; ///< Pedestal subtraction enable switch + bool iIsContinuousReadout; ///< Continous simulation switch + std::shared_ptr> iNoiseObject; ///< Pointer to noise object + std::shared_ptr> iPedestalObject; ///< Pointer to pedestal object }; + /// Processing the digits, made static to allow for multithreading + /// \param digits Reference to digit container + /// \param clusterFinder Reference to container holding all cluster finder instances + /// \param cluster Reference to container for found clusters + /// \param label Reference to container for MC labels of found clusters + /// \param config Configuration for the cluster finding static void processDigits( - const std::vector>>>& digits, - const std::vector>& clusterFinder, + const std::vector>>& digits, + const std::vector>>& clusterFinder, std::vector& cluster, + std::vector>>& label, CfConfig config); - void ProcessTimeBins(int iTimeBinMin, int iTimeBinMax); - - std::vector>> mClusterFinder; - std::vector>>>> mDigitContainer; - - std::vector> mClusterStorage; - - Processing mProcessingType; - - int mGlobalTime; - int mCRUMin; - int mCRUMax; - float mMinQDiff; - bool mAssignChargeUnique; - bool mEnableNoiseSim; - bool mEnablePedestalSubtraction; - bool mIsContinuousReadout; ///< Switch for continuous readout - int mPadsPerCF; - int mTimebinsPerCF; - int mCfPerRow; - int mLastTimebin; - - std::vector* mClusterArray; ///< Internal cluster storage - - CalDet* mNoiseObject; - CalDet* mPedestalObject; + /// Handling of the parallel cluster finder threads + /// \param iTimeBinMin Minimum time bin to be processed + /// \param iTimeBinMax Maximum time bin to be processed + /// \param mcDigitTruth MC Truth information associated with the digits + /// \param eventCount Event counter + void ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, + MCLabelContainer const* mcDigitTruth, int eventCount); + + /* + * class members + */ + Processing mProcessingType; ///< Processing type for cluster finding + bool mAssignChargeUnique; ///< Setting for CF to use charge only for one cluster + bool mEnableNoiseSim; ///< Switch for noise simulation + bool mEnablePedestalSubtraction; ///< Switch for pedestal subtraction + bool mIsContinuousReadout; ///< Switch for continuous readout + int mCRUMin; ///< Minimum CRU ID to be processed + int mCRUMax; ///< Maximum CRU ID to be processed + int mPadsPerCF; ///< Number of pads per cluster finder instance + int mTimebinsPerCF; ///< Number of time bins per cluster finder instance + int mLastTimebin; ///< Last time bin of previous event + float mMinQDiff; ///< Minimum charge difference between neighboring pads / time bins + + std::vector>>> mClusterFinder; ///< Cluster finder container for each row in each CRU + std::vector>>> mDigitContainer; ///< Sorted digit container for each row in each CRU. Tuple consists of pointer to digit, original digit index and event count + + std::vector> mClusterStorage; ///< Cluster storage for each CRU + std::vector>>> mClusterDigitIndexStorage; ///< Container for digit indices, used in found clusters. Pair consists of original digit index and event count + + std::vector* mClusterArray; ///< Pointer to output cluster storage + std::unique_ptr& mClusterMcLabelArray; ///< Internal MC Label storage + + std::shared_ptr> mNoiseObject; ///< Pointer to the CalDet object for noise simulation + std::shared_ptr> mPedestalObject; ///< Pointer to the CalDet object for the pedestal subtraction + + std::map> mLastMcDigitTruth; }; } } diff --git a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C index f4bb0214c6d45..649ecc100fa2a 100644 --- a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C +++ b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C @@ -176,7 +176,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt // HW cluster finder std::unique_ptr cl; if (clustererType == ClustererType::HW) { - HwClusterer *hwCl = new HwClusterer(&arrCluster, HwClusterer::Processing::Parallel, 0, 0, 4); + HwClusterer *hwCl = new HwClusterer(&arrCluster, HwClusterer::Processing::Parallel, 0, 3); hwCl->setContinuousReadout(false); hwCl->setPedestalObject(pedestal); cl = std::unique_ptr(hwCl); diff --git a/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C new file mode 100644 index 0000000000000..2eb01afa416c6 --- /dev/null +++ b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C @@ -0,0 +1,43 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 readCLusterMCtruth.cxx +/// \brief This macro demonstrates how to extract the MC truth information from +/// the clusters +/// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de +#include + +void readClusterMCtruth(std::string filename) +{ + TFile *clusterFile = TFile::Open(filename.data()); + TTree *clusterTree = (TTree*)clusterFile->Get("o2sim"); + + std::vector* clusters = nullptr; + clusterTree->SetBranchAddress("TPCClusterHW",&clusters); + + o2::dataformats::MCTruthContainer mMCTruthArray; + o2::dataformats::MCTruthContainer *mcTruthArray(&mMCTruthArray); + clusterTree->SetBranchAddress("TPCClusterHWMCTruth", &mcTruthArray); + + for(int iEvent=0; iEventGetEntriesFast(); ++iEvent) { + int cluster = 0; + clusterTree->GetEntry(iEvent); + for(auto& inputcluster : *clusters) { + gsl::span mcArray = mMCTruthArray.getLabels(cluster); + for(int j=0; j(mcArray.size()); ++j) { + std::cout << "Cluster " << cluster << " from Event " << + mMCTruthArray.getElement(mMCTruthArray.getMCTruthHeader(cluster).index+j).getEventID() << " with Track ID " << + mMCTruthArray.getElement(mMCTruthArray.getMCTruthHeader(cluster).index+j).getTrackID() << "\n"; + } + ++cluster; + } + } + +} diff --git a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx index dc5e1b1b58cc2..7f043a38f451a 100644 --- a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx @@ -110,23 +110,6 @@ BoxClusterer::BoxClusterer(std::vector *output): mAllSigBins(nullptr), mAllNSigBins(nullptr), mPedestals(nullptr) -{ -} - -//________________________________________________________________________ -BoxClusterer::~BoxClusterer() -{ - for (Int_t iRow = 0; iRow < mRowsMax; iRow++) { - delete [] mAllBins[iRow]; - delete [] mAllSigBins[iRow]; - } - delete [] mAllBins; - delete [] mAllSigBins; - delete [] mAllNSigBins; -} - -//________________________________________________________________________ -void BoxClusterer::Init() { mAllBins = new Float_t*[mRowsMax]; mAllSigBins = new Int_t*[mRowsMax]; @@ -144,7 +127,19 @@ void BoxClusterer::Init() } //________________________________________________________________________ -void BoxClusterer::Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) +BoxClusterer::~BoxClusterer() +{ + for (Int_t iRow = 0; iRow < mRowsMax; iRow++) { + delete [] mAllBins[iRow]; + delete [] mAllSigBins[iRow]; + } + delete [] mAllBins; + delete [] mAllSigBins; + delete [] mAllNSigBins; +} + +//________________________________________________________________________ +void BoxClusterer::Process(std::vector const &digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); // check this @@ -180,7 +175,7 @@ void BoxClusterer::Process(std::vector const &digits, MCLabel co } //________________________________________________________________________ -void BoxClusterer::Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel& mcClusterTruth) +void BoxClusterer::Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); diff --git a/Detectors/TPC/reconstruction/src/Clusterer.cxx b/Detectors/TPC/reconstruction/src/Clusterer.cxx index 11ea740f4dfe6..2ded5216110c0 100644 --- a/Detectors/TPC/reconstruction/src/Clusterer.cxx +++ b/Detectors/TPC/reconstruction/src/Clusterer.cxx @@ -8,29 +8,9 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. - /// \file Clusterer.cxx +/// \file Clusterer.cxx /// \brief Base class for TPC Clusterer #include "TPCReconstruction/Clusterer.h" - - using namespace o2::TPC; - -Clusterer::Clusterer() - : Clusterer(18, 138, 1024, 5, true, true) -{ -} - -//________________________________________________________________________ -Clusterer::Clusterer(int rowsMax, int padsMax, int timeBinsMax, int minQMax, - bool requirePositiveCharge, bool requireNeighbouringPad) - : mRowsMax(rowsMax) - , mPadsMax(padsMax) - , mTimeBinsMax(timeBinsMax) - , mMinQMax(minQMax) - , mRequirePositiveCharge(requirePositiveCharge) - , mRequireNeighbouringPad(requireNeighbouringPad) -{ -} - diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 07c1d78a575b0..e4540c515f78b 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -17,7 +17,7 @@ // #include "TPCReconstruction/ClustererTask.h" -#include "TPCReconstruction/ClusterContainer.h" // for ClusterContainer +#include "TPCReconstruction/Cluster.h" #include "TPCBase/Digit.h" #include "FairLogger.h" // for LOG #include "FairRootManager.h" // for FairRootManager @@ -32,14 +32,15 @@ ClustererTask::ClustererTask() , mBoxClustererEnable(false) , mHwClustererEnable(false) , mIsContinuousReadout(true) + , mEventCount(0) , mBoxClusterer(nullptr) , mHwClusterer(nullptr) , mDigitsArray(nullptr) , mDigitMCTruthArray(nullptr) , mClustersArray(nullptr) , mHwClustersArray(nullptr) - , mClustersMCTruthArray() - , mHwClustersMCTruthArray() + , mClustersMCTruthArray(nullptr) + , mHwClustersMCTruthArray(nullptr) { } @@ -48,13 +49,11 @@ ClustererTask::~ClustererTask() { LOG(DEBUG) << "Enter Destructor of ClustererTask" << FairLogger::endl; - if (mBoxClustererEnable) delete mBoxClusterer; - if (mHwClustererEnable) delete mHwClusterer; - if (mClustersArray) delete mClustersArray; if (mHwClustersArray) delete mHwClustersArray; + } //_____________________________________________________________________ @@ -79,7 +78,7 @@ InitStatus ClustererTask::Init() mDigitMCTruthArray = mgr->InitObjectAs("TPCDigitMCTruth"); if( !mDigitMCTruthArray ) { - LOG(ERROR) << "TPC MC truth not registered in the FairRootManager. Exiting ..." << FairLogger::endl; + LOG(ERROR) << "TPC MC Truth not registered in the FairRootManager. Exiting ..." << FairLogger::endl; return kERROR; } @@ -89,11 +88,12 @@ InitStatus ClustererTask::Init() mClustersArray = new std::vector; mgr->RegisterAny("TPCCluster", mClustersArray, kTRUE); - mgr->Register("TPCClusterMCTruth", "TPC", &mClustersMCTruthArray, kTRUE); + // Register MC Truth output container + mClustersMCTruthArray = std::make_unique(); + mgr->Register("TPCClusterMCTruth", "TPC", mClustersMCTruthArray.get(), kTRUE); // create clusterer and pass output pointer - mBoxClusterer = new BoxClusterer(mClustersArray); - mBoxClusterer->Init(); + mBoxClusterer = std::make_unique(mClustersArray); } if (mHwClustererEnable) { @@ -101,12 +101,13 @@ InitStatus ClustererTask::Init() mHwClustersArray = new std::vector; mgr->RegisterAny("TPCClusterHW", mHwClustersArray, kTRUE); - mgr->Register("TPCClusterHWMCTruth", "TPC", &mHwClustersMCTruthArray, kTRUE); + // Register MC Truth output container + mHwClustersMCTruthArray = std::make_unique(); + mgr->Register("TPCClusterHWMCTruth", "TPC", mHwClustersMCTruthArray.get(), kTRUE); // create clusterer and pass output pointer - mHwClusterer = new HwClusterer(mHwClustersArray); + mHwClusterer = std::make_unique(mHwClustersArray,mHwClustersMCTruthArray);//,HwClusterer::Processing::Parallel,0,359); mHwClusterer->setContinuousReadout(mIsContinuousReadout); - mHwClusterer->Init(); // TODO: implement noise/pedestal objecta // mHwClusterer->setNoiseObject(); // mHwClusterer->setPedestalObject(); @@ -122,13 +123,16 @@ void ClustererTask::Exec(Option_t *option) if (mBoxClustererEnable) { mClustersArray->clear(); - mBoxClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mClustersMCTruthArray); + if(mClustersMCTruthArray) mClustersMCTruthArray->clear(); + mBoxClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); + LOG(DEBUG) << "Box clusterer found " << mClustersArray->size() << " clusters" << FairLogger::endl; } if (mHwClustererEnable) { - mHwClustersArray->clear(); - mHwClustersMCTruthArray.clear(); - mHwClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mHwClustersMCTruthArray); + if(mHwClustersArray) mHwClustersArray->clear(); + if(mHwClustersMCTruthArray) mHwClustersMCTruthArray->clear(); + mHwClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl; } + ++mEventCount; } diff --git a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx index f89725844671b..8c0f1bf3e2107 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx @@ -66,17 +66,17 @@ HwClusterFinder::HwClusterFinder( } short t,p; - mData = new float*[mTimebins]; + mData = new MiniDigit*[mTimebins]; for (t = 0; t < mTimebins; ++t){ - mData[t] = new float[mPads]; - for (p = 0; p < mPads; ++p){ - mData[t][p] = 0; - } + mData[t] = new MiniDigit[mPads]; +// for (p = 0; p < mPads; ++p){ +// mData[t][p] = MiniDigit(); +// } } - tmpCluster = new float*[mClusterSizeTime]; + tmpCluster = new MiniDigit*[mClusterSizeTime]; for (t=0; t= mData[t][p]) continue; - if (mData[t+1][p ] > mData[t][p]) continue; - if (mData[t ][p-1] >= mData[t][p]) continue; - if (mData[t ][p+1] > mData[t][p]) continue; - if (mData[t-1][p-1] >= mData[t][p]) continue; - if (mData[t+1][p+1] > mData[t][p]) continue; - if (mData[t+1][p-1] > mData[t][p]) continue; - if (mData[t-1][p+1] >= mData[t][p]) continue; + if (mData[t-1][p ].charge >= mData[t][p].charge) continue; + if (mData[t+1][p ].charge > mData[t][p].charge) continue; + if (mData[t ][p-1].charge >= mData[t][p].charge) continue; + if (mData[t ][p+1].charge > mData[t][p].charge) continue; + if (mData[t-1][p-1].charge >= mData[t][p].charge) continue; + if (mData[t+1][p+1].charge > mData[t][p].charge) continue; + if (mData[t+1][p-1].charge > mData[t][p].charge) continue; + if (mData[t-1][p+1].charge >= mData[t][p].charge) continue; // printf("##\n"); // printf("## cluster found at t=%d, p=%d (in CF %d in row %d of CRU %d)\n",t,p,mId,mRow,mCRU); // printf("##\n"); @@ -262,7 +263,7 @@ bool HwClusterFinder::findCluster() // prepare temp storage for (tt=0; tt 0 || tmpCluster[tt][pp].event >= 0) + clusterDigitIndices.back().emplace_back(std::make_pair(tmpCluster[tt][pp].index, tmpCluster[tt][pp].event)); qTot += charge; @@ -360,7 +363,6 @@ bool HwClusterFinder::findCluster() minP = std::min(minP,pp); maxP = std::max(maxP,pp); minT = std::min(minT,tt); maxT = std::max(maxT,tt); } - } } @@ -394,9 +396,11 @@ bool HwClusterFinder::findCluster() // // subtract found cluster from storage // + // TODO: really nexessary?? or just set to 0 for (tt=0; tt<5; ++tt) { for (pp=0; pp<5; ++pp) { - mData[t+(tt-2)][p+(pp-2)] -= tmpCluster[tt][pp]; + //mData[t+(tt-2)][p+(pp-2)].charge -= tmpCluster[tt][pp]; + mData[t+(tt-2)][p+(pp-2)].clear();// = MiniDigit(); } } } @@ -408,16 +412,16 @@ bool HwClusterFinder::findCluster() } //________________________________________________________________________ -float HwClusterFinder::chargeForCluster(float* charge, float* toCompare) +HwClusterFinder::MiniDigit HwClusterFinder::chargeForCluster(MiniDigit* charge, MiniDigit* toCompare) { //printf("%.2f - %.2f = %.f compared to %.2f)?\n",toCompare,*charge,toCompare-*charge,-mDiffThreshold); - if ((mRequirePositiveCharge && (*charge > 0)) & - (*toCompare > mDiffThreshold)) {//mChargeThreshold)) { + if ((mRequirePositiveCharge && (charge->charge > 0)) & + (toCompare->charge > mDiffThreshold)) {//mChargeThreshold)) { //printf("\tyes\n"); return *charge; } else { //printf("\tno\n"); - return 0; + return MiniDigit(); } } @@ -443,7 +447,8 @@ void HwClusterFinder::setNextCF(HwClusterFinder* nextCF) } //________________________________________________________________________ -void HwClusterFinder::clusterAlreadyUsed(short time, short pad, float** cluster) +// TODO: really nexessary?? or just set to 0 +void HwClusterFinder::clusterAlreadyUsed(short time, short pad, MiniDigit** cluster) { short localPad = pad - mPadOffset; @@ -453,7 +458,8 @@ void HwClusterFinder::clusterAlreadyUsed(short time, short pad, float** cluster) for (p=localPad-2; p<=localPad+2; ++p){ if (p < 0 || p >= mPads) continue; - mData[t][p] -= cluster[t-time+2][p-localPad+2]; + //mData[t][p].charge -= cluster[t-time+2][p-localPad+2].charge; + mData[t][p].clear();// = MiniDigit(); } } } @@ -464,7 +470,7 @@ void HwClusterFinder::reset(unsigned globalTimeAfterReset) short t,p; for (t = 0; t < mTimebins; ++t){ for (p = 0; p < mPads; ++p){ - mData[t][p] = 0; + mData[t][p].clear();// = MiniDigit(); } } @@ -478,7 +484,7 @@ void HwClusterFinder::printCluster(short time, short pad) for (t = time-2; t <= time+2; ++t) { printf("%d\t\t",t); for (p = pad-2; p <= pad+2; ++p) { - printf("%.2f\t", mData[t][p]); + printf("%.2f\t", mData[t][p].charge); } printf("\n"); } diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 2956a5ae3fc7f..840019fdd738c 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -21,22 +21,23 @@ #include "FairLogger.h" #include "TMath.h" -#include + #include #include -#include std::mutex g_display_mutex; using namespace o2::TPC; //________________________________________________________________________ -HwClusterer::HwClusterer(std::vector *output, Processing processingType, int globalTime, int cruMin, int cruMax, float minQDiff, - bool assignChargeUnique, bool enableNoiseSim, bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF, int cfPerRow) +HwClusterer::HwClusterer(std::vector *clusterOutput, + std::unique_ptr &labelOutput, + Processing processingType, int cruMin, int cruMax, float minQDiff, + bool assignChargeUnique, bool enableNoiseSim, bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF) : Clusterer() - , mClusterArray(output) + , mClusterArray(clusterOutput) + , mClusterMcLabelArray(labelOutput) , mProcessingType(processingType) - , mGlobalTime(globalTime) , mCRUMin(cruMin) , mCRUMax(cruMax) , mMinQDiff(minQDiff) @@ -46,45 +47,24 @@ HwClusterer::HwClusterer(std::vector *output, Processing proce , mIsContinuousReadout(true) , mPadsPerCF(padsPerCF) , mTimebinsPerCF(timebinsPerCF) - , mCfPerRow(cfPerRow) , mLastTimebin(-1) , mNoiseObject(nullptr) , mPedestalObject(nullptr) + , mLastMcDigitTruth()//new MCLabelContainer) { -} - -//________________________________________________________________________ -HwClusterer::~HwClusterer() -{ - LOG(DEBUG) << "Enter Destructor of HwClusterer" << FairLogger::endl; - - for (auto it : mClusterFinder) { - for (auto itt : it) { - for (auto ittt : itt) { - delete ittt; - } - } - } -} - -//________________________________________________________________________ -void HwClusterer::Init() -{ - LOG(DEBUG) << "Enter Initializer of HwClusterer" << FairLogger::endl; - /* * initalize all cluster finder */ - mCfPerRow = (int)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); - mClusterFinder.resize(mCRUMax); + int iCfPerRow = (int)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); + mClusterFinder.resize(mCRUMax+1); const Mapper& mapper = Mapper::instance(); - for (int iCRU = mCRUMin; iCRU < mCRUMax; iCRU++){ + for (int iCRU = mCRUMin; iCRU <= mCRUMax; iCRU++){ mClusterFinder[iCRU].resize(mapper.getNumberOfRowsPartition(iCRU)); for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); iRow++){ - mClusterFinder[iCRU][iRow].resize(mCfPerRow); - for (int iCF = 0; iCF < mCfPerRow; iCF++){ + mClusterFinder[iCRU][iRow].resize(iCfPerRow); + for (int iCF = 0; iCF < iCfPerRow; iCF++){ int padOffset = iCF*(mPadsPerCF-2-2)-2; - mClusterFinder[iCRU][iRow][iCF] = new HwClusterFinder(iCRU,iRow,iCF,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); + mClusterFinder[iCRU][iRow][iCF] = std::make_unique(iCRU,iRow,iCF,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); mClusterFinder[iCRU][iRow][iCF]->setAssignChargeUnique(mAssignChargeUnique); @@ -94,7 +74,7 @@ void HwClusterer::Init() * already used for a cluster. */ if (iCF != 0) { - mClusterFinder[iCRU][iRow][iCF]->setNextCF(mClusterFinder[iCRU][iRow][iCF-1]); + mClusterFinder[iCRU][iRow][iCF]->setNextCF(mClusterFinder[iCRU][iRow][iCF-1].get()); } } } @@ -105,7 +85,8 @@ void HwClusterer::Init() * vector of HwCluster vectors, one vector for each CRU (possible thread) * to store the clusters found there */ - mClusterStorage.resize(mCRUMax); + mClusterStorage.resize(mCRUMax+1); + mClusterDigitIndexStorage.resize(mCRUMax+1); /* @@ -113,17 +94,26 @@ void HwClusterer::Init() * store there only those digits which are relevant for this particular * CRU (thread) */ - mDigitContainer.resize(mCRUMax); - for (int iCRU = mCRUMin; iCRU < mCRUMax; iCRU++) + mDigitContainer.resize(mCRUMax+1); + for (int iCRU = mCRUMin; iCRU <= mCRUMax; iCRU++) mDigitContainer[iCRU].resize(mapper.getNumberOfRowsPartition(iCRU)); } +//________________________________________________________________________ +HwClusterer::~HwClusterer() +{ + LOG(DEBUG) << "Enter Destructor of HwClusterer" << FairLogger::endl; + +// delete mLastMcDigitTruth; +} + //________________________________________________________________________ void HwClusterer::processDigits( - const std::vector>>>& digits, - const std::vector>& clusterFinder, + const std::vector>>& digits, + const std::vector>>& clusterFinder, std::vector& cluster, + std::vector>>& label, const CfConfig config) { int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; @@ -133,7 +123,7 @@ void HwClusterer::processDigits( // std::cout << "thread " << this_id << " started.\n"; // g_display_mutex.unlock(); - float iAllBins[timeDiff][config.iMaxPads]; + HwClusterFinder::MiniDigit iAllBins[timeDiff][config.iMaxPads]; for (int iRow = 0; iRow < config.iMaxRows; iRow++){ @@ -145,27 +135,31 @@ void HwClusterer::processDigits( if (config.iEnableNoiseSim && config.iNoiseObject != nullptr) { for (t=timeDiff; t--;) { for (p=config.iMaxPads; p--;) { - iAllBins[t][p] = config.iNoiseObject->getValue(CRU(config.iCRU),iRow,p); + iAllBins[t][p].charge = config.iNoiseObject->getValue(CRU(config.iCRU),iRow,p); + iAllBins[t][p].index = -1; + iAllBins[t][p].event = -1; } } } else { - std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, 0.f); + std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, HwClusterFinder::MiniDigit()); } /* * fill in digits */ for (auto& digit : digits[iRow]){ - const Int_t iTime = digit.first->getTimeStamp(); - const Int_t iPad = digit.first->getPad() + 2; // offset to have 2 empty pads on the "left side" - const Float_t charge = digit.first->getChargeFloat(); + const Int_t iTime = std::get<0>(digit)->getTimeStamp(); + const Int_t iPad = std::get<0>(digit)->getPad() + 2; // offset to have 2 empty pads on the "left side" + const Float_t charge = std::get<0>(digit)->getChargeFloat(); // std::cout << iCRU << " " << iRow << " " << iPad << " " << iTime << " (" << iTime-minTime << "," << timeDiff << ") " << charge << std::endl; - iAllBins[iTime-config.iMinTimeBin][iPad] += charge; + iAllBins[iTime-config.iMinTimeBin][iPad].charge += charge; + iAllBins[iTime-config.iMinTimeBin][iPad].index = std::get<1>(digit); + iAllBins[iTime-config.iMinTimeBin][iPad].event = std::get<2>(digit); if (config.iEnablePedestalSubtraction && config.iPedestalObject != nullptr) { const float pedestal = config.iPedestalObject->getValue(CRU(config.iCRU),iRow,iPad-2); //printf("digit: %.2f, pedestal: %.2f\n", iAllBins[iTime-config.iMinTimeBin][iPad], pedestal); - iAllBins[iTime-config.iMinTimeBin][iPad] -= pedestal; + iAllBins[iTime-config.iMinTimeBin][iPad].charge -= pedestal; } } @@ -174,7 +168,7 @@ void HwClusterer::processDigits( */ const Short_t iPadsPerCF = clusterFinder[iRow][0]->getNpads(); const Short_t iTimebinsPerCF = clusterFinder[iRow][0]->getNtimebins(); - std::vector::const_reverse_iterator> cfWithCluster; + std::vector>::const_reverse_iterator> cfWithCluster; unsigned time,pad; for (time = 0; time < timeDiff; ++time){ // ordering important!! for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { @@ -229,12 +223,16 @@ void HwClusterer::processDigits( /* * collect found cluster */ - for (std::vector::const_reverse_iterator &cf_it : cfWithCluster) { - std::vector* cc = (*cf_it)->getClusterContainer(); - for (Cluster& c : *cc){ - cluster.push_back(c); + for (auto &cf_rit : cfWithCluster) { + auto cc = (*cf_rit)->getClusterContainer(); + for (auto& c : *cc) cluster.push_back(c); + + auto ll = (*cf_rit)->getClusterDigitIndices(); + for (auto& l : *ll) { + label.push_back(l); } - (*cf_it)->clearClusterContainer(); + + (*cf_rit)->clearClusterContainer(); } } @@ -245,15 +243,17 @@ void HwClusterer::processDigits( } //________________________________________________________________________ -void HwClusterer::Process(std::vector const &digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth - ) +void HwClusterer::Process(std::vector const &digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); + mClusterMcLabelArray->clear(); + /* * clear old storages */ for (auto& cs : mClusterStorage) cs.clear(); + for (auto& cdis : mClusterDigitIndexStorage) cdis.clear(); for (auto& dc : mDigitContainer ) { for (auto& dcc : dc) dcc.clear(); } @@ -263,6 +263,9 @@ void HwClusterer::Process(std::vector const &digits, MCLabel con //int iTimeBinMin = mLastTimebin + 1; int iTimeBinMax = mLastTimebin; + if (mLastMcDigitTruth.find(eventCount) == mLastMcDigitTruth.end()) + mLastMcDigitTruth[eventCount] = std::make_unique(); + gsl::span mcArray; /* * Loop over digits */ @@ -274,26 +277,46 @@ void HwClusterer::Process(std::vector const &digits, MCLabel con */ iTimeBin = digit.getTimeStamp(); - if (iTimeBin < iTimeBinMin) continue; + if (digit.getCRU() < mCRUMin || digit.getCRU() > mCRUMax) { + LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit.getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; + if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); + ++digitIndex; + continue; + } + if (iTimeBin < iTimeBinMin) { + LOG(DEBUG) << "Digit [" << digitIndex << "] time stamp too small (" << iTimeBin << " < " << iTimeBinMin << ")" << FairLogger::endl; + if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); + ++digitIndex; + continue; + } + iTimeBinMax = std::max(iTimeBinMax,iTimeBin); if (mcDigitTruth == nullptr) - mDigitContainer[digit.getCRU()][digit.getRow()].push_back(std::make_pair(&digit,nullptr)); - else - mDigitContainer[digit.getCRU()][digit.getRow()].push_back(std::make_pair(&digit,mcDigitTruth->getLabels(digitIndex))); + mDigitContainer[digit.getCRU()][digit.getRow()].emplace_back(std::make_tuple(&digit,-1,eventCount)); + else { + mDigitContainer[digit.getCRU()][digit.getRow()].emplace_back(std::make_tuple(&digit,digitIndex,eventCount)); + mcArray = mcDigitTruth->getLabels(digitIndex); + for (auto &l : mcArray) mLastMcDigitTruth[eventCount]->addElement(digitIndex,l); + } ++digitIndex; } - ProcessTimeBins(iTimeBinMin, iTimeBinMax); + + ProcessTimeBins(iTimeBinMin, iTimeBinMax, mcDigitTruth, eventCount); + + mLastMcDigitTruth.erase(eventCount-mTimebinsPerCF); } //________________________________________________________________________ -void HwClusterer::Process(std::vector>& digits, MCLabel const* mcDigitTruth, MCLabel &mcClusterTruth) +void HwClusterer::Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); + mClusterMcLabelArray->clear(); /* * clear old storages */ for (auto& cs : mClusterStorage) cs.clear(); + for (auto& cdis : mClusterDigitIndexStorage) cdis.clear(); for (auto& dc : mDigitContainer ) { for (auto& dcc : dc) dcc.clear(); } @@ -302,6 +325,10 @@ void HwClusterer::Process(std::vector>& digits, MCLabel c int iTimeBinMin = (mIsContinuousReadout)?mLastTimebin + 1 : 0; int iTimeBinMax = mLastTimebin; + if (mLastMcDigitTruth.find(eventCount) == mLastMcDigitTruth.end()) + mLastMcDigitTruth[eventCount] = std::make_unique(); + gsl::span mcArray; + /* * Loop over digits */ @@ -313,19 +340,37 @@ void HwClusterer::Process(std::vector>& digits, MCLabel c * add current digit to storage */ iTimeBin = digit->getTimeStamp(); - if (iTimeBin < iTimeBinMin) continue; + if (digit->getCRU() < mCRUMin || digit->getCRU() > mCRUMax) { + LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit->getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; + if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); + ++digitIndex; + continue; + } + if (iTimeBin < iTimeBinMin) { + LOG(DEBUG) << "Digit [" << digitIndex << "] time stamp too small (" << iTimeBin << " < " << iTimeBinMin << ")" << FairLogger::endl; + if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); + ++digitIndex; + continue; + } + iTimeBinMax = std::max(iTimeBinMax,iTimeBin); if (mcDigitTruth == nullptr) - mDigitContainer[digit->getCRU()][digit->getRow()].push_back(std::make_pair(digit,nullptr)); - else - mDigitContainer[digit->getCRU()][digit->getRow()].push_back(std::make_pair(digit,mcDigitTruth->getLabels(digitIndex))); + mDigitContainer[digit->getCRU()][digit->getRow()].emplace_back(std::make_tuple(digit,-1,eventCount)); + else { + mDigitContainer[digit->getCRU()][digit->getRow()].emplace_back(std::make_tuple(digit,digitIndex,eventCount)); + mcArray = mcDigitTruth->getLabels(digitIndex); + for (auto &l : mcArray) mLastMcDigitTruth[eventCount]->addElement(digitIndex,l); + } ++digitIndex; } - ProcessTimeBins(iTimeBinMin, iTimeBinMax); + ProcessTimeBins(iTimeBinMin, iTimeBinMax, mcDigitTruth, eventCount); + + mLastMcDigitTruth.erase(eventCount-mTimebinsPerCF); + } -void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) +void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelContainer const* mcDigitTruth, int eventCount) { /* @@ -343,7 +388,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) */ const Mapper& mapper = Mapper::instance(); - for (int iCRU = mCRUMin; iCRU < mCRUMax; ++iCRU) { + for (int iCRU = mCRUMin; iCRU <= mCRUMax; ++iCRU) { struct CfConfig cfConfig = { iCRU, mapper.getNumberOfRowsPartition(iCRU), @@ -356,13 +401,14 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) mNoiseObject, mPedestalObject }; -// std::cout << "starting CRU " << iCRU << std::endl; +// LOG(DEBUG) << "Start processing CRU " << iCRU << FairLogger::endl; if (mProcessingType == Processing::Parallel) thread_vector.emplace_back( processDigits, // function name std::ref(mDigitContainer[iCRU]), // digit container for individual CRUs std::ref(mClusterFinder[iCRU]), // cluster finder for individual CRUs std::ref(mClusterStorage[iCRU]), // container to store found clusters + std::ref(mClusterDigitIndexStorage[iCRU]), // container to store found cluster MC Labels cfConfig ); else { @@ -370,6 +416,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) std::ref(mDigitContainer[iCRU]), std::ref(mClusterFinder[iCRU]), std::ref(mClusterStorage[iCRU]), + std::ref(mClusterDigitIndexStorage[iCRU]), cfConfig); } } @@ -385,10 +432,17 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) /* * collect clusters from individual cluster finder */ - for (auto& cc : mClusterStorage) { - if (cc.size() != 0) { - for (auto& c : cc) { - mClusterArray->emplace_back(c); + for (int cru = 0; cru < mClusterStorage.size(); ++cru) { + std::vector* clustersFromCRU = &mClusterStorage[cru]; + std::vector>>* labelsFromCRU = &mClusterDigitIndexStorage[cru]; + + for(int c = 0; c < clustersFromCRU->size(); ++c) { + const auto clusterPos = mClusterArray->size(); + mClusterArray->emplace_back(clustersFromCRU->at(c)); + for (auto &digitIndex : labelsFromCRU->at(c)) { + if (digitIndex.first < 0) continue; + for (auto &l : mLastMcDigitTruth[digitIndex.second]->getLabels(digitIndex.first)) + mClusterMcLabelArray->addElement(clusterPos,l); } } } From fddec3a9ba975da17c0dcdf5a9a080de600087e5 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Wed, 15 Nov 2017 11:30:24 +0100 Subject: [PATCH 06/12] more efficient thread spawning + num threads configurable at runtime --- .../include/TPCReconstruction/HwClusterer.h | 39 ++- .../reconstruction/macro/RawClusterFinder.C | 2 +- .../TPC/reconstruction/src/ClustererTask.cxx | 2 +- .../TPC/reconstruction/src/HwClusterer.cxx | 261 +++++++++--------- macro/run_clus_tpc.C | 5 +- run/runTPC.cxx | 8 +- 6 files changed, 159 insertions(+), 158 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index 4face4f5e87a3..101409de86507 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -41,12 +41,10 @@ class HwClusterer : public Clusterer { using MCLabelContainer = o2::dataformats::MCTruthContainer; public: - enum class Processing : int { Sequential, Parallel}; /// Constructor /// \param clusterOutput is pointer to vector to be filled with clusters /// \param labelOutput is reference to storage to be filled with MC labels - /// \param processingType parallel or sequential /// \param cru Number of CRUs to process /// \param minQDiff Min charge difference /// \param assignChargeUnique Avoid using same charge for multiple nearby clusters @@ -57,7 +55,6 @@ class HwClusterer : public Clusterer { HwClusterer( std::vector *clusterOutput, std::unique_ptr &labelOutput, - Processing processingType = Processing::Parallel, int cruMin = 0, int cruMax = 359, float minQDiff = 0, @@ -80,10 +77,6 @@ class HwClusterer : public Clusterer { void Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) override; - /// Switch processing type between parallel and sequential - /// \param type Type to be used - void setProcessingType(Processing type) { mProcessingType = type; }; - /// Setter for noise object, noise will be added before cluster finding /// \param noiseObject CalDet object, containing noise simulation void setNoiseObject(std::shared_ptr> noiseObject) { mNoiseObject = noiseObject; }; @@ -97,10 +90,14 @@ class HwClusterer : public Clusterer { void setContinuousReadout(bool isContinuous) { mIsContinuousReadout = isContinuous; }; /// Setters for CRU range to be processed - /// param cru ID of min/max CRU to be processed + /// \param cru ID of min/max CRU to be processed void setCRUMin(int cru) { mCRUMin = cru; }; void setCRUMax(int cru) { mCRUMax = cru; }; + /// Set number of parallel threads + /// \param threads Number to be set, if 0 hardware default value is used + void setNumThreads(unsigned threads); + private: /* @@ -109,8 +106,10 @@ class HwClusterer : public Clusterer { /// Configuration struct for the processDigits function struct CfConfig { - int iCRU; ///< CRU ID - int iMaxRows; ///< Maximum row number to be processed + unsigned iThreadID; ///< Index of thread + unsigned iThreadMax; ///< Total number of started threads + unsigned iCRUMin; ///< Minimum CRU number to process + unsigned iCRUMax; ///< Maximum CRU number to process int iMaxPads; ///< Maximum number of pads per row int iMinTimeBin; ///< Minumum digit time bin int iMaxTimeBin; ///< Maximum digit time bin @@ -128,10 +127,10 @@ class HwClusterer : public Clusterer { /// \param label Reference to container for MC labels of found clusters /// \param config Configuration for the cluster finding static void processDigits( - const std::vector>>& digits, - const std::vector>>& clusterFinder, - std::vector& cluster, - std::vector>>& label, + const std::vector>>>& digits, + const std::vector>>>& clusterFinder, + std::vector>& cluster, + std::vector>>>& label, CfConfig config); /// Handling of the parallel cluster finder threads @@ -145,16 +144,16 @@ class HwClusterer : public Clusterer { /* * class members */ - Processing mProcessingType; ///< Processing type for cluster finding bool mAssignChargeUnique; ///< Setting for CF to use charge only for one cluster bool mEnableNoiseSim; ///< Switch for noise simulation bool mEnablePedestalSubtraction; ///< Switch for pedestal subtraction bool mIsContinuousReadout; ///< Switch for continuous readout - int mCRUMin; ///< Minimum CRU ID to be processed - int mCRUMax; ///< Maximum CRU ID to be processed - int mPadsPerCF; ///< Number of pads per cluster finder instance - int mTimebinsPerCF; ///< Number of time bins per cluster finder instance int mLastTimebin; ///< Last time bin of previous event + unsigned mCRUMin; ///< Minimum CRU ID to be processed + unsigned mCRUMax; ///< Maximum CRU ID to be processed + unsigned mPadsPerCF; ///< Number of pads per cluster finder instance + unsigned mTimebinsPerCF; ///< Number of time bins per cluster finder instance + unsigned mNumThreads; ///< Number of parallel processing threads float mMinQDiff; ///< Minimum charge difference between neighboring pads / time bins std::vector>>> mClusterFinder; ///< Cluster finder container for each row in each CRU @@ -169,7 +168,7 @@ class HwClusterer : public Clusterer { std::shared_ptr> mNoiseObject; ///< Pointer to the CalDet object for noise simulation std::shared_ptr> mPedestalObject; ///< Pointer to the CalDet object for the pedestal subtraction - std::map> mLastMcDigitTruth; + std::map> mLastMcDigitTruth; ///< Buffer for digit MC truth information }; } } diff --git a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C index 649ecc100fa2a..be85f5c36b677 100644 --- a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C +++ b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C @@ -176,7 +176,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt // HW cluster finder std::unique_ptr cl; if (clustererType == ClustererType::HW) { - HwClusterer *hwCl = new HwClusterer(&arrCluster, HwClusterer::Processing::Parallel, 0, 3); + HwClusterer *hwCl = new HwClusterer(&arrCluster, 0, 3); hwCl->setContinuousReadout(false); hwCl->setPedestalObject(pedestal); cl = std::unique_ptr(hwCl); diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index e4540c515f78b..5f499b8154b74 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -106,7 +106,7 @@ InitStatus ClustererTask::Init() mgr->Register("TPCClusterHWMCTruth", "TPC", mHwClustersMCTruthArray.get(), kTRUE); // create clusterer and pass output pointer - mHwClusterer = std::make_unique(mHwClustersArray,mHwClustersMCTruthArray);//,HwClusterer::Processing::Parallel,0,359); + mHwClusterer = std::make_unique(mHwClustersArray,mHwClustersMCTruthArray);//,0,359); mHwClusterer->setContinuousReadout(mIsContinuousReadout); // TODO: implement noise/pedestal objecta // mHwClusterer->setNoiseObject(); diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 840019fdd738c..02128cb88d5ad 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -31,29 +31,33 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterer::HwClusterer(std::vector *clusterOutput, - std::unique_ptr &labelOutput, - Processing processingType, int cruMin, int cruMax, float minQDiff, - bool assignChargeUnique, bool enableNoiseSim, bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF) + std::unique_ptr &labelOutput, int cruMin, int cruMax, + float minQDiff, bool assignChargeUnique, bool enableNoiseSim, + bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF) : Clusterer() - , mClusterArray(clusterOutput) - , mClusterMcLabelArray(labelOutput) - , mProcessingType(processingType) - , mCRUMin(cruMin) - , mCRUMax(cruMax) - , mMinQDiff(minQDiff) , mAssignChargeUnique(assignChargeUnique) , mEnableNoiseSim(enableNoiseSim) , mEnablePedestalSubtraction(enablePedestalSubtraction) , mIsContinuousReadout(true) + , mLastTimebin(-1) + , mCRUMin(cruMin) + , mCRUMax(cruMax) , mPadsPerCF(padsPerCF) , mTimebinsPerCF(timebinsPerCF) - , mLastTimebin(-1) + , mNumThreads(std::thread::hardware_concurrency()) + , mMinQDiff(minQDiff) + , mClusterFinder() + , mDigitContainer() + , mClusterStorage() + , mClusterDigitIndexStorage() + , mClusterArray(clusterOutput) + , mClusterMcLabelArray(labelOutput) , mNoiseObject(nullptr) , mPedestalObject(nullptr) - , mLastMcDigitTruth()//new MCLabelContainer) + , mLastMcDigitTruth() { /* - * initalize all cluster finder + * initialize all cluster finder */ int iCfPerRow = (int)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); mClusterFinder.resize(mCRUMax+1); @@ -69,7 +73,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, /* - * Connect always two CFs to be able to comunicate found clusters. So + * Connect always two CFs to be able to communicate found clusters. So * the "right" one can tell the one "on the left" which pads were * already used for a cluster. */ @@ -110,131 +114,136 @@ HwClusterer::~HwClusterer() //________________________________________________________________________ void HwClusterer::processDigits( - const std::vector>>& digits, - const std::vector>>& clusterFinder, - std::vector& cluster, - std::vector>>& label, + const std::vector>>>& digits, + const std::vector>>>& clusterFinder, + std::vector>& cluster, + std::vector>>>& label, const CfConfig config) { - int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; - if (timeDiff < 0) return; // std::thread::id this_id = std::this_thread::get_id(); // g_display_mutex.lock(); // std::cout << "thread " << this_id << " started.\n"; // g_display_mutex.unlock(); + int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; + if (timeDiff < 0) return; + const Mapper& mapper = Mapper::instance(); HwClusterFinder::MiniDigit iAllBins[timeDiff][config.iMaxPads]; - for (int iRow = 0; iRow < config.iMaxRows; iRow++){ + for (int iCRU = config.iCRUMin; iCRU <= config.iCRUMax; ++iCRU) { + if (iCRU % config.iThreadMax != config.iThreadID) continue; - /* - * prepare local storage - */ - short t,p; - float noise; - if (config.iEnableNoiseSim && config.iNoiseObject != nullptr) { - for (t=timeDiff; t--;) { - for (p=config.iMaxPads; p--;) { - iAllBins[t][p].charge = config.iNoiseObject->getValue(CRU(config.iCRU),iRow,p); - iAllBins[t][p].index = -1; - iAllBins[t][p].event = -1; - } - } - } else { - std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, HwClusterFinder::MiniDigit()); - } - - /* - * fill in digits - */ - for (auto& digit : digits[iRow]){ - const Int_t iTime = std::get<0>(digit)->getTimeStamp(); - const Int_t iPad = std::get<0>(digit)->getPad() + 2; // offset to have 2 empty pads on the "left side" - const Float_t charge = std::get<0>(digit)->getChargeFloat(); - -// std::cout << iCRU << " " << iRow << " " << iPad << " " << iTime << " (" << iTime-minTime << "," << timeDiff << ") " << charge << std::endl; - iAllBins[iTime-config.iMinTimeBin][iPad].charge += charge; - iAllBins[iTime-config.iMinTimeBin][iPad].index = std::get<1>(digit); - iAllBins[iTime-config.iMinTimeBin][iPad].event = std::get<2>(digit); - if (config.iEnablePedestalSubtraction && config.iPedestalObject != nullptr) { - const float pedestal = config.iPedestalObject->getValue(CRU(config.iCRU),iRow,iPad-2); - //printf("digit: %.2f, pedestal: %.2f\n", iAllBins[iTime-config.iMinTimeBin][iPad], pedestal); - iAllBins[iTime-config.iMinTimeBin][iPad].charge -= pedestal; - } - } + for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); iRow++){ - /* - * copy data to cluster finders - */ - const Short_t iPadsPerCF = clusterFinder[iRow][0]->getNpads(); - const Short_t iTimebinsPerCF = clusterFinder[iRow][0]->getNtimebins(); - std::vector>::const_reverse_iterator> cfWithCluster; - unsigned time,pad; - for (time = 0; time < timeDiff; ++time){ // ordering important!! - for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { - const Short_t cf = pad / (iPadsPerCF-2-2); - clusterFinder[iRow][cf]->AddTimebin(&iAllBins[time][pad],time+config.iMinTimeBin,(config.iMaxPads-pad)>=iPadsPerCF?iPadsPerCF:(config.iMaxPads-pad)); + /* + * prepare local storage + */ + short t,p; + float noise; + if (config.iEnableNoiseSim && config.iNoiseObject != nullptr) { + for (t=timeDiff; t--;) { + for (p=config.iMaxPads; p--;) { + iAllBins[t][p].charge = config.iNoiseObject->getValue(CRU(iCRU),iRow,p); + iAllBins[t][p].index = -1; + iAllBins[t][p].event = -1; + } + } + } else { + std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, HwClusterFinder::MiniDigit()); } /* - * search for clusters and store reference to CF if one was found + * fill in digits */ - if (clusterFinder[iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { - /* - * ordering is important: from right to left, so that the CFs could inform each other if cluster was found - */ - for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { - if ((*rit)->findCluster()) { - cfWithCluster.push_back(rit); - } + for (auto& digit : digits[iCRU][iRow]){ + const Int_t iTime = std::get<0>(digit)->getTimeStamp(); + const Int_t iPad = std::get<0>(digit)->getPad() + 2; // offset to have 2 empty pads on the "left side" + const Float_t charge = std::get<0>(digit)->getChargeFloat(); + + // std::cout << iCRU << " " << iRow << " " << iPad << " " << iTime << " (" << iTime-minTime << "," << timeDiff << ") " << charge << std::endl; + iAllBins[iTime-config.iMinTimeBin][iPad].charge += charge; + iAllBins[iTime-config.iMinTimeBin][iPad].index = std::get<1>(digit); + iAllBins[iTime-config.iMinTimeBin][iPad].event = std::get<2>(digit); + if (config.iEnablePedestalSubtraction && config.iPedestalObject != nullptr) { + const float pedestal = config.iPedestalObject->getValue(CRU(iCRU),iRow,iPad-2); + //printf("digit: %.2f, pedestal: %.2f\n", iAllBins[iTime-config.iMinTimeBin][iPad], pedestal); + iAllBins[iTime-config.iMinTimeBin][iPad].charge -= pedestal; } } - } - /* - * add empty timebins to find last clusters - */ - if (!config.iIsContinuousReadout) { - // +2 so that for sure all data is processed - for (time = 0; time < clusterFinder[iRow][0]->getNtimebins()+2; ++time){ - for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { - (*rit)->AddZeroTimebin(time+timeDiff+config.iMinTimeBin,iPadsPerCF); + /* + * copy data to cluster finders + */ + const Short_t iPadsPerCF = clusterFinder[iCRU][iRow][0]->getNpads(); + const Short_t iTimebinsPerCF = clusterFinder[iCRU][iRow][0]->getNtimebins(); + std::vector>::const_reverse_iterator> cfWithCluster; + unsigned time,pad; + for (time = 0; time < timeDiff; ++time){ // ordering important!! + for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { + const Short_t cf = pad / (iPadsPerCF-2-2); + clusterFinder[iCRU][iRow][cf]->AddTimebin(&iAllBins[time][pad],time+config.iMinTimeBin,(config.iMaxPads-pad)>=iPadsPerCF?iPadsPerCF:(config.iMaxPads-pad)); } /* * search for clusters and store reference to CF if one was found */ - if (clusterFinder[iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { + if (clusterFinder[iCRU][iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { /* * ordering is important: from right to left, so that the CFs could inform each other if cluster was found */ - for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { + for (auto rit = clusterFinder[iCRU][iRow].crbegin(); rit != clusterFinder[iCRU][iRow].crend(); ++rit) { if ((*rit)->findCluster()) { cfWithCluster.push_back(rit); } } } } - for (auto rit = clusterFinder[iRow].crbegin(); rit != clusterFinder[iRow].crend(); ++rit) { - (*rit)->setTimebinsAfterLastProcessing(0); + + /* + * add empty timebins to find last clusters + */ + if (!config.iIsContinuousReadout) { + // +2 so that for sure all data is processed + for (time = 0; time < clusterFinder[iCRU][iRow][0]->getNtimebins()+2; ++time){ + for (auto rit = clusterFinder[iCRU][iRow].crbegin(); rit != clusterFinder[iCRU][iRow].crend(); ++rit) { + (*rit)->AddZeroTimebin(time+timeDiff+config.iMinTimeBin,iPadsPerCF); + } + + /* + * search for clusters and store reference to CF if one was found + */ + if (clusterFinder[iCRU][iRow][0]->getTimebinsAfterLastProcessing() == iTimebinsPerCF-2 -2) { + /* + * ordering is important: from right to left, so that the CFs could inform each other if cluster was found + */ + for (auto rit = clusterFinder[iCRU][iRow].crbegin(); rit != clusterFinder[iCRU][iRow].crend(); ++rit) { + if ((*rit)->findCluster()) { + cfWithCluster.push_back(rit); + } + } + } + } + for (auto rit = clusterFinder[iCRU][iRow].crbegin(); rit != clusterFinder[iCRU][iRow].crend(); ++rit) { + (*rit)->setTimebinsAfterLastProcessing(0); + } } - } - /* - * collect found cluster - */ - for (auto &cf_rit : cfWithCluster) { - auto cc = (*cf_rit)->getClusterContainer(); - for (auto& c : *cc) cluster.push_back(c); + /* + * collect found cluster + */ + for (auto &cf_rit : cfWithCluster) { + auto cc = (*cf_rit)->getClusterContainer(); + for (auto& c : *cc) cluster[iCRU].push_back(c); + + auto ll = (*cf_rit)->getClusterDigitIndices(); + for (auto& l : *ll) { + label[iCRU].push_back(l); + } - auto ll = (*cf_rit)->getClusterDigitIndices(); - for (auto& l : *ll) { - label.push_back(l); + (*cf_rit)->clearClusterContainer(); } - (*cf_rit)->clearClusterContainer(); } - } // g_display_mutex.lock(); @@ -279,12 +288,14 @@ void HwClusterer::Process(std::vector const &digits, MCLabelCont iTimeBin = digit.getTimeStamp(); if (digit.getCRU() < mCRUMin || digit.getCRU() > mCRUMax) { LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit.getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; + // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); ++digitIndex; continue; } if (iTimeBin < iTimeBinMin) { LOG(DEBUG) << "Digit [" << digitIndex << "] time stamp too small (" << iTimeBin << " < " << iTimeBinMin << ")" << FairLogger::endl; + // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); ++digitIndex; continue; @@ -342,12 +353,14 @@ void HwClusterer::Process(std::vector>& digits, MCLabelCo iTimeBin = digit->getTimeStamp(); if (digit->getCRU() < mCRUMin || digit->getCRU() > mCRUMax) { LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit->getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; + // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); ++digitIndex; continue; } if (iTimeBin < iTimeBinMin) { LOG(DEBUG) << "Digit [" << digitIndex << "] time stamp too small (" << iTimeBin << " < " << iTimeBinMin << ")" << FairLogger::endl; + // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); ++digitIndex; continue; @@ -374,24 +387,18 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta { /* - * vector to store all threads for parallel processing - * one thread per CRU (360 in total) + * vector to store threads for parallel processing */ std::vector thread_vector; - if (mProcessingType == Processing::Parallel) - LOG(DEBUG) << std::thread::hardware_concurrency() << " concurrent threads are supported." << FairLogger::endl; - /* - * if CRU number of current digit changes, start processing (either - * sequential or parallel) of all CRUs in between the last processed - * one and the current one. - */ + LOG(DEBUG) << "Starting " << mNumThreads << " threads, hardware supports " << std::thread::hardware_concurrency() << " parallel threads." << FairLogger::endl; - const Mapper& mapper = Mapper::instance(); - for (int iCRU = mCRUMin; iCRU <= mCRUMax; ++iCRU) { + for (unsigned threadId = 0; threadId < std::min(mNumThreads,mCRUMax); ++threadId) { struct CfConfig cfConfig = { - iCRU, - mapper.getNumberOfRowsPartition(iCRU), + threadId, + mNumThreads, + mCRUMin, + mCRUMax, mPadsMax+2+2, iTimeBinMin, iTimeBinMax, @@ -401,24 +408,14 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta mNoiseObject, mPedestalObject }; -// LOG(DEBUG) << "Start processing CRU " << iCRU << FairLogger::endl; - if (mProcessingType == Processing::Parallel) - thread_vector.emplace_back( - processDigits, // function name - std::ref(mDigitContainer[iCRU]), // digit container for individual CRUs - std::ref(mClusterFinder[iCRU]), // cluster finder for individual CRUs - std::ref(mClusterStorage[iCRU]), // container to store found clusters - std::ref(mClusterDigitIndexStorage[iCRU]), // container to store found cluster MC Labels - cfConfig + thread_vector.emplace_back( + processDigits, // function name + std::ref(mDigitContainer), // digit container for individual CRUs + std::ref(mClusterFinder), // cluster finder for individual CRUs + std::ref(mClusterStorage), // container to store found clusters + std::ref(mClusterDigitIndexStorage), // container to store found cluster MC Labels + cfConfig // configuration ); - else { - processDigits( - std::ref(mDigitContainer[iCRU]), - std::ref(mClusterFinder[iCRU]), - std::ref(mClusterStorage[iCRU]), - std::ref(mClusterDigitIndexStorage[iCRU]), - cfConfig); - } } @@ -450,3 +447,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta mLastTimebin = iTimeBinMax; } +void HwClusterer::setNumThreads(unsigned threads) { + mNumThreads = (threads == 0) ? std::thread::hardware_concurrency() : threads; +} + diff --git a/macro/run_clus_tpc.C b/macro/run_clus_tpc.C index c04a4872fa870..8be966c2e600c 100644 --- a/macro/run_clus_tpc.C +++ b/macro/run_clus_tpc.C @@ -15,7 +15,7 @@ #include "TPCReconstruction/ClustererTask.h" #endif -void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isContinuous=true) +void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isContinuous=true, unsigned threads = 0) { // Initialize logger FairLogger *logger = FairLogger::GetLogger(); @@ -56,8 +56,7 @@ void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isConti // Initialize everything run->Init(); - clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Parallel); -// clustTPC->getHwClusterer()->setProcessingType(o2::TPC::HwClusterer::Processing::Sequential); + clustTPC->getHwClusterer()->setNumThreads(threads); // Start simulation timer.Start(); diff --git a/run/runTPC.cxx b/run/runTPC.cxx index d4b7b73d249a8..2ed2e0caee292 100644 --- a/run/runTPC.cxx +++ b/run/runTPC.cxx @@ -34,7 +34,8 @@ int main(int argc, char *argv[]) ("mode,m", bpo::value()->default_value("sim"), R"(mode of processing, "sim", "digi", "clus", "track" or "all".)") ("nEvents,n", bpo::value()->default_value(2), "number of events to simulate.") ("mcEngine,e", bpo::value()->default_value("TGeant3"), "MC generator to be used.") - ("continuous,c", bpo::value()->default_value(1), "Running in continuous mode 1 - Triggered mode 0"); + ("continuous,c", bpo::value()->default_value(1), "Running in continuous mode 1 - Triggered mode 0") + ("threads,j", bpo::value()->default_value(0), "Parallel processing threads"); bpo::store(parse_command_line(argc, argv, desc), vm); bpo::notify(vm); @@ -49,6 +50,7 @@ int main(int argc, char *argv[]) const std::string engine = vm["mcEngine"].as(); const std::string mode = vm["mode"].as(); const int isContinuous = vm["continuous"].as(); + const unsigned threads = vm["threads"].as(); std::cout << "####" << std::endl; std::cout << "#### Starting TPC simulation tool for" << std::endl; @@ -62,7 +64,7 @@ int main(int argc, char *argv[]) } else if (mode == "digi") { run_digi_tpc(events,engine, isContinuous); } else if (mode == "clus") { - run_clus_tpc(events,engine, isContinuous); + run_clus_tpc(events,engine, isContinuous, threads); } else if (mode == "track") { std::stringstream inputfile, outputfile; inputfile << "AliceO2_" << engine << ".tpc.clusters_" << events << "_event.root"; @@ -82,7 +84,7 @@ int main(int argc, char *argv[]) PID = fork(); if (PID == -1) { std::cout << "ERROR" << std::endl; return EXIT_FAILURE;} - if (PID == 0) { run_clus_tpc(events,engine,isContinuous); return EXIT_SUCCESS;} + if (PID == 0) { run_clus_tpc(events,engine,isContinuous,threads); return EXIT_SUCCESS;} else waitpid(PID,&status,0); PID = fork(); From 0637d572c6e83c59ba80686ac1f64cae21c8ec7f Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Wed, 15 Nov 2017 17:32:53 +0100 Subject: [PATCH 07/12] continued Cluster MCLabels + cleaning up WIP: Cluster MCLabels, cleaning up continued finished Cluster MCLabels + cleaning up --- .../TPCReconstruction/HwClusterFinder.h | 221 ++++++------- .../include/TPCReconstruction/HwClusterer.h | 4 +- .../TPC/reconstruction/src/ClustererTask.cxx | 4 +- .../reconstruction/src/HwClusterFinder.cxx | 295 ++++++------------ .../TPC/reconstruction/src/HwClusterer.cxx | 20 +- macro/run_clus_tpc.C | 4 +- 6 files changed, 220 insertions(+), 328 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index 2d04d5ded74b9..59ba90994011c 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -14,11 +14,10 @@ #ifndef ALICEO2_TPC_HWClusterFinder_H_ #define ALICEO2_TPC_HWClusterFinder_H_ -#include "SimulationDataFormat/MCTruthContainer.h" -#include "SimulationDataFormat/MCCompLabel.h" - #include +#include #include +#include #include namespace o2{ @@ -30,6 +29,8 @@ class Cluster; /// \brief Class for TPC HW cluster finder class HwClusterFinder { public: + /// Mini digit struct, consisting of charge, event number and digit index. + /// The last two are needed to retrieve MC truth information. struct MiniDigit { float charge; int event; @@ -43,93 +44,91 @@ class HwClusterFinder { /// Default Constructor /// \param cru CRU of this cluster finder /// \param row Row of this cluster finder - /// \param id ID for the cluster finder /// \param padOffset Offset in pad direction of the cluster finder /// \param pad Number of pads - /// \param timebins Number of timebins + /// \param timebins Number of time bins /// \param diffThreshold Minimum charge difference at neighboring pads /// \param chargeThreshold Minimum charge of cluster peak /// \param requirePositiveCharge Charge >0 required - HwClusterFinder(short cru, short row, short id, + HwClusterFinder(short cru, short row, short padOffset, short pads=8, short timebins=8, float diffThreshold=0, float chargeThreshold=5, bool requirePositiveCharge=true); /// Destructor - ~HwClusterFinder(); - - /// Copy constructor - HwClusterFinder(const HwClusterFinder& other); - - /// Add a new timebin to cluster finder - /// \param timebin Array of size "length" with new charges - /// \param globalTime Global time of this timebin - /// \param length Size of array "timebin" - bool AddTimebin(MiniDigit* timebin, unsigned globalTime, int length = 8); + ~HwClusterFinder() = default; -// /// Add multiple timebins at once -// /// \param nBins Number of timebins -// /// \param timebins 2D array with new charges -// /// \param globalTime Global time of this timebin -// /// \param length Size of array "timebin" -// bool AddTimebins(int nBins, float** timebins, unsigned globalTimeOfLast, int length = 8); + /// Add a new time bin to cluster finder + /// \param timebin Iterator to a vector with data + /// \param globalTime Global time of this time bin + /// \param length Number of pads to be used after iterator starts + /// \param zeroBin Switch to fill timebin with zero's instead + void addTimebin(std::vector::iterator timebin, unsigned globalTime, int length = 8, bool zeroBin = false); /// Add a timebin with charges of 0 /// \param globalTime Global time of this timebin /// \param length Size of array "timebin" - void AddZeroTimebin(unsigned globalTime = 0, int lengt = 8); + void addZeroTimebin(unsigned globalTime = 0, int lengt = 8); - /// Print the local storagae of charges - void PrintLocalStorage(); + /// Print the local storage of charges + void printLocalStorage(); /// Resets the local storage to zeros - /// \param globalTimeAfterReset Global time of the first timebin after reset + /// \param globalTimeAfterReset Global time of the first time bin after reset void reset(unsigned globalTimeAfterReset); - // Getter functions - int getGlobalTimeOfLast() const { return mGlobalTimeOfLast; } - int getTimebinsAfterLastProcessing() const { return mTimebinsAfterLastProcessing; }; - short getCRU() const { return mCRU; } - short getRow() const { return mRow; } - short getId() const { return mId; } - short getPadOffset() const { return mPadOffset; } - short getNpads() const { return mPads; } - short getNtimebins() const { return mTimebins; } - short getClusterSizeP() const { return mClusterSizePads; } - short getClusterSizeT() const { return mClusterSizeTime; } - float getDiffThreshold() const { return mDiffThreshold; } - float getChargeThreshold() const { return mChargeThreshold; } - bool getRequirePositiveCharge() const { return mRequirePositiveCharge; } - bool getRequireNeighbouringPad() const { return mRequireNeighbouringPad; } - bool getRequireNeighbouringTimebin() const { return mRequireNeighbouringTimebin; } - bool getAutoProcessing() const { return mAutoProcessing; } - bool getmAssignChargeUnique() const { return mAssignChargeUnique; } - HwClusterFinder* getNextCF() const { return mNextCF; } - std::vector* getClusterContainer() { return &clusterContainer; } - std::vector>>* getClusterDigitIndices() { return &clusterDigitIndices; } - - // Setter functions - void setTimebinsAfterLastProcessing(int val) { mTimebinsAfterLastProcessing = val; }; - void setCRU(short val) { mCRU = val; } - void setRow(short val) { mRow = val; } - void setId(short val) { mId = val; } - void setPadOffset(short val) { mPadOffset = val; } - void setNpads(short val) { mPads = val; } - void setNtimebins(short val) { mTimebins = val; } - void setClusterSizeP(short val) { mClusterSizePads = val; } - void setClusterSizeT(short val) { mClusterSizeTime = val; } - void setDiffThreshold(float val) { mDiffThreshold = val; } - void setChargeThreshold(float val) { mChargeThreshold = val; } - void setRequirePositiveCharge(bool val) { mRequirePositiveCharge = val; } - void setRequireNeighbouringPad(bool val) { mRequireNeighbouringPad = val; } - void setRequireNeighbouringTimebin(bool val) { mRequireNeighbouringTimebin = val; } - void setAutoProcessing(bool val) { mAutoProcessing = val; } - void setAssignChargeUnique(bool val) { mAssignChargeUnique = val; } - void setNextCF(HwClusterFinder* nextCF); + /// Getter function + /// \return Time of last inserted time bin + unsigned getGlobalTimeOfLast() const { return mGlobalTimeOfLast; } + + /// Getter function + /// \return Number of inserted time bins since last processing + unsigned getTimebinsAfterLastProcessing() const { return mTimebinsAfterLastProcessing; }; + + /// Getter function + /// \return Configured CRU number + short getCRU() const { return mCRU; } + + /// Getter function + /// \return Configured row number + short getRow() const { return mRow; } + + /// Getter function + /// \return Configured pad offset of this CF + short getPadOffset() const { return mPadOffset; } + + /// Getter function + /// \return Width of CF in pad direction + short getNpads() const { return mPads; } + + /// Getter function + /// \return Width of CF in time direction + short getNtimebins() const { return mTimebins; } + + /// Getter function + /// \return Pointer to cluster container with all found clusters + std::vector* getClusterContainer() { return &mClusterContainer; } + + /// Getter function + /// \return Pointer to container with digits indices used for the clusters + std::vector>>* getClusterDigitIndices() { return &mClusterDigitIndices; } + + + /// Setter function + /// \param val Number of time bins since last processing + void setTimebinsAfterLastProcessing(unsigned val) { mTimebinsAfterLastProcessing = val; }; + + /// Setter function + /// \param val Switch whether charge should be used only for one cluster ("charge splitting", TODO: not yet properly implemented) + void setAssignChargeUnique(bool val) { mAssignChargeUnique = val; } + + /// Setter function + /// \param nextCF Pointer to neighboring CF instance (on the "left" side) + void setNextCF(std::shared_ptr nextCF) { mNextCF = nextCF; }; /// Clears the local cluster storage - void clearClusterContainer() { clusterContainer.clear(); clusterDigitIndices.clear(); } + void clearClusterContainer() { mClusterContainer.clear(); mClusterDigitIndices.clear(); } /// Process the cluster finding bool findCluster(); @@ -137,66 +136,68 @@ class HwClusterFinder { /// Neighboring cluster finder can inform about already used charges /// \param time Time bin of found cluster peak /// \param pad Pad of found cluster peak - /// \param cluster Cluster charges - void clusterAlreadyUsed(short time, short pad, MiniDigit** cluster); + void clusterAlreadyUsed(short time, short pad); private: - MiniDigit chargeForCluster(MiniDigit* charge, MiniDigit* toCompare); + /// Comparator helper function for two MiniDigits. Checks if outerCharge is positiv (if required) and if innerCharge is above threshold + /// \param outerCharge Charge of the "outer" pad + /// \param innerCharge Charge of the "inner" pad + /// \return MiniDigit to be used for cluster. Will be outerCharge, if requirements fulfilled, otherwise "0-MiniDigit" + bool chargeForCluster(float outerCharge, float innerCharge); + + /// Prints 5x5 matrix of internal storage around given parameters + /// \param time relative time bin of center + /// \param pad relative pad number of center void printCluster(short time, short pad); - // local variables - std::vector clusterContainer; - std::vector>> clusterDigitIndices; - int mTimebinsAfterLastProcessing; - MiniDigit** mData; - MiniDigit** tmpCluster; - MiniDigit* mZeroTimebin; - - - // configuration - int mGlobalTimeOfLast; - short mCRU; - short mRow; - short mId; - short mPadOffset; - short mPads; - short mTimebins; - short mClusterSizePads; - short mClusterSizeTime; - float mDiffThreshold; - float mChargeThreshold; - bool mRequirePositiveCharge; - bool mRequireNeighbouringPad; - bool mRequireNeighbouringTimebin; - bool mAutoProcessing; - bool mAssignChargeUnique; - - HwClusterFinder* mNextCF; + /* + * Class members + */ + bool mRequirePositiveCharge; ///< Switch if positive charge is required for individual pad + bool mRequireNeighbouringPad; ///< Switch if at least one neighboring pad needs charge > 0 + bool mRequireNeighbouringTimebin; ///< Switch if at least one neighboring time bin needs charge > 0 + bool mAssignChargeUnique; ///< Switch for "charge splitting", TODO: not yet properly implemented + short mCRU; ///< CRU number + short mRow; ///< Row number + short mPadOffset; ///< Pad number in row of leftmost pad + short mPads; ///< Size of CF in pad direction + short mTimebins; ///< Size of CF in time direction + short mClusterSizePads; ///< Size of cluster in pad direction + short mClusterSizeTime; ///< Size of cluster in time direction + float mDiffThreshold; ///< Charge difference threshold, not yet used + float mChargeThreshold; ///< Charge threshold + unsigned mGlobalTimeOfLast; ///< Global time of last added time bin + unsigned mTimebinsAfterLastProcessing; ///< Number of time bins added after last processing + std::weak_ptr mNextCF; ///< Not owning pointer to neighboring cluster finder (on the "left" side) + + std::vector>> mData; ///< local data storage + std::vector> mTmpCluster; ///< local temporary cluster data storage + std::vector mClusterContainer; ///< Container for found clusters + std::vector>> mClusterDigitIndices; ///< Container for digit indices associated with found clusters }; //________________________________________________________________________ -inline bool HwClusterFinder::AddTimebin(MiniDigit* timebin, unsigned globalTime, int length) +inline void HwClusterFinder::addTimebin(std::vector::iterator timebin, unsigned globalTime, int length, bool zeroBin) { mGlobalTimeOfLast = globalTime; ++mTimebinsAfterLastProcessing; // - // reordering of the local arrays + // reordering of the local array // - MiniDigit* data0 = mData[0]; - std::memmove(mData,mData+1,(mTimebins-1)*sizeof(mData[0])); - mData[mTimebins-1] = data0; - if (length < mPads) { - std::memset(*(mData+mTimebins-1)+length,0 ,(mPads-length)*sizeof(mData[mTimebins-1][0])); - std::memcpy(*(mData+mTimebins-1),timebin,length*sizeof(timebin[0])); + std::rotate(mData.begin(), mData.begin() + 1, mData.end()); + + // + // fillin with data + // + if (zeroBin) { + for (auto &digi : *mData.back()) digi.clear(); } else { - std::memcpy(*(mData+mTimebins-1),timebin,mPads*sizeof(timebin[0])); + std::copy(timebin,timebin+length,mData.back()->begin()); + mData.back()->resize(mPads,MiniDigit()); } - - if (mAutoProcessing & (mTimebinsAfterLastProcessing >= (mTimebins -2 -2))) findCluster(); - return true; } } diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index 101409de86507..ec7639fcc9276 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -128,7 +128,7 @@ class HwClusterer : public Clusterer { /// \param config Configuration for the cluster finding static void processDigits( const std::vector>>>& digits, - const std::vector>>>& clusterFinder, + const std::vector>>>& clusterFinder, std::vector>& cluster, std::vector>>>& label, CfConfig config); @@ -156,7 +156,7 @@ class HwClusterer : public Clusterer { unsigned mNumThreads; ///< Number of parallel processing threads float mMinQDiff; ///< Minimum charge difference between neighboring pads / time bins - std::vector>>> mClusterFinder; ///< Cluster finder container for each row in each CRU + std::vector>>> mClusterFinder; ///< Cluster finder container for each row in each CRU std::vector>>> mDigitContainer; ///< Sorted digit container for each row in each CRU. Tuple consists of pointer to digit, original digit index and event count std::vector> mClusterStorage; ///< Cluster storage for each CRU diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 5f499b8154b74..8309e5595baf3 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -125,14 +125,14 @@ void ClustererTask::Exec(Option_t *option) mClustersArray->clear(); if(mClustersMCTruthArray) mClustersMCTruthArray->clear(); mBoxClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); - LOG(DEBUG) << "Box clusterer found " << mClustersArray->size() << " clusters" << FairLogger::endl; + LOG(DEBUG) << "Box clusterer found " << mClustersArray->size() << " clusters" << FairLogger::endl << FairLogger::endl; } if (mHwClustererEnable) { if(mHwClustersArray) mHwClustersArray->clear(); if(mHwClustersMCTruthArray) mHwClustersMCTruthArray->clear(); mHwClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); - LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl; + LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl << FairLogger::endl; } ++mEventCount; } diff --git a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx index 8c0f1bf3e2107..8d456f4e1a271 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx @@ -24,150 +24,63 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterFinder::HwClusterFinder( - short cru, short row, short id, + short cru, short row, short padOffset, short pads, short timebins, float diffThreshold, float chargeThreshold, bool requirePositiveCharge) - : mGlobalTimeOfLast(0) - , mTimebinsAfterLastProcessing(0) + : mRequirePositiveCharge(requirePositiveCharge) + , mRequireNeighbouringPad(false)//true) + , mRequireNeighbouringTimebin(true) + , mAssignChargeUnique(false)//true) , mCRU(cru) , mRow(row) - , mId(id) , mPadOffset(padOffset) , mPads(pads) , mTimebins(timebins) , mClusterSizePads(5) , mClusterSizeTime(5) - , mDiffThreshold(diffThreshold) + , mDiffThreshold(diffThreshold) // not yet used , mChargeThreshold(chargeThreshold) - , mRequirePositiveCharge(requirePositiveCharge) - , mRequireNeighbouringPad(false)//true) - , mRequireNeighbouringTimebin(true) - , mAutoProcessing(false) - , mAssignChargeUnique(false)//true) - , mData(nullptr) - , tmpCluster(nullptr) - , mZeroTimebin(nullptr) - , mNextCF(nullptr) + , mGlobalTimeOfLast(0) + , mTimebinsAfterLastProcessing(0) + , mNextCF() + , mData() // couldn't find a way to initialize a vector of unique_ptr's here + , mTmpCluster(5, std::vector(5, MiniDigit())) + , mClusterContainer() + , mClusterDigitIndices() { if (mPads < mClusterSizePads) { LOG(ERROR) << "Given width in pad direction is smaller than cluster size in pad direction." - << " width in pad direction was increased to cluster size." << FairLogger::endl; + << " Width in pad direction was increased to cluster size " << mClusterSizePads << FairLogger::endl; mPads = mClusterSizePads; } if (mTimebins < mClusterSizeTime) { LOG(ERROR) << "Given width in time direction is smaller than cluster size in time direction." - << " width in time direction was increased to cluster size." << FairLogger::endl; + << " Width in time direction was increased to cluster size " << mClusterSizeTime << FairLogger::endl; mTimebins = mClusterSizeTime; } - if (mTimebins*mPads > 64) { - LOG(WARNING) << "Bins in pad direction X bins in time direction is larger than 64." << FairLogger::endl; - } - - short t,p; - mData = new MiniDigit*[mTimebins]; - for (t = 0; t < mTimebins; ++t){ - mData[t] = new MiniDigit[mPads]; -// for (p = 0; p < mPads; ++p){ -// mData[t][p] = MiniDigit(); -// } - } + for (short t = 0; t < mTimebins; ++t) + mData.emplace_back(std::make_unique>(pads, MiniDigit())); - tmpCluster = new MiniDigit*[mClusterSizeTime]; - for (t=0; tbegin() /* some iterator, is not used */,globalTime,length,true); } //________________________________________________________________________ -HwClusterFinder::~HwClusterFinder() +void HwClusterFinder::printLocalStorage() { - short t; - for (t = 0; t < mTimebins; ++t){ - delete [] mData[t]; - } - delete [] mData; - delete [] mZeroTimebin; - - for (t=0; tat(p).charge < mChargeThreshold) continue; // Require at least one neighboring time bin with signal - if (mRequireNeighbouringTimebin && (mData[t-1][p ].charge + mData[t+1][p ].charge <= 0)) continue; + if (mRequireNeighbouringTimebin && (mData[t-1]->at(p ).charge <= 0 && mData[t+1]->at(p ).charge <= 0)) continue; // Require at least one neighboring pad with signal - if (mRequireNeighbouringPad && (mData[t ][p-1].charge + mData[t ][p+1].charge <= 0)) continue; + if (mRequireNeighbouringPad && (mData[t ]->at(p-1).charge <= 0 && mData[t ]->at(p+1).charge <= 0)) continue; // check for local maximum - if (mData[t-1][p ].charge >= mData[t][p].charge) continue; - if (mData[t+1][p ].charge > mData[t][p].charge) continue; - if (mData[t ][p-1].charge >= mData[t][p].charge) continue; - if (mData[t ][p+1].charge > mData[t][p].charge) continue; - if (mData[t-1][p-1].charge >= mData[t][p].charge) continue; - if (mData[t+1][p+1].charge > mData[t][p].charge) continue; - if (mData[t+1][p-1].charge > mData[t][p].charge) continue; - if (mData[t-1][p+1].charge >= mData[t][p].charge) continue; + if (mData[t-1]->at(p ).charge >= mData[t]->at(p).charge) continue; + if (mData[t+1]->at(p ).charge > mData[t]->at(p).charge) continue; + if (mData[t ]->at(p-1).charge >= mData[t]->at(p).charge) continue; + if (mData[t ]->at(p+1).charge > mData[t]->at(p).charge) continue; + if (mData[t-1]->at(p-1).charge >= mData[t]->at(p).charge) continue; + if (mData[t+1]->at(p+1).charge > mData[t]->at(p).charge) continue; + if (mData[t+1]->at(p-1).charge > mData[t]->at(p).charge) continue; + if (mData[t-1]->at(p+1).charge >= mData[t]->at(p).charge) continue; // printf("##\n"); -// printf("## cluster found at t=%d, p=%d (in CF %d in row %d of CRU %d)\n",t,p,mId,mRow,mCRU); +// printf("## cluster found at t=%d, p=%d (in row %d of CRU %d)\n",t,p,mRow,mCRU); // printf("##\n"); // printCluster(t,p); @@ -263,7 +176,7 @@ bool HwClusterFinder::findCluster() // prepare temp storage for (tt=0; ttat(p+(pp-2)).charge < 0) continue; + mTmpCluster[tt][pp] = mData[t+(tt-2)]->at(p+(pp-2)); // mData[t+(tt-2)][p+(pp-2)] = 0; } } @@ -295,11 +208,11 @@ bool HwClusterFinder::findCluster() // // [p] 0 1 2 3 4 - //tmpCluster[t][p] - tmpCluster[0][2] = chargeForCluster(&mData[t-2][p ],&mData[t-1][p ]); // t-X -> older - tmpCluster[4][2] = chargeForCluster(&mData[t+2][p ],&mData[t+1][p ]); // t+X -> newer - tmpCluster[2][0] = chargeForCluster(&mData[t ][p-2],&mData[t ][p-1]); - tmpCluster[2][4] = chargeForCluster(&mData[t ][p+2],&mData[t ][p+1]); + // o i mTmpCluster[t][p] + if(chargeForCluster(mData[t-2]->at(p ).charge, mData[t-1]->at(p ).charge)) mTmpCluster[0][2] = mData[t-2]->at(p ); // t-X -> older + if(chargeForCluster(mData[t+2]->at(p ).charge, mData[t+1]->at(p ).charge)) mTmpCluster[4][2] = mData[t+2]->at(p ); // t+X -> newer + if(chargeForCluster(mData[t ]->at(p-2).charge, mData[t ]->at(p-1).charge)) mTmpCluster[2][0] = mData[t ]->at(p-2); + if(chargeForCluster(mData[t ]->at(p+2).charge, mData[t ]->at(p+1).charge)) mTmpCluster[2][4] = mData[t ]->at(p+2); // The cells of the corners have 3 neighbours. @@ -310,27 +223,27 @@ bool HwClusterFinder::findCluster() // o o o o // bottom left - tmpCluster[3][0] = chargeForCluster(&mData[t+1][p-2],&mData[t+1][p-1]); - tmpCluster[4][0] = chargeForCluster(&mData[t+2][p-2],&mData[t+1][p-1]); - tmpCluster[4][1] = chargeForCluster(&mData[t+2][p-1],&mData[t+1][p-1]); + if(chargeForCluster(mData[t+1]->at(p-2).charge,mData[t+1]->at(p-1).charge)) mTmpCluster[3][0] = mData[t+1]->at(p-2); + if(chargeForCluster(mData[t+2]->at(p-2).charge,mData[t+1]->at(p-1).charge)) mTmpCluster[4][0] = mData[t+2]->at(p-2); + if(chargeForCluster(mData[t+2]->at(p-1).charge,mData[t+1]->at(p-1).charge)) mTmpCluster[4][1] = mData[t+2]->at(p-1); // bottom right - tmpCluster[4][3] = chargeForCluster(&mData[t+2][p+1],&mData[t+1][p+1]); - tmpCluster[4][4] = chargeForCluster(&mData[t+2][p+2],&mData[t+1][p+1]); - tmpCluster[3][4] = chargeForCluster(&mData[t+1][p+2],&mData[t+1][p+1]); + if(chargeForCluster(mData[t+2]->at(p+1).charge,mData[t+1]->at(p+1).charge)) mTmpCluster[4][3] = mData[t+2]->at(p+1); + if(chargeForCluster(mData[t+2]->at(p+2).charge,mData[t+1]->at(p+1).charge)) mTmpCluster[4][4] = mData[t+2]->at(p+2); + if(chargeForCluster(mData[t+1]->at(p+2).charge,mData[t+1]->at(p+1).charge)) mTmpCluster[3][4] = mData[t+1]->at(p+2); // top right - tmpCluster[1][4] = chargeForCluster(&mData[t-1][p+2],&mData[t-1][p+1]); - tmpCluster[0][4] = chargeForCluster(&mData[t-2][p+2],&mData[t-1][p+1]); - tmpCluster[0][3] = chargeForCluster(&mData[t-2][p+1],&mData[t-1][p+1]); + if(chargeForCluster(mData[t-1]->at(p+2).charge,mData[t-1]->at(p+1).charge)) mTmpCluster[1][4] = mData[t-1]->at(p+2); + if(chargeForCluster(mData[t-2]->at(p+2).charge,mData[t-1]->at(p+1).charge)) mTmpCluster[0][4] = mData[t-2]->at(p+2); + if(chargeForCluster(mData[t-2]->at(p+1).charge,mData[t-1]->at(p+1).charge)) mTmpCluster[0][3] = mData[t-2]->at(p+1); // top left - tmpCluster[0][1] = chargeForCluster(&mData[t-2][p-1],&mData[t-1][p-1]); - tmpCluster[0][0] = chargeForCluster(&mData[t-2][p-2],&mData[t-1][p-1]); - tmpCluster[1][0] = chargeForCluster(&mData[t-1][p-2],&mData[t-1][p-1]); + if(chargeForCluster(mData[t-2]->at(p-1).charge,mData[t-1]->at(p-1).charge)) mTmpCluster[0][1] = mData[t-2]->at(p-1); + if(chargeForCluster(mData[t-2]->at(p-2).charge,mData[t-1]->at(p-1).charge)) mTmpCluster[0][0] = mData[t-2]->at(p-2); + if(chargeForCluster(mData[t-1]->at(p-2).charge,mData[t-1]->at(p-1).charge)) mTmpCluster[1][0] = mData[t-1]->at(p-2); // // calculate cluster Properties // - qMax = tmpCluster[2][2].charge; + qMax = mTmpCluster[2][2].charge; qTot = 0; meanP = 0; meanT = 0; @@ -341,23 +254,23 @@ bool HwClusterFinder::findCluster() maxP = 0; maxT = 0; clusterSize = 0; - clusterDigitIndices.emplace_back(); + mClusterDigitIndices.emplace_back(); for (tt = 0; tt < mClusterSizeTime; ++tt) { deltaT = tt - mClusterSizeTime/2; for (pp = 0; pp < mClusterSizePads; ++pp) { deltaP = pp - mClusterSizePads/2; - charge = tmpCluster[tt][pp].charge; - if (charge > 0 || tmpCluster[tt][pp].event >= 0) - clusterDigitIndices.back().emplace_back(std::make_pair(tmpCluster[tt][pp].index, tmpCluster[tt][pp].event)); + charge = mTmpCluster[tt][pp].charge; + if (charge > 0 || mTmpCluster[tt][pp].event >= 0) + mClusterDigitIndices.back().emplace_back(std::make_pair(mTmpCluster[tt][pp].index, mTmpCluster[tt][pp].event)); qTot += charge; - meanP += charge * deltaP; - meanT += charge * deltaT; + meanP += charge * static_cast(deltaP); + meanT += charge * static_cast(deltaT); - sigmaP += charge * deltaP*deltaP; - sigmaT += charge * deltaT*deltaT; + sigmaP += charge * static_cast(deltaP)*static_cast(deltaP); + sigmaT += charge * static_cast(deltaT)*static_cast(deltaT); if (charge > 0) { minP = std::min(minP,pp); maxP = std::max(maxP,pp); @@ -382,14 +295,15 @@ bool HwClusterFinder::findCluster() } ++foundNclusters; -// clusterContainer.emplace_back(mCRU, mRow, mClusterSizePads, mClusterSizeTime, tmpCluster,p+mPadOffset,mGlobalTimeOfLast-(mTimebins-1)+t); - clusterContainer.emplace_back(mCRU, mRow, qTot, qMax, meanP, sigmaP, meanT, sigmaT); + mClusterContainer.emplace_back(mCRU, mRow, qTot, qMax, meanP, sigmaP, meanT, sigmaT); if (mAssignChargeUnique) { if (p < (pMin+4)) { // If the cluster peak is in one of the 6 leftmost pads, the Cluster Finder // on the left has to know about it to ignore the already used pads. - if (mNextCF != nullptr) mNextCF->clusterAlreadyUsed(t,p+mPadOffset,tmpCluster); + if (auto next = mNextCF.lock()) { + next->clusterAlreadyUsed(t,p+mPadOffset);//,mTmpCluster); + } } @@ -399,8 +313,8 @@ bool HwClusterFinder::findCluster() // TODO: really nexessary?? or just set to 0 for (tt=0; tt<5; ++tt) { for (pp=0; pp<5; ++pp) { - //mData[t+(tt-2)][p+(pp-2)].charge -= tmpCluster[tt][pp]; - mData[t+(tt-2)][p+(pp-2)].clear();// = MiniDigit(); + //mData[t+(tt-2)][p+(pp-2)].charge -= mTmpCluster[tt][pp]; + mData[t+(tt-2)]->at(p+(pp-2)).clear(); } } } @@ -412,43 +326,19 @@ bool HwClusterFinder::findCluster() } //________________________________________________________________________ -HwClusterFinder::MiniDigit HwClusterFinder::chargeForCluster(MiniDigit* charge, MiniDigit* toCompare) +bool HwClusterFinder::chargeForCluster(float outerCharge, float innerCharge) { //printf("%.2f - %.2f = %.f compared to %.2f)?\n",toCompare,*charge,toCompare-*charge,-mDiffThreshold); - if ((mRequirePositiveCharge && (charge->charge > 0)) & - (toCompare->charge > mDiffThreshold)) {//mChargeThreshold)) { - //printf("\tyes\n"); - return *charge; - } else { - //printf("\tno\n"); - return MiniDigit(); + if ((mRequirePositiveCharge && (outerCharge > 0)) && + (innerCharge > mChargeThreshold)) { + return true; } -} - -//________________________________________________________________________ -void HwClusterFinder::setNextCF(HwClusterFinder* nextCF) -{ - if (nextCF == nullptr) { - LOG(WARNING) << R"(Got "nullptrd" as neighboring Cluster Finder.)" << FairLogger::endl; - return; - } - if (mNextCF != nullptr) { - LOG(WARNING) << "This Cluster Finder (" - << "CRU " << mCRU << ", row " << mRow << ",pad offset " << mPadOffset - << ") had already a neighboring CF set (" - << "CRU " << mNextCF->getCRU() << ", row " << mNextCF->getRow() << ",pad offset " << mNextCF->getPadOffset() - << "). It will be replaced with a new one (" - << "CRU " << nextCF->getCRU() << ", row " << nextCF->getRow() << ",pad offset " << nextCF->getPadOffset() - << ")." << FairLogger::endl; - } - mNextCF = nextCF; - - return; + return false; } //________________________________________________________________________ // TODO: really nexessary?? or just set to 0 -void HwClusterFinder::clusterAlreadyUsed(short time, short pad, MiniDigit** cluster) +void HwClusterFinder::clusterAlreadyUsed(short time, short pad) { short localPad = pad - mPadOffset; @@ -459,7 +349,7 @@ void HwClusterFinder::clusterAlreadyUsed(short time, short pad, MiniDigit** clus if (p < 0 || p >= mPads) continue; //mData[t][p].charge -= cluster[t-time+2][p-localPad+2].charge; - mData[t][p].clear();// = MiniDigit(); + mData[t]->at(p).clear(); } } } @@ -467,12 +357,8 @@ void HwClusterFinder::clusterAlreadyUsed(short time, short pad, MiniDigit** clus //________________________________________________________________________ void HwClusterFinder::reset(unsigned globalTimeAfterReset) { - short t,p; - for (t = 0; t < mTimebins; ++t){ - for (p = 0; p < mPads; ++p){ - mData[t][p].clear();// = MiniDigit(); - } - } + for(auto &tb : mData) + for (auto &digi : *tb) digi.clear(); mGlobalTimeOfLast = globalTimeAfterReset; } @@ -482,10 +368,11 @@ void HwClusterFinder::printCluster(short time, short pad) { short t,p; for (t = time-2; t <= time+2; ++t) { - printf("%d\t\t",t); + LOG(DEBUG) << "t " << t << ":\t"; for (p = pad-2; p <= pad+2; ++p) { - printf("%.2f\t", mData[t][p].charge); + LOG(DEBUG) << mData[t]->at(p).charge << "\t"; } - printf("\n"); + LOG(DEBUG) << FairLogger::endl; } + LOG(DEBUG) << FairLogger::endl; } diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 02128cb88d5ad..8b84dbffda13d 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -68,7 +68,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, mClusterFinder[iCRU][iRow].resize(iCfPerRow); for (int iCF = 0; iCF < iCfPerRow; iCF++){ int padOffset = iCF*(mPadsPerCF-2-2)-2; - mClusterFinder[iCRU][iRow][iCF] = std::make_unique(iCRU,iRow,iCF,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); + mClusterFinder[iCRU][iRow][iCF] = std::make_shared(iCRU,iRow,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); mClusterFinder[iCRU][iRow][iCF]->setAssignChargeUnique(mAssignChargeUnique); @@ -78,7 +78,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, * already used for a cluster. */ if (iCF != 0) { - mClusterFinder[iCRU][iRow][iCF]->setNextCF(mClusterFinder[iCRU][iRow][iCF-1].get()); + mClusterFinder[iCRU][iRow][iCF]->setNextCF(mClusterFinder[iCRU][iRow][iCF-1]); } } } @@ -115,7 +115,7 @@ HwClusterer::~HwClusterer() //________________________________________________________________________ void HwClusterer::processDigits( const std::vector>>>& digits, - const std::vector>>>& clusterFinder, + const std::vector>>>& clusterFinder, std::vector>& cluster, std::vector>>>& label, const CfConfig config) @@ -128,7 +128,7 @@ void HwClusterer::processDigits( int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; if (timeDiff < 0) return; const Mapper& mapper = Mapper::instance(); - HwClusterFinder::MiniDigit iAllBins[timeDiff][config.iMaxPads]; + std::vector> iAllBins(timeDiff,std::vector(config.iMaxPads)); for (int iCRU = config.iCRUMin; iCRU <= config.iCRUMax; ++iCRU) { if (iCRU % config.iThreadMax != config.iThreadID) continue; @@ -149,7 +149,8 @@ void HwClusterer::processDigits( } } } else { - std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, HwClusterFinder::MiniDigit()); + for (auto &bins : iAllBins) std::fill(bins.begin(),bins.end(),HwClusterFinder::MiniDigit()); +// std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, HwClusterFinder::MiniDigit()); } /* @@ -176,12 +177,12 @@ void HwClusterer::processDigits( */ const Short_t iPadsPerCF = clusterFinder[iCRU][iRow][0]->getNpads(); const Short_t iTimebinsPerCF = clusterFinder[iCRU][iRow][0]->getNtimebins(); - std::vector>::const_reverse_iterator> cfWithCluster; + std::vector>::const_reverse_iterator> cfWithCluster; unsigned time,pad; for (time = 0; time < timeDiff; ++time){ // ordering important!! for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { const Short_t cf = pad / (iPadsPerCF-2-2); - clusterFinder[iCRU][iRow][cf]->AddTimebin(&iAllBins[time][pad],time+config.iMinTimeBin,(config.iMaxPads-pad)>=iPadsPerCF?iPadsPerCF:(config.iMaxPads-pad)); + clusterFinder[iCRU][iRow][cf]->addTimebin(iAllBins[time].begin()+pad,time+config.iMinTimeBin,(config.iMaxPads-pad)>=iPadsPerCF?iPadsPerCF:(config.iMaxPads-pad)); } /* @@ -206,7 +207,7 @@ void HwClusterer::processDigits( // +2 so that for sure all data is processed for (time = 0; time < clusterFinder[iCRU][iRow][0]->getNtimebins()+2; ++time){ for (auto rit = clusterFinder[iCRU][iRow].crbegin(); rit != clusterFinder[iCRU][iRow].crend(); ++rit) { - (*rit)->AddZeroTimebin(time+timeDiff+config.iMinTimeBin,iPadsPerCF); + (*rit)->addZeroTimebin(time+timeDiff+config.iMinTimeBin,iPadsPerCF); } /* @@ -315,6 +316,8 @@ void HwClusterer::Process(std::vector const &digits, MCLabelCont ProcessTimeBins(iTimeBinMin, iTimeBinMax, mcDigitTruth, eventCount); mLastMcDigitTruth.erase(eventCount-mTimebinsPerCF); + + LOG(DEBUG) << "Event ranged from time bin " << iTimeBinMin << " to " << iTimeBinMax << "." << FairLogger::endl; } //________________________________________________________________________ @@ -381,6 +384,7 @@ void HwClusterer::Process(std::vector>& digits, MCLabelCo mLastMcDigitTruth.erase(eventCount-mTimebinsPerCF); + LOG(DEBUG) << "Event ranged from time bin " << iTimeBinMin << " to " << iTimeBinMax << "." << FairLogger::endl; } void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelContainer const* mcDigitTruth, int eventCount) diff --git a/macro/run_clus_tpc.C b/macro/run_clus_tpc.C index 8be966c2e600c..24f3c977c63a1 100644 --- a/macro/run_clus_tpc.C +++ b/macro/run_clus_tpc.C @@ -32,7 +32,7 @@ void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isConti TStopwatch timer; // Setup FairRoot analysis manager - FairRunAna * run = new FairRunAna; + FairRunAna * run = new FairRunAna(); FairFileSource *fFileSource = new FairFileSource(inputfile.str().c_str()); run->SetSource(fFileSource); run->SetOutputFile(outputfile.str().c_str()); @@ -46,7 +46,7 @@ void run_clus_tpc(Int_t nEvents = 10, TString mcEngine = "TGeant3", bool isConti TGeoManager::Import("geofile_full.root"); // Setup clusterer - o2::TPC::ClustererTask *clustTPC = new o2::TPC::ClustererTask; + o2::TPC::ClustererTask *clustTPC = new o2::TPC::ClustererTask(); clustTPC->setContinuousReadout(isContinuous); clustTPC->setClustererEnable(o2::TPC::ClustererTask::ClustererType::Box,false); clustTPC->setClustererEnable(o2::TPC::ClustererTask::ClustererType::HW,true); From 93cb1b1b9e04c4bb1be9e963bf21ec44b5459706 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Fri, 17 Nov 2017 09:51:50 +0100 Subject: [PATCH 08/12] fixed some compiler warnings, comparison signed/unsigned --- .../TPCReconstruction/HwClusterFinder.h | 28 ++++++++--------- .../include/TPCReconstruction/HwClusterer.h | 2 +- .../include/TPCReconstruction/RawReader.h | 6 ++-- .../reconstruction/src/HwClusterFinder.cxx | 12 ++------ .../TPC/reconstruction/src/HwClusterer.cxx | 30 +++++++++---------- .../TPC/reconstruction/src/RawReader.cxx | 8 ++--- 6 files changed, 39 insertions(+), 47 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index 59ba90994011c..3a3a00728b978 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -50,8 +50,8 @@ class HwClusterFinder { /// \param diffThreshold Minimum charge difference at neighboring pads /// \param chargeThreshold Minimum charge of cluster peak /// \param requirePositiveCharge Charge >0 required - HwClusterFinder(short cru, short row, - short padOffset, short pads=8, short timebins=8, + HwClusterFinder(unsigned short cru, unsigned short row, + unsigned short padOffset, unsigned short pads=8, unsigned short timebins=8, float diffThreshold=0, float chargeThreshold=5, bool requirePositiveCharge=true); /// Destructor @@ -87,23 +87,23 @@ class HwClusterFinder { /// Getter function /// \return Configured CRU number - short getCRU() const { return mCRU; } + unsigned short getCRU() const { return mCRU; } /// Getter function /// \return Configured row number - short getRow() const { return mRow; } + unsigned short getRow() const { return mRow; } /// Getter function /// \return Configured pad offset of this CF - short getPadOffset() const { return mPadOffset; } + unsigned short getPadOffset() const { return mPadOffset; } /// Getter function /// \return Width of CF in pad direction - short getNpads() const { return mPads; } + unsigned short getNpads() const { return mPads; } /// Getter function /// \return Width of CF in time direction - short getNtimebins() const { return mTimebins; } + unsigned short getNtimebins() const { return mTimebins; } /// Getter function /// \return Pointer to cluster container with all found clusters @@ -158,13 +158,13 @@ class HwClusterFinder { bool mRequireNeighbouringPad; ///< Switch if at least one neighboring pad needs charge > 0 bool mRequireNeighbouringTimebin; ///< Switch if at least one neighboring time bin needs charge > 0 bool mAssignChargeUnique; ///< Switch for "charge splitting", TODO: not yet properly implemented - short mCRU; ///< CRU number - short mRow; ///< Row number - short mPadOffset; ///< Pad number in row of leftmost pad - short mPads; ///< Size of CF in pad direction - short mTimebins; ///< Size of CF in time direction - short mClusterSizePads; ///< Size of cluster in pad direction - short mClusterSizeTime; ///< Size of cluster in time direction + unsigned short mCRU; ///< CRU number + unsigned short mRow; ///< Row number + unsigned short mPadOffset; ///< Pad number in row of leftmost pad + unsigned short mPads; ///< Size of CF in pad direction + unsigned short mTimebins; ///< Size of CF in time direction + unsigned short mClusterSizePads; ///< Size of cluster in pad direction + unsigned short mClusterSizeTime; ///< Size of cluster in time direction float mDiffThreshold; ///< Charge difference threshold, not yet used float mChargeThreshold; ///< Charge threshold unsigned mGlobalTimeOfLast; ///< Global time of last added time bin diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index ec7639fcc9276..ca4119e17fbb5 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -110,7 +110,7 @@ class HwClusterer : public Clusterer { unsigned iThreadMax; ///< Total number of started threads unsigned iCRUMin; ///< Minimum CRU number to process unsigned iCRUMax; ///< Maximum CRU number to process - int iMaxPads; ///< Maximum number of pads per row + unsigned iMaxPads; ///< Maximum number of pads per row int iMinTimeBin; ///< Minumum digit time bin int iMaxTimeBin; ///< Maximum digit time bin bool iEnableNoiseSim; ///< Noise simulation enable switch diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReader.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReader.h index f2f3495d92a78..906b012b6d069 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReader.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReader.h @@ -198,12 +198,12 @@ class RawReader { bool decodeRawGBTFrames(EventInfo eventInfo); bool decodePreprocessedData(EventInfo eventInfo); - int mRegion; ///< Region of the data - int mLink; ///< FEC of the data - int mRun; ///< Run number bool mUseRawInMode3; ///< in readout mode 3 decode GBT frames bool mApplyChannelMask; ///< apply channel mask bool mCheckAdcClock; ///< check the ADC clock + int mRegion; ///< Region of the data + int mLink; ///< FEC of the data + int mRun; ///< Run number int64_t mLastEvent; ///< Number of last loaded event std::array mTimestampOfFirstData; ///< Time stamp of first decoded ADC value, individually for each half SAMPA std::map>> mEvents; ///< all "event data" - headers, file path, etc. NOT actual data diff --git a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx index 8d456f4e1a271..65a36b56a46ff 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx @@ -24,8 +24,8 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterFinder::HwClusterFinder( - short cru, short row, - short padOffset, short pads, short timebins, + unsigned short cru, unsigned short row, + unsigned short padOffset, unsigned short pads, unsigned short timebins, float diffThreshold, float chargeThreshold, bool requirePositiveCharge) : mRequirePositiveCharge(requirePositiveCharge) @@ -111,16 +111,12 @@ bool HwClusterFinder::findCluster() // short tMax = (mTimebins-1) - 2; short tMin = 2; - short delTm = -2; - short delTp = 2; // // In pad direction // int pMin = 2; int pMax = (mPads-1)-2; - short delPm = -2; - short delPp = 2; float qMax; float qTot; @@ -130,7 +126,6 @@ bool HwClusterFinder::findCluster() short minP, minT; short maxP, maxT; short deltaP, deltaT; - short clusterSize; // // peak finding // @@ -253,7 +248,6 @@ bool HwClusterFinder::findCluster() minT = mClusterSizeTime; maxP = 0; maxT = 0; - clusterSize = 0; mClusterDigitIndices.emplace_back(); for (tt = 0; tt < mClusterSizeTime; ++tt) { deltaT = tt - mClusterSizeTime/2; @@ -279,8 +273,6 @@ bool HwClusterFinder::findCluster() } } - clusterSize = (maxP-minP+1)*10 + (maxT-minT+1); - if (qTot > 0) { meanP /= qTot; meanT /= qTot; diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 8b84dbffda13d..b711bc27ca406 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -59,14 +59,14 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, /* * initialize all cluster finder */ - int iCfPerRow = (int)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); + unsigned iCfPerRow = (unsigned)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); mClusterFinder.resize(mCRUMax+1); const Mapper& mapper = Mapper::instance(); - for (int iCRU = mCRUMin; iCRU <= mCRUMax; iCRU++){ + for (unsigned iCRU = mCRUMin; iCRU <= mCRUMax; ++iCRU){ mClusterFinder[iCRU].resize(mapper.getNumberOfRowsPartition(iCRU)); - for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); iRow++){ + for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); ++iRow){ mClusterFinder[iCRU][iRow].resize(iCfPerRow); - for (int iCF = 0; iCF < iCfPerRow; iCF++){ + for (unsigned iCF = 0; iCF < iCfPerRow; ++iCF){ int padOffset = iCF*(mPadsPerCF-2-2)-2; mClusterFinder[iCRU][iRow][iCF] = std::make_shared(iCRU,iRow,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); mClusterFinder[iCRU][iRow][iCF]->setAssignChargeUnique(mAssignChargeUnique); @@ -99,7 +99,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, * CRU (thread) */ mDigitContainer.resize(mCRUMax+1); - for (int iCRU = mCRUMin; iCRU <= mCRUMax; iCRU++) + for (unsigned iCRU = mCRUMin; iCRU <= mCRUMax; ++iCRU) mDigitContainer[iCRU].resize(mapper.getNumberOfRowsPartition(iCRU)); } @@ -130,7 +130,7 @@ void HwClusterer::processDigits( const Mapper& mapper = Mapper::instance(); std::vector> iAllBins(timeDiff,std::vector(config.iMaxPads)); - for (int iCRU = config.iCRUMin; iCRU <= config.iCRUMax; ++iCRU) { + for (unsigned iCRU = config.iCRUMin; iCRU <= config.iCRUMax; ++iCRU) { if (iCRU % config.iThreadMax != config.iThreadID) continue; for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); iRow++){ @@ -139,7 +139,6 @@ void HwClusterer::processDigits( * prepare local storage */ short t,p; - float noise; if (config.iEnableNoiseSim && config.iNoiseObject != nullptr) { for (t=timeDiff; t--;) { for (p=config.iMaxPads; p--;) { @@ -175,10 +174,11 @@ void HwClusterer::processDigits( /* * copy data to cluster finders */ - const Short_t iPadsPerCF = clusterFinder[iCRU][iRow][0]->getNpads(); - const Short_t iTimebinsPerCF = clusterFinder[iCRU][iRow][0]->getNtimebins(); + const unsigned iPadsPerCF = static_cast(clusterFinder[iCRU][iRow][0]->getNpads()); + const unsigned iTimebinsPerCF = static_cast(clusterFinder[iCRU][iRow][0]->getNtimebins()); std::vector>::const_reverse_iterator> cfWithCluster; - unsigned time,pad; + int time; + unsigned pad; for (time = 0; time < timeDiff; ++time){ // ordering important!! for (pad = 0; pad < config.iMaxPads; pad = pad + (iPadsPerCF -2 -2 )) { const Short_t cf = pad / (iPadsPerCF-2-2); @@ -287,7 +287,7 @@ void HwClusterer::Process(std::vector const &digits, MCLabelCont */ iTimeBin = digit.getTimeStamp(); - if (digit.getCRU() < mCRUMin || digit.getCRU() > mCRUMax) { + if (digit.getCRU() < static_cast(mCRUMin) || digit.getCRU() > static_cast(mCRUMax)) { LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit.getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); @@ -354,7 +354,7 @@ void HwClusterer::Process(std::vector>& digits, MCLabelCo * add current digit to storage */ iTimeBin = digit->getTimeStamp(); - if (digit->getCRU() < mCRUMin || digit->getCRU() > mCRUMax) { + if (digit->getCRU() < static_cast(mCRUMin) || digit->getCRU() > static_cast(mCRUMax)) { LOG(DEBUG) << "Digit [" << digitIndex << "] is out of CRU range (" << digit->getCRU() << " < " << mCRUMin << " or > " << mCRUMax << ")" << FairLogger::endl; // Necessary because MCTruthContainer requires continuous indexing if (mcDigitTruth != nullptr) mLastMcDigitTruth[eventCount]->addElement(digitIndex,o2::MCCompLabel(-1,-1,-1)); @@ -403,7 +403,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta mNumThreads, mCRUMin, mCRUMax, - mPadsMax+2+2, + static_cast(mPadsMax)+2+2, iTimeBinMin, iTimeBinMax, mEnableNoiseSim, @@ -433,11 +433,11 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta /* * collect clusters from individual cluster finder */ - for (int cru = 0; cru < mClusterStorage.size(); ++cru) { + for (unsigned cru = 0; cru < mClusterStorage.size(); ++cru) { std::vector* clustersFromCRU = &mClusterStorage[cru]; std::vector>>* labelsFromCRU = &mClusterDigitIndexStorage[cru]; - for(int c = 0; c < clustersFromCRU->size(); ++c) { + for(unsigned c = 0; c < clustersFromCRU->size(); ++c) { const auto clusterPos = mClusterArray->size(); mClusterArray->emplace_back(clustersFromCRU->at(c)); for (auto &digitIndex : labelsFromCRU->at(c)) { diff --git a/Detectors/TPC/reconstruction/src/RawReader.cxx b/Detectors/TPC/reconstruction/src/RawReader.cxx index 70a0eb06ccf44..326dccbb42684 100644 --- a/Detectors/TPC/reconstruction/src/RawReader.cxx +++ b/Detectors/TPC/reconstruction/src/RawReader.cxx @@ -27,13 +27,13 @@ using namespace o2::TPC; RawReader::RawReader(int region, int link, int run) - : mRegion(region) + : mUseRawInMode3(true) + , mApplyChannelMask(false) + , mCheckAdcClock(false) + , mRegion(region) , mLink(link) , mRun(run) , mLastEvent(-1) - , mUseRawInMode3(true) - , mApplyChannelMask(false) - , mCheckAdcClock(false) , mTimestampOfFirstData({0,0,0,0,0}) , mEvents() , mData() From 444caa1473cdce5648058e3fb70b7e68c71afd10 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Fri, 17 Nov 2017 13:28:50 +0100 Subject: [PATCH 09/12] add MC label sorting --- .../TPC/reconstruction/macro/addInclude.C | 2 +- .../reconstruction/macro/readClusterMCtruth.C | 10 ++++++++- .../TPC/reconstruction/src/HwClusterer.cxx | 22 +++++++++++++++++-- 3 files changed, 30 insertions(+), 4 deletions(-) diff --git a/Detectors/TPC/reconstruction/macro/addInclude.C b/Detectors/TPC/reconstruction/macro/addInclude.C index 490cb1990ae88..ce2e0ad8daa86 100644 --- a/Detectors/TPC/reconstruction/macro/addInclude.C +++ b/Detectors/TPC/reconstruction/macro/addInclude.C @@ -10,5 +10,5 @@ void addInclude() { - gSystem->AddIncludePath("-I$O2_ROOT/include -I$FAIRROOT_ROOT/include -I$VC_ROOT/include -I$BOOST_ROOT/include"); + gSystem->AddIncludePath("-I$O2_ROOT/include -I$FAIRROOT_ROOT/include -I$VC_ROOT/include -I$BOOST_ROOT/include -I$MS_GSL_ROOT/include"); } diff --git a/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C index 2eb01afa416c6..529d728893538 100644 --- a/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C +++ b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C @@ -12,8 +12,16 @@ /// \brief This macro demonstrates how to extract the MC truth information from /// the clusters /// \author Andi Mathis, TU München, andreas.mathis@ph.tum.de -#include +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "TFile.h" +#include "TTree.h" + +#include "TPCReconstruction/Cluster.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include +#endif void readClusterMCtruth(std::string filename) { TFile *clusterFile = TFile::Open(filename.data()); diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index b711bc27ca406..4506fbec95c74 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -433,18 +433,36 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta /* * collect clusters from individual cluster finder */ + + // map to count unique MC labels + std::map labelCount; + + // multiset to sort labels according to occurrence + auto mcComp = [](const std::pair& a, const std::pair& b) { return a.second > b.second;}; + std::multiset,decltype(mcComp)> labelSort(mcComp); + + // for each CRU for (unsigned cru = 0; cru < mClusterStorage.size(); ++cru) { std::vector* clustersFromCRU = &mClusterStorage[cru]; std::vector>>* labelsFromCRU = &mClusterDigitIndexStorage[cru]; + // for each found cluster for(unsigned c = 0; c < clustersFromCRU->size(); ++c) { const auto clusterPos = mClusterArray->size(); mClusterArray->emplace_back(clustersFromCRU->at(c)); + labelCount.clear(); + labelSort.clear(); + + // for each used digit for (auto &digitIndex : labelsFromCRU->at(c)) { if (digitIndex.first < 0) continue; - for (auto &l : mLastMcDigitTruth[digitIndex.second]->getLabels(digitIndex.first)) - mClusterMcLabelArray->addElement(clusterPos,l); + for (auto &l : mLastMcDigitTruth[digitIndex.second]->getLabels(digitIndex.first)) { + try { labelCount.at(l) = labelCount.at(l)+1;} // increment counter + catch (std::out_of_range &e) { labelCount[l] = 1;} // label didn't exist yet, set to 1 + } } + for (auto &l : labelCount) labelSort.insert(l); + for (auto &l : labelSort) mClusterMcLabelArray->addElement(clusterPos,l.first); } } From f3a0f20fa8eefbdaa57ff1bbbe972842e7f77271 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Fri, 17 Nov 2017 14:02:11 +0100 Subject: [PATCH 10/12] fix macro compilation error also wrong include paths, introduced by to moving Cluster related classes into reconstruction --- .../include/TPCReconstruction/HwClusterer.h | 12 ++++++++---- .../TPC/reconstruction/macro/RawClusterFinder.C | 13 ++++++------- .../reconstruction/macro/convertClusters.C_backup | 2 +- Detectors/TPC/reconstruction/macro/convertTracks.C | 2 +- Detectors/TPC/reconstruction/macro/dEdxRes.C | 2 +- .../TPC/reconstruction/macro/testClustererData.C | 9 ++++----- Detectors/TPC/reconstruction/macro/testRawRead.C | 8 ++++++++ Detectors/TPC/reconstruction/macro/testTracks.C | 3 +-- Detectors/TPC/reconstruction/src/ClustererTask.cxx | 2 +- Detectors/TPC/reconstruction/src/HwClusterer.cxx | 7 ++++--- .../TPC/simulation/macro/testHitsDigitsClusters.C | 3 +-- 11 files changed, 36 insertions(+), 27 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index ca4119e17fbb5..c70a9738b3dca 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -44,7 +44,7 @@ class HwClusterer : public Clusterer { /// Constructor /// \param clusterOutput is pointer to vector to be filled with clusters - /// \param labelOutput is reference to storage to be filled with MC labels + /// \param labelOutput is pointer to storage to be filled with MC labels /// \param cru Number of CRUs to process /// \param minQDiff Min charge difference /// \param assignChargeUnique Avoid using same charge for multiple nearby clusters @@ -54,7 +54,7 @@ class HwClusterer : public Clusterer { /// \param timebinsPerCF Time bins per cluster finder HwClusterer( std::vector *clusterOutput, - std::unique_ptr &labelOutput, + MCLabelContainer *labelOutput, int cruMin = 0, int cruMax = 359, float minQDiff = 0, @@ -84,6 +84,10 @@ class HwClusterer : public Clusterer { /// Setter for pedestal object, pedestal value will be subtracted before cluster finding /// \param pedestalObject CalDet object, containing pedestals for each pad void setPedestalObject(std::shared_ptr> pedestalObject) { mPedestalObject = pedestalObject; }; + void setPedestalObject(CalDet* pedestalObject) { + LOG(DEBUG) << "Consider using std::shared_ptr for the pedestal object." << FairLogger::endl; + mPedestalObject = std::shared_ptr>(pedestalObject); + }; /// Switch for triggered / continuous readout /// \param isContinuous - false for triggered readout, true for continuous readout @@ -162,8 +166,8 @@ class HwClusterer : public Clusterer { std::vector> mClusterStorage; ///< Cluster storage for each CRU std::vector>>> mClusterDigitIndexStorage; ///< Container for digit indices, used in found clusters. Pair consists of original digit index and event count - std::vector* mClusterArray; ///< Pointer to output cluster storage - std::unique_ptr& mClusterMcLabelArray; ///< Internal MC Label storage + std::vector* mClusterArray; ///< Pointer to output cluster storage + MCLabelContainer* mClusterMcLabelArray; ///< Pointer to MC Label storage std::shared_ptr> mNoiseObject; ///< Pointer to the CalDet object for noise simulation std::shared_ptr> mPedestalObject; ///< Pointer to the CalDet object for the pedestal subtraction diff --git a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C index be85f5c36b677..436da7ce6e0b6 100644 --- a/Detectors/TPC/reconstruction/macro/RawClusterFinder.C +++ b/Detectors/TPC/reconstruction/macro/RawClusterFinder.C @@ -29,13 +29,13 @@ #include "TPCBase/CalDet.h" #include "TPCBase/CRU.h" #include "TPCCalibration/CalibRawBase.h" -#include "TPCSimulation/HwClusterer.h" -#include "TPCSimulation/BoxClusterer.h" -#include "TPCSimulation/ClusterContainer.h" +#include "TPCReconstruction/HwClusterer.h" +#include "TPCReconstruction/BoxClusterer.h" +#include "TPCReconstruction/ClusterContainer.h" #include "TPCBase/Digit.h" #include "TPCBase/Mapper.h" #include "TPCReconstruction/TrackTPC.h" -#include "TPCSimulation/Cluster.h" +#include "TPCReconstruction/Cluster.h" #include "TCanvas.h" #endif @@ -176,7 +176,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt // HW cluster finder std::unique_ptr cl; if (clustererType == ClustererType::HW) { - HwClusterer *hwCl = new HwClusterer(&arrCluster, 0, 3); + HwClusterer *hwCl = new HwClusterer(&arrCluster, nullptr, 0, 3); hwCl->setContinuousReadout(false); hwCl->setPedestalObject(pedestal); cl = std::unique_ptr(hwCl); @@ -191,7 +191,6 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt } //cl->setRequirePositiveCharge(false); - cl->Init(); // Box cluster finder @@ -212,7 +211,7 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt if (!arr.size()) {++events; continue;} //printf("Converted digits: %zu %f\n", arr.size(), arr.at(0)->getChargeFloat()); - cl->Process(arr); + cl->Process(arr,nullptr,events); // ---| set cherenkov value|--------------------------------------------- if (istr.is_open()) { diff --git a/Detectors/TPC/reconstruction/macro/convertClusters.C_backup b/Detectors/TPC/reconstruction/macro/convertClusters.C_backup index 98fadf2c1add7..1532a9ac963d6 100644 --- a/Detectors/TPC/reconstruction/macro/convertClusters.C_backup +++ b/Detectors/TPC/reconstruction/macro/convertClusters.C_backup @@ -23,7 +23,7 @@ #include "TPCBase/Sector.h" #include "TPCBase/Mapper.h" #include "TPCBase/PadRegionInfo.h" -#include "TPCSimulation/Cluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCSimulation/Constants.h" /* diff --git a/Detectors/TPC/reconstruction/macro/convertTracks.C b/Detectors/TPC/reconstruction/macro/convertTracks.C index b28294f3ebc0a..3b12ebcbb6cf5 100644 --- a/Detectors/TPC/reconstruction/macro/convertTracks.C +++ b/Detectors/TPC/reconstruction/macro/convertTracks.C @@ -22,7 +22,7 @@ #include "TPCReconstruction/TrackTPC.h" #include "DetectorsBase/Track.h" -#include "TPCSimulation/Cluster.h" +#include "TPCReconstruction/Cluster.h" #endif using namespace o2::TPC; diff --git a/Detectors/TPC/reconstruction/macro/dEdxRes.C b/Detectors/TPC/reconstruction/macro/dEdxRes.C index 3ad379a1b0386..1b9021218aae7 100644 --- a/Detectors/TPC/reconstruction/macro/dEdxRes.C +++ b/Detectors/TPC/reconstruction/macro/dEdxRes.C @@ -9,7 +9,7 @@ // or submit itself to any jurisdiction. #if !defined(__CLING__) || defined(__ROOTCLING__) -#include "TPCSimulation/Cluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCBase/Mapper.h" #include "TPCReconstruction/TrackTPC.h" #include "TPCBase/CRU.h" diff --git a/Detectors/TPC/reconstruction/macro/testClustererData.C b/Detectors/TPC/reconstruction/macro/testClustererData.C index ec7305140ca76..54fdd6675c58f 100644 --- a/Detectors/TPC/reconstruction/macro/testClustererData.C +++ b/Detectors/TPC/reconstruction/macro/testClustererData.C @@ -17,8 +17,8 @@ #include "TClonesArray.h" #include "TString.h" -#include "TPCSimulation/HwClusterer.h" -#include "TPCSimulation/ClusterContainer.h" +#include "TPCReconstruction/HwClusterer.h" +#include "TPCReconstruction/ClusterContainer.h" #include "TPCBase/Digit.h" #include "TPCReconstruction/GBTFrameContainer.h" #endif @@ -48,9 +48,8 @@ void testClustererData(Int_t maxEvents=50, TString fileInfo="GBTx0_Run005:0:0;GB t.Branch("cl", &arrCluster); // ===| cluster finder |====================================================== - HwClusterer cl(&arrCluster); + HwClusterer cl(&arrCluster, nullptr); cl.setPedestalObject(pedestal); - cl.Init(); // ===| loop over all data |================================================== int events = 0; @@ -84,7 +83,7 @@ void testClustererData(Int_t maxEvents=50, TString fileInfo="GBTx0_Run005:0:0;GB //} //printf("\n"); - cl.Process(arr); + cl.Process(arr,nullptr,events); t.Fill(); printf("Found clusters: %lu\n", arrCluster.size()); diff --git a/Detectors/TPC/reconstruction/macro/testRawRead.C b/Detectors/TPC/reconstruction/macro/testRawRead.C index 877c1ac0e54fa..de6c7a4a78425 100644 --- a/Detectors/TPC/reconstruction/macro/testRawRead.C +++ b/Detectors/TPC/reconstruction/macro/testRawRead.C @@ -8,6 +8,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "TH1F.h" +#include "TH2F.h" +#include "TCanvas.h" + +#include "TPCReconstruction/GBTFrameContainer.h" +#endif + using namespace o2::TPC; void testRawRead(std::string filename) { diff --git a/Detectors/TPC/reconstruction/macro/testTracks.C b/Detectors/TPC/reconstruction/macro/testTracks.C index ca7969d0f4c0d..ca70d593b999e 100644 --- a/Detectors/TPC/reconstruction/macro/testTracks.C +++ b/Detectors/TPC/reconstruction/macro/testTracks.C @@ -21,8 +21,7 @@ #include "TPCSimulation/Digitizer.h" #include "TPCReconstruction/TrackTPC.h" #include "DetectorsBase/Track.h" -#include "TPCSimulation/Cluster.h" -#include "TPCSimulation/Cluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCBase/Mapper.h" #endif diff --git a/Detectors/TPC/reconstruction/src/ClustererTask.cxx b/Detectors/TPC/reconstruction/src/ClustererTask.cxx index 8309e5595baf3..da75496119a15 100644 --- a/Detectors/TPC/reconstruction/src/ClustererTask.cxx +++ b/Detectors/TPC/reconstruction/src/ClustererTask.cxx @@ -106,7 +106,7 @@ InitStatus ClustererTask::Init() mgr->Register("TPCClusterHWMCTruth", "TPC", mHwClustersMCTruthArray.get(), kTRUE); // create clusterer and pass output pointer - mHwClusterer = std::make_unique(mHwClustersArray,mHwClustersMCTruthArray);//,0,359); + mHwClusterer = std::make_unique(mHwClustersArray,mHwClustersMCTruthArray.get());//,0,359); mHwClusterer->setContinuousReadout(mIsContinuousReadout); // TODO: implement noise/pedestal objecta // mHwClusterer->setNoiseObject(); diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 4506fbec95c74..4185d259ca1f8 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -31,7 +31,7 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterer::HwClusterer(std::vector *clusterOutput, - std::unique_ptr &labelOutput, int cruMin, int cruMax, + MCLabelContainer *labelOutput, int cruMin, int cruMax, float minQDiff, bool assignChargeUnique, bool enableNoiseSim, bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF) : Clusterer() @@ -256,7 +256,7 @@ void HwClusterer::processDigits( void HwClusterer::Process(std::vector const &digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); - mClusterMcLabelArray->clear(); + if(mClusterMcLabelArray) mClusterMcLabelArray->clear(); /* @@ -324,7 +324,7 @@ void HwClusterer::Process(std::vector const &digits, MCLabelCont void HwClusterer::Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); - mClusterMcLabelArray->clear(); + if(mClusterMcLabelArray) mClusterMcLabelArray->clear(); /* * clear old storages @@ -450,6 +450,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta for(unsigned c = 0; c < clustersFromCRU->size(); ++c) { const auto clusterPos = mClusterArray->size(); mClusterArray->emplace_back(clustersFromCRU->at(c)); + if (mClusterMcLabelArray == nullptr) continue; labelCount.clear(); labelSort.clear(); diff --git a/Detectors/TPC/simulation/macro/testHitsDigitsClusters.C b/Detectors/TPC/simulation/macro/testHitsDigitsClusters.C index d2683a1c11891..712c929af4d77 100644 --- a/Detectors/TPC/simulation/macro/testHitsDigitsClusters.C +++ b/Detectors/TPC/simulation/macro/testHitsDigitsClusters.C @@ -30,7 +30,6 @@ #include "TPCReconstruction/TrackTPC.h" #include "DetectorsBase/Track.h" #include "TPCReconstruction/Cluster.h" -#include "TPCReconstruction/HwCluster.h" #include "TPCBase/Mapper.h" #endif @@ -156,7 +155,7 @@ void testHitsDigitsClusters(int iEv=0, TFile *clusterFile = TFile::Open(clusFile.data()); TTree *clusterTree = (TTree *)gDirectory->Get("o2sim"); - std::vector *clustersArray = nullptr; + std::vector *clustersArray = nullptr; clusterTree->SetBranchAddress("TPCClusterHW", &clustersArray); TGraph *grClustersA = new TGraph(); From e9b736090b9c3f27394bb9ee9e3abb4fc5a0a42c Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Fri, 17 Nov 2017 15:58:55 +0100 Subject: [PATCH 11/12] removed try/catch construct --- .../include/TPCReconstruction/RawReaderEventSync.h | 12 +++--------- Detectors/TPC/reconstruction/src/HwClusterer.cxx | 3 +-- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReaderEventSync.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReaderEventSync.h index 255af905d5a44..99fa280de8d84 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReaderEventSync.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/RawReaderEventSync.h @@ -68,15 +68,9 @@ inline void RawReaderEventSync::addEventOffset(int event, long time) { // set time offset to next multiple of mFrames (+ mFrames) - long iTimeOffset = time + (mFrames - (time%mFrames)) + mFrames; - - try { - // using at() function to provoke std::out_of_range exception if "event" does not yet exist in map - mMaxOffset.at(event) = std::max(mMaxOffset.at(event),iTimeOffset); - } catch (const std::out_of_range& e) { - // using [] operator to insert new element - mMaxOffset[event] = iTimeOffset; // First element, therefore has to be the maximum and use it - } + const long iTimeOffset = time + (mFrames - (time%mFrames)) + mFrames; + + mMaxOffset[event] = std::max(mMaxOffset[event],iTimeOffset); } } diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index 4185d259ca1f8..a1722ec78140a 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -458,8 +458,7 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelConta for (auto &digitIndex : labelsFromCRU->at(c)) { if (digitIndex.first < 0) continue; for (auto &l : mLastMcDigitTruth[digitIndex.second]->getLabels(digitIndex.first)) { - try { labelCount.at(l) = labelCount.at(l)+1;} // increment counter - catch (std::out_of_range &e) { labelCount[l] = 1;} // label didn't exist yet, set to 1 + labelCount[l]++; } } for (auto &l : labelCount) labelSort.insert(l); From 92d2d3f1dadabe647d3b0a01fb66d9e6bbb76b82 Mon Sep 17 00:00:00 2001 From: Sebastian Klewin Date: Mon, 20 Nov 2017 15:01:42 +0100 Subject: [PATCH 12/12] implemented Matthias' comments + fix padMean --- .../include/TPCReconstruction/BoxClusterer.h | 17 ++++++++++++++++- .../include/TPCReconstruction/Clusterer.h | 11 ++++++++--- .../include/TPCReconstruction/ClustererTask.h | 2 +- .../include/TPCReconstruction/HwClusterFinder.h | 6 +++--- .../include/TPCReconstruction/HwClusterer.h | 16 ++++++++++++++-- .../TPC/reconstruction/src/BoxClusterer.cxx | 6 ++++-- .../TPC/reconstruction/src/HwClusterFinder.cxx | 4 ++-- .../TPC/reconstruction/src/HwClusterer.cxx | 10 ++++++---- 8 files changed, 54 insertions(+), 18 deletions(-) diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h index c5b0e4f487ae8..8701300c8cd40 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h @@ -32,7 +32,22 @@ namespace o2{ class BoxClusterer : public Clusterer { public: - BoxClusterer(std::vector *output); + + /// Constructor + /// \param output is pointer to vector to be filled with clusters + /// \param rowsMax Max number of rows to process + /// \param padsMax Max number of pads to process + /// \param timeBinsMax Max number of timebins to process + /// \param minQMax Minimum peak charge for cluster + /// \param requirePositiveCharge Positive charge is required + /// \param requireNeighbouringPad Requires at least 2 adjecent pads with charge above threshold + BoxClusterer(std::vector *output, + int rowsMax = 18, + int padsMax = 138, + int timeBinsMax = 1024, + int minQMax = 5, + bool requirePositiveCharge = true, + bool requireNeighbouringPad = true); /// Destructor ~BoxClusterer() override; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h index 8b23f4dae72d6..5c9f901d59a95 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h @@ -34,7 +34,7 @@ class Clusterer { public: /// Default Constructor - Clusterer() : Clusterer(18, 138, 1024, 5, true, true) {}; + Clusterer() = delete; /// Constructor /// \param rowsMax Max number of rows to process @@ -43,8 +43,13 @@ class Clusterer { /// \param minQMax Minimum peak charge for cluster /// \param requirePositiveCharge Positive charge is required /// \param requireNeighbouringPad Requires at least 2 adjecent pads with charge above threshold - Clusterer(int rowsMax, int padsMax, int timeBinsMax, int minQMax, - bool requirePositiveCharge, bool requireNeighbouringPad); + Clusterer( + int rowsMax = 18, + int padsMax = 138, + int timeBinsMax = 1024, + int minQMax = 5, + bool requirePositiveCharge = true, + bool requireNeighbouringPad = true); /// Destructor virtual ~Clusterer() = default; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index 6206bf1df0c4c..7fc7bce734315 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -78,7 +78,7 @@ class ClustererTask : public FairTask{ /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() BoxClusterer* getBoxClusterer() const { return mBoxClusterer.get(); }; - /// Returns pointer to Box Clusterer + /// Returns pointer to Hw Clusterer /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() HwClusterer* getHwClusterer() const { return mHwClusterer.get(); }; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index 3a3a00728b978..f3446187d8ac0 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -51,7 +51,7 @@ class HwClusterFinder { /// \param chargeThreshold Minimum charge of cluster peak /// \param requirePositiveCharge Charge >0 required HwClusterFinder(unsigned short cru, unsigned short row, - unsigned short padOffset, unsigned short pads=8, unsigned short timebins=8, + short padOffset, unsigned short pads=8, unsigned short timebins=8, float diffThreshold=0, float chargeThreshold=5, bool requirePositiveCharge=true); /// Destructor @@ -95,7 +95,7 @@ class HwClusterFinder { /// Getter function /// \return Configured pad offset of this CF - unsigned short getPadOffset() const { return mPadOffset; } + short getPadOffset() const { return mPadOffset; } /// Getter function /// \return Width of CF in pad direction @@ -158,9 +158,9 @@ class HwClusterFinder { bool mRequireNeighbouringPad; ///< Switch if at least one neighboring pad needs charge > 0 bool mRequireNeighbouringTimebin; ///< Switch if at least one neighboring time bin needs charge > 0 bool mAssignChargeUnique; ///< Switch for "charge splitting", TODO: not yet properly implemented + short mPadOffset; ///< Pad number in row of leftmost pad (can be negative for leftmost CF) unsigned short mCRU; ///< CRU number unsigned short mRow; ///< Row number - unsigned short mPadOffset; ///< Pad number in row of leftmost pad unsigned short mPads; ///< Size of CF in pad direction unsigned short mTimebins; ///< Size of CF in time direction unsigned short mClusterSizePads; ///< Size of cluster in pad direction diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index c70a9738b3dca..7c75d2bd918f5 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -52,17 +52,29 @@ class HwClusterer : public Clusterer { /// \param enablePedestalSubtraction Enables the Pedestal subtraction (pedestal object has to be set) /// \param padsPerCF Pads per cluster finder /// \param timebinsPerCF Time bins per cluster finder + /// \param rowsMax Max number of rows to process + /// \param padsMax Max number of pads to process + /// \param timeBinsMax Max number of timebins to process + /// \param minQMax Minimum peak charge for cluster + /// \param requirePositiveCharge Positive charge is required + /// \param requireNeighbouringPad Requires at least 2 adjecent pads with charge above threshold HwClusterer( std::vector *clusterOutput, MCLabelContainer *labelOutput, int cruMin = 0, int cruMax = 359, float minQDiff = 0, - bool assignChargeUnique = false, + bool assignChargeUnique = true,//false, bool enableNoiseSim = true, bool enablePedestalSubtraction = true, int padsPerCF = 8, - int timebinsPerCF = 8); + int timebinsPerCF = 8, + int rowsMax = 18, + int padsMax = 138, + int timeBinsMax = 1024, + int minQMax = 5, + bool requirePositiveCharge = true, + bool requireNeighbouringPad = true); /// Destructor ~HwClusterer() override; diff --git a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx index 7f043a38f451a..d8e1dcf65891e 100644 --- a/Detectors/TPC/reconstruction/src/BoxClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx @@ -103,8 +103,10 @@ ClassImp(o2::TPC::BoxClusterer) using namespace o2::TPC; //________________________________________________________________________ -BoxClusterer::BoxClusterer(std::vector *output): - Clusterer(), +BoxClusterer::BoxClusterer(std::vector *output, + int rowsMax, int padsMax, int timeBinsMax, int minQMax, + bool requirePositiveCharge, bool requireNeighbouringPad): + Clusterer(rowsMax,padsMax,timeBinsMax,minQMax,requirePositiveCharge,requireNeighbouringPad), mClusterArray(output), mAllBins(nullptr), mAllSigBins(nullptr), diff --git a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx index 65a36b56a46ff..a11545963f397 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterFinder.cxx @@ -25,16 +25,16 @@ using namespace o2::TPC; //________________________________________________________________________ HwClusterFinder::HwClusterFinder( unsigned short cru, unsigned short row, - unsigned short padOffset, unsigned short pads, unsigned short timebins, + short padOffset, unsigned short pads, unsigned short timebins, float diffThreshold, float chargeThreshold, bool requirePositiveCharge) : mRequirePositiveCharge(requirePositiveCharge) , mRequireNeighbouringPad(false)//true) , mRequireNeighbouringTimebin(true) , mAssignChargeUnique(false)//true) + , mPadOffset(padOffset) , mCRU(cru) , mRow(row) - , mPadOffset(padOffset) , mPads(pads) , mTimebins(timebins) , mClusterSizePads(5) diff --git a/Detectors/TPC/reconstruction/src/HwClusterer.cxx b/Detectors/TPC/reconstruction/src/HwClusterer.cxx index a1722ec78140a..0bd11fb1ff40e 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -33,8 +33,10 @@ using namespace o2::TPC; HwClusterer::HwClusterer(std::vector *clusterOutput, MCLabelContainer *labelOutput, int cruMin, int cruMax, float minQDiff, bool assignChargeUnique, bool enableNoiseSim, - bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF) - : Clusterer() + bool enablePedestalSubtraction, int padsPerCF, int timebinsPerCF, + int rowsMax, int padsMax, int timeBinsMax, int minQMax, + bool requirePositiveCharge, bool requireNeighbouringPad) + : Clusterer(rowsMax,padsMax,timeBinsMax,minQMax,requirePositiveCharge,requireNeighbouringPad) , mAssignChargeUnique(assignChargeUnique) , mEnableNoiseSim(enableNoiseSim) , mEnablePedestalSubtraction(enablePedestalSubtraction) @@ -59,7 +61,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, /* * initialize all cluster finder */ - unsigned iCfPerRow = (unsigned)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); + unsigned iCfPerRow = (unsigned)ceil((double)(mPadsMax+2+2)/(static_cast(mPadsPerCF)-2-2)); mClusterFinder.resize(mCRUMax+1); const Mapper& mapper = Mapper::instance(); for (unsigned iCRU = mCRUMin; iCRU <= mCRUMax; ++iCRU){ @@ -67,7 +69,7 @@ HwClusterer::HwClusterer(std::vector *clusterOutput, for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); ++iRow){ mClusterFinder[iCRU][iRow].resize(iCfPerRow); for (unsigned iCF = 0; iCF < iCfPerRow; ++iCF){ - int padOffset = iCF*(mPadsPerCF-2-2)-2; + int padOffset = iCF*(static_cast(mPadsPerCF)-2-2)-2; mClusterFinder[iCRU][iRow][iCF] = std::make_shared(iCRU,iRow,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); mClusterFinder[iCRU][iRow][iCF]->setAssignChargeUnique(mAssignChargeUnique);