icarus::opdet::PMTsimulationAlg Class Reference

Algorithm class for the full simulation of PMT channels. More...

#include <PMTsimulationAlg.h>

Classes
struct	ConfigurationParameters_t
	Type holding all configuration parameters for this algorithm. More...
class	GainFluctuator
class	TimeToTickAndSubtickConverter
	Functor to convert tick point into a tick number and a subsample index. More...

Public Types
using	microseconds = util::quantities::microsecond
using	nanoseconds = util::quantities::nanosecond
using	hertz = util::quantities::hertz
using	megahertz = util::quantities::megahertz
using	picocoulomb = util::quantities::picocoulomb
using	tick = util::quantities::tick
using	ADCcount = DiscretePhotoelectronPulse::ADCcount
using	time_interval = detinfo::timescales::time_interval
using	optical_tick = detinfo::timescales::optical_tick

Public Member Functions
	PMTsimulationAlg (ConfigurationParameters_t const &config)
	Constructor. More...
std::tuple< std::vector < raw::OpDetWaveform > , std::optional < sim::SimPhotons > >	simulate (sim::SimPhotons const &photons, sim::SimPhotonsLite const &lite_photons)
	Returns the waveforms originating from simulated photons. More...
template<typename Stream >
void	printConfiguration (Stream &&out, std::string indent="") const
	Prints the configuration into the specified output stream. More...

Private Types
using	OpDetWaveformMaker_t = icarus::opdet::OpDetWaveformMakerClass< ADCcount >
using	Waveform_t = OpDetWaveformMaker_t::WaveformData_t
	Type internally used for storing waveforms. More...
using	WaveformValue_t = ADCcount::value_t
	Numeric type in waveforms. More...
using	PulseSampling_t = DiscretePhotoelectronPulse::Subsample_t
	Type of sampled pulse shape: sequence of samples, one per tick. More...
using	NoiseAdderFunc_t = void(PMTsimulationAlg::*)(Waveform_t &) const
	Type of member function to add electronics noise. More...

Private Member Functions
auto	makeGainFluctuator () const
	Returns a configured gain fluctuator object. More...
Waveform_t	CreateFullWaveform (sim::SimPhotons const &photons, sim::SimPhotonsLite const &lite_photons, std::optional< sim::SimPhotons > &photons_used) const
	Creates raw::OpDetWaveform objects from simulated photoelectrons. More...
std::vector< raw::OpDetWaveform >	CreateFixedSizeOpDetWaveforms (raw::Channel_t opChannel, Waveform_t const &waveform) const
	Creates raw::OpDetWaveform objects from a waveform data. More...
template<typename Combine >
void	AddPulseShape (PulseSampling_t const &pulse, Waveform_t &wave, tick const time_bin, Combine combination) const
	Adds a pulse to a waveform, starting at a given tick. More...
void	AddPhotoelectrons (PulseSampling_t const &pulse, Waveform_t &wave, tick const time_bin, WaveformValue_t const n) const
	Adds a number of pulses to a waveform, starting at a given tick. More...
void	AddNoise (Waveform_t &wave) const
void	AddNoise_faster (Waveform_t &wave) const
	Same as AddNoise() but using an alternative generator. More...
void	AddDarkNoise (Waveform_t &wave) const
std::vector< optical_tick >	FindTriggers (Waveform_t const &wvfm) const
	Ticks in the specified waveform where some signal activity starts. More...
std::vector< optical_tick >	CreateBeamGateTriggers () const
	Generate periodic interest points regardless the actual activity. More...
bool	KicksPhotoelectron () const
	Returns a random response whether a photon generates a photoelectron. More...
std::pair< ADCcount, ADCcount >	saturationRange () const
	Returns the ADC range allowed for photoelectron saturation. More...
void	ApplySaturation (Waveform_t &waveform, std::pair< ADCcount, ADCcount > const &range) const
void	ApplySaturation (Waveform_t &waveform) const
	Applies the configured photoelectron saturation on the waveform. More...

Static Private Member Functions
static void	ClipWaveform (Waveform_t &waveform, ADCcount min, ADCcount max)
	Forces waveform ADC within the min to max range (max included). More...

Private Attributes
ConfigurationParameters_t	fParams
	Complete algorithm configuration. More...
double	fQE
	PMT quantum efficiency. More...
megahertz	fSampling
	Wave sampling frequency [MHz]. More...
std::size_t	fNsamples
	Samples per waveform. More...
DiscretePhotoelectronPulse	wsp
NoiseAdderFunc_t const	fNoiseAdder
	Single photon pulse (sampled). More...

Static Private Attributes
static util::FastAndPoorGauss< 32768U, float > const	fFastGauss

Detailed Description

Algorithm class for the full simulation of PMT channels.

The algorithm creates simulated PMT waveforms as read out by ICARUS, including the generation of trigger primitives. Contributions to the waveforms include:

• physical photons
• dark noise
• electronics noise

The algorithm processes an optical channel at a time, independently and uncorrelated from the other channels. For each channel, multiple waveforms may be generated according to the readout parameters.

Activity sources

Physical photons

Photons are read from sim::SimPhotons data objects, each one pertaining a single optical detector channel. Each photon on the channel is assumed to have successfully reached the external surface of the photocathode, with the wavelength shifter. Depending on the upstream simulation, and in particular on the photon visibility library settings, the photon might have also already passed the wavelength shifting and even triggered the conversion to a detectable photoelectron.

Quantum efficiency is simulated to determine if each photon converts into a photoelectron on the internal side of the photocathode. The target quantum efficiency is specified via the QE configuration parameter. It is assumed that some level of quantum efficiency has already been simulated upstream: more precisely, that the quantum efficiency already applied is in the amount returned by detinfo::LArProperties::ScintPreScale(). Therefore:

1. the quantum efficiency applied here is only the residual one to go from detinfo::LArProperties::ScintPreScale() to the value in QE
2. there is no implement here to increase quantum efficiency, i.e. QE must not exceed detinfo::LArProperties::ScintPreScale()
3. if the configuration specifies a target quantum efficiency QE larger than the one applied upstream (detinfo::LArProperties::ScintPreScale()), a warning message is printed, and no change to quantum efficiency is performed

Note that if the upstream code has not applied any quantum efficiency, the configuration should give a detinfo::LArProperties::ScintPreScale() of 1.0.

Note

If the photon visibility library already includes the probability of the photon converting to a photoelectron, the quantum efficiency check here should be skipped by setting the efficiency to 1.

For each converting photon, a photoelectron is added to the channel by placing a template waveform shape into the channel waveform.

The timestamp of each waveform is based on the same scale as the trigger time, as defined by detinfo::DetectorClocks::TriggerTime(). On that scale, the timestamp pins down the time of the first sample of the waveform. Note that this is typically earlier than when the actual signal starts. More precisely, the signal is defined to start at an interest point (see FindTriggers() for their definition), and the waveform starts (at tick #0) earlier than that by a fraction PreTrigFraction of the readout window size ReadoutWindowSize (both are configuration parameters of the algorithm), allowing for that amount of pre-trigger data.

The configuration parameter TriggerOffsetPMT describes how much earlier than the trigger time the optical readout has started. Note that if an interest point (see above) happens early after optical readout has started, there might be not enough data to fill the pre-trigger data. In such case, the interest point will just be located earlier than usual within the final waveform. This situation may be caused for example by asynchronous physics events like scintillation light from cosmic rays or radioactive decay of the detector materials, or from a fluctuation of the noise.

Photoelectrons

The response of the PMT to a single photoelectron is passed to the algorithm as a full blown function of type SinglePhotonResponseFunc_t. The function needs to be valid for the lifetime of the algorithm, since the algorithm refers to without owning it, and it is expected not to change during that time. See icarus::opdet::SimPMTIcarus

To account for gain fluctuations, that shape is considered to correspond to a nominal gain (PMTspecs.gain configuration parameter), which is then fluctuated to obtain the effective gain. This feature can be disabled by setting configuration parameter FluctuateGain to false. The approximation used here is that the fluctuation is entirely due to the first stage of multiplication. The gain on the first stage is described as a random variable with Poisson distribution around the mean gain. The gain on a single photoelectron at the first stage is, in fact, an integral number in each case. The time spread of the signal may be increased by the difference in time of the different branches of the multiplication avalanche. Therefore, increasing or decreasing the number of branches, as it is done by changing the gain of the first stage, the time evolution of the signal will also be likewise affected. At this time we do not take this aspect into account in the simulation. For the nominal settings of a Hamamatsu 5912 photomultiplier (gain 10 ^7^, high multiplication on the first stage) the gain at the first stage is around 20, causing a fluctuation of about 20%. If the multiplication were equally distributed across the stages, that fluctuation would be almost 45%.

The first stage gain is computed by icarus::opdet::PMTsimulationAlg::ConfigurationParameters_t::PMTspecs_t::multiplicationStageGain().

Dark noise

Dark noise, i.e. the noise originating by "spontaneous" emission of a photoelectron in the photocathode without any external stimulation, is simulated by randomly extracting the time such emission happens. Each emission causes a photoelectron template waveform to be added at the extracted time. The rate of dark noise emission is set by configuration with DarkNoiseRate parameter.

Electronics noise

Electronics noise is described by Gaussian fluctuations of a given standard deviation, controlled by the configuration parameter AmpNoise. No noise correlation is simulated neither in time nor in space.

Configuration

PMT specifications

PMT specifications are used to evaluate the variance of the gain. The details of the calculation are documented in icarus::opdet::PMTsimulationAlg::ConfigurationParameters_t::PMTspecs_t::multiplicationStageGain().

The available parameters include:

• gain (default: 1e7): the nominal gain of the PMT; this is just a reference value.
• voltageDistribution is a sequence of values, one for each stage of the multiplication chain of the PMT. Each number represents the relative size of the resistance that determines the fall of the potential on that stage. Only the stages that contribute to the gain need to be included. The absolute value of each element is inconsequential. For example, a 10-stage PMT with the first stage having twice the resistance of all the other would be represented by the setting [ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1 ].
• dynodeK (default: 0.75) represents the dependence of the gain of a stage on the potential drop: (with the gain for stage , the drop of potential of that stage and the parameter set by dynodeK.

Random number generators

Three independent random engines are currently used in the simulation:

• "main" random engine:
  • residual quantum efficiency;
  • gain fluctuations;
• "dark noise" engine: dark current noise only;
• "electronics noise" engine: electronics noise only.

Structure of the algorithm

This section needs completion.

The algorithm is serviceable immediately after construction. Construction relies on a custom parameter data structure.

An utility, PMTsimulationAlgMaker, splits the set up in two parts:

1. configuration, where the full set of parameters is learned;
2. set up, where service providers, random number engines and the external single photon response function are acquired.

This is supposed to make the creation of the filling of parameter structure and the creation of the algorithm easier. At that point, a single sim::SimPhotons can be processed (simulate()) at a time, or multiple at the same time (see the multithreading notes below).

The function used to describe the single particle response is customizable and it must in fact be specified by the caller, since there is no default form. The function must implement the PulseFunction_t interface.

Multithreading notes

The algorithm processes one channel at a time, and it does not depend on event-level information; therefore, the same algorithm object can be used to process many events in sequence. On the other hand, multithreading is impaired by the random number generation, in the sense that multithreading will break reproducibility if the random engine is not magically thread-resistant.

If the set up is event-dependent, then this object can't be used for multiple events at the same time. There is no global state, so at least different instances of the algorithm can be run at the same time. In fact, the creation of an algorithm is expected to take negligible time compared to its run time on a single event, and it is conceivable to create a new algorithm instance for each event.

Definition at line 341 of file PMTsimulationAlg.h.

Member Typedef Documentation

using icarus::opdet::PMTsimulationAlg::ADCcount = DiscretePhotoelectronPulse::ADCcount

Definition at line 350 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::hertz = util::quantities::hertz

Definition at line 346 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::megahertz = util::quantities::megahertz

Definition at line 347 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::microseconds = util::quantities::microsecond

Definition at line 344 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::nanoseconds = util::quantities::nanosecond

Definition at line 345 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::NoiseAdderFunc_t = void (PMTsimulationAlg::*)(Waveform_t&) const

private

Type of member function to add electronics noise.

Definition at line 531 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::OpDetWaveformMaker_t = icarus::opdet::OpDetWaveformMakerClass<ADCcount>

private

Definition at line 521 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::optical_tick = detinfo::timescales::optical_tick

Definition at line 353 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::picocoulomb = util::quantities::picocoulomb

Definition at line 348 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::PulseSampling_t = DiscretePhotoelectronPulse::Subsample_t

private

Type of sampled pulse shape: sequence of samples, one per tick.

Definition at line 528 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::tick = util::quantities::tick

Definition at line 349 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::time_interval = detinfo::timescales::time_interval

Definition at line 352 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::Waveform_t = OpDetWaveformMaker_t::WaveformData_t

private

Type internally used for storing waveforms.

Definition at line 524 of file PMTsimulationAlg.h.

using icarus::opdet::PMTsimulationAlg::WaveformValue_t = ADCcount::value_t

private

Numeric type in waveforms.

Definition at line 525 of file PMTsimulationAlg.h.

Constructor & Destructor Documentation

icarus::opdet::PMTsimulationAlg::PMTsimulationAlg

(

ConfigurationParameters_t const &

config

)

Constructor.

Definition at line 102 of file PMTsimulationAlg.cxx.

  : fParams(config)
  , fQE(fParams.QEbase / fParams.larProp->ScintPreScale())
  , fSampling(fParams.clockData->OpticalClock().Frequency())
  , fNsamples(fParams.readoutEnablePeriod * fSampling) // us * MHz cancels out
  , wsp( // NOTE: wsp amplitude already includes sign from polarity
    *(fParams.pulseFunction),
    fSampling,
    fParams.pulseSubsamples, // tick subsampling
    1.0e-4_ADCf // stop sampling when ADC counts are below this value
    )
  , fNoiseAdder(fParams.useFastElectronicsNoise
      ? &icarus::opdet::PMTsimulationAlg::AddNoise_faster
      : &icarus::opdet::PMTsimulationAlg::AddNoise
    )
{
  using namespace util::quantities::electronics_literals;

  //  mf::LogDebug("PMTsimulationAlg") << "Sampling = " << fSampling << std::endl;

  // shape of single pulse

  // Correction due to scalling factor applied during simulation if any
  mf::LogDebug("PMTsimulationAlg") << "PMT corrected efficiency = " << fQE;

  if (fQE >= 1.0001) {
    mf::LogWarning("PMTsimulationAlg")
      << "WARNING: Quantum efficiency set in the configuration (QE="
        << fParams.QEbase << ") seems to be too large!"
      "\nThe photon visibility library is assumed to already include a"
        " quantum efficiency of " << fParams.larProp->ScintPreScale() <<
        " (`ScintPreScale` setting of `LArProperties` service)"
        " and here we are requesting a higher one."
      "\nThis configuration is not supported by `PMTsimulationAlg`,"
        " and this simulation will effectively apply a quantum efficiency of "
        << fParams.larProp->ScintPreScale();
      ;
  }

  // check that the sampled waveform has a sufficiently large range, so that
  // tails are below 10^-3 ADC counts (absolute value);
  // if this test fails, it's better to reduce the threshold in wsp constructor
  // (for analytical pulses converging to 0) or have a longer sampling that does
  // converge to 0 (10^-3 ADC is quite low though).
  wsp.checkRange(1.0e-3_ADCf, "PMTsimulationAlg");

} // icarus::opdet::PMTsimulationAlg::PMTsimulationAlg()

Member Function Documentation

void icarus::opdet::PMTsimulationAlg::AddDarkNoise ( Waveform_t &  wave ) const

private

Definition at line 595 of file PMTsimulationAlg.cxx.

                                                                     {
  /*
   * We assume leakage current ("dark noise") is completely stochastic and
   * distributed uniformly in time with a fixed and known rate.
   *
   * In these condition, the time between two consecutive events follows an
   * exponential distribution with as decay constant the same rate.
   * We extact all the leakage events (including the first one) according to
   * that distribution.
   *
   * We follow the "standard" approach of subsampling of tick as for the
   * photoelectron.
   *
   */
  using namespace util::quantities::frequency_literals;

  if (fParams.darkNoiseRate <= 0.0_Hz) return; // no dark noise

  // CLHEP random objects do not understand quantities, so we use scalars;
  // we choose to work with nanosecond
  CLHEP::RandExponential random(*(fParams.darkNoiseRandomEngine),
    (1.0 / fParams.darkNoiseRate).convertInto<nanoseconds>().value());

  // time to stop at: full duration of the waveform
  nanoseconds const maxTime = static_cast<double>(wave.size()) / fSampling;

  // the time of first leakage event:
  nanoseconds darkNoiseTime { random.fire() };

  TimeToTickAndSubtickConverter const toTickAndSubtick(wsp.nSubsamples());

  auto gainFluctuation = makeGainFluctuator();

  MF_LOG_TRACE("PMTsimulationAlg")
    << "Adding dark noise (" << fParams.darkNoiseRate << ") up to " << maxTime;

  while (darkNoiseTime < maxTime) {

    auto const [ tick, subtick ] = toTickAndSubtick(darkNoiseTime * fSampling);

    double const n = gainFluctuation(1.0); // leakage is one photoelectron
    MF_LOG_TRACE("PMTsimulationAlg")
      << " * at " << darkNoiseTime << " (" << tick << ", subsample " << subtick
      << ") x" << n;

    AddPhotoelectrons
      (wsp.subsample(subtick), wave, tick, static_cast<WaveformValue_t>(n));

    // time of the next leakage event:
    darkNoiseTime += nanoseconds{ random.fire() };

  } // while

} // icarus::opdet::PMTsimulationAlg::AddDarkNoise()

void icarus::opdet::PMTsimulationAlg::AddNoise ( Waveform_t &  wave ) const

private

Definition at line 558 of file PMTsimulationAlg.cxx.

                                                                 {

  CLHEP::RandGaussQ random
    (*fParams.elecNoiseRandomEngine, 0.0, fParams.ampNoise.value());
  for(auto& sample: wave) {
    ADCcount const noise { static_cast<float>(random.fire()) }; // Gaussian noise
    sample += noise;
  } // for sample

} // PMTsimulationAlg::AddNoise()

void icarus::opdet::PMTsimulationAlg::AddNoise_faster ( Waveform_t &  wave ) const

private

Same as AddNoise() but using an alternative generator.

Definition at line 571 of file PMTsimulationAlg.cxx.

                                                                        {

  /*
    * Compared to AddNoise(), we use a somehow faster random generator;
    * to squeeze the CPU cycles, we avoid the CLHEP interface as much as
    * possible; the random number from the engine is immediately converted
    * to single precision, and the rest of the math happens in there as well.
    * No virtual interfaces nor indirection is involved within this function
    * (except for CLHEP random engine). We generate a normal variable _z_
    * (standard deviation 1, mean 0) and we just scale it to the desired
    * standard deviation, not bothering to add the mean offset of 0.
    * Note that unless the random engine is multi-thread safe, this function
    * won't gain anything from multi-threading.
    */
  auto& engine = *fParams.elecNoiseRandomEngine;

  for(auto& sample: wave) {
    sample += fParams.ampNoise * fFastGauss(engine.flat()); // Gaussian noise
  } // for sample

} // PMTsimulationAlg::AddNoise_faster()

void icarus::opdet::PMTsimulationAlg::AddPhotoelectrons ( PulseSampling_t const &  pulse, Waveform_t &  wave, tick const  time_bin, WaveformValue_t const  n  ) const

private

Adds a number of pulses to a waveform, starting at a given tick.

Parameters

pulse	the sampling to add, scaled, to the waveform
wave	the waveform the pulses will be added to
time_bin	the tick of the waveform where the pulses start being added
n	the number of pulses added (it may be fractional)

All the samples of pulse are scaled by the factor n and then added to the sampling waveform wave, starting from the time_bin sample of this waveform.

The pulse samples are a sequence of ADC counts describing the single photoelectron pulse shape. The waveform is also a sequence of samples representing a optical detector channel digitized waveform, starting at tick #0.

Definition at line 517 of file PMTsimulationAlg.cxx.

        {

  if (n == 0.0) return;

  if (n == 1.0) {
    // simple addition
    AddPulseShape(pulse, wave, time_bin, std::plus<ADCcount>());
  }
  else {
    // multiply each `pulse` sample by `n`:
    AddPulseShape(pulse, wave, time_bin, [n](auto a, auto b) { return a+n*b; });
  }

} // icarus::opdet::PMTsimulationAlg::AddPhotoelectrons()

template<typename Combine >

void icarus::opdet::PMTsimulationAlg::AddPulseShape ( PulseSampling_t const &  pulse, Waveform_t &  wave, tick const  time_bin, Combine  combination  ) const

private

Adds a pulse to a waveform, starting at a given tick.

Template Parameters

Combine

binary operation combining two ADC counts into one

Parameters

pulse	the sampling to add to the waveform
wave	the waveform the pulse will be added to
time_bin	the tick of the waveform where the pulse starts
combination	how to combine the pulse and the waveform

This is the internal implementation of AddPhotoelectrons().

The combination functor behaves as a binary function taking the existing wave sample and the sample from the pulse at the same time and returning their combination as a new sample value.

Definition at line 538 of file PMTsimulationAlg.cxx.

{
  std::size_t const min = time_bin.value();
  std::size_t const max = std::min(min + pulse.size(), fNsamples);
  if (min >= max) return;

  std::transform(
    wave.begin() + min, wave.begin() + max,
    pulse.begin(),
    wave.begin() + min,
    combination
    );

} // icarus::opdet::PMTsimulationAlg::AddPulseShape()

void icarus::opdet::PMTsimulationAlg::ApplySaturation ( Waveform_t &  waveform, std::pair< ADCcount, ADCcount > const &  range  ) const

private

Applies the configured photoelectron saturation on the waveform, only if the saturation is cutting into the digitisation range.

Definition at line 681 of file PMTsimulationAlg.cxx.

{
  //
  // simple sharp capping/cupping of ADC values
  //
  auto const boundaries = saturationRange();
  
  // saturation is out of boundaries: do not apply it.
  if ((boundaries.first <= range.first) && (boundaries.second >= range.second))
    return;
  
  ClipWaveform(waveform, boundaries.first, boundaries.second);
  
} // icarus::opdet::PMTsimulationAlg::ApplySaturation(range)

void icarus::opdet::PMTsimulationAlg::ApplySaturation ( Waveform_t &  waveform ) const

private

Applies the configured photoelectron saturation on the waveform.

Definition at line 668 of file PMTsimulationAlg.cxx.

{
  //
  // simple sharp capping/cupping of ADC values
  //
  auto const boundaries = saturationRange();
  ClipWaveform(waveform, boundaries.first, boundaries.second);
  
} // icarus::opdet::PMTsimulationAlg::ApplySaturation()

void icarus::opdet::PMTsimulationAlg::ClipWaveform ( Waveform_t &  waveform, ADCcount  min, ADCcount  max  )

staticprivate

Forces waveform ADC within the min to max range (max included).

Definition at line 700 of file PMTsimulationAlg.cxx.

{
  using ClamperFunc_t = std::function<ADCcount(ADCcount)>;
  ClamperFunc_t const clamper =
    (min == std::numeric_limits<ADCcount>::min())
      ? ((max == std::numeric_limits<ADCcount>::max())
        ? ClamperFunc_t{ [](ADCcount s){ return s; } }
        : ClamperFunc_t{ [max](ADCcount s){ return std::min(s, max); } }
        )
      : ((max == std::numeric_limits<ADCcount>::max())
        ? ClamperFunc_t{ [min](ADCcount s){ return std::max(s, min); } }
        : ClamperFunc_t{ [min,max](ADCcount s){ return std::clamp(s, min, max); } }
        )
      ;
  
  std::transform(waveform.cbegin(), waveform.cend(), waveform.begin(), clamper);
  
} // icarus::opdet::PMTsimulationAlg::ClipWaveform()

auto icarus::opdet::PMTsimulationAlg::CreateBeamGateTriggers ( ) const

private

Generate periodic interest points regardless the actual activity.

Returns

a collection of ticks where we pretend interesting activity to be

See Also

FindTriggers()

This methods emits a list of interest points according to the algorithm configuration. More precisely, if CreateBeamGateTriggers is configured true, BeamGateTriggerNReps interest points are generated at BeamGateTriggerRepPeriod intervals, starting from the beam gate time as defined by detinfo::DetectorClocks::BeamGateTime().

See FindTriggers() for the meaning of "interest point".

Note

It is assumed that tick 0 happens at a time defined by triggerOffsetPMT configuration parameter after the trigger (but since the value of that parameter is expected to be negative, tick 0 effectively happens before the trigger).

Definition at line 341 of file PMTsimulationAlg.cxx.

  {
    using namespace util::quantities::time_literals;
    using detinfo::timescales::trigger_time;

    detinfo::DetectorTimings const& timings
      = detinfo::makeDetectorTimings(fParams.clockData);

    std::vector<optical_tick> trigger_locations;
    trigger_locations.reserve(fParams.beamGateTriggerNReps);
    trigger_time const beamTime = timings.toTriggerTime(timings.BeamGateTime());

    for(unsigned int i_trig=0; i_trig<fParams.beamGateTriggerNReps; ++i_trig) {
      trigger_time const trig_time
        = beamTime
        + i_trig * fParams.beamGateTriggerRepPeriod
        - fParams.triggerOffsetPMT
        ;
      if (trig_time < 0_us) continue;
      if (trig_time > fParams.readoutEnablePeriod) break;
      trigger_locations.push_back
        (optical_tick::castFrom(trig_time.quantity()*fSampling));
    }
    return trigger_locations;
  }

std::vector< raw::OpDetWaveform > icarus::opdet::PMTsimulationAlg::CreateFixedSizeOpDetWaveforms ( raw::Channel_t  opChannel, Waveform_t const &  waveform  ) const

private

Creates raw::OpDetWaveform objects from a waveform data.

Parameters

opChannel	number of optical detector channel the data belongs to
waveform	the waveform data

Returns

a collection of raw::OpDetWaveform

The waveform data is a sequence of samples on the optical detector channel, starting at the beginning of the optical time clock, that is set by the algorithm configuration as the global hardware trigger time (configured in detinfo::DetectorClocks) and an offset. It may be obtained from CreateFullWaveform().

Multiple waveforms may be returned for a single channel, since a form of zero suppression is applied by the simulated hardware. The waveforms are guaranteed to be non-overlapping, non-contiguous and sorted with increasing timestamp.

The procedure attempts to mimic the working mode of CAEN V1730B boards as described on page 32 of the manual "UM2792_V1730_V1725_rev2.pdf" in SBN DocDB 15024.

Algorithm details

In the board readout language, the "trigger" is the per-channel information that the signal on the channel has passed the threshold.

It is critical to understand the details of the generation of the triggers. At the time of writing, the algorithm in FindTriggers() emits a trigger every time the signal passes from under threshold to beyond threshold. Assuming that the threshold is at less than one photoelectron, this kind of behavior implies that every scintillation photons arriving on the tail of the first (large?) signal will add a new waveform in tail. It is not clear what happens if two more triggers happen while the window from the first trigger is still open.

The assumptions in the code include:

1. overlapped triggers are not discarded (p. 32);
2. board is set for self-trigger (p. 38);
3. each waveform has a fixed duration (except the ones overlapping the previous one);
4. when a trigger starts the recording window on a channel, only the last trigger happening during that window ("overlapping"), if any, is kept.
5. at decoding time, contiguous buffers are merged in a single waveform

The fixed duration of the waveform if the sum of pre- and post-trigger window length. If during this window another trigger is emitted, a new window will follow the current one and it will be long enough to contain all the post-trigger data.

Definition at line 409 of file PMTsimulationAlg.cxx.

{
  /*
   * Plan:
   * 
   * 1. set up
   * 2. get the trigger points of the waveform
   * 3. define the size of data around each trigger to commit to waveforms
   *    (also merge contiguous and overlapping intervals)
   * 4. create the actual `raw::OpDetWaveform` objects
   * 
   */
  
  //
  // parameters check and setup
  //
  
  // not a big deal if this assertion fails, but a bit more care needs to be
  // taken in comparisons and subtractions
  static_assert(
    std::is_signed_v<optical_tick::value_t>,
    "This algorithm requires tick type to be signed."
    );
  
  using namespace detinfo::timescales; // electronics_time, time_interval, ...

  auto const pretrigSize = optical_time_ticks::castFrom(fParams.pretrigSize());
  auto const posttrigSize = optical_time_ticks::castFrom(fParams.posttrigSize());
  
  detinfo::DetectorTimings const& timings
    = detinfo::makeDetectorTimings(fParams.clockData);

  // first viable tick number: since this is the item index in `wvfm`, it's 0
  optical_tick const firstTick { 0 };

  // use hardware trigger time plus the configured offset as waveform start time
  OpDetWaveformMaker_t createOpDetWaveform {
    waveform,
    timings.TriggerTime() + time_interval{ fParams.triggerOffsetPMT },
    1.0 / fSampling
    };
  
  //
  // get the PMT channel triggers to consider
  //
  
  // prepare the set of triggers
  std::vector<optical_tick> const trigger_locations = FindTriggers(waveform);
  auto const tend = trigger_locations.end();
  
  // find the first viable trigger
  auto tooEarlyTrigger = [earliest = firstTick + pretrigSize](optical_tick t)
    { return t < earliest; };
  auto iNextTrigger
    = std::find_if_not(trigger_locations.begin(), tend, tooEarlyTrigger);
  
  //
  // collect all buffer ranges and merge them
  //
  using BufferRange_t = OpDetWaveformMaker_t::BufferRange_t;
  auto makeBuffer
    = [pretrigSize, posttrigSize](optical_tick triggerTime) -> BufferRange_t
    { return { triggerTime - pretrigSize, triggerTime + posttrigSize }; }
    ;
  
  std::vector<BufferRange_t> buffers;
  buffers.reserve(std::distance(iNextTrigger, tend)); // worst case
  
  auto earliestBufferStart { firstTick };
  while (iNextTrigger != tend) {
    
    BufferRange_t const buffer = makeBuffer(*iNextTrigger);
    
    if (buffer.first < earliestBufferStart) { // extend the previous buffer
      assert(!buffers.empty()); // guaranteed because we skipped early triggers
      buffers.back().second = buffer.second;
    }
    else buffers.emplace_back(buffer);
    
    earliestBufferStart = buffer.second;
    
    ++iNextTrigger;
  } // while
  
  //
  // turn each buffer into a waveform
  //
  MF_LOG_TRACE("PMTsimulationAlg")
    << "Channel #" << opChannel << ": " << buffers.size() << " waveforms for "
    << trigger_locations.size() << " triggers"
    ;
  std::vector<raw::OpDetWaveform> output_opdets;
  for (BufferRange_t const& buffer: buffers) {
    
    output_opdets.push_back(createOpDetWaveform(opChannel, buffer));
    
  } // for buffers
  
  return output_opdets;
} // icarus::opdet::PMTsimulationAlg::CreateFixedSizeOpDetWaveforms()

auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform ( sim::SimPhotons const &  photons, sim::SimPhotonsLite const &  lite_photons, std::optional< sim::SimPhotons > &  photons_used  ) const

private

Creates raw::OpDetWaveform objects from simulated photoelectrons.

Parameters

photons	the simulated list of photoelectrons
photons_used	(output) list of used photoelectrons

Returns

a collection of digitised raw::OpDetWaveform objects

This function performs the digitization of a optical detector channel which is collecting the photoelectrons in the photons list.

The photoelectrons are already screened for quantum efficiency: all of them are considered for use. It is still possible, if a reduction of the quantum efficiency is requested, that some of them are discarded.

The photons_used output argument is constructed and filled only if the configuration of the algorithm requires the creation of a list of used photons.

Definition at line 185 of file PMTsimulationAlg.cxx.

{

    using namespace util::quantities::time_literals;
    using namespace util::quantities::frequency_literals;
    using namespace util::quantities::electronics_literals;
    using namespace detinfo::timescales;

    detinfo::DetectorTimings const& timings
      = detinfo::makeDetectorTimings(fParams.clockData);

    tick const endSample = tick::castFrom(fNsamples);

    //
    // collect the amount of photoelectrons arriving at each subtick;
    // the waveform is split in groups of photons at the same relative subtick
    // (i.e. first all the photons on the first subtick of a tick, then
    // all the photons on the second subtick of a tick, and so on);
    // storage is by subtick group (vector index is the subtick), then by
    // tick (unordered map index is the tick number).
    //
    std::vector<std::unordered_map<tick, unsigned int>> peMaps
      (wsp.nSubsamples());

    // returns tick and relative subtick number
    TimeToTickAndSubtickConverter const toTickAndSubtick(peMaps.size());

//     auto start = std::chrono::high_resolution_clock::now();
    
    if (photons_used) {
      photons_used->clear();
      photons_used->SetChannel(photons.OpChannel());
    }
    for(auto const& ph : photons) {
      if (!KicksPhotoelectron()) continue;

      if (photons_used) photons_used->push_back(ph); // copy

      simulation_time const photonTime { ph.Time };

      trigger_time const mytime
        = timings.toTriggerTime(photonTime)
        - fParams.triggerOffsetPMT
        ;
      if ((mytime < 0.0_us) || (mytime >= fParams.readoutEnablePeriod)) continue;

      auto const [ tick, subtick ]
        = toTickAndSubtick(mytime.quantity() * fSampling);
      /*
      mf::LogTrace("PMTsimulationAlg")
        << "Photon at " << photonTime << ", optical time " << mytime
        << " => tick " << tick_d
        << " => sample " << tick << " subsample " << subtick
        ;
      */
      if (tick >= endSample) continue;
      ++peMaps[subtick][tick];
    } // for photons

//     auto end = std::chrono::high_resolution_clock::now();
//     std::chrono::duration<double> diff = end-start;
//     std::cout << "\tcollected pes... " << photons.OpChannel() << " " << diff.count() << std::endl;
//     start=std::chrono::high_resolution_clock::now();

    // Add SimPhotonsLite.  Loop over geant ticks (=ns).

    for(auto const& [ time_ns, nphotons ]: lite_photons.DetectedPhotons) {

      // Count photoelectrons.

      unsigned int nPE = 0;
      for(int i=0; i<nphotons; ++i) {
        if(KicksPhotoelectron())
          ++nPE;
      }

      // Convert photon time bin to ticks.

      simulation_time const photonTime { time_ns + 0.5 };
      trigger_time const mytime
        = timings.toTriggerTime(photonTime)
        - fParams.triggerOffsetPMT;
      if ((mytime < 0.0_us) || (mytime >= fParams.readoutEnablePeriod)) continue;

      auto const [ tick, subtick ]
        = toTickAndSubtick(mytime.quantity() * fSampling);
      if (tick < endSample)
        peMaps[subtick][tick] += nPE;
    }

    //
    // add the collected photoelectrons to the waveform
    //
    Waveform_t waveform(fNsamples, fParams.baseline);
    
    unsigned int nTotalPE [[gnu::unused]] = 0U; // unused if not in `debug` mode
    double nTotalEffectivePE [[gnu::unused]] = 0U; // unused if not in `debug` mode

    auto gainFluctuation = makeGainFluctuator();

    // go though all subsamples (starting each at a fraction of a tick)
    for (auto const& [ iSubsample, peMap ]: util::enumerate(peMaps)) {

      // this is the waveform sampling for the selected subsample:
      auto const& subsample = wsp.subsample(iSubsample);

      for (auto const& [ startTick, nPE ]: peMap) {
        nTotalPE += nPE;

        double const nEffectivePE = gainFluctuation(nPE);
        nTotalEffectivePE += nEffectivePE;

        AddPhotoelectrons(
          subsample, waveform, startTick,
          static_cast<WaveformValue_t>(nEffectivePE)
          );
        //        std::cout << "Channel=" << photons.OpChannel() << ", subsample=" << iSubsample << ", tick=" << startTick << ", nPE=" << nPE << ", ePE=" << nEffectivePE << std::endl;

      } // for sample
    } // for subsamples
    MF_LOG_TRACE("PMTsimulationAlg")
      << nTotalPE << " photoelectrons at "
      << std::accumulate(
           peMaps.begin(), peMaps.end(), 0U,
           [](unsigned int n, auto const& map){ return n + map.size(); }
         )
      << " times in channel " << photons.OpChannel()
      ;

//       end=std::chrono::high_resolution_clock::now(); diff = end-start;
//       std::cout << "\tadded pes... " << photons.OpChannel() << " " << diff.count() << std::endl;
//       start=std::chrono::high_resolution_clock::now();

      if(fParams.ampNoise > 0.0_ADCf) (this->*fNoiseAdder)(waveform);
      if(fParams.darkNoiseRate > 0.0_Hz) AddDarkNoise(waveform);

//       end=std::chrono::high_resolution_clock::now(); diff = end-start;
//       std::cout << "\tadded noise... " << photons.OpChannel() << " " << diff.count() << std::endl;
//       start=std::chrono::high_resolution_clock::now();
      
    // saturation in terms of photoelectrons (sharp);
    auto const ADCrange = fParams.ADCrange();
    ApplySaturation(waveform, ADCrange);
    
    // clip to the ADC range, 0 -- 2
    ClipWaveform(waveform, ADCrange.first, ADCrange.second);
    
//       end=std::chrono::high_resolution_clock::now(); diff = end-start;
//       std::cout << "\tadded saturation... " << photons.OpChannel() << " " << diff.count() << std::endl;
    
    return waveform;
  } // CreateFullWaveform()

auto icarus::opdet::PMTsimulationAlg::FindTriggers ( Waveform_t const &  wvfm ) const

private

Ticks in the specified waveform where some signal activity starts.

Returns

a collection of ticks with interesting activity, sorted

See Also

CreateBeamGateTriggers()

We define an "interest point" a time when some activity in the waveform is considered interesting enough to be recorded. This method returns a list of interest points, in the form of the index they are located at in the waveform wvfm.

In general (but with the exception noted below), a time becomes an interest point if the sample recorded at that time is above the threshold set by the configuration parameter ThresholdADC.

These interest points are local readout triggers that drive zero suppression on the optical readout channel and that are not necessarily causing any level of event trigger.

This method also adds the mandatory beam gate interest points as explained in CreateBeamGateTriggers() documentation. These are additional interest points that are added independently of whether there is actual interesting activity in there.

Definition at line 368 of file PMTsimulationAlg.cxx.

  {
    std::vector<optical_tick> trigger_locations;

    // find all ticks at which we would trigger readout
    bool above_thresh=false;
    for(size_t i_t=0; i_t<wvfm.size(); ++i_t){

      auto const val { fParams.pulsePolarity* (wvfm[i_t]-fParams.baseline) };

      if(!above_thresh && val>=fParams.thresholdADC){
        above_thresh=true;
        trigger_locations.push_back(optical_tick::castFrom(i_t));
      }
      else if(above_thresh && val<fParams.thresholdADC){
        above_thresh=false;
      }

    }//end loop over waveform

    // next, add the triggers injected at beam gate time
    if (fParams.createBeamGateTriggers) {
      auto beamGateTriggers = CreateBeamGateTriggers();

      // insert the new triggers and sort them
      trigger_locations.insert(trigger_locations.end(),
        beamGateTriggers.begin(), beamGateTriggers.end());
      std::inplace_merge(
        trigger_locations.begin(),
        trigger_locations.end() - beamGateTriggers.size(),
        trigger_locations.end()
        );
    }

    return trigger_locations;
  }

bool icarus::opdet::PMTsimulationAlg::KicksPhotoelectron ( ) const

private

Returns a random response whether a photon generates a photoelectron.

Definition at line 512 of file PMTsimulationAlg.cxx.

  { return CLHEP::RandFlat::shoot(fParams.randomEngine) < fQE; }

auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator ( ) const

private

Returns a configured gain fluctuator object.

Definition at line 169 of file PMTsimulationAlg.cxx.

                                                           {

  using Fluctuator_t = GainFluctuator<CLHEP::RandPoisson>;

  if (fParams.doGainFluctuations) {
    double const refGain = fParams.PMTspecs.firstStageGain();
    return Fluctuator_t
      { refGain, CLHEP::RandPoisson{ *fParams.gainRandomEngine, refGain } };
  }
  else return Fluctuator_t{}; // default-constructed does not fluctuate anything

} // icarus::opdet::PMTsimulationAlg::makeGainFluctuator()

template<typename Stream >

void icarus::opdet::PMTsimulationAlg::printConfiguration	(	Stream &&	out,
		std::string	indent = ""
	)		const

Prints the configuration into the specified output stream.

Definition at line 1064 of file PMTsimulationAlg.h.

{
  using namespace util::quantities::electronics_literals;

  out
            << indent << "ADC bits:            " << fParams.ADCbits
      << " (" << fParams.ADCrange().first << " -- " << fParams.ADCrange().second
      << ")"
    << '\n' << indent << "Baseline:            " << fParams.baseline
    << '\n' << indent << "ReadoutWindowSize:   " << fParams.readoutWindowSize << " ticks"
    << '\n' << indent << "PreTrigFraction:     " << fParams.pretrigFraction
    << '\n' << indent << "ThresholdADC:        " << fParams.thresholdADC
    << '\n' << indent << "Saturation:          " << fParams.saturation << " p.e."
    << '\n' << indent << "doGainFluctuations:  "
      << std::boolalpha << fParams.doGainFluctuations
    << '\n' << indent << "PulsePolarity:       " << ((fParams.pulsePolarity == 1)? "positive": "negative") << " (=" << fParams.pulsePolarity << ")"
    << '\n' << indent << "Sampling:            " << fSampling;
  if (fParams.pulseSubsamples > 1U)
    out << " (subsampling: x" << fParams.pulseSubsamples << ")";
  out
    << '\n' << indent << "Samples/waveform:    " << fNsamples << " ticks"
    << '\n' << indent << "Gain at first stage: " << fParams.PMTspecs.firstStageGain()
    ;

  out << '\n' << indent << "Electronics noise:   ";
  if (fParams.ampNoise > 0_ADCf) {
    out << fParams.ampNoise << " RMS ("
      << (fParams.useFastElectronicsNoise? "faster": "slower") << " algorithm)";
  }
  else out << "none";

  if (fParams.createBeamGateTriggers) {
    out << '\n' << indent << "Create " << fParams.beamGateTriggerNReps
      << " beam gate triggers, one every " << fParams.beamGateTriggerRepPeriod << ".";
  }
  else out << '\n' << indent << "Do not create beam gate triggers.";

  out << '\n' << indent << "... and more.";

  out << '\n' << indent << "Template photoelectron waveform settings:"
    << '\n';
  wsp.dump(std::forward<Stream>(out), indent + "  ");
  
  out << '\n' << indent << "Track used photons:  "
    << std::boolalpha << fParams.trackSelectedPhotons
    << '\n';
} // icarus::opdet::PMTsimulationAlg::printConfiguration()

auto icarus::opdet::PMTsimulationAlg::saturationRange ( ) const

private

Returns the ADC range allowed for photoelectron saturation.

Definition at line 653 of file PMTsimulationAlg.cxx.

{
  std::pair<ADCcount, ADCcount> range = fParams.ADCrange();
  if (fParams.saturation > 0) {
    ADCcount const saturationLevel
      = fParams.baseline + fParams.saturation*wsp.peakAmplitude();
    ((fParams.pulsePolarity > 0)? range.second: range.first) = saturationLevel;
  }
  return range;
} // icarus::opdet::PMTsimulationAlg::saturationRange()

std::tuple< std::vector< raw::OpDetWaveform >, std::optional< sim::SimPhotons > > icarus::opdet::PMTsimulationAlg::simulate	(	sim::SimPhotons const &	photons,
		sim::SimPhotonsLite const &	lite_photons
	)

Returns the waveforms originating from simulated photons.

Parameters

photons

all the photons simulated to land on the channel

Returns

a list of optical waveforms, response to those photons, and which photons were used (if requested)

Due to threshold readout, a single channel may result in multiple waveforms, which are all on the same channel but disjunct in time.

The second element of the return value is optional and filled only if the trackSelectedPhotons configuration parameter is set to true. In that case, the returned sim::SimPhotons contains a copy of each of the photons contributing to any of the waveforms.

Definition at line 153 of file PMTsimulationAlg.cxx.

{
  std::optional<sim::SimPhotons> photons_used;

  Waveform_t const waveform = CreateFullWaveform(photons, lite_photons, photons_used);

  return {
    CreateFixedSizeOpDetWaveforms(photons.OpChannel(), waveform),
    std::move(photons_used)
    };
  
} // icarus::opdet::PMTsimulationAlg::simulate()

Member Data Documentation

util::FastAndPoorGauss< 32768U, float > const icarus::opdet::PMTsimulationAlg::fFastGauss

staticprivate

Definition at line 589 of file PMTsimulationAlg.h.

NoiseAdderFunc_t const icarus::opdet::PMTsimulationAlg::fNoiseAdder

private

Single photon pulse (sampled).

Selected electronics noise method. Transformation uniform to Gaussian for electronics noise.

Definition at line 586 of file PMTsimulationAlg.h.

std::size_t icarus::opdet::PMTsimulationAlg::fNsamples

private

Samples per waveform.

Definition at line 582 of file PMTsimulationAlg.h.

ConfigurationParameters_t icarus::opdet::PMTsimulationAlg::fParams

private

Complete algorithm configuration.

Definition at line 578 of file PMTsimulationAlg.h.

double icarus::opdet::PMTsimulationAlg::fQE

private

PMT quantum efficiency.

Definition at line 580 of file PMTsimulationAlg.h.

megahertz icarus::opdet::PMTsimulationAlg::fSampling

private

Wave sampling frequency [MHz].

Definition at line 581 of file PMTsimulationAlg.h.

DiscretePhotoelectronPulse icarus::opdet::PMTsimulationAlg::wsp

private

Definition at line 584 of file PMTsimulationAlg.h.

---

The documentation for this class was generated from the following files:

• PMTsimulationAlg.h
• PMTsimulationAlg.cxx