11 #include <top/modules/TOPDigitizer/TimeDigitizer.h>
12 #include <top/geometry/TOPGeometryPar.h>
13 #include <mdst/dataobjects/MCParticle.h>
14 #include <framework/logging/Logger.h>
26 unsigned TimeDigitizer::s_storageDepth = 0;
27 unsigned TimeDigitizer::s_readoutWindows = 0;
28 int TimeDigitizer::s_offsetWindows = 0;
29 unsigned TimeDigitizer::s_window = 0;
30 bool TimeDigitizer::s_maskSamples =
false;
31 DBObjPtr<TOPCalTimeWalk>* TimeDigitizer::s_timeWalk = 0;
33 TimeDigitizer::TimeDigitizer(
int moduleID,
int pixelID,
double timeOffset,
34 double calErrorsSq,
int shift,
double rmsNoise,
36 m_moduleID(moduleID), m_pixelID(pixelID), m_timeOffset(timeOffset),
37 m_calErrorsSq(calErrorsSq), m_rmsNoise(rmsNoise), m_sampleTimes(&sampleTimes),
41 m_channel = channelMapper.getChannel(pixelID);
42 if (!channelMapper.isChannelValid(
m_channel)) {
43 B2ERROR(
"TimeDigitizer::TimeDigitizer: invalid channel");
53 B2ERROR(
"TimeDigitizer::TimeDigitizer: no valid frontend map found");
66 hit.pulseHeight = pulseHeight;
70 case static_cast<int>(c_Hit):
73 case static_cast<int>(c_ChargeShare):
76 case static_cast<int>(c_CrossTalk):
78 B2ERROR(
"TOP::TimeDigitizer: waveform shape of cross-talk not yet available");
80 case static_cast<int>(c_CalPulse):
86 m_times.insert(std::pair<double, const Hit>(t, hit));
93 typedef std::multimap<double, const TimeDigitizer::Hit>::const_iterator
Iterator;
99 double timeJitter)
const
108 double pileupTime = tdc.getPileupTime();
111 std::vector<Iterator> ranges;
112 ranges.push_back(
m_times.begin());
113 double prevTime =
m_times.begin()->first;
115 if (it->first - prevTime > hitResolveTime) ranges.push_back(it);
116 prevTime = it->first;
118 ranges.push_back(
m_times.end());
121 for (
unsigned k = 0; k < ranges.size() - 1; k++) {
124 std::vector<double> times;
125 std::vector<double> weights;
126 std::vector<const TOPSimHit*> simHits;
129 prevTime = ranges[k]->first;
130 for (Iterator it = ranges[k]; it != ranges[k + 1]; ++it) {
131 if (it->first - prevTime > pileupTime)
break;
132 times.push_back(it->first);
133 double pulseHeight = gRandom->Gaus(it->second.pulseHeight,
m_rmsNoise);
134 weights.push_back(pulseHeight);
135 simHits.push_back(it->second.simHit);
136 prevTime = it->first;
141 for (
const auto& weight : weights) height += weight;
144 for (
auto& weight : weights) weight /= height;
147 double width = gRandom->Gaus(2.3, 0.4);
148 if (width < 0.5) width = 0.5;
152 for (
unsigned j = 0; j < times.size(); j++) time += times[j] * weights[j];
155 if (timeJitter > 0) time += gRandom->Gaus(0., timeJitter);
159 double mean = time + width / 2;
160 double sigma = width / 2.35482;
162 int samplePeak = tdc.getSample(mean);
163 if (mean - tdc.getSampleTime(samplePeak) > tdc.getSampleWidth() / 2) samplePeak++;
164 short vPeak = height *
gauss(tdc.getSampleTime(samplePeak), mean, sigma);
165 if (vPeak < threshold)
continue;
166 if (vPeak > 3000) vPeak = 3000;
168 double halfWid = sigma * sqrt(-2 * log(0.5 * vPeak / height));
169 int sampleRise = tdc.getSample(mean - halfWid);
170 if (!tdc.isSampleValid(sampleRise))
continue;
171 short vRise0 = height *
gauss(tdc.getSampleTime(sampleRise), mean, sigma);
172 short vRise1 = height *
gauss(tdc.getSampleTime(sampleRise + 1), mean, sigma);
174 int sampleFall = tdc.getSample(mean + halfWid);
175 if (!tdc.isSampleValid(sampleFall + 1))
continue;
176 short vFall0 = height *
gauss(tdc.getSampleTime(sampleFall), mean, sigma);
177 short vFall1 = height *
gauss(tdc.getSampleTime(sampleFall + 1), mean, sigma);
180 double halfDt = sigma * sqrt(-2 * log(threshold / height));
181 int overThr = tdc.getSample(mean + halfDt) - tdc.getSample(mean - halfDt);
182 if (overThr < thresholdCount)
continue;
186 int numSamples = (sampleFall - sampleRise) * 4;
187 double sigmaIntegral = 0;
188 if (numSamples > 1) sigmaIntegral =
m_rmsNoise * sqrt(numSamples - 1);
189 double integral = gRandom->Gaus(height * 7.0, sigmaIntegral);
194 rawDigit->setASICNumber(
m_asic);
195 rawDigit->setASICChannel(
m_chan);
197 rawDigit->setSampleRise(sampleRise);
198 rawDigit->setDeltaSamplePeak(samplePeak - sampleRise);
199 rawDigit->setDeltaSampleFall(sampleFall - sampleRise);
200 rawDigit->setValueRise0(vRise0);
201 rawDigit->setValueRise1(vRise1);
202 rawDigit->setValueFall0(vFall0);
203 rawDigit->setValueFall1(vFall1);
204 rawDigit->setValuePeak(vPeak);
205 rawDigit->setIntegral(integral);
207 double rawTime = rawDigit->getCFDLeadingTime();
208 double cfdTime = rawTime * tdc.getSampleWidth() - tdc.getOffset();
209 double cfdWidth = rawDigit->getFWHM() * tdc.getSampleWidth();
210 int sampleDivisions = 0x1 << tdc.getSubBits();
211 int tfine = int(rawTime * sampleDivisions) % sampleDivisions;
212 if (tfine < 0) tfine += sampleDivisions;
213 rawDigit->setTFine(tfine);
218 digit->setTime(cfdTime);
219 digit->setTimeError(timeJitter);
220 digit->setPulseHeight(rawDigit->getValuePeak());
221 digit->setIntegral(rawDigit->getIntegral());
222 digit->setPulseWidth(cfdWidth);
224 digit->setFirstWindow(rawDigit->getASICWindow());
225 digit->setStatus(TOPDigit::c_OffsetSubtracted);
226 digit->addRelationTo(rawDigit);
230 for (
unsigned j = 0; j < simHits.size(); j++) {
231 const auto* simHit = simHits[j];
232 const double& weight = weights[j];
234 digit->addRelationTo(simHit, weight);
236 for (
unsigned i = 0; i < particles.size(); ++i) {
237 digit->
addRelationTo(particles[i], particles.weight(i) * weight);
253 int thresholdCount)
const
263 std::vector<unsigned short> windowNumbers;
265 int winnum =
s_window + offsetWindows;
267 windowNumbers.push_back(winnum);
269 windowNumbers.push_back((windowNumbers.back() + 1) %
s_storageDepth);
271 std::vector<double> baselines(windowNumbers.size(), 0);
272 std::vector<double> rmsNoises(windowNumbers.size(),
m_rmsNoise);
273 double averagePedestal = tdc.getAveragePedestal();
274 std::vector<double> pedestals(windowNumbers.size(), averagePedestal);
275 int adcRange = tdc.getADCRange();
278 auto wfData =
generateWaveform(startSample, baselines, rmsNoises, pedestals,
284 windowNumbers[0], 0, wfData);
285 waveform->setStorageWindows(windowNumbers);
286 waveform->setPedestalSubtractedFlag(
true);
290 waveform->featureExtraction(threshold, hysteresis, thresholdCount);
294 for (
const auto& feature : waveform->getFeatureExtractionData()) {
295 int sampleRise = feature.sampleRise;
296 if (sampleRise < 0)
continue;
298 if (iwin >= windowNumbers.size())
continue;
302 unsigned DsamplePeak = feature.samplePeak - feature.sampleRise;
303 if ((DsamplePeak & 0x0F) != DsamplePeak)
continue;
305 unsigned DsampleFall = feature.sampleFall - feature.sampleRise;
306 if ((DsampleFall & 0x0F) != DsampleFall)
continue;
308 std::vector<unsigned short> winNumbers;
309 for (
unsigned i = iwin; i < windowNumbers.size(); i++) {
310 winNumbers.push_back(windowNumbers[i]);
316 rawDigit->setASICNumber(
m_asic);
317 rawDigit->setASICChannel(
m_chan);
318 rawDigit->setASICWindow(windowNumbers[iwin]);
319 rawDigit->setStorageWindows(winNumbers);
320 rawDigit->setSampleRise(sampleRise);
321 rawDigit->setDeltaSamplePeak(DsamplePeak);
322 rawDigit->setDeltaSampleFall(DsampleFall);
323 rawDigit->setValueRise0(feature.vRise0);
324 rawDigit->setValueRise1(feature.vRise1);
325 rawDigit->setValueFall0(feature.vFall0);
326 rawDigit->setValueFall1(feature.vFall1);
327 rawDigit->setValuePeak(feature.vPeak);
328 rawDigit->setIntegral(feature.integral);
329 double rawTimeLeading = rawDigit->getCFDLeadingTime();
330 int sampleDivisions = 0x1 << tdc.getSubBits();
331 int tfine = int(rawTimeLeading * sampleDivisions) % sampleDivisions;
332 if (tfine < 0) tfine += sampleDivisions;
333 rawDigit->setTFine(tfine);
334 rawDigit->addRelationTo(waveform);
336 double rawTimeFalling = rawDigit->getCFDFallingTime();
337 double rawTimeErr = rawDigit->getCFDLeadingTimeError(
m_rmsNoise);
340 int nwin = iwin + offsetWindows;
348 int sample =
static_cast<int>(rawTimeLeading);
349 if (rawTimeLeading < 0) sample--;
352 int pulseHeight = rawDigit->getValuePeak();
355 cfdTime -= (*s_timeWalk)->getTimeWalk(pulseHeight);
356 timeErrorSq += (*s_timeWalk)->getSigmaSq(pulseHeight);
361 digit->setTime(cfdTime);
362 digit->setTimeError(sqrt(timeErrorSq));
363 digit->setPulseHeight(pulseHeight);
364 digit->setIntegral(rawDigit->getIntegral());
365 digit->setPulseWidth(width);
369 digit->addStatus(TOPDigit::c_TimeBaseCalibrated);
371 digit->addRelationTo(rawDigit);
375 if (!rawDigit->isFEValid() or rawDigit->isPedestalJump()) {
376 digit->setHitQuality(TOPDigit::c_Junk);
381 std::multimap<double, const Hit*, std::greater<double>> weights;
382 for (
const auto& hit :
m_times) {
383 double hitTime = hit.first;
385 const auto* pulseShape = hit.second.shape;
387 weight = fabs(pulseShape->getValue(cfdTime - hitTime));
390 weight *= hit.second.pulseHeight;
391 weights.insert(std::pair<double, const Hit*>(weight, &hit.second));
395 for (
const auto& w : weights) sum += w.first;
396 if (sum == 0)
continue;
397 for (
const auto& w : weights) {
398 auto weight = w.first / sum;
399 const auto* simHit = w.second->simHit;
400 if (simHit and weight > 0) {
401 digit->addRelationTo(simHit, weight);
403 for (
unsigned i = 0; i < particles.size(); ++i) {
404 digit->
addRelationTo(particles[i], particles.weight(i) * weight);
406 }
else if (w.second->type == c_CalPulse and weight > 0.90) {
407 digit->setHitQuality(TOPDigit::c_CalPulse);
416 const vector<double>& baselines,
417 const vector<double>& rmsNoises,
418 const vector<double>& pedestals,
422 if (baselines.empty() or baselines.size() != rmsNoises.size() or
423 baselines.size() != pedestals.size())
424 B2FATAL(
"TOP::TimeDigitizer: inconsistent vector sizes");
430 const auto& signalShape = geo->getSignalShape();
432 double dt = tdc.getSampleWidth();
433 double fmax = 1 / (2 * dt);
434 double tau1 = 1 / (2 * M_PI * signalShape.getPole1());
435 double tau2 = 1 / (2 * M_PI * signalShape.getPole2());
436 double bw1 = 1 / (4 * tau1);
437 double bw2 = 1 / (4 * (tau1 + tau2));
441 vector<double> waveform;
442 for (
auto rms : rmsNoises) {
443 double rms0 = rms * sqrt(fmax / bw2);
445 waveform.push_back(gRandom->Gaus(0, rms0));
451 double a1 = exp(-dt / tau1);
452 double v0 = gRandom->Gaus(0, rmsNoises[0] * sqrt(bw1 / bw2));
453 double prevValue = v0;
454 for (
auto& value : waveform) {
455 value = a1 * prevValue + (1 - a1) * value;
458 double a2 = exp(-dt / tau2);
460 for (
auto& value : waveform) {
461 value = a2 * prevValue + (1 - a2) * value;
468 for (
auto baseline : baselines) {
470 waveform[k] += baseline;
477 for (
const auto& hit :
m_times) {
478 double hitTime = hit.first;
479 double pulseHeight = hit.second.pulseHeight;
480 const auto* pulseShape = hit.second.shape;
481 if (not pulseShape)
continue;
484 double t = (*s_timeWalk)->getTimeWalk(pulseHeight);
485 double sig = (*s_timeWalk)->getSigma(pulseHeight);
486 timeWalk = gRandom->Gaus(t, sig);
488 int size = waveform.size();
489 for (
int sample = 0; sample < size; sample++) {
491 waveform[sample] += pulseHeight * pulseShape->getValue(t - timeWalk - hitTime);
500 for (
auto pedestal : pedestals) {
502 int adc = waveform[k] + pedestal;
503 if (adc < 0) adc = 0;
504 if (adc > adcRange) adc = adcRange;