Extracts a baseline from PMT waveforms. More...

Inheritance diagram for icarus::PMTWaveformBaselinesFromChannelData:

Classes
struct	BaselineInfo_t
	Record of baseline information to be written. More...

struct	Config

Public Types
using	Parameters = art::SharedProducer::Table< Config >

Public Member Functions
	PMTWaveformBaselinesFromChannelData (Parameters const &config, art::ProcessingFrame const &)

virtual void	beginJob (art::ProcessingFrame const &frame) override
	Prepares the plots to be filled. More...

virtual void	beginRun (art::Run &run, art::ProcessingFrame const &) override
	Reads needed PMT configuration. More...

virtual void	produce (art::Event &event, art::ProcessingFrame const &) override
	Creates the data products. More...

virtual void	endJob (art::ProcessingFrame const &frame) override
	Remove empty plots. More...

Private Member Functions
unsigned int	removeWaveformsAround (std::vector< raw::OpDetWaveform const * > &waveforms, double time) const
	Removes `waveforms` containing `time`, retuning how many were removed. More...

unsigned int	getPretriggerBuffer (sbn::PMTconfiguration const &PMTconfig) const
	Returns the smallest pre-trigger buffer size available among all boards. More...

void	setupPlots (art::ProcessingFrame const &frame)
	Creates all the plots to be filled by the module. More...

void	buildBaselineGraphs (art::ProcessingFrame const &frame)
	Removes the empty plots. More...

Static Private Member Functions
static std::vector < std::vector < raw::OpDetWaveform const * > >	groupByChannel (std::vector< raw::OpDetWaveform > const &waveforms)

Private Attributes
art::InputTag const	fOpDetWaveformTag
	Input optical waveform tag. More...

unsigned int	fExcludeSpillTimeIfMoreThan
	Minimum number of waveforms when excluding trigger one is allowed. More...

unsigned int	fPretriggerSamples
	Samples in the pre-trigger buffer. More...

std::optional< art::InputTag >	fPMTconfigTag
	Tag for PMT configuration data product. More...

float const	fSampleFraction
	Fraction of pretrigger buffer to use. More...

opdet::SharedWaveformBaseline::Params_t	fAlgoParams
	Parameters for the baseline algorithm. More...

bool const	fPlotBaselines
	Whether to produce plots. More...

double const	fBaselineTimeAverage { 0.0 }
	Width of baseline time profile binning [s]. More...

std::string const	fLogCategory
	Category name for the console output stream. More...

util::quantities::intervals::nanoseconds	fOpticalTick
	Duration of a PMT digitizer tick. More...

std::size_t	fNPlotChannels = 0U
	Number of plotted channels. More...

TH2 *	fHBaselines = nullptr
	All baselines, per channel. More...

std::vector< std::vector < std::pair< double, double > > >	fBaselinesVsTime
	For each channel, all event times and their baselines. More...

Detailed Description

Extracts a baseline from PMT waveforms.

This module produces a baseline data product for each optical detector waveform. The algorithm extracts a baseline per channel per event, considering all the waveforms from one channel together, under the assumptions that:

the baseline is not going to change during the few milliseconds of readout
the beginning of the waveform is on its baseline level

The algorithm

The core of the algorithm is described in its own class, opdet::SharedWaveformBaseline.

On each event independently, the waveforms are grouped by channel. For each waveform, opdet::SharedWaveformBaseline considers their first part for baseline calculation. The number of samples of that part is specified as a fraction of the pre-trigger buffer. The pre-trigger buffer includes the samples that were collected before the physics activity that causes the data acquisition happened, and as such is expected to be almost completely free of physics signal and to be made of just electronics noise. The size of this prebuffer is determined in data by the readout configuration, which can be accessed by specifying the parameter PMTconfigurationTag, and in simulation by a digitization module parameter, which can be replicated here by the parameter PretriggerBufferSize. The fraction of such buffer used for the baseline estimation is also specified by a configuration parameter (PretriggerBufferFractionForBaseline). Depending on the value of the configuration parameter ExcludeSpillTimeIfMoreThan, the one waveform that covers the trigger time may be excluded and not used; the rationale is that there is a conspicuous number of events triggered by the late light of a cosmic ray happening before the beam spill, in which case the activity may contaminate the pre-trigger buffer and could bias the estimation of the baseline.

The result of the opdet::SharedWaveformBaseline is currently used directly as the baseline for all the waveforms on the channel on that event.

Todo:: Add run-level information and checks.

Output data products

a data product of type std::vector<icarus::WaveformBaseline>, with one baseline per input waveform; the baselines are guaranteed to be in the same order as the waveforms in the input collection;
a data product of type std::vector<icarus::WaveformRMS>, with one baseline RMS per input waveform; the baselines are guaranteed to be in the same order as the waveforms in the input collection;
two art associations, with one baseline and one RMS per input optical waveform; normally this is not needed as long as the original PMT waveform data product is available; the type of the associations are art::Assns<icarus::WaveformBaseline, raw::OpDetWaveform> and art::Assns<icarus::WaveformRMS, raw::OpDetWaveform>.

Output plots

Baselines: baseline distribution, per channel; average baseline [ADC] per event per channel; all waveforms on the same channels in a single event contribute to the average, and channels with no waveforms in an event do not contribute an entry for that event.

Input data products

std::vector<raw::OpDetWaveform>: a single waveform for each recorded optical detector activity; the activity belongs to a single channel, but there may be multiple waveforms on the same channel.

Service requirements

DetectorClocksService for the determination of the optical tick duration and the trigger time (if the exclusion of waveform with trigger is enabled).
TFileService and Geometry if PlotBaselines is enabled.

Configuration parameters

A terse description of the parameters is printed by running lar --print-description PMTWaveformBaselinesFromChannelData.

OpticalWaveforms (input tag, mandatory): the data product containing all optical detector waveforms;
PretriggerBufferSize (positive integer, optional): the number of samples collected by the PMT readout before the PMT signal. This number is used as base for the size of the initial part of the waveform to be used in baseline calculation. This and PMTconfigurationTag parameters are exclusive. This one is expected to be preferred for simulated samples.
PMTconfigurationTag (input tag): the (run) data product containing the PMT readout configuration. This is extracted from the run configuration for data events. The pre-trigger buffer size is read from this data product, and used as base for the size of the initial part of the waveform to be used in baseline calculation. This and PretriggerBufferSize parameters are exclusive. This one is expected to be preferred for data samples.
PretriggerBufferFractionForBaseline (real, default: 0.5): the fraction of the pre-trigger buffer (see PretriggerBufferSize and PMTconfigurationTag parameters) to be used to extract the baseline.
ExcludeSpillTimeIfMoreThan (integer, default: disable): the minimum number of PMT waveforms on a single channel of a single event, in order for the exclusion of the waveform containing the global trigger to be enabled. By default (huge number) the feature is disabled.
AlgoParams (table): configuration of the core algorithm extracting the baseline from a set of waveforms on the same channel. The configuration is as follows:
- AcceptedSampleRangeRMS (real): for a waveform to be considered for the baseline, the values of the samples must stay within a certain range, which is defined as this number of RMS away in either directions from a central value that is a rougher estimation of the baseline.
- ExcessSampleLimit (integer): for a waveform to be considered for the baseline, there must be less than this number of samples in a row that are outside of the range defined by AcceptedSampleRangeRMS parameter.
PlotBaselines (flag, default: true): whether to produce distributions of the extracted baselines.
BaselineTimeAverage (real number, default: 600.0): binning of the baseline profile vs. time, in seconds. Requires PlotBaselines to be set.
OutputCategory (string, default: "PMTWaveformBaselinesFromChannelData"): label for the category of messages in the console output; this is the label that can be used for filtering messages via MessageFacility service configuration.

Definition at line 206 of file PMTWaveformBaselinesFromChannelData_module.cc.

Member Typedef Documentation

using icarus::PMTWaveformBaselinesFromChannelData::Parameters = art::SharedProducer::Table<Config>

Definition at line 297 of file PMTWaveformBaselinesFromChannelData_module.cc.

Constructor & Destructor Documentation

icarus::PMTWaveformBaselinesFromChannelData::PMTWaveformBaselinesFromChannelData	(	Parameters const &	config,
		art::ProcessingFrame const &	frame
	)

explicit

Definition at line 462 of file PMTWaveformBaselinesFromChannelData_module.cc.

   : art::SharedProducer(config)
   // configuration
   , fOpDetWaveformTag(config().OpticalWaveforms())
   , fExcludeSpillTimeIfMoreThan(config().ExcludeSpillTimeIfMoreThan())
   , fPretriggerSamples(config().PretriggerBufferSize().value_or(0U))
   , fPMTconfigTag(config().PMTconfigurationTag())
   , fSampleFraction(config().PretriggerBufferFractionForBaseline())
   , fAlgoParams(config().AlgoParams())
   , fPlotBaselines(config().PlotBaselines())
   , fBaselineTimeAverage(config().BaselineTimeAverage())
   , fLogCategory(config().OutputCategory())
   // service caching
   , fOpticalTick(
       detinfo::DetectorTimings
         { frame.serviceHandle<detinfo::DetectorClocksService>()->DataForJob() }
         .OpticalClockPeriod()
     )
   // algorithms
 {
   
   if (fPlotBaselines)
 //     serialize<art::InEvent>(art::TFileService::resource_name()); // TODO isn't art supposed to provide this method?
       serializeExternal<art::InEvent>(std::string{ "TFileService" });
   else
     async<art::InEvent>();
   
   
   //
   // optional configuration parameters
   //
   if (fPMTconfigTag.has_value() == (fPretriggerSamples != 0)) {
     throw art::Exception(art::errors::Configuration)
       << "Exactly one between parameters '"
       << config().PMTconfigurationTag.name() << "' and '"
       << config().PretriggerBufferFractionForBaseline.name()
       << "' (with a positive value) must be specified!\n";
   }
   
   //
   // configuration report
   //
   {
     mf::LogInfo log{ fLogCategory };
     log << "Using the standard (median) algorithm, waveform by waveform, on '"
       << fOpDetWaveformTag.encode() << "'";
   }
   
   //
   // declaration of input
   //
   consumes<std::vector<raw::OpDetWaveform>>(fOpDetWaveformTag);
   if (fPMTconfigTag) consumes<sbn::PMTconfiguration>(*fPMTconfigTag);
   
   //
   // declaration of output
   //
   produces<std::vector<icarus::WaveformBaseline>>();
   produces<std::vector<icarus::WaveformRMS>>();
   produces<art::Assns<icarus::WaveformBaseline, raw::OpDetWaveform>>();
   produces<art::Assns<icarus::WaveformRMS, raw::OpDetWaveform>>();
   
 } // icarus::PMTWaveformBaselinesFromChannelData::PMTWaveformBaselinesFromChannelData()

Member Function Documentation

void icarus::PMTWaveformBaselinesFromChannelData::beginJob ( art::ProcessingFrame const & frame )

overridevirtual

Prepares the plots to be filled.

Definition at line 529 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   //
   // set up the plots, if needed
   //
   if (fPlotBaselines) setupPlots(frame);
   
 } // icarus::PMTWaveformBaselinesFromChannelData::beginJob()

void icarus::PMTWaveformBaselinesFromChannelData::beginRun	(	art::Run &	run,
		art::ProcessingFrame const &
	)

overridevirtual

Reads needed PMT configuration.

Definition at line 542 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   if (fPMTconfigTag) { // update the pretrigger samples number
     fPretriggerSamples = getPretriggerBuffer
       (run.getProduct<sbn::PMTconfiguration>(*fPMTconfigTag));
   }
   assert(fPretriggerSamples > 0U);
   
   fAlgoParams.nSample = std::max(
     1U,
     static_cast<unsigned int>(std::round(fSampleFraction * fPretriggerSamples))
     );
   
   {
     // not clear whether this is debug or info level
     mf::LogInfo log{ fLogCategory };
     log << run.id() << ": baseline algorithm configuration:\n";
     fAlgoParams.dump(log, " - ");
   }
   
 } // icarus::PMTWaveformBaselinesFromChannelData::beginRun()

void icarus::PMTWaveformBaselinesFromChannelData::buildBaselineGraphs ( art::ProcessingFrame const & frame )

private

Removes the empty plots.

Definition at line 799 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   auto& tfs = *(frame.serviceHandle<art::TFileService>());
   
   art::TFileDirectory graphDir = tfs.mkdir("graphs", "Baseline vs. time");
   art::TFileDirectory profileDir
     = tfs.mkdir("profiles", "Baseline profiles vs. time");
   
   for (auto const channel: util::counter(fBaselinesVsTime.size())) {
     
     auto& timeAndBaselines = fBaselinesVsTime[channel];
     if (timeAndBaselines.empty()) continue;
     
     // sort by time (entries with the same time would be sorted by baseline,
     // but that does not really happen nor it mattered if it happened)
     std::sort(timeAndBaselines.begin(), timeAndBaselines.end());
     
     // graph, one point per event
     auto* const graph = graphDir.makeAndRegister<TGraph>(
       "BaselineCh" + std::to_string(channel),
       "PMT channel #" + std::to_string(channel) + ": baseline vs. run time",
       timeAndBaselines.size()
       );
     assert(graph->GetXaxis());
     assert(graph->GetYaxis());
     graph->GetXaxis()->SetTitle("event time");
     graph->GetYaxis()->SetTitle("baseline  [ ADC ]");
     
     
     // profile, one point every 10 minutes (or best offer)
     double const startTime = timeAndBaselines.front().first;
     double const totalTime = timeAndBaselines.back().first - startTime;
     auto const nBins = std::max(1U,
       static_cast<unsigned int>(std::ceil(totalTime / fBaselineTimeAverage))
       );
     double const endTime = startTime + nBins * fBaselineTimeAverage;
     
     auto* const profile = profileDir.make<TProfile>(
       ("BaselineCh" + std::to_string(channel) + "profile").c_str(),
       ("PMT channel #" + std::to_string(channel) + ": baseline vs. run time")
         .c_str(),
       nBins, startTime, endTime
       );
     assert(profile->GetXaxis());
     assert(profile->GetYaxis());
     profile->GetXaxis()->SetTitle(
       ("event time  [ / " + std::to_string(fBaselineTimeAverage) + " s ]")
       .c_str()
       );
     profile->GetYaxis()->SetTitle("average baseline  [ ADC ]");
     
     for (auto const& [ i, data ]: util::enumerate(timeAndBaselines)) {
       auto const [ time, baseline ] = data;
       graph->SetPoint(i, time, baseline);
       profile->Fill(time, baseline);
     } // for
     
   } // for
   
 } // icarus::PMTWaveformBaselinesFromChannelData::buildBaselineGraphs()

void icarus::PMTWaveformBaselinesFromChannelData::endJob ( art::ProcessingFrame const & frame )

overridevirtual

Remove empty plots.

Definition at line 698 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   if (fPlotBaselines) buildBaselineGraphs(frame);
   
 } // icarus::PMTWaveformBaselinesFromChannelData::endJob()

unsigned int icarus::PMTWaveformBaselinesFromChannelData::getPretriggerBuffer ( sbn::PMTconfiguration const & PMTconfig ) const

private

Returns the smallest pre-trigger buffer size available among all boards.

Definition at line 734 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   lar::util::MinMaxCollector<unsigned int> prebuffer;
   for (sbn::V1730Configuration const& board: PMTconfig.boards)
     prebuffer.add(board.preTriggerTicks());
   if (!prebuffer.has_data()) {
     throw art::Exception{ art::errors::Unknown }
       << "No boards found in the PMT configuration!\n";
   }
   else if (prebuffer.max() > prebuffer.min()) {
     throw art::Exception{ art::errors::Unknown }
       << "Found different pre-trigger readout buffer sizes between "
       << prebuffer.min() << " and " << prebuffer.max() << " ticks.\n";
   }
   return prebuffer.min();
 } // icarus::PMTWaveformBaselinesFromChannelData::getPretriggerBuffer()

std::vector< std::vector< raw::OpDetWaveform const * > > icarus::PMTWaveformBaselinesFromChannelData::groupByChannel ( std::vector< raw::OpDetWaveform > const & waveforms )

staticprivate

Returns a common baseline for all the specified waveforms. Returns a vector of waveforms per channel (which is the index).

Definition at line 755 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   std::vector<std::vector<raw::OpDetWaveform const*>> groups;
   VectorExpander groupForChannel { groups };
   
   for (raw::OpDetWaveform const& waveform: waveforms) {
     groupForChannel(waveform.ChannelNumber()).push_back(&waveform);
   }
   
   return groups;
 } // icarus::PMTWaveformBaselinesFromChannelData::groupByChannel()

void icarus::PMTWaveformBaselinesFromChannelData::produce	(	art::Event &	event,
		art::ProcessingFrame const &	frame
	)

overridevirtual

Creates the data products.

Definition at line 568 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   detinfo::DetectorTimings const detTimings{
     frame.serviceHandle<detinfo::DetectorClocksService const>()->DataFor(event)
     };
   
   detinfo::timescales::electronics_time const triggerTime
      = detTimings.TriggerTime();
   
   mf::LogDebug{ fLogCategory }
     << event.id() << " trigger time: " << triggerTime;
   
   //
   // fetch input
   //
   auto const& waveformHandle
     = event.getValidHandle<std::vector<raw::OpDetWaveform>>(fOpDetWaveformTag);
   auto const& waveforms = *waveformHandle;
   
   //
   // compute all the baselines (and plot them)
   //
   opdet::SharedWaveformBaseline const sharedWaveformBaselineAlgo
     { fAlgoParams, fLogCategory };
   
   double const eventTime = static_cast<double>(event.time().timeHigh())
     + static_cast<double>(event.time().timeHigh()) * 1e-9;
   
   std::vector<BaselineInfo_t> channelBaselines;
   VectorExpander baselineForChannel { channelBaselines };
   
   auto waveformsByChannel = groupByChannel(waveforms);
   
   for (auto const& [ channel, waveforms ]: util::enumerate(waveformsByChannel))
   {
     if (waveforms.empty()) continue;
     
     mf::LogTrace{ fLogCategory }
       << "Processing " << waveforms.size() << " waveforms for channel "
       << channel;
     
     //
     // remove global trigger waveform
     //
     if (waveforms.size() >= fExcludeSpillTimeIfMoreThan) {
       
       unsigned int const nExcluded
         = removeWaveformsAround(waveforms, triggerTime.value());
       if (nExcluded > 0U) {
         mf::LogTrace{ fLogCategory }
           << "Removed " << nExcluded << "/" << (waveforms.size() + nExcluded)
           << " waveforms at trigger time " << triggerTime;
       }
       
     } // if many waveforms
     
     //
     // extract baseline
     //
     opdet::SharedWaveformBaseline::BaselineInfo_t const baseline
       = sharedWaveformBaselineAlgo(waveforms);
 
     mf::LogTrace{ fLogCategory }
       << "Channel " << channel << ": baseline " << baseline.baseline
       << " ADC# from " << baseline.nSamples << " samples in "
       << baseline.nWaveforms << "/" << waveforms.size()
       << " waveforms; found RMS=" << baseline.RMS << " ADC#";
     
     auto const chIndex = static_cast<std::size_t>(channel);
     baselineForChannel(chIndex) = { baseline.baseline, baseline.RMS };
     
     if (fHBaselines) fHBaselines->Fill(double(channel), baseline.baseline);
     if (chIndex < fBaselinesVsTime.size())
       fBaselinesVsTime[chIndex].emplace_back(eventTime, baseline.baseline);
     
   } // for grouped waveforms
   
   //
   // assign baselines to waveforms
   //
   
   std::vector<icarus::WaveformBaseline> baselines;
   baselines.reserve(waveforms.size());
   art::Assns<icarus::WaveformBaseline, raw::OpDetWaveform> baselineToWaveforms;
   art::PtrMaker<icarus::WaveformBaseline> const makeBaselinePtr(event);
   
   std::vector<icarus::WaveformRMS> RMSs;
   RMSs.reserve(waveforms.size());
   art::Assns<icarus::WaveformRMS, raw::OpDetWaveform> RMStoWaveforms;
   art::PtrMaker<icarus::WaveformRMS> const makeRMSPtr(event);
   
   for (auto const& [ iWaveform, waveform ]: util::enumerate(waveforms)) {
     
     BaselineInfo_t baselineInfo
       = channelBaselines[static_cast<std::size_t>(waveform.ChannelNumber())];
     
     art::Ptr<raw::OpDetWaveform> const waveformPtr{ waveformHandle, iWaveform };
     
     baselines.push_back(baselineInfo.baseline);
     RMSs.push_back(baselineInfo.RMS);
     baselineToWaveforms.addSingle(makeBaselinePtr(iWaveform), waveformPtr);
     RMStoWaveforms.addSingle(makeRMSPtr(iWaveform), waveformPtr);
     
   } // for
   
   //
   // output
   //
   event.put(
     std::make_unique<std::vector<icarus::WaveformBaseline>>
       (std::move(baselines))
     );
   event.put
     (std::make_unique<std::vector<icarus::WaveformRMS>>(std::move(RMSs)));
   event.put(
     std::make_unique<art::Assns<icarus::WaveformBaseline, raw::OpDetWaveform>>
       (std::move(baselineToWaveforms))
     );
   event.put(
     std::make_unique<art::Assns<icarus::WaveformRMS, raw::OpDetWaveform>>
       (std::move(RMStoWaveforms))
     );
   
   
 } // icarus::PMTWaveformBaselinesFromChannelData::produce()

unsigned int icarus::PMTWaveformBaselinesFromChannelData::removeWaveformsAround	(	std::vector< raw::OpDetWaveform const * > &	waveforms,
		double	time
	)		const

private

Removes waveforms containing time, retuning how many were removed.

Definition at line 771 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   // tick duration in the same unit as the electronics time scale times
   // (that is microseconds, actually)
   double const tickDuration
     = detinfo::timescales::electronics_time::interval_t{ fOpticalTick }.value();
   
   auto const hasTriggerTime
     = [tickDuration,time](raw::OpDetWaveform const* waveform)
     {
       double const relTime = time - waveform->TimeStamp();
       return (relTime >= 0.0) && (relTime < (waveform->size() * tickDuration));
     };
   
   auto const iLast
     = std::remove_if(waveforms.begin(), waveforms.end(), hasTriggerTime);
   
   unsigned int const nRemoved = std::distance(iLast, waveforms.end());
   waveforms.erase(iLast, waveforms.end());
   
   return nRemoved;
   
 } // icarus::PMTWaveformBaselinesFromChannelData::removeWaveformsAround()

void icarus::PMTWaveformBaselinesFromChannelData::setupPlots ( art::ProcessingFrame const & frame )

private

Creates all the plots to be filled by the module.

Definition at line 708 of file PMTWaveformBaselinesFromChannelData_module.cc.

 {
   
   auto const& tfs = *(frame.serviceHandle<art::TFileService>());
   auto const& geom = *(frame.serviceHandle<geo::Geometry>()->provider());
   
   fNPlotChannels = geom.NOpChannels();
   
   fHBaselines = tfs.make<TH2F>(
     "Baselines",
     "PMT baseline;channel;baseline per channel [ / 8 ADC ]",
     fNPlotChannels, 0.0, double(fNPlotChannels),
     256, 13312.0, 15360.0
     );
   
   // these are graphs, and it is more convenient to carry around their data
   // in a vector than carrying around the graphs themselves;
   // `buildBaselineGraphs()` will turn that data into graph at end of the job;
   // here we just declare which channels we are going to plot.
   fBaselinesVsTime.resize(fNPlotChannels);
   
 } // icarus::PMTWaveformBaselinesFromChannelData::setupPlots()