diff --git a/Detectors/TPC/reconstruction/CMakeLists.txt b/Detectors/TPC/reconstruction/CMakeLists.txt index eef0425f832ac..27c2808652fec 100644 --- a/Detectors/TPC/reconstruction/CMakeLists.txt +++ b/Detectors/TPC/reconstruction/CMakeLists.txt @@ -6,8 +6,6 @@ link_directories( ${LINK_DIRECTORIES}) set(SRCS src/AdcClockMonitor.cxx - - src/BoxCluster.cxx src/BoxClusterer.cxx src/Cluster.cxx src/Clusterer.cxx @@ -16,7 +14,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 +28,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 +36,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/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/BoxClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h index d7d6d8fd440b8..8701300c8cd40 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/BoxClusterer.h @@ -18,32 +18,48 @@ #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); - + + /// 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; - - // 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 eventCount event counter /// @return Container with clusters - void Process(std::vector const & digits) override; - void Process(std::vector>& digits) 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; } @@ -51,21 +67,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 // @@ -74,12 +90,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/Clusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h index b56771755a203..5c9f901d59a95 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/Clusterer.h @@ -17,20 +17,24 @@ #include #include -#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 MCLabelContainer = o2::dataformats::MCTruthContainer; + public: /// Default Constructor - Clusterer(); + Clusterer() = delete; /// Constructor /// \param rowsMax Max number of rows to process @@ -39,20 +43,24 @@ 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; - - /// Initialization function for clusterer - virtual void Init() = 0; - + /// Processing all digits /// \param digits Container with TPC digits + /// \param mcDigitTruth MC Digit Truth container + /// \param eventCount event counter /// \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, 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; }; @@ -67,20 +75,31 @@ 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 - - }; + +}; + +//________________________________________________________________________ +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) +{} + } } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h index 54a08b04a1343..7fc7bce734315 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h @@ -18,26 +18,35 @@ #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; - + InitStatus Init() override; 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; @@ -45,41 +54,57 @@ class ClustererTask : public FairTask{ }; }; - bool isClustererEnable(ClustererType type) const { + /// 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; case ClustererType::Box: return mBoxClustererEnable; }; }; - Clusterer* getClusterer(ClustererType type) { + /// 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(); }; }; + + /// Returns pointer to Box Clusterer + /// \return Pointer to Clusterer, nullptr if Clusterer was not enabled during Init() + BoxClusterer* getBoxClusterer() const { return mBoxClusterer.get(); }; - BoxClusterer* getBoxClusterer() const { return mBoxClusterer; }; - HwClusterer* getHwClusterer() const { return mHwClusterer; }; - // Clusterer *GetClusterer() const { return fClusterer; } - + /// Returns pointer to Hw 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 + 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 - BoxClusterer *mBoxClusterer; - HwClusterer *mHwClusterer; - - std::vector const *mDigitsArray; - // produced data containers - std::vector *mClustersArray; - std::vector *mHwClustersArray; - ClassDefOverride(ClustererTask, 1) }; diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwCluster.h deleted file mode 100644 index 7100df03c1540..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/include/TPCReconstruction/HwClusterFinder.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h index 6ab8aac820712..f3446187d8ac0 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterFinder.h @@ -15,106 +15,120 @@ #define ALICEO2_TPC_HWClusterFinder_H_ #include +#include +#include +#include #include namespace o2{ namespace TPC { -class HwCluster; - +class Cluster; + /// \class HwClusterFinder /// \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; + 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 - /// \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, - short padOffset, short pads=8, short timebins=8, + HwClusterFinder(unsigned short cru, unsigned short row, + short padOffset, unsigned short pads=8, unsigned short timebins=8, float diffThreshold=0, float chargeThreshold=5, bool requirePositiveCharge=true); /// Destructor - ~HwClusterFinder(); + ~HwClusterFinder() = default; - /// 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(float* 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 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; } - - // 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 + unsigned short getCRU() const { return mCRU; } + + /// Getter function + /// \return Configured row number + unsigned 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 + unsigned short getNpads() const { return mPads; } + + /// Getter function + /// \return Width of CF in time direction + unsigned 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(); } + void clearClusterContainer() { mClusterContainer.clear(); mClusterDigitIndices.clear(); } /// Process the cluster finding bool findCluster(); @@ -122,69 +136,72 @@ 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, float** cluster); + void clusterAlreadyUsed(short time, short pad); private: - float chargeForCluster(float* charge, float* 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; - int mTimebinsAfterLastProcessing; - float** mData; - float** tmpCluster; - float* 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 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 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 + 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(float* 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 // - float* 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; } } } -#endif +#endif diff --git a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h index 37a0fd491c57c..7c75d2bd918f5 100644 --- a/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h +++ b/Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h @@ -15,135 +15,179 @@ #define ALICEO2_TPC_HWClusterer_H_ #include "TPCReconstruction/Clusterer.h" -#include "TPCReconstruction/HwCluster.h" -#include "TPCBase/CalDet.h" +#include "TPCReconstruction/Cluster.h" +#include "TPCBase/CalDet.h" + +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" #include +#include +#include +#include +#include namespace o2{ namespace TPC { - -class ClusterContainer; + class ClustererTask; class HwClusterFinder; -class HwCluster; class Digit; /// \class HwClusterer /// \brief Class for TPC HW cluster finding class HwClusterer : public Clusterer { + + 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 processingType parallel or sequential - /// \param globalTime value of first timebin + /// \param clusterOutput is pointer to vector to be filled with clusters + /// \param labelOutput is pointer to storage to be filled with MC labels /// \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 + /// \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 *output, - Processing processingType = Processing::Parallel, - int globalTime = 0, + std::vector *clusterOutput, + MCLabelContainer *labelOutput, 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); - + bool assignChargeUnique = true,//false, + bool enableNoiseSim = true, + bool enablePedestalSubtraction = true, + int padsPerCF = 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; - - // 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 eventCount event counter /// @return Container with clusters - void Process(std::vector const &digits) override; - void Process(std::vector>& digits) override; - - void setProcessingType(Processing processing) { mProcessingType = processing; }; - - void setNoiseObject(CalDet* noiseObject) { mNoiseObject = noiseObject; }; - void setPedestalObject(CalDet* pedestalObject) { mPedestalObject = pedestalObject; }; + void Process(std::vector const &digits, + MCLabelContainer const* mcDigitTruth, int eventCount) override; + void Process(std::vector>& digits, + MCLabelContainer const* mcDigitTruth, int eventCount) override; + + /// 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; }; + 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 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; }; - + + /// Set number of parallel threads + /// \param threads Number to be set, if 0 hardware default value is used + void setNumThreads(unsigned threads); + 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; + 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 + 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 + 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, - std::vector& cluster, + const std::vector>>>& digits, + const std::vector>>>& clusterFinder, + std::vector>& cluster, + std::vector>>>& label, 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> 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 + */ + 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 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 + 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 + 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 + + std::map> mLastMcDigitTruth; ///< Buffer for digit MC truth information }; } } -#endif +#endif 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/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/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 9d7e335504728..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 @@ -143,8 +143,8 @@ void RawClusterFinder::processEvents(TString fileInfo, TString pedestalFile, TSt } // ===| output file and container |=========================================== - std::vector arrCluster; - std::vector *arrClusterBox = nullptr; + std::vector arrCluster; + std::vector *arrClusterBox = nullptr; float cherenkovValue = 0.; int runNumber = 0; @@ -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, 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/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/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 95c0badb5635b..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; @@ -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/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/readClusterMCtruth.C b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C new file mode 100644 index 0000000000000..529d728893538 --- /dev/null +++ b/Detectors/TPC/reconstruction/macro/readClusterMCtruth.C @@ -0,0 +1,51 @@ +// 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 +#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()); + 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/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..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 @@ -42,15 +42,14 @@ 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); // ===| 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 33e88aa9cf185..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/HwCluster.h" +#include "TPCReconstruction/Cluster.h" #include "TPCBase/Mapper.h" #endif @@ -44,7 +43,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/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/BoxClusterer.cxx b/Detectors/TPC/reconstruction/src/BoxClusterer.cxx index e56477347ece2..d8e1dcf65891e 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,30 +103,15 @@ 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), 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]; @@ -145,7 +129,19 @@ void BoxClusterer::Init() } //________________________________________________________________________ -void BoxClusterer::Process(std::vector const &digits) +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 @@ -181,7 +177,7 @@ void BoxClusterer::Process(std::vector const &digits) } //________________________________________________________________________ -void BoxClusterer::Process(std::vector>& digits) +void BoxClusterer::Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); @@ -352,8 +348,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/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 184d7ddc151b2..da75496119a15 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,11 +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(nullptr) + , mHwClustersMCTruthArray(nullptr) { } @@ -45,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; + } //_____________________________________________________________________ @@ -74,26 +76,38 @@ 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; + mClustersArray = new std::vector; mgr->RegisterAny("TPCCluster", mClustersArray, 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) { // Register output container - mHwClustersArray = new std::vector; + mHwClustersArray = new std::vector; mgr->RegisterAny("TPCClusterHW", mHwClustersArray, 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.get());//,0,359); mHwClusterer->setContinuousReadout(mIsContinuousReadout); - mHwClusterer->Init(); // TODO: implement noise/pedestal objecta // mHwClusterer->setNoiseObject(); // mHwClusterer->setPedestalObject(); @@ -109,12 +123,16 @@ void ClustererTask::Exec(Option_t *option) if (mBoxClustererEnable) { mClustersArray->clear(); - mBoxClusterer->Process(*mDigitsArray); + if(mClustersMCTruthArray) mClustersMCTruthArray->clear(); + mBoxClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); + LOG(DEBUG) << "Box clusterer found " << mClustersArray->size() << " clusters" << FairLogger::endl << FairLogger::endl; } if (mHwClustererEnable) { - mHwClustersArray->clear(); - mHwClusterer->Process(*mDigitsArray); - LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl; + if(mHwClustersArray) mHwClustersArray->clear(); + if(mHwClustersMCTruthArray) mHwClustersMCTruthArray->clear(); + mHwClusterer->Process(*mDigitsArray,mDigitMCTruthArray,mEventCount); + LOG(DEBUG) << "Hw clusterer found " << mHwClustersArray->size() << " clusters" << FairLogger::endl << FairLogger::endl; } + ++mEventCount; } 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(5, MiniDigit())) + , mClusterContainer() + , mClusterDigitIndices() { - 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; + << " 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 float*[mTimebins]; - for (t = 0; t < mTimebins; ++t){ - mData[t] = new float[mPads]; - for (p = 0; p < mPads; ++p){ - mData[t][p] = 0; - } - } - - tmpCluster = new float*[mClusterSizeTime]; - for (t=0; t>(pads, MiniDigit())); -//________________________________________________________________________ -HwClusterFinder::~HwClusterFinder() -{ - short t; - for (t = 0; t < mTimebins; ++t){ - delete [] mData[t]; - } - delete [] mData; - delete [] mZeroTimebin; - - for (t=0; tbegin() /* some iterator, is not used */,globalTime,length,true); } //________________________________________________________________________ -void HwClusterFinder::AddZeroTimebin(unsigned globalTime, int length) +void HwClusterFinder::printLocalStorage() { - if (mZeroTimebin == nullptr) { - mZeroTimebin = new float[length]; - for (short i = 0; i < length; ++i) mZeroTimebin[i] = 0; + for (short t = 0; t < mTimebins; ++t){ + LOG(DEBUG) << "t " << t << ":\t"; + for (auto &digi : *mData[t]) + LOG(DEBUG) << digi.charge << "\t"; + LOG(DEBUG) << FairLogger::endl; } - bool ret = AddTimebin(mZeroTimebin,globalTime,length); -} - -//________________________________________________________________________ -void HwClusterFinder::PrintLocalStorage() -{ - short t,p; - for (t = 0; t < mTimebins; ++t){ - printf("t %d:\t",t); - for (p = 0; p < mPads; ++p){ - printf("%.2f\t", mData[t][p]); - } - printf("\n"); - } - printf("\n"); + LOG(DEBUG) << FairLogger::endl; } //________________________________________________________________________ @@ -189,25 +103,29 @@ 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 // 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; + float charge; + float meanP, meanT; + float sigmaP, sigmaT; + short minP, minT; + short maxP, maxT; + short deltaP, deltaT; // // peak finding // @@ -225,35 +143,35 @@ bool HwClusterFinder::findCluster() // o i i i o // o o o o o // - if (mData[t ][p ] < mChargeThreshold) continue; + if (mData[t ]->at(p).charge < mChargeThreshold) continue; // Require at least one neighboring time bin with signal - if (mRequireNeighbouringTimebin && (mData[t-1][p ] + mData[t+1][p ] <= 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] + mData[t ][p+1] <= 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 ] >= 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]->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); - + // // cluster was found!! // - + // 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; } } - + // // The outer cells of the 5x5 matrix (o) are taken only if the // neighboring inner cell (i) has a signal above threshold. // - + // // The cells of the "inner cross" have here only 1 neighbour. - // [t] - // 0 o - // 1 i - // 2 o i C i o - // 3 i - // 4 o + // [t] + // 0 o + // 1 i + // 2 o i C i o + // 3 i + // 4 o // // [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. // 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]); - 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); -// 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 = mTmpCluster[2][2].charge; + qTot = 0; + meanP = 0; + meanT = 0; + sigmaP = 0; + sigmaT = 0; + minP = mClusterSizePads; + minT = mClusterSizeTime; + maxP = 0; + maxT = 0; + mClusterDigitIndices.emplace_back(); + for (tt = 0; tt < mClusterSizeTime; ++tt) { + deltaT = tt - mClusterSizeTime/2; + for (pp = 0; pp < mClusterSizePads; ++pp) { + deltaP = pp - mClusterSizePads/2; + + 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 * static_cast(deltaP); + meanT += charge * static_cast(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); + minT = std::min(minT,tt); maxT = std::max(maxT,tt); + } + } + } + + 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); + mClusterContainer.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); + if (auto next = mNextCF.lock()) { + next->clusterAlreadyUsed(t,p+mPadOffset);//,mTmpCluster); + } } - + // // 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 -= mTmpCluster[tt][pp]; + mData[t+(tt-2)]->at(p+(pp-2)).clear(); } } } @@ -351,42 +318,19 @@ bool HwClusterFinder::findCluster() } //________________________________________________________________________ -float HwClusterFinder::chargeForCluster(float* charge, float* toCompare) +bool HwClusterFinder::chargeForCluster(float outerCharge, float innerCharge) { //printf("%.2f - %.2f = %.f compared to %.2f)?\n",toCompare,*charge,toCompare-*charge,-mDiffThreshold); - if ((mRequirePositiveCharge && (*charge > 0)) & - (*toCompare > mDiffThreshold)) {//mChargeThreshold)) { - //printf("\tyes\n"); - return *charge; - } else { - //printf("\tno\n"); - return 0; - } -} - -//________________________________________________________________________ -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; + if ((mRequirePositiveCharge && (outerCharge > 0)) && + (innerCharge > mChargeThreshold)) { + return true; } - mNextCF = nextCF; - - return; + return false; } //________________________________________________________________________ -void HwClusterFinder::clusterAlreadyUsed(short time, short pad, float** cluster) +// TODO: really nexessary?? or just set to 0 +void HwClusterFinder::clusterAlreadyUsed(short time, short pad) { short localPad = pad - mPadOffset; @@ -395,8 +339,9 @@ 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]; + + //mData[t][p].charge -= cluster[t-time+2][p-localPad+2].charge; + mData[t]->at(p).clear(); } } } @@ -404,12 +349,8 @@ void HwClusterFinder::clusterAlreadyUsed(short time, short pad, float** cluster) //________________________________________________________________________ void HwClusterFinder::reset(unsigned globalTimeAfterReset) { - short t,p; - for (t = 0; t < mTimebins; ++t){ - for (p = 0; p < mPads; ++p){ - mData[t][p] = 0; - } - } + for(auto &tb : mData) + for (auto &digi : *tb) digi.clear(); mGlobalTimeOfLast = globalTimeAfterReset; } @@ -419,10 +360,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]); + 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 0bc3491eedeb2..0bd11fb1ff40e 100644 --- a/Detectors/TPC/reconstruction/src/HwClusterer.cxx +++ b/Detectors/TPC/reconstruction/src/HwClusterer.cxx @@ -11,87 +11,71 @@ /// \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" #include "TPCBase/PadSecPos.h" -#include "TPCBase/CalArray.h" +#include "TPCBase/CalArray.h" #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) - : Clusterer() - , mClusterArray(output) - , mProcessingType(processingType) - , mGlobalTime(globalTime) - , mCRUMin(cruMin) - , mCRUMax(cruMax) - , mMinQDiff(minQDiff) +HwClusterer::HwClusterer(std::vector *clusterOutput, + MCLabelContainer *labelOutput, int cruMin, int cruMax, + float minQDiff, bool assignChargeUnique, bool enableNoiseSim, + 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) , mIsContinuousReadout(true) + , mLastTimebin(-1) + , mCRUMin(cruMin) + , mCRUMax(cruMax) , mPadsPerCF(padsPerCF) , mTimebinsPerCF(timebinsPerCF) - , mCfPerRow(cfPerRow) - , mLastTimebin(-1) + , mNumThreads(std::thread::hardware_concurrency()) + , mMinQDiff(minQDiff) + , mClusterFinder() + , mDigitContainer() + , mClusterStorage() + , mClusterDigitIndexStorage() + , mClusterArray(clusterOutput) + , mClusterMcLabelArray(labelOutput) , mNoiseObject(nullptr) , mPedestalObject(nullptr) + , mLastMcDigitTruth() { -} - -//________________________________________________________________________ -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 + * initialize all cluster finder */ - mCfPerRow = (int)ceil((double)(mPadsMax+2+2)/(mPadsPerCF-2-2)); - mClusterFinder.resize(mCRUMax); + unsigned iCfPerRow = (unsigned)ceil((double)(mPadsMax+2+2)/(static_cast(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++){ - mClusterFinder[iCRU][iRow].resize(mCfPerRow); - for (int iCF = 0; iCF < mCfPerRow; iCF++){ - int padOffset = iCF*(mPadsPerCF-2-2)-2; - mClusterFinder[iCRU][iRow][iCF] = new HwClusterFinder(iCRU,iRow,iCF,padOffset,mPadsPerCF,mTimebinsPerCF,mMinQDiff,mMinQMax,mRequirePositiveCharge); + for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); ++iRow){ + mClusterFinder[iCRU][iRow].resize(iCfPerRow); + for (unsigned iCF = 0; iCF < iCfPerRow; ++iCF){ + 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); /* - * 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. */ @@ -107,159 +91,162 @@ 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); - /* + /* * 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); - for (int iCRU = mCRUMin; iCRU < mCRUMax; iCRU++) + mDigitContainer.resize(mCRUMax+1); + for (unsigned 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, - std::vector& cluster, + const std::vector>>>& digits, + const std::vector>>>& clusterFinder, + std::vector>& cluster, + std::vector>>>& label, const CfConfig config) -// int iCRU, -// int maxRows, -// int maxPads, -// unsigned minTime, -// unsigned maxTime) { - 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(); -// float iAllBins[maxTime][maxPads]; -// for (float (&timebin)[maxPads] : iAllBins) { -// for (float &pad : timebin) { -// pad = 0.0; -// } -// } + int timeDiff = (config.iMaxTimeBin+1) - config.iMinTimeBin; + if (timeDiff < 0) return; + const Mapper& mapper = Mapper::instance(); + std::vector> iAllBins(timeDiff,std::vector(config.iMaxPads)); - float iAllBins[timeDiff][config.iMaxPads]; + for (unsigned iCRU = config.iCRUMin; iCRU <= config.iCRUMax; ++iCRU) { + if (iCRU % config.iThreadMax != config.iThreadID) continue; - for (int iRow = 0; iRow < config.iMaxRows; iRow++){ + for (int iRow = 0; iRow < mapper.getNumberOfRowsPartition(iCRU); iRow++){ - /* - * 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] = config.iNoiseObject->getValue(CRU(config.iCRU),iRow,p); + /* + * prepare local storage + */ + short t,p; + 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 { + 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()); } - } else { - std::fill(&iAllBins[0][0], &iAllBins[0][0]+timeDiff*config.iMaxPads, 0.f); - } - /* - * 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(); - -// std::cout << iCRU << " " << iRow << " " << iPad << " " << iTime << " (" << iTime-minTime << "," << timeDiff << ") " << charge << std::endl; - iAllBins[iTime-config.iMinTimeBin][iPad] += charge; - 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; - } - } - - /* - * 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){ - 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 + * 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 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; + 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); + clusterFinder[iCRU][iRow][cf]->addTimebin(iAllBins[time].begin()+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 (std::vector::const_reverse_iterator &cf_it : cfWithCluster) { - std::vector* cc = (*cf_it)->getClusterContainer(); - for (HwCluster& 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); + } + + (*cf_rit)->clearClusterContainer(); } - (*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,108 +255,157 @@ void HwClusterer::processDigits( } //________________________________________________________________________ -void HwClusterer::Process(std::vector const &digits) +void HwClusterer::Process(std::vector const &digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); + if(mClusterMcLabelArray) 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(); } -// int iCRU; -// int iRow; -// int iPad; -// float charge; int iTimeBin; int iTimeBinMin = (mIsContinuousReadout)?mLastTimebin + 1 : 0; //int iTimeBinMin = mLastTimebin + 1; int iTimeBinMax = mLastTimebin; - /* + if (mLastMcDigitTruth.find(eventCount) == mLastMcDigitTruth.end()) + mLastMcDigitTruth[eventCount] = std::make_unique(); + gsl::span mcArray; + /* * 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; + 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)); + ++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; + } + iTimeBinMax = std::max(iTimeBinMax,iTimeBin); - mDigitContainer[digit.getCRU()][digit.getRow()].push_back(&digit); + if (mcDigitTruth == nullptr) + 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); + + LOG(DEBUG) << "Event ranged from time bin " << iTimeBinMin << " to " << iTimeBinMax << "." << FairLogger::endl; } //________________________________________________________________________ -void HwClusterer::Process(std::vector>& digits) +void HwClusterer::Process(std::vector>& digits, MCLabelContainer const* mcDigitTruth, int eventCount) { mClusterArray->clear(); + if(mClusterMcLabelArray) 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(); + for (auto& dcc : dc) dcc.clear(); } int iTimeBin; 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 */ + 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; + 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)); + ++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; + } + iTimeBinMax = std::max(iTimeBinMax,iTimeBin); - mDigitContainer[digit->getCRU()][digit->getRow()].push_back(digit); + if (mcDigitTruth == nullptr) + 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); + + LOG(DEBUG) << "Event ranged from time bin " << iTimeBinMin << " to " << iTimeBinMax << "." << FairLogger::endl; } -void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) +void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax, MCLabelContainer const* mcDigitTruth, int eventCount) { /* - * 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), - mPadsMax+2+2, + threadId, + mNumThreads, + mCRUMin, + mCRUMax, + static_cast(mPadsMax)+2+2, iTimeBinMin, iTimeBinMax, mEnableNoiseSim, @@ -378,33 +414,14 @@ void HwClusterer::ProcessTimeBins(int iTimeBinMin, int iTimeBinMax) mNoiseObject, mPedestalObject }; -// std::cout << "starting CRU " << iCRU << std::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 - 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 - + 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]), - cfConfig); -// iCRU, -// mRowsMax, -// mPadsMax+2+2, -// iTimeBinMin, -// iTimeBinMax); - } } @@ -418,24 +435,43 @@ 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){ - 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(); + + // 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)); + if (mClusterMcLabelArray == nullptr) continue; + 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)) { + labelCount[l]++; + } } + for (auto &l : labelCount) labelSort.insert(l); + for (auto &l : labelSort) mClusterMcLabelArray->addElement(clusterPos,l.first); } } mLastTimebin = iTimeBinMax; } - + +void HwClusterer::setNumThreads(unsigned threads) { + mNumThreads = (threads == 0) ? std::thread::hardware_concurrency() : threads; +} + 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() 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+; 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(); 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..24f3c977c63a1 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(); @@ -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); @@ -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();