Skip to content

Optimization of the Digitizer #3047

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sawenzel merged 2 commits into from
Mar 2, 2020
Merged

Optimization of the Digitizer #3047

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
940d4f2
Optimization of the Digitizer
mbroz84 Feb 28, 2020
96b952d
Update Digitizer.cxx
sawenzel Mar 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 18 additions & 15 deletions Detectors/FIT/FDD/simulation/include/FDDSimulation/DigitizationParameters.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
48 changes: 31 additions & 17 deletions Detectors/FIT/FDD/simulation/include/FDDSimulation/Digitizer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
156 changes: 108 additions & 48 deletions Detectors/FIT/FDD/simulation/src/Digitizer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);

Expand All @@ -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
Expand All @@ -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;
}
}
//_____________________________________________________________________________
Expand All @@ -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)
{
}
Expand All @@ -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));
}