algorithms/digi/PulseNoise.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Simon Gardner, Minho Kim
0003 //
0004 // Adds noise to a time series pulse
0005 //
0006
0007 #include <DDDigi/noise/FalphaNoise.h>
0008 #include <edm4hep/MCParticle.h>
0009 #include <edm4hep/SimCalorimeterHit.h>
0010 #include <edm4hep/SimTrackerHit.h>
0011 #include <podio/RelationRange.h>
0012 #include <cstddef>
0013 #include <gsl/pointers>
0014 #include <random>
0015 #include <vector>
0016
0017 #include "PulseNoise.h"
0018
0019 namespace eicrecon {
0020
0021 void PulseNoise::init() {}
0022
0023 void PulseNoise::process(const PulseNoise::Input& input, const PulseNoise::Output& output) const {
0024   const auto [headers, inPulses] = input;
0025   auto [outPulses]               = output;
0026
0027   // local random generator
0028   auto seed = m_uid.getUniqueID(*headers, name());
0029   std::default_random_engine generator(seed);
0030   dd4hep::detail::FalphaNoise falpha(m_cfg.poles, m_cfg.alpha, m_cfg.variance);
0031
0032   for (const auto& pulse : *inPulses) {
0033
0034     //Clone input pulse to a mutable output pulse
0035     auto out_pulse = outPulses->create();
0036     out_pulse.setCellID(pulse.getCellID());
0037     out_pulse.setInterval(pulse.getInterval());
0038     out_pulse.setTime(pulse.getTime());
0039
0040     float integral = 0;
0041     //Add noise to the pulse
0042     for (std::size_t i = 0; i < pulse.getAmplitude().size(); i++) {
0043       double noise     = falpha(generator) * m_cfg.scale + m_cfg.pedestal;
0044       double amplitude = pulse.getAmplitude()[i] + noise;
0045       out_pulse.addToAmplitude(amplitude);
0046       integral += amplitude;
0047     }
0048
0049 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0050     out_pulse.setIntegral(integral);
0051     out_pulse.setPosition(pulse.getPosition());
0052     out_pulse.addToPulses(pulse);
0053
0054     for (auto particle : pulse.getParticles()) {
0055       out_pulse.addToParticles(particle);
0056     }
0057     // Not sure if we want/need to keep the hits themselves at this point?
0058     for (auto hit : pulse.getTrackerHits()) {
0059       out_pulse.addToTrackerHits(hit);
0060     }
0061     for (auto hit : pulse.getCalorimeterHits()) {
0062       out_pulse.addToCalorimeterHits(hit);
0063     }
0064
0065 #endif
0066   }
0067
0068 } // PulseNoise:process
0069 } // namespace eicrecon