[AFIT-77] FDD reconstruction refactor

DataFormats/Detectors/FIT/FDD/include/DataFormatsFDD/FDDRecoConfig.h

@@ -13,13 +13,60 @@
 #define ALICEO2_FDD_RECO_PARAM
 
 #include "CommonUtils/ConfigurableParamHelper.h"
+#include "DataFormatsFDD/ChannelData.h"
 
 namespace o2::fdd
 {
 struct FDDRecoConfig : o2::conf::ConfigurableParamHelper<FDDRecoConfig> {
-  float AmplitudeCutOnCollisionTimeWeights = 3;
+  int16_t ChannelAmplitudeUpper = 4096;
+  int16_t ChannelAmplitudeLower = 0;
+  int16_t ChannelTimeLower = -2050;
+  int16_t ChannelTimeUpper = 2050;
+
+  uint8_t ChannelValidPmInputFlagMask = 0;
+  uint8_t ChannelValidPmInputFlags = 0;
+  bool isChannelDataGood(uint8_t flags) const
+  {
+    return (flags & ChannelValidPmInputFlagMask) == ChannelValidPmInputFlags;
+  }
+
+  uint8_t TriggerBitsValidationMask = 0;
+  uint8_t TriggerBitsValid = 0;
+  bool areTriggerBitsGood(uint8_t flags) const
+  {
+    return (flags & TriggerBitsValidationMask) == TriggerBitsValid;
+  }
+
+  int16_t CollisionTimeAmpltudeUpper = 4096;
+  int16_t CollisionTimeAmplitudeLower = 0;
+  int16_t CollisionTimeTimeLower = -500;
+  int16_t CollisionTimeTimeUpper = 500;
+
+  uint8_t CollisionTimeValidPmInputFlagMask = static_cast<uint8_t>(~(1u << ChannelData::kNumberADC));
+  uint8_t CollisionTimeValidPmInputFlags = static_cast<uint8_t>((1u << ChannelData::kIsCFDinADCgate) | (1u << ChannelData::kIsEventInTVDC));
+  bool isChannelDataGoodForCollisionTime(uint8_t flags) const
+  {
+    return (flags & CollisionTimeValidPmInputFlagMask) == CollisionTimeValidPmInputFlags;
+  }
+
+  int16_t SwOrGate = 153;
+  int16_t SwOrTriggerAmplitudeLowThresholdSideA = 6;
+  int16_t SwOrTriggerAmplitudeLowThresholdSideC = 0;
+  int16_t SwOrTriggerAmplitudeHighThresholdSideA = 4095;
+  int16_t SwOrTriggerAmplitudeHighThresholdSideC = 4095;
+  int16_t SwVtxTriggerTimeWindow = 100;
+  bool SwTriggerCoincidenceMode = true;
+  bool SwTriggerEnabled = false;
+
+  uint8_t SwTriggerValidPmInputFlagMask = static_cast<uint8_t>(~((1u << ChannelData::kNumberADC) | (1u << ChannelData::kIsEventInTVDC)));
+  uint8_t SwTriggerValidPmInputFlags = static_cast<uint8_t>(1u << ChannelData::kIsCFDinADCgate);
+  bool isChannelDataGoodForSwTrigger(uint8_t flags) const
+  {
+    return (flags & SwTriggerValidPmInputFlagMask) == SwTriggerValidPmInputFlags;
+  }
+
   O2ParamDef(FDDRecoConfig, "FDDRecoConfig");
-  ClassDefNV(FDDRecoConfig, 1);
+  ClassDefNV(FDDRecoConfig, 2);
 };
 } // namespace o2::fdd
-#endif
+#endif

DataFormats/Detectors/FIT/FDD/include/DataFormatsFDD/RecPoint.h

@@ -53,15 +53,19 @@ class RecPoint
  public:
   enum TimeTypeIndex : int { TimeA,
                              TimeC };
+  enum ETriggerBits { kSwOrTrgSideA = 5,
+                      kSwOrTrgSideC = 6,
+                      kSwVtxTrg = 7 };
   RecPoint() = default;
   RecPoint(const std::array<int, 2>& collisiontime, int first, int ne,
-           o2::InteractionRecord iRec, o2::fdd::Triggers triggers)
+           o2::InteractionRecord iRec, o2::fdd::Triggers triggers, bool swOrTrgSideA, bool swOrTrgSideC, bool swVtxTrg)
     : mCollisionTimePs(collisiontime)
   {
     mRef.setFirstEntry(first);
     mRef.setEntries(ne);
     mIntRecord = iRec;
     mTriggers = triggers;
+    updateTriggers(swOrTrgSideA, swOrTrgSideC, swVtxTrg);
   }
   ~RecPoint() = default;
 
@@ -86,12 +90,20 @@ class RecPoint
   bool operator==(const RecPoint&) const = default;
 
  private:
+  void updateTriggers(bool swOrTrgSideA, bool swOrTrgSideC, bool swVertexTrigger)
+  {
+    auto [triggerBits, technicalWords] = o2::fit::Triggers::parseDigitTriggerWord(mTriggers.getTriggersignals(), true);
+    triggerBits |= (uint8_t(swVertexTrigger) << kSwVtxTrg) | (uint8_t(swOrTrgSideA) << kSwOrTrgSideA) | (uint8_t(swOrTrgSideC) << kSwOrTrgSideC);
+    mTriggers.setTriggers(triggerBits);
+    mTechnicalWords = technicalWords;
+  }
   o2::dataformats::RangeReference<int, int> mRef;
   o2::InteractionRecord mIntRecord;
-  o2::fdd::Triggers mTriggers;                                                        // pattern of triggers  in this BC
+  o2::fdd::Triggers mTriggers; // pattern of triggers  in this BC
+  uint8_t mTechnicalWords{0};
   std::array<int, 2> mCollisionTimePs = {sDummyCollissionTime, sDummyCollissionTime}; // in picoseconds
 
-  ClassDefNV(RecPoint, 3);
+  ClassDefNV(RecPoint, 4);
 };
 } // namespace fdd
 } // namespace o2

Detectors/FIT/FDD/reconstruction/include/FDDReconstruction/Reconstructor.h

@@ -18,6 +18,8 @@
 #include "DataFormatsFDD/Digit.h"
 #include "DataFormatsFDD/RecPoint.h"
 #include "DataFormatsFIT/DeadChannelMap.h"
+#include "DataFormatsFDD/FDDRecoConfig.h"
+#include "FDDBase/Constants.h"
 namespace o2
 {
 namespace fdd
@@ -40,6 +42,52 @@ class Reconstructor
   }
 
  private:
+  static constexpr uint8_t coincidencePairs[Nchannels] = {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11};
+
+  bool isInAmplitudeWindow(int16_t amplitude)
+  {
+    return amplitude > FDDRecoConfig::Instance().ChannelAmplitudeLower && amplitude < FDDRecoConfig::Instance().ChannelAmplitudeUpper;
+  }
+  bool isInTimeWindow(int16_t time)
+  {
+    return time > FDDRecoConfig::Instance().ChannelTimeLower && time < FDDRecoConfig::Instance().ChannelTimeUpper;
+  }
+  bool isChannelDataGood(uint8_t channelBits)
+  {
+    return FDDRecoConfig::Instance().isChannelDataGood(channelBits);
+  }
+
+  bool isInCollisionTimeAmplitudeWindow(int16_t amplitude)
+  {
+    return amplitude > FDDRecoConfig::Instance().CollisionTimeAmplitudeLower && amplitude < FDDRecoConfig::Instance().CollisionTimeAmpltudeUpper;
+  }
+  bool isInCollisionTimeTimeWindow(int16_t time)
+  {
+    return time > FDDRecoConfig::Instance().CollisionTimeTimeLower && time < FDDRecoConfig::Instance().CollisionTimeTimeUpper;
+  }
+  bool isChannelDataGoodForCollisionTime(uint8_t channelBits)
+  {
+    return FDDRecoConfig::Instance().isChannelDataGoodForCollisionTime(channelBits);
+  }
+
+  bool checkSwOrGateConditions(const ChannelData& channelData)
+  {
+    bool isAmplitudeInRange = false;
+    if (channelData.mPMNumber >= sSideAStartId) {
+      isAmplitudeInRange = channelData.mChargeADC > FDDRecoConfig::Instance().SwOrTriggerAmplitudeLowThresholdSideA && channelData.mChargeADC < FDDRecoConfig::Instance().SwOrTriggerAmplitudeHighThresholdSideA;
+    } else {
+      isAmplitudeInRange = channelData.mChargeADC > FDDRecoConfig::Instance().SwOrTriggerAmplitudeLowThresholdSideC && channelData.mChargeADC < FDDRecoConfig::Instance().SwOrTriggerAmplitudeHighThresholdSideC;
+    }
+    bool isInTimeWindow = std::abs(channelData.mTime) < FDDRecoConfig::Instance().SwOrGate;
+    return isAmplitudeInRange && isInTimeWindow && FDDRecoConfig::Instance().isChannelDataGoodForSwTrigger(channelData.mFEEBits);
+  }
+
+  bool calculateSwTrigger(float meanTimeSideA, float meanTimeSideC, bool orGateSideA, bool orGateSideC)
+  {
+    return (std::abs(meanTimeSideC - meanTimeSideA) < o2::fdd::timePerTDC * FDDRecoConfig::Instance().SwVtxTriggerTimeWindow) && orGateSideA && orGateSideC;
+  }
+  static constexpr uint8_t sSideAStartId = 8;
+
   o2::fit::DeadChannelMap const* mDeadChannelMap = nullptr;
   ClassDefNV(Reconstructor, 3);
 };

Detectors/FIT/FDD/reconstruction/src/Reconstructor.cxx

@@ -27,66 +27,99 @@ void Reconstructor::process(o2::fdd::Digit const& digitBC, gsl::span<const o2::f
                             std::vector<o2::fdd::RecPoint>& RecPoints, std::vector<o2::fdd::ChannelDataFloat>& outChData)
 // gsl::span<o2::fdd::ChannelDataFloat> outChData)
 {
-  // Compute charge weighted average time
-  Double_t timeFDA = 0, timeFDC = 0;
-  Double_t weightFDA = 0.0, weightFDC = 0.0;
-  Int_t nInTimeA = 0, nInTimeC = 0;
-  int firstEntry = outChData.size();
-  int nStored = 0;
-  int nch = inChData.size();
-
-  for (int ich = 0; ich < nch; ich++) {
-    if (mDeadChannelMap && !mDeadChannelMap->isChannelAlive(inChData[ich].mPMNumber)) {
-      LOG(debug) << "Channel " << ich << " is dead - discarding data";
+  const int firstEntry = outChData.size();
+
+  bool atLeastOneVertexChannel = false;
+  bool atLeastOneChannelInSwAmplitudeRange = false;
+
+  int64_t timeSideA = 0;
+  int64_t timeSideC = 0;
+  int64_t entriesNumSideA = 0;
+  int64_t entriesNumSideC = 0;
+  uint32_t storedEntries = 0;
+
+  bool swOrGateSideA = false;
+  bool swOrGateSideC = false;
+  int64_t swTimeSideA = 0;
+  int64_t swTimeSideC = 0;
+  int64_t swEntriesNumSideA = 0;
+  int64_t swEntriesNumSideC = 0;
+
+  bool coincidencePairsOrTrg[Nchannels] = {false};
+  int64_t coincidencePairsTime[Nchannels] = {0};
+
+  for (const auto& channelData : inChData) {
+    if (mDeadChannelMap && !mDeadChannelMap->isChannelAlive(channelData.mPMNumber)) {
+      LOG(debug) << "Channel " << channelData.mPMNumber << " is dead - discarding data";
       continue;
     }
-    bool inTime = inChData[ich].getFlag(ChannelData::EEventDataBit::kIsEventInTVDC);
-    bool inAdcGate = inChData[ich].getFlag(ChannelData::EEventDataBit::kIsCFDinADCgate);
-    if (inAdcGate) {
-      outChData.emplace_back((int)inChData[ich].mPMNumber, (inChData[ich].mTime) * timePerTDC,
-                             (double)inChData[ich].mChargeADC, inChData[ich].mFEEBits);
-      nStored++;
+    if (!isInAmplitudeWindow(channelData.mChargeADC) || !isInTimeWindow(channelData.mTime) || !isChannelDataGood(channelData.mFEEBits)) {
+      continue;
     }
+    const auto& chNum = channelData.mPMNumber;
 
-    Float_t adc = inChData[ich].mChargeADC;
-    Float_t time = (inChData[ich].mTime) * timePerTDC;
-    if (time == o2::InteractionRecord::DummyTime) {
-      continue;
+    outChData.emplace_back((int)chNum, (channelData.mTime) * timePerTDC,
+                           (double)channelData.mChargeADC, channelData.mFEEBits);
+    storedEntries++;
+
+    if (FDDRecoConfig::Instance().SwTriggerEnabled) {
+      coincidencePairsOrTrg[chNum] |= checkSwOrGateConditions(channelData);
+      coincidencePairsTime[chNum] = channelData.mTime;
+      if (FDDRecoConfig::Instance().SwTriggerCoincidenceMode && coincidencePairsOrTrg[chNum] && coincidencePairsOrTrg[coincidencePairs[chNum]]) {
+        if (chNum >= sSideAStartId) {
+          swOrGateSideA = true;
+          swTimeSideA += channelData.mTime + coincidencePairsTime[coincidencePairs[chNum]];
+          swEntriesNumSideA += 2;
+        } else {
+          swOrGateSideC = true;
+          swTimeSideC += channelData.mTime + coincidencePairsTime[coincidencePairs[chNum]];
+          swEntriesNumSideC += 2;
+        }
+      } else if (!FDDRecoConfig::Instance().SwTriggerCoincidenceMode && coincidencePairsOrTrg[chNum]) {
+        if (chNum >= sSideAStartId) {
+          swOrGateSideA = true;
+          swTimeSideA += channelData.mTime;
+          swEntriesNumSideA++;
+        } else {
+          swOrGateSideC = true;
+          swTimeSideC += channelData.mTime;
+          swEntriesNumSideC++;
+        }
+      }
     }
-    Float_t timeErr = 1;
-    if (adc > FDDRecoConfig::Instance().AmplitudeCutOnCollisionTimeWeights) {
-      timeErr = 1. / adc;
+
+    if (!isChannelDataGoodForCollisionTime(channelData.mFEEBits) || !isInCollisionTimeAmplitudeWindow(channelData.mChargeADC) || !isInCollisionTimeTimeWindow(channelData.mTime)) {
+      continue;
     }
-    if ((int)inChData[ich].mPMNumber < 8) {
-      nInTimeC += inTime;
-      timeFDC += time / (timeErr * timeErr);
-      weightFDC += 1. / (timeErr * timeErr);
+
+    if (channelData.mPMNumber >= sSideAStartId) {
+      timeSideA += channelData.mTime;
+      entriesNumSideA++;
     } else {
-      nInTimeA += inTime;
-      timeFDA += time / (timeErr * timeErr);
-      weightFDA += 1. / (timeErr * timeErr);
+      timeSideC += channelData.mTime;
+      entriesNumSideC++;
     }
   }
-  const int nsToPs = 1e3;
+
+  constexpr float nsToPs = 1e3;
+  float meanTimeSideA = entriesNumSideA > 0 ? float(timeSideA) * timePerTDC / float(entriesNumSideA) : 0;
+  float meanTimeSideC = entriesNumSideC > 0 ? float(timeSideC) * timePerTDC / float(entriesNumSideC) : 0;
+
+  bool swVertexTrg = false;
+  if (FDDRecoConfig::Instance().SwTriggerEnabled) {
+    float swMeanTimeSideA = swEntriesNumSideA > 0 ? float(swTimeSideA) * timePerTDC / float(swEntriesNumSideA) : 0;
+    float swMeanTimeSideC = swEntriesNumSideC > 0 ? float(swTimeSideC) * timePerTDC / float(swEntriesNumSideC) : 0;
+    swVertexTrg = calculateSwTrigger(swMeanTimeSideA, swMeanTimeSideC, swOrGateSideA, swOrGateSideC);
+  }
+
   std::array<int, 2> mCollisionTime = {o2::fdd::RecPoint::sDummyCollissionTime, o2::fdd::RecPoint::sDummyCollissionTime};
-  /// Avg time for each side, only if one channel satisfy the TVDC condition (if not, also avg time is propagated for background study using AO2D)
-  if (nInTimeA > 0) {
-    mCollisionTime[o2::fdd::RecPoint::TimeA] = (weightFDA > 1) ? round(timeFDA / weightFDA * nsToPs)
-                                                               : o2::fdd::RecPoint::sDummyCollissionTime;
-  } else {
-    if (weightFDA > 0) {
-      mCollisionTime[o2::fdd::RecPoint::TimeA] = round(timeFDA / weightFDA * nsToPs);
-    }
+  if (entriesNumSideA > 0) {
+    mCollisionTime[o2::fdd::RecPoint::TimeA] = std::round(meanTimeSideA * nsToPs);
   }
-  if (nInTimeC > 0) {
-    mCollisionTime[o2::fdd::RecPoint::TimeC] = (weightFDC > 1) ? round(timeFDC / weightFDC * nsToPs)
-                                                               : o2::fdd::RecPoint::sDummyCollissionTime;
-  } else {
-    if (weightFDC > 0) {
-      mCollisionTime[o2::fdd::RecPoint::TimeC] = round(timeFDC / weightFDC * nsToPs);
-    }
+  if (entriesNumSideC > 0) {
+    mCollisionTime[o2::fdd::RecPoint::TimeC] = std::round(meanTimeSideC * nsToPs);
   }
-  RecPoints.emplace_back(mCollisionTime, firstEntry, nStored, digitBC.getIntRecord(), digitBC.mTriggers);
+  RecPoints.emplace_back(mCollisionTime, firstEntry, storedEntries, digitBC.getIntRecord(), digitBC.mTriggers, swOrGateSideA, swOrGateSideC, swVertexTrg);
 }
 //________________________________________________________
 void Reconstructor::finish()