[EMCAL-565] Add cut on pre and post trigger pile-up (#13018)

- The EMCal QC observed a few runs with cells, which have a significant
  contribution of pre-trigger pile-up. These come from FECs that loose
the correct timeing signal during a run and then fire too early or too
late, resulting in a second peak
- A new cut is introduced that calculates the fraction of pre trigger
  (-500 - -25ns before the main time peak) and post trigger (25 - 500ns)
after the main timing peak
- The distribution of this post and pre-trigger pile-up fraction is then
  used to calculate the mean and sigma.
- The cut is chosen quite loose (defined in the CalibParams) in order to only capture the extreme
  cases

Author: jokonig

Common/Utils/include/CommonUtils/BoostHistogramUtils.h

@@ -714,6 +714,29 @@ auto ReduceBoostHistoFastSliceByValue(boost::histogram::histogram<axes...>& hist
   return ReduceBoostHistoFastSlice(hist2d, binXLow, binXHigh, binYLow, binYHigh, includeOverflowUnderflow);
 }
 
+/// \brief Function to integrate 1d boost histogram in specified range
+/// \param hist 1d boost histogram
+/// \param min lower integration range
+/// \param max upper integration range
+/// \return sum of bin contents in specified range
+template <typename... axes>
+double getIntegralBoostHist(boost::histogram::histogram<axes...>& hist, double min, double max)
+{
+  // find bins for min and max values
+  std::array<int, 2> axisLimitsIndex = {hist.axis(0).index(min), hist.axis(0).index(max)};
+  // over/underflow bin have to be protected
+  for (auto& bin : axisLimitsIndex) {
+    if (bin < 0) {
+      bin = 0;
+    } else if (bin >= hist.axis(0).size()) {
+      bin = hist.axis(0).size() - 1;
+    }
+  }
+  // Reduce histogram to desired range
+  auto slicedHist = ReduceBoostHistoFastSlice1D(hist, axisLimitsIndex[0], axisLimitsIndex[1], false);
+  return boost::histogram::algorithm::sum(slicedHist);
+}
+
 } // namespace utils
 } // end namespace o2
 

Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibExtractor.h

@@ -58,8 +58,12 @@ class EMCALCalibExtractor
   };
 
   struct BadChannelCalibTimeInfo {
-    std::array<double, 17664> sigmaCell; // sigma value of time distribution for single cells
-    double goodCellWindow;               // cut value for good cells
+    std::array<double, 17664> sigmaCell;         // sigma value of time distribution for single cells
+    double goodCellWindow;                       // cut value for good cells
+    std::array<double, 17664> fracHitsPreTrigg;  // fraction of hits before the main time peak (pre-trigger pile-up)
+    double goodCellWindowFracHitsPreTrigg;       // cut value for good cells for pre-trigger pile-up
+    std::array<double, 17664> fracHitsPostTrigg; // fraction of hits after the main time peak (post-trigger pile-up)
+    double goodCellWindowFracHitsPostTrigg;      // cut value for good cells for post-trigger pile-up
   };
 
  public:
@@ -191,6 +195,12 @@ class EMCALCalibExtractor
             if (calibrationTimeInfo.sigmaCell[cellID] > calibrationTimeInfo.goodCellWindow) {
               LOG(debug) << "Cell " << cellID << " is flagged due to time distribution";
               failed = true;
+            } else if (calibrationTimeInfo.fracHitsPreTrigg[cellID] > calibrationTimeInfo.goodCellWindowFracHitsPreTrigg) {
+              LOG(debug) << "Cell " << cellID << " is flagged due to time distribution (pre-trigger)";
+              failed = true;
+            } else if (calibrationTimeInfo.fracHitsPostTrigg[cellID] > calibrationTimeInfo.goodCellWindowFracHitsPostTrigg) {
+              LOG(debug) << "Cell " << cellID << " is flagged due to time distribution (post-trigger)";
+              failed = true;
             }
           }
         }
@@ -321,7 +331,7 @@ class EMCALCalibExtractor
   template <typename... axes>
   BadChannelCalibTimeInfo buildTimeMeanAndSigma(const boost::histogram::histogram<axes...>& histCellTime)
   {
-    std::array<double, 17664> meanSigma;
+    BadChannelCalibTimeInfo timeInfo;
     for (int i = 0; i < mNcells; ++i) {
       // calculate sigma per cell
       const int indexLow = histCellTime.axis(1).index(i);
@@ -333,24 +343,45 @@ class EMCALCalibExtractor
         maxElementIndex = 0;
       }
       float maxElementCenter = 0.5 * (boostHistCellSlice.axis(0).bin(maxElementIndex).upper() + boostHistCellSlice.axis(0).bin(maxElementIndex).lower());
-      meanSigma[i] = std::sqrt(o2::utils::getVarianceBoost1D(boostHistCellSlice, -999999, maxElementCenter - 50, maxElementCenter + 50));
+      timeInfo.sigmaCell[i] = std::sqrt(o2::utils::getVarianceBoost1D(boostHistCellSlice, -999999, maxElementCenter - 50, maxElementCenter + 50));
+
+      // get number of hits within mean+-25ns (trigger bunch), from -500ns to -25ns before trigger bunch (pre-trigger), and for 25ns to 500ns (post-trigger)
+      double sumTrigg = o2::utils::getIntegralBoostHist(boostHistCellSlice, maxElementCenter - 25, maxElementCenter + 25);
+      double sumPreTrigg = o2::utils::getIntegralBoostHist(boostHistCellSlice, maxElementCenter - 500, maxElementCenter - 25);
+      double sumPostTrigg = o2::utils::getIntegralBoostHist(boostHistCellSlice, maxElementCenter + 25, maxElementCenter + 500);
+
+      // calculate fraction of hits of post and pre-trigger to main trigger bunch
+      timeInfo.fracHitsPreTrigg[i] = sumTrigg == 0 ? 0. : sumPreTrigg / sumTrigg;
+      timeInfo.fracHitsPostTrigg[i] = sumTrigg == 0 ? 0. : sumPostTrigg / sumTrigg;
     }
 
     // get the mean sigma and the std. deviation of the sigma distribution
     // those will be the values we cut on
     double avMean = 0, avSigma = 0;
     TRobustEstimator robustEstimator;
-    robustEstimator.EvaluateUni(meanSigma.size(), meanSigma.data(), avMean, avSigma, 0.5 * meanSigma.size());
+    robustEstimator.EvaluateUni(timeInfo.sigmaCell.size(), timeInfo.sigmaCell.data(), avMean, avSigma, 0.5 * timeInfo.sigmaCell.size());
     // protection for the following case: For low statistics cases, it can happen that more than half of the cells is in one bin
     // in that case the sigma will be close to zero. In that case, we take 95% of the data to calculate the truncated mean
     if (std::abs(avMean) < 0.001 && std::abs(avSigma) < 0.001) {
-      robustEstimator.EvaluateUni(meanSigma.size(), meanSigma.data(), avMean, avSigma, 0.95 * meanSigma.size());
+      robustEstimator.EvaluateUni(timeInfo.sigmaCell.size(), timeInfo.sigmaCell.data(), avMean, avSigma, 0.95 * timeInfo.sigmaCell.size());
     }
-
-    BadChannelCalibTimeInfo timeInfo;
-    timeInfo.sigmaCell = meanSigma;
+    // timeInfo.sigmaCell = meanSigma;
     timeInfo.goodCellWindow = avMean + (avSigma * o2::emcal::EMCALCalibParams::Instance().sigmaTime_bc); // only upper limit needed
 
+    double avMeanPre = 0, avSigmaPre = 0;
+    robustEstimator.EvaluateUni(timeInfo.fracHitsPreTrigg.size(), timeInfo.fracHitsPreTrigg.data(), avMeanPre, avSigmaPre, 0.5 * timeInfo.fracHitsPreTrigg.size());
+    if (std::abs(avMeanPre) < 0.001 && std::abs(avSigmaPre) < 0.001) {
+      robustEstimator.EvaluateUni(timeInfo.fracHitsPreTrigg.size(), timeInfo.fracHitsPreTrigg.data(), avMeanPre, avSigmaPre, 0.95 * timeInfo.fracHitsPreTrigg.size());
+    }
+    timeInfo.goodCellWindowFracHitsPreTrigg = avMeanPre + (avSigmaPre * o2::emcal::EMCALCalibParams::Instance().sigmaTimePreTrigg_bc); // only upper limit needed
+
+    double avMeanPost = 0, avSigmaPost = 0;
+    robustEstimator.EvaluateUni(timeInfo.fracHitsPostTrigg.size(), timeInfo.fracHitsPostTrigg.data(), avMeanPost, avSigmaPost, 0.5 * timeInfo.fracHitsPostTrigg.size());
+    if (std::abs(avMeanPost) < 0.001 && std::abs(avSigmaPost) < 0.001) {
+      robustEstimator.EvaluateUni(timeInfo.fracHitsPostTrigg.size(), timeInfo.fracHitsPostTrigg.data(), avMeanPost, avSigmaPost, 0.95 * timeInfo.fracHitsPostTrigg.size());
+    }
+    timeInfo.goodCellWindowFracHitsPostTrigg = avMeanPost + (avSigmaPost * o2::emcal::EMCALCalibParams::Instance().sigmaTimePostTrigg_bc); // only upper limit needed
+
     return timeInfo;
   }
 

Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALCalibParams.h

@@ -43,6 +43,8 @@ struct EMCALCalibParams : public o2::conf::ConfigurableParamHelper<EMCALCalibPar
   float rangeTimeAxisHigh_bc = 500;     ///< maximum value of time for histogram range
   float minCellEnergyTime_bc = 0.1;     ///< minimum energy needed to fill the time histogram
   float sigmaTime_bc = 5;               ///< sigma value for the upper cut on the time-variance distribution
+  float sigmaTimePreTrigg_bc = 10.;     ///< sigma value for the upper cut on the fraction of cells in the pre-trigger region
+  float sigmaTimePostTrigg_bc = 10.;    ///< sigma value for the upper cut on the fraction of cells in the post-trigger region
   unsigned int slotLength_bc = 0;       ///< Lenght of the slot before calibration is triggered. If set to 0 calibration is triggered when hasEnoughData returns true
   bool UpdateAtEndOfRunOnly_bc = false; ///< switch to enable trigger of calibration only at end of run
   float minNHitsForMeanEnergyCut = 100; ///< mean number of hits per cell that is needed to cut on the mean energy per hit. Needed for high energy intervals as outliers can distort the distribution