FT0: continue crosstalk handling in digitizer

Author: Szymon Pulawski

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

@@ -50,6 +50,8 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper<FT0DigParam> {
   float mMV_2_Nchannels = 2.;         // amplitude channel 7 mV ->14channels
   float mMV_2_NchannelsInverse = 0.5; // inverse amplitude channel 7 mV ->14channels (nowhere used)
 
+  float Cross_Talk_Frac = 0.05f;      // Crosstalk between channels
+
   O2ParamDef(FT0DigParam, "FT0DigParam");
 };
 } // namespace o2::ft0

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

@@ -20,6 +20,7 @@
 #include "FT0Base/Constants.h"
 #include <map>
 #include <array>
+#include <unordered_map>
 #include <regex>
 #include <string>
 
@@ -376,19 +377,53 @@ void Digitizer::storeBC(BCCache& bc,
 
   int vertex_time;
   const auto& params = FT0DigParam::Instance();
+
+  static bool pmGroupsInitialized = false;
+  static std::vector<std::array<int, 4>> pmtChannelGroups;
+  if (!pmGroupsInitialized) {
+    std::unordered_map<uint8_t, std::vector<int>> tmpGroups;
+    for (int ch = 0; ch < o2::ft0::Constants::sNCHANNELS_PM; ++ch) {
+      tmpGroups[mChID2PMhash[static_cast<uint8_t>(ch)]].push_back(ch);
+    }
+
+    for (auto& [pmHash, chVec] : tmpGroups) {
+      std::sort(chVec.begin(), chVec.end());
+      if (chVec.size() % 4 != 0) {
+        LOG(fatal) << "PM hash " << int(pmHash) << " has " << chVec.size()
+                   << " channels in LUT, expected multiplicity of 4";
+      }
+      for (size_t i = 0; i < chVec.size(); i += 4) {
+        std::array<int, 4> arr = {chVec[i + 0], chVec[i + 1], chVec[i + 2], chVec[i + 3]};
+        pmtChannelGroups.push_back(arr);
+      }
+    }
+    pmGroupsInitialized = true;
+  }
+
   int first = digitsCh.size(), nStored = 0;
   auto& particles = bc.hits;
   std::sort(std::begin(particles), std::end(particles));
   auto channel_end = particles.begin();
   std::vector<float> channel_times;
-//  std::set<int> disabledChannels = {40, 41, 42, 43, 88, 89, 90, 91, 56, 57, 58, 59, 60, 61, 62, 63, 72, 73, 74, 75, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 164, 165, 166, 167, 184, 185, 186, 187, 160, 161, 162, 163, 188, 189, 190, 191, 156, 157, 158, 159, 192, 193, 194, 195, 152, 153, 154, 155, 196, 197, 198, 199, 148, 149, 150, 151, 144, 145, 146, 147, 204, 205, 206, 207, 200, 201, 202, 203};  // przykładowe kanały 
+  std::vector<float> baseAmp(params.mMCPs, 0.f);
+  std::vector<float> finalAmp(params.mMCPs, 0.f);
+  std::vector<int> chTime(params.mMCPs, -5000);
+  std::vector<int> chChain(params.mMCPs, 0);
+  std::vector<bool> chValid(params.mMCPs, false);
+
+  static const std::array<std::array<int, 3>, 4> localNeighbours = {{
+    {{1, 2, 3}},
+    {{0, 3, 2}},
+    {{0, 3, 1}},
+    {{1, 2, 0}}
+  }};
+
+//  std::set<int> disabledChannels = {40, 41, 42, 43, 88, 89, 90, 91, 56, 57, 58, 59, 60, 61, 62, 63, 72, 73, 74, 75, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 164, 165, 166, 167, 184, 185, 186, 187, 160, 161, 162, 163, 188, 189, 190, 191, 156, 157, 158, 159, 192, 193, 194, 195, 152, 153, 154, 155, 196, 197, 198, 199, 148, 149, 150, 151, 144, 145, 146, 147, 204, 205, 206, 207, 200, 201, 202, 203};  // przykładowe kanały
   for (Int_t ipmt = 0; ipmt < params.mMCPs; ++ipmt) {
     auto channel_begin = channel_end;
     channel_end = std::find_if(channel_begin, particles.end(),
                                [ipmt](BCCache::particle const& p) { return p.hit_ch != ipmt; });
 
-    // The hits between 'channel_begin' and 'channel_end' now contains all hits for channel 'ipmt'
-
     if (channel_end - channel_begin < params.mAmp_trsh) {
       continue;
     }
@@ -409,31 +444,86 @@ void Digitizer::storeBC(BCCache& bc,
     if (mCalibOffset) {
       miscalib = mCalibOffset->mTimeOffsets[ipmt];
     }
-    int smeared_time = 1000. * (*cfd.particle - params.mCfdShift) * params.mChannelWidthInverse + miscalib; // + int(1000. * mIntRecord.getTimeOffsetWrtBC() * params.mChannelWidthInverse);
+    int smeared_time = 1000. * (*cfd.particle - params.mCfdShift) * params.mChannelWidthInverse + miscalib;
     bool is_time_in_signal_gate = (smeared_time > -params.mTime_trg_gate && smeared_time < params.mTime_trg_gate);
     float charge = measure_amplitude(channel_times) * params.mCharge2amp;
-    float amp = is_time_in_signal_gate ? params.mMV_2_Nchannels * charge : 0;
-    if (amp > 4095) {
-      amp = 4095;
+    float amp = is_time_in_signal_gate ? params.mMV_2_Nchannels * charge : 0.f;
+    if (amp > 4095.f) {
+      amp = 4095.f;
     }
 //    if (!disabledChannels.count(ipmt)) {
 //      continue;
 //    }
 
-    LOG(debug) << mEventID << " bc " << firstBCinDeque.bc << " orbit " << firstBCinDeque.orbit << ", ipmt " << ipmt << ", smeared_time " << smeared_time << " nStored " << nStored << " offset " << miscalib;
+    LOG(debug) << mEventID << " bc " << firstBCinDeque.bc << " orbit " << firstBCinDeque.orbit
+               << ", ipmt " << ipmt << ", smeared_time " << smeared_time
+               << " nStored " << nStored << " offset " << miscalib
+               << " base amp " << amp;
     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);
+    }
+
+    baseAmp[ipmt] = amp;
+    finalAmp[ipmt] = amp;
+    chTime[ipmt] = smeared_time;
+    chChain[ipmt] = chain;
+    chValid[ipmt] = true;
+  }
+
+  for (const auto& channels : pmtChannelGroups) {
+    for (int localIdx = 0; localIdx < 4; ++localIdx) {
+      const int src = channels[localIdx];
+      if (!chValid[src] || baseAmp[src] <= 0.f) {
+        continue;
       }
+
+      const int nb1 = channels[localNeighbours[localIdx][0]];
+      const int nb2 = channels[localNeighbours[localIdx][1]];
+      const int diag = channels[localNeighbours[localIdx][2]];
+
+      const float directXtalk = baseAmp[src] * params.Cross_Talk_Frac;
+      const float diagXtalk = baseAmp[src] * (params.Cross_Talk_Frac / 3.f);
+
+      finalAmp[nb1] += directXtalk;
+      finalAmp[nb2] += directXtalk;
+      finalAmp[diag] += diagXtalk;
+
+      auto propagateIfEmpty = [&](int dst) {
+        if (!chValid[dst]) {
+          chValid[dst] = true;
+          chTime[dst] = chTime[src];
+          chChain[dst] = chChain[src];
+        }
+      };
+
+      propagateIfEmpty(nb1);
+      propagateIfEmpty(nb2);
+      propagateIfEmpty(diag);
+    }
+  }
+
+  for (Int_t ipmt = 0; ipmt < params.mMCPs; ++ipmt) {
+    if (!chValid[ipmt]) {
+      continue;
     }
+
+    float amp = finalAmp[ipmt];
+    if (amp > 4095.f) {
+      amp = 4095.f;
+    }
+
+    const int smeared_time = chTime[ipmt];
+    const int chain = chChain[ipmt];
+    const bool is_time_in_signal_gate = (smeared_time > -params.mTime_trg_gate && smeared_time < params.mTime_trg_gate);
+
+    if (is_time_in_signal_gate && 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++;
 
-    // fill triggers
-
     Bool_t is_A_side = (ipmt < 4 * mGeometry.NCellsA);
     if (!is_time_in_signal_gate) {
       continue;