add TPC cluster reject options (#1376)

* change output container in task only in one place

* add cluster reject options

* fix HwClusterer test case 6

* add option to not clear output container first

* fix clang-format

* update ClustererSpec.cxx

Author: KlewinS

Detectors/TPC/reconstruction/include/TPCReconstruction/ClustererTask.h

@@ -34,44 +34,46 @@ namespace TPC{
 class ClustererTask : public FairTask{
 
   using MCLabelContainer = o2::dataformats::MCTruthContainer<o2::MCCompLabel>;
+  //  using OutputType          = Cluster;
+  using OutputType = ClusterHardwareContainer8kb;
 
-  public:
-   /// Default constructor
-   /// \param sectorid Sector to be processed
-   ClustererTask(int sectorid = -1);
+ public:
+  /// Default constructor
+  /// \param sectorid Sector to be processed
+  ClustererTask(int sectorid = -1);
 
-   /// Destructor
-   ~ClustererTask() override = default;
+  /// Destructor
+  ~ClustererTask() override = default;
 
-   /// Initializes the clusterer and connects input and output container
-   InitStatus Init() override;
+  /// Initializes the clusterer and connects input and output container
+  InitStatus Init() override;
 
-   /// Clusterization
-   void Exec(Option_t* option) override;
+  /// Clusterization
+  void Exec(Option_t* option) override;
 
-   /// Complete Clusterization
-   void FinishTask() override;
+  /// Complete Clusterization
+  void FinishTask() override;
 
-   /// Switch for triggered / continuous readout
-   /// \param isContinuous - false for triggered readout, true for continuous readout
-   void setContinuousReadout(bool isContinuous);
+  /// Switch for triggered / continuous readout
+  /// \param isContinuous - false for triggered readout, true for continuous readout
+  void setContinuousReadout(bool isContinuous);
 
-  private:
-   bool mIsContinuousReadout; ///< Switch for continuous readout
-   int mEventCount;           ///< Event counter
-   int mClusterSector;        ///< Sector to be processed
+ private:
+  bool mIsContinuousReadout; ///< Switch for continuous readout
+  int mEventCount;           ///< Event counter
+  int mClusterSector;        ///< Sector to be processed
 
-   std::unique_ptr<HwClusterer> mHwClusterer; ///< Hw Clusterfinder instance
+  std::unique_ptr<HwClusterer> mHwClusterer; ///< Hw Clusterfinder instance
 
-   // Digit arrays
-   std::unique_ptr<const std::vector<Digit>> mDigitsArray;     ///< Array of TPC digits
-   std::unique_ptr<const MCLabelContainer> mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray
+  // Digit arrays
+  std::unique_ptr<const std::vector<Digit>> mDigitsArray;     ///< Array of TPC digits
+  std::unique_ptr<const MCLabelContainer> mDigitMCTruthArray; ///< Array for MCTruth information associated to digits in mDigitsArrray
 
-   // Cluster arrays
-   std::unique_ptr<std::vector<ClusterHardwareContainer8kb>> mHwClustersArray; ///< Array of clusters found by Hw Clusterfinder
-   std::unique_ptr<MCLabelContainer> mHwClustersMCTruthArray;                  ///< Array for MCTruth information associated to cluster in mHwClustersArrays
+  // Cluster arrays
+  std::unique_ptr<std::vector<OutputType>> mHwClustersArray; ///< Array of clusters found by Hw Clusterfinder
+  std::unique_ptr<MCLabelContainer> mHwClustersMCTruthArray; ///< Array for MCTruth information associated to cluster in mHwClustersArrays
 
-   ClassDefOverride(ClustererTask, 1)
+  ClassDefOverride(ClustererTask, 1)
 };
 
 inline

Detectors/TPC/reconstruction/include/TPCReconstruction/HwClusterer.h

@@ -76,12 +76,16 @@ class HwClusterer : public Clusterer
   /// Process digits
   /// \param digits Container with TPC digits
   /// \param mcDigitTruth MC Digit Truth container
+  /// \param clearContainerFirst Clears the outpcontainer for clusters and MC labels first, before processing
   void process(std::vector<o2::TPC::Digit> const& digits, MCLabelContainer const* mcDigitTruth) override;
+  void process(std::vector<o2::TPC::Digit> const& digits, MCLabelContainer const* mcDigitTruth, bool clearContainerFirst);
 
   /// Finish processing digits
   /// \param digits Container with TPC digits
   /// \param mcDigitTruth MC Digit Truth container
+  /// \param clearContainerFirst Clears the outpcontainer for clusters and MC labels first, before processing
   void finishProcess(std::vector<o2::TPC::Digit> const& digits, MCLabelContainer const* mcDigitTruth) override;
+  void finishProcess(std::vector<o2::TPC::Digit> const& digits, MCLabelContainer const* mcDigitTruth, bool clearContainerFirst);
 
   /// Switch for triggered / continuous readout
   /// \param isContinuous - false for triggered readout, true for continuous readout
@@ -95,6 +99,24 @@ class HwClusterer : public Clusterer
   /// \param charge Threshold which will be used
   void setContributionChargeThreshold(unsigned charge);
 
+  /// Switch to reject single pad clusters
+  /// \param doReject - true to reject clusters with sigmaPad2Pre == 0
+  void setRejectSinglePadClusters(bool doReject);
+
+  /// Switch to reject single time clusters
+  /// \param doReject - true to reject clusters with sigmaTime2Pre == 0
+  void setRejectSingleTimeClusters(bool doReject);
+
+  /// Switch to reject peaks in following timebins
+  /// \param doReject - true to reject peak, false to NOT reject peaks in following bins
+  void setRejectLaterTimebin(bool doReject);
+
+  /// Switch for mode, how the charge should be shared among nearby clusters
+  /// \param mode   0 for no splitting, charge is used for all found peaks,
+  ///               1 for minimum contributes half to all peaks
+  ///               2 for minimum contributes only to left/older peak
+  void setSplittingMode(short mode);
+
  private:
   /*
    * Helper functions
@@ -103,23 +125,23 @@ class HwClusterer : public Clusterer
   /// HW Cluster Processor
   /// \param peakMask       VC-mask with only peaks enabled
   /// \param qMaxIndex      Buffer index of center pad
-  /// \param center_pad     Pad number to be checked for cluster
-  /// \param center_time    Time to be checked for cluster
+  /// \param centerPad     Pad number to be checked for cluster
+  /// \param centerTime    Time to be checked for cluster
   /// \param row            Row number for cluster properties
-  void hwClusterProcessor(Vc::uint_m peakMask, unsigned qMaxIndex, short center_pad, int center_time, unsigned short row);
+  void hwClusterProcessor(Vc::uint_m peakMask, unsigned qMaxIndex, short centerPad, int centerTime, unsigned short row);
 
   /// HW Peak Finder
   /// \param qMaxIndex        Buffer index of central pad
-  /// \param center_pad       Pad number to be checked for cluster
-  /// \param center_time      Time to be checked for cluster
+  /// \param centerPad        Pad number to be checked for cluster
+  /// \param mappedCenterTime Time to be checked for cluster, mapped in available time space
   /// \param row              Row number for cluster properties
-  void hwPeakFinder(unsigned qMaxIndex, short center_pad, int center_time, unsigned short row);
+  void hwPeakFinder(unsigned qMaxIndex, short centerPad, int mappedCenterTime, unsigned short row);
 
   /// Helper function to update cluster properties and MC labels
   /// \param selectionMask  VC-mask with slected pads enabled
   /// \param row            Current row
-  /// \param center_pad     Pad of peak
-  /// \param center_time    Timebin of peak
+  /// \param centerPad     Pad of peak
+  /// \param centerTime    Timebin of peak
   /// \param dp             delta pad
   /// \param dt             delta time
   /// \param qTot           Total charge
@@ -128,7 +150,7 @@ class HwClusterer : public Clusterer
   /// \param sigmaPad2      Weighted sigma pad ^2 parameter
   /// \param sigmaTime2     Weighted sigma time ^2 parameter
   /// \param mcLabel        Vector with MClabel-counter-pairs
-  void updateCluster(const Vc::uint_m selectionMask, int row, short center_pad, int center_time, short dp, short dt, Vc::uint_v& qTot, Vc::int_v& pad, Vc::int_v& time, Vc::int_v& sigmaPad2, Vc::int_v& sigmaTime2, std::vector<std::unique_ptr<std::vector<std::pair<MCCompLabel, unsigned>>>>& mcLabels);
+  void updateCluster(const Vc::uint_m selectionMask, int row, short centerPad, int centerTime, short dp, short dt, Vc::uint_v& qTot, Vc::int_v& pad, Vc::int_v& time, Vc::int_v& sigmaPad2, Vc::int_v& sigmaTime2, std::vector<std::unique_ptr<std::vector<std::pair<MCCompLabel, unsigned>>>>& mcLabels, Vc::uint_m splitMask = Vc::Mask<uint>(false));
 
   /// Writes clusters from temporary storage to cluster output
   /// \param timeOffset   Time offset of cluster container
@@ -146,7 +168,7 @@ class HwClusterer : public Clusterer
   /// \param clear    Clears data buffer afterwards (for not continuous readout)
   void finishFrame(bool clear = false);
 
-  /// Clears the buffer at given timebin TODO: and fills timebin with noise + pedestal
+  /// Clears the buffer at given timebin
   /// \param timebin  Timebin to be cleared
   void clearBuffer(int timebin);
 
@@ -163,29 +185,25 @@ class HwClusterer : public Clusterer
   /// \param value  some value
   Vc::uint_v getFpOfADC(const Vc::uint_v value);
 
-  /// Does the comparison between two pads and sets the bits accordingly
-  /// \param qMaxIndex      Buffer index of central pad
-  /// \param compareIndex   Buffer index of pad to compare with
-  /// \param bitMax         Bit to be set for peak finding
-  /// \param bitMin         Bit to be set for minimum finding
-  /// \param row            Row number
-  void compareForPeak(const unsigned qMaxIndex, const unsigned compareIndex, const unsigned bitMax, const unsigned bitMin, const unsigned short row);
-
   /*
    * class members
    */
-  static const int mTimebinsInBuffer = 5;
+  static const int mTimebinsInBuffer = 6;
 
   unsigned short mNumRows;               ///< Number of rows in this sector
   unsigned short mNumRowSets;            ///< Number of row sets (Number of rows / Vc::Size) in this sector
   short mCurrentMcContainerInBuffer;     ///< Bit field, where to find the current MC container in buffer
+  short mSplittingMode;                  ///< Cluster splitting mode, 0 no splitting, 1 for minimum contributes half to both, 2 for miminum corresponds to left/older cluster
   int mClusterSector;                    ///< Sector to be processed
   int mLastTimebin;                      ///< Last time bin of previous event
   unsigned mLastHB;                      ///< Last HB bin of previous event
   unsigned mPeakChargeThreshold;         ///< Charge threshold for the central peak in ADC counts
   unsigned mContributionChargeThreshold; ///< Charge threshold for the contributing pads in ADC counts
   unsigned mClusterCounter;              ///< Cluster counter in output container for MC truth matching
   bool mIsContinuousReadout;             ///< Switch for continuous readout
+  bool mRejectSinglePadClusters;         ///< Switch to reject single pad clusters, sigmaPad2Pre == 0
+  bool mRejectSingleTimeClusters;        ///< Switch to reject single time clusters, sigmaTime2Pre == 0
+  bool mRejectLaterTimebin;              ///< Switch to reject peaks in later timebins of the same pad
 
   std::vector<unsigned short> mPadsPerRow;                       ///< Number of pads for given row (offset of 2 pads on both sides is already added)
   std::vector<unsigned short> mPadsPerRowSet;                    ///< Number of pads for given row set (offset of 2 pads on both sides is already added), a row set combines rows for parallel SIMD processing
@@ -205,6 +223,16 @@ class HwClusterer : public Clusterer
   MCLabelContainer* mClusterMcLabelArray;                  ///< Pointer to MC Label container
 };
 
+inline void HwClusterer::process(std::vector<o2::TPC::Digit> const& digits, o2::dataformats::MCTruthContainer<o2::MCCompLabel> const* mcDigitTruth)
+{
+  process(digits, mcDigitTruth, true);
+}
+
+inline void HwClusterer::finishProcess(std::vector<o2::TPC::Digit> const& digits, o2::dataformats::MCTruthContainer<o2::MCCompLabel> const* mcDigitTruth)
+{
+  finishProcess(digits, mcDigitTruth, true);
+}
+
 inline void HwClusterer::setContinuousReadout(bool isContinuous)
 {
   mIsContinuousReadout = isContinuous;
@@ -220,6 +248,26 @@ inline void HwClusterer::setContributionChargeThreshold(unsigned charge)
   mContributionChargeThreshold = charge;
 }
 
+inline void HwClusterer::setRejectSinglePadClusters(bool doReject)
+{
+  mRejectSinglePadClusters = doReject;
+}
+
+inline void HwClusterer::setRejectSingleTimeClusters(bool doReject)
+{
+  mRejectSingleTimeClusters = doReject;
+}
+
+inline void HwClusterer::setRejectLaterTimebin(bool doReject)
+{
+  mRejectLaterTimebin = doReject;
+}
+
+inline void HwClusterer::setSplittingMode(short mode)
+{
+  mSplittingMode = mode;
+}
+
 inline int HwClusterer::mapTimeInRange(int time)
 {
   return (mTimebinsInBuffer + (time % mTimebinsInBuffer)) % mTimebinsInBuffer;
@@ -239,57 +287,6 @@ inline short HwClusterer::getFirstSetBitOfField()
   return -1;
 }
 
-inline void HwClusterer::compareForPeak(const unsigned qMaxIndex, const unsigned compareIndex, const unsigned bitMax, const unsigned bitMin, const unsigned short row)
-{
-
-  const auto tmpMask = getFpOfADC(mDataBuffer[row][qMaxIndex]) >= getFpOfADC(mDataBuffer[row][compareIndex]);
-
-  // current center could be peak in one direction
-  where(tmpMask) | mDataBuffer[row][qMaxIndex] |= (0x1 << bitMax);
-
-  // other is smaller than center
-  // and if other one was not a peak candidate before (bit is 0), it is a minimum in one direction,
-  // so bitMin has to be set to inverse of bitMax of pad to compare
-  where(tmpMask) | mDataBuffer[row][compareIndex] |= (~mDataBuffer[row][compareIndex] & (0x1 << bitMax)) >> (bitMax - bitMin);
-
-  // other is not peak
-  where(tmpMask) | mDataBuffer[row][compareIndex] &= ~(0x1 << bitMax);
-
-  // other is peak if bit was already set
-  where(!tmpMask) | mDataBuffer[row][compareIndex] |= (mDataBuffer[row][compareIndex] & (0x1 << bitMax));
-}
-
-inline void HwClusterer::updateCluster(
-  const Vc::uint_m selectionMask, int row, short center_pad, int center_time, short dp, short dt,
-  Vc::uint_v& qTot, Vc::int_v& pad, Vc::int_v& time, Vc::int_v& sigmaPad2, Vc::int_v& sigmaTime2,
-  std::vector<std::unique_ptr<std::vector<std::pair<MCCompLabel, unsigned>>>>& mcLabels)
-{
-  const int mappedTime = mapTimeInRange(center_time + dt);
-  const int index = mappedTime * mPadsPerRowSet[row] + center_pad + dp;
-
-  where(selectionMask) | qTot += getFpOfADC(mDataBuffer[row][index]);
-  where(selectionMask) | pad += getFpOfADC(mDataBuffer[row][index]) * dp;
-  where(selectionMask) | time += getFpOfADC(mDataBuffer[row][index]) * dt;
-  where(selectionMask) | sigmaPad2 += getFpOfADC(mDataBuffer[row][index]) * dp * dp;
-  where(selectionMask) | sigmaTime2 += getFpOfADC(mDataBuffer[row][index]) * dt * dt;
-
-  for (int i = 0; i < Vc::uint_v::Size; ++i) {
-    if (selectionMask[i] && mMCtruth[mappedTime] != nullptr) {
-      for (auto& label : mMCtruth[mappedTime]->getLabels(mIndexBuffer[row][index][i])) {
-        bool isKnown = false;
-        for (auto& vecLabel : *mcLabels[i]) {
-          if (label == vecLabel.first) {
-            ++vecLabel.second;
-            isKnown = true;
-          }
-        }
-        if (!isKnown) {
-          mcLabels[i]->emplace_back(label, 1);
-        }
-      }
-    }
-  }
-}
 }
 }
 

Detectors/TPC/reconstruction/src/ClustererTask.cxx

@@ -72,7 +72,7 @@ InitStatus ClustererTask::Init()
   }
 
   // Register output container
-  mHwClustersArray = std::make_unique<std::vector<ClusterHardwareContainer8kb>>();
+  mHwClustersArray = std::make_unique<std::vector<OutputType>>();
   // a trick to register the unique pointer with FairRootManager
   static auto clusterArrayTmpPtr = mHwClustersArray.get();
   mgr->RegisterAny(Form("TPCClusterHW%i", mClusterSector), clusterArrayTmpPtr, kTRUE);