diff --git a/CMakeLists.txt b/CMakeLists.txt index a6abcc2..411d0b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,8 @@ find_package(TBB REQUIRED) find_package(spdlog REQUIRED) find_package(phlex REQUIRED) +add_subdirectory(gaus_hit_finder) + # Proxies for framework-agnostic experiment libraries # N.B. The specified library type should be SHARED. add_library(my_add SHARED my_add.cpp) diff --git a/gaus_hit_finder/CMakeLists.txt b/gaus_hit_finder/CMakeLists.txt new file mode 100644 index 0000000..b5b3ee0 --- /dev/null +++ b/gaus_hit_finder/CMakeLists.txt @@ -0,0 +1,23 @@ +# Proxies for framework-agnostic experiment libraries +# N.B. The specified library type should be SHARED. +add_library(find_hits_with_gaussians SHARED + find_hits_with_gaussians.cpp + copied_from_larsoft_minor_edits/Hit.cxx + copied_from_larsoft_minor_edits/geo_types.cxx + copied_from_larsoft_minor_edits/CandHitStandard.cxx + copied_from_larsoft_minor_edits/PeakFitterMrqdt.cxx + copied_from_larsoft_minor_edits/HitFilterAlg.cxx + copied_from_larsoft_minor_edits/MarqFitAlg.cxx + copied_from_larsoft_minor_edits/Wire.cxx + wire_serialization.cpp + print_hits.cpp +) +target_link_libraries(find_hits_with_gaussians PRIVATE TBB::tbb fmt::fmt) + +# Register providers for data products +add_library(wires_source MODULE wires_source.cpp) +target_link_libraries(wires_source PRIVATE phlex::module find_hits_with_gaussians) + +# Register experiment libraries with Phlex +add_library(find_hits_with_gaussians_hof MODULE register_find_hits_with_gaussians.cpp) +target_link_libraries(find_hits_with_gaussians_hof PRIVATE phlex::module find_hits_with_gaussians) diff --git a/gaus_hit_finder/README.md b/gaus_hit_finder/README.md new file mode 100644 index 0000000..7dced91 --- /dev/null +++ b/gaus_hit_finder/README.md @@ -0,0 +1,203 @@ +# Purpose of this directory + +This directory holds an example showing how to convert a +module that works with the `art` framework into an +algorithm that works with the `phlex` framework. This +particular example converts the module defined in +the file GausHitFinder_module.cc from the LArSoft larreco +repository. + +This is only an example and a prototype to be used in +testing. It is not intended to be the version of this +module used by DUNE or any other experiment +for production or physics. That belongs in some other +repository. + +This example is based on version v0.1.0 of `phlex`. It will need +some minor modifications to work with the newest version of `phlex`. +The migration to phlex started with version of GausHitFinder_module.cc +in v10_05_00 of larreco which was the version in use with DUNE software +at the time the migration work was started. + +# Work in Progress + +This is work in progress. It is not final. It is intended +that this example will change as we learn more about using +`phlex` and as we receive feedback from DUNE. Currently, this +is more like a rough first draft of our first attempts +to use `phlex`. Please don't expect stability. Many features +of `phlex` have not been implemented and the design of `phlex` +itself is changing. Those changes will also cause this example +to change. + +As we study this example, it will guide how phlex is used. +As we run tests, it may point to problems in phlex and this +may guide which areas developers should focus on to fix +issues and remove performance bottlenecks. + +# Files in this directory + +## Interesting files that are the heart of the example + +1. find_hits_with_gaussians.hpp +2. find_hits_with_gaussians.cpp +3. register_find_hits_with_gaussians.cpp +4. test_find_hits_with_gaussians.jsonnet + +## These files exist only for test purposes + +There is one function that will print out reconstructed hits. +This printout can be used to compare results between an +`art` process and a `phlex` process running GausHitFinder. +There are functions that can be used to persistently store +the input from an `art` process (a vector of `Wire`) so +it can be used in the `phlex` process. And there is a `phlex` provider +to read this input. This persistence mechanism is necessary +because there is not an alternative yet. These files are all +temporary and not part of the example. + +1. print_hits.hpp +2. print_hits.cpp +3. wire_serialization.hpp +4. wire_serialization.cpp +5. wires_source.cpp + +## Files copied from LArSoft + +The infrastructure does not yet exist to link the +code here with LArSoft code or include headers. +Necessary files are copied in with minimal changes. +Some changes are necessary because the files are in +a different location and some remove unwanted dependences, +but the modifications are small. When it is possible +to depend on LArSoft code in some other repository, +these files should be deleted to remove duplication. + +These files are in the subdirectory: + + ```phlex-examples/gaus_hit_finder/copied_from_larsoft_minor_edits``` + +# Where more work is needed + +At the moment, this is implemented as a single `phlex` +transform. Does it make sense to split it with unfolds +and folds? Does it make sense to have multiple transforms? +We don't know the answers to these questions yet. +Getting the single transform to work is a good first step. + +The next thing we plan to work on is creating a second version +of the migrated module. We will try replacing the +outer `parallel_for` with an `unfold`, `transform`, and `fold` +sequence of algorithms. After that we will create other versions +pushing the division into separate algorithms to lower levels. +Then we plan to run tests and compare the different versions. + +`Phlex` will not support the `Tools` feature that existed in +`art`. One possibility is that algorithm nodes scheduled +by `tbb::flow_graph` offer sufficient configurability that we +don't need a separate plugin system like `Tools`. +Another possibility is that we eventually implement a +plugin system for `phlex`. For now, I replaced the two `Tool` +types used in `GausHitFinder` with direct instantiations of one +of the plugin types. I instantiated a vector of objects of type +`CandHitStandard` and one object of type `PeakFitterMrqdt`. They +needed some modification but are as close as reasonably possible +to the original versions. + +I do not know yet how to deal with the fact that `GausHitFinder` +can be configured to produce 1 or 2 output data products. +For now, the example can only produce one `std::vector`. +It's configurable whether to filter hits or not, but it cannot +produce both outputs at the same time. + +For now, this ignores the fact that `GausHitFinder` also can +produce `art::Assns` and +`art::Assns`. In the `art` version, +it is configurable whether to produce those or not. +We are ignoring this because `phlex` does not yet support `art::Assns` +or whatever will be used to replace that type. For now, the +example can only produce one `std::vector`. + +Here is a list of the files that were used as a basis +for this example, but significantly modified from the +LArSoft version: + +``` +https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/GausHitFinder_module.cc +``` + +Some files were copied with minimal modifications. The "tool" +classes were split into .h and .cxx files and are no longer implemented +as `art Tools`, but otherwise the changes are minimal. +Here is a list of the files that were copied from LArSoft +with minimal modifications: + +``` + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFinderTools/CandHitStandard_tool.cc + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFinderTools/PeakFitterMrqdt_tool.cc + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFilterAlg.h + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFilterAlg.cxx + https://github.com/LArSoft/lardataobj/blob/develop/lardataobj/RecoBase/Wire.h + https://github.com/LArSoft/lardataobj/blob/develop/lardataobj/RecoBase/Wire.cxx + https://github.com/LArSoft/lardataobj/blob/develop/lardataobj/RecoBase/Hit.h + https://github.com/LArSoft/lardataobj/blob/develop/lardataobj/RecoBase/Hit.cxx + https://github.com/LArSoft/larcoreobj/blob/develop/larcoreobj/SimpleTypesAndConstants/geo_types.h + https://github.com/LArSoft/larcoreobj/blob/develop/larcoreobj/SimpleTypesAndConstants/geo_types.cxx + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFinderTools/ICandidateHitFinder.h + https://github.com/LArSoft/larcoreobj/blob/develop/larcoreobj/SimpleTypesAndConstants/RawTypes.h + https://github.com/LArSoft/lardataobj/blob/develop/lardataobj/Utilities/sparse_vector.h + https://github.com/LArSoft/larvecutils/blob/develop/larvecutils/MarqFitAlg/MarqFitAlg.h + https://github.com/LArSoft/larvecutils/blob/develop/larvecutils/MarqFitAlg/MarqFitAlg.cxx + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFinderTools/ICandidateHitFinder.h + https://github.com/LArSoft/larreco/blob/develop/larreco/HitFinder/HitFinderTools/IPeakFitter.h +``` + +For now, the part of `GausHitFinder` that fills 2 histograms is removed. +Why? Because `phlex` does not yet have a histogramming facility like `TFileService` +yet. Also those histograms were booked in `beginJob`. That entire method +was removed. Booking the histograms was the only thing that was done in `beginJob`. +I'm not sure how to support `beginJob` activity in `phlex` either. + +Dropped `fEventCount`. It was not used. Seems like not a good thing +to keep track of in a `phlex` algorithm anyway. Reconstruction algorithms +should not depend on event count. Also removed the related count argument +from the findHitCandidates method of ICandidateHitFinder and its implementations. +It was not used in CandHitStandard anyway. + +MessageLogger calls are commented out for now but not deleted. Eventually +we will want to replace those with calls to whatever logging system `phlex` +eventually uses or delete those lines entirely. + +For now, calls to the Geometry service are commented out but not deleted. +The current plan is to replace those with a `phlex` provider that provides the geometry +data as data products and these providers will be scheduled by the flow graph. +In the main algorithm the geometry is used to get the plane number for each +wire. For now, I just hard coded the plane number to 0. The signal type in +the `recob::Hit` objects is also from the geometry and for now always set +to 0. Similarly, the `WireID` in the `recob::Hit` objects is from the geometry +and for now set to a default value. `FillOutHitParameterVector` depends +on the geometry also, but it only checks the size of a configuration +vector against the number of planes. If this is skipped but the +configuration vector is sized to match the number of planes, then things +should still work ok. For now, I just skip FillOutHitParameterVector. + +I copied in the content of the `TMath::Gaus` function from ROOT to +avoid depending on ROOT in `phlex` code for now (slightly edited). +Are we allowed to depend on ROOT in `phlex` algorithm code? That could +easily be restored if the build system allows it. + +# Testing + +This was tested by running `GausHitFinder` with the `art` framework and running the modified `GausHitFinder` with the `phlex` framework. The two processes used identical input (the `std::vector`). The output `std::vector` objects were compared using a text file printed during each process. The output was identical except for the following: + +Two output data members of `Hit` were ignored because they depend on the `Geometry` and that is not implemented yet for `phlex`. In the `phlex` version, these data members are filled with default values. +``` + geo::SigType_t fSignalType; ///< signal type for the plane of the hit + geo::WireID fWireID; ///< WireID for the hit (Cryostat, TPC, Plane, Wire) +``` + +The `phlex` process was run with multithreading and that causes the order of `Hit` objects to vary from one execution to the next and also the order of events to vary. In the comparison the `Hit` objects were sorted and we had to be careful to compare matching events. + +Contact the `phlex` group if you are interested in details related +to running `art` side of this test or want the input files (or +to regenerate them). diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.cxx new file mode 100644 index 0000000..7738bf6 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.cxx @@ -0,0 +1,124 @@ +#include +#include + +#include "CandHitStandard.h" + +namespace examples { + + CandHitStandard::CandHitStandard(const CandHitStandardCfg& cfg) + : fRoiThreshold(cfg.fRoiThreshold) {} + + void CandHitStandard::findHitCandidates( + const recob::Wire::RegionsOfInterest_t::datarange_t& dataRange, + const size_t roiStartTick, + const size_t channel, + HitCandidateVec& hitCandidateVec) const + { + // Recover the actual waveform + const Waveform& waveform = dataRange.data(); + + // Use the recursive version to find the candidate hits + findHitCandidates(waveform.begin(), waveform.end(), roiStartTick, hitCandidateVec); + } + + void CandHitStandard::findHitCandidates(std::vector::const_iterator startItr, + std::vector::const_iterator stopItr, + const size_t roiStartTick, + HitCandidateVec& hitCandidateVec) const + { + // Need a minimum number of ticks to do any work here + if (std::distance(startItr, stopItr) > 4) { + // Find the highest peak in the range given + auto maxItr = std::max_element(startItr, stopItr); + + float maxValue = *maxItr; + int maxTime = std::distance(startItr, maxItr); + + if (maxValue > fRoiThreshold) { + // backwards to find first bin for this candidate hit + auto firstItr = std::distance(startItr, maxItr) > 2 ? maxItr - 1 : startItr; + + while (firstItr != startItr) { + // Check for pathology where waveform goes too negative + if (*firstItr < -fRoiThreshold) break; + + // Check both sides of firstItr and look for min/inflection point + if (*firstItr < *(firstItr + 1) && *firstItr <= *(firstItr - 1)) break; + + firstItr--; + } + + int firstTime = std::distance(startItr, firstItr); + + // Recursive call to find all candidate hits earlier than this peak + findHitCandidates(startItr, firstItr + 1, roiStartTick, hitCandidateVec); + + // forwards to find last bin for this candidate hit + auto lastItr = std::distance(maxItr, stopItr) > 2 ? maxItr + 1 : stopItr - 1; + + while (lastItr != stopItr - 1) { + // Check for pathology where waveform goes too negative + if (*lastItr < -fRoiThreshold) break; + + // Check both sides of firstItr and look for min/inflection point + if (*lastItr <= *(lastItr + 1) && *lastItr < *(lastItr - 1)) break; + + lastItr++; + } + + int lastTime = std::distance(startItr, lastItr); + + // Now save this candidate's start and max time info + HitCandidate hitCandidate; + hitCandidate.startTick = roiStartTick + firstTime; + hitCandidate.stopTick = roiStartTick + lastTime; + hitCandidate.maxTick = roiStartTick + firstTime; + hitCandidate.minTick = roiStartTick + lastTime; + hitCandidate.maxDerivative = *(startItr + firstTime); + hitCandidate.minDerivative = *(startItr + lastTime); + hitCandidate.hitCenter = roiStartTick + maxTime; + hitCandidate.hitSigma = std::max(2., float(lastTime - firstTime) / 6.); + hitCandidate.hitHeight = maxValue; + + hitCandidateVec.push_back(hitCandidate); + + // Recursive call to find all candidate hits later than this peak + findHitCandidates(lastItr + 1, + stopItr, + roiStartTick + std::distance(startItr, lastItr + 1), + hitCandidateVec); + } + } + } + + void CandHitStandard::MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&, + const HitCandidateVec& hitCandidateVec, + MergeHitCandidateVec& mergedHitsVec) const + { + // If no hits then nothing to do here + if (hitCandidateVec.empty()) return; + + // The idea is to group hits that "touch" so they can be part of common fit, those that + // don't "touch" are fit independently. So here we build the output vector to achieve that + HitCandidateVec groupedHitVec; + int lastTick = hitCandidateVec.front().stopTick; + + // Step through the input hit candidates and group them by proximity + for (const auto& hitCandidate : hitCandidateVec) { + // Check condition that we have a new grouping + if (int(hitCandidate.startTick) - lastTick > 1) { + mergedHitsVec.emplace_back(groupedHitVec); + + groupedHitVec.clear(); + } + + // Add the current hit to the current group + groupedHitVec.emplace_back(hitCandidate); + + lastTick = hitCandidate.stopTick; + } + + // Check end condition + if (!groupedHitVec.empty()) mergedHitsVec.emplace_back(groupedHitVec); + } +} diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.h new file mode 100644 index 0000000..683528b --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/CandHitStandard.h @@ -0,0 +1,44 @@ +#ifndef CANDHITSTANDARD_H +#define CANDHITSTANDARD_H + +//////////////////////////////////////////////////////////////////////// +/// \file CandHitStandard.cc +/// \author T. Usher +//////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "ICandidateHitFinder.h" +#include "Wire.h" + +namespace examples { + + struct CandHitStandardCfg { + float fRoiThreshold; + }; + + class CandHitStandard : public ICandidateHitFinder { + public: + explicit CandHitStandard(const CandHitStandardCfg& cfg); + + void findHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&, + const size_t, + const size_t, + HitCandidateVec&) const override; + + void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&, + const HitCandidateVec&, + MergeHitCandidateVec&) const override; + + private: + void findHitCandidates(std::vector::const_iterator, + std::vector::const_iterator, + const size_t, + HitCandidateVec&) const; + + const float fRoiThreshold; ///< minimum maximum to minimum peak distance + }; + +} // namespace examples +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.cxx new file mode 100644 index 0000000..6e6f9c5 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.cxx @@ -0,0 +1,120 @@ +/** **************************************************************************** + * @file lardataobj/RecoBase/Hit.cxx + * @brief Definition of signal hit object. + * @author mitchell.soderberg@yale.edu + * @see lardataobj/RecoBase/Hit.h + * + * ****************************************************************************/ + +// Hit header +#include "Hit.h" + +// C/C++ standard library +#include +#include + +namespace recob { + + //---------------------------------------------------------------------- + Hit::Hit() + : fChannel(raw::InvalidChannelID) + , fStartTick(0) + , fEndTick(0) + , fPeakTime(0.) + , fSigmaPeakTime(-1.) + , fRMS(0.) + , fPeakAmplitude(0.) + , fSigmaPeakAmplitude(-1.) + , fROISummedADC(0.) + , fHitSummedADC(0.) + , fIntegral(0.) + , fSigmaIntegral(-1.) + , fMultiplicity(0) + , fLocalIndex(-1) + , fGoodnessOfFit(0.) + , fNDF(-1) + , fView(geo::kUnknown) + , fSignalType(geo::kMysteryType) + , fWireID() // invalid + {} + + //---------------------------------------------------------------------- + Hit::Hit(raw::ChannelID_t channel, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float peak_time, + float sigma_peak_time, + float rms, + float peak_amplitude, + float sigma_peak_amplitude, + float ROIsummedADC, + float HitsummedADC, + float hit_integral, + float hit_sigma_integral, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + geo::View_t view, + geo::SigType_t signal_type, + geo::WireID wireID) + : fChannel(channel) + , fStartTick(start_tick) + , fEndTick(end_tick) + , fPeakTime(peak_time) + , fSigmaPeakTime(sigma_peak_time) + , fRMS(rms) + , fPeakAmplitude(peak_amplitude) + , fSigmaPeakAmplitude(sigma_peak_amplitude) + , fROISummedADC(ROIsummedADC) + , fHitSummedADC(HitsummedADC) + , fIntegral(hit_integral) + , fSigmaIntegral(hit_sigma_integral) + , fMultiplicity(multiplicity) + , fLocalIndex(local_index) + , fGoodnessOfFit(goodness_of_fit) + , fNDF(dof) + , fView(view) + , fSignalType(signal_type) + , fWireID(wireID) + {} + + //---------------------------------------------------------------------- + // ostream operator. + // + std::ostream& operator<<(std::ostream& o, Hit const& hit) + { + o << std::setiosflags(std::ios::fixed) << std::setprecision(2); + o << " Channel " << std::setw(5) << std::right << hit.Channel() << " View = " << std::setw(3) + << std::right << hit.View() << " Signal type = " << std::setw(3) << std::right + << hit.SignalType() << " Wire = " << std::setw(3) << std::right << hit.WireID() + << "\n\tStartTick = " << std::setw(7) << std::right << hit.StartTick() + << "\tEndTick = " << std::setw(7) << std::right << hit.EndTick() + << "\tPeakTime = " << std::setw(7) << std::right << hit.PeakTime() << " +/- " << std::setw(7) + << std::right << hit.SigmaPeakTime() << "\tRMS = " << std::setw(7) << std::right << hit.RMS() + << "\n\tAmplitude = " << std::setw(7) << std::right << hit.PeakAmplitude() << " +/- " + << std::setw(7) << std::right << hit.SigmaPeakAmplitude() << "\tIntegral = " << std::setw(7) + << std::right << hit.Integral() << " +/- " << std::setw(7) << std::right + << hit.SigmaIntegral() << "\tROIADCsum = " << std::setw(7) << std::right << hit.ROISummedADC() + << "\tHitADCsum = " << std::setw(7) << std::right << hit.HitSummedADC() + << "\tMultiplicity = " << std::setw(5) << std::right << hit.LocalIndex() << " of " + << hit.Multiplicity() << "\tGoodnessOfFit = " << std::setw(7) << std::right + << hit.GoodnessOfFit() << " DoF = " << std::setw(7) << std::right << hit.DegreesOfFreedom() + << std::endl; + return o; + } // operator<< (std::ostream, Hit) + + //---------------------------------------------------------------------- + // < operator. + // + bool operator<(const Hit& a, const Hit& b) + { + if (a.Channel() != b.Channel()) return a.Channel() < b.Channel(); + if (a.View() != b.View()) return a.View() < b.View(); + if (a.PeakTime() != b.PeakTime()) return a.PeakTime() < b.PeakTime(); + + return false; //They are equal + } // operator< (Hit, Hit) + + //---------------------------------------------------------------------- +} // namespace recob diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.h new file mode 100644 index 0000000..aea4102 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/Hit.h @@ -0,0 +1,310 @@ +/** **************************************************************************** + * @file lardataobj/RecoBase/Hit.h + * @brief Declaration of signal hit object. + * @author mitchell.soderberg@yale.edu + * @see lardataobj/RecoBase/Hit.cxx + * + * Changes: + * 20141212 Gianluca Petrillo (petrillo@fnal.gov) + * data architecture revision changes (v13 -> v14): + * - fRawDigit and RawDigit() removed + * - fWire and Wire() removed + * - fSignalType and SignalType() removed + * - fChannel and Channel() added + * - constructors now take directly a RawDigit, not its art::Ptr + * 20150129 Gianluca Petrillo (petrillo@fnal.gov) + * data architecture revision changes (v14 -> v15): + * - removed fHitSignal + * + * ****************************************************************************/ + +#ifndef LARDATAOBJ_RECOBASE_HIT_H +#define LARDATAOBJ_RECOBASE_HIT_H + +// C/C++ standard librraies +#include + +// LArSoft libraries +#include "RawTypes.h" // raw::ChannelID_t +#include "geo_types.h" // geo::View_t, geo::SignalType, geo::WireID + +namespace recob { + + /** + * @brief 2D representation of charge deposited in the TDC/wire plane + * + * Hits are assumed to be made from deconvoluted, unipolar, calibrated + * signals. + * They identify a charge deposit in a specific location and time; + * the location is absolute and unique in the detector, while the time is + * relative to the start of sampling (tick time). + * + * The version 14 of recob::Hit introduces the following changes: + * - StartTime() becomes StartTick(), StopTime() becomes StopTick() + * - Charge(true) is now PeakAmplitude(), Charge(false) (default) is Integral() + */ + class Hit { + public: + /// Default constructor: a hit with no signal + Hit(); + + private: + raw::ChannelID_t fChannel; ///< ID of the readout channel the hit was extracted from + raw::TDCtick_t fStartTick; ///< initial tdc tick for hit + raw::TDCtick_t fEndTick; ///< final tdc tick for hit + float fPeakTime; ///< time of the signal peak, in tick units + float fSigmaPeakTime; ///< uncertainty for the signal peak, in tick units + float fRMS; ///< RMS of the hit shape, in tick units + float fPeakAmplitude; ///< the estimated amplitude of the hit at its peak, in ADC units + float + fSigmaPeakAmplitude; ///< uncertainty on estimated amplitude of the hit at its peak, in ADC units + float fROISummedADC; ///< the sum of calibrated ADC counts of the ROI + float fHitSummedADC; ///< the sum of calibrated ADC counts of the hit + float + fIntegral; ///< the integral under the calibrated signal waveform of the hit, in tick x ADC units + float + fSigmaIntegral; ///< the uncertainty of integral under the calibrated signal waveform of the hit, in ADC units + short int fMultiplicity; ///< how many hits could this one be shared with + short int fLocalIndex; ///< index of this hit among the Multiplicity() hits in the signal window + float fGoodnessOfFit; ///< how well do we believe we know this hit? + int fNDF; ///< degrees of freedom in the determination of the hit shape + geo::View_t fView; ///< view for the plane of the hit + geo::SigType_t fSignalType; ///< signal type for the plane of the hit + geo::WireID fWireID; ///< WireID for the hit (Cryostat, TPC, Plane, Wire) + + friend class HitCreator; // helper to create hits + + public: + /** + * @brief Constructor: directly sets all the fields + * @param channel ID of the readout channel the hit was extracted from + * @param start_tick initial tdc tick for hit + * @param end_tick final tdc tick for hit + * @param peak_time tdc for the peak charge deposition + * @param sigma_peak_time uncertainty for tdc of the peak + * @param rms RMS of the hit shape + * @param peak_amplitude the estimated amplitude of the hit at its peak + * @param sigma_peak_amplitude the uncertainty on the estimated amplitude of the hit at its peak + * @param ROIsummedADC the sum of calibrated ADC counts of the ROI + * @param HitsummedADC the sum of calibrated ADC counts of the hit + * @param hit_integral the integral under the calibrated signal waveform of the hit + * @param hit_sigma_integral uncertainty on the integral under the calibrated signal waveform of the hit + * @param multiplicity how many hits could this one be shared with + * @param local_index index of this hit among the Multiplicity() hits in the signal window + * @param goodness_of_fit how well do we believe we know this hit? + * @param dof number of degrees of freedom in the determination of hit shape + * @param view view for the plane of the hit + * @param signal_type signal type for the plane of the hit + * @param wireID WireID for the hit (Cryostat TPC Plane Wire) + * + * The tick parameters are real numbers, since they can in principle come + * from some processing. + */ + Hit(raw::ChannelID_t channel, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float peak_time, + float sigma_peak_time, + float rms, + float peak_amplitude, + float sigma_peak_amplitude, + float ROIsummedADC, + float HitsummedADC, + float hit_integral, + float hit_sigma_integral, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + geo::View_t view, + geo::SigType_t signal_type, + geo::WireID wireID); + + /// @{ + /// @name Accessors + + /// Initial tdc tick for hit + raw::TDCtick_t StartTick() const; + + /// Final tdc tick for hit + raw::TDCtick_t EndTick() const; + + /// Time of the signal peak, in tick units + float PeakTime() const; + + /// Uncertainty for the signal peak, in tick units + float SigmaPeakTime() const; + + /// RMS of the hit shape, in tick units + float RMS() const; + + /// The estimated amplitude of the hit at its peak, in ADC units + float PeakAmplitude() const; + + /// Uncertainty on estimated amplitude of the hit at its peak, in ADC units + float SigmaPeakAmplitude() const; + + /// The sum of calibrated ADC counts of the ROI (0. by default) + float ROISummedADC() const; + + /// The sum of calibrated ADC counts of the hit (0. by default) + float HitSummedADC() const; + + /// Integral under the calibrated signal waveform of the hit, in tick x ADC units + float Integral() const; + + ///< Uncertainty of integral under the calibrated signal waveform of the hit, in ADC units + float SigmaIntegral() const; + + /// How many hits could this one be shared with + short int Multiplicity() const; + + ///< Index of this hit among the Multiplicity() hits in the signal window (-1 by default) + short int LocalIndex() const; + + ///< How well do we believe we know this hit? + float GoodnessOfFit() const; + + ///< Degrees of freedom in the determination of the hit signal shape (-1 by default) + int DegreesOfFreedom() const; + + /// ID of the readout channel the hit was extracted from + raw::ChannelID_t Channel() const; + + /// View for the plane of the hit + geo::View_t View() const; + + /// Signal type for the plane of the hit + geo::SigType_t SignalType() const; + + ///< ID of the wire the hit is on (Cryostat, TPC, Plane, Wire) + geo::WireID const& WireID() const; + + /// @} + + //@{ + /** + * @brief Returns a time sigmas RMS away from the peak time + * @param sigmas the number of RMS units to move away + * @return the shifted time in TDC ticks + * + * PeakTimePlusRMS() returns PeakTime() + sigmas x RMS(); + * PeakTimeMinusRMS() returns PeakTime() - sigmas x RMS(). + * + * @note StartTime() of recob::Hit version <=13 was defined by + * GausHitFinder to be PeakTimePlusRMS(-1.), and EndTime() was + * PeakTimePlusRMS(+1.). + */ + float PeakTimePlusRMS(float sigmas = +1.) const; + float PeakTimeMinusRMS(float sigmas = +1.) const; + //@} + + /** + * @brief Returns the distance of the specified time from peak, in RMS units + * @param time the time, in TDC tick units + * @return the distance of the specified time from peak, in RMS units + * + * This returns (ticks - PeakTime()) / RMS(). + * There is no protection in case RMS is 0! + */ + float TimeDistanceAsRMS(float time) const; + + friend std::ostream& operator<<(std::ostream& o, const Hit& a); + friend bool operator<(const Hit& a, const Hit& b); + + }; // class Hit +} // namespace recob + +inline raw::TDCtick_t recob::Hit::StartTick() const +{ + return fStartTick; +} +inline raw::TDCtick_t recob::Hit::EndTick() const +{ + return fEndTick; +} +inline float recob::Hit::PeakTime() const +{ + return fPeakTime; +} +inline float recob::Hit::SigmaPeakTime() const +{ + return fSigmaPeakTime; +} +inline float recob::Hit::RMS() const +{ + return fRMS; +} +inline float recob::Hit::PeakAmplitude() const +{ + return fPeakAmplitude; +} +inline float recob::Hit::SigmaPeakAmplitude() const +{ + return fSigmaPeakAmplitude; +} +inline float recob::Hit::ROISummedADC() const +{ + return fROISummedADC; +} +inline float recob::Hit::HitSummedADC() const +{ + return fHitSummedADC; +} +inline float recob::Hit::Integral() const +{ + return fIntegral; +} +inline float recob::Hit::SigmaIntegral() const +{ + return fSigmaIntegral; +} +inline short int recob::Hit::Multiplicity() const +{ + return fMultiplicity; +} +inline short int recob::Hit::LocalIndex() const +{ + return fLocalIndex; +} +inline float recob::Hit::GoodnessOfFit() const +{ + return fGoodnessOfFit; +} +inline int recob::Hit::DegreesOfFreedom() const +{ + return fNDF; +} +inline raw::ChannelID_t recob::Hit::Channel() const +{ + return fChannel; +} +inline geo::SigType_t recob::Hit::SignalType() const +{ + return fSignalType; +} +inline geo::View_t recob::Hit::View() const +{ + return fView; +} +inline geo::WireID const& recob::Hit::WireID() const +{ + return fWireID; +} + +inline float recob::Hit::PeakTimePlusRMS(float sigmas /* = +1. */) const +{ + return PeakTime() + sigmas * RMS(); +} + +inline float recob::Hit::PeakTimeMinusRMS(float sigmas /* = +1. */) const +{ + return PeakTimePlusRMS(-sigmas); +} + +inline float recob::Hit::TimeDistanceAsRMS(float time) const +{ + return (time - PeakTime()) / RMS(); +} + +#endif // LARDATAOBJ_RECOBASE_HIT_H diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.cxx new file mode 100644 index 0000000..4d30d8a --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.cxx @@ -0,0 +1,33 @@ +#include "HitFilterAlg.h" + +#include "geo_types.h" +#include "Hit.h" + +//#include "messagefacility/MessageLogger/MessageLogger.h" + +namespace examples { + + HitFilterAlg::HitFilterAlg(HitFilterAlgCfg const& cfg) + : fMinPulseHeight(cfg.fMinPulseHeight) + , fMinPulseSigma(cfg.fMinPulseSigma) {} + + bool HitFilterAlg::IsGoodHit(const recob::Hit& hit) const + { + + const float hitPH = hit.PeakAmplitude(); + const float hitSigma = hit.RMS(); + + const geo::WireID& wireID = hit.WireID(); + const size_t view = wireID.Plane; + + if (view >= fMinPulseSigma.size() || view >= fMinPulseHeight.size()) { + // mf::LogError("HitFilterAlg") << "Filtering settings not configured for all views! Will not " + // "filter hits in unconfigured views!"; + return true; + } + + if (hitPH > fMinPulseHeight[view] && hitSigma > fMinPulseSigma[view]) { return true; } + else + return false; + } +} //end namespace examples diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.h new file mode 100644 index 0000000..be9c045 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/HitFilterAlg.h @@ -0,0 +1,38 @@ +//////////////////////////////////////////////////////////////////////// +// Class: HitFilterAlg +// Purpose: Provide configurable hit filtering functions for identifying +// noise and other undesireable hits +// +// Original code by Tracy Usher, converted to a larsoft algorithm by Brandon Eberly +//////////////////////////////////////////////////////////////////////// + +#ifndef HITFILTERALG_H +#define HITFILTERALG_H + +#include + +namespace recob { + class Hit; +} + +namespace examples { + + struct HitFilterAlgCfg { + std::vector fMinPulseHeight; + std::vector fMinPulseSigma; + }; + + class HitFilterAlg { + public: + explicit HitFilterAlg(HitFilterAlgCfg const& cfg); + virtual ~HitFilterAlg() {} + + bool IsGoodHit(const recob::Hit& hit) const; + + private: + const std::vector fMinPulseHeight; + const std::vector fMinPulseSigma; + }; + +} //end namespace examples +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/ICandidateHitFinder.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/ICandidateHitFinder.h new file mode 100644 index 0000000..61d79a6 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/ICandidateHitFinder.h @@ -0,0 +1,56 @@ +/////////////////////////////////////////////////////////////////////// +/// +/// \file ICandidateHitFinder.h +/// +/// \brief This provides an interface for tools which are tasked with +/// finding candidate hits on input waveforms +/// +/// \author T. Usher +/// +//////////////////////////////////////////////////////////////////////// + +#ifndef ICandidateHitFinder_H +#define ICandidateHitFinder_H + +#include +#include + +#include "Wire.h" + +namespace examples { + class ICandidateHitFinder { + public: + virtual ~ICandidateHitFinder() noexcept = default; + + // Define a structure to contain hits + struct HitCandidate { + size_t startTick; + size_t stopTick; + size_t maxTick; + size_t minTick; + float maxDerivative; + float minDerivative; + float hitCenter; + float hitSigma; + float hitHeight; + }; + + using HitCandidateVec = std::vector; + using MergeHitCandidateVec = std::vector; + + using Waveform = std::vector; + + // Search for candidate hits on the input waveform + virtual void findHitCandidates( + const recob::Wire::RegionsOfInterest_t::datarange_t&, // Waveform (with range info) to analyze + const size_t, // waveform start tick + const size_t, // channel # + HitCandidateVec&) const = 0; // output candidate hits + + virtual void MergeHitCandidates(const recob::Wire::RegionsOfInterest_t::datarange_t&, + const HitCandidateVec&, + MergeHitCandidateVec&) const = 0; + }; +} + +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/IPeakFitter.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/IPeakFitter.h new file mode 100644 index 0000000..828a020 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/IPeakFitter.h @@ -0,0 +1,45 @@ +/////////////////////////////////////////////////////////////////////// +/// +/// \file IPeakFitter.h +/// +/// \brief This provides an interface for tools which are tasked with +/// fitting peaks on input waveforms +/// +/// \author T. Usher +/// +//////////////////////////////////////////////////////////////////////// + +#ifndef IPeakFitter_H +#define IPeakFitter_H + +#include "ICandidateHitFinder.h" + +#include + +namespace examples { + class IPeakFitter { + public: + // Define standard art tool interface + + // Define a structure to contain hits + struct PeakFitParams_t { + float peakCenter; + float peakCenterError; + float peakSigma; + float peakSigmaError; + float peakAmplitude; + float peakAmplitudeError; + }; + + using PeakParamsVec = std::vector; + virtual ~IPeakFitter() = default; + // Get parameters for input candidate peaks + virtual void findPeakParameters(const std::vector&, + const ICandidateHitFinder::HitCandidateVec&, + PeakParamsVec&, + double&, + int&) const = 0; + }; +} + +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.cxx new file mode 100644 index 0000000..aee277a --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.cxx @@ -0,0 +1,393 @@ +#include "MarqFitAlg.h" +#include + +namespace gshf { + + /* multi-Gaussian function, number of Gaussians is npar divided by 3 */ + void MarqFitAlg::fgauss(const float yd[], + const float p[], + const int npar, + const int ndat, + std::vector& res) + { +#if defined WITH_OPENMP +#pragma omp simd +#endif + for (int i = 0; i < ndat; i++) { + float yf = 0.; + for (int j = 0; j < npar; j += 3) { + yf = yf + p[j] * std::exp(-0.5 * std::pow((float(i) - p[j + 1]) / p[j + 2], 2)); + } + res[i] = yd[i] - yf; + } + } + + /* analytic derivatives for multi-Gaussian function in fgauss */ + void MarqFitAlg::dgauss(const float p[], const int npar, const int ndat, std::vector& dydp) + { +#if defined WITH_OPENMP +#pragma omp simd +#ifdef __INTEL_COMPILER +#pragma ivdep +#endif +#endif + for (int i = 0; i < ndat; i++) { + for (int j = 0; j < npar; j += 3) { + const float xmu = float(i) - p[j + 1]; + const float xmu_sg = xmu / p[j + 2]; + const float xmu_sg2 = xmu_sg * xmu_sg; + dydp[i * npar + j] = std::exp(-0.5 * xmu_sg2); + dydp[i * npar + j + 1] = p[j] * dydp[i * npar + j] * xmu_sg / p[j + 2]; + dydp[i * npar + j + 2] = dydp[i * npar + j + 1] * xmu_sg; + } + } + } + + /* calculate ChiSquared */ + float MarqFitAlg::cal_xi2(const std::vector& res, const int ndat) + { + int i; + float xi2; + xi2 = 0.; + for (i = 0; i < ndat; i++) { + xi2 += res[i] * res[i]; + } + return xi2; + } + + /* setup the beta and (curvature) matrices */ + void MarqFitAlg::setup_matrix(const std::vector& res, + const std::vector& dydp, + const int npar, + const int ndat, + std::vector& beta, + std::vector& alpha) + { + int i, j, k; + +/* ... Calculate beta */ +#if defined WITH_OPENMP +#pragma omp simd +#ifdef __INTEL_COMPILER +#pragma ivdep +#endif +#endif + for (j = 0; j < npar; j++) { + beta[j] = 0.0; + for (i = 0; i < ndat; i++) { + beta[j] += res[i] * dydp[i * npar + j]; + } + } + +/* ... Calculate alpha */ +#if defined WITH_OPENMP +#pragma omp simd +#ifdef __INTEL_COMPILER +#pragma ivdep +#endif +#endif + for (j = 0; j < npar; j++) { + for (k = j; k < npar; k++) { + alpha[j * npar + k] = 0.0; + for (i = 0; i < ndat; i++) { + alpha[j * npar + k] += dydp[i * npar + j] * dydp[i * npar + k]; + } + if (k != j) alpha[k * npar + j] = alpha[j * npar + k]; + } + } + } + + /* solve system of linear equations */ + void MarqFitAlg::solve_matrix(const std::vector& beta, + const std::vector& alpha, + const int npar, + std::vector& dp) + { + int i, j, k, imax; + float hmax, hsav; + + std::vector> h(npar, std::vector(npar + 1, 0)); + + /* ... set up augmented N x N+1 matrix */ + for (i = 0; i < npar; i++) { + h[i][npar] = beta[i]; + for (j = 0; j < npar; j++) { + h[i][j] = alpha[i * npar + j]; + } + } + + /* ... diagonalize N x N matrix but do only terms required for solution */ + for (i = 0; i < npar; i++) { + hmax = h[i][i]; + imax = i; + for (j = i + 1; j < npar; j++) { + if (h[j][i] > hmax) { + hmax = h[j][i]; + imax = j; + } + } + if (imax != i) { + for (k = 0; k <= npar; k++) { + hsav = h[i][k]; + h[i][k] = h[imax][k]; + h[imax][k] = hsav; + } + } + for (j = 0; j < npar; j++) { + if (j == i) continue; + for (k = i; k < npar; k++) { + h[j][k + 1] -= h[i][k + 1] * h[j][i] / h[i][i]; + } + } + } + /* ... scale (N+1)'th column with factor which normalizes the diagonal */ + for (i = 0; i < npar; i++) { + dp[i] = h[i][npar] / h[i][i]; + } + } + + float MarqFitAlg::invrt_matrix(std::vector& alphaf, const int npar) + { + /* + Inverts the curvature matrix alpha using Gauss-Jordan elimination and + returns the determinant. This is based on the implementation in "Data + Reduction and Error Analysis for the Physical Sciences" by P. R. Bevington. + That implementation, in turn, uses the algorithm of the subroutine MINV, + described on page 44 of the IBM System/360 Scientific Subroutine Package + Program Description Manual (GH20-0586-0). This only needs to be called + once after the fit has converged if the parameter errors are desired. + */ + + //turn input alphas into doubles + std::vector alpha(npar * npar); + + int i, j, k; + std::vector ik(npar); + std::vector jk(npar); + double aMax, save, det; + float detf; + + for (i = 0; i < npar * npar; i++) { + alpha[i] = alphaf[i]; + } + + det = 0; + /* ... search for the largest element which we will then put in the diagonal */ + for (k = 0; k < npar; k++) { + aMax = 0; + for (i = k; i < npar; i++) { + for (j = k; j < npar; j++) { + + // alpha[i*npar+j]=alphaf[i*npar+j]; + + if (fabs(alpha[i * npar + j]) > fabs(aMax)) { + aMax = alpha[i * npar + j]; + ik[k] = i; + jk[k] = j; + } + } + } + if (aMax == 0) return (det); /* return 0 determinant to signal problem */ + det = 1; + /* ... interchange rows if necessary to put aMax in diag */ + i = ik[k]; + if (i > k) { + for (j = 0; j < npar; j++) { + save = alpha[k * npar + j]; + alpha[k * npar + j] = alpha[i * npar + j]; + alpha[i * npar + j] = -save; + } + } + /* ... interchange columns if necessary to put aMax in diag */ + j = jk[k]; + if (j > k) { + for (i = 0; i < npar; i++) { + save = alpha[i * npar + k]; + alpha[i * npar + k] = alpha[i * npar + j]; + alpha[i * npar + j] = -save; + } + } + /* ... accumulate elements of inverse matrix */ + for (i = 0; i < npar; i++) { + if (i != k) alpha[i * npar + k] = -alpha[i * npar + k] / aMax; + } +#if defined WITH_OPENMP +#pragma omp simd +#ifdef __INTEL_COMPILER +#pragma ivdep +#endif +#endif + for (i = 0; i < npar; i++) { + for (j = 0; j < npar; j++) { + if ((i != k) && (j != k)) + alpha[i * npar + j] = alpha[i * npar + j] + alpha[i * npar + k] * alpha[k * npar + j]; + } + } + for (j = 0; j < npar; j++) { + if (j != k) alpha[k * npar + j] = alpha[k * npar + j] / aMax; + } + alpha[k * npar + k] = 1 / aMax; + det = det * aMax; + } + + /* ... restore ordering of matrix */ + for (k = npar - 1; k >= 0; k--) { + j = ik[k]; + if (j > k) { + for (i = 0; i < npar; i++) { + save = alpha[i * npar + k]; + alpha[i * npar + k] = -alpha[i * npar + j]; + alpha[i * npar + j] = save; + } + } + i = jk[k]; + if (i > k) { + for (j = 0; j < npar; j++) { + save = alpha[k * npar + j]; + alpha[k * npar + j] = -alpha[i * npar + j]; + alpha[i * npar + j] = save; + } + } + } + + for (i = 0; i < npar * npar; i++) { + alphaf[i] = alpha[i]; + } + + detf = det; + return (detf); + } + + /* Calculate parameter errors */ + int MarqFitAlg::cal_perr(float p[], float y[], const int nParam, const int nData, float perr[]) + { + int i, j; + float det; + + std::vector res(nData); + std::vector dydp(nData * nParam); + std::vector beta(nParam); + std::vector alpha(nParam * nParam); + std::vector> alpsav(nParam, std::vector(nParam)); + + fgauss(y, p, nParam, nData, res); + dgauss(p, nParam, nData, dydp); + setup_matrix(res, dydp, nParam, nData, beta, alpha); + for (i = 0; i < nParam; i++) { + for (j = 0; j < nParam; j++) { + alpsav[i][j] = alpha[i * nParam + j]; + } + } + det = invrt_matrix(alpha, nParam); + + if (det == 0) return 1; + for (i = 0; i < nParam; i++) { + if (alpha[i * nParam + i] >= 0.) { perr[i] = sqrt(alpha[i * nParam + i]); } + else { + perr[i] = alpha[i * nParam + i]; + } + } + + return 0; + } + + int MarqFitAlg::mrqdtfit(float& lambda, + float p[], + float y[], + const int nParam, + const int nData, + float& chiSqr, + float& dchiSqr) + { + std::vector plimmin(nParam, std::numeric_limits::lowest()); + std::vector plimmax(nParam, std::numeric_limits::max()); + return mrqdtfit(lambda, p, &plimmin[0], &plimmax[0], y, nParam, nData, chiSqr, dchiSqr); + } + + int MarqFitAlg::mrqdtfit(float& lambda, + float p[], + float plimmin[], + float plimmax[], + float y[], + const int nParam, + const int nData, + float& chiSqr, + float& dchiSqr) + { + int j; + float nu, rho, lzmlh, amax, chiSq0; + + std::vector res(nData); + std::vector beta(nParam); + std::vector dp(nParam); + std::vector alpsav(nParam); + std::vector psav(nParam); + std::vector dydp(nData * nParam); + std::vector alpha(nParam * nParam); + + bool haslimits = false; + for (j = 0; j < nParam; j++) { + if (plimmin[j] > std::numeric_limits::lowest() || + plimmax[j] < std::numeric_limits::max()) { + haslimits = true; + break; + } + } + + fgauss(y, p, nParam, nData, res); + chiSq0 = cal_xi2(res, nData); + dgauss(p, nParam, nData, dydp); + setup_matrix(res, dydp, nParam, nData, beta, alpha); + if (lambda < 0.) { + amax = -999.; + for (j = 0; j < nParam; j++) { + if (alpha[j * nParam + j] > amax) amax = alpha[j * nParam + j]; + } + lambda = 0.001 * amax; + } + for (j = 0; j < nParam; j++) { + alpsav[j] = alpha[j * nParam + j]; + alpha[j * nParam + j] = alpsav[j] + lambda; + } + solve_matrix(beta, alpha, nParam, dp); + + nu = 2.; + rho = -1.; + + do { + for (j = 0; j < nParam; j++) { + psav[j] = p[j]; + p[j] = p[j] + dp[j]; + } + fgauss(y, p, nParam, nData, res); + chiSqr = cal_xi2(res, nData); + if (haslimits) { + for (j = 0; j < nParam; j++) { + if (p[j] <= plimmin[j] || p[j] >= plimmax[j]) + chiSqr *= 10000; //penalty for going out of limits! + } + } + + lzmlh = 0.; + for (j = 0; j < nParam; j++) { + lzmlh += dp[j] * (lambda * dp[j] + beta[j]); + } + rho = 2. * (chiSq0 - chiSqr) / lzmlh; + if (rho < 0.) { + for (j = 0; j < nParam; j++) + p[j] = psav[j]; + chiSqr = chiSq0; + lambda = nu * lambda; + nu = 2. * nu; + for (j = 0; j < nParam; j++) { + alpha[j * nParam + j] = alpsav[j] + lambda; + } + solve_matrix(beta, alpha, nParam, dp); + } + } while (rho < 0.); + lambda = lambda * fmax(0.333333, 1. - pow(2. * rho - 1., 3)); + dchiSqr = chiSqr - chiSq0; + return 0; + } + +} //end namespace gshf diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.h new file mode 100644 index 0000000..7b825f8 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/MarqFitAlg.h @@ -0,0 +1,67 @@ +//////////////////////////////////////////////////////////////////////// +// Class: MarqFitAlg +// Purpose: Fit gaussians +// +// +// Original code by Mike Wang, converted to a larsoft algorithm by S. Berkman +//////////////////////////////////////////////////////////////////////// + +#ifndef MARQFITALG_H +#define MARQFITALG_H + +#include +#include +#include + +//#include "fhiclcpp/ParameterSet.h" +//#include "lardataobj/RecoBase/Hit.h" + +namespace gshf { + + class MarqFitAlg { + public: + explicit MarqFitAlg(); + virtual ~MarqFitAlg() {} + + int cal_perr(float p[], float y[], const int nParam, const int nData, float perr[]); + int mrqdtfit(float& lambda, + float p[], + float y[], + const int nParam, + const int nData, + float& chiSqr, + float& dchiSqr); + int mrqdtfit(float& lambda, + float p[], + float plimmin[], + float plimmax[], + float y[], + const int nParam, + const int nData, + float& chiSqr, + float& dchiSqr); + + private: + //these functions are called by the public functions + void fgauss(const float yd[], + const float p[], + const int npar, + const int ndat, + std::vector& res); + void dgauss(const float p[], const int npar, const int ndat, std::vector& dydp); + float cal_xi2(const std::vector& res, const int ndat); + void setup_matrix(const std::vector& res, + const std::vector& dydp, + const int npar, + const int ndat, + std::vector& beta, + std::vector& alpha); + void solve_matrix(const std::vector& beta, + const std::vector& alpha, + const int npar, + std::vector& dp); + float invrt_matrix(std::vector& alphaf, const int npar); + }; + +} //end namespace gshf +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.cxx new file mode 100644 index 0000000..f16d57a --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.cxx @@ -0,0 +1,123 @@ +#include +#include +#include + +#include "PeakFitterMrqdt.h" + +namespace examples { + +//-------------------------- + //constructor + + PeakFitterMrqdt::PeakFitterMrqdt(const PeakFitterMrqdtCfg& cfg) + : fMinWidth(cfg.fMinWidth) + , fMaxWidthMult(cfg.fMaxWidthMult) + , fPeakRange(cfg.fPeakRange) + , fAmpRange(cfg.fAmpRange) {} + + //------------------------ + //output parameters should be replaced with a returned value + void PeakFitterMrqdt::findPeakParameters(const std::vector& signal, + const ICandidateHitFinder::HitCandidateVec& fhc_vec, + PeakParamsVec& mhpp_vec, + double& chi2PerNDF, + int& NDF) const + { + if (fhc_vec.empty()) return; + + const float chiCut = 1e-3; + float lambda = 0.001; /* Marquardt damping parameter */ + float chiSqr = std::numeric_limits::max(), dchiSqr = std::numeric_limits::max(); + int nParams = 0; + + int startTime = fhc_vec.front().startTick; + int endTime = fhc_vec.back().stopTick; + + int roiSize = endTime - startTime; + + // init to a large number in case the fit fails + chi2PerNDF = chiSqr; + + std::vector y(roiSize); + std::vector p(3 * fhc_vec.size()); + std::vector plimmin(3 * fhc_vec.size()); + std::vector plimmax(3 * fhc_vec.size()); + std::vector perr(3 * fhc_vec.size()); + + /* choose the fit function and set the parameters */ + nParams = 0; + + for (size_t ih = 0; ih < fhc_vec.size(); ih++) { + float const peakMean = + fhc_vec[ih].hitCenter - 0.5 - (float)startTime; //shift by 0.5 to account for bin width + float const peakWidth = fhc_vec[ih].hitSigma; + float const amplitude = fhc_vec[ih].hitHeight; + float const meanLowLim = fmax(peakMean - fPeakRange * peakWidth, 0.); + float const meanHiLim = fmin(peakMean + fPeakRange * peakWidth, (float)roiSize); + p[0 + nParams] = amplitude; + p[1 + nParams] = peakMean; + p[2 + nParams] = peakWidth; + + plimmin[0 + nParams] = amplitude * 0.1; + plimmin[1 + nParams] = meanLowLim; + plimmin[2 + nParams] = std::max(fMinWidth, 0.1 * peakWidth); + + plimmax[0 + nParams] = amplitude * fAmpRange; + plimmax[1 + nParams] = meanHiLim; + plimmax[2 + nParams] = fMaxWidthMult * peakWidth; + + nParams += 3; + } + int fitResult = -1; + + for (size_t idx = 0; idx < size_t(roiSize); idx++) { + y[idx] = signal[startTime + idx]; + } + + int trial = 0; + lambda = -1.; /* initialize lambda on first call */ + do { + fitResult = fMarqFitAlg->gshf::MarqFitAlg::mrqdtfit( + lambda, &p[0], &plimmin[0], &plimmax[0], &y[0], nParams, roiSize, chiSqr, dchiSqr); + trial++; + if (fitResult || (trial > 100)) break; + } while (fabs(dchiSqr) >= chiCut); + + if (!fitResult) { + int fitResult2 = + fMarqFitAlg->gshf::MarqFitAlg::cal_perr(&p[0], &y[0], nParams, roiSize, &perr[0]); + if (!fitResult2) { + NDF = roiSize - nParams; + chi2PerNDF = chiSqr / NDF; + int parIdx = 0; + for (size_t i = 0; i < fhc_vec.size(); i++) { + + /* stand alone method + mhpp_vec[i].peakAmplitude = p[parIdx + 0]; + mhpp_vec[i].peakAmplitudeError = perr[parIdx + 0]; + mhpp_vec[i].peakCenter = p[parIdx + 1] + 0.5 + float(startTime); + mhpp_vec[i].peakCenterError = perr[parIdx + 1]; + mhpp_vec[i].peakSigma = p[parIdx + 2]; + mhpp_vec[i].peakSigmaError = perr[parIdx + 2]; + */ + + PeakFitParams_t mhpp; + mhpp.peakAmplitude = p[parIdx + 0]; + mhpp.peakAmplitudeError = perr[parIdx + 0]; + mhpp.peakCenter = + p[parIdx + 1] + 0.5 + float(startTime); //shift by 0.5 to account for bin width + mhpp.peakCenterError = perr[parIdx + 1]; + mhpp.peakSigma = std::abs(p[parIdx + 2]); + mhpp.peakSigmaError = perr[parIdx + 2]; + + mhpp_vec.emplace_back(mhpp); + + parIdx += 3; + } + + // fitStat=0; + } + } + return; + } +} diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.h new file mode 100644 index 0000000..f3e1901 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/PeakFitterMrqdt.h @@ -0,0 +1,39 @@ +#ifndef PEAKFITTERMRQDT_H +#define PEAKFITTERMRQDT_H + +#include +#include + +#include "ICandidateHitFinder.h" +#include "IPeakFitter.h" +#include "MarqFitAlg.h" //marqfit functions + +namespace examples { + + struct PeakFitterMrqdtCfg { + double fMinWidth; + double fMaxWidthMult; + double fPeakRange; + double fAmpRange; + }; + + class PeakFitterMrqdt : public IPeakFitter { + public: + explicit PeakFitterMrqdt(const PeakFitterMrqdtCfg& cfg); + + void findPeakParameters(const std::vector&, + const ICandidateHitFinder::HitCandidateVec&, + PeakParamsVec&, + double&, + int&) const override; + + private: + const double fMinWidth; + const double fMaxWidthMult; + const double fPeakRange; + const double fAmpRange; + + std::unique_ptr fMarqFitAlg; + }; +} +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/RawTypes.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/RawTypes.h new file mode 100644 index 0000000..4a06234 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/RawTypes.h @@ -0,0 +1,42 @@ +#ifndef RAW_TYPES_H +#define RAW_TYPES_H + +#include // std::numeric_limits<> + +namespace raw { + + typedef enum _compress { + kNone, ///< no compression + kHuffman, ///< Huffman Encoding + kZeroSuppression, ///< Zero Suppression algorithm + kZeroHuffman, ///< Zero Suppression followed by Huffman Encoding + kDynamicDec, ///< Dynamic decimation + kFibonacci ///< Fibonacci encoding + } Compress_t; + + typedef enum _auxdettype { + kUnknownAuxDet, ///< no idea + kScintillator, ///< Scintillator paddle + kTimeOfFlight, ///< Time of flight + kCherenkov ///< Cherenkov counter + } AuxDetType_t; + + /// Type representing a TDC tick + typedef int TDCtick_t; + + /// Type representing the ID of a readout channel + typedef unsigned int ChannelID_t; + + /// ID of an invalid channel + constexpr ChannelID_t InvalidChannelID = std::numeric_limits::max(); + + /// Returns whether the specified channel ID is valid + /// @note This does not mean that channel exists in the current geometry. + inline constexpr bool isValidChannelID(raw::ChannelID_t channel) + { + return channel != InvalidChannelID; + } + +} // namespace raw + +#endif diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.cxx new file mode 100644 index 0000000..650664e --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.cxx @@ -0,0 +1,36 @@ +/** **************************************************************************** + * @file Wire.cxx + * @brief Definition of basic channel signal object. + * @author brebel@fnal.gov + * @see Wire.h + * + * ****************************************************************************/ + +#include "Wire.h" + +// C/C++ standard libraries +#include // std::move() + +namespace recob { + + //---------------------------------------------------------------------- + Wire::Wire() : fChannel(raw::InvalidChannelID), fView(geo::kUnknown), fSignalROI() {} + + //---------------------------------------------------------------------- + Wire::Wire(RegionsOfInterest_t const& sigROIlist, raw::ChannelID_t channel, geo::View_t view) + : fChannel(channel), fView(view), fSignalROI(sigROIlist) + {} + + //---------------------------------------------------------------------- + Wire::Wire(RegionsOfInterest_t&& sigROIlist, raw::ChannelID_t channel, geo::View_t view) + : fChannel(channel), fView(view), fSignalROI(std::move(sigROIlist)) + {} + + //---------------------------------------------------------------------- + std::vector Wire::Signal() const + { + return {fSignalROI.begin(), fSignalROI.end()}; + } // Wire::Signal() + +} +//////////////////////////////////////////////////////////////////////// diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.h new file mode 100644 index 0000000..2cee4a0 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/Wire.h @@ -0,0 +1,236 @@ +/** **************************************************************************** + * @file lardataobj/RecoBase/Wire.h + * @brief Declaration of basic channel signal object. + * @author brebel@fnal.gov + * @see lardataobj/RecoBase/Wire.cxx + */ + +/* + * Changes + * 20190510 Gianluca Petrillo (petrillo@slac.stanford.edu) + * updated documentation + * 20141211 Gianluca Petrillo (petrillo@fnal.gov) + * data architecture revision changes: + * - fSignalType and SignalType() removed + * - fRawDigit and RawDigit() removed + * - fChannel and Channel() added + * - constructors now take directly a RawDigit, not its art::Ptr + * + * ****************************************************************************/ + +#ifndef LARDATAOBJ_RECOBASE_WIRE_H +#define LARDATAOBJ_RECOBASE_WIRE_H + +// LArSoft libraries +#include "RawTypes.h" // raw::ChannelID_t +#include "geo_types.h" +#include "sparse_vector.h" + +// C/C++ standard libraries +#include // std::size_t +#include + +namespace recob { + + /** + * @brief Class holding the regions of interest of signal from a channel. + * @note Despite the name, this class is associated to a readout channel, not + * just a wire + * + * The channel content is expected to have been filtered from noise and + * corrected for electronics response, a process often called in jargon + * "deconvolution". + * + * The content is presented as calibrated ADC counts, pedestal removed, as + * function of time in discrete TDC units. The time is expected to be the same + * as for the `raw::RawDigit` that originates it, i.e. starting from + * @ref DetectorClocksTPCelectronicsStartTime "TPC electronics start time" + * (use `detinfo::DetectorClocks` to discover the exact extent of each tick). + * + * The content is organized as time intervals where some signal is present + * ("regions of interest", RoI), outside which we assume no signal. + * By principle, the definition of the regions of interest is a negative one: + * we determine time intervals where we are confident no signal is present; + * the rest will constitute regions of interest where there _might_ be signal. + * The identification of such regions is responsibility of the algorithm + * creating the `Wire` object. In the simplest approach, the whole readout + * window is stored in a single region of interest, meaning that we don't + * claim any of the channel signal to be definitely signal free. + * More generally, the first tick of the waveform is #0, and if the first + * region of interest starts after that tick, it implies that the region + * between tick #0 and the start of that first region lacks signal. + * Likewise, any samples in the end of the covered time window (defined above) + * which lack signal are indicated by the overall size of the content, while + * the last region of interest ends earlier. + * + * Algorithms using the regions of interest can access the channel signal + * information either ignoring the regions of interest, and being potentially + * flooded by zeroes from the non-signal regions: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * for (float ADCcount: wire.Signal()) ... + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * or they can analyze region by region: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * for (auto iROI = wire.begin_range(); iROI != wire.end_range(); ++iROI) { + * const datarange_t& ROI = *iROI; + * const int FirstTick = ROI.begin_index(); + * const int EndTick = ROI.end_index(); + * const float FirstADC = ROI[FirstTick]; // index access is by absolute tick number + * for (float ADC: ROI) // ... or you can just iterate through + * // ... + * } // for + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * An alternative to the first form is: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * for (float ADCcount: wire.SignalROI()) ... + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * which does not create a temporary dense vector, as `Signal()` does instead. + * + * Note that the indexed access is always by absolute tick number. + * More examples of the use of `SignalROI()` return value are documented in + * `lar::sparse_vector`. + * + * Each channel is associated with a `raw::RawDigit` object. These + * associations should be stored together with `recob::Wire` by the producer + * in a `art::Assns` data product. + * + * + * Creating `recob::Wire` objects + * =============================== + * + * LArSoft "protocol" prescribes: + * + * * creation of an association of each `recob::Wire` with the `raw::RawDigit` + * it comes from; the association should be created by the same _art_ module + * producing the `recob::Wire` collection + * * `recob::Wire` time span should be the same as its `raw::RawDigit` + * + * Two patterns can be followed for creating a `recob::Wire`: + * + * 1. direct construction: see the constructor documentation + * 2. with `recob::WireCreator`, which extracts some relevant information from + * the source `raw::RawDigit` but _does not help with their association_ + * + * In both cases, please read the documentation of `recob::Wire` constructors. + */ + class Wire { + public: + /// a region of interest is a pair (TDC offset, readings) + typedef lar::sparse_vector RegionsOfInterest_t; + + /// Default constructor: a wire with no signal information + Wire(); + + private: + raw::ChannelID_t fChannel; ///< ID of the associated channel. + geo::View_t fView; ///< View corresponding to the plane of this wire. + RegionsOfInterest_t fSignalROI; ///< Signal on the channel as function of time tick. + + friend class WireCreator; // helper to create wires in art + + public: + // --- BEGIN -- Constructors --------------------------------------------- + /** + * @brief Constructor: uses specified signal in regions of interest. + * @param sigROIlist signal organized in regions of interest + * @param channel the ID of the channel + * @param view the view the channel belongs to + * + * Signal is copied into the `recob::Wire` object, including the sparse + * region of interest structure within `sigROIlist`. + * If possible, use the other constructor that moves the data instead. + * + * For more details, see the other constructor documentation. + */ + Wire(RegionsOfInterest_t const& sigROIlist, raw::ChannelID_t channel, geo::View_t view); + + /** + * @brief Constructor: uses specified signal in regions of interest. + * @param sigROIlist signal organized in regions of interest + * @param channel the ID of the channel + * @param view the view the channel belongs to + * + * The `recob::Wire` object is constructed with the waveform information + * in `sigROIlist` and assigned the specified `channel` and `view`. + * + * The signal is stored in a sparse vector, each entry corresponding to a + * tick in the calibrated waveform. The tick range of the sparse vector + * reflects the one in the wire, i.e. the first sample in `sigROIlist` + * becomes the sample #0 of the `recob::Wire` waveform. + * The total length of the waveform (that is, its duration in ticks) is + * also learned from the (nominal) size of `sigROIlist` (see also + * `lar::sparse_vector::resize()`), which can and should extend beyond + * the last region of interest. + * + * This constructor moves the signal information is moved `sigROIlist`, + * that becomes invalid. + * This also preserves the sparse region of interest structure within + * `sigROIlist`. + */ + Wire(RegionsOfInterest_t&& sigROIlist, raw::ChannelID_t channel, geo::View_t view); + // --- END -- Constructors ----------------------------------------------- + + // --- BEGIN -- Accessors ------------------------------------------------ + ///@name Accessors + ///@{ + + /// Return a zero-padded full length vector filled with RoI signal + std::vector Signal() const; + + /// Returns the list of regions of interest + const RegionsOfInterest_t& SignalROI() const; + + /// Returns the number of time ticks, or samples, in the channel + std::size_t NSignal() const; + + /// Returns the view the channel belongs to + geo::View_t View() const; + + /// Returns the ID of the channel (or InvalidChannelID) + raw::ChannelID_t Channel() const; + + ///@} + // --- END -- Accessors -------------------------------------------------- + + // --- BEGIN -- Sorting and comparison operations ------------------------ + /// @name Sorting and comparison operations + /// @{ + + /// Returns whether this channel ID is smaller than the other + bool operator<(const Wire& than) const; + + // --- END -- Sorting and comparison operations -------------------------- + + }; // class Wire + +} // namespace recob + +//------------------------------------------------------------------------------ +//--- inline implementation +//------------------------------------------------------------------------------ +inline const recob::Wire::RegionsOfInterest_t& recob::Wire::SignalROI() const +{ + return fSignalROI; +} +inline std::size_t recob::Wire::NSignal() const +{ + return fSignalROI.size(); +} +inline geo::View_t recob::Wire::View() const +{ + return fView; +} +inline raw::ChannelID_t recob::Wire::Channel() const +{ + return fChannel; +} +inline bool recob::Wire::operator<(const Wire& than) const +{ + return Channel() < than.Channel(); +} + +//------------------------------------------------------------------------------ + +#endif // LARDATAOBJ_RECOBASE_WIRE_H + +//////////////////////////////////////////////////////////////////////// diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.cxx b/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.cxx new file mode 100644 index 0000000..4972b75 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.cxx @@ -0,0 +1,64 @@ +/** + * @file larcoreobj/SimpleTypesAndConstants/geo_types.cxx + * @brief Definition of data types for geometry description (implementation). + * @see larcoreobj/SimpleTypesAndConstants/geo_types.h + * @ingroup Geometry + */ + +// header library +#include "geo_types.h" + +// C++ standard libraries +#include +#include // std::logic_error + +namespace geo { + + std::ostream& operator<<(std::ostream& os, Coordinate const coordinate) + { + switch (coordinate) { + case Coordinate::X: os << 'X'; break; + case Coordinate::Y: os << 'Y'; break; + case Coordinate::Z: os << 'Z'; break; + } + return os; + } + + std::ostream& operator<<(std::ostream& os, DriftSign const sign) + { + switch (sign) { + case DriftSign::Positive: os << '+'; break; + case DriftSign::Negative: os << '-'; break; + case DriftSign::Unknown: os << '?'; break; + } + return os; + } + + bool operator==(DriftAxis const a, DriftAxis const b) + { + return a.coordinate == b.coordinate and a.sign == b.sign; + } + + bool operator!=(DriftAxis const a, DriftAxis const b) + { + return not(a == b); + } + + std::ostream& operator<<(std::ostream& os, DriftAxis const driftAxis) + { + return os << driftAxis.sign << driftAxis.coordinate; + } + + std::string SignalTypeName(SigType_t sigType) + { + switch (sigType) { + case kInduction: return "induction"; + case kCollection: return "collection"; + case kMysteryType: return "unknown"; + } // switch + throw std::logic_error("geo::SignalTypeName(): unexpected signal type #" + + std::to_string(static_cast(sigType))); + } +} + +// ----------------------------------------------------------------------------- diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.h new file mode 100644 index 0000000..a321313 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/geo_types.h @@ -0,0 +1,682 @@ +/** + * @defgroup Geometry Detector geometry information + * @file larcoreobj/SimpleTypesAndConstants/geo_types.h + * @brief Definition of data types for geometry description. + * @author Brian Rebel (brebel@fnal.gov) + * @see larcoreobj/SimpleTypesAndConstants/geo_types.cxx + * @ingroup Geometry + * + * This library depends only on standard C++. + * + */ + +#ifndef LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_TYPES_H +#define LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_TYPES_H + +// C++ standard libraries +#include +#include +#include // std::ostream +#include // std::numeric_limits<> +#include +#include +#include +#include + +namespace geo::details { + /// Write the argument into a string + template + std::string writeToString(T const& value); + + template + struct has_parent : std::false_type {}; + + template + struct has_parent> : std::true_type {}; + + /// Whether `ID` represents an element on top of the hierarchy. + template + constexpr bool isTopGeoElementID = !has_parent::value; + + template + constexpr auto index_impl(std::size_t& index) + { + if constexpr (has_parent::value) { index_impl(++index); } + } + + template + constexpr auto const index_for() + { + std::size_t index{}; + index_impl(index); + return index; + } + + template + struct AbsIDtypeStruct; + + template + using AbsIDtype = typename AbsIDtypeStruct::type; + + template + constexpr auto getAbsIDindex(ID const& id) + { + static_assert(Index <= ID::Level, "Index not available for this type."); + if constexpr (Index == ID::Level) + return id.deepestIndex(); + else + return getAbsIDindex(id.parentID()); + } + + template + auto& getAbsIDindex(ID& id) + { + static_assert(Index <= ID::Level, "Index not available for this type."); + if constexpr (Index == ID::Level) + return id.deepestIndex(); + else + return getAbsIDindex(id.parentID()); + } + +} // namespace geo::details + +// BEGIN Geometry -------------------------------------------------------------- +/** + * @addtogroup Geometry + * @see `geo::GeometryCore` + * + * The description of the detector as seen by LArSoft is accessed via LArSoft + * "geometry system". + * + * A high level explanation of the geometry description is present in the + * official [LArSoft collaboration web site](http://larsoft.org), at the bottom + * of which further sources are detailed. + * + * The geometry system revolves around the `geo::GeometryCore` service provider, + * which is the hub to access all the geometry information. + * The information accessible via this interface include: + * * geometric description of the TPC components (wires, cryostats...) + * * geometric description of the optical detectors + * * geometric description of the "auxiliary" detectors (i.e. every detector + * which is not the two above) + * * a ROOT representation of the detector geometry + * * relation between readout channels and physical components (especially TPC + * readout channels vs. wires) + * * tools to navigate this information + * * ... and a load more stuff. + * + * The service provider can be obtained including `larcore/Geometry/Geometry.h` + * header and calling: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * auto const& geom = *(lar::providerFrom()); + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + */ +/// @{ +namespace geo { + + // --- BEGIN -- Geometry enumerators ----------------------------------------- + /// @name Geometry enumerators + /// @{ + + enum class Coordinate { X, Y, Z }; + std::ostream& operator<<(std::ostream& os, Coordinate); + constexpr int to_int(Coordinate const coord) noexcept + { + return static_cast(coord); + } + + /// Enumerate the possible plane projections + typedef enum _plane_proj { + kU, ///< Planes which measure U. + kV, ///< Planes which measure V. + kW, ///< Planes which measure W (third view for Bo, MicroBooNE, etc). + kZ = kW, ///< Planes which measure Z direction. + kY, ///< Planes which measure Y direction. + kX, ///< Planes which measure X direction. + k3D, ///< 3-dimensional objects, potentially hits, clusters, prongs, etc. + kUnknown ///< Unknown view. + } View_t; + + typedef enum _plane_orient { + kHorizontal, ///< Planes that are in the horizontal plane. + kVertical ///< Planes that are in the vertical plane (e.g. ArgoNeuT). + } Orient_t; + + typedef enum _plane_sigtype { + kInduction, ///< Signal from induction planes. + kCollection, ///< Signal from collection planes. + kMysteryType ///< Who knows? + } SigType_t; + + /** + * @brief Drift sign: positive or negative + */ + enum class DriftSign { + Unknown, ///< Drift direction is unknown. + Positive, ///< Drift towards positive values. + Negative ///< Drift towards negative values. + }; + std::ostream& operator<<(std::ostream& os, DriftSign); + + constexpr int to_int(DriftSign const sign) + { + switch (sign) { + case DriftSign::Positive: return 1; + case DriftSign::Negative: return -1; + case DriftSign::Unknown: break; + } + return 0; + } + + struct DriftAxis { + Coordinate coordinate; + DriftSign sign; + }; + bool operator==(DriftAxis a, DriftAxis b); + bool operator!=(DriftAxis a, DriftAxis b); + std::ostream& operator<<(std::ostream& os, DriftAxis); + + /// @} + // --- END -- Geometry enumerators ------------------------------------------- + + /// @{ + /// @name Geometry element IDs + + /// The data type to uniquely identify a cryostat. + struct CryostatID { + using CryostatID_t = unsigned int; ///< Type for the ID number. + + // not constexpr because we would need an implementation file to define it + /// Special code for an invalid ID. + static constexpr CryostatID_t InvalidID = std::numeric_limits::max(); + + bool isValid = false; ///< Whether this ID points to a valid element. + CryostatID_t Cryostat = InvalidID; ///< Index of cryostat. + + /// Default constructor: an invalid cryostat. + constexpr CryostatID() = default; + + /// Constructor: valid ID of cryostat with index c + explicit constexpr CryostatID(CryostatID_t c) : isValid(true), Cryostat(c) {} + + /// Constructor: valid ID of cryostat with index c + constexpr CryostatID(CryostatID_t c, bool valid) : isValid(valid), Cryostat(c) {} + + static constexpr auto first() { return CryostatID{0, true}; } + /// @{ + /// @name ID validity + + /// Returns true if the ID is valid + explicit constexpr operator bool() const { return isValid; } + + /// Sets the validity of the ID. + void setValidity(bool valid) { isValid = valid; } + + /// Sets the ID as valid. + void markValid() { setValidity(true); } + + /// Sets the ID as invalid. + void markInvalid() { setValidity(false); } + /// @} + + // comparison operators are out of class + + //@{ + /// Human-readable representation of the cryostat ID. + std::string toString() const { return details::writeToString(*this); } + explicit operator std::string() const { return toString(); } + //@} + + // the following four methods are useful for (templated) abstraction + /// Returns the value of the deepest ID available (cryostat's). + constexpr auto const& deepestIndex() const { return Cryostat; } + /// Returns the deepest ID available (cryostat's). + auto& deepestIndex() { return Cryostat; } + /// Returns the index level `Index` of this type. + template + constexpr auto getIndex() const; + + /// Level of this element. + static constexpr auto Level = details::index_for(); + + /// Return the value of the invalid ID as a r-value + static constexpr CryostatID_t getInvalidID() { return CryostatID::InvalidID; } + + }; // struct CryostatID + + /// The data type to uniquely identify a optical detector. + struct OpDetID : public CryostatID { + using OpDetID_t = unsigned int; ///< Type for the ID number. + + using ParentID_t = CryostatID; ///< Type of the parent ID. + + // not constexpr because we would need an implementation file to define it + /// Special code for an invalid ID. + static constexpr OpDetID_t InvalidID = std::numeric_limits::max(); + + /// Index of the optical detector within its cryostat. + OpDetID_t OpDet = InvalidID; + + /// Default constructor: an invalid optical detector ID. + constexpr OpDetID() = default; + + /// Constructor: optical detector with index `o` in the cryostat identified by + /// `cryoid` + constexpr OpDetID(CryostatID const& cryoid, OpDetID_t o) : CryostatID(cryoid), OpDet(o) {} + + /// Constructor: opdtical detector with index `o` in the cryostat index `c` + constexpr OpDetID(CryostatID_t c, OpDetID_t o) : CryostatID(c), OpDet(o) {} + + static constexpr OpDetID first() { return OpDetID{CryostatID::first(), 0}; } + + // comparison operators are out of class + + //@{ + /// Human-readable representation of the optical detector ID. + std::string toString() const { return details::writeToString(*this); } + explicit operator std::string() const { return toString(); } + //@} + + // the following two methods are useful for (templated) abstraction + /// Returns the value of the deepest ID available (OpDet's). + constexpr auto const& deepestIndex() const { return OpDet; } + /// Returns the deepest ID available (OpDet's). + auto& deepestIndex() { return OpDet; } + /// Return the parent ID of this one (a cryostat ID). + constexpr ParentID_t const& parentID() const { return *this; } + constexpr ParentID_t& parentID() { return *this; } + /// Returns the index level `Index` of this type. + template + constexpr auto getIndex() const; + + /// Conversion to CryostatID (for convenience of notation). + constexpr CryostatID const& asCryostatID() const { return parentID(); } + constexpr CryostatID& asCryostatID() { return parentID(); } + + /// Level of this element. + static constexpr auto Level = details::index_for(); + + /// Return the value of the invalid optical detector ID as a r-value + static constexpr OpDetID_t getInvalidID() { return OpDetID::InvalidID; } + + }; // struct OpDetID + + /// The data type to uniquely identify a TPC. + struct TPCID : public CryostatID { + using TPCID_t = unsigned int; ///< Type for the ID number. + + using ParentID_t = CryostatID; ///< Type of the parent ID. + + // not constexpr because we would need an implementation file to define it + /// Special code for an invalid ID. + static constexpr TPCID_t InvalidID = std::numeric_limits::max(); + + TPCID_t TPC = InvalidID; ///< Index of the TPC within its cryostat. + + /// Default constructor: an invalid TPC ID. + constexpr TPCID() = default; + + /// Constructor: TPC with index t in the cryostat identified by cryoid + constexpr TPCID(CryostatID const& cryoid, TPCID_t t) : CryostatID(cryoid), TPC(t) {} + + /// Constructor: TPC with index t in the cryostat index c + constexpr TPCID(CryostatID_t c, TPCID_t t) : CryostatID(c), TPC(t) {} + + TPCID next() const { return {Cryostat, TPC + 1}; } + + static constexpr auto first() { return TPCID{CryostatID::first(), 0}; } + static constexpr auto first(CryostatID const& id) { return TPCID{id, 0}; } + + // comparison operators are out of class + + //@{ + /// Human-readable representation of the TPC ID. + std::string toString() const { return details::writeToString(*this); } + explicit operator std::string() const { return toString(); } + //@} + + // the following two methods are useful for (templated) abstraction + /// Returns the value of the deepest ID available (TPC's). + constexpr auto const& deepestIndex() const { return TPC; } + /// Returns the deepest ID available (TPC's). + auto& deepestIndex() { return TPC; } + /// Return the parent ID of this one (a cryostat ID). + constexpr ParentID_t const& parentID() const { return *this; } + constexpr ParentID_t& parentID() { return *this; } + /// Returns the index level `Index` of this type. + template + constexpr auto getIndex() const; + + /// Conversion to TPCID (for convenience of notation). + constexpr CryostatID const& asCryostatID() const { return parentID(); } + constexpr CryostatID& asCryostatID() { return parentID(); } + + /// Level of this element. + static constexpr auto Level = details::index_for(); + + /// Return the value of the invalid TPC ID as a r-value + static constexpr TPCID_t getInvalidID() { return TPCID::InvalidID; } + + }; // struct TPCID + + /// The data type to uniquely identify a Plane. + struct PlaneID : public TPCID { + using PlaneID_t = unsigned int; ///< Type for the ID number. + + using ParentID_t = TPCID; ///< Type of the parent ID. + + // not constexpr because we would need an implementation file to define it + /// Special code for an invalid ID. + static constexpr PlaneID_t InvalidID = std::numeric_limits::max(); + + PlaneID_t Plane = InvalidID; ///< Index of the plane within its TPC. + + /// Default constructor: an invalid plane ID. + constexpr PlaneID() = default; + + /// Constructor: plane with index p in the TPC identified by tpcid + constexpr PlaneID(TPCID const& tpcid, PlaneID_t p) : TPCID(tpcid), Plane(p) {} + + /// Constructor: plane with index p in the cryostat index c, TPC index t + constexpr PlaneID(CryostatID_t c, TPCID_t t, PlaneID_t p) : TPCID(c, t), Plane(p) {} + + static constexpr auto first() { return PlaneID{TPCID::first(), 0}; } + static constexpr auto first(CryostatID const& id) { return PlaneID{TPCID::first(id), 0}; } + static constexpr auto first(TPCID const& id) { return PlaneID{id, 0}; } + + // comparison operators are out of class + + //@{ + /// Human-readable representation of the plane ID. + std::string toString() const { return details::writeToString(*this); } + explicit operator std::string() const { return toString(); } + //@} + + // the following two methods are useful for (templated) abstraction + /// Returns the value of the deepest ID available (plane's). + constexpr auto const& deepestIndex() const { return Plane; } + /// Returns the deepest ID available (plane's). + auto& deepestIndex() { return Plane; } + /// Return the parent ID of this one (a TPC ID). + constexpr ParentID_t const& parentID() const { return *this; } + constexpr ParentID_t& parentID() { return *this; } + /// Returns the index level `Index` of this type. + template + constexpr auto getIndex() const; + + /// Conversion to PlaneID (for convenience of notation). + constexpr TPCID const& asTPCID() const { return parentID(); } + constexpr TPCID& asTPCID() { return parentID(); } + + /// Level of this element. + static constexpr auto Level = details::index_for(); + + /// Return the value of the invalid plane ID as a r-value + static constexpr PlaneID_t getInvalidID() { return PlaneID::InvalidID; } + + }; // struct PlaneID + + // The data type to uniquely identify a code wire segment. + struct WireID : public PlaneID { + using WireID_t = unsigned int; ///< Type for the ID number. + + using ParentID_t = PlaneID; ///< Type of the parent ID. + + // not constexpr because we would need an implementation file to define it + /// Special code for an invalid ID. + static constexpr WireID_t InvalidID = std::numeric_limits::max(); + + WireID_t Wire = InvalidID; ///< Index of the wire within its plane. + + /// Default constructor: an invalid TPC ID. + constexpr WireID() = default; + + /// Constructor: wire with index w in the plane identified by planeid + constexpr WireID(PlaneID const& planeid, WireID_t w) : PlaneID(planeid), Wire(w) {} + + /// Constructor: wire with index w in cryostat index c, TPC index t, plane index p + constexpr WireID(CryostatID_t c, TPCID_t t, PlaneID_t p, WireID_t w) : PlaneID(c, t, p), Wire(w) + {} + + static constexpr auto first() { return WireID{PlaneID::first(), 0}; } + static constexpr auto first(CryostatID const& id) { return WireID{PlaneID::first(id), 0}; } + static constexpr auto first(TPCID const& id) { return WireID{PlaneID::first(id), 0}; } + static constexpr auto first(PlaneID const& id) { return WireID{id, 0}; } + + //@{ + /// Human-readable representation of the wire ID. + std::string toString() const { return details::writeToString(*this); } + explicit operator std::string() const { return toString(); } + //@} + + // the following two methods are useful for (templated) abstraction + /// Returns the value of the deepest ID available (wire's). + constexpr auto const& deepestIndex() const { return Wire; } + /// Returns the deepest ID available (wire's). + auto& deepestIndex() { return Wire; } + /// Return the parent ID of this one (a plane ID). + constexpr ParentID_t const& parentID() const { return *this; } + constexpr ParentID_t& parentID() { return *this; } + /// Returns the index level `Index` of this type. + template + constexpr auto getIndex() const; + + /// Conversion to WireID (for convenience of notation). + constexpr PlaneID const& asPlaneID() const { return parentID(); } + constexpr PlaneID& asPlaneID() { return parentID(); } + + /// Backward compatibility; use the wire directly or a explicit cast instead + /// @todo Remove the instances of geo::WireID::planeID() in the code + constexpr PlaneID const& planeID() const { return *this; } + + /// Level of this element. + static constexpr auto Level = details::index_for(); + + /// Return the value of the invalid wire ID as a r-value + static constexpr WireID_t getInvalidID() { return WireID::InvalidID; } + + }; // struct WireID + + //---------------------------------------------------------------------------- + //--- ID output operators + //--- + /// Generic output of CryostatID to stream + inline std::ostream& operator<<(std::ostream& out, CryostatID const& cid) + { + out << "C:" << cid.Cryostat; + return out; + } + + /// Generic output of OpDetID to stream. + inline std::ostream& operator<<(std::ostream& out, OpDetID const& oid) + { + out << oid.parentID() << " O:" << oid.OpDet; + return out; + } + + /// Generic output of TPCID to stream + inline std::ostream& operator<<(std::ostream& out, TPCID const& tid) + { + out << tid.parentID() << " T:" << tid.TPC; + return out; + } + + /// Generic output of PlaneID to stream + inline std::ostream& operator<<(std::ostream& out, PlaneID const& pid) + { + out << pid.parentID() << " P:" << pid.Plane; + return out; + } + + /// Generic output of WireID to stream + inline std::ostream& operator<<(std::ostream& out, WireID const& wid) + { + out << wid.parentID() << " W:" << wid.Wire; + return out; + } + + /// @} + // Geometry element IDs + + //---------------------------------------------------------------------------- + //--- ID comparison operators + //--- + + /// @{ + /// @name ID comparison operators + /// @details The result of comparison with invalid IDs is undefined. + + /// Comparison: the IDs point to the same cryostat (validity is ignored) + inline constexpr bool operator==(CryostatID const& a, CryostatID const& b) + { + return a.Cryostat == b.Cryostat; + } + + /// Comparison: the IDs point to different cryostats (validity is ignored) + inline constexpr bool operator!=(CryostatID const& a, CryostatID const& b) + { + return !(a == b); + } + + /// Order cryostats with increasing ID + inline constexpr bool operator<(CryostatID const& a, CryostatID const& b) + { + return a.Cryostat < b.Cryostat; + } + + template + static constexpr bool is_base_of_strict{std::is_base_of{} && + !std::is_same{}}; + + /// Comparison: the IDs point to same optical detector (validity is ignored) + template + inline constexpr std::enable_if_t, bool> operator==( + GeoID const& a, + GeoID const& b) + { + return std::tie(a.parentID(), a.deepestIndex()) == std::tie(b.parentID(), b.deepestIndex()); + } + + /// Comparison: IDs point to different optical detectors (validity is ignored) + template + inline constexpr std::enable_if_t, bool> operator!=( + GeoID const& a, + GeoID const& b) + { + return !(a == b); + } + + /// Order OpDetID in increasing Cryo, then OpDet + template + inline constexpr std::enable_if_t, bool> operator<( + GeoID const& a, + GeoID const& b) + { + return std::tie(a.parentID(), a.deepestIndex()) < std::tie(b.parentID(), b.deepestIndex()); + } + + /// @} + + //---------------------------------------------------------------------------- + struct WireIDIntersection { + double y; ///< y position of intersection + double z; ///< z position of intersection + unsigned int TPC; ///< TPC of intersection + + // In APAs, we want this to increase in the direction wireID + // index increases in: moving inward vertically towards y=0 + bool operator<(WireIDIntersection const& otherIntersect) const + { + return std::abs(y) > std::abs(otherIntersect.y); + } + + static constexpr WireIDIntersection invalid() + { + return { + std::numeric_limits::infinity(), std::numeric_limits::infinity(), -1u}; + } + }; + + //---------------------------------------------------------------------------- + /// Returns the name of the specified signal type. + std::string SignalTypeName(geo::SigType_t sigType); + + //---------------------------------------------------------------------------- + +} // namespace geo +/// @} +// END Geometry ---------------------------------------------------------------- + +namespace geo::details { + + //-------------------------------------------------------------------------- + template + struct AbsIDtypeStruct { + static_assert(Index <= ID::Level, "Requested ID index is not available."); + using type = typename AbsIDtypeStruct::type; + }; + + template + struct AbsIDtypeStruct> { + using type = ID; + }; + + //-------------------------------------------------------------------------- + template + inline std::string writeToString(T const& value) + { + std::ostringstream sstr; + sstr << value; + return sstr.str(); + } + + //-------------------------------------------------------------------------- + +} // namespace geo::details + +//------------------------------------------------------------------------------ +//--- template implementation +//------------------------------------------------------------------------------ +template +constexpr auto geo::CryostatID::getIndex() const +{ + static_assert(Index <= Level, "This ID type does not have the requested Index level."); + return details::getAbsIDindex(*this); +} + +//------------------------------------------------------------------------------ +template +constexpr auto geo::OpDetID::getIndex() const +{ + static_assert(Index <= Level, "This ID type does not have the requested Index level."); + return details::getAbsIDindex(*this); +} + +//------------------------------------------------------------------------------ +template +constexpr auto geo::TPCID::getIndex() const +{ + static_assert(Index <= Level, "This ID type does not have the requested Index level."); + return details::getAbsIDindex(*this); +} + +//------------------------------------------------------------------------------ +template +constexpr auto geo::PlaneID::getIndex() const +{ + static_assert(Index <= Level, "This ID type does not have the requested Index level."); + return details::getAbsIDindex(*this); +} + +//------------------------------------------------------------------------------ +template +constexpr auto geo::WireID::getIndex() const +{ + static_assert(Index <= Level, "This ID type does not have the requested Index level."); + return details::getAbsIDindex(*this); +} + +//------------------------------------------------------------------------------ + +#endif // LARCOREOBJ_SIMPLETYPESANDCONSTANTS_GEO_TYPES_H diff --git a/gaus_hit_finder/copied_from_larsoft_minor_edits/sparse_vector.h b/gaus_hit_finder/copied_from_larsoft_minor_edits/sparse_vector.h new file mode 100644 index 0000000..d0ec8b5 --- /dev/null +++ b/gaus_hit_finder/copied_from_larsoft_minor_edits/sparse_vector.h @@ -0,0 +1,2512 @@ +/** + * @file lardataobj/Utilities/sparse_vector.h + * @brief Class defining a sparse vector (holes are zeroes) + * @author Gianluca Petrillo (petrillo@fnal.gov) + * @date April 22, 2014 + * @version 1.0 + * + */ + +#ifndef LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H +#define LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H + +// C/C++ standard library +#include // std::upper_bound(), std::max() +#include // std::ptrdiff_t +#include // std::distance() +#include // std::accumulate +#include +#include // std::out_of_range() +#include // std::is_integral +#include + +/// Namespace for generic LArSoft-related utilities. +namespace lar { + + // ----------------------------------------------------------------------------- + // --- utility classes for sparse_vector + // --- + + /** + * @brief Little class storing a value. + * @tparam T type of the stored value + * + * This class stores a constant value and returns it as conversion to type T. + * It also acts as a left-value of type T, except that the assigned value + * is ignored. + */ + template + class const_value_box { + typedef const_value_box this_t; + + public: + typedef T value_type; ///< type of the value stored + + /// Default constructor: stores default value + const_value_box() {} + + /// Constructor: stores the specified value + explicit const_value_box(value_type new_value) : value(new_value) {} + + /// Assignment: the assigned value is ignored + this_t& operator=(value_type) { return *this; } + + //@{ + /// Conversion to the basic type: always returns the stored value + operator value_type() const { return value; } + operator const value_type&() const { return value; } + //@} + + protected: + value_type value; ///< the value stored for delivery + }; // class const_value_box + + /// @brief A constant iterator returning always the same value + /// @tparam T type of the value returned by dereferenciation + template + class value_const_iterator { + typedef value_const_iterator this_t; ///< alias for this type + + public: + using iterator_category = std::random_access_iterator_tag; + using value_type = T; + using difference_type = std::ptrdiff_t; + using pointer = T*; + using reference = T&; + + /// Default constructor: use the default value + value_const_iterator() : value() {} + + /// Constructor: value that will be returned + value_const_iterator(value_type new_value) : value(new_value) {} + + /// Constructor: value to be returned and current iterator "position" + value_const_iterator(value_type new_value, difference_type offset) + : index(offset), value(new_value) + {} + + /// Returns a copy of the stored value + value_type operator*() const { return value; } + + /// Returns a copy of the stored value + value_type operator[](difference_type) const { return value; } + + /// Increments the position of the iterator, returns the new position + this_t& operator++() + { + ++index; + return *this; + } + + /// Increments the position of the iterator, returns the old position + this_t operator++(int) + { + this_t copy(*this); + ++index; + return copy; + } + + /// Decrements the position of the iterator, returns the new position + this_t& operator--() + { + --index; + return *this; + } + + /// Decrements the position of the iterator, returns the old position + this_t operator--(int) + { + this_t copy(*this); + --index; + return copy; + } + + /// Increments the position of the iterator by the specified steps + this_t& operator+=(difference_type ofs) + { + index += ofs; + return *this; + } + + /// Decrements the position of the iterator by the specified steps + this_t& operator-=(difference_type ofs) + { + index -= ofs; + return *this; + } + + /// Returns an iterator pointing ahead of this one by the specified steps + this_t operator+(difference_type ofs) const { return this_t(value, index + ofs); } + + /// Returns an iterator pointing behind this one by the specified steps + this_t operator-(difference_type ofs) const { return this_t(value, index - ofs); } + + /// @brief Returns the distance between this iterator and the other + /// @return how many steps this iterator is ahead of the other + difference_type operator-(const this_t& iter) const { return index - iter.index; } + + //@{ + /// Comparison operators: determined by the position pointed by the iterators + bool operator==(const this_t& as) const { return index == as.index; } + bool operator!=(const this_t& as) const { return index != as.index; } + + bool operator<(const this_t& as) const { return index < as.index; } + bool operator<=(const this_t& as) const { return index <= as.index; } + bool operator>(const this_t& as) const { return index > as.index; } + bool operator>=(const this_t& as) const { return index >= as.index; } + //@} + + protected: + difference_type index{0}; ///< (arbitrary) position pointed by the iterator + value_type value; ///< value to be returned when dereferencing + }; // class value_const_iterator<> + + /** + * @brief Returns an iterator pointing ahead of this one by the specified steps + * @param ofs number of steps ahead + * @param iter base iterator + */ + template + value_const_iterator operator+(typename value_const_iterator::difference_type ofs, + value_const_iterator& iter) + { + return {iter.value, iter.index + ofs}; + } + + /* // not ready yet... +template +class value_iterator: public value_const_iterator { + using base_t = value_const_iterator; + using this_t = value_iterator; + using const_value_box = value_box; + public: + + using typename base_t::value_type; + using typename base_t::reference; + + value_box operator* () const { return value_box(value); } + value_box operator[] (difference_type) const { return value_box(value); } + +}; // class value_iterator<> +*/ + + //------------------------------------------------------------------------------ + //--- lar::range_t + //--- + /** + * @brief A range (interval) of integers + * @tparam SIZE type of the indices (expected integral) + * + * Includes a first and an after-the-last value, and some relation metods. + */ + template + class range_t { + public: + typedef SIZE size_type; ///< type for the indices in the range + typedef std::ptrdiff_t difference_type; ///< type for index difference + + static_assert(std::is_integral::value, "range_t needs an integral type"); + + size_type offset; ///< offset (absolute index) of the first element + size_type last; ///< offset (absolute index) after the last element + + /// Default constructor: empty range + range_t() : offset(0), last(0) {} + + /// Constructor from first and last index + range_t(size_type from, size_type to) : offset(from), last(std::max(from, to)) {} + + /// Sets the borders of the range + void set(size_type from, size_type to) + { + offset = from; + last = std::max(from, to); + } + + /// Returns the first absolute index included in the range + size_type begin_index() const { return offset; } + + /// Returns the first absolute index not included in the range + size_type end_index() const { return last; } + + /// Returns the position within the range of the absolute index specified + size_type relative_index(size_type index) const { return index - offset; } + + /// Returns the size of the range + size_type size() const { return last - offset; } + + /// Moves the end of the range to fit the specified size + void resize(size_type new_size) { last = offset + new_size; } + + /// Moves the begin of the range by the specified amount + void move_head(difference_type shift) { offset += shift; } + + /// Moves the end of the range by the specified amount + void move_tail(difference_type shift) { last += shift; } + + /// Returns whether the range is empty + bool empty() const { return last <= offset; } + + /// Returns whether the specified absolute index is included in this range + bool includes(size_type index) const { return (index >= offset) && (index < last); } + + /// Returns whether the specified range is completely included in this one + bool includes(const range_t& r) const + { + return includes(r.begin_index()) && includes(r.end_index()); + } + + /// Returns if this and the specified range overlap + bool overlap(const range_t& r) const + { + return (begin_index() < r.end_index()) && (end_index() > r.begin_index()); + } + + /// Returns if there are elements in between this and the specified range + bool separate(const range_t& r) const + { + return (begin_index() > r.end_index()) || (end_index() < r.begin_index()); + } + + /// @brief Returns whether an index is within or immediately after this range + /// + /// Returns whether the specified absolute index is included in this range + /// or is immediately after it (not before it!) + bool borders(size_type index) const { return (index >= offset) && (index <= last); } + + //@{ + /// Sort: this range is smaller if its offset is smaller + bool operator<(const range_t& than) const { return less(*this, than); } + //@} + + /// Returns whether the specified range has our same offset and size + bool operator==(const range_t& as) const { return (offset == as.offset) && (last == as.last); } + + /// Returns whether the range is valid (that is, non-negative size) + bool is_valid() const { return last >= offset; } + + //@{ + /// Returns if a is "less" than b + static bool less(const range_t& a, const range_t& b) { return a.offset < b.offset; } + static bool less(const range_t& a, size_type b) { return a.offset < b; } + static bool less(size_type a, const range_t& b) { return a < b.offset; } + //@} + + /// Helper type to be used for binary searches + typedef bool (*less_int_range)(size_type, const range_t& b); + }; // class range_t + + // ----------------------------------------------------------------------------- + // --- lar::sparse_vector + // --- + + // the price of having used subclasses extensively: + template + class sparse_vector; + + namespace details { + template + decltype(auto) make_const_datarange_t(typename sparse_vector::datarange_t& r); + } // namespace details + + /** **************************************************************************** + * @brief A sparse vector + * @tparam T type of data stored in the vector + * @todo backward iteration; reverse iterators; iterator on non-void elements + * only; iterator on non-void elements only, returning a pair (index;value) + * + * A sparse_vector is a container of items marked by consecutive indices + * (that is, a vector like std::vector), where only non-zero elements are + * actually stored. + * The implementation is a container of ranges of non-zero consecutive values; + * the zero elements are effectively not stored in the object, and a zero is + * returned whenever they are accessed. + * In the following, the regions of zeros between the non-zero ranges are + * collectively called "the void". + * + * See sparse_vector_test.cc for examples of usage. + * + * Although some level of dynamic assignment is present, the class is not very + * flexible and it is best assigned just once, by adding ranges (add_range(); or + * by push_back(), which is less efficient). + * While this class mimics a good deal of the STL vector interface, it is + * *not* a std::vector and it does not support all the tricks of it. + * + * For the following, let's assume: + * ~~~ + * lar::sparse_vector sv; + * std::vector buffer; + * ~~~ + * + * Supported usage + * --------------- + * + * ~~~ + * for (const auto& value: sv) std::cout << " " << value; + * ~~~ + * ~~~ + * lar::sparse_vector::const_iterator iValue = sv.cbegin(), vend = sv.cend(); + * while (iSV != sv.end()) std::cout << " " << *(iSV++); + * ~~~ + * ~~~ + * for (auto value: sv) { ... } + * ~~~ + * Common iteration on all elements. Better to do a constant one. + * The first two examples print the full content of the sparse vector, void + * included. + * + * --- + * + * ~~~ + * sv[10] = 3.; + * ~~~ + * Assign to an *existing* (not void) element. Assigning to void elements + * is not supported, and there is no way to make an existing element to become + * void (assigning 0 will yield to an existing element with value 0). + * + * --- + * + * ~~~ + * sv.set_at(10, 3.); + * ~~~ + * Assign a value to an element. The element could be in the void; in any case, + * after this call the element will not be in the void anymore (even if the + * assigned value is zero; to cast a cell into the void, use unset_at()). + * + * --- + * + * ~~~ + * sv.add_range(20, buffer) + * ~~~ + * ~~~ + * sv.add_range(20, buffer.begin(), buffer.end()) + * ~~~ + * ~~~ + * sv.add_range(20, std::move(buffer)) + * ~~~ + * Add the content of buffer starting at the specified position. + * The last line will attempt to use the buffer directly (only happens if the + * start of the new range -- 20 in the example -- is in the void and therefore a + * new range will be added). The new range is merged with the existing ones when + * needed, and it overwrites their content in case of overlap. + * If the specified position is beyond the current end of the sparse vector, + * the gap will be filled by void. + * + * --- + * + * ~~~ + * sv.append(buffer) + * ~~~ + * ~~~ + * sv.append(buffer.begin(), buffer.end()) + * ~~~ + * ~~~ + * sv.append(std::move(buffer)) + * ~~~ + * Add the content of buffer at the end of the sparse vector. + * The last line will attempt to use the buffer directly (only happens if the + * end of the sparse vector is in the void and therefore a new range will be + * added). + * + * --- + * + * ~~~ + * sv.resize(30) + * ~~~ + * Resizes the sparse vector to the specified size. Truncation may occur, in + * which case the data beyond the new size is removed. If an extension occurs + * instead, the new area is void. + * + * --- + * + * ~~~ + * for (const lar::sparse_vector::datarange_t& range: sv.get_ranges()) { + * size_t first_item = range.begin_index(); // index of the first item + * size_t nItems = range.size(); // number of items in this range + * for (double value: range) { ... } + * } + * ~~~ + * A sparse vector can be parsed range by range, skipping the void. Each range + * object itself supports iteration. + * Neither the content nor the shape of the ranges can be changed this way. + * + * + * Possible future improvements + * ---------------------------- + * + * ~~~ + * sv[10] = 3.; + * ~~~ + * is not supported if the vector is shorter than 11 (similarly to a + * std::vector too), and not even if the item #10 is currently in the void. + * This latter could be changed to create a new element/range; this would + * require to ship the pointer to the container with the reference (the return + * value of "sv[10]"), just in case an assignment will occur, or to create the + * element regardless, like for std::map (but this would provoke a disaster + * if the caller uses the non-constant operator[] for iteration). + * So far, set_at() is the closest thing to it. + * + * + * Non-supported usage + * ------------------- + * + * ~~~ + * sv.reserve(20); + * ~~~ + * This has no clear meaning. A usage analogous to STL would precede a sequence + * of push_back's. In that case, we should create a new empty range at the end + * of the vector, and reserve that. Empty ranges are currently not allowed. + * A replacement of this pattern is to create a new `std::vector`, reserve space + * for it and fill it, and finally use `sparse_vector::append()`. If the end + * of the vector is void, there will be no performance penalty, otherwise a + * reserve + copy will happen. + * + * --- + * + * ~~~ + * for (auto& value: sv) { ... } + * ~~~ + * In order to allow for assignment in an item which is currently not void, + * the non-constant iterators do not dereference directly into a reference to + * the vector element (which would be non-existent if in the void), but instead + * into a lightweight object (still called "reference"). These objects are + * semantically working as references, but they are formally rvalues (i.e., just + * values, not C++ references), so they can't be assigned to references (like + * "auto&"). Nevertheless they work as references and assigning to them does + * change the original value. Currently assigning to void cells is not supported + * (see above). + * + */ + template + class sparse_vector { + typedef sparse_vector this_t; + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // - - - public interface + public: + // - - - types + typedef T value_type; ///< type of the stored values + typedef std::vector vector_t; + ///< type of STL vector holding this data + typedef typename vector_t::size_type size_type; ///< size type + typedef typename vector_t::difference_type difference_type; + ///< index difference type + typedef typename vector_t::pointer pointer; + typedef typename vector_t::const_pointer const_pointer; + // typedef typename vector_t::reference reference; + // typedef typename vector_t::const_reference const_reference; + + // --- declarations only --- + class reference; + class const_reference; + class iterator; + class const_iterator; + + class datarange_t; + class const_datarange_t; + // --- ----------------- --- + + typedef std::vector range_list_t; + ///< type of sparse vector data + typedef typename range_list_t::iterator range_iterator; + ///< type of iterator over ranges + typedef typename range_list_t::const_iterator range_const_iterator; + ///< type of constant iterator over ranges + + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // - - - public methods + /// Default constructor: an empty vector + sparse_vector() : nominal_size(0), ranges() {} + + /// Constructor: a vector with new_size elements in the void + sparse_vector(size_type new_size) : nominal_size(0), ranges() { resize(new_size); } + + /** + * @brief Constructor: a solid vector from an existing STL vector + * @param from vector to copy data from + * @param offset (default: 0) index the data starts from (preceeded by void) + */ + sparse_vector(const vector_t& from, size_type offset = 0) : nominal_size(0), ranges() + { + add_range(offset, from.begin(), from.end()); + } + + /// Copy constructor: default + sparse_vector(sparse_vector const&) = default; + + /// Move constructor + sparse_vector(sparse_vector&& from) + : nominal_size(from.nominal_size), ranges(std::move(from.ranges)) + { + from.nominal_size = 0; + } + + /// Copy assignment: default + sparse_vector& operator=(sparse_vector const&) = default; + + /// Move assignment + sparse_vector& operator=(sparse_vector&& from) + { + ranges = std::move(from.ranges); + nominal_size = from.nominal_size; + from.nominal_size = 0; + return *this; + } // operator= (sparse_vector&&) + + /** + * @brief Constructor: a solid vector from an existing STL vector + * @param from vector to move data from + * @param offset (default: 0) index the data starts from (preceeded by void) + */ + sparse_vector(vector_t&& from, size_type offset = 0) : nominal_size(0), ranges() + { + add_range(offset, std::move(from)); + } + + /// Destructor: default + ~sparse_vector() = default; + + // - - - STL-like interface + /// Removes all the data, making the vector empty + void clear() + { + ranges.clear(); + nominal_size = 0; + } + + /// Returns the size of the vector + size_type size() const { return nominal_size; } + + /// Returns whether the vector is empty + bool empty() const { return size() == 0; } + + /// Returns the capacity of the vector (compatibility only) + size_type capacity() const { return nominal_size; } + + //@{ + /// Resizes the vector to the specified size, adding void + void resize(size_type new_size); + /// Resizes the vector to the specified size, adding def_value + void resize(size_type new_size, value_type def_value); + //@} + + //@{ + /// Standard iterators interface + iterator begin(); + iterator end(); + const_iterator begin() const; + const_iterator end() const; + const_iterator cbegin() const { return begin(); } + const_iterator cend() const { return end(); } + //@} + + /// Access to an element (read only) + value_type operator[](size_type index) const; + + /// Access to an element (read/write for non-void elements only!) + reference operator[](size_type index); + + // - - - special interface + + ///@{ @name Cell test + /// + /// The methods in this group test the single vector cells. + + // --- element test + /** + * @brief Returns whether the specified position is void + * @param index position of the cell to be tested + * @throw out_of_range if index is not in the vector + * @see is_back_void() + */ + bool is_void(size_type index) const; + + /// @brief Returns whether the sparse vector ends with void + /// @see is_void() + bool back_is_void() const { return ranges.empty() || (ranges.back().end_index() < size()); } + + /// Returns the number of non-void cells + size_type count() const; + ///@} + + ///@{ @name Cell set + /// + /// The methods in this group access and/or change the cell values. + /** + * @brief Writes into an element (creating or expanding a range if needed) + * @param index the index of the element to set + * @param value the value to be set + * @return a reference to the changed value + * @see unset_at() + * + * Note that setting the value to zero will not cast the element into void. + * Use unset_at for that. + */ + value_type& set_at(size_type index, value_type value); + + /// Casts the element with the specified index into the void + void unset_at(size_type index); + + //@{ + /// Adds one element to the end of the vector (zero values too) + /// @param value value to be added + void push_back(value_type value) { resize(size() + 1, value); } + + /** + * @brief Adds one element to the end of the vector (if zero, just adds void) + * @param value value to be added + * @param thr threshold below which the value is considered zero + * + * If the threshold is strictly negative, all values are pushed back + */ + void push_back(value_type value, value_type thr) + { + if (is_zero(value, thr)) + resize(size() + 1); + else + push_back(value); + } // push_back() + //@} + + //@{ + /** + * @brief Copies data from a sequence between two iterators + * @tparam ITER type of iterator + * @param first iterator pointing to the first element to be copied + * @param last iterator pointing after the last element to be copied + * + * The previous content of the sparse vector is lost. + */ + template + void assign(ITER first, ITER last) + { + clear(); + append(first, last); + } + + /** + * @brief Copies data from a container + * @tparam CONT type of container supporting the standard begin/end interface + * @param new_data container with the data to be copied + * + * The previous content of the sparse vector is lost. + */ + template + void assign(const CONT& new_data) + { + clear(); + append(new_data); + } + + /** + * @brief Moves data from a vector + * @param new_data vector with the data to be moved + * + * The previous content of the sparse vector is lost. + */ + void assign(vector_t&& new_data) + { + clear(); + append(std::move(new_data)); + } + //@} + + ///@} + + // --- BEGIN Ranges --------------------------------------------------------- + /** + * @name Ranges + * + * A range is a contiguous region of the sparse vector which contains all + * non-void values. + * + * A sparse vector is effectively a sorted collection of ranges. + * This interface allows: + * * iteration through all ranges (read only) + * * random access to a range by index (read only) + * * access to the data of a selected range (read/write) + * * look-up of the range containing a specified index + * (read/write; use write with care!) + * * addition (and more generically, combination with existing data) of a new + * range + * * extension of the sparse vector by appending a range at the end of it + * * remove the data of a range, making it void + */ + /// @{ + + // --- range location + + /// Returns the internal list of non-void ranges + const range_list_t& get_ranges() const { return ranges; } + auto iterate_ranges() -> decltype(auto); + + /// Returns the internal list of non-void ranges + size_type n_ranges() const { return ranges.size(); } + + /// Returns the i-th non-void range (zero-based) + const datarange_t& range(size_t i) const { return ranges[i]; } + + //@{ + /** + * @brief Provides direct access to data of i-th non-void range (zero-based) + * @param i index of the range + * @return an object suitable for ranged-for iteration + * + * No information about the positioning of the range itself is provided, + * which can be obtained with other means (e.g. `range(i).begin_index()`). + * The returned object can be used in a ranged-for loop: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * for (std::size_t iRange = 0; iRange < sv.n_ranges(); ++iRange) { + * for (auto& value: sv.range_data(iRange)) { + * v *= 2.0; + * } + * } + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * (with `sv` a `lar::sparse_vector` instance). + * + * While this is a somehow clumsier interface than `get_ranges()`, it allows, + * using the non-`const` version, write access to the data elements. + * It intentionally provides no write access to the location of each range, + * though. + */ + auto range_data(std::size_t i); + auto range_data(std::size_t const i) const { return range_const_data(i); } + //@} + + /// Like `range_data()` but with explicitly read-only access to data. + auto range_const_data(std::size_t i) const; + + /// Returns a constant iterator to the first data range + range_const_iterator begin_range() const { return ranges.begin(); } + + /// Returns a constant iterator to after the last data range + range_const_iterator end_range() const { return ranges.end(); } + + //@{ + /** + * @brief Returns an iterator to the range containing the specified index + * @param index absolute index of the element to be sought + * @return iterator to containing range, or get_ranges().end() if in void + * @throw std::out_of_range if index is not in the vector + * @see is_void() + */ + range_const_iterator find_range_iterator(size_type index) const; + range_iterator find_range_iterator(size_type index) + { + return ranges.begin() + find_range_number(index); + } + //@} + + /** + * @brief Returns the number (0-based) of range containing `index`. + * @param index absolute index of the element to be sought + * @return index of containing range, or `n_ranges()` if in void + * @throw std::out_of_range if index is not in the vector + * @see is_void() + */ + std::size_t find_range_number(size_type index) const + { + return find_range_iterator(index) - begin_range(); + } + + //@{ + /** + * @brief Returns the range containing the specified index + * @param index absolute index of the element to be sought + * @return the containing range + * @throw std::out_of_range if index is in no range (how appropriate!) + * @see is_void() + */ + const datarange_t& find_range(size_type index) const; + datarange_t& find_range(size_type index); + //@} + + /** + * @brief Casts the whole range with the specified item into the void + * @param index absolute index of the element whose range is cast to void + * @return the range just voided + * @throw std::out_of_range if index is not in the vector + * @see unset(), make_void() + */ + datarange_t make_void_around(size_type index); + + //@{ + /** + * @brief Adds a sequence of elements as a range with specified offset + * @tparam ITER type of iterator + * @param offset where to add the elements + * @param first iterator to the first element to be added + * @param last iterator after the last element to be added + * @return the range where the new data was added + * + * If the offset is beyond the current end of the sparse vector, void is + * added before the new range. + * + * Existing ranges can be merged with the new data if they overlap. + */ + template + const datarange_t& add_range(size_type offset, ITER first, ITER last); + + /** + * @brief Copies the elements in container to a range with specified offset + * @tparam CONT type of container supporting the standard begin/end interface + * @param offset where to add the elements + * @param new_data container holding the data to be copied + * @return the range where the new data was added + * + * If the offset is beyond the current end of the sparse vector, void is + * added before the new range. + * + * Existing ranges can be merged with the new data if they overlap. + */ + template + const datarange_t& add_range(size_type offset, const CONT& new_data) + { + return add_range(offset, new_data.begin(), new_data.end()); + } + + /** + * @brief Adds a sequence of elements as a range with specified offset + * @param offset where to add the elements + * @param new_data container holding the data to be moved + * @return the range where the new data was added + * + * If the offset is beyond the current end of the sparse vector, void is + * added before the new range. + * + * Existing ranges can be merged with the new data if they overlap. + * If no merging happens, new_data vector is directly used as the new range + * added; otherwise, it is just copied. + */ + const datarange_t& add_range(size_type offset, vector_t&& new_data); + //@} + + /// @{ + /** + * @brief Combines a sequence of elements as a range with data at `offset`. + * @tparam ITER type of iterator + * @tparam OP combination operation + * @param offset where to add the elements + * @param first iterator to the first element to be added + * @param last iterator after the last element to be added + * @param op operation to be executed element by element + * @param void_value (default: `value_zero`) the value to use for void cells + * @return the range where the new data landed + * + * This is a more generic version of `add_range()`, where instead of + * replacing the target data with the data in [ `first`, `last` [, the + * existing data is combined with the one in that interval. + * The operation `op` is a binary operation with signature equivalent to + * `Data_t op(Data_t, Data_t)`, and the operation is equivalent to + * `v[i + offset] = op(v[i + offset], *(first + offset))`: `op` is a binary + * operation whose first operand is the existing value and the second one + * is the one being provided. + * If the cell `i + offset` is currently void, it will be created and the + * value used in the combination will be `void_value`. + * + * If the offset is beyond the current end of the sparse vector, void is + * added before the new range. + * + * Existing ranges can be merged with the new data if they overlap. + */ + template + const datarange_t& combine_range(size_type offset, + ITER first, + ITER last, + OP&& op, + value_type void_value = value_zero); + + /** + * @brief Combines the elements in container with the data at `offset`. + * @tparam CONT type of container supporting the standard begin/end interface + * @tparam OP combination operation + * @param offset where to add the elements + * @param other container holding the data to be combined + * @param op operation to be executed element by element + * @param void_value (default: `value_zero`) the value to use for void cells + * @return the range where the new data was added + * @see `combine_range()` + * + * This is equivalent to `combine_range(size_type, ITER, ITER, OP, Data_t)` + * using as the new data range the full content of `other` container. + */ + template + const datarange_t& combine_range(size_type offset, + const CONT& other, + OP&& op, + value_type void_value = value_zero) + { + return combine_range(offset, other.begin(), other.end(), std::forward(op), void_value); + } + /// @} + + //@{ + /** + * @brief Adds a sequence of elements as a range at the end of the vector. + * @param first iterator to the first element to be added + * @param last iterator after the last element to be added + * @return the range where the new data was added + * + * The input range is copied at the end of the sparse vector. + * If the end of the sparse vector was the end of a range, that range is + * expanded, otherwise a new one is created. + */ + template + const datarange_t& append(ITER first, ITER last) + { + return add_range(size(), first, last); + } + + /** + * @brief Adds a sequence of elements as a range at the end of the vector. + * @param range_data contained holding the data to be copied or moved + * @return the range where the new data was added + * + * The input data is copied at the end of the sparse vector. + * If the end of the sparse vector was the end of a range, that range is + * expanded, otherwise a new one is created. + */ + template + const datarange_t& append(const CONT& range_data) + { + return add_range(size(), range_data); + } + + /** + * @brief Adds a sequence of elements as a range at the end of the vector. + * @param range_data contained holding the data to be copied or moved + * @return the range where the new data was added + * + * If there is a range at the end of the sparse vector, it will be expanded + * with the new data. + * Otherwise, this method will use the data vector directly as a new range + * added at the end of the sparse vector. + */ + const datarange_t& append(vector_t&& range_data) + { + return add_range(size(), std::move(range_data)); + } + //@} + + //@{ + /** + * @brief Turns the specified range into void + * @param iRange iterator or index of range to be deleted + * @return the range just voided + * @see make_void(), unset_at() + * + * The range is effectively removed from the sparse vector, rendering void + * the interval it previously covered. + * The range object itself is returned (no copy is performed). + * + * The specified range must be valid. Trying to void an invalid range + * (including `end_range()`) yields undefined behavior. + */ + datarange_t void_range(range_iterator const iRange); + datarange_t void_range(std::size_t const iRange) { return void_range(ranges.begin() + iRange); } + //@} + ///@} + + /** + * @brief Returns if the vector is in a valid state + * + * The vector is in a valid state if: + * - no ranges overlap or touch each other (a void gap must exist) + * - no range is empty + * - all ranges are sorted + * - the size of the vector is not smaller than the sum of the size of + * the ranges plus the internal gaps + * An invalid state can be the result of “too smart” use of this + * class, or of a bug, which should be reported to the author. + */ + bool is_valid() const; + + /// Makes all the elements from first and before last void + void make_void(iterator first, iterator last); + + //@{ + /// Performs internal optimization, returns whether the object was changed + bool optimize() { return optimize(min_gap()); } + bool optimize(size_t) { return false; /* no optimization implemented yet */ } + //@} + + ///@{ @name Static members for dealing with this type of value + + static constexpr value_type value_zero{0}; ///< a representation of 0 + + /// Returns the module of the specified value + static value_type abs(value_type v) { return (v < value_zero) ? -v : v; } + + /// Returns whether the value is exactly zero + static value_type is_zero(value_type v) { return v == value_zero; } + + /// Returns whether the value is zero below a given threshold + static value_type is_zero(value_type v, value_type thr) { return abs(v - value_zero) <= thr; } + + /// Returns whether two values are the same + static value_type is_equal(value_type a, value_type b) { return is_zero(abs(a - b)); } + + /// Returns whether two values are the same below a given threshold + static value_type is_equal(value_type a, value_type b, value_type thr) + { + return is_zero(abs(a - b), thr); + } + ///@} + + ///@{ @name Static members related to data size and optimization + + /// Returns the expected size taken by a vector of specified size + static size_t expected_vector_size(size_t size); + + /// Minimum optimal gap between ranges (a guess) + static size_t min_gap(); + + /// Returns if merging the two specified ranges would save memory + static bool should_merge(const typename datarange_t::base_t& a, + const typename datarange_t::base_t& b); + ///@} + + // -------------------------------------------------------------------------- + protected: + size_type nominal_size; ///< current size + range_list_t ranges; ///< list of ranges + + //@{ + /** + * @brief Returns an iterator to the range after `index`. + * @param index the absolute index + * @return iterator to the next range not including index, or ranges.end() + * if none + */ + range_iterator find_next_range_iter(size_type index) + { + return find_next_range_iter(index, ranges.begin()); + } + range_const_iterator find_next_range_iter(size_type index) const + { + return find_next_range_iter(index, ranges.cbegin()); + } + //@} + //@{ + /** + * @brief Returns an iterator to the range after `index`, + * or `ranges.end()` if none. + * @param index the absolute index + * @param rbegin consider only from this range on + * @return iterator to the next range not including index, or `ranges.end()` + * if none + */ + range_iterator find_next_range_iter(size_type index, range_iterator rbegin); + range_const_iterator find_next_range_iter(size_type index, range_const_iterator rbegin) const; + //@} + + //@{ + /** + * @brief Returns an iterator to the range at or after `index`. + * @param index the absolute index + * @return iterator to the next range including index or after it, + * or ranges.end() if none + * @see `find_next_range_iter()` + * + * If `index` is in a range, an iterator to that range is returned. + * If `index` is in the void, an iterator to the next range after `index` is + * returned. If there is no such range after `index`, `ranges.end()` is + * returned. + */ + range_iterator find_range_iter_at_or_after(size_type index); + range_const_iterator find_range_iter_at_or_after(size_type index) const; + //@} + + //@{ + /** + * @brief Returns an iterator to the range no earlier than `index`, + or `end()` if none. + * @param index the absolute index + * @return iterator to the range including index, or the next range if none + * + * The returned iterator points to a range that "borders" the specified index, + * meaning that the cell at `index` is either within the range, or it is the + * one immediately after that range. + * If `index` is in the middle of the void, though (i.e. if the previous cell + * is void), the next range is returned instead. + * Finally, if there is no next range, `end_range()` is returned. + * + * The result may be also interpreted as follow: if the start of the returned + * range is lower than `index`, then the cell at `index` belongs to that + * range. Otherwise, it initiates its own range (but that range might end up + * being contiguous to the next(. + */ + range_iterator find_extending_range_iter(size_type index) + { + return find_extending_range_iter(index, ranges.begin()); + } + range_const_iterator find_extending_range_iter(size_type index) const + { + return find_extending_range_iter(index, ranges.cbegin()); + } + //@} + + //@{ + /** + * @brief Returns an iterator to the range that contains the first non-void + * element after `index`, or `end()` if none. + * @param index the absolute index + * @param rbegin consider only from this range on + * @return iterator to the next range not including index, or `ranges.end()` + * if none + * + * Note that the returned range can contain `index` as well. + */ + range_iterator find_extending_range_iter(size_type index, range_iterator rbegin); + range_const_iterator find_extending_range_iter(size_type index, + range_const_iterator rbegin) const; + //@} + + /// Returns the size determined by the ranges already present + size_type minimum_size() const { return ranges.empty() ? 0 : ranges.back().end_index(); } + + //@{ + /// Plug a new data range in the specified position; no check performed + range_iterator insert_range(range_iterator iInsert, const datarange_t& data) + { + return data.empty() ? iInsert : ranges.insert(iInsert, data); + } + range_iterator insert_range(range_iterator iInsert, datarange_t&& data) + { + return data.empty() ? iInsert : ranges.insert(iInsert, std::move(data)); + } + //@} + + /// Implementation detail of `add_range()`, with where to add the range. + const datarange_t& add_range_before(size_type offset, + vector_t&& new_data, + range_iterator nextRange); + + /// Voids the starting elements up to index (excluded) of a given range + range_iterator eat_range_head(range_iterator iRange, size_t index); + + /** + * @brief Merges all the following contiguous ranges. + * @param iRange iterator to the range to merge the following ones into + * @return iterator to the merged range + * + * Starting from the range next to `iRange`, if that range is contiguous to + * `iRange` the two are merged. The merging continues while there are + * ranges contiguous to `iRange`. In the end, an iterator is returned pointing + * to a range that has the same starting point that `iRange` had, and that is + * not followed by a contiuguous range, since all contiguous ranges have been + * merged into it. + */ + datarange_t& merge_ranges(range_iterator iRange); + + /// Extends the vector size according to the last range + size_type fix_size(); + + }; // class sparse_vector<> + +} // namespace lar + +/** + * @brief Prints a sparse vector into a stream + * @tparam T template type of the sparse vector + * @param out output stream + * @param v the sparse vector to be written + * @return the output stream (out) + * + * The output is in the form: + * 
+ * Sparse vector of size ## with ## ranges:
+ *   [min1 - max1] (size1) { elements of the first range }
+ *   [min2 - max2] (size2) { elements of the second range }
+ *   ...
+ *
+ */ +template +std::ostream& operator<<(std::ostream& out, const lar::sparse_vector& v); + +// ----------------------------------------------------------------------------- +// --- sparse_vector::datarange_t definition +// --- + +/// Range class, with range and data +template +class lar::sparse_vector::datarange_t : public range_t { +public: + typedef range_t base_t; ///< base class + + typedef typename vector_t::iterator iterator; + typedef typename vector_t::const_iterator const_iterator; + + /// Default constructor: an empty range + datarange_t() : base_t(), values() {} + + /// Constructor: range initialized with 0 + datarange_t(const base_t& range) : base_t(range), values(range.size()) {} + + /// Constructor: offset and data + template + datarange_t(size_type offset, ITER first, ITER last) + : base_t(offset, offset + std::distance(first, last)), values(first, last) + {} + + /// Constructor: offset and data as a vector (which will be used directly) + datarange_t(size_type offset, vector_t&& data) + : base_t(offset, offset + data.size()), values(data) + {} + + //@{ + /// Returns an iterator to the specified absolute value (no check!) + iterator get_iterator(size_type index) { return values.begin() + index - base_t::begin_index(); } + const_iterator get_iterator(size_type index) const { return get_const_iterator(index); } + const_iterator get_const_iterator(size_type index) const + { + return values.begin() + index - base_t::begin_index(); + } + //@} + + //@{ + /// begin and end iterators + iterator begin() { return values.begin(); } + iterator end() { return values.end(); } + const_iterator begin() const { return values.begin(); } + const_iterator end() const { return values.end(); } + const_iterator cbegin() const { return values.cbegin(); } + const_iterator cend() const { return values.cend(); } + //@} + + //@{ + /// Resizes the range (optionally filling the new elements with def_value) + void resize(size_t new_size) + { + values.resize(new_size); + fit_size_from_data(); + } + void resize(size_t new_size, value_type def_value) + { + values.resize(new_size, def_value); + fit_size_from_data(); + } + //@} + + //@{ + /// Returns the value at the specified absolute index + value_type& operator[](size_type index) { return values[base_t::relative_index(index)]; } + const value_type& operator[](size_type index) const + { + return values[base_t::relative_index(index)]; + } + //@} + + //@{ + /// Return the vector of data values + const vector_t& data() const { return values; } + // vector_t& data() { return values; } + //@} + + /** + * @brief Appends the specified elements to this range. + * @tparam ITER type of iterator of the range + * @param index the starting point + * @param first iterator to the first object to copy + * @param last iterator after the last object to copy + * @return the extended range + */ + template + datarange_t& extend(size_type index, ITER first, ITER last); + + /** + * @brief Moves the begin of this range + * @param to_index absolute index to move the head to + * @param def_value value to be inserted in case of expansion of the range + */ + void move_head(size_type to_index, value_type def_value = value_zero); + + /** + * @brief Moves the end of this range + * @param to_index absolute index to move the tail to + * @param def_value value to be inserted in case of expansion of the range + */ + void move_tail(size_type to_index, value_type def_value = value_zero) + { + resize(base_t::relative_index(to_index), def_value); + } + + /** + * @brief Dumps the content of this data range into a stream. + * @tparam Stream type of stream to sent the dump to + * @param out stream to sent the dump to + * + * The output format is: + * + * [min - max] (size) { values... } + * + * Output is on a single line, which is not terminated. + */ + template + void dump(Stream&& out) const; + +protected: + vector_t values; ///< data in the range + + void fit_size_from_data() { base_t::resize(values.size()); } + +}; // class datarange_t + +/** + * @brief A constant reference to a data range. + * + * Values in the range can be modified, but their position and number can not. + */ +template +class lar::sparse_vector::const_datarange_t : private datarange_t { +public: + using iterator = typename datarange_t::iterator; + using const_iterator = typename datarange_t::const_iterator; + + // `range_t` interface + using datarange_t::begin_index; + using datarange_t::borders; + using datarange_t::empty; + using datarange_t::end_index; + using datarange_t::includes; + using datarange_t::last; + using datarange_t::offset; + using datarange_t::overlap; + using datarange_t::relative_index; + using datarange_t::separate; + using datarange_t::size; + using datarange_t::operator<; + using datarange_t::operator==; + using datarange_t::is_valid; + using datarange_t::less; + using less_int_range = typename datarange_t::less_int_range; + + // `datarange_t` interface + using datarange_t::get_iterator; + decltype(auto) begin() { return datarange_t::begin(); } + decltype(auto) end() { return datarange_t::end(); } + // using datarange_t::begin; + // using datarange_t::end; + // using datarange_t::cbegin; + // using datarange_t::cend; + using datarange_t::operator[]; + using datarange_t::dump; + + /// Return the vector of data values (only constant access). + const vector_t& data() const { return base().values; } + + ~const_datarange_t() = delete; // can't destroy; can cast into it though + +protected: + friend decltype(auto) details::make_const_datarange_t(datarange_t&); + + datarange_t const& base() const { return static_cast(*this); } + datarange_t& base() { return static_cast(*this); } + + static void static_check() { static_assert(sizeof(const_datarange_t) == sizeof(datarange_t)); } + +}; // lar::sparse_vector::const_datarange_t + +// ----------------------------------------------------------------------------- +// --- sparse_vector iterators definition +// --- + +/// Special little box to allow void elements to be treated as references. +template +class lar::sparse_vector::const_reference { +protected: + const value_type* ptr; + +public: + const_reference(const value_type* pValue = 0) : ptr(pValue) {} + const_reference(const value_type& value) : const_reference(&value) {} + + explicit operator value_type() const { return ptr ? *ptr : value_zero; } + operator const value_type&() const { return ptr ? *ptr : value_zero; } +}; // lar::sparse_vector::const_reference + +/** + * @brief A class representing a cell in a sparse vector. + * + * This class is a little box allowing assignment of values into it; + * if the internal pointer is invalid (as in case of void cell), + * dereferencing or assigning will provoke a segmentation fault. + */ +template +class lar::sparse_vector::reference : public const_reference { + friend class iterator; + + // This object "disappears" when assigned to: either the assignment is not + // possible, and then a segmentation fault will occur, or the return value + // is an actual C++ reference to the assigned value. + // The same is true when explicitly converting it to a reference. +public: + reference(value_type* pValue = 0) : const_reference(pValue) {} + reference(value_type& value) : reference(&value) {} + + reference& operator=(const reference&) = default; + value_type& operator=(value_type v) { return const_cast(*const_reference::ptr) = v; } + + explicit operator value_type&() { return const_cast(*const_reference::ptr); } + +protected: + explicit reference(const const_reference& from) : const_reference(from) {} + +}; // lar::sparse_vector::reference + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +/// Iterator to the sparse vector values +template +class lar::sparse_vector::const_iterator { + // + // This iterator fulfils the traits of an immutable forward iterator. + // + +protected: + friend class container_t; + + typedef sparse_vector container_t; + typedef typename container_t::size_type size_type; + typedef typename container_t::range_list_t::const_iterator ranges_const_iterator; + +public: + // so far, only forward interface is implemented; + // backward is tricky... + typedef std::forward_iterator_tag iterator_category; + + /// Namespace for special initialization + struct special { + class begin {}; + class end {}; + }; + + // types + typedef typename container_t::value_type value_type; + typedef typename container_t::difference_type difference_type; + typedef typename container_t::pointer pointer; + typedef typename container_t::reference reference; + // note that this is not out little box, it's the normal C++ constant + // reference from the underlying vector + typedef typename vector_t::const_reference const_reference; + + /// Default constructor, does not iterate anywhere + const_iterator() : cont(nullptr), index(0), currentRange() {} + + /// Constructor from a container and a offset + const_iterator(const container_t& c, size_type offset) + : cont(&c), index(std::min(offset, c.size())), currentRange() + { + refresh_state(); + } + + /// Special constructor: initializes at the beginning of the container + const_iterator(const container_t& c, const typename special::begin) + : cont(&c), index(0), currentRange(c.get_ranges().begin()) + {} + + /// Special constructor: initializes at the end of the container + const_iterator(const container_t& c, const typename special::end) + : cont(&c), index(c.size()), currentRange(c.get_ranges().end()) + {} + + /// Random access + const_reference operator[](size_type offset) const { return (*cont)[index + offset]; } + + /// Constant dereferenciation operator + const_reference operator*() const; + + //@{ + /// Increment and decrement operators + const_iterator& operator++(); + const_iterator operator++(int) + { + const_iterator copy(*this); + const_iterator::operator++(); + return copy; + } + //@} + + //@{ + /// Increment and decrement operators + const_iterator& operator+=(difference_type delta); + const_iterator& operator-=(difference_type delta); + const_iterator operator+(difference_type delta) const; + const_iterator operator-(difference_type delta) const; + //@} + + /// Distance operator + difference_type operator-(const const_iterator& iter) const; + + //@{ + /// Iterator comparisons + bool operator==(const const_iterator& as) const + { + return (cont == as.cont) && (index == as.index); + } + bool operator!=(const const_iterator& as) const + { + return (cont != as.cont) || (index != as.index); + } + bool operator<(const const_iterator& than) const + { + return (cont == than.cont) && (index < than.index); + } + bool operator>(const const_iterator& than) const + { + return (cont == than.cont) && (index > than.index); + } + bool operator<=(const const_iterator& than) const + { + return (cont == than.cont) && (index <= than.index); + } + bool operator>=(const const_iterator& than) const + { + return (cont == than.cont) && (index >= than.index); + } + //@} + + /// Returns the current range internal value; use it at your own risk!! + range_const_iterator get_current_range() const { return currentRange; } + +protected: + // + // Logic of the internals: + // - cont provides the list of ranges + // - index is the absolute index in the container + // - currentRange is a pointer to the "current" range, cached for speed + // + // An iterator past-the-end has the index equal to the size of the + // container it points to. Operations on it are undefined, except that + // going backward by one decrement goes back to the container + // (if not empty) TODO currently backward iteration is not supported + // + // When the iterator is pointing to an element within a range, + // currentRange points to that range. + // When the iterator is pointing into the void of the vector, currentRange + // points to the range next to that void area, or end() if there is none. + // When the iterator is past-the-end, currentRange is not defined. + // + // Conversely, when the index is equal (or greater) to the size of the + // container, the iterator is not dereferenciable, the iterator points + // past-the-end of the container. + // When currentRange points to a valid range and the index + // is before that range, the iterator is pointing to the void. + // When currentRange points to a valid range and the index + // is inside that range, the iterator is pointing to an item in that + // range. + // When currentRange points to a valid range, the index can't point beyond + // that range. + // When currentRange points to the past-to-last range, the iterator is + // pointing to the void (past the last range). + // + const container_t* cont; ///< pointer to the container + size_type index; ///< pointer to the current value, as absolute index + ranges_const_iterator currentRange; ///< pointer to the current (or next) range + + /// Reassigns the internal state according to the index + void refresh_state(); + +}; // class lar::sparse_vector::const_iterator + +/** + * @brief Iterator to the sparse vector values + * + * This iterator respects the traits of an immutable forward iterator, + * EXCEPT that the iterator can be non-dereferenciable even when it's a + * "past the end" iterator. + * That is due to the fact that currently dereferencing (and assigning) + * to a cell which is not in a range already is not supported yet + * (it can be done with some complicate mechanism). + */ +template +class lar::sparse_vector::iterator : public const_iterator { + typedef typename const_iterator::container_t container_t; + friend typename const_iterator::container_t; + +public: + typedef typename const_iterator::reference reference; + typedef typename const_iterator::const_reference const_reference; + typedef typename const_iterator::special special; + + /// Default constructor, does not iterate anywhere + iterator() : const_iterator() {} + + /// Constructor from a container and an offset + iterator(container_t& c, size_type offset = 0) : const_iterator(c, offset) {} + + /// Special constructor: initializes at the beginning of the container + iterator(const container_t& c, const typename special::begin _) : const_iterator(c, _) {} + + /// Special constructor: initializes at the end of the container + iterator(const container_t& c, const typename special::end _) : const_iterator(c, _) {} + + /// Random access + reference operator[](size_type offset) const + { + return (*const_iterator::cont)[const_iterator::index + offset]; + } + + //@{ + /// Increment and decrement operators + iterator& operator++() + { + const_iterator::operator++(); + return *this; + } + iterator operator++(int _) { return const_iterator::operator++(_); } + //@} + + //@{ + /// Increment and decrement operators + iterator& operator+=(difference_type delta) + { + const_iterator::operator+=(delta); + return *this; + } + iterator& operator-=(difference_type delta) + { + const_iterator::operator+=(delta); + return *this; + } + iterator operator+(difference_type delta) const { return const_iterator::operator+(delta); } + iterator operator-(difference_type delta) const { return const_iterator::operator-(delta); } + //@} + + /// Dereferenciation operator (can't write non-empty elements!) + reference operator*() const { return reference(const_iterator::operator*()); } + + // /// Returns the current range internal value; use it at your own risk!! + // range_iterator get_current_range() const + // { return const_iterator::currentRange; } + +protected: + // also worth conversion + iterator(const_iterator from) : const_iterator(from) {} + +}; // class lar::sparse_vector::iterator + +// ----------------------------------------------------------------------------- +// --- implementation -------------------------------------------------------- +// ----------------------------------------------------------------------------- +namespace lar::details { + + // -------------------------------------------------------------------------- + /// Enclosure to use two iterators representing a range in a range-for loop. + template + class iteratorRange { + BITER b; + EITER e; + + public: + iteratorRange(BITER const& b, EITER const& e) : b(b), e(e) {} + auto const& begin() const { return b; } + auto const& end() const { return e; } + auto const& size() const { return std::distance(begin(), end()); } + }; // iteratorRange() + + // -------------------------------------------------------------------------- + template + decltype(auto) make_const_datarange_t(typename sparse_vector::datarange_t& r) + { + return static_cast::const_datarange_t&>(r); + } + + // -------------------------------------------------------------------------- + template + class const_datarange_iterator { + using const_datarange_t = typename sparse_vector::const_datarange_t; + using base_iterator = typename sparse_vector::range_iterator; + base_iterator it; + + public: + // minimal set of features for ranged-for loops + const_datarange_iterator() = default; + const_datarange_iterator(base_iterator it) : it(it) {} + + const_datarange_iterator& operator++() + { + ++it; + return *this; + } + const_datarange_t& operator*() const { return make_const_datarange_t(*it); } + bool operator!=(const_datarange_iterator const& other) const { return it != other.it; } + }; + + // -------------------------------------------------------------------------- + +} // namespace lar::details + +//------------------------------------------------------------------------------ +//--- sparse_vector implementation + +template +constexpr typename lar::sparse_vector::value_type lar::sparse_vector::value_zero; + +template +decltype(auto) lar::sparse_vector::iterate_ranges() +{ + return details::iteratorRange(details::const_datarange_iterator(ranges.begin()), + details::const_datarange_iterator(ranges.end())); +} // lar::sparse_vector::iterate_ranges() + +template +void lar::sparse_vector::resize(size_type new_size) +{ + if (new_size >= size()) { + nominal_size = new_size; + return; + } + + // truncating... + range_iterator iLastRange = find_next_range_iter(new_size); + ranges.erase(iLastRange, ranges.end()); + if (!ranges.empty()) { + range_iterator iLastRange = ranges.end() - 1; + if (new_size == iLastRange->begin_index()) + ranges.erase(iLastRange); + else if (new_size < iLastRange->end_index()) + iLastRange->resize(new_size - iLastRange->begin_index()); + } // if we have ranges + + // formally resize + nominal_size = new_size; +} // lar::sparse_vector::resize() + +template +void lar::sparse_vector::resize(size_type new_size, value_type def_value) +{ + if (new_size == size()) return; + if (new_size > size()) { + + if (back_is_void()) // add a new range + append(vector_t(new_size - size(), def_value)); + else // extend the last range + ranges.back().resize(new_size - ranges.back().begin_index(), def_value); + + nominal_size = new_size; + return; + } + // truncating is the same whether there is a default value or not + resize(new_size); +} // lar::sparse_vector::resize() + +template +inline typename lar::sparse_vector::iterator lar::sparse_vector::begin() +{ + return iterator(*this, typename iterator::special::begin()); +} + +template +inline typename lar::sparse_vector::iterator lar::sparse_vector::end() +{ + return iterator(*this, typename iterator::special::end()); +} + +template +inline typename lar::sparse_vector::const_iterator lar::sparse_vector::begin() const +{ + return const_iterator(*this, typename const_iterator::special::begin()); +} + +template +inline typename lar::sparse_vector::const_iterator lar::sparse_vector::end() const +{ + return const_iterator(*this, typename const_iterator::special::end()); +} + +template +typename lar::sparse_vector::value_type lar::sparse_vector::operator[](size_type index) const +{ + // first range not including the index + range_const_iterator iNextRange = find_next_range_iter(index); + + // if not even the first range includes the index, we are in the void + // (or in some negative index space, if index signedness allowed); + if (iNextRange == ranges.begin()) return value_zero; + + // otherwise, let's take the previous range; + // if it includes the index, we return its value; + // or it precedes it, and our index is in the void: return zero + const datarange_t& range(*--iNextRange); + return (index < range.end_index()) ? range[index] : value_zero; +} // lar::sparse_vector::operator[] + +template +typename lar::sparse_vector::reference lar::sparse_vector::operator[](size_type index) +{ + // first range not including the index + range_iterator iNextRange = find_next_range_iter(index); + + // if not even the first range includes the index, we are in the void + // (or in some negative index space, if index signedness allowed); + if (iNextRange == ranges.begin()) return reference(); + + // otherwise, let's take the previous range; + // if it includes the index, we return its value; + // or it precedes it, and our index is in the void: return zero + datarange_t& range(*--iNextRange); + return (index < range.end_index()) ? reference(range[index]) : reference(); +} // lar::sparse_vector::operator[] + +template +bool lar::sparse_vector::is_void(size_type index) const +{ + if (ranges.empty() || (index >= size())) throw std::out_of_range("empty sparse vector"); + // range after the index: + range_const_iterator iNextRange = find_next_range_iter(index); + return ((iNextRange == ranges.begin()) || ((--iNextRange)->end_index() <= index)); +} // lar::sparse_vector::is_void() + +template +inline typename lar::sparse_vector::size_type lar::sparse_vector::count() const +{ + return std::accumulate(begin_range(), + end_range(), + size_type(0), + [](size_type s, const datarange_t& rng) { return s + rng.size(); }); +} // count() + +template +typename lar::sparse_vector::value_type& lar::sparse_vector::set_at(size_type const index, + value_type value) +{ + // first range not including the index + range_iterator iNextRange = find_next_range_iter(index); + + // if not even the first range includes the index, we are in the void + // (or in some negative index space, if index signedness allowed); + if (iNextRange != ranges.begin()) { + // otherwise, let's take the previous range; + // if it includes the index, we set the existing value; + // or it precedes it, and our index is in the void, add the value as a + // range + datarange_t& range(*std::prev(iNextRange)); + if (index < range.end_index()) return range[index] = value; + } + // so we are in the void; add the value as a new range + return const_cast(add_range(index, {value}))[index]; +} // lar::sparse_vector::set_at() + +template +void lar::sparse_vector::unset_at(size_type index) +{ + // first range not including the index + range_iterator iNextRange = find_next_range_iter(index); + + // if not even the first range includes the index, we are in the void + // (or in some negative index space, if index signedness allowed); + if (iNextRange == ranges.begin()) return; + + // otherwise, let's take the previous range; + // or it precedes the index, and our index is in the void + datarange_t& range(*--iNextRange); + if (index >= range.end_index()) return; // void already + + // it includes the index: + // - if it's a one-element range, remove the range + if (range.size() == 1) ranges.erase(iNextRange); + // - if it's on a border, resize the range + else if (index == range.begin_index()) + range.move_head(index + 1); + else if (index == range.end_index() - 1) + range.move_tail(index); + // - if it's inside, break the range in two + else if (index < range.end_index()) { + // we are going to put a hole in a range; + // this effectively creates two ranges; + // first create the rightmost range, and add it to the list + ranges.emplace( + ++iNextRange, index + 1, range.begin() + range.relative_index(index + 1), range.end()); + // then cut the existing one + range.move_tail(index); + } +} // lar::sparse_vector::unset_at() + +template +auto lar::sparse_vector::range_data(std::size_t const i) +{ + auto& r = ranges[i]; + return details::iteratorRange(r.begin(), r.end()); +} + +template +auto lar::sparse_vector::range_const_data(std::size_t const i) const +{ + auto& r = ranges[i]; + return details::iteratorRange(r.cbegin(), r.cend()); +} + +template +typename lar::sparse_vector::range_const_iterator lar::sparse_vector::find_range_iterator( + size_type index) const +{ + if (ranges.empty()) throw std::out_of_range("empty sparse vector"); + // range after the index: + range_const_iterator iNextRange = find_next_range_iter(index); + return ((iNextRange == ranges.begin()) || (index >= (--iNextRange)->end_index())) ? ranges.end() : + iNextRange; +} // lar::sparse_vector::find_range_iterator() const + +template +const typename lar::sparse_vector::datarange_t& lar::sparse_vector::find_range( + size_type index) const +{ + if (ranges.empty()) throw std::out_of_range("empty sparse vector"); + // range on the index: + range_const_iterator iNextRange = find_range_iterator(index); + if (iNextRange == ranges.end()) throw std::out_of_range("index in no range of the sparse vector"); + return *iNextRange; +} // lar::sparse_vector::find_range() const + +template +inline typename lar::sparse_vector::datarange_t& lar::sparse_vector::find_range( + size_type index) +{ + return const_cast(const_cast(this)->find_range(index)); +} // lar::sparse_vector::find_range() + +template +auto lar::sparse_vector::make_void_around(size_type index) -> datarange_t +{ + if (ranges.empty() || (index >= size())) throw std::out_of_range("empty sparse vector"); + // range after the index: + range_iterator iNextRange = find_next_range_iter(index); + if ((iNextRange == ranges.begin()) || ((--iNextRange)->end_index() <= index)) { return {}; } + return void_range(iNextRange); +} // lar::sparse_vector::make_void_around() + +template +template +const typename lar::sparse_vector::datarange_t& +lar::sparse_vector::add_range(size_type offset, ITER first, ITER last) +{ + // insert the new range before the existing range which starts after offset + range_iterator iInsert = find_next_range_iter(offset); + + // is there a range before this, which includes the offset? + if ((iInsert != ranges.begin()) && std::prev(iInsert)->borders(offset)) { + // then we should extend it + (--iInsert)->extend(offset, first, last); + } + else { + // no range before the insertion one includes the offset of the new range; + // ... we need to add it as a new range + iInsert = insert_range(iInsert, {offset, first, last}); + } + return merge_ranges(iInsert); +} // lar::sparse_vector::add_range() + +template +const typename lar::sparse_vector::datarange_t& lar::sparse_vector::add_range( + size_type offset, + vector_t&& new_data) +{ + // insert the new range before the existing range which starts after offset + return add_range_before( + offset, + std::move(new_data), + std::upper_bound(ranges.begin(), + ranges.end(), + offset, + typename datarange_t::less_int_range(datarange_t::less))); + +} // lar::sparse_vector::add_range(vector) + +template +template +auto lar::sparse_vector::combine_range(size_type offset, + ITER first, + ITER last, + OP&& op, + value_type void_value /* = value_zero */ + ) -> const datarange_t& +{ + /* + * This is a complicate enough task, that we go brute force: + * 1) combine all the input elements within the datarange where offset falls + * in + * 2) create a new datarange combining void with the remaining input elements + * 3) if the void area is over before the input elements are, repeat steps + * (1) and (2) with the updated offset + * 4) at this point we'll have a train of contiguous ranges, result of + * combination of the input elements alternating with existing elements + * and with void cells: apply the regular merge algorithm + * + */ + + auto src = first; + auto const insertionPoint = offset; // saved for later + auto destRange = find_range_iter_at_or_after(offset); + while (src != last) { + + // + // (1) combine all the input elements within the datarange including offset + // (NOT extending the range) + // + if ((destRange != end_range()) && destRange->includes(offset)) { + // combine input data until this range is over (or input data is over) + auto dest = destRange->get_iterator(offset); + + auto const end = destRange->end(); + while (src != last) { + *dest = op(*dest, *src); + ++src; + ++offset; + if (++dest == end) break; + } // while + if (src == last) break; + offset = destRange->end_index(); + ++destRange; + } // if + + // + // (2) create a new datarange combining void with input elements + // + // at this point, offset is in the void, we do have more input data, + // and `destRange` does _not_ contain the current insertion offset; + // we fill as much void as we can with data, creating a new range. + // When to stop? at the beginning of the next range, or when data is over + // + size_type const newRangeSize = + (destRange == end_range()) ? std::distance(src, last) : (destRange->begin_index() - offset); + + // prepare the data (we'll plug it in directly) + vector_t combinedData; + combinedData.reserve(newRangeSize); + size_type i = 0; + while (i++ < newRangeSize) { + combinedData.push_back(op(void_value, *src)); + if (++src == last) break; // no more data + } + // move the data as a new range inserted before the next range we just found + // return value is the iterator to the inserted range + destRange = insert_range(destRange, {offset, std::move(combinedData)}); + + // + // (3) if there is more input, repeat steps (1) and (2) with updated offset + // + offset = destRange->end_index(); + ++destRange; + } // while + + // + // (4) apply the regular merge algorithm; + // since we did not extend existing ranges, it may happen that we have + // created a new range at `insertionPoint` contiguous to the previous one; + // so we take care of checking that too + // + return merge_ranges(find_extending_range_iter(insertionPoint == 0 ? 0 : insertionPoint - 1)); + +} // lar::sparse_vector::combine_range() + +template +void lar::sparse_vector::make_void(iterator first, iterator last) +{ + // iterators have in "currentRange" either the range they point into, + // or the range next to the void they point to + + if ((first.cont != this) || (last.cont != this)) { + throw std::runtime_error("lar::sparse_vector::make_void(): iterators from alien container"); + } + // if first is after next, no range is identified + if (first >= last) return; + + range_iterator first_range = ranges.begin() + (first.get_current_range() - ranges.begin()), + last_range = ranges.begin() + (last.get_current_range() - ranges.begin()); + + //--- "first": void up to this iterator --- + // if first in the last void region, there is nothing to erase + if (first_range == ranges.end()) return; + + // if first is in the middle of a valid range instead, we have to resize it + if (first.index > first_range->begin_index()) { + if (first_range == last_range) { + // we are going to erase a subset of a range; + // this effectively creates two ranges; + // first create the rightmost (last) range, and add it to the list + last_range = ranges.emplace(++last_range, + last.index, + first_range->begin() + first_range->relative_index(last.index), + first_range->end()); + // then cut the existing one + first_range->move_tail(first.index); + return; + } + first_range->move_tail(first.index); + ++first_range; // from next range on, start voiding + } + + //--- "last" + // if the index is inside a range, remove up to it + if ((last_range != ranges.end()) && (last.index > last_range->begin_index())) + eat_range_head(last_range, last.index); + + // finally, remove entirely the ranges in between + ranges.erase(first_range, last_range); +} // lar::sparse_vector::make_void() + +template +auto lar::sparse_vector::void_range(range_iterator iRange) -> datarange_t +{ + auto r{std::move(*iRange)}; // triggering move constructor + ranges.erase(iRange); // the emptied range is removed from vector + return r; // returning it as a temporary avoids copies +} // lar::sparse_vector::void_range() + +template +bool lar::sparse_vector::is_valid() const +{ + // a sparse vector with no non-null elements can't be detected invalid + if (ranges.empty()) return true; + + range_const_iterator iNext = ranges.begin(), rend = ranges.end(); + while (iNext != rend) { + range_const_iterator iRange = iNext++; + if (iRange->empty()) return false; + if (iNext != rend) { + if (!(*iRange < *iNext)) return false; + if (!iRange->separate(*iNext)) return false; + } + } // while + if (nominal_size < ranges.back().end_index()) return false; + return true; +} // lar::sparse_vector::is_valid() + +// --- private methods + +template +typename lar::sparse_vector::range_iterator lar::sparse_vector::find_next_range_iter( + size_type index, + range_iterator rbegin) +{ + // this range has the offset (first index) above the index argument: + return std::upper_bound( + rbegin, ranges.end(), index, typename datarange_t::less_int_range(datarange_t::less)); +} // lar::sparse_vector::find_next_range_iter() + +template +typename lar::sparse_vector::range_const_iterator lar::sparse_vector::find_next_range_iter( + size_type index, + range_const_iterator rbegin) const +{ + // this range has the offset (first index) above the index argument: + return std::upper_bound( + rbegin, ranges.end(), index, typename datarange_t::less_int_range(datarange_t::less)); +} // lar::sparse_vector::find_next_range_iter() const + +template +typename lar::sparse_vector::range_const_iterator +lar::sparse_vector::find_range_iter_at_or_after(size_type index) const +{ + // this range has the offset (first index) above the index argument: + auto after = find_next_range_iter(index); + // if the previous range exists and contains index, that's the one we want + return ((after != ranges.begin()) && (index < std::prev(after)->end_index())) ? std::prev(after) : + after; +} // lar::sparse_vector::find_range_iter_at_or_after() const + +template +typename lar::sparse_vector::range_iterator lar::sparse_vector::find_range_iter_at_or_after( + size_type index) +{ + // this range has the offset (first index) above the index argument: + auto after = find_next_range_iter(index); + // if the previous range exists and contains index, that's the one we want + return ((after != ranges.begin()) && (index < std::prev(after)->end_index())) ? std::prev(after) : + after; +} // lar::sparse_vector::find_range_iter_at_or_after() + +template +typename lar::sparse_vector::range_iterator lar::sparse_vector::find_extending_range_iter( + size_type index, + range_iterator rbegin) +{ + // this range has the offset (first index) above the index argument: + auto it = find_next_range_iter(index, rbegin); + // if index were not void, would it belong to the previous range? + // if so, the previous range is the one we want + return ((it != rbegin) && std::prev(it)->borders(index)) ? std::prev(it) : it; +} // lar::sparse_vector::find_extending_range_iter() + +template +typename lar::sparse_vector::range_const_iterator +lar::sparse_vector::find_extending_range_iter(size_type index, range_const_iterator rbegin) const +{ + // this range has the offset (first index) above the index argument: + auto it = find_next_range_iter(index, rbegin); + // if index were not void, would it belong to the previous range? + // if so, the previus range is the one we want + return ((it != rbegin) && std::prev(it)->borders(index)) ? std::prev(it) : it; +} // lar::sparse_vector::find_extending_range_iter() const + +template +const typename lar::sparse_vector::datarange_t& lar::sparse_vector::add_range_before( + size_type offset, + vector_t&& new_data, + range_iterator nextRange) +{ + // insert the new range before the existing range which starts after offset + range_iterator iInsert = nextRange; + + // is there a range before this, which includes the offset? + if ((iInsert != ranges.begin()) && (iInsert - 1)->borders(offset)) { + // then we should extend it + (--iInsert)->extend(offset, new_data.begin(), new_data.end()); + } + else { + // no range before the insertion one includes the offset of the new range; + // ... we need to add it as a new range; + // it is not really clear to me [GP] why I need a std::move here, since + // new_data is a rvalue already; in doubt, I have painted all this kind + // of constructs with move()s, just in case + iInsert = insert_range(iInsert, {offset, std::move(new_data)}); + } + return merge_ranges(iInsert); +} // lar::sparse_vector::add_range_before(vector, iterator) + +template +typename lar::sparse_vector::datarange_t& lar::sparse_vector::merge_ranges( + range_iterator iRange) +{ + range_iterator iNext = iRange + 1; + while (iNext != ranges.end()) { + if (!iRange->borders(iNext->begin_index())) break; + // iRange content dominates, but if iNext has data beyond it, + // then we copy it + if (iNext->end_index() > iRange->end_index()) { + iRange->extend(iRange->end_index(), iNext->get_iterator(iRange->end_index()), iNext->end()); + } + iNext = ranges.erase(iNext); + } // while + fix_size(); + return *iRange; +} // lar::sparse_vector::merge_ranges() + +template +typename lar::sparse_vector::range_iterator lar::sparse_vector::eat_range_head( + range_iterator iRange, + size_t index) +{ + if (index <= iRange->begin_index()) return iRange; + if (index >= iRange->end_index()) return ranges.erase(iRange); + iRange->move_head(index); + return iRange; +} // lar::sparse_vector::eat_range_head() + +template +typename lar::sparse_vector::size_type lar::sparse_vector::fix_size() +{ + if (!ranges.empty()) nominal_size = std::max(nominal_size, ranges.back().end_index()); + return nominal_size; +} // lar::sparse_vector::fix_size() + +// --- static methods + +template +inline size_t lar::sparse_vector::expected_vector_size(size_t size) +{ + // apparently, a chunk of heap memory takes at least 32 bytes; + // that means that a vector of 1 or 5 32-bit integers takes the same + // space; the overhead appears to be 8 bytes, which can be allocated + return sizeof(vector_t) + std::max(size_t(32), (alignof(datarange_t) * size + 8)); +} // lar::sparse_vector::expected_vector_size() + +template +inline size_t lar::sparse_vector::min_gap() +{ + // we assume here that there is no additional overhead by alignment; + // the gap adds the space of another datarange_t, including the vector, + // its data and overhead from heap (apparently, 8 bytes); + // + return (sizeof(datarange_t) + 8) / sizeof(value_type) + 1; // round up +} // lar::sparse_vector::min_gap() + +template +inline bool lar::sparse_vector::should_merge(const typename datarange_t::base_t& a, + const typename datarange_t::base_t& b) +{ + size_type gap_size = (a < b) ? b.begin_index() - a.begin_index() - a.size() : + a.begin_index() - b.begin_index() - b.size(); + return expected_vector_size(a.size() + b.size() + gap_size) <= + expected_vector_size(a.size()) + expected_vector_size(b.size()); +} // lar::sparse_vector::should_merge() + +// --- non-member functions +template +std::ostream& operator<<(std::ostream& out, const lar::sparse_vector& v) +{ + + out << "Sparse vector of size " << v.size() << " with " << v.get_ranges().size() << " ranges:"; + typename lar::sparse_vector::range_const_iterator iRange = v.begin_range(), + rend = v.end_range(); + while (iRange != rend) { + out << "\n "; + iRange->dump(out); + /* + out << "\n [" << iRange->begin_index() << " - " << iRange->end_index() + << "] (" << iRange->size() << "):"; + typename lar::sparse_vector::datarange_t::const_iterator + iValue = iRange->begin(), vend = iRange->end(); + while (iValue != vend) out << " " << (*(iValue++)); + */ + ++iRange; + } // for + return out << std::endl; +} // operator<< (ostream, sparse_vector) + +// ----------------------------------------------------------------------------- +// --- lar::sparse_vector::datarange_t implementation +// --- +template +template +typename lar::sparse_vector::datarange_t& +lar::sparse_vector::datarange_t::extend(size_type index, ITER first, ITER last) +{ + size_type new_size = + std::max(base_t::relative_index(index) + std::distance(first, last), base_t::size()); + base_t::resize(new_size); + values.resize(new_size); + std::copy(first, last, get_iterator(index)); + return *this; +} // lar::sparse_vector::datarange_t::extend() + +template +void lar::sparse_vector::datarange_t::move_head(size_type to_index, + value_type def_value /* = value_zero */) +{ + difference_type delta = to_index - base_t::begin_index(); + if (delta == 0) return; + base_t::move_head(delta); + if (delta > 0) // shrink + values.erase(values.begin(), values.begin() + delta); + else { // extend + values.insert(values.begin(), + value_const_iterator(def_value), + value_const_iterator(def_value) - delta); + } +} // lar::sparse_vector::datarange_t::move_head() + +template +template +void lar::sparse_vector::datarange_t::dump(Stream&& out) const +{ + out << "[" << this->begin_index() << " - " << this->end_index() << "] (" << this->size() + << "): {"; + for (auto const& v : this->values) + out << " " << v; + out << " }"; +} // lar::sparse_vector::datarange_t::dump() + +// ----------------------------------------------------------------------------- +// --- lar::sparse_vector::const_iterator implementation +// --- +template +typename lar::sparse_vector::const_iterator& lar::sparse_vector::const_iterator::operator++() +{ + // no container, not doing anything; + // index beyond the end: stays there + if (!cont || (index >= cont->size())) return *this; + + // increment the internal index + ++index; + + // were we in a valid range? + if (currentRange != cont->ranges.end()) { + // if the new index is beyond the current range, pick the next + if (currentRange->end_index() <= index) ++currentRange; + } + // if we have no valid range, we are forever in the void + return *this; +} // lar::sparse_vector::iterator::operator++() + +template +typename lar::sparse_vector::const_iterator::const_reference +lar::sparse_vector::const_iterator::operator*() const +{ + // no container, no idea what to do + if (!cont) throw std::out_of_range("iterator to no sparse vector"); + + // index beyond the end: can't return any reference + if (index >= cont->size()) return value_zero; + + // are we in a valid range? if not, we are past the last range + if (currentRange == cont->ranges.end()) return value_zero; + + // if the index is beyond the current range, we are in the void + if (index < currentRange->begin_index()) return value_zero; + + return (*currentRange)[index]; +} // lar::sparse_vector::const_iterator::operator*() + +template +typename lar::sparse_vector::const_iterator& lar::sparse_vector::const_iterator::operator+=( + difference_type delta) +{ + if (delta == 1) return this->operator++(); + index += delta; + if ((currentRange == cont->ranges.end()) || !currentRange->includes(index)) refresh_state(); + return *this; +} // lar::sparse_vector::const_iterator::operator+=() + +template +inline typename lar::sparse_vector::const_iterator& +lar::sparse_vector::const_iterator::operator-=(difference_type delta) +{ + return this->operator+=(-delta); +} + +template +typename lar::sparse_vector::const_iterator lar::sparse_vector::const_iterator::operator+( + difference_type delta) const +{ + if ((currentRange == cont->ranges.end()) || !currentRange->includes(index + delta)) + return const_iterator(*cont, index + delta); + const_iterator iter(*this); + iter.index += delta; + return iter; +} // lar::sparse_vector::const_iterator::operator+() + +template +inline typename lar::sparse_vector::const_iterator +lar::sparse_vector::const_iterator::operator-(difference_type delta) const +{ + return this->operator+(-delta); +} + +/// distance operator +template +inline typename lar::sparse_vector::const_iterator::difference_type +lar::sparse_vector::const_iterator::operator-(const const_iterator& iter) const +{ + if (cont != iter.cont) { + throw std::runtime_error("lar::sparse_vector::const_iterator:" + " difference with alien iterator"); + } + return index - iter.index; +} // lar::sparse_vector::const_iterator::operator-(const_iterator) + +template +void lar::sparse_vector::const_iterator::refresh_state() +{ + // update the currentRange + // currentRange is the range including the current item, or next to it + if (cont) { + // the following is actually the range after the index: + currentRange = cont->find_next_range_iter(index); + // if the index is inside the previous index, then we want to move there: + if (currentRange != cont->ranges.begin()) { + if ((currentRange - 1)->end_index() > index) --currentRange; + } + } + else { + currentRange = {}; + } +} // lar::sparse_vector::const_iterator::refresh_state() + +// ----------------------------------------------------------------------------- +// --- lar::sparse_vector::iterator implementation +// --- +// +// nothing new so far +// + +#endif // LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H diff --git a/gaus_hit_finder/find_hits_with_gaussians.cpp b/gaus_hit_finder/find_hits_with_gaussians.cpp new file mode 100644 index 0000000..4c2df9d --- /dev/null +++ b/gaus_hit_finder/find_hits_with_gaussians.cpp @@ -0,0 +1,478 @@ +// See REAMDME.md for some general comments about this example. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "tbb/concurrent_vector.h" +#include "tbb/parallel_for.h" + +#include "find_hits_with_gaussians.hpp" +#include "copied_from_larsoft_minor_edits/ICandidateHitFinder.h" +#include "print_hits.hpp" + +namespace { + + // This is an edited copy of the TMath::Gaus function from ROOT, since we + // don't want to depend on ROOT in this example. + double Gaus(double x, double mean, double sigma, bool norm) { + if (sigma == 0) return 1.e30; + double arg = (x-mean)/sigma; + // for |arg| > 39 result is zero in double precision + if (arg < -39.0 || arg > 39.0) return 0.0; + double res = std::exp(-0.5*arg*arg); + if (!norm) return res; + return res/(2.50662827463100024*sigma); //sqrt(2*Pi)=2.50662827463100024 + } +} + +namespace examples { + std::vector find_hits_with_gaussians(find_hits_with_gaussians_cfg const& cfg, + std::vector const& wires, + std::vector> const& cand_hit_standard, + PeakFitterMrqdt const& peak_fitter_mrqdt, + HitFilterAlg const& hit_filter_alg) { + + std::cout << "Finding hits with Gaussians." << std::endl; + + // auto const& wireReadoutGeom = art::ServiceHandle()->Get(); + + //store in a thread safe way + struct hitstruct { + recob::Hit hit_tbb; + }; + + tbb::concurrent_vector hitstruct_vec; + tbb::concurrent_vector filthitstruct_vec; + + //################################################# + //### Set the charge determination method ### + //### Default is to compute the normalized area ### + //################################################# + std::function chargeFunc = + [](double /* peakMean */, + double peakAmp, + double peakWidth, + double areaNorm, + int /* low */, + int /* hi */) { return std::sqrt(2 * std::numbers::pi) * peakAmp * peakWidth / areaNorm; }; + + //############################################## + //### Alternative is to integrate over pulse ### + //############################################## + if (cfg.area_method == 0) + chargeFunc = [](double peakMean, + double peakAmp, + double peakWidth, + double /* areaNorm */, + int low, + int hi) { + double charge(0); + for (int sigPos = low; sigPos < hi; sigPos++) + charge += peakAmp * Gaus(sigPos, peakMean, peakWidth, false); + return charge; + }; + + //############################## + //### Looping over the wires ### + //############################## + tbb::parallel_for( + static_cast(0), + wires.size(), + [&](size_t& wireIter) { + // #################################### + // ### Getting this particular wire ### + // #################################### + recob::Wire const* wire = &wires[wireIter]; + + // --- Setting Channel Number and Signal type --- + + raw::ChannelID_t channel = wire->Channel(); + + // TODO. The plane number should come from the Geometry. + // That is not implemented yet... Below is the commented out + // version of the code that worked with the art framework. + geo::PlaneID::PlaneID_t plane = 0; + + // get the WireID for this hit + // std::vector wids = wireReadoutGeom.ChannelToWire(channel); + // for now, just take the first option returned from ChannelToWire + // geo::WireID wid = wids[0]; + // We need to know the plane to look up parameters + // geo::PlaneID::PlaneID_t plane = wid.Plane; + + // ---------------------------------------------------------- + // -- Setting the appropriate signal widths and thresholds -- + // -- for the right plane. -- + // ---------------------------------------------------------- + + // ################################################# + // ### Set up to loop over ROI's for this wire ### + // ################################################# + const recob::Wire::RegionsOfInterest_t& signalROI = wire->SignalROI(); + + tbb::parallel_for( + static_cast(0), + signalROI.n_ranges(), + [&](size_t& rangeIter) { + const auto& range = signalROI.range(rangeIter); + // ROI start time + raw::TDCtick_t roiFirstBinTick = range.begin_index(); + + // ########################################################### + // ### Scan the waveform and find candidate peaks + merge ### + // ########################################################### + + examples::ICandidateHitFinder::HitCandidateVec hitCandidateVec; + examples::ICandidateHitFinder::MergeHitCandidateVec mergedCandidateHitVec; + + cand_hit_standard.at(plane)->findHitCandidates( + range, 0, channel, hitCandidateVec); + cand_hit_standard.at(plane)->MergeHitCandidates( + range, hitCandidateVec, mergedCandidateHitVec); + + // ####################################################### + // ### Lets loop over the pulses we found on this wire ### + // ####################################################### + + for (auto& mergedCands : mergedCandidateHitVec) { + int startT = mergedCands.front().startTick; + int endT = mergedCands.back().stopTick; + + // ### Putting in a protection in case things went wrong ### + // ### In the end, this primarily catches the case where ### + // ### a fake pulse is at the start of the ROI ### + if (endT - startT < 5) continue; + + // ####################################################### + // ### Clearing the parameter vector for the new Pulse ### + // ####################################################### + + // === Setting The Number Of Gaussians to try === + int nGausForFit = mergedCands.size(); + + // ################################################## + // ### Calling the function for fitting Gaussians ### + // ################################################## + double chi2PerNDF(0.); + int NDF(1); + /*stand alone + reco_tool::IPeakFitter::PeakParamsVec peakParamsVec(nGausForFit); + */ + examples::IPeakFitter::PeakParamsVec peakParamsVec; + + // ####################################################### + // ### If # requested Gaussians is too large then punt ### + // ####################################################### + if (mergedCands.size() <= cfg.max_multi_hit) { + peak_fitter_mrqdt.findPeakParameters( + range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF); + + // If the chi2 is infinite then there is a real problem so we bail + if (!(chi2PerNDF < std::numeric_limits::infinity())) { + chi2PerNDF = 2. * cfg.chi2_ndf; + NDF = 2; + } + } + + // ####################################################### + // ### If too large then force alternate solution ### + // ### - Make n hits from pulse train where n will ### + // ### depend on the fhicl parameter fLongPulseWidth ### + // ### Also do this if chi^2 is too large ### + // ####################################################### + if (mergedCands.size() > cfg.max_multi_hit || nGausForFit * chi2PerNDF > cfg.chi2_ndf) { + int longPulseWidth = cfg.long_pulse_width_vec.at(plane); + int nHitsThisPulse = (endT - startT) / longPulseWidth; + + if (nHitsThisPulse > cfg.long_max_hits_vec.at(plane)) { + nHitsThisPulse = cfg.long_max_hits_vec.at(plane); + longPulseWidth = (endT - startT) / nHitsThisPulse; + } + + if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++; + + int firstTick = startT; + int lastTick = std::min(firstTick + longPulseWidth, endT); + + peakParamsVec.clear(); + nGausForFit = nHitsThisPulse; + NDF = 1.; + chi2PerNDF = chi2PerNDF > cfg.chi2_ndf ? chi2PerNDF : -1.; + + for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) { + // This hit parameters + double ROIsumADC = + std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.); + double peakSigma = (lastTick - firstTick) / 3.; // Set the width... + double peakAmp = 0.3989 * ROIsumADC / peakSigma; // Use gaussian formulation + double peakMean = (firstTick + lastTick) / 2.; + + // Store hit params + examples::IPeakFitter::PeakFitParams_t peakParams; + + peakParams.peakCenter = peakMean; + peakParams.peakCenterError = 0.1 * peakMean; + peakParams.peakSigma = peakSigma; + peakParams.peakSigmaError = 0.1 * peakSigma; + peakParams.peakAmplitude = peakAmp; + peakParams.peakAmplitudeError = 0.1 * peakAmp; + + peakParamsVec.push_back(peakParams); + + // set for next loop + firstTick = lastTick; + lastTick = std::min(lastTick + longPulseWidth, endT); + } + } + + // ####################################################### + // ### Loop through returned peaks and make recob hits ### + // ####################################################### + + int numHits(0); + + // Make a container for what will be the filtered collection + std::vector filteredHitVec; + + float nextpeak(0); + float prevpeak(0); + float nextpeakSig(0); + float prevpeakSig(0); + float nsigmaADC(2.0); + float newright(0); + float newleft(0); + for (const auto& peakParams : peakParamsVec) { + // Extract values for this hit + float peakAmp = peakParams.peakAmplitude; + float peakMean = peakParams.peakCenter; + float peakWidth = peakParams.peakSigma; + + //std::cout<<" ans hits "< 0) { + prevpeak = (peakParamsVec.at(numHits - 1)).peakCenter; + prevpeakSig = (peakParamsVec.at(numHits - 1)).peakSigma; + //std::cout<<" ans size "<::const_iterator sumStartItr = range.begin() + startT; + std::vector::const_iterator sumEndItr = range.begin() + endT; + + //### limits for the sum of the Hit based on the gaussian peak and sigma + std::vector::const_iterator HitsumStartItr = + range.begin() + peakMean - nsigmaADC * peakWidth; + std::vector::const_iterator HitsumEndItr = + range.begin() + peakMean + nsigmaADC * peakWidth; + + if (nGausForFit > 1) { + if (numHits > 0) { + if ((peakMean - nsigmaADC * peakWidth) < (prevpeak + nsigmaADC * prevpeakSig)) { + float difPeak = peakMean - prevpeak; + float weightpeak = prevpeakSig / (prevpeakSig + peakWidth); + HitsumStartItr = range.begin() + prevpeak + difPeak * weightpeak; + newleft = prevpeak + difPeak * weightpeak; + } + } + + if (numHits < nGausForFit - 1) { + if ((peakMean + nsigmaADC * peakWidth) > (nextpeak - nsigmaADC * nextpeakSig)) { + float difPeak = nextpeak - peakMean; + float weightpeak = peakWidth / (nextpeakSig + peakWidth); + HitsumEndItr = range.begin() + peakMean + difPeak * weightpeak; + newright = peakMean + difPeak * weightpeak; + } + } + } + + //protection to avoid negative ranges + if (newright - newleft < 0) continue; + + //avoid ranges out of ROI if it happens + if (HitsumStartItr < sumStartItr) HitsumStartItr = sumStartItr; + + if (HitsumEndItr > sumEndItr) HitsumEndItr = sumEndItr; + + if (HitsumStartItr > HitsumEndItr) continue; + + // ### Sum of ADC counts + double ROIsumADC = std::accumulate(sumStartItr, sumEndItr, 0.); + double HitsumADC = std::accumulate(HitsumStartItr, HitsumEndItr, 0.); + + recob::Hit hit(wire->Channel(), + startT + roiFirstBinTick, + endT + roiFirstBinTick, + peakMean + roiFirstBinTick, + peakMeanErr, + peakWidth, + peakAmp, + peakAmpErr, + ROIsumADC, + HitsumADC, + charge, + chargeErr, + nGausForFit, + numHits, + chi2PerNDF, + NDF, + wire->View(), + // Geometry system for Phlex is not yet implemented, + // so we just set signal type to 0 (kInduction) for now. + geo::kInduction, + // art::ServiceHandle()->Get().SignalType(wire.Channel()), + // wid also comes from the Geometry, which is not yet implemented, so we just set + // it to a default value for now. + geo::WireID()); + // wid); + + if (cfg.filter_hits) filteredHitVec.push_back(hit); + + // This loop will store ALL hits + hitstruct tmp{std::move(hit)}; + hitstruct_vec.push_back(std::move(tmp)); + + numHits++; + } // <---End loop over gaussians + + // Should we filter hits? + if (cfg.filter_hits && !filteredHitVec.empty()) { + // ####################################################################### + // Is all this sorting really necessary? Would it be faster to just loop + // through the hits and perform simple cuts on amplitude and width on a + // hit-by-hit basis, either here in the module (using fPulseHeightCuts and + // fPulseWidthCuts) or in HitFilterAlg? + // ####################################################################### + + // Sort in ascending peak height + // (I believe the preceding comment is incorrect. The sort below + // is in descending order of peak height, not ascending.) + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](const auto& left, const auto& right) { + return left.PeakAmplitude() > right.PeakAmplitude(); + }); + + // Reject if the first hit fails the PH/wid cuts + if (filteredHitVec.front().PeakAmplitude() < cfg.pulse_height_cuts.at(plane) || + filteredHitVec.front().RMS() < cfg.pulse_width_cuts.at(plane)) + filteredHitVec.clear(); + + // Now check other hits in the snippet + if (filteredHitVec.size() > 1) { + // The largest pulse height will now be at the front... + float largestPH = filteredHitVec.front().PeakAmplitude(); + + // Find where the pulse heights drop below threshold + float threshold(cfg.pulse_ratio_cuts.at(plane)); + + std::vector::iterator smallHitItr = + std::find_if(filteredHitVec.begin(), + filteredHitVec.end(), + [largestPH, threshold](const auto& hit) { + return hit.PeakAmplitude() < 8. && + hit.PeakAmplitude() / largestPH < threshold; + }); + + // Shrink to fit + if (smallHitItr != filteredHitVec.end()) + filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr)); + + // Resort in time order + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](const auto& left, const auto& right) { + return left.PeakTime() < right.PeakTime(); + }); + } + + // Copy the hits we want to keep to the filtered hit collection + for (const auto& filteredHit : filteredHitVec) { + if (!cfg.filter_hits || hit_filter_alg.IsGoodHit(filteredHit)) { + hitstruct tmp{std::move(filteredHit)}; + filthitstruct_vec.push_back(std::move(tmp)); + } + } + } + } //<---End loop over merged candidate hits + } //<---End looping over ROI's + ); //end tbb parallel for + } //<---End looping over all the wires + ); //end tbb parallel for + + std::vector hits; + + if (cfg.filter_hits) { + for (size_t i = 0; i < filthitstruct_vec.size(); i++) { + hits.emplace_back(filthitstruct_vec[i].hit_tbb); + } + } else { + for (size_t i = 0; i < hitstruct_vec.size(); i++) { + hits.emplace_back(hitstruct_vec[i].hit_tbb); + } + } + + // Sort hits by PeakTime ascending, then by PeakAmplitude ascending + // I added the sorting when trying to compare the output of this + // code with output from the art version of GausHitFinder. The parallel_for + // loops above scramble the order and make comparisons difficult. + // This should be temporary and removed. + std::sort(hits.begin(), hits.end(), [](const recob::Hit& a, const recob::Hit& b) { + if (a.PeakTime() != b.PeakTime()) { + return a.PeakTime() < b.PeakTime(); + } + return a.PeakAmplitude() < b.PeakAmplitude(); + }); + + // Temporary, only for testing purposes, print + // the hits to a file here. This should also be + // removed eventually. + static std::atomic i = 0; + int current_value = i.fetch_add(1); + std::string filename = std::string("hits_") + + std::to_string(current_value) + ".txt"; + print_hits_to_file(hits, filename); + + return hits; + } // End of find_hits_with_gaussians +} // end of examples namespace diff --git a/gaus_hit_finder/find_hits_with_gaussians.hpp b/gaus_hit_finder/find_hits_with_gaussians.hpp new file mode 100644 index 0000000..f1eefd5 --- /dev/null +++ b/gaus_hit_finder/find_hits_with_gaussians.hpp @@ -0,0 +1,65 @@ +#ifndef PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_HPP +#define PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_HPP + +// See REAMDME.md for some general comments about this example. + +//////////////////////////////////////////////////////////////////////// +// +// GaussHitFinder class +// +// jaasaadi@syr.edu +// +// This algorithm is designed to find hits on wires after deconvolution. +// ----------------------------------- +// This algorithm is based on the FFTHitFinder written by Brian Page, +// Michigan State University, for the ArgoNeuT experiment. +// +// +// The algorithm walks along the wire and looks for pulses above threshold +// The algorithm then attempts to fit n-gaussians to these pulses where n +// is set by the number of peaks found in the pulse +// If the Chi2/NDF returned is "bad" it attempts to fit n+1 gaussians to +// the pulse. If this is a better fit it then uses the parameters of the +// Gaussian fit to characterize the "hit" object +// +// To use this simply include the following in your producers: +// gaushit: @local::microboone_gaushitfinder +// gaushit: @local::argoneut_gaushitfinder +//////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "copied_from_larsoft_minor_edits/CandHitStandard.h" +#include "copied_from_larsoft_minor_edits/Hit.h" +#include "copied_from_larsoft_minor_edits/HitFilterAlg.h" +#include "copied_from_larsoft_minor_edits/PeakFitterMrqdt.h" +#include "copied_from_larsoft_minor_edits/Wire.h" + + +namespace examples { + + struct find_hits_with_gaussians_cfg { + bool filter_hits; + + std::vector long_max_hits_vec; /// long_pulse_width_vec; /// + area_norms_vec; /// pulse_height_cuts; + std::vector pulse_width_cuts; + std::vector pulse_ratio_cuts; + }; + + std::vector find_hits_with_gaussians(find_hits_with_gaussians_cfg const& cfg, + std::vector const& wires, + std::vector> const& cand_hit_standard, + PeakFitterMrqdt const& peak_fitter_mrqdt, + HitFilterAlg const& hit_filter_alg); + +} +#endif // PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_HPP diff --git a/gaus_hit_finder/print_hits.cpp b/gaus_hit_finder/print_hits.cpp new file mode 100644 index 0000000..ec4849a --- /dev/null +++ b/gaus_hit_finder/print_hits.cpp @@ -0,0 +1,48 @@ +#include +#include +#include +#include "print_hits.hpp" + +void print_hits_to_file(std::vector const& hits, + std::string const& filename) { + std::FILE* file = std::fopen(filename.c_str(), "w"); + if (!file) { + throw std::runtime_error("Failed to open file for writing: " + filename); + return; + } + + fmt::print(file, "Contents of a std::vector\n\n"); + + int default_precision = 6; + + for (int i = 0; auto const& hit : hits) { + fmt::print(file, "i = {}\n", i++); + fmt::print(file, "{}\n", hit.Channel()); + fmt::print(file, "{}\n", hit.StartTick()); + fmt::print(file, "{}\n", hit.EndTick()); + fmt::print(file, "{:.{}g}\n", hit.PeakTime(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.SigmaPeakTime(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.RMS(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.PeakAmplitude(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.SigmaPeakAmplitude(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.ROISummedADC(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.HitSummedADC(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.Integral(), default_precision); + fmt::print(file, "{:.{}g}\n", hit.SigmaIntegral(), default_precision); + fmt::print(file, "{}\n", hit.Multiplicity()); + fmt::print(file, "{}\n", hit.LocalIndex()); + fmt::print(file, "{:.{}g}\n", hit.GoodnessOfFit(), default_precision); + fmt::print(file, "{}\n", hit.DegreesOfFreedom()); + fmt::print(file, "{}\n", static_cast(hit.View())); + + // These are commented out for now because I am + // using this print out to make a comparison to + // the output of the art version of GausHitFinder. + // These fields depend on Geometry and that is + // not implemented yet in phlex. + // + // hit.SignalType() + // hit.WireID() + } + std::fclose(file); +} diff --git a/gaus_hit_finder/print_hits.hpp b/gaus_hit_finder/print_hits.hpp new file mode 100644 index 0000000..bbe52c1 --- /dev/null +++ b/gaus_hit_finder/print_hits.hpp @@ -0,0 +1,11 @@ +#ifndef PHLEX_EXAMPLES_PRINT_HITS_HPP +#define PHLEX_EXAMPLES_PRINT_HITS_HPP + +#include +#include +#include "copied_from_larsoft_minor_edits/Hit.h" + +void print_hits_to_file(std::vector const& hits, + std::string const& filename); + +#endif // PHLEX_EXAMPLES_PRINT_HITS_HPP diff --git a/gaus_hit_finder/register_find_hits_with_gaussians.cpp b/gaus_hit_finder/register_find_hits_with_gaussians.cpp new file mode 100644 index 0000000..3c593e0 --- /dev/null +++ b/gaus_hit_finder/register_find_hits_with_gaussians.cpp @@ -0,0 +1,89 @@ +// See REAMDME.md for some general comments about this example. + +#include +#include +#include +#include +#include +#include + +#include "copied_from_larsoft_minor_edits/CandHitStandard.h" +#include "copied_from_larsoft_minor_edits/HitFilterAlg.h" +#include "copied_from_larsoft_minor_edits/PeakFitterMrqdt.h" +#include "copied_from_larsoft_minor_edits/Wire.h" + +#include "phlex/module.hpp" +#include "find_hits_with_gaussians.hpp" + +using namespace phlex; + +PHLEX_REGISTER_ALGORITHMS(m, config) +{ + auto const layer = config.get("layer"); + + examples::find_hits_with_gaussians_cfg cfg = { + .filter_hits = config.get("filter_hits"), + .long_max_hits_vec = config.get>("long_max_hits_vec"), + .long_pulse_width_vec = config.get>("long_pulse_width_vec"), + .max_multi_hit = config.get("max_multi_hit"), + .area_method = config.get("area_method"), + .area_norms_vec = config.get>("area_norms_vec"), + .chi2_ndf = config.get("chi2_ndf"), + .pulse_height_cuts = config.get>("pulse_height_cuts"), + .pulse_width_cuts = config.get>("pulse_width_cuts"), + .pulse_ratio_cuts = config.get>("pulse_ratio_cuts") + }; + + std::vector> cand_hit_standard_vec; + + auto finder_configs = config.get("cand_hit_standard_configs"); + cand_hit_standard_vec.resize(finder_configs.keys().size()); + for (auto const& key : finder_configs.keys()) { + auto finder_config = finder_configs.get(key); + auto hit_finder = std::make_shared( + examples::CandHitStandardCfg{ + .fRoiThreshold = finder_config.get("roiThreshold") + }); + unsigned int plane = finder_config.get("Plane"); + if (plane >= cand_hit_standard_vec.size()) { + std::cerr << "Error: plane number " << plane << " is out of range for cand_hit_standard_vec of size " << cand_hit_standard_vec.size() << std::endl; + throw std::runtime_error("Invalid plane number"); + } + cand_hit_standard_vec[plane] = std::move(hit_finder); + } + + std::shared_ptr peak_fitter_mrqdt; + auto fitter_config = config.get("peak_fitter_mrqdt_config"); + peak_fitter_mrqdt = std::make_shared( + examples::PeakFitterMrqdtCfg{ + .fMinWidth = fitter_config.get("min_width"), + .fMaxWidthMult = fitter_config.get("max_width_mult"), + .fPeakRange = fitter_config.get("peak_range_fact"), + .fAmpRange = fitter_config.get("peak_amp_range") + }); + + + std::shared_ptr hit_filter_alg; + auto filter_config = config.get("hit_filter_alg_config"); + hit_filter_alg = std::make_shared( + examples::HitFilterAlgCfg{ + .fMinPulseHeight = filter_config.get>("min_pulse_height"), + .fMinPulseSigma = filter_config.get>("min_pulse_sigma") + }); + + m.transform("find_hits_with_gaussians", + [cfg = std::move(cfg), + cand_hit_standard_vec = std::move(cand_hit_standard_vec), + peak_fitter_mrqdt = std::move(peak_fitter_mrqdt), + hit_filter_alg = std::move(hit_filter_alg)] + (std::vector const& wires) { + return examples::find_hits_with_gaussians(cfg, + wires, + cand_hit_standard_vec, + *peak_fitter_mrqdt, + *hit_filter_alg); + }, + concurrency::unlimited) + .input_family("wires"_in(layer)) + .output_products("hits"); +} diff --git a/gaus_hit_finder/test_find_hits_with_gaussians.jsonnet b/gaus_hit_finder/test_find_hits_with_gaussians.jsonnet new file mode 100644 index 0000000..114e8ef --- /dev/null +++ b/gaus_hit_finder/test_find_hits_with_gaussians.jsonnet @@ -0,0 +1,61 @@ +{ + driver: { + cpp: 'generate_layers', + layers: { + spill: { parent: 'job', total: 5, starting_number: 1 }, + }, + }, + sources: { + wires_source: { + cpp: 'wires_source', + layer: 'spill', + }, + }, + modules: { + find_hits_with_gaussians_cpp: { + cpp: 'find_hits_with_gaussians_hof', + layer: 'spill', + // Copied in the configuration parameters that + // I obtained from "Far Detector Simulation/Reconstruction conveners": + // Dominic Brailsford and Laura Paulucci (in cc). See Dom and Kyle + // 2/17/2026 on SLACK and also email around the same time. These were + // copied out of the expanded fcl configuration. Note that the ones + // that have different values for different planes are not used when + // filter_hits is false. + filter_hits: false, + long_max_hits_vec: [1, 1, 1], + long_pulse_width_vec: [16, 16, 16], + max_multi_hit: 4, + area_method: 0, + area_norms_vec: [13.25, 13.25, 13.25], + chi2_ndf: 50.0, + pulse_height_cuts: [3.0, 3.0, 3.0], + pulse_width_cuts: [2.0, 1.5, 1.0], + pulse_ratio_cuts: [0.35, 0.40, 0.20], + cand_hit_standard_configs: { + CandidateHitsPlane0: { + Plane: 0, + roiThreshold: 6.0, + }, + CandidateHitsPlane1: { + Plane: 1, + roiThreshold: 6.0, + }, + CandidateHitsPlane2: { + Plane: 2, + roiThreshold: 6.0, + }, + }, + peak_fitter_mrqdt_config: { + min_width: 0.5, + max_width_mult: 3.0, + peak_range_fact: 2.0, + peak_amp_range: 2.0, + }, + hit_filter_alg_config: { + min_pulse_height: [5.0, 5.0, 5.0], + min_pulse_sigma: [1.0, 1.0, 1.0], + }, + }, + }, +} diff --git a/gaus_hit_finder/wire_serialization.cpp b/gaus_hit_finder/wire_serialization.cpp new file mode 100644 index 0000000..ff67128 --- /dev/null +++ b/gaus_hit_finder/wire_serialization.cpp @@ -0,0 +1,134 @@ +#include +#include +#include "copied_from_larsoft_minor_edits/sparse_vector.h" +#include "wire_serialization.hpp" + +bool write_wires_to_file(std::vector const& wires, + std::string const& filename) { + std::ofstream out(filename, std::ios::binary); + if (!out) { + throw std::runtime_error("Failed to open file for writing: " + filename); + return false; + } + + // Write number of wires + size_t num_wires = wires.size(); + out.write(reinterpret_cast(&num_wires), sizeof(num_wires)); + + for (const auto& wire : wires) { + // Write channel ID + raw::ChannelID_t channel = wire.Channel(); + out.write(reinterpret_cast(&channel), sizeof(channel)); + + // Write view + geo::View_t view = wire.View(); + out.write(reinterpret_cast(&view), sizeof(view)); + + // Write sparse vector data + const auto& signal_roi = wire.SignalROI(); + size_t nominal_size = signal_roi.size(); + out.write(reinterpret_cast(&nominal_size), sizeof(nominal_size)); + + // Write number of ranges + size_t num_ranges = std::distance(signal_roi.begin_range(), signal_roi.end_range()); + out.write(reinterpret_cast(&num_ranges), sizeof(num_ranges)); + + // Write each range + for (auto range_it = signal_roi.begin_range(); range_it != signal_roi.end_range(); ++range_it) { + // Write range offset and last + size_t offset = range_it->begin_index(); + size_t last = range_it->end_index(); + out.write(reinterpret_cast(&offset), sizeof(offset)); + out.write(reinterpret_cast(&last), sizeof(last)); + + // Write number of values + size_t num_values = range_it->size(); + out.write(reinterpret_cast(&num_values), sizeof(num_values)); + + // Write values + if (num_values > 0) { + // Create a temporary vector to hold the values + std::vector values(range_it->begin(), range_it->end()); + out.write(reinterpret_cast(values.data()), + num_values * sizeof(float)); + } + } + } + out.close(); + return out.good(); +} + +bool read_wires_from_file(std::string const& filename, + std::vector& wires) { + std::ifstream in(filename, std::ios::binary); + if (!in) { + throw std::runtime_error("Failed to open file for reading: " + filename); + return false; + } + + wires.clear(); + + // Read number of wires + size_t num_wires; + in.read(reinterpret_cast(&num_wires), sizeof(num_wires)); + if (!in) return false; + + wires.reserve(num_wires); + + for (size_t i = 0; i < num_wires; ++i) { + // Read channel ID + raw::ChannelID_t channel; + in.read(reinterpret_cast(&channel), sizeof(channel)); + if (!in) return false; + + // Read view + geo::View_t view; + in.read(reinterpret_cast(&view), sizeof(view)); + if (!in) return false; + + // Read sparse vector data + size_t nominal_size; + in.read(reinterpret_cast(&nominal_size), sizeof(nominal_size)); + if (!in) return false; + + // Read number of ranges + size_t num_ranges; + in.read(reinterpret_cast(&num_ranges), sizeof(num_ranges)); + if (!in) return false; + + // Create sparse vector + lar::sparse_vector signal_roi; + signal_roi.resize(nominal_size); + + // Read each range + for (size_t j = 0; j < num_ranges; ++j) { + // Read range offset and last + size_t offset, last; + in.read(reinterpret_cast(&offset), sizeof(offset)); + in.read(reinterpret_cast(&last), sizeof(last)); + if (!in) return false; + + // Read number of values + size_t num_values; + in.read(reinterpret_cast(&num_values), sizeof(num_values)); + if (!in) return false; + + // Read values + std::vector values(num_values); + if (num_values > 0) { + in.read(reinterpret_cast(values.data()), + num_values * sizeof(float)); + if (!in) return false; + } + + // Add range to sparse vector + signal_roi.add_range(offset, values); + } + + // Create wire + wires.emplace_back(std::move(signal_roi), channel, view); + } + + in.close(); + return in.good(); +} diff --git a/gaus_hit_finder/wire_serialization.hpp b/gaus_hit_finder/wire_serialization.hpp new file mode 100644 index 0000000..9b104a1 --- /dev/null +++ b/gaus_hit_finder/wire_serialization.hpp @@ -0,0 +1,15 @@ +#ifndef PHLEX_EXAMPLES_WIRE_SERIALIZATION_H +#define PHLEX_EXAMPLES_WIRE_SERIALIZATION_H + +#include +#include +#include +#include "copied_from_larsoft_minor_edits/Wire.h" + +bool write_wires_to_file(std::vector const& wires, + std::string const& filename); + +bool read_wires_from_file(std::string const& filename, + std::vector& wires); + +#endif // PHLEX_EXAMPLES_WIRE_SERIALIZATION_H diff --git a/gaus_hit_finder/wires_source.cpp b/gaus_hit_finder/wires_source.cpp new file mode 100644 index 0000000..a4d1e76 --- /dev/null +++ b/gaus_hit_finder/wires_source.cpp @@ -0,0 +1,40 @@ +#include +#include +#include +#include + +#include "copied_from_larsoft_minor_edits/Wire.h" + +#include "phlex/configuration.hpp" +#include "phlex/model/data_cell_index.hpp" +#include "phlex/source.hpp" +#include "wire_serialization.hpp" + +using namespace phlex; + +static std::atomic fileCounter{0}; + +// This provider is for temporary testing purposes only. This +// is NOT intended to be a long term solution for persistence. +// The purpose of this provider is to read input for +// GausHitFinder from files created in an "art" framework +// job, without using the art framework or ROOT. +// One must run an art job that produces these data +// files before running the phlex job that uses this provider. + +PHLEX_REGISTER_PROVIDERS(m, config) +{ + auto const layer = config.get("layer"); + + m.provide("provide_wires", [](data_cell_index const& id) -> std::vector { + int current_value = fileCounter.fetch_add(1); + std::string filename = std::string("wires_") + std::to_string(current_value) + ".dat"; + std::vector wires; + bool status = read_wires_from_file(filename, wires); + if (!status) { + throw std::runtime_error("Failure while reading from file: " + filename); + } + return wires; + }) + .output_product("wires"_in(layer)); +}