FDD: Optimization of the Digitizer (#3047)

* Optimization of the Digitizer (applying techniques already done for FV0)

Author: mbroz84

Detectors/FIT/FDD/simulation/include/FDDSimulation/DigitizationParameters.h

@@ -14,21 +14,24 @@
 namespace o2::fdd
 {
 struct DigitizationParameters {
-  const UShort_t mNchannels = 16;
-  const UShort_t mNtriggers = 5;
-  const Float_t mIntTimeRes = 0.4;
-  const Float_t mPhotoCathodeEfficiency = 0.18;
-  const Float_t mLightYield = 0.01;
-  const Float_t mPmGain = 1e6;
-  const Float_t mChargePerADC = 0.6e-12;
-  const Float_t mPMTransitTime = 6.0;   // PM response time (corresponds to 1.9 ns rise time)
-  const Float_t mPMTransparency = 0.25; // Transparency of the first dynode of the PM
-  const Float_t mPMNbOfSecElec = 6.0;   // Number of secondary electrons emitted from first dynode (per ph.e.)
-  const Float_t mShapeAlpha = -0.445;
-  const Float_t mShapeN = 2.65;
-  const Float_t mShapeSigma = 3.25;
-  //const Float_t mPedestal = 0;
-  const Float_t mTimeShiftCFD = 1.42;
+  static constexpr UShort_t mNchannels = 16;
+  static constexpr UShort_t mNtriggers = 5;
+  static constexpr Float_t mIntTimeRes = 0.4;
+  static constexpr Float_t mPhotoCathodeEfficiency = 0.18;
+  static constexpr Float_t mLightYield = 0.01;
+  static constexpr Float_t mPmGain = 1e6;
+  static constexpr Float_t mChargePerADC = 0.6e-12;
+  static constexpr Float_t mPMTransitTime = 6.0;   // PM response time (corresponds to 1.9 ns rise time)
+  static constexpr Float_t mPMTransparency = 0.25; // Transparency of the first dynode of the PM
+  static constexpr Float_t mPMNbOfSecElec = 6.0;   // Number of secondary electrons emitted from first dynode (per ph.e.)
+  static constexpr Float_t mShapeAlpha = -0.445;
+  static constexpr Float_t mShapeN = 2.65;
+  static constexpr Float_t mShapeSigma = 3.25;
+  //static constexpr Float_t mPedestal = 0;
+  static constexpr Float_t mTimeShiftCFD = 1.42;
+  static constexpr int mPheRRSize = 1e5;     // size of random ring to be used inside photoelectron loop
+  static constexpr int mHitRRSize = 1e4;     // size of random ring to be used inside hit loop
+  static constexpr int mNResponseTables = 9; // number of PMT response tables
 };
 } // namespace o2::fdd
 #endif

Detectors/FIT/FDD/simulation/include/FDDSimulation/Digitizer.h

@@ -18,6 +18,7 @@
 #include "SimulationDataFormat/MCTruthContainer.h"
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "FDDSimulation/DigitizationParameters.h"
+#include "MathUtils/RandomRing.h"
 #include "MathUtils/CachingTF1.h"
 
 namespace o2
@@ -26,14 +27,19 @@ namespace fdd
 {
 class Digitizer
 {
+
+ private:
+  typedef math_utils::RandomRing<float_v::size() * DigitizationParameters::mPheRRSize> HitRandomRingType;
+  typedef math_utils::RandomRing<float_v::size() * DigitizationParameters::mHitRRSize> PheRandomRingType;
+
  public:
-  Digitizer(const DigitizationParameters& params, Int_t mode = 0) : mMode(mode), parameters(params), mTime(16) { init(); };
+  Digitizer(const DigitizationParameters& params, Int_t mode = 0) : mEventTime(0), mIntRecord(), mEventID(-1), mSrcID(-1), mMCLabels(), parameters(params), mTime(), mRndScintDelay(HitRandomRingType::RandomType::CustomTF1), mRndGainVar(PheRandomRingType::RandomType::CustomTF1), mRndSignalShape(PheRandomRingType::RandomType::CustomTF1), mPMResponseTables() { init(); };
   ~Digitizer() = default;
 
   void process(const std::vector<o2::fdd::Hit>* hits, o2::fdd::Digit* digit);
 
   void initParameters();
-  void SetEventTime(double value) { mEventTime = value; }
+  void SetEventTime(long value) { mEventTime = value; }
   void SetEventID(Int_t id) { mEventID = id; }
   void SetSrcID(Int_t id) { mSrcID = id; }
   void SetInteractionRecord(uint16_t bc, uint32_t orbit)
@@ -50,33 +56,41 @@ class Digitizer
   void SetTriggers(o2::fdd::Digit* digit);
   Int_t SimulateLightYield(Int_t pmt, Int_t nPhot);
   Float_t SimulateTimeCFD(Int_t channel);
-  Double_t PMResponse(Double_t* x, Double_t* par);
-  Double_t SinglePhESpectrum(Double_t* x, Double_t* par);
 
   void init();
   void finish();
 
   void setMCLabels(o2::dataformats::MCTruthContainer<o2::fdd::MCLabel>* mclb) { mMCLabels = mclb; }
 
  private:
-  Int_t mMode;                      //triggered or continuos
-  o2::InteractionRecord mIntRecord; // Interaction record (orbit, bc)
-  Int_t mEventID;
-  Int_t mSrcID;        // signal, background or QED
-  Double_t mEventTime; // timestamp
+  long mEventTime;              // TF (run) timestamp
+  InteractionRecord mIntRecord; // Interaction record (orbit, bc) -> InteractionTimeRecord
+  Int_t mEventID;               // ID of the current event
+  Int_t mSrcID;                 // signal, background or QED
 
   DigitizationParameters parameters;
   o2::dataformats::MCTruthContainer<o2::fdd::MCLabel>* mMCLabels = nullptr;
 
-  std::vector<std::vector<Float_t>> mTime;                  // Charge time series aka analogue signal pulse from PM
-  std::vector<Float_t> mTimeCFD;                            // Time series for CFD measurement
-  UInt_t mNBins;                                            // Number of bins in pulse series
-  Float_t mBinSize;                                         // Time width of the pulse bin - HPTDC resolution
-  std::unique_ptr<o2::base::CachingTF1> mPMResponse;        // function which describes the PM time response
-  std::unique_ptr<o2::base::CachingTF1> mSinglePhESpectrum; // function which describes the single ph.e. PM response
-  std::unique_ptr<o2::base::CachingTF1> mSignalShape;
+  std::array<std::vector<Float_t>, DigitizationParameters::mNchannels> mTime; // Charge time series aka analogue signal pulse from PM
+  std::vector<Float_t> mTimeCFD;                                              // Time series for CFD measurement
+  UInt_t mNBins;                                                              // Number of bins in pulse series
+  Float_t mBinSize;                                                           // Time width of the pulse bin - HPTDC resolution
+  Float_t mPmtTimeIntegral;
+
+  // Random rings
+  HitRandomRingType mRndScintDelay;
+  PheRandomRingType mRndGainVar;
+  PheRandomRingType mRndSignalShape;
+
+  // 8 tables starting at different sub-bin positions, i.e, [-4:4] / 8 * mBinSize
+  // wit each table containg values for start + [-2:2:mBinSize] * DigitizationParameters::mPmtTransitTime
+  std::array<std::vector<Float_t>, DigitizationParameters::mNResponseTables> mPMResponseTables;
+
+  static Double_t PMResponse(Double_t x);
+  static Double_t PMResponse(Double_t* x, Double_t*);
+  static Double_t SinglePhESpectrum(Double_t* x, Double_t* par);
 
-  ClassDefNV(Digitizer, 2);
+  ClassDefNV(Digitizer, 3);
 };
 } // namespace fdd
 } // namespace o2

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

@@ -14,11 +14,11 @@
 
 #include "TMath.h"
 #include "TRandom.h"
-#include "MathUtils/CachingTF1.h"
 #include <algorithm>
 #include <cassert>
 #include <iostream>
 
+using namespace o2::math_utils;
 using namespace o2::fdd;
 
 ClassImp(Digitizer);
@@ -39,9 +39,11 @@ void Digitizer::process(const std::vector<o2::fdd::Hit>* hits, o2::fdd::Digit* d
     for (int i = 0; i < parameters.mNchannels; ++i)
       channel_data.emplace_back(o2::fdd::ChannelData{i, o2::InteractionRecord::DummyTime, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0});
   }
+
+  auto const roundVc = [&](int i) -> int {
+    return (i / Vc::float_v::Size) * Vc::float_v::Size;
+  };
   Int_t parent = -10;
-  Float_t integral = mPMResponse->Integral(-parameters.mPMTransitTime, 2. * parameters.mPMTransitTime);
-  Float_t meansPhE = mSinglePhESpectrum->Mean(0, 20);
   for (Int_t i = 0; i < parameters.mNchannels; i++)
     std::fill(mTime[i].begin(), mTime[i].end(), 0);
 
@@ -60,22 +62,43 @@ void Digitizer::process(const std::vector<o2::fdd::Hit>* hits, o2::fdd::Digit* d
     Int_t pmt = hit.GetDetectorID();
     Int_t nPhE = SimulateLightYield(pmt, hit.GetNphot());
 
-    Float_t dt_scintillator = gRandom->Gaus(0, parameters.mIntTimeRes);
+    Float_t dt_scintillator = mRndScintDelay.getNextValue();
     Float_t t = dt_scintillator + hit.GetTime();
+    // LOG(INFO) << "Nphot = " << hit.GetNphot() << " time =" << hit.GetTime();
+    // LOG(INFO) << "NphE = " << nPhE;
+    Float_t charge = TMath::Qe() * parameters.mPmGain * mBinSize / mPmtTimeIntegral;
+
+    auto& analogSignal = mTime[pmt];
 
-    //LOG(INFO) << "Nphot = "<<hit.GetNphot()<<" time ="<<hit.GetTime()<<FairLogger::endl;
-    //LOG(INFO) << "NphE = " << nPhE << FairLogger::endl;
-    Float_t charge = TMath::Qe() * parameters.mPmGain * mBinSize / integral;
     for (Int_t iPhE = 0; iPhE < nPhE; ++iPhE) {
-      Float_t tPhE = t + mSignalShape->GetRandom(0, mBinSize * Float_t(mNBins));
+      Float_t tPhE = t + mRndSignalShape.getNextValue();
       //LOG(INFO) <<"t = "<<t<<"tPhE = "<<tPhE;
-      Float_t gainVar = mSinglePhESpectrum->GetRandom(0, 20) / meansPhE;
-      Int_t firstBin = TMath::Max((UInt_t)0, (UInt_t)((tPhE - parameters.mPMTransitTime) / mBinSize));
-      Int_t lastBin = TMath::Min(mNBins - 1, (UInt_t)((tPhE + 2. * parameters.mPMTransitTime) / mBinSize));
-      //LOG(INFO) << "firstBin = "<<firstBin<<" lastbin "<<lastBin<<FairLogger::endl;
-      for (Int_t iBin = firstBin; iBin <= lastBin; ++iBin) {
-        Float_t tempT = mBinSize * (0.5 + iBin) - tPhE;
-        mTime[pmt][iBin] += gainVar * charge * mPMResponse->Eval(tempT);
+      Int_t const firstBin = roundVc(TMath::Max((Int_t)0, (Int_t)((tPhE - parameters.mPMTransitTime) / mBinSize)));
+      Int_t const lastBin = TMath::Min((Int_t)mNBins - 1, (Int_t)((tPhE + 2. * parameters.mPMTransitTime) / mBinSize));
+      //LOG(INFO) << "firstBin = "<<firstBin<<" lastbin "<<lastBin;
+      Float_t const tempT = mBinSize * (0.5f + firstBin) - tPhE;
+      Float_t* p = analogSignal.data() + firstBin;
+      long iStart = std::lround((tempT + 2.0f * parameters.mPMTransitTime) / mBinSize);
+      float const offset = tempT + 2.0f * parameters.mPMTransitTime - Float_t(iStart) * mBinSize;
+      long const iOffset = std::lround(offset / mBinSize * Float_t(parameters.mNResponseTables - 1));
+      if (iStart < 0) { // this should not happen
+        LOG(ERROR) << "FDDDigitizer: table lookup failure";
+      }
+      iStart = roundVc(std::max(long(0), iStart));
+
+      Vc::float_v workVc;
+      Vc::float_v pmtVc;
+      Float_t const* q = mPMResponseTables[parameters.mNResponseTables / 2 + iOffset].data() + iStart;
+      Float_t const* qEnd = &mPMResponseTables[parameters.mNResponseTables / 2 + iOffset].back();
+      for (Int_t i = firstBin, iEnd = roundVc(lastBin); q < qEnd && i < iEnd; i += Vc::float_v::Size) {
+        pmtVc.load(q);
+        q += Vc::float_v::Size;
+        Vc::prefetchForOneRead(q);
+        workVc.load(p);
+        workVc += mRndGainVar.getNextValueVc() * charge * pmtVc;
+        workVc.store(p);
+        p += Vc::float_v::Size;
+        Vc::prefetchForOneRead(p);
       }
     }
     //MCLabels
@@ -93,7 +116,7 @@ void Digitizer::process(const std::vector<o2::fdd::Hit>* hits, o2::fdd::Digit* d
     channel_data[ipmt].mTime = SimulateTimeCFD(ipmt);
     for (Int_t iBin = 0; iBin < mNBins; ++iBin)
       channel_data[ipmt].mChargeADC += mTime[ipmt][iBin] / parameters.mChargePerADC;
-    //LOG(INFO) << "ADC " << channel_data[ipmt].mChargeADC << " Time " << channel_data[ipmt].mTime << FairLogger::endl;
+    // LOG(INFO) << "ADC " << channel_data[ipmt].mChargeADC << " Time " << channel_data[ipmt].mTime;
   }
 }
 //_____________________________________________________________________________
@@ -120,31 +143,6 @@ Float_t Digitizer::SimulateTimeCFD(Int_t channel)
   return timeCFD;
 }
 //_____________________________________________________________________________
-Int_t Digitizer::SimulateLightYield(Int_t pmt, Int_t nPhot)
-{
-  const Float_t p = parameters.mLightYield * parameters.mPhotoCathodeEfficiency;
-  if (p == 1.0f || nPhot == 0)
-    return nPhot;
-  const Int_t n = Int_t(nPhot < 100 ? gRandom->Binomial(nPhot, p) : gRandom->Gaus(p * nPhot + 0.5, TMath::Sqrt(p * (1 - p) * nPhot)));
-  return n;
-}
-//_____________________________________________________________________________
-Double_t Digitizer::PMResponse(Double_t* x, Double_t*)
-{
-  // this function describes the PM time response to a single photoelectron
-  Double_t y = x[0] + parameters.mPMTransitTime;
-  return y * y * TMath::Exp(-y * y / (parameters.mPMTransitTime * parameters.mPMTransitTime));
-}
-//_____________________________________________________________________________
-Double_t Digitizer::SinglePhESpectrum(Double_t* x, Double_t*)
-{
-  // this function describes the PM amplitude response to a single photoelectron
-  Double_t y = x[0];
-  if (y < 0)
-    return 0;
-  return (TMath::Poisson(y, parameters.mPMNbOfSecElec) + parameters.mPMTransparency * TMath::Poisson(y, 1.0));
-}
-//_____________________________________________________________________________
 void Digitizer::SetTriggers(o2::fdd::Digit* digit)
 {
 }
@@ -156,21 +154,83 @@ void Digitizer::initParameters()
 //_______________________________________________________________________
 void Digitizer::init()
 {
+  mEventTime = 0;
 
   mNBins = 2000;           //Will be computed using detector set-up from CDB
   mBinSize = 25.0 / 256.0; //Will be set-up from CDB
   for (Int_t i = 0; i < parameters.mNchannels; i++)
     mTime[i].resize(mNBins);
   mTimeCFD.resize(mNBins);
 
-  if (!mPMResponse)
-    mPMResponse = std::make_unique<o2::base::CachingTF1>("mPMResponse", this, &Digitizer::PMResponse, -parameters.mPMTransitTime, 2. * parameters.mPMTransitTime, 0);
-  if (!mSinglePhESpectrum)
-    mSinglePhESpectrum = std::make_unique<o2::base::CachingTF1>("mSinglePhESpectrum", this, &Digitizer::SinglePhESpectrum, 0, 20, 0);
-  if (!mSignalShape) {
-    mSignalShape = std::make_unique<o2::base::CachingTF1>("mSignalShape", "crystalball", 0, 300);
-    mSignalShape->SetParameters(1, parameters.mShapeSigma, parameters.mShapeSigma, parameters.mShapeAlpha, parameters.mShapeN);
+  auto const roundVc = [&](int i) -> int {
+    return (i / Vc::float_v::Size) * Vc::float_v::Size;
+  };
+  // set up PMT response tables
+  Float_t offset = -0.5f * mBinSize; // offset \in [-0.5..0.5] * mBinSize
+  Int_t const nBins = roundVc(std::lround(4.0f * parameters.mPMTransitTime / mBinSize));
+  for (auto& table : mPMResponseTables) {
+    table.resize(nBins);
+    Float_t t = -2.0f * parameters.mPMTransitTime + offset; // t \in offset + [-2 2] * DP::mPmtTransitTime
+    for (Int_t j = 0; j < nBins; ++j) {
+      table[j] = Digitizer::PMResponse(t);
+      t += mBinSize;
+    }
+    offset += mBinSize / Float_t(parameters.mNResponseTables - 1);
   }
+
+  TF1 scintDelayFn("fScintDelay", "gaus", -6.0f * parameters.mIntTimeRes, +6.0f * parameters.mIntTimeRes);
+  scintDelayFn.SetParameters(1, 0, parameters.mIntTimeRes);
+  mRndScintDelay.initialize(scintDelayFn);
+
+  // Initialize function describing the PMT time response
+  TF1 pmtResponseFn("mPmtResponseFn", &Digitizer::PMResponse, -1.0f * parameters.mPMTransitTime, +2.0f * parameters.mPMTransitTime, 0);
+  pmtResponseFn.SetNpx(100);
+  mPmtTimeIntegral = pmtResponseFn.Integral(-1.0f * parameters.mPMTransitTime, +2.0f * parameters.mPMTransitTime);
+
+  // Initialize function describing PMT response to the single photoelectron
+  TF1 singlePhESpectrumFn("mSinglePhESpectrum",
+                          &Digitizer::SinglePhESpectrum, 0, 30, 0);
+  Float_t const meansPhE = singlePhESpectrumFn.Mean(0, 30);
+  mRndGainVar.initialize([&]() -> float {
+    return singlePhESpectrumFn.GetRandom(0, 30) / meansPhE;
+  });
+
+  TF1 signalShapeFn("signalShape", "crystalball", 0, 300);
+  signalShapeFn.SetParameters(1, parameters.mShapeSigma, parameters.mShapeSigma, parameters.mShapeAlpha, parameters.mShapeN);
+  mRndSignalShape.initialize([&]() -> float {
+    return signalShapeFn.GetRandom(0, mBinSize * Float_t(mNBins));
+  });
 }
 //_______________________________________________________________________
 void Digitizer::finish() {}
+
+//_____________________________________________________________________________
+Int_t Digitizer::SimulateLightYield(Int_t pmt, Int_t nPhot)
+{
+  const Float_t p = parameters.mLightYield * parameters.mPhotoCathodeEfficiency;
+  if (p == 1.0f || nPhot == 0)
+    return nPhot;
+  const Int_t n = Int_t(nPhot < 100 ? gRandom->Binomial(nPhot, p) : gRandom->Gaus(p * nPhot + 0.5, TMath::Sqrt(p * (1 - p) * nPhot)));
+  return n;
+}
+//_____________________________________________________________________________
+Double_t Digitizer::PMResponse(Double_t* x, Double_t*)
+{
+  return Digitizer::PMResponse(x[0]);
+}
+//_____________________________________________________________________________
+Double_t Digitizer::PMResponse(Double_t x)
+{
+  // this function describes the PM time response to a single photoelectron
+  Double_t y = x + DigitizationParameters::mPMTransitTime;
+  return y * y * TMath::Exp(-y * y / (DigitizationParameters::mPMTransitTime * DigitizationParameters::mPMTransitTime));
+}
+//_____________________________________________________________________________
+Double_t Digitizer::SinglePhESpectrum(Double_t* x, Double_t*)
+{
+  // this function describes the PM amplitude response to a single photoelectron
+  Double_t y = x[0];
+  if (y < 0)
+    return 0;
+  return (TMath::Poisson(y, DigitizationParameters::mPMNbOfSecElec) + DigitizationParameters::mPMTransparency * TMath::Poisson(y, 1.0));
+}