FT0: update Digitizer signal shape and trigger logic; FV0: update trigger logic in digitizer

Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h

@@ -31,8 +31,8 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper<FT0DigParam> {
   float mAmpRecordUp = 15;           // to [ns]
   float hitTimeOffsetA = 0;          ///< hit time offset on the A side [ns]
   float hitTimeOffsetC = 0;          ///< hit time offset on the C side [ns]
-  int mtrg_central_trh = 600.;       // channels
-  int mtrg_semicentral_trh = 300.;   // channels
+  int mtrg_central_trh = 40;         // Tclu units (40 for pp and 1433 for PbPb in Run3)
+  int mtrg_semicentral_trh = 20;     // Tclu units (20 for pp and 35 for PbPb in Run3)
 
   float mMip_in_V = 7;       // MIP to mV
   float mPe_in_mip = 0.004;  // invserse Np.e. in MIP 1./250.
@@ -43,11 +43,12 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper<FT0DigParam> {
   float mNoiseVar = 0.1;              // noise level
   float mNoisePeriod = 1 / 0.9;       // GHz low frequency noise period;
   short mTime_trg_gate = 153;         // #channels as in TCM as in Pilot beams ('OR gate' setting in TCM tab in ControlServer)
+  short mTime_trg_vertex_gate = 100;  // #channels as in TCM as in Pilot beams ('OR gate' setting in TCM tab in ControlServer)
   float mAmpThresholdForReco = 5;     // only channels with amplitude higher will participate in calibration and collision time: 0.3 MIP
   short mTimeThresholdForReco = 1000; // only channels with time below will participate in calibration and collision time
 
-  float mMV_2_Nchannels = 2.2857143;          // amplitude channel 7 mV ->16channels
-  float mMV_2_NchannelsInverse = 0.437499997; // inverse amplitude channel 7 mV ->16channels
+  float mMV_2_Nchannels = 2.;         // amplitude channel 7 mV ->14channels
+  float mMV_2_NchannelsInverse = 0.5; // inverse amplitude channel 7 mV ->14channels (nowhere used)
 
   O2ParamDef(FT0DigParam, "FT0DigParam");
 };

Detectors/FIT/FT0/simulation/src/Digitizer.cxx

@@ -16,6 +16,13 @@
 #include "CommonConstants/PhysicsConstants.h"
 #include "CommonDataFormat/InteractionRecord.h"
 
+#include "DataFormatsFT0/LookUpTable.h"
+#include "FT0Base/Constants.h"
+#include <map>
+#include <array>
+#include <regex>
+#include <string>
+
 #include "TMath.h"
 #include "TRandom.h"
 #include <TH1F.h>
@@ -35,24 +42,84 @@ namespace o2::ft0
 template <typename Float>
 Float signalForm_i(Float x)
 {
-  using namespace std;
-  Float const a = -0.45458;
-  Float const b = -0.83344945;
-  return x > Float(0) ? -(exp(b * x) - exp(a * x)) / Float(7.8446501) : Float(0);
+  float p0, p1, p2, p3, p4, p5, p6, p7;
+  p0 = 1.30853;
+  p1 = 0.516807;
+  p2 = 3.36714;
+  p3 = -1.01206;
+  p4 = 1.42832;
+  p5 = 1.1589;
+  p6 = 1.22019;
+  p7 = 0.426818;
+
+  Double_t val = 0;
+
+  if (x > p3) {
+    Double_t x1 = x - p3;
+    Double_t arg1 = (log(x1) - p0) / p1;
+    val += p2 * (1.0 / (x1 * p1 * sqrt(2 * TMath::Pi()))) * exp(-0.5 * arg1 * arg1);
+  }
+
+  if (x > p7) {
+    Double_t x2 = x - p7;
+    Double_t arg2 = (log(x2) - p4) / p5;
+    val += p6 * (1.0 / (x2 * p5 * sqrt(2 * TMath::Pi()))) * exp(-0.5 * arg2 * arg2);
+  }
+
+  return val;
 };
 
 // integrated signal shape function
 inline float signalForm_integral(float x)
 {
-  using namespace std;
-  double const a = -0.45458;
-  double const b = -0.83344945;
-  if (x < 0) {
-    x = 0;
+  float p0, p1, p2, p3, p4, p5, p6, p7;
+  p0 = 1.30853;
+  p1 = 0.516807;
+  p2 = 3.36714;
+  p3 = -1.01206;
+  p4 = 1.42832;
+  p5 = 1.1589;
+  p6 = 1.22019;
+  p7 = 0.426818;
+  Double_t val = 0;
+
+  if (x > p3) {
+    Double_t x1 = x - p3;
+    Double_t z1 = (log(x1) - p0) / (sqrt(2) * p1);
+    val += p2 * 0.5 * (1 + TMath::Erf(z1)); // norm1 * CDF1
   }
-  return -(exp(b * x) / b - exp(a * x) / a) / 7.8446501;
+
+  if (x > p7) {
+    Double_t x2 = x - p7;
+    Double_t z2 = (log(x2) - p4) / (sqrt(2) * p5);
+    val += p6 * 0.5 * (1 + TMath::Erf(z2)); // norm2 * CDF2
+  }
+
+  return val;
+};
+/*
+// signal shape function
+template <typename Float>
+Float signalForm_i(Float x)
+{
+using namespace std;
+Float const a = -0.45458;
+Float const b = -0.83344945;
+return x > Float(0) ? -(exp(b * x) - exp(a * x)) / Float(7.8446501) : Float(0);
 };
 
+// integrated signal shape function
+inline float signalForm_integral(float x)
+{
+using namespace std;
+double const a = -0.45458;
+double const b = -0.83344945;
+if (x < 0) {
+  x = 0;
+}
+return -(exp(b * x) / b - exp(a * x) / a) / 7.8446501;
+};
+*/
 // SIMD version of the integrated signal shape function
 inline Vc::float_v signalForm_integralVc(Vc::float_v x)
 {
@@ -249,8 +316,64 @@ void Digitizer::storeBC(BCCache& bc,
   if (bc.hits.empty()) {
     return;
   }
+  // Initialize mapping channelID -> PM hash and PM side (A/C) using FT0 LUT
+  static bool pmLutInitialized = false;
+  static std::array<uint8_t, o2::ft0::Constants::sNCHANNELS_PM> mChID2PMhash{};
+  static std::map<uint8_t, bool> mMapPMhash2isAside; // hashed PM -> is A side
+
+  if (!pmLutInitialized) {
+    std::map<std::string, uint8_t> mapFEE2hash; // module name -> hashed PM id
+    uint8_t tcmHash = 0;
+
+    const auto& lut = o2::ft0::SingleLUT::Instance().getVecMetadataFEE();
+    auto lutSorted = lut;
+    std::sort(lutSorted.begin(), lutSorted.end(),
+              [](const auto& first, const auto& second) { return first.mModuleName < second.mModuleName; });
+
+    uint8_t binPos = 0;
+    for (const auto& lutEntry : lutSorted) {
+      const auto& moduleName = lutEntry.mModuleName;
+      const auto& moduleType = lutEntry.mModuleType;
+      const auto& strChID = lutEntry.mChannelID;
+
+      auto [it, inserted] = mapFEE2hash.insert({moduleName, binPos});
+      if (inserted) {
+        if (moduleName.find("PMA") != std::string::npos) {
+          mMapPMhash2isAside.insert({binPos, true});
+        } else if (moduleName.find("PMC") != std::string::npos) {
+          mMapPMhash2isAside.insert({binPos, false});
+        }
+        ++binPos;
+      }
+
+      if (std::regex_match(strChID, std::regex("^[0-9]{1,3}$"))) {
+        int chID = std::stoi(strChID);
+        if (chID < o2::ft0::Constants::sNCHANNELS_PM) {
+          mChID2PMhash[chID] = mapFEE2hash[moduleName];
+        } else {
+          LOG(fatal) << "Incorrect LUT entry: chID " << strChID << " | " << moduleName;
+        }
+      } else if (moduleType != "TCM") {
+        LOG(fatal) << "Non-TCM module w/o numerical chID: chID " << strChID << " | " << moduleName;
+      } else { // TCM
+        tcmHash = mapFEE2hash[moduleName];
+      }
+    }
+
+    pmLutInitialized = true;
+  }
+
   int n_hit_A = 0, n_hit_C = 0, mean_time_A = 0, mean_time_C = 0;
   int summ_ampl_A = 0, summ_ampl_C = 0;
+  int sum_A_ampl = 0, sum_C_ampl = 0;
+  int nPMTs = mGeometry.NCellsA * 4 + mGeometry.NCellsC * 4;
+  std::vector<int> sum_ampl_ipmt(nPMTs, 0);
+  // Per-PM summed charge (like in digits2trgFT0)
+  std::map<uint8_t, int> mapPMhash2sumAmpl;
+  for (const auto& entry : mMapPMhash2isAside) {
+    mapPMhash2sumAmpl.insert({entry.first, 0});
+  }
+
   int vertex_time;
   const auto& params = FT0DigParam::Instance();
   int first = digitsCh.size(), nStored = 0;
@@ -297,6 +420,10 @@ void Digitizer::storeBC(BCCache& bc,
     if (is_time_in_signal_gate) {
       chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsCFDinADCgate);
       chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsEventInTVDC);
+      // Sum channel charge per PM (similar logic as in digits2trgFT0)
+      if (ipmt < o2::ft0::Constants::sNCHANNELS_PM) {
+        mapPMhash2sumAmpl[mChID2PMhash[static_cast<uint8_t>(ipmt)]] += static_cast<int>(amp);
+      }
     }
     digitsCh.emplace_back(ipmt, smeared_time, int(amp), chain);
     nStored++;
@@ -308,6 +435,8 @@ void Digitizer::storeBC(BCCache& bc,
       continue;
     }
 
+    sum_ampl_ipmt[ipmt] += amp;
+
     if (is_A_side) {
       n_hit_A++;
       summ_ampl_A += amp;
@@ -318,17 +447,47 @@ void Digitizer::storeBC(BCCache& bc,
       mean_time_C += smeared_time;
     }
   }
+
+  for (size_t i = 0; i < sum_ampl_ipmt.size(); i++) {
+    sum_ampl_ipmt[i] = sum_ampl_ipmt[i] >> 3;
+    if (i < 4 * mGeometry.NCellsA) {
+      sum_A_ampl += sum_ampl_ipmt[i];
+    } else {
+      sum_C_ampl += sum_ampl_ipmt[i];
+    }
+  }
+
+  // Sum over PMs (using per-PM map) for debug/monitoring
+  int sum_PM_ampl_debug = 0;
+  int sum_PM_ampl_A_debug = 0;
+  int sum_PM_ampl_C_debug = 0;
+  for (const auto& entry : mapPMhash2sumAmpl) {
+    int pmAmpl = (entry.second >> 3);
+    sum_PM_ampl_debug += pmAmpl;
+    auto itSide = mMapPMhash2isAside.find(entry.first);
+    if (itSide != mMapPMhash2isAside.end()) {
+      if (itSide->second) {
+        sum_PM_ampl_A_debug += pmAmpl;
+      } else {
+        sum_PM_ampl_C_debug += pmAmpl;
+      }
+    }
+  }
+  LOG(debug) << "Sum PM amplitude (LUT-based): total=" << sum_PM_ampl_debug
+             << " A-side=" << sum_PM_ampl_A_debug
+             << " C-side=" << sum_PM_ampl_C_debug;
+
   Bool_t is_A, is_C, isVertex, is_Central, is_SemiCentral = 0;
   is_A = n_hit_A > 0;
   is_C = n_hit_C > 0;
-  is_Central = summ_ampl_A + summ_ampl_C >= params.mtrg_central_trh;
-  is_SemiCentral = summ_ampl_A + summ_ampl_C >= params.mtrg_semicentral_trh;
+  is_Central = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 * params.mtrg_central_trh;
+  is_SemiCentral = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 * params.mtrg_semicentral_trh && !is_Central;
   uint32_t amplA = is_A ? summ_ampl_A * 0.125 : -5000; // sum amplitude A side / 8 (hardware)
   uint32_t amplC = is_C ? summ_ampl_C * 0.125 : -5000; // sum amplitude C side / 8 (hardware)
   int timeA = is_A ? mean_time_A / n_hit_A : -5000;    // average time A side
   int timeC = is_C ? mean_time_C / n_hit_C : -5000;    // average time C side
   vertex_time = (timeC - timeA) * 0.5;
-  isVertex = is_A && is_C && (vertex_time > -params.mTime_trg_gate && vertex_time < params.mTime_trg_gate);
+  isVertex = is_A && is_C && (vertex_time > -params.mTime_trg_vertex_gate && vertex_time < params.mTime_trg_vertex_gate);
   LOG(debug) << " A " << is_A << " timeA " << timeA << " mean_time_A " << mean_time_A << "  n_hit_A " << n_hit_A << " C " << is_C << " timeC " << timeC << " mean_time_C " << mean_time_C << "  n_hit_C " << n_hit_C << " vertex_time " << vertex_time;
   Triggers triggers;
   bool isLaser = false;

Detectors/FIT/FV0/simulation/include/FV0Simulation/FV0DigParam.h

@@ -69,6 +69,10 @@ struct FV0DigParam : o2::conf::ConfigurableParamHelper<FV0DigParam> {
   uint8_t defaultChainQtc = 0x48;     // only 2 flags are set by default in simulation: kIsCFDinADCgate and kIsEventInTVDC
   float mAmpThresholdForReco = 24;    // only channels with amplitude higher will participate in calibration and collision time
   short mTimeThresholdForReco = 1000; // only channels with time below will participate in calibration and collision time
+  int NchannelsLevel = 2;             // trigger Nchannels
+  float InnerChargeLevel = 4;         // InnerRingsChargeLevel
+  float OuterChargeLevel = 4;         // OuterRingsChargeLevel
+  float ChargeLevel = 8;              // ChargeLevel
 
   O2ParamDef(FV0DigParam, "FV0DigParam");
 };

Detectors/FIT/FV0/simulation/src/Digitizer.cxx

@@ -38,8 +38,8 @@ void Digitizer::clear()
 void Digitizer::init()
 {
   LOG(info) << "init";
-  mNBins = FV0DigParam::Instance().waveformNbins;      //Will be computed using detector set-up from CDB
-  mBinSize = FV0DigParam::Instance().waveformBinWidth; //Will be set-up from CDB
+  mNBins = FV0DigParam::Instance().waveformNbins;                                  // Will be computed using detector set-up from CDB
+  mBinSize = FV0DigParam::Instance().waveformBinWidth;                             // Will be set-up from CDB
   mNTimeBinsPerBC = std::lround(o2::constants::lhc::LHCBunchSpacingNS / mBinSize); // 1920 bins/BC
 
   for (Int_t detID = 0; detID < Constants::nFv0Channels; detID++) {
@@ -149,8 +149,8 @@ void Digitizer::process(const std::vector<o2::fv0::Hit>& hits,
 
       createPulse(mipFraction, hit.GetTrackID(), hitTime, hit.GetPos().R(), cachedIR, nCachedIR, detId);
 
-    } //while loop
-  }   //hitloop
+    } // while loop
+  } // hitloop
 }
 
 void Digitizer::createPulse(float mipFraction, int parID, const double hitTime, const float hitR,
@@ -200,7 +200,7 @@ void Digitizer::createPulse(float mipFraction, int parID, const double hitTime,
     }
     added[ir] = true;
   }
-  ///Add MC labels to BCs for those contributed to the PMT signal
+  /// Add MC labels to BCs for those contributed to the PMT signal
   for (int ir = 0; ir < nCachedIR; ir++) {
     if (added[ir]) {
       auto bcCache = getBCCache(cachedIR[ir]);
@@ -238,6 +238,8 @@ void Digitizer::storeBC(const BCCache& bc,
   int8_t nTotFiredCells = 0;
   int8_t nTrgFiredCells = 0; // number of fired cells, that follow additional trigger conditions (time gate)
   int totalChargeAllRing = 0;
+  int totalChargeInnerRing = 0;
+  int totalChargeOuterRing = 0;
   int32_t avgTime = 0;
   double nSignalInner = 0;
   double nSignalOuter = 0;
@@ -285,8 +287,10 @@ void Digitizer::storeBC(const BCCache& bc,
       avgTime += iCfdZero;
       if (iPmt < 24) {
         nSignalInner++;
+        totalChargeInnerRing += iTotalCharge;
       } else {
         nSignalOuter++;
+        totalChargeOuterRing += iTotalCharge;
       }
     }
   }
@@ -300,24 +304,26 @@ void Digitizer::storeBC(const BCCache& bc,
   } else {
     avgTime = o2::fit::Triggers::DEFAULT_TIME;
   }
-  ///Triggers for FV0
-  bool isA, isAIn, isAOut, isCen, isSCen;
+  /// Triggers for FV0
+  bool isA, isNchannels, isAIn, isAOut, isTotalCharge;
   isA = nTrgFiredCells > 0;
-  isAIn = nSignalInner > 0;  // ring 1,2 and 3
-  isAOut = nSignalOuter > 0; // ring 4 and 5
-  isCen = totalChargeAllRing > FV0DigParam::Instance().adcChargeCenThr;
-  isSCen = totalChargeAllRing > FV0DigParam::Instance().adcChargeSCenThr;
+  isNchannels = nTrgFiredCells > FV0DigParam::Instance().NchannelsLevel;
+  // isAIn = nSignalInner > FV0DigParam::Instance().NchannelsLevel;  // ring 1,2 and 3
+  isAIn = 0.125 * totalChargeInnerRing > 2 * FV0DigParam::Instance().InnerChargeLevel; // ring 1,2 and 3
+  // isAOut = nSignalOuter > FV0DigParam::Instance().NchannelsLevel; // ring 4 and 5
+  isAOut = 0.125 * totalChargeOuterRing > 2 * FV0DigParam::Instance().OuterChargeLevel; // ring 4 and 5
+  isTotalCharge = 0.125 * totalChargeAllRing > 2 * FV0DigParam::Instance().ChargeLevel;
 
   Triggers triggers;
   const int unusedCharge = o2::fit::Triggers::DEFAULT_AMP;
   const int unusedTime = o2::fit::Triggers::DEFAULT_TIME;
   const int unusedZero = o2::fit::Triggers::DEFAULT_ZERO;
   const bool unusedBitsInSim = false; // bits related to laser and data validity
   const bool bitDataIsValid = true;
-  triggers.setTriggers(isA, isAIn, isAOut, isCen, isSCen, nTrgFiredCells, (int8_t)unusedZero,
+  triggers.setTriggers(isA, isAIn, isAOut, isTotalCharge, isNchannels, nTrgFiredCells, (int8_t)unusedZero,
                        (int32_t)(0.125 * totalChargeAllRing), (int32_t)unusedCharge, (int16_t)avgTime, (int16_t)unusedTime, unusedBitsInSim, unusedBitsInSim, bitDataIsValid);
   digitsBC.emplace_back(first, nTotFiredCells, bc, triggers, mEventId - 1);
-  digitsTrig.emplace_back(bc, isA, isAIn, isAOut, isCen, isSCen);
+  digitsTrig.emplace_back(bc, isA, isAIn, isAOut, isTotalCharge, isNchannels);
   for (auto const& lbl : bc.labels) {
     labels.addElement(nBC, lbl);
   }
@@ -342,8 +348,8 @@ Int_t Digitizer::SimulateLightYield(Int_t pmt, Int_t nPhot) const
 //---------------------------------------------------------------------------
 Float_t Digitizer::IntegrateCharge(const ChannelDigitF& pulse) const
 {
-  int const chargeIntMin = FV0DigParam::Instance().isIntegrateFull ? 0 : (FV0DigParam::Instance().avgCfdTimeForMip - 6.0) / mBinSize;                //Charge integration offset (cfd mean time - 6 ns)
-  int const chargeIntMax = FV0DigParam::Instance().isIntegrateFull ? mNTimeBinsPerBC : (FV0DigParam::Instance().avgCfdTimeForMip + 14.0) / mBinSize; //Charge integration offset (cfd mean time + 14 ns)
+  int const chargeIntMin = FV0DigParam::Instance().isIntegrateFull ? 0 : (FV0DigParam::Instance().avgCfdTimeForMip - 6.0) / mBinSize;                // Charge integration offset (cfd mean time - 6 ns)
+  int const chargeIntMax = FV0DigParam::Instance().isIntegrateFull ? mNTimeBinsPerBC : (FV0DigParam::Instance().avgCfdTimeForMip + 14.0) / mBinSize; // Charge integration offset (cfd mean time + 14 ns)
   if (chargeIntMin < 0 || chargeIntMin > mNTimeBinsPerBC || chargeIntMax > mNTimeBinsPerBC) {
     LOG(fatal) << "invalid indicess: chargeInMin=" << chargeIntMin << " chargeIntMax=" << chargeIntMax;
   }
@@ -400,7 +406,7 @@ float Digitizer::getDistFromCellCenter(UInt_t cellId, double hitx, double hity)
   double a = -(y0 - pCell->y) / (x0 - pCell->x);
   double b = 1;
   double c = -(y0 - a * x0);
-  //Return the distance from hit to this line
+  // Return the distance from hit to this line
   return (a * hitx + b * hity + c) / TMath::Sqrt(a * a + b * b);
 }