TPC: Add possibility for MC anchoring of basic calibrations

Detectors/TPC/base/include/TPCBase/CDBInterface.h

@@ -135,6 +135,9 @@ class CDBInterface
   /// \return noise object
   const CalPad& getNoise();
 
+  /// Return the zero suppression threshold map
+  const CalPad& getZeroSuppressionThreshold();
+
   /// Return the gain map object
   ///
   /// The function checks if the object is already loaded and returns it
@@ -201,6 +204,9 @@ class CDBInterface
   /// \param default switch if to use default values
   void setUseDefaults(bool defaults = true) { mUseDefaults = defaults; }
 
+  /// return defaults usage
+  bool getUseDefaults() const { return mUseDefaults; }
+
   /// set CDB time stamp for object retrieval
   void setTimeStamp(long time)
   {
@@ -220,16 +226,18 @@ class CDBInterface
   {
     mPedestals.reset();
     mNoise.reset();
+    mZeroSuppression.reset();
     mGainMap.reset();
   }
 
  private:
   CDBInterface() = default;
 
   // ===| Pedestal and noise |==================================================
-  std::unique_ptr<CalPad> mPedestals; ///< Pedestal object
-  std::unique_ptr<CalPad> mNoise;     ///< Noise object
-  std::unique_ptr<CalPad> mGainMap;   ///< Gain map object
+  std::unique_ptr<CalPad> mPedestals;       ///< Pedestal object
+  std::unique_ptr<CalPad> mNoise;           ///< Noise object
+  std::unique_ptr<CalPad> mZeroSuppression; ///< Noise object
+  std::unique_ptr<CalPad> mGainMap;         ///< Gain map object
 
   // ===| switches and parameters |=============================================
   bool mUseDefaults = false; ///< use defaults instead of CCDB
@@ -244,6 +252,7 @@ class CDBInterface
   void loadGainMapFromFile();          ///< load gain map from mGainmapFileName
   void createDefaultPedestals();       ///< creation of default pedestals if requested
   void createDefaultNoise();           ///< creation of default noise if requested
+  void createDefaultZeroSuppression(); ///< creation of default noise if requested
   void createDefaultGainMap();         ///< creation of default gain map if requested
 
   template <typename T>

Detectors/TPC/base/include/TPCBase/ParameterElectronics.h

@@ -28,23 +28,24 @@ namespace tpc
 
 enum class DigitzationMode : char {
   FullMode = 0,         ///< Apply noise, pedestal and saturation
-  SubtractPedestal = 1, ///< Apply noise, pedestal and saturation and then from that subtract the pedestal
-  NoSaturation = 2,     ///< Apply only noise and pedestal
-  PropagateADC = 3      ///< Just propagate the bare ADC value
+  ZeroSuppression = 1,  ///< Apply noise, pedestal and saturation and then from that subtract the pedestal
+  SubtractPedestal = 2, ///< Apply noise, pedestal and saturation and then from that subtract the pedestal
+  NoSaturation = 3,     ///< Apply only noise and pedestal
+  PropagateADC = 4      ///< Just propagate the bare ADC value
 };
 
 struct ParameterElectronics : public o2::conf::ConfigurableParamHelper<ParameterElectronics> {
   static constexpr int TIMEBININBC = 8;
 
-  int NShapedPoints = 8;                                        ///< Number of ADC samples which are taken into account for a given, shaped signal (should fit
-                                                                /// into SSE registers)
-  float PeakingTime = 160e-3f;                                  ///< Peaking time of the SAMPA [us]
-  float ChipGain = 20.f;                                        ///< Gain of the SAMPA [mV/fC] - may be either 20 or 30
-  float ADCdynamicRange = 2200.f;                               ///< Dynamic range of the ADC [mV]
-  float ADCsaturation = 1024.f;                                 ///< ADC saturation [ADC counts]
+  int NShapedPoints = 8;                                                        ///< Number of ADC samples which are taken into account for a given, shaped signal (should fit
+                                                                                /// into SSE registers)
+  float PeakingTime = 160e-3f;                                                  ///< Peaking time of the SAMPA [us]
+  float ChipGain = 20.f;                                                        ///< Gain of the SAMPA [mV/fC] - may be either 20 or 30
+  float ADCdynamicRange = 2200.f;                                               ///< Dynamic range of the ADC [mV]
+  float ADCsaturation = 1024.f;                                                 ///< ADC saturation [ADC counts]
   float ZbinWidth = TIMEBININBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3; ///< Width of a z bin [us]
-  float ElectronCharge = 1.602e-19f;                            ///< Electron charge [C]
-  DigitzationMode DigiMode = DigitzationMode::SubtractPedestal; ///< Digitization mode [full / ... ]
+  float ElectronCharge = 1.602e-19f;                                            ///< Electron charge [C]
+  DigitzationMode DigiMode = DigitzationMode::ZeroSuppression;                  ///< Digitization mode [full / ... ]
 
   /// Average time from the start of the signal shaping to the COG of the sampled distribution
   ///

Detectors/TPC/base/src/CDBInterface.cxx

@@ -16,6 +16,7 @@
 // system includes
 #include <cxxabi.h>
 #include <ctime>
+#include <memory>
 #include <fmt/format.h>
 #include <fmt/chrono.h>
 
@@ -94,6 +95,21 @@ const CalPad& CDBInterface::getNoise()
   return *mNoise;
 }
 
+//______________________________________________________________________________
+const CalPad& CDBInterface::getZeroSuppressionThreshold()
+{
+  if (mUseDefaults) {
+    if (!mZeroSuppression) {
+      createDefaultZeroSuppression();
+    }
+    return *mZeroSuppression;
+  } else {
+    // return from CDB, assume that check for object existence are done there
+    return getObjectFromCDB<CalPadMapType>(CDBTypeMap.at(CDBType::ConfigFEEPad)).at("ThresholdMap");
+    ;
+  }
+}
+
 //______________________________________________________________________________
 const CalPad& CDBInterface::getGainMap()
 {
@@ -262,6 +278,19 @@ void CDBInterface::createDefaultNoise()
   }
 }
 
+//______________________________________________________________________________
+void CDBInterface::createDefaultZeroSuppression()
+{
+  // default map is 3*noise
+  mZeroSuppression = std::unique_ptr<CalPad>(new CalPad(getNoise()));
+  mZeroSuppression->setName("ThresholdMap");
+
+  for (auto& calArray : mZeroSuppression->getData()) {
+    auto& data = calArray.getData();
+    std::transform(data.begin(), data.end(), data.begin(), [](const auto value) { return 3.f * value; });
+  }
+}
+
 //______________________________________________________________________________
 void CDBInterface::createDefaultGainMap()
 {

Detectors/TPC/base/test/testTPCParameters.cxx

@@ -72,7 +72,7 @@ BOOST_AUTO_TEST_CASE(ParameterElectronics_test1)
   BOOST_CHECK_CLOSE(ParameterElectronics::Instance().ADCsaturation, 1024, 1e-3);
   BOOST_CHECK_CLOSE(ParameterElectronics::Instance().ZbinWidth, NominalTimeBin, 1e-3);
   BOOST_CHECK_CLOSE(ParameterElectronics::Instance().ElectronCharge, 1.602e-19, 1e-3);
-  BOOST_CHECK(ParameterElectronics::Instance().DigiMode == DigitzationMode::SubtractPedestal);
+  BOOST_CHECK(ParameterElectronics::Instance().DigiMode == DigitzationMode::ZeroSuppression);
 
   BOOST_CHECK(o2::conf::ConfigurableParam::getValueAs<int>("TPCEleParam.NShapedPoints") == 8);
   BOOST_CHECK_CLOSE(o2::conf::ConfigurableParam::getValueAs<float>("TPCEleParam.PeakingTime"), 160e-3, 1e-3);

Detectors/TPC/simulation/include/TPCSimulation/SAMPAProcessing.h

@@ -118,6 +118,9 @@ class SAMPAProcessing
   /// \return Noise on the channel of interest
   float getNoise(const int sector, const int globalPadInSector);
 
+  /// Get the zero suppression threshold for a given channel
+  float getZeroSuppression(const int sector, const int globalPadInSector) const;
+
   /// Get the pedestal for a given channel
   /// \param cru CRU of the channel of interest
   /// \param padPos PadPos of the channel of interest
@@ -127,11 +130,12 @@ class SAMPAProcessing
  private:
   SAMPAProcessing();
 
-  const ParameterGas* mGasParam;         ///< Caching of the parameter class to avoid multiple CDB calls
-  const ParameterDetector* mDetParam;    ///< Caching of the parameter class to avoid multiple CDB calls
-  const ParameterElectronics* mEleParam; ///< Caching of the parameter class to avoid multiple CDB calls
-  const CalPad* mNoiseMap;               ///< Caching of the parameter class to avoid multiple CDB calls
-  const CalPad* mPedestalMap;            ///< Caching of the parameter class to avoid multiple CDB calls
+  const ParameterGas* mGasParam;             ///< Caching of the parameter class to avoid multiple CDB calls
+  const ParameterDetector* mDetParam;        ///< Caching of the parameter class to avoid multiple CDB calls
+  const ParameterElectronics* mEleParam;     ///< Caching of the parameter class to avoid multiple CDB calls
+  const CalPad* mNoiseMap;                   ///< Caching of the parameter class to avoid multiple CDB calls
+  const CalPad* mPedestalMap;                ///< Caching of the parameter class to avoid multiple CDB calls
+  const CalPad* mZeroSuppression;            ///< Caching of the parameter class to avoid multiple CDB calls
   math_utils::RandomRing<> mRandomNoiseRing; ///< Ring with random number for noise
 };
 
@@ -158,11 +162,23 @@ inline float SAMPAProcessing::makeSignal(float ADCcounts, const int sector, cons
       return getADCSaturation(signal);
       break;
     }
+    case DigitzationMode::ZeroSuppression: {
+      signal -= commonMode;
+      signal += noise;
+      signal += pedestal;
+      const float signalSubtractPedestal = getADCSaturation(signal) - pedestal;
+      const float zeroSuppression = getZeroSuppression(sector, globalPadInSector);
+      if (signalSubtractPedestal < zeroSuppression) {
+        return 0.f;
+      }
+      return signalSubtractPedestal;
+      break;
+    }
     case DigitzationMode::SubtractPedestal: {
       signal -= commonMode;
       signal += noise;
       signal += pedestal;
-      float signalSubtractPedestal = getADCSaturation(signal) - pedestal;
+      const float signalSubtractPedestal = getADCSaturation(signal) - pedestal;
       return signalSubtractPedestal;
       break;
     }
@@ -239,6 +255,11 @@ inline float SAMPAProcessing::getNoise(const int sector, const int globalPadInSe
   return mRandomNoiseRing.getNextValue() * mNoiseMap->getValue(sector, globalPadInSector);
 }
 
+inline float SAMPAProcessing::getZeroSuppression(const int sector, const int globalPadInSector) const
+{
+  return mZeroSuppression->getValue(sector, globalPadInSector);
+}
+
 inline float SAMPAProcessing::getPedestal(const int sector, const int globalPadInSector) const
 {
   return mPedestalMap->getValue(sector, globalPadInSector);

Detectors/TPC/simulation/src/DigitContainer.cxx

@@ -42,6 +42,10 @@ void DigitContainer::fillOutputContainer(std::vector<Digit>& output,
           time.fillOutputContainer<DigitzationMode::FullMode>(output, mcTruth, commonModeOutput, sector, timeBin);
           break;
         }
+        case DigitzationMode::ZeroSuppression: {
+          time.fillOutputContainer<DigitzationMode::ZeroSuppression>(output, mcTruth, commonModeOutput, sector, timeBin);
+          break;
+        }
         case DigitzationMode::SubtractPedestal: {
           time.fillOutputContainer<DigitzationMode::SubtractPedestal>(output, mcTruth, commonModeOutput, sector, timeBin);
           break;

Detectors/TPC/simulation/src/GEMAmplification.cxx

@@ -19,7 +19,7 @@
 #include <TFile.h>
 #include "TPCBase/CDBInterface.h"
 #include <fstream>
-#include "FairLogger.h"
+#include "Framework/Logger.h"
 #include <filesystem>
 
 using namespace o2::tpc;
@@ -161,7 +161,7 @@ int GEMAmplification::getGEMMultiplication(int nElectrons, int GEM)
     return 0;
   } else if (nElectrons > 500) {
     /// For this condition the central limit theorem holds and we can approximate the amplification fluctuations by
-    ///a Gaussian for all electrons
+    /// a Gaussian for all electrons
     /// The mean is given by nElectrons * G_abs and the width by sqrt(nElectrons) * Sigma/Mu (Polya) * G_abs
     return ((mRandomGaus.getNextValue() * std::sqrt(static_cast<float>(nElectrons)) *
              mGasParam->SigmaOverMu) +

Detectors/TPC/simulation/src/SAMPAProcessing.cxx

@@ -20,7 +20,7 @@
 #include <iostream>
 #include <sstream>
 #include <string>
-#include "FairLogger.h"
+#include "Framework/Logger.h"
 
 using namespace o2::tpc;
 
@@ -37,6 +37,7 @@ void SAMPAProcessing::updateParameters()
   auto& cdb = CDBInterface::instance();
   mPedestalMap = &(cdb.getPedestals());
   mNoiseMap = &(cdb.getNoise());
+  mZeroSuppression = &(cdb.getZeroSuppressionThreshold());
 }
 
 void SAMPAProcessing::getShapedSignal(float ADCsignal, float driftTime, std::vector<float>& signalArray) const

Detectors/TPC/simulation/test/testTPCDigitContainer.cxx

@@ -19,6 +19,7 @@
 #include <boost/test/unit_test.hpp>
 #include <memory>
 #include <vector>
+#include <fmt/format.h>
 #include "DataFormatsTPC/Digit.h"
 #include "TPCSimulation/DigitContainer.h"
 #include "TPCSimulation/SAMPAProcessing.h"
@@ -36,7 +37,7 @@ BOOST_AUTO_TEST_CASE(DigitContainer_test1)
 {
   auto& cdb = CDBInterface::instance();
   cdb.setUseDefaults();
-  o2::conf::ConfigurableParam::updateFromString("TPCEleParam.DigiMode=3"); // propagate the ADC values, otherwise the computation get complicated
+  o2::conf::ConfigurableParam::updateFromString(fmt::format("TPCEleParam.DigiMode={}", (int)o2::tpc::DigitzationMode::PropagateADC)); // propagate the ADC values, otherwise the computation get complicated
   const Mapper& mapper = Mapper::instance();
   const SAMPAProcessing& sampa = SAMPAProcessing::instance();
   DigitContainer digitContainer;
@@ -96,7 +97,7 @@ BOOST_AUTO_TEST_CASE(DigitContainer_test2)
 {
   auto& cdb = CDBInterface::instance();
   cdb.setUseDefaults();
-  o2::conf::ConfigurableParam::updateFromString("TPCEleParam.DigiMode=3"); // propagate the ADC values, otherwise the computation get complicated
+  o2::conf::ConfigurableParam::updateFromString(fmt::format("TPCEleParam.DigiMode={}", (int)o2::tpc::DigitzationMode::PropagateADC)); // propagate the ADC values, otherwise the computation get complicated
   const Mapper& mapper = Mapper::instance();
   const SAMPAProcessing& sampa = SAMPAProcessing::instance();
   DigitContainer digitContainer;

Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx

@@ -119,6 +119,8 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
     mWithMCTruth = o2::conf::DigiParams::Instance().mctruth;
     auto useDistortions = ic.options().get<int>("distortionType");
     auto triggeredMode = ic.options().get<bool>("TPCtriggered");
+    mUseCalibrationsFromCCDB = ic.options().get<bool>("TPCuseCCDB");
+    LOG(info) << "TPC calibrations from CCDB: " << mUseCalibrationsFromCCDB;
 
     if (useDistortions > 0) {
       if (useDistortions == 1) {
@@ -233,7 +235,7 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
 
     /// For the time being use the defaults for the CDB
     auto& cdb = o2::tpc::CDBInterface::instance();
-    cdb.setUseDefaults();
+    cdb.setUseDefaults(!mUseCalibrationsFromCCDB);
     if (std::filesystem::exists("GainMap.root")) {
       LOG(info) << "TPC: Using gain map from 'GainMap.root'";
       cdb.setGainMapFromFile("GainMap.root");
@@ -250,7 +252,7 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
           delete mInternalROOTFlushFile;
           mInternalROOTFlushFile = nullptr;
         }
-        //TODO: make generic reset method?
+        // TODO: make generic reset method?
         mFlushCounter = 0;
         mDigitCounter = 0;
       }
@@ -466,6 +468,7 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer
   bool mWriteGRP = false;
   bool mWithMCTruth = true;
   bool mInternalWriter = false;
+  bool mUseCalibrationsFromCCDB = false;
 };
 
 o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter)
@@ -498,10 +501,13 @@ o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP,
     Inputs{InputSpec{"collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast<SubSpecificationType>(channel), Lifetime::Timeframe}},
     outputs,
     AlgorithmSpec{adaptFromTask<TPCDPLDigitizerTask>(internalwriter)},
-    Options{{"distortionType", VariantType::Int, 0, {"Distortion type to be used. 0 = no distortions (default), 1 = realistic distortions (not implemented yet), 2 = constant distortions"}},
-            {"initialSpaceChargeDensity", VariantType::String, "", {"Path to root file containing TH3 with initial space-charge density and name of the TH3 (comma separated)"}},
-            {"readSpaceCharge", VariantType::String, "", {"Path to root file containing pre-calculated space-charge object and name of the object (comma separated)"}},
-            {"TPCtriggered", VariantType::Bool, false, {"Impose triggered RO mode (default: continuous)"}}}};
+    Options{
+      {"distortionType", VariantType::Int, 0, {"Distortion type to be used. 0 = no distortions (default), 1 = realistic distortions (not implemented yet), 2 = constant distortions"}},
+      {"initialSpaceChargeDensity", VariantType::String, "", {"Path to root file containing TH3 with initial space-charge density and name of the TH3 (comma separated)"}},
+      {"readSpaceCharge", VariantType::String, "", {"Path to root file containing pre-calculated space-charge object and name of the object (comma separated)"}},
+      {"TPCtriggered", VariantType::Bool, false, {"Impose triggered RO mode (default: continuous)"}},
+      {"TPCuseCCDB", VariantType::Bool, false, {"true: load calibrations from CCDB; false: use random calibratoins"}},
+    }};
 }
 
 o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector<int> const& sectors, bool mctruth, bool internalwriter)