39template <
typename Float>
42 float p0,
p1,
p2, p3, p4, p5, p6, p7;
56 Double_t arg1 = (log(
x1) - p0) /
p1;
57 val +=
p2 * (1.0 / (
x1 *
p1 * sqrt(2 * TMath::Pi()))) *
exp(-0.5 * arg1 * arg1);
62 Double_t arg2 = (log(x2) - p4) / p5;
63 val += p6 * (1.0 / (x2 * p5 * sqrt(2 * TMath::Pi()))) *
exp(-0.5 * arg2 * arg2);
72 float p0,
p1,
p2, p3, p4, p5, p6, p7;
85 Double_t z1 = (log(
x1) - p0) / (sqrt(2) *
p1);
86 val +=
p2 * 0.5 * (1 + TMath::Erf(z1));
91 Double_t z2 = (log(x2) - p4) / (sqrt(2) * p5);
92 val += p6 * 0.5 * (1 + TMath::Erf(z2));
123 auto const mask = (
x >= 0.0f);
126 Vc::float_v
const a(-0.45458f);
127 Vc::float_v
const b(-0.83344945f);
128 Vc::float_v
result = -(Vc::exp(
b * arg) /
b - Vc::exp(
a * arg) /
a) / 7.8446501f;
135 assert(std::is_sorted(std::begin(
times), std::end(
times)));
138 for (
float&
n : mNoiseSamples) {
139 n = mRndGaus.getNextValue();
143 auto value_at = [&](
float time) {
147 Vc::float_v tableVal(0);
148 const float* tp =
times.data();
149 size_t m =
times.size() / Vc::float_v::size();
150 for (
size_t i = 0;
i <
m; ++
i) {
152 tp += Vc::float_v::size();
153 Vc::prefetchForOneRead(tp);
158 for (
size_t i = Vc::float_v::size() *
m;
i <
times.size(); ++
i, ++tp) {
164 int timeOffset = timeIndex / mSincTable.size();
165 timeIndex %= mSincTable.size();
166 if (timeOffset >= mNumNoiseSamples) {
167 timeOffset = mNumNoiseSamples - 1;
168 LOG(
debug) <<
"timeOffset >= mNumNoiseSamples";
170 if (timeOffset <= -mNumNoiseSamples) {
171 timeOffset = -mNumNoiseSamples + 1;
172 LOG(
debug) <<
"timeOffset <= -mNumNoiseSamples";
174 Vc::float_v noiseVal(0);
175 const float* np = mNoiseSamples.data();
176 tp = mSincTable[timeIndex].data() + mNumNoiseSamples - timeOffset;
178 m = mNumNoiseSamples / Vc::float_v::Size;
179 for (
size_t i = 0;
i <
m; ++
i) {
181 tp += Vc::float_v::Size;
182 Vc::prefetchForOneRead(tp);
184 np += Vc::float_v::Size;
185 Vc::prefetchForOneRead(np);
186 acc += noiseVal * tableVal;
190 for (
size_t i = Vc::float_v::Size *
m;
i < mNumNoiseSamples; ++
i, ++tp, ++np) {
191 val += (*np) * (*tp);
195 auto const min_time = std::max(deadTime, *std::min_element(std::begin(
times),
198 bool is_positive =
true;
201 std::fill_n(std::begin(mSignalCache), std::size(mSignalCache), -1.0f);
205 float const val = value_at(
time);
207 if (
index >= 0 &&
index < mSignalCache.size()) {
211 float val_prev = 0.0f;
213 val_prev = ((index_prev < 0 || index_prev >= mSignalCache.size() || mSignalCache[index_prev] < 0.0f)
215 : mSignalCache[index_prev]);
216 float const cfd_val = 5.0f * val_prev -
val;
217 if (std::abs(
val) >
params.mCFD_trsh && !is_positive && cfd_val > 0.0f) {
225 is_positive = cfd_val > 0.0f;
238 const float* tp =
times.data();
239 size_t const m =
times.size() / Vc::float_v::Size;
240 for (
size_t i = 0;
i <
m; ++
i) {
242 tp += Vc::float_v::Size;
243 Vc::prefetchForOneRead(tp);
248 for (
size_t i = Vc::float_v::Size *
m;
i <
times.size(); ++
i, ++tp) {
255 std::vector<o2::ft0::Digit>& digitsBC,
256 std::vector<o2::ft0::ChannelData>& digitsCh,
257 std::vector<o2::ft0::DetTrigInput>& digitsTrig,
262 LOG(
debug) <<
" process firstBCinDeque " << firstBCinDeque <<
" mIntRecord " << mIntRecord;
263 if (firstBCinDeque != mIntRecord) {
268 for (
auto const& hit : *hits) {
269 if (hit.GetEnergyLoss() > 0) {
273 Int_t hit_ch = hit.GetDetectorID();
276 if (mDeadChannelMap && !mDeadChannelMap->
isChannelAlive(hit_ch)) {
282 Bool_t is_A_side = (hit_ch < 4 * mGeometry.
NCellsA);
285 const Float_t timeOfFlight = hit.GetPos().R() / o2::constants::physics::LightSpeedCm2NS;
287 Double_t hit_time = hit.GetTime() - timeOfFlight + timeOffset + mIntRecord.getTimeOffsetWrtBC();
289 if (hit_time > 150) {
294 if (mCache.size() <= relBC.bc) {
295 mCache.resize(relBC.bc + 1);
297 mCache[relBC.bc].hits.emplace_back(
BCCache::particle{hit_ch, hit_time - relBC.bc2ns()});
299 Int_t parentID = hit.GetTrackID();
300 if (parentID != parent) {
301 mCache[relBC.bc].labels.emplace(parentID, mEventID, mSrcID, hit_ch);
307void Digitizer::storeBC(BCCache&
bc,
308 std::vector<o2::ft0::Digit>& digitsBC,
309 std::vector<o2::ft0::ChannelData>& digitsCh,
310 std::vector<o2::ft0::DetTrigInput>& digitsTrig,
313 if (
bc.hits.empty()) {
318 int n_hit_A = 0, n_hit_C = 0, mean_time_A = 0, mean_time_C = 0;
319 int summ_ampl_A = 0, summ_ampl_C = 0;
320 int sum_A_ampl = 0, sum_C_ampl = 0;
322 std::vector<int> sum_ampl_ipmt(nPMTs, 0);
324 std::map<uint8_t, int> mapPMhash2sumAmpl;
325 for (
const auto&
entry : pmLookupTable.getPMs()) {
326 mapPMhash2sumAmpl.insert({
entry.first, 0});
331 int first = digitsCh.size(), nStored = 0;
332 auto& particles =
bc.hits;
333 std::sort(std::begin(particles), std::end(particles));
334 auto channel_end = particles.begin();
335 std::vector<float> channel_times;
337 for (Int_t ipmt = 0; ipmt <
params.mMCPs; ++ipmt) {
338 auto channel_begin = channel_end;
339 channel_end = std::find_if(channel_begin, particles.end(),
340 [ipmt](BCCache::particle
const& p) { return p.hit_ch != ipmt; });
344 if (channel_end - channel_begin <
params.mAmp_trsh) {
347 channel_times.resize(channel_end - channel_begin);
348 std::transform(channel_begin, channel_end, channel_times.begin(), [](BCCache::particle
const& p) { return p.hit_time; });
349 int chain = (std::rand() % 2) ? 1 : 0;
350 auto cfd =
get_time(channel_times, mDeadTimes[ipmt].intrec.bc2ns() -
351 firstBCinDeque.
bc2ns() +
352 mDeadTimes[ipmt].deadTime);
353 mDeadTimes[ipmt].intrec = firstBCinDeque;
354 mDeadTimes[ipmt].deadTime = cfd.deadTime;
364 int smeared_time = 1000. * (*cfd.particle -
params.mCfdShift) *
params.mChannelWidthInverse + miscalib;
365 bool is_time_in_signal_gate = (smeared_time > -
params.mTime_trg_gate && smeared_time <
params.mTime_trg_gate);
367 float amp = is_time_in_signal_gate ?
params.mMV_2_Nchannels *
charge : 0;
375 LOG(
debug) << mEventID <<
" bc " << firstBCinDeque.
bc <<
" orbit " << firstBCinDeque.
orbit <<
", ipmt " << ipmt <<
", smeared_time " << smeared_time <<
" nStored " << nStored <<
" offset " << miscalib;
376 if (is_time_in_signal_gate) {
384 digitsCh.emplace_back(ipmt, smeared_time,
int(amp),
chain);
389 Bool_t is_A_side = (ipmt < 4 * mGeometry.
NCellsA);
390 if (!is_time_in_signal_gate) {
394 sum_ampl_ipmt[ipmt] += amp;
399 mean_time_A += smeared_time;
403 mean_time_C += smeared_time;
407 for (
size_t i = 0;
i < sum_ampl_ipmt.size();
i++) {
408 sum_ampl_ipmt[
i] = sum_ampl_ipmt[
i] >> 3;
410 sum_A_ampl += sum_ampl_ipmt[
i];
412 sum_C_ampl += sum_ampl_ipmt[
i];
417 int sum_PM_ampl_debug = 0;
418 int sum_PM_ampl_A_debug = 0;
419 int sum_PM_ampl_C_debug = 0;
420 for (
const auto&
entry : mapPMhash2sumAmpl) {
421 int pmAmpl = (
entry.second >> 3);
422 sum_PM_ampl_debug += pmAmpl;
423 if (pmLookupTable.isASide(
entry.first)) {
424 sum_PM_ampl_A_debug += pmAmpl;
426 sum_PM_ampl_C_debug += pmAmpl;
429 LOG(
debug) <<
"Sum PM amplitude (LUT-based): total=" << sum_PM_ampl_debug
430 <<
" A-side=" << sum_PM_ampl_A_debug
431 <<
" C-side=" << sum_PM_ampl_C_debug;
433 Bool_t is_A, is_C, isVertex, is_Central, is_SemiCentral = 0;
436 is_Central = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 *
params.mtrg_central_trh;
437 is_SemiCentral = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 *
params.mtrg_semicentral_trh && !is_Central;
440 int timeA = is_A ? mean_time_A / n_hit_A : -5000;
441 int timeC = is_C ? mean_time_C / n_hit_C : -5000;
443 isVertex = is_A && is_C && (vertex_time > -
params.mTime_trg_vertex_gate && vertex_time <
params.mTime_trg_vertex_gate);
444 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;
446 bool isLaser =
false;
447 bool isOutputsAreBlocked =
false;
448 bool isDataValid =
true;
452 digitsBC.emplace_back(
first, nStored, firstBCinDeque, triggers, mEventID - 1);
453 digitsTrig.emplace_back(firstBCinDeque, is_A, is_C, isVertex, is_Central, is_SemiCentral);
454 size_t const nBC = digitsBC.size();
456 labels.addElement(nBC - 1, lbl);
461 LOG(
debug) <<
"Event ID: " << mEventID <<
", bc " << firstBCinDeque.
bc <<
", N hit " <<
bc.hits.size();
464 <<
" mean time C: " << triggers.
getTimeC() <<
" nStored " << nStored;
471 std::vector<o2::ft0::ChannelData>& digitsCh,
472 std::vector<o2::ft0::DetTrigInput>& digitsTrig,
476 assert(firstBCinDeque <= mIntRecord);
478 while (firstBCinDeque < mIntRecord && !mCache.empty()) {
479 storeBC(mCache.front(), digitsBC, digitsCh, digitsTrig,
labels);
483 firstBCinDeque = mIntRecord;
487 std::vector<o2::ft0::ChannelData>& digitsCh,
488 std::vector<o2::ft0::DetTrigInput>& digitsTrig,
492 assert(firstBCinDeque <= mIntRecord);
494 while (!mCache.empty()) {
495 storeBC(mCache.front(), digitsBC, digitsCh, digitsTrig,
labels);
503 auto const sinc = [](
double x) {
x *= TMath::Pi();
return (std::abs(
x) < 1
e-12) ? 1.0 : std::sin(
x) /
x; };
507 mNumNoiseSamples = std::ceil(
params.mBunchWidth /
params.mNoisePeriod);
508 mNoiseSamples.resize(mNumNoiseSamples);
511 for (
size_t i = 0,
n = mSincTable.size();
i <
n; ++
i) {
515 mSincTable[
i].resize(3 * mNumNoiseSamples);
516 for (
int j = -mNumNoiseSamples;
j < 2 * mNumNoiseSamples; ++
j) {
517 mSincTable[
i][mNumNoiseSamples +
j] =
params.mNoiseVar * sinc((
time + 0.5f *
params.mBunchWidth) /
params.mNoisePeriod -
j);
532 LOG(info) <<
" @@@ Digitizer::init " << std::endl;
546 LOG(info) <<
" Run Digitzation with parametrs: \n"
547 <<
" CFD amplitude threshold \n " <<
params.mCFD_trsh <<
" CFD signal gate in ps \n"
548 <<
params.mTime_trg_gate <<
"shift to have signal around zero after CFD trancformation \n"
549 <<
params.mCfdShift <<
"CFD distance between 0.3 of max amplitude to max \n"
550 <<
params.mCFDShiftPos <<
"MIP -> mV " <<
params.mMip_in_V <<
" Pe in MIP \n"
551 <<
params.mPe_in_mip <<
"noise level " <<
params.mNoiseVar <<
" noise frequency \n"
std::vector< unsigned long > times
#define O2ParamImpl(classname)
uint64_t exp(uint64_t base, uint8_t exp) noexcept
General constants in FT0.
Configurable digitization parameters.
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char e
configStandalone Reconstruction rec_tpc float
constexpr int p1()
constexpr to accelerate the coordinates changing
std::vector< o2::MCCompLabel > labels
Definition of a container to keep Monte Carlo truth external to simulation objects.
Header to collect physics constants.
static const FT0DigParam & Instance()
uint8_t getNChanC() const
void setTriggers(uint8_t trgsig, uint8_t chanA, uint8_t chanC, int32_t aamplA, int32_t aamplC, int16_t atimeA, int16_t atimeC)
uint8_t getNChanA() const
void 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 > &label)
void printParameters() const
CFDOutput get_time(const std::vector< float > ×, float deadTime)
VcType signalFormVc(VcType x) const
void 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)
void 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 > &label)
double measure_amplitude(const std::vector< float > ×) const
float signalForm(float x) const
static constexpr int NCellsA
static constexpr int NCellsC
static const PMLookupTable & Instance()
Return the process-wide immutable lookup table.
GLuint GLfloat GLfloat GLfloat x1
GLboolean GLboolean GLboolean b
GLenum const GLfloat * params
GLuint GLsizei const GLchar * label
GLboolean GLboolean GLboolean GLboolean a
uint8_t int statusCode int
Float signalForm_i(Float x)
Vc::float_v signalForm_integralVc(Vc::float_v x)
float signalForm_integral(float x)
uint16_t bc
bunch crossing ID of interaction
static double bc2ns(int bc, unsigned int orbit)
const bool isChannelAlive(const uint8_t &chId) const
static constexpr std::size_t sNCHANNELS_PM
static constexpr double SIGNAL_CACHE_DT
static constexpr int SIGNAL_TABLE_SIZE
std::array< int16_t, o2::ft0::Geometry::Nchannels > mTimeOffsets
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"