[EMCAL-686] Implement error handling in Raw fitter

Error handling is now exception-based, where
the errors thrown are defined as error codes in
enum class, and not inherit from std::exception.
Raw fitter errors are caught in the
RawToCellConverterSpec and pushed to the
error container which is pushed to the QC.

Futhermore a faulty 0 check was found in the
CaloRawFitterGamma2 in the parabola fit
where the numeric check for 0 is replaced by
DBL_EPSILON.

Author: mfasDa

Detectors/EMCAL/reconstruction/include/EMCALReconstruction/AltroDecoder.h

@@ -64,6 +64,10 @@ class AltroDecoderError : public std::exception
   /// \return the error type
   static ErrorType_t intToErrorType(int errornumber);
 
+  /// \brief Get the number of error types handled by the AltroDecoderError
+  /// \return Number of error types
+  static constexpr int getNumberOfErrorTypes() noexcept { return 8; }
+
   /// \brief Access to the error type connected to the erro
   /// \return Error type
   const ErrorType_t getErrorType() const noexcept { return mErrorType; }

Detectors/EMCAL/reconstruction/include/EMCALReconstruction/CaloRawFitter.h

@@ -39,13 +39,38 @@ class CaloRawFitter
 {
 
  public:
+  /**
+   * \enum RawFitterError_t
+   * \brief Error codes for failures in raw fitter procedure
+   */
+  enum class RawFitterError_t {
+    SAMPLE_UNINITIALIZED, ///< Samples not initialized or length is 0
+    FIT_ERROR,            ///< Fit procedure failed
+    CHI2_ERROR,           ///< Chi2 cannot be determined (usually due to insufficient amount of samples)
+    BUNCH_NOT_OK          ///< Bunch selection failed
+  };
+
+  /// \brief Create error message for a given error type
+  /// \param fiterror Fit error type
+  /// \return Error message connected to the error type
+  static std::string createErrorMessage(RawFitterError_t fiterror);
+
+  /// \brief Convert error type to numeric representation
+  /// \param fiterror Fit error type
+  /// \return Numeric representation of the raw fitter error
+  static int getErrorNumber(RawFitterError_t fiterror);
+
+  /// \brief Get the number of raw fit error types supported
+  /// \return Number of error types (4)
+  static constexpr int getNumberOfErrorTypes() noexcept { return 4; }
+
   /// \brief Constructor
   CaloRawFitter(const char* name, const char* nameshort);
 
   /// \brief Destructor
   virtual ~CaloRawFitter() = default;
 
-  virtual CaloFitResults evaluate(const std::vector<Bunch>& bunchvector,
+  virtual CaloFitResults evaluate(const gsl::span<const Bunch> bunchvector,
                                   std::optional<unsigned int> altrocfg1,
                                   std::optional<unsigned int> altrocfg2) = 0;
 
@@ -58,7 +83,7 @@ class CaloRawFitter
   /// \return pedestal,
   /// \return first time bin,
   /// \return last time bin,
-  std::tuple<int, int, float, short, short, float, int, int> preFitEvaluateSamples(const std::vector<Bunch>& bunchvector,
+  std::tuple<int, int, float, short, short, float, int, int> preFitEvaluateSamples(const gsl::span<const Bunch> bunchvector,
                                                                                    std::optional<unsigned int> altrocfg1, std::optional<unsigned int> altrocfg2, int acut);
 
   /// \brief The require time range if the maximum ADC value is between min and max (timebin)
@@ -103,7 +128,7 @@ class CaloRawFitter
   /// \return The index of the array with the maximum aplitude
   /// \return The bin where we have a maximum amp
   /// \return The maximum ADC signal
-  std::tuple<short, short, short> selectBunch(const std::vector<Bunch>& bunchvector);
+  std::tuple<short, short, short> selectBunch(const gsl::span<const Bunch>& bunchvector);
 
   /// \brief Selection of subset of data from one bunch that will be used for fitting or Peak finding.
   /// Go to the left and right of index of the maximum time bin
@@ -113,7 +138,7 @@ class CaloRawFitter
   std::tuple<int, int> selectSubarray(const gsl::span<double> data, short maxindex, int cut) const;
 
   /// \brief Pedestal evaluation if not zero suppressed
-  float evaluatePedestal(const std::vector<uint16_t>& data, std::optional<int> length) const;
+  float evaluatePedestal(const gsl::span<const uint16_t> data, std::optional<int> length) const;
 
   /// \brief Calculates the chi2 of the fit
   ///

Detectors/EMCAL/reconstruction/include/EMCALReconstruction/CaloRawFitterGamma2.h

@@ -52,8 +52,10 @@ class CaloRawFitterGamma2 final : public CaloRawFitter
   int getNiterationsMax() { return mNiterationsMax; }
 
   /// \brief Evaluation Amplitude and TOF
-  /// return Container with the fit results (amp, time, chi2, ...)
-  CaloFitResults evaluate(const std::vector<Bunch>& bunchvector,
+  /// \param
+  /// \throw RawFitterError_t::FIT_ERROR in case the peak fit failed
+  /// \return Container with the fit results (amp, time, chi2, ...)
+  CaloFitResults evaluate(const gsl::span<const Bunch> bunchvector,
                           std::optional<unsigned int> altrocfg1,
                           std::optional<unsigned int> altrocfg2) final;
 
@@ -62,12 +64,22 @@ class CaloRawFitterGamma2 final : public CaloRawFitter
   int mNiterationsMax = 15; ///< max number of iteraions
 
   /// \brief Fits the raw signal time distribution
-  /// \return chi2, fit status.
-  std::tuple<float, bool> doFit_1peak(int firstTimeBin, int nSamples, float& ampl, float& time);
+  /// \param firstTimeBin First timebin in the ALTRO bunch
+  /// \param nSamples Number of time samples of the ALTRO bunch
+  /// \param[in] ampl Initial guess of the amplitude for the fit
+  /// \param[out] ampl Amplitude result of the peak fit
+  /// \param[in] time Initial guess of the time for the fit
+  /// \param[out] time Time result of the peak fit
+  /// \return chi2 of the fit
+  /// \throw RawFitterError_t::FIT_ERROR in case of fit errors (insufficient number of time samples, matrix diagonalization error, ...)
+  float doFit_1peak(int firstTimeBin, int nSamples, float& ampl, float& time);
 
   /// \brief Fits the raw signal time distribution
+  /// \param maxTimeBin Time bin of the max. amplitude
   /// \return the fit parameters: amplitude, time.
-  std::tuple<float, float> doParabolaFit(int x) const;
+  ///
+  /// Fit performed as parabola fit to the signal
+  std::tuple<float, float> doParabolaFit(int maxTimeBin) const;
 
   ClassDefNV(CaloRawFitterGamma2, 1);
 }; // End of CaloRawFitterGamma2

Detectors/EMCAL/reconstruction/include/EMCALReconstruction/CaloRawFitterStandard.h

@@ -50,20 +50,27 @@ class CaloRawFitterStandard final : public CaloRawFitter
   ~CaloRawFitterStandard() final = default;
 
   /// \brief Approximate response function of the EMCal electronics.
-  /// \param x: bin
-  /// \param par: function parameters
+  /// \param x bin
+  /// \param par function parameters
   /// \return double with signal for a given time bin
   static double rawResponseFunction(double* x, double* par);
 
   /// \brief Evaluation Amplitude and TOF
-  /// return Container with the fit results (amp, time, chi2, ...)
-  CaloFitResults evaluate(const std::vector<Bunch>& bunchvector,
+  /// \param bunchvector Calo bunches for the tower and event
+  /// \param altrocfg1 ALTRO config register 1 from RCU trailer
+  /// \param altrocfg2 ALTRO config register 2 from RCU trailer
+  /// \return Container with the fit results (amp, time, chi2, ...)
+  /// \throw RawFitterError_t in case the fit failed (including all possible errors from upstream)
+  CaloFitResults evaluate(const gsl::span<const Bunch> bunchvector,
                           std::optional<unsigned int> altrocfg1,
                           std::optional<unsigned int> altrocfg2) final;
 
-  /// \brief Fits the raw signal time distribution
-  /// \return the fit parameters: amplitude, time, chi2, fit status.
-  std::tuple<float, float, float, bool> fitRaw(int firstTimeBin, int lastTimeBin) const;
+  /// \brief Fits the raw signal time distribution using TMinuit
+  /// \param firstTimeBin First timebin of the ALTRO bunch
+  /// \param lastTimeBin Last timebin of the ALTRO bunch
+  /// \return the fit parameters: amplitude, time, chi2
+  /// \throw RawFitter_t::FIT_ERROR in case the fit failed (insufficient number of samples or fit error from MINUIT)
+  std::tuple<float, float, float> fitRaw(int firstTimeBin, int lastTimeBin) const;
 
  private:
   ClassDefNV(CaloRawFitterStandard, 1);

Detectors/EMCAL/reconstruction/macros/RawFitterTESTs.C

@@ -97,10 +97,14 @@ void RawFitterTESTs(const char* filename = "")
           std::cout << "processing next channel idx " << chan.getChannelIndex() << ", " << chan.getHardwareAddress() << std::endl;
           // define the conatiner for the fit results, and perform the raw fitting using the stadnard raw fitter
           continue;
-          o2::emcal::CaloFitResults fitResults = RawFitter.evaluate(chan.getBunches(), 0, 0);
+          try {
+            o2::emcal::CaloFitResults fitResults = RawFitter.evaluate(chan.getBunches(), 0, 0);
 
-          // print the fit output
-          std::cout << "The Time is : " << fitResults.getTime() << " And the Amplitude is : " << fitResults.getAmp() << std::endl;
+            // print the fit output
+            std::cout << "The Time is : " << fitResults.getTime() << " And the Amplitude is : " << fitResults.getAmp() << std::endl;
+          } catch (o2::emcal::CaloRawFitter::RawFitterError_t& fiterror) {
+            std::cerr << "Error processing raw fit: " << o2::emcal::CaloRawFitter::createErrorMessage(fiterror) << std::endl;
+          }
         }
       }
     }

Detectors/EMCAL/reconstruction/src/CaloRawFitter.cxx

@@ -25,6 +25,40 @@
 
 using namespace o2::emcal;
 
+std::string CaloRawFitter::createErrorMessage(CaloRawFitter::RawFitterError_t errorcode)
+{
+  switch (errorcode) {
+    case RawFitterError_t::SAMPLE_UNINITIALIZED:
+      return "Sample for fit not initialzied or bunch length is 0";
+    case RawFitterError_t::FIT_ERROR:
+      return "Fit of the raw bunch was not successful";
+    case RawFitterError_t::CHI2_ERROR:
+      return "Chi2 of the fit could not be determined";
+    case RawFitterError_t::BUNCH_NOT_OK:
+      return "Calo bunch could not be selected";
+  };
+  // Silence compiler warnings for false positives
+  // can never enter here due to usage of enum class
+  return "Unknown error code";
+}
+
+int CaloRawFitter::getErrorNumber(CaloRawFitter::RawFitterError_t fiterror)
+{
+  switch (fiterror) {
+    case RawFitterError_t::SAMPLE_UNINITIALIZED:
+      return 0;
+    case RawFitterError_t::FIT_ERROR:
+      return 1;
+    case RawFitterError_t::CHI2_ERROR:
+      return 2;
+    case RawFitterError_t::BUNCH_NOT_OK:
+      return 3;
+  };
+  // Silence compiler warnings for false positives
+  // can never enter here due to usage of enum class
+  return -1;
+}
+
 //Default constructor
 CaloRawFitter::CaloRawFitter(const char* name, const char* nameshort) : mMinTimeIndex(-1),
                                                                         mMaxTimeIndex(-1),
@@ -108,7 +142,7 @@ std::tuple<float, std::array<double, constants::EMCAL_MAXTIMEBINS>> CaloRawFitte
 {
   std::array<double, constants::EMCAL_MAXTIMEBINS> outarray;
   int length = bunch.getBunchLength();
-  const std::vector<uint16_t>& sig = bunch.getADC();
+  const gsl::span<const uint16_t> sig(bunch.getADC());
 
   double ped = evaluatePedestal(sig, length);
 
@@ -119,8 +153,11 @@ std::tuple<float, std::array<double, constants::EMCAL_MAXTIMEBINS>> CaloRawFitte
   return std::make_tuple(ped, outarray);
 }
 
-float CaloRawFitter::evaluatePedestal(const std::vector<uint16_t>& data, std::optional<int> length) const
+float CaloRawFitter::evaluatePedestal(const gsl::span<const uint16_t> data, std::optional<int> length) const
 {
+  if (!mNsamplePed) {
+    throw RawFitterError_t::SAMPLE_UNINITIALIZED;
+  }
   double tmp = 0;
 
   if (mIsZerosupressed == false) {
@@ -173,15 +210,15 @@ bool CaloRawFitter::checkBunchEdgesForMax(const Bunch& bunch) const
   return bunchOK;
 }
 
-std::tuple<short, short, short> CaloRawFitter::selectBunch(const std::vector<Bunch>& bunchvector)
+std::tuple<short, short, short> CaloRawFitter::selectBunch(const gsl::span<const Bunch>& bunchvector)
 {
   short bunchindex = -1;
   short maxall = -1;
   int indx = -1;
   short maxampbin(0), maxamplitude(0);
 
   for (unsigned int i = 0; i < bunchvector.size(); i++) {
-    short max = maxAmp(bunchvector.at(i), indx); // CRAP PTH, bug fix, trouble if more than one bunches
+    short max = maxAmp(bunchvector[i], indx); // CRAP PTH, bug fix, trouble if more than one bunches
     if (isInTimeRange(indx, mMaxTimeIndex, mMinTimeIndex)) {
       if (max > maxall) {
         maxall = max;
@@ -193,9 +230,9 @@ std::tuple<short, short, short> CaloRawFitter::selectBunch(const std::vector<Bun
   }
 
   if (bunchindex >= 0) {
-    bool bunchOK = checkBunchEdgesForMax(bunchvector.at(bunchindex));
+    bool bunchOK = checkBunchEdgesForMax(bunchvector[bunchindex]);
     if (!bunchOK) {
-      bunchindex = -1;
+      throw RawFitterError_t::BUNCH_NOT_OK;
     }
   }
 
@@ -217,7 +254,7 @@ double CaloRawFitter::calculateChi2(double amp, double time,
 {
   if (first == last || first < 0) { // signal consists of single sample, chi2 estimate (0) not too well defined..
                                     // or, first is negative, the indices are not valid
-    return -1;
+    throw RawFitterError_t::CHI2_ERROR;
   }
 
   int nsamples = last - first + 1;
@@ -243,7 +280,7 @@ double CaloRawFitter::calculateChi2(double amp, double time,
 
   return chi2;
 }
-std::tuple<int, int, float, short, short, float, int, int> CaloRawFitter::preFitEvaluateSamples(const std::vector<Bunch>& bunchvector,
+std::tuple<int, int, float, short, short, float, int, int> CaloRawFitter::preFitEvaluateSamples(const gsl::span<const Bunch> bunchvector,
                                                                                                 std::optional<unsigned int> altrocfg1, std::optional<unsigned int> altrocfg2, int acut)
 {
 
@@ -261,8 +298,8 @@ std::tuple<int, int, float, short, short, float, int, int> CaloRawFitter::preFit
     //std::tie(ped, mReversed) = reverseAndSubtractPed((bunchvector.at(index)), altrocfg1, altrocfg2);
     //maxf = (float)*std::max_element(mReversed.begin(), mReversed.end());
 
-    int length = bunchvector.at(index).getBunchLength();
-    const std::vector<uint16_t>& sig = bunchvector.at(index).getADC();
+    int length = bunchvector[index].getBunchLength();
+    const std::vector<uint16_t>& sig = bunchvector[index].getADC();
 
     ped = evaluatePedestal(sig, length);
 
@@ -276,8 +313,8 @@ std::tuple<int, int, float, short, short, float, int, int> CaloRawFitter::preFit
     if (maxf >= acut) // possibly significant signal
     {
       // select array around max to possibly be used in fit
-      maxrev = maxampindex - bunchvector.at(index).getStartTime();
-      std::tie(first, last) = selectSubarray(gsl::span<double>(&mReversed[0], bunchvector.at(index).getBunchLength()), maxrev, acut);
+      maxrev = maxampindex - bunchvector[index].getStartTime();
+      std::tie(first, last) = selectSubarray(gsl::span<double>(&mReversed[0], bunchvector[index].getBunchLength()), maxrev, acut);
 
       // sanity check: maximum should not be in first or last bin
       // if we should do a fit