Digitizer.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#include "FT0Simulation/Digitizer.h"
#include "FT0Simulation/DigitizationConstants.h"
#include "FT0Base/FT0DigParam.h"
#include "SimulationDataFormat/MCTruthContainer.h"
#include "CommonConstants/PhysicsConstants.h"
#include "CommonDataFormat/InteractionRecord.h"
 
#include "DataFormatsFT0/LookUpTable.h"
#include "FT0Base/Constants.h"
#include <map>
#include <array>
#include <unordered_map>
#include <regex>
#include <string>
 
#include "TMath.h"
#include "TRandom.h"
#include <TH1F.h>
#include <algorithm>
#include <cassert>
#include <iostream>
#include <optional>
#include <Vc/Vc>
 
using namespace o2::ft0;
 
ClassImp(Digitizer);
 
namespace o2::ft0
{
// signal shape function
template <typename Float>
Float signalForm_i(Float x)
{
  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)
{
  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
  }
 
  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)
{
  auto const mask = (x >= 0.0f);
  Vc::float_v arg(0);
  arg.assign(x, mask); // branchless if
  Vc::float_v const a(-0.45458f);
  Vc::float_v const b(-0.83344945f);
  Vc::float_v result = -(Vc::exp(b * arg) / b - Vc::exp(a * arg) / a) / 7.8446501f;
  return result;
};
} // namespace o2::ft0
 
Digitizer::CFDOutput Digitizer::get_time(const std::vector<float>& times, float deadTime)
{
  assert(std::is_sorted(std::begin(times), std::end(times)));
 
  // get a new batch of random values
  for (float& n : mNoiseSamples) {
    n = mRndGaus.getNextValue();
  }
 
  // this lambda function evaluates the signal shape at a give time
  auto value_at = [&](float time) {
    float val = 0;
    // (1) sum over individual hits
    Vc::float_v acc(0);
    Vc::float_v tableVal(0);
    const float* tp = times.data();
    size_t m = times.size() / Vc::float_v::size();
    for (size_t i = 0; i < m; ++i) {
      tableVal.load(tp);
      tp += Vc::float_v::size();
      Vc::prefetchForOneRead(tp);
      acc += signalFormVc(time - tableVal);
    }
    val += acc.sum();
    // non-SIMD tail
    for (size_t i = Vc::float_v::size() * m; i < times.size(); ++i, ++tp) {
      val += signalForm(time - (*tp));
    }
    // (2) add noise
    // find the right indices into the sinc table
    int timeIndex = std::lround(time / FT0DigParam::Instance().mNoisePeriod * mSincTable.size());
    int timeOffset = timeIndex / mSincTable.size();
    timeIndex %= mSincTable.size();
    if (timeOffset >= mNumNoiseSamples) { // this happens when time >= 25 ns
      timeOffset = mNumNoiseSamples - 1;
      LOG(debug) << "timeOffset >= mNumNoiseSamples";
    }
    if (timeOffset <= -mNumNoiseSamples) { // this happens when time <= -25 ns
      timeOffset = -mNumNoiseSamples + 1;
      LOG(debug) << "timeOffset <= -mNumNoiseSamples";
    }
    Vc::float_v noiseVal(0);
    const float* np = mNoiseSamples.data();
    tp = mSincTable[timeIndex].data() + mNumNoiseSamples - timeOffset;
    acc = 0.0f;
    m = mNumNoiseSamples / Vc::float_v::Size;
    for (size_t i = 0; i < m; ++i) {
      tableVal.load(tp);
      tp += Vc::float_v::Size;
      Vc::prefetchForOneRead(tp);
      noiseVal.load(np);
      np += Vc::float_v::Size;
      Vc::prefetchForOneRead(np);
      acc += noiseVal * tableVal;
    }
    val += acc.sum(); // horizontal sum
    // non-SIMD tail
    for (size_t i = Vc::float_v::Size * m; i < mNumNoiseSamples; ++i, ++tp, ++np) {
      val += (*np) * (*tp);
    }
    return val;
  };
  auto const min_time = std::max(deadTime, *std::min_element(std::begin(times),
                                                             std::end(times)));
  CFDOutput result{std::nullopt, -0.5f * FT0DigParam::Instance().mBunchWidth};
  bool is_positive = true;
 
  // reset the chache
  std::fill_n(std::begin(mSignalCache), std::size(mSignalCache), -1.0f);
  const auto& params = FT0DigParam::Instance();
  // we need double precision for time in order to match previous behaviour
  for (double time = min_time; time < 0.5 * params.mBunchWidth; time += DP::SIGNAL_CACHE_DT) {
    float const val = value_at(time);
    int const index = std::lround((time + 0.5 * params.mBunchWidth) / DP::SIGNAL_CACHE_DT);
    if (index >= 0 && index < mSignalCache.size()) { // save the value for later use
      mSignalCache[index] = val;
    }
    // look up the time-shifted signal value from the past
    float val_prev = 0.0f;
    int const index_prev = std::lround((time - params.mCFDShiftPos + 0.5f * params.mBunchWidth) / DP::SIGNAL_CACHE_DT);
    val_prev = ((index_prev < 0 || index_prev >= mSignalCache.size() || mSignalCache[index_prev] < 0.0f)
                  ? value_at(time - params.mCFDShiftPos) //  was not computed before
                  : mSignalCache[index_prev]);           //  is available in the cache
    float const cfd_val = 5.0f * val_prev - val;
    if (std::abs(val) > params.mCFD_trsh && !is_positive && cfd_val > 0.0f) {
      if (!result.particle) {
        result.particle = time;
      }
      result.deadTime = time + params.mCFDdeadTime;
      time += params.mCFDdeadTime - DP::SIGNAL_CACHE_DT;
      is_positive = true;
    } else {
      is_positive = cfd_val > 0.0f;
    }
  }
  return result;
}
 
double Digitizer::measure_amplitude(const std::vector<float>& times) const
{
  float const from = FT0DigParam::Instance().mAmpRecordLow;
  float const to = from + FT0DigParam::Instance().mAmpRecordUp;
  // SIMD version has a negligible effect on the total wall time
  Vc::float_v acc(0);
  Vc::float_v tv(0);
  const float* tp = times.data();
  size_t const m = times.size() / Vc::float_v::Size;
  for (size_t i = 0; i < m; ++i) {
    tv.load(tp);
    tp += Vc::float_v::Size;
    Vc::prefetchForOneRead(tp);
    acc += signalForm_integralVc(to - tv) - signalForm_integralVc(from - tv);
  }
  float result = acc.sum(); // horizontal sum
  // non-SIMD tail
  for (size_t i = Vc::float_v::Size * m; i < times.size(); ++i, ++tp) {
    result += signalForm_integral(to - (*tp)) - signalForm_integral(from - (*tp));
  }
  return result;
}
 
void Digitizer::process(const std::vector<o2::ft0::HitType>* hits,
                        std::vector<o2::ft0::Digit>& digitsBC,
                        std::vector<o2::ft0::ChannelData>& digitsCh,
                        std::vector<o2::ft0::DetTrigInput>& digitsTrig,
                        o2::dataformats::MCTruthContainer<o2::ft0::MCLabel>& label)
{
  ;
  // Calculating signal time, amplitude in mean_time +- time_gate --------------
  LOG(debug) << " process firstBCinDeque " << firstBCinDeque << " mIntRecord " << mIntRecord;
  if (firstBCinDeque != mIntRecord) {
    flush(digitsBC, digitsCh, digitsTrig, label);
  }
 
  Int_t parent = -10;
  for (auto const& hit : *hits) {
    if (hit.GetEnergyLoss() > 0) {
      continue;
    }
 
    Int_t hit_ch = hit.GetDetectorID();
 
    // If the dead channel map is used, and the channel with ID 'hit_ch' is dead, don't process this hit.
    if (mDeadChannelMap && !mDeadChannelMap->isChannelAlive(hit_ch)) {
      continue;
    }
 
    const auto& params = FT0DigParam::Instance();
 
    Bool_t is_A_side = (hit_ch < 4 * mGeometry.NCellsA);
 
    // Subtract time-of-flight from hit time
    const Float_t timeOfFlight = hit.GetPos().R() / o2::constants::physics::LightSpeedCm2NS;
    const Float_t timeOffset = is_A_side ? params.hitTimeOffsetA : params.hitTimeOffsetC;
    Double_t hit_time = hit.GetTime() - timeOfFlight + timeOffset + mIntRecord.getTimeOffsetWrtBC();
 
    if (hit_time > 150) {
      continue; // not collect very slow particles
    }
 
    auto relBC = o2::InteractionRecord{hit_time};
    if (mCache.size() <= relBC.bc) {
      mCache.resize(relBC.bc + 1);
    }
    mCache[relBC.bc].hits.emplace_back(BCCache::particle{hit_ch, hit_time - relBC.bc2ns()});
    // charge particles in MCLabel
    Int_t parentID = hit.GetTrackID();
    if (parentID != parent) {
      mCache[relBC.bc].labels.emplace(parentID, mEventID, mSrcID, hit_ch);
      parent = parentID;
    }
  }
}
 
void Digitizer::storeBC(BCCache& bc,
                        std::vector<o2::ft0::Digit>& digitsBC,
                        std::vector<o2::ft0::ChannelData>& digitsCh,
                        std::vector<o2::ft0::DetTrigInput>& digitsTrig,
                        o2::dataformats::MCTruthContainer<o2::ft0::MCLabel>& labels)
{
  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();
 
  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::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; });
 
    if (channel_end - channel_begin < params.mAmp_trsh) {
      continue;
    }
    channel_times.resize(channel_end - channel_begin);
    std::transform(channel_begin, channel_end, channel_times.begin(), [](BCCache::particle const& p) { return p.hit_time; });
    int chain = (std::rand() % 2) ? 1 : 0;
    auto cfd = get_time(channel_times, mDeadTimes[ipmt].intrec.bc2ns() -
                                         firstBCinDeque.bc2ns() +
                                         mDeadTimes[ipmt].deadTime);
    mDeadTimes[ipmt].intrec = firstBCinDeque;
    mDeadTimes[ipmt].deadTime = cfd.deadTime;
 
    if (!cfd.particle) {
      continue;
    }
    // miscalibrate CFD with cahnnel offsets
    int miscalib = 0;
    if (mCalibOffset) {
      miscalib = mCalibOffset->mTimeOffsets[ipmt];
    }
    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.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
               << " base amp " << amp;
    if (is_time_in_signal_gate) {
      chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsCFDinADCgate);
      chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsEventInTVDC);
    }
 
    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;
 
      if (!chValid[nb1] && directXtalk >= params.mAmpThresholdForCrossTalkDigit) {
        chValid[nb1] = true;
        chTime[nb1] = chTime[src];
        chChain[nb1] = chChain[src];
      }
 
      if (!chValid[nb2] && directXtalk >= params.mAmpThresholdForCrossTalkDigit) {
        chValid[nb2] = true;
        chTime[nb2] = chTime[src];
        chChain[nb2] = chChain[src];
      }
 
      if (!chValid[diag] && diagXtalk >= params.mAmpThresholdForCrossTalkDigit) {
        chValid[diag] = true;
        chTime[diag] = chTime[src];
        chChain[diag] = chChain[src];
      }
    }
  }
 
  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 bool hasPrimarySignal = (baseAmp[ipmt] > 0.f);
    const bool isCrossTalkOnly = (!hasPrimarySignal && amp > 0.f);
 
    if (isCrossTalkOnly && amp < params.mAmpThresholdForCrossTalkDigit) {
      continue;
    }
 
    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++;
 
    Bool_t is_A_side = (ipmt < 4 * mGeometry.NCellsA);
    if (!is_time_in_signal_gate) {
      continue;
    }
 
    sum_ampl_ipmt[ipmt] += amp;
 
    if (is_A_side) {
      n_hit_A++;
      summ_ampl_A += amp;
      mean_time_A += smeared_time;
    } else {
      n_hit_C++;
      summ_ampl_C += amp;
      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 = 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_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;
  bool isOutputsAreBlocked = false;
  bool isDataValid = true;
  if (nStored > 0) {
    triggers.setTriggers(is_A, is_C, isVertex, is_Central, is_SemiCentral, int8_t(n_hit_A), int8_t(n_hit_C),
                         amplA, amplC, timeA, timeC, isLaser, isOutputsAreBlocked, isDataValid);
    digitsBC.emplace_back(first, nStored, firstBCinDeque, triggers, mEventID - 1);
    digitsTrig.emplace_back(firstBCinDeque, is_A, is_C, isVertex, is_Central, is_SemiCentral);
    size_t const nBC = digitsBC.size();
    for (auto const& lbl : bc.labels) {
      labels.addElement(nBC - 1, lbl);
    }
  }
  // Debug output -------------------------------------------------------------
 
  LOG(debug) << "Event ID: " << mEventID << ", bc " << firstBCinDeque.bc << ", N hit " << bc.hits.size();
  LOG(debug) << "N hit A: " << int(triggers.getNChanA()) << " N hit C: " << int(triggers.getNChanC()) << " summ ampl A: " << int(triggers.getAmplA())
             << " summ ampl C: " << int(triggers.getAmplC()) << " mean time A: " << triggers.getTimeA()
             << " mean time C: " << triggers.getTimeC() << " nStored " << nStored;
 
  LOG(debug) << "IS A " << triggers.getOrA() << " IsC " << triggers.getOrC() << " vertex " << triggers.getVertex() << " is Central " << triggers.getCen() << " is SemiCentral " << triggers.getSCen();
}
 
//------------------------------------------------------------------------
void Digitizer::flush(std::vector<o2::ft0::Digit>& digitsBC,
                      std::vector<o2::ft0::ChannelData>& digitsCh,
                      std::vector<o2::ft0::DetTrigInput>& digitsTrig,
                      o2::dataformats::MCTruthContainer<o2::ft0::MCLabel>& labels)
{
 
  assert(firstBCinDeque <= mIntRecord);
 
  while (firstBCinDeque < mIntRecord && !mCache.empty()) {
    storeBC(mCache.front(), digitsBC, digitsCh, digitsTrig, labels);
    mCache.pop_front();
    ++firstBCinDeque;
  }
  firstBCinDeque = mIntRecord;
}
 
void Digitizer::flush_all(std::vector<o2::ft0::Digit>& digitsBC,
                          std::vector<o2::ft0::ChannelData>& digitsCh,
                          std::vector<o2::ft0::DetTrigInput>& digitsTrig,
                          o2::dataformats::MCTruthContainer<o2::ft0::MCLabel>& labels)
{
 
  assert(firstBCinDeque <= mIntRecord);
  ++mEventID;
  while (!mCache.empty()) {
    storeBC(mCache.front(), digitsBC, digitsCh, digitsTrig, labels);
    mCache.pop_front();
    ++firstBCinDeque;
  }
}
 
void Digitizer::initParameters()
{
  auto const sinc = [](double x) { x *= TMath::Pi(); return (std::abs(x) < 1e-12) ? 1.0 : std::sin(x) / x; };
 
  // number of noise samples in one BC
  const auto& params = FT0DigParam::Instance();
  mNumNoiseSamples = std::ceil(params.mBunchWidth / params.mNoisePeriod);
  mNoiseSamples.resize(mNumNoiseSamples);
 
  // set up tables with sinc function values (times noiseVar)
  for (size_t i = 0, n = mSincTable.size(); i < n; ++i) {
    float const time = i / float(n) * params.mNoisePeriod; // [0 .. 1/params.mNoisePeriod)
    LOG(debug) << "initParameters " << i << "/" << n << " " << time;
    // we make a table of sinc values between -num_noise_samples and 2*num_noise_samples
    mSincTable[i].resize(3 * mNumNoiseSamples);
    for (int j = -mNumNoiseSamples; j < 2 * mNumNoiseSamples; ++j) {
      mSincTable[i][mNumNoiseSamples + j] = params.mNoiseVar * sinc((time + 0.5f * params.mBunchWidth) / params.mNoisePeriod - j);
    }
  }
  // set up the lookup table for the signal form
  for (size_t i = 0; i < DP::SIGNAL_TABLE_SIZE; ++i) {
    float const x = float(i) / float(DP::SIGNAL_TABLE_SIZE) * params.mBunchWidth;
    mSignalTable[i] = signalForm_i(x);
  }
 
  // cache for signal time series used by the CFD -BC/2 .. +3BC/2
  mSignalCache.resize(std::lround(params.mBunchWidth / DP::SIGNAL_CACHE_DT));
}
//_______________________________________________________________________
void Digitizer::init()
{
  LOG(info) << " @@@ Digitizer::init " << std::endl;
  mDeadTimes.fill({InteractionRecord(0), -100.});
  printParameters();
}
//_______________________________________________________________________
void Digitizer::finish()
{
  printParameters();
}
 
//_______________________________________________________________________
void Digitizer::printParameters() const
{
  const auto& params = FT0DigParam::Instance();
  LOG(info) << " Run Digitzation with parametrs: \n"
            << " CFD amplitude threshold \n " << params.mCFD_trsh << " CFD signal gate in ps \n"
            << params.mTime_trg_gate << "shift to have signal around zero after CFD trancformation  \n"
            << params.mCfdShift << "CFD distance between 0.3 of max amplitude  to max \n"
            << params.mCFDShiftPos << "MIP -> mV " << params.mMip_in_V << " Pe in MIP \n"
            << params.mPe_in_mip << "noise level " << params.mNoiseVar << " noise frequency \n"
            << params.mNoisePeriod << " mMCPs " << params.mMCPs;
}
 
O2ParamImpl(FT0DigParam);