MCH: add digit noise filtering options

aphecetche · shahor02 · commit 8cad5f1eab5b · 2022-05-13T00:29:36.000+04:00
diff --git a/Detectors/MUON/MCH/DigitFiltering/CMakeLists.txt b/Detectors/MUON/MCH/DigitFiltering/CMakeLists.txt
@@ -10,7 +10,7 @@
 # or submit itself to any jurisdiction.
 
 o2_add_library(MCHDigitFiltering
-        SOURCES src/DigitFilterParam.cxx src/DigitFilteringSpec.cxx
+        SOURCES src/DigitFilter.cxx src/DigitFilterParam.cxx src/DigitFilteringSpec.cxx
         PUBLIC_LINK_LIBRARIES O2::MCHBase O2::Framework O2::SimulationDataFormat)
 
 o2_add_executable(
diff --git a/Detectors/MUON/MCH/DigitFiltering/README.md b/Detectors/MUON/MCH/DigitFiltering/README.md
@@ -5,4 +5,18 @@
 
 # MCH Digit Filtering
 
+```shell
+o2-mch-digits-filtering-workflow
+```
 
+Filter out some digits.
+
+Inputs :
+- list of all digits ([Digit](/DataFormats/Detectors/MUON/MCH/include/DataFormatsMCH/Digit.h)) in the current time frame, with the (default) data description `DIGITS` (can be changed with `--input-digits-data-description` option)
+- the list of ROF records ([ROFRecord](/DataFormats/Detectors/MUON/MCH/include/DataFormatsMCH/ROFRecord.h)) pointing to the digits associated to each interaction, with the (default) data description `DIGITROFS` (can be changed with `--input-digit-rofs-data-description` option)
+
+Outputs :
+- list of digits that pass the filtering criteria (for the moment ADC>0), with the (default) data description `F-DIGITS`  (can be changed with `--output-digits-data-description` option)
+- list of ROF records corresponding to the digits above, with a (default) data description of `F-DIGITROFS` (can be changed with `--output-digit-rofs-data-description` option)
+
+The exact behavior of the filtering is governed by the [MCHDigitFilterParam](/Detectors/MUON/MCH/DigitFiltering/include/MCHDigitFiltering/DigitFilterParam.h) configurable param, where you can select the minimum ADC value to consider, and whether to select signal (i.e. killing as much background as possible, possibly killing some signal as well) and/or to reject background (while not killing signal).
diff --git a/Detectors/MUON/MCH/DigitFiltering/include/MCHDigitFiltering/DigitFilter.h b/Detectors/MUON/MCH/DigitFiltering/include/MCHDigitFiltering/DigitFilter.h
@@ -0,0 +1,27 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef O2_MCH_DIGITFILTERING_DIGITFILTER_H_
+#define O2_MCH_DIGITFILTERING_DIGITFILTER_H_
+
+#include <functional>
+#include "DataFormatsMCH/Digit.h"
+
+namespace o2::mch
+{
+
+typedef std::function<bool(const Digit&)> DigitFilter;
+
+DigitFilter createDigitFilter(int minADC, bool rejectBackground, bool selectSignal);
+
+} // namespace o2::mch
+
+#endif
diff --git a/Detectors/MUON/MCH/DigitFiltering/include/MCHDigitFiltering/DigitFilterParam.h b/Detectors/MUON/MCH/DigitFiltering/include/MCHDigitFiltering/DigitFilterParam.h
@@ -24,9 +24,10 @@ namespace o2::mch
  */
 struct DigitFilterParam : public o2::conf::ConfigurableParamHelper<DigitFilterParam> {
 
-  bool sanityCheck = false; ///< whether or not to perform some sanity checks on the input digits
-  uint32_t minADC = 1;      ///< digits with an ADC below this value are discarded
-
+  bool sanityCheck = false;      ///< whether or not to perform some sanity checks on the input digits
+  uint32_t minADC = 1;           ///< digits with an ADC below this value are discarded
+  bool rejectBackground = false; ///< attempts to reject background (loose background selection, don't kill signal)
+  bool selectSignal = false;     ///< attempts to select only signal (strict background selection, might loose signal)
   O2ParamDef(DigitFilterParam, "MCHDigitFilter");
 };
 
diff --git a/Detectors/MUON/MCH/DigitFiltering/src/DigitFilter.cxx b/Detectors/MUON/MCH/DigitFiltering/src/DigitFilter.cxx
@@ -0,0 +1,120 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "MCHDigitFiltering/DigitFilter.h"
+
+#include "DataFormatsMCH/Digit.h"
+#include "MCHDigitFiltering/DigitFilterParam.h"
+#include <vector>
+#include <functional>
+#include <gsl/span>
+
+namespace
+{
+/** function used to select the signal.
+ * might cut some signal, the focus here is to kill as
+ * much background as possible.
+ */
+double signalCut(double* x, const double* p)
+{
+  double x0 = pow(p[0] / p[2], 1. / p[3]) + p[1];
+  if (x[0] < x0) {
+    return p[0];
+  } else {
+    return p[2] * pow(x[0] - p[1], p[3]);
+  }
+}
+
+/** function used to reject the background.
+ * might not kill all background, focus here is ensure
+ * we're not killing some signal along the way.
+ */
+double backgroundCut(double* x, const double* p)
+{
+
+  double x0 = (p[3] * p[2] - p[1] * p[0]) / (p[3] - p[1]);
+  if (x[0] < x0) {
+    return p[1] * (x[0] - p[0]);
+  } else {
+    return p[3] * (x[0] - p[2]);
+  }
+}
+
+o2::mch::DigitFilter createMinAdcCut(int minADC)
+{
+  return [minADC](const o2::mch::Digit& digit) -> bool {
+    if (digit.getADC() < minADC) {
+      return false;
+    }
+    return true;
+  };
+}
+
+o2::mch::DigitFilter createRejectBackground()
+{
+  int minNSamplesBackground = 14;
+  double backgroundParam[4] = {18., 24., -20., 7.0};
+
+  auto backgroundCut = [backgroundParam](double* x) {
+    return ::backgroundCut(x, backgroundParam);
+  };
+
+  return [backgroundCut, minNSamplesBackground](const o2::mch::Digit& digit) -> bool {
+    double nSample = digit.getNofSamples();
+    if (digit.getNofSamples() < minNSamplesBackground || digit.getADC() < backgroundCut(&nSample)) {
+      return false;
+    }
+    return true;
+  };
+}
+
+o2::mch::DigitFilter createSelectSignal()
+{
+  int minNSamplesSignal = 17;
+  double signalParam[4] = {80., 16., 12., 1.2};
+
+  auto signalCut = [signalParam](double* x) {
+    return ::signalCut(x, signalParam);
+  };
+  return [signalCut, minNSamplesSignal](const o2::mch::Digit& digit) -> bool {
+    double nSample = digit.getNofSamples();
+    if (digit.getNofSamples() < minNSamplesSignal || digit.getADC() < signalCut(&nSample)) {
+      return false;
+    }
+    return true;
+  };
+}
+} // namespace
+
+namespace o2::mch
+{
+DigitFilter createDigitFilter(int minADC, bool rejectBackground, bool selectSignal)
+{
+  std::vector<DigitFilter> parts;
+
+  parts.emplace_back(createMinAdcCut(minADC));
+  if (rejectBackground) {
+    parts.emplace_back(createRejectBackground());
+  }
+  if (selectSignal) {
+    parts.emplace_back(createSelectSignal());
+  }
+  return [parts](const Digit& digit) {
+    for (const auto& p : parts) {
+      if (!p(digit)) {
+        return false;
+      }
+    }
+    return true;
+  };
+}
+
+} // namespace o2::mch
diff --git a/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx b/Detectors/MUON/MCH/DigitFiltering/src/DigitFilteringSpec.cxx
@@ -20,10 +20,12 @@
 #include "Framework/Task.h"
 #include "Framework/WorkflowSpec.h"
 #include "MCHBase/SanityCheck.h"
+#include "MCHDigitFiltering/DigitFilter.h"
 #include "MCHDigitFiltering/DigitFilterParam.h"
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
 #include <fmt/format.h>
+#include <functional>
 #include <iostream>
 #include <string>
 #include <vector>
@@ -41,12 +43,15 @@ class DigitFilteringTask
   void init(InitContext& ic)
   {
     mSanityCheck = DigitFilterParam::Instance().sanityCheck;
-    mMinADC = DigitFilterParam::Instance().minADC;
-  }
-
-  bool isGoodDigit(const Digit& digit) const
-  {
-    return digit.getADC() >= mMinADC;
+    int minADC = DigitFilterParam::Instance().minADC;
+    bool rejectBackground = DigitFilterParam::Instance().rejectBackground;
+    mIsGoodDigit = createDigitFilter(minADC, rejectBackground, false);
+    // at digit filtering stage it is important to keep the 3rd parameter
+    // to false in the call above : the idea is to not cut too much
+    // on the tails of the charge distributions otherwise the clustering
+    // resolution will suffer.
+    // That's why we only apply the "reject background" filter, which
+    // is a loose background cut that does not penalize the signal
   }
 
   void run(ProcessingContext& pc)
@@ -84,7 +89,7 @@ class DigitFilteringTask
         // filter the digits from the current ROF
         for (auto i = 0; i < digits.size(); i++) {
           const auto& d = digits[i];
-          if (isGoodDigit(d)) {
+          if (mIsGoodDigit(d)) {
             oDigits.emplace_back(d);
             if (iLabels) {
               oLabels->addElements(oLabels->getIndexedSize(), iLabels->getLabels(i + irof.getFirstIdx()));
@@ -119,7 +124,7 @@ class DigitFilteringTask
  private:
   bool mSanityCheck;
   bool mUseMC;
-  int mMinADC;
+  DigitFilter mIsGoodDigit;
 };
 
 framework::DataProcessorSpec
diff --git a/Detectors/MUON/MCH/Workflow/README.md b/Detectors/MUON/MCH/Workflow/README.md
@@ -81,20 +81,7 @@ filePath = /home/data/data-de819-ped-raw.raw
 
 ## Digit filtering
 
-```shell
-o2-mch-digits-filtering-workflow
-```
-
-Filter out some digits. For the moment only removes digits that have a null ADC.
-
-Inputs :
-- list of all digits ([Digit](/DataFormats/Detectors/MUON/MCH/include/DataFormatsMCH/Digit.h)) in the current time frame, with the (default) data description `DIGITS` (can be changed with `--input-digits-data-description` option)
-- the list of ROF records ([ROFRecord](../../../../DataFormats/Detectors/MUON/MCH/include/DataFormatsMCH/ROFRecord.h)) pointing to the digits associated to each interaction, with the (default) data description `DIGITROFS` (can be changed with `--input-digit-rofs-data-description` option)
-
-Outputs :
-- list of digits that pass the filtering criteria (for the moment ADC>0), with the (default) data description `F-DIGITS`  (can be changed with `--output-digits-data-description` option)
-- list of ROF records corresponding to the digits above, with a (default) data description of `F-DIGITROFS` (can be changed with `--output-digit-rofs-data-description` option)
-
+Filter out (i.e. remove) some digits [more...](/Detectors/MUON/MCH/DigitFiltering/README.md)
 
 ## Time clustering
 
