Do full waveform time digitization.
As a result, the digitized hits are appended to TOPRawDigits, then they are converted to TOPDigits and the relations to TOPSimHits and MCParticles are set with proper weights.
257 {
258
259
260
261 const auto* geo = TOPGeometryPar::Instance()->getGeometry();
262 const auto& tdc = geo->getNominalTDC();
263
264
265
266 std::vector<unsigned short> windowNumbers;
267 int offsetWindows = s_offsetWindows + m_windowShift;
268 int winnum = s_window + offsetWindows;
269 if (winnum < 0) winnum += s_storageDepth;
270 windowNumbers.push_back(winnum % s_storageDepth);
271 for (unsigned i = 1; i < s_readoutWindows; i++) {
272 windowNumbers.push_back((windowNumbers.back() + 1) % s_storageDepth);
273 }
274 std::vector<double> baselines(windowNumbers.size(), 0);
275 std::vector<double> rmsNoises(windowNumbers.size(), m_rmsNoise);
276 double averagePedestal = tdc.getAveragePedestal();
277 std::vector<double> pedestals(windowNumbers.size(), averagePedestal);
278 int adcRange = tdc.getADCRange();
279
280 int startSample = -offsetWindows * TOPRawWaveform::c_WindowSize;
281 auto wfData = generateWaveform(startSample, baselines, rmsNoises, pedestals,
282 adcRange);
283
284
285
286 auto* waveform = waveforms.appendNew(m_moduleID, m_pixelID, m_channel, m_scrodID,
287 windowNumbers[0], 0, wfData);
288 waveform->setStorageWindows(windowNumbers);
289 waveform->setPedestalSubtractedFlag(true);
290
291
292
293 waveform->featureExtraction(threshold, hysteresis, thresholdCount);
294
295
296
297 for (const auto& feature : waveform->getFeatureExtractionData()) {
298 int sampleRise = feature.sampleRise;
299 if (sampleRise < 0) continue;
300 unsigned iwin = sampleRise / TOPRawWaveform::c_WindowSize;
301 if (iwin >= windowNumbers.size()) continue;
302 sampleRise -= iwin * TOPRawWaveform::c_WindowSize;
303 if (s_maskSamples and (sampleRise % 64) > 55) continue;
304
305 unsigned DsamplePeak = feature.samplePeak - feature.sampleRise;
306 if ((DsamplePeak & 0x0F) != DsamplePeak) continue;
307
308 unsigned DsampleFall = feature.sampleFall - feature.sampleRise;
309 if ((DsampleFall & 0x0F) != DsampleFall) continue;
310
311 std::vector<unsigned short> winNumbers;
312 for (unsigned i = iwin; i < windowNumbers.size(); i++) {
313 winNumbers.push_back(windowNumbers[i]);
314 }
315
316
317 auto* rawDigit = rawDigits.appendNew(m_scrodID, TOPRawDigit::c_ProductionDebug);
318 rawDigit->setCarrierNumber(m_carrier);
319 rawDigit->setASICNumber(m_asic);
320 rawDigit->setASICChannel(m_chan);
321 rawDigit->setASICWindow(windowNumbers[iwin]);
322 rawDigit->setStorageWindows(winNumbers);
323 rawDigit->setSampleRise(sampleRise);
324 rawDigit->setDeltaSamplePeak(DsamplePeak);
325 rawDigit->setDeltaSampleFall(DsampleFall);
326 rawDigit->setValueRise0(feature.vRise0);
327 rawDigit->setValueRise1(feature.vRise1);
328 rawDigit->setValueFall0(feature.vFall0);
329 rawDigit->setValueFall1(feature.vFall1);
330 rawDigit->setValuePeak(feature.vPeak);
331 rawDigit->setIntegral(feature.integral);
332 double rawTimeLeading = rawDigit->getCFDLeadingTime();
333 int sampleDivisions = 0x1 << tdc.getSubBits();
334 int tfine = int(rawTimeLeading * sampleDivisions) % sampleDivisions;
335 if (tfine < 0) tfine += sampleDivisions;
336 rawDigit->setTFine(tfine);
337 rawDigit->addRelationTo(waveform);
338
339 double rawTimeFalling = rawDigit->getCFDFallingTime();
340 double rawTimeErr = rawDigit->getCFDLeadingTimeError(m_rmsNoise);
341
342
343 int nwin = iwin + offsetWindows;
344 rawTimeLeading += nwin * TOPRawDigit::c_WindowSize;
345 rawTimeFalling += nwin * TOPRawDigit::c_WindowSize;
346
347 double cfdTime = m_sampleTimes->getTime(s_window, rawTimeLeading) - m_timeOffset;
348 double width = m_sampleTimes->getDeltaTime(s_window,
349 rawTimeFalling,
350 rawTimeLeading);
351 int sample = static_cast<int>(rawTimeLeading);
352 if (rawTimeLeading < 0) sample--;
353 double timeError = rawTimeErr * m_sampleTimes->getTimeBin(s_window, sample);
354 double timeErrorSq = timeError * timeError + m_calErrorsSq;
355 int pulseHeight = rawDigit->getValuePeak();
356
357 if (s_timeWalk) {
358 cfdTime -= (*s_timeWalk)->getTimeWalk(pulseHeight);
359 timeErrorSq += (*s_timeWalk)->getSigmaSq(pulseHeight);
360 }
361
362
363 auto* digit = digits.appendNew(m_moduleID, m_pixelID, rawTimeLeading);
364 digit->setTime(cfdTime);
365 digit->setTimeError(
sqrt(timeErrorSq));
366 digit->setPulseHeight(pulseHeight);
367 digit->setIntegral(rawDigit->getIntegral());
368 digit->setPulseWidth(width);
369 digit->setChannel(m_channel);
370 digit->setFirstWindow(s_window);
371 if (m_sampleTimes->isCalibrated()) {
372 digit->addStatus(TOPDigit::c_TimeBaseCalibrated);
373 }
374 digit->addRelationTo(rawDigit);
375
376
377
378 if (!rawDigit->isFEValid() or rawDigit->isPedestalJump()) {
379 digit->setHitQuality(TOPDigit::c_Junk);
380 }
381
382
383
384 std::multimap<double, const Hit*, std::greater<double>> weights;
385 for (const auto& hit : m_times) {
386 double hitTime = hit.first;
387 double weight = 0;
388 const auto* pulseShape = hit.second.shape;
389 if (pulseShape) {
390 weight = fabs(pulseShape->getValue(cfdTime - hitTime));
391 }
392 if (weight > 0.01) {
393 weight *= hit.second.pulseHeight;
394 weights.insert(std::pair<double, const Hit*>(weight, &hit.second));
395 }
396 }
397 double sum = 0;
398 for (const auto& w : weights) sum += w.first;
399 if (sum == 0) continue;
400
401 std::map<MCParticle*, double> relatedMCParticles;
402 for (const auto& w : weights) {
403 auto weight = w.first / sum;
404 const auto* simHit = w.second->simHit;
405 if (simHit and weight > 0) {
406 digit->addRelationTo(simHit, weight);
407 RelationVector<MCParticle> particles = simHit->getRelationsFrom<MCParticle>();
408 for (unsigned i = 0; i < particles.size(); ++i) {
409 relatedMCParticles[particles[i]] += particles.weight(i) * weight;
410 }
411 } else if (w.second->type == c_CalPulse and weight > 0.90) {
412 digit->setHitQuality(TOPDigit::c_CalPulse);
413 }
414 }
415 for (const auto& x : relatedMCParticles) digit->addRelationTo(x.first, x.second);
416 }
417
418 }
1.13.2