From 940d4f23e059997f5cd5ed58aa9705458c584080 Mon Sep 17 00:00:00 2001 From: Michal Broz Date: Fri, 28 Feb 2020 09:36:18 -0500 Subject: [PATCH 1/2] Optimization of the Digitizer --- .../FDDSimulation/DigitizationParameters.h | 33 ++-- .../include/FDDSimulation/Digitizer.h | 48 ++++-- .../FIT/FDD/simulation/src/Digitizer.cxx | 156 ++++++++++++------ 3 files changed, 157 insertions(+), 80 deletions(-) diff --git a/Detectors/FIT/FDD/simulation/include/FDDSimulation/DigitizationParameters.h b/Detectors/FIT/FDD/simulation/include/FDDSimulation/DigitizationParameters.h index be77c52abc4ef..9bdcd0da09af3 100644 --- a/Detectors/FIT/FDD/simulation/include/FDDSimulation/DigitizationParameters.h +++ b/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 diff --git a/Detectors/FIT/FDD/simulation/include/FDDSimulation/Digitizer.h b/Detectors/FIT/FDD/simulation/include/FDDSimulation/Digitizer.h index ec38e6a5bcf8d..62a61348952dc 100644 --- a/Detectors/FIT/FDD/simulation/include/FDDSimulation/Digitizer.h +++ b/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 HitRandomRingType; + typedef math_utils::RandomRing 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* 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,8 +56,6 @@ 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(); @@ -59,24 +63,34 @@ class Digitizer void setMCLabels(o2::dataformats::MCTruthContainer* 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* mMCLabels = nullptr; - std::vector> mTime; // Charge time series aka analogue signal pulse from PM - std::vector 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 mPMResponse; // function which describes the PM time response - std::unique_ptr mSinglePhESpectrum; // function which describes the single ph.e. PM response - std::unique_ptr mSignalShape; + std::array, DigitizationParameters::mNchannels> mTime; // Charge time series aka analogue signal pulse from PM + std::vector 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, 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 diff --git a/Detectors/FIT/FDD/simulation/src/Digitizer.cxx b/Detectors/FIT/FDD/simulation/src/Digitizer.cxx index e233937e20387..a601481aaffc4 100644 --- a/Detectors/FIT/FDD/simulation/src/Digitizer.cxx +++ b/Detectors/FIT/FDD/simulation/src/Digitizer.cxx @@ -14,11 +14,11 @@ #include "TMath.h" #include "TRandom.h" -#include "MathUtils/CachingTF1.h" #include #include #include +using namespace o2::math_utils; using namespace o2::fdd; ClassImp(Digitizer); @@ -39,9 +39,11 @@ void Digitizer::process(const std::vector* 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* 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 = "<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 = "<* 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,6 +154,7 @@ 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 @@ -163,14 +162,75 @@ void Digitizer::init() mTime[i].resize(mNBins); mTimeCFD.resize(mNBins); - if (!mPMResponse) - mPMResponse = std::make_unique("mPMResponse", this, &Digitizer::PMResponse, -parameters.mPMTransitTime, 2. * parameters.mPMTransitTime, 0); - if (!mSinglePhESpectrum) - mSinglePhESpectrum = std::make_unique("mSinglePhESpectrum", this, &Digitizer::SinglePhESpectrum, 0, 20, 0); - if (!mSignalShape) { - mSignalShape = std::make_unique("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)); +} From 96b952dd8d05d57e52f5fc26c9ea7947a48dc5e9 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 2 Mar 2020 09:15:47 +0100 Subject: [PATCH 2/2] Update Digitizer.cxx --- Detectors/FIT/FDD/simulation/src/Digitizer.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Detectors/FIT/FDD/simulation/src/Digitizer.cxx b/Detectors/FIT/FDD/simulation/src/Digitizer.cxx index a601481aaffc4..549fdc1129764 100644 --- a/Detectors/FIT/FDD/simulation/src/Digitizer.cxx +++ b/Detectors/FIT/FDD/simulation/src/Digitizer.cxx @@ -64,8 +64,8 @@ void Digitizer::process(const std::vector* hits, o2::fdd::Digit* d 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; + // 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]; @@ -116,7 +116,7 @@ void Digitizer::process(const std::vector* 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; + // LOG(INFO) << "ADC " << channel_data[ipmt].mChargeADC << " Time " << channel_data[ipmt].mTime; } } //_____________________________________________________________________________