All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
caldata::RawDigitFFTAlg< T > Class Template Reference

#include <RawDigitFFTAlg.h>

Public Member Functions

 RawDigitFFTAlg (fhicl::ParameterSet const &pset)
 
 ~RawDigitFFTAlg ()
 Destructor. More...
 
void reconfigure (fhicl::ParameterSet const &pset)
 
void initializeHists (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::ServiceHandle< art::TFileService > &)
 Begin job method. More...
 
void getFFTCorrection (std::vector< T > &, double) const
 
void getFFTCorrection (std::vector< T > &, size_t) const
 
void filterFFT (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< short > &, size_t, size_t, float pedestal=0.)
 

Private Attributes

std::vector< bool > fTransformViewVec
 apply FFT transform to this view More...
 
bool fFillHistograms
 if true then will fill diagnostic hists More...
 
std::string fHistDirName
 If writing histograms, the folder name. More...
 
std::vector< size_t > fLoWireByPlane
 Low wire for individual wire histograms. More...
 
std::vector< size_t > fHiWireByPlane
 Hi wire for individual wire histograms. More...
 
std::map< size_t, std::vector
< std::complex< float > > > 
fFilterVecMap
 
std::vector< float > fFFTInputVec
 
std::vector< std::complex
< float > > 
fFFTOutputVec
 
std::vector< float > fPowerVec
 
std::vector< std::vector
< TProfile * > > 
fFFTPowerVec
 
std::vector< TProfile * > fAveFFTPowerVec
 
std::vector< TProfile * > fConvFFTPowerVec
 
std::vector< TProfile * > fFilterFuncVec
 
icarus_signal_processing::WaveformTools
< T > 
fWaveformTool
 
std::map< size_t,
std::unique_ptr
< icarus_tool::IFilter > > 
fFilterToolMap
 
std::unique_ptr< Eigen::FFT
< float > > 
fEigenFFT
 

Detailed Description

template<class T>
class caldata::RawDigitFFTAlg< T >

Definition at line 44 of file RawDigitFFTAlg.h.

Constructor & Destructor Documentation

template<class T >
caldata::RawDigitFFTAlg< T >::RawDigitFFTAlg ( fhicl::ParameterSet const &  pset)

Constructor.

Arguments:

pset - Fcl parameters.

Definition at line 29 of file RawDigitFFTAlg.cxx.

30 {
31  reconfigure(pset);
32 
33  // Report.
34  mf::LogInfo("RawDigitFFTAlg") << "RawDigitFFTAlg configured\n";
35 }
void reconfigure(fhicl::ParameterSet const &pset)
template<class T >
caldata::RawDigitFFTAlg< T >::~RawDigitFFTAlg ( )

Destructor.

Definition at line 39 of file RawDigitFFTAlg.cxx.

40 {}

Member Function Documentation

template<class T >
void caldata::RawDigitFFTAlg< T >::filterFFT ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< short > &  rawadc,
size_t  plane,
size_t  wire,
float  pedestal = 0. 
)

Definition at line 210 of file RawDigitFFTAlg.cxx.

213 {
214  // Check there is something to do
215  if (!fTransformViewVec.at(plane)) return;
216 
217  // Step one is to setup and then get the FFT transform of the input waveform
218  size_t const fftDataSize = rawadc.size();
219  double sampleRate = sampling_rate(clockData);
220  double readOutSize = detProp.ReadOutWindowSize();
221 
222  if (fFFTInputVec.size() != fftDataSize) fFFTInputVec.resize(fftDataSize, 0.);
223 
224  std::transform(rawadc.begin(),rawadc.end(),fFFTInputVec.begin(),[pedestal](const auto& val){return float(float(val) - pedestal);});
225 
226  if (fFFTOutputVec.size() != fftDataSize) fFFTOutputVec.resize(fftDataSize);
227 
229 
230  size_t halfFFTDataSize(fftDataSize/2 + 1);
231 
232  if (fPowerVec.size() != halfFFTDataSize + 1) fPowerVec.resize(halfFFTDataSize + 1, 0.);
233 
234  std::transform(fFFTOutputVec.begin(), fFFTOutputVec.begin() + halfFFTDataSize, fPowerVec.begin(), [](const auto& complex){return std::abs(complex);});
235 
236  // Recover the filter function we are using...
237  const icarusutil::FrequencyVec& filter = fFilterToolMap.at(plane)->getResponseVec();
238  const std::vector<std::complex<float>>& filterVec = fFilterVecMap.at(plane);;
239 
240  // Make sure the filter has been correctly initialized
241  if (filter.size() != halfFFTDataSize)
242  {
243  fFilterToolMap.at(plane)->setResponse(fftDataSize,1.,1.);
244 
245  // Set up the internal FFT vector
246  fFilterVecMap.at(plane).reserve(halfFFTDataSize);
247  for(auto& complex : filter) fFilterVecMap.at(plane).emplace_back(complex);
248  }
249 
250  //std::transform(fFFTOutputVec.begin(), fFFTOutputVec.begin() + fFFTOutputVec.size()/2, filterVec.begin(), fFFTOutputVec.begin(), std::multiplies<std::complex<float>>());
251 
252  //for(size_t idx = 0; idx < fFFTOutputVec.size()/2; idx++) fFFTOutputVec[fFFTOutputVec.size() - idx - 1] = fFFTOutputVec[idx];
253 
254  std::transform(fFFTOutputVec.begin(), fFFTOutputVec.begin() + halfFFTDataSize, filterVec.begin(), fFFTOutputVec.begin(), std::multiplies<std::complex<float>>());
255 
256  for(size_t idx = 1; idx < fFFTOutputVec.size()/2; idx++) fFFTOutputVec[fFFTOutputVec.size() - idx] = fFFTOutputVec[idx];
257 
259 
260  // Fill hists
261  if (fFillHistograms)
262  {
263  // Fill any individual wire histograms we want to look at
264  if (wire >= fLoWireByPlane[plane] && wire < fHiWireByPlane[plane])
265  {
266  // Fill the power spectrum histogram
267  for(size_t idx = 0; idx < halfFFTDataSize; idx++)
268  {
269  float freq = 1.e6 * float(idx + 1)/ (sampleRate * readOutSize);
270  fFFTPowerVec[plane][wire-fLoWireByPlane[plane]]->Fill(freq, std::min(fPowerVec[idx],float(999.)), 1.);
271  }
272  }
273 
274  for(size_t idx = 0; idx < halfFFTDataSize; idx++)
275  {
276  float freq = 1.e6 * float(idx + 1)/ (sampleRate * readOutSize);
277  fAveFFTPowerVec[plane]->Fill(freq, std::min(fPowerVec[idx],float(999.)), 1.);
278  }
279 
280  // Get the filter power vec
281  std::transform(fFFTOutputVec.begin(), fFFTOutputVec.begin() + halfFFTDataSize, fPowerVec.begin(), [](const auto& val){return std::abs(val);});
282 
283  for(size_t idx = 0; idx < halfFFTDataSize; idx++)
284  {
285  float freq = 1.e6 * float(idx)/ (sampleRate * readOutSize);
286  fConvFFTPowerVec[plane]->Fill(freq, std::min(fPowerVec[idx],float(999.)), 1.);
287  fFilterFuncVec[plane]->Fill(freq, double(std::abs(filter[idx])), 1.);
288  }
289  }
290 
291  // Finally, we invert the resulting time domain values to recover the new waveform
292  std::transform(fFFTInputVec.begin(), fFFTInputVec.end(), rawadc.begin(), [pedestal](const float& adc){return std::round(adc + pedestal);});
293 
294  return;
295 }
static constexpr Sample_t transform(Sample_t sample)
std::map< size_t, std::vector< std::complex< float > > > fFilterVecMap
std::vector< size_t > fLoWireByPlane
Low wire for individual wire histograms.
std::vector< ComplexVal > FrequencyVec
std::vector< TProfile * > fAveFFTPowerVec
std::vector< TProfile * > fFilterFuncVec
std::vector< size_t > fHiWireByPlane
Hi wire for individual wire histograms.
T abs(T value)
std::vector< float > fPowerVec
bool fFillHistograms
if true then will fill diagnostic hists
std::vector< std::complex< float > > fFFTOutputVec
std::vector< std::vector< TProfile * > > fFFTPowerVec
std::unique_ptr< Eigen::FFT< float > > fEigenFFT
physics filters filter
std::vector< bool > fTransformViewVec
apply FFT transform to this view
std::vector< TProfile * > fConvFFTPowerVec
std::map< size_t, std::unique_ptr< icarus_tool::IFilter > > fFilterToolMap
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
std::vector< float > fFFTInputVec
auto const detProp
template<class T >
void caldata::RawDigitFFTAlg< T >::getFFTCorrection ( std::vector< T > &  corValVec,
double  minPowerThreshold 
) const

Definition at line 129 of file RawDigitFFTAlg.cxx.

130 {
131  // This version will take FFT of input waveform and then remove bins in the time domain with a power less
132  // than the threshold input above.
133  size_t const fftDataSize = corValVec.size();
134 
135  Eigen::FFT<T> eigenFFT;
136 
137  std::vector<std::complex<T>> fftOutputVec(corValVec.size());
138 
139  eigenFFT.fwd(fftOutputVec, corValVec);
140 
141  size_t halfFFTDataSize(fftDataSize/2 + 1);
142 
143  std::vector<T> powerVec(halfFFTDataSize);
144 
145  std::transform(fftOutputVec.begin(), fftOutputVec.begin() + halfFFTDataSize, powerVec.begin(), [](const auto& val){return std::abs(val);});
146 
147  // Want the first derivative
148  std::vector<T> firstDerivVec(powerVec.size());
149 
150  fWaveformTool.firstDerivative(powerVec, firstDerivVec);
151 
152  // Find the peaks
153  icarus_signal_processing::WaveformTools<float>::PeakTupleVec peakTupleVec;
154 
155  fWaveformTool.findPeaks(firstDerivVec.begin(),firstDerivVec.end(),peakTupleVec,minPowerThreshold,0);
156 
157  if (!peakTupleVec.empty())
158  {
159  for(const auto& peakTuple : peakTupleVec)
160  {
161  size_t startTick = std::get<0>(peakTuple);
162  size_t stopTick = std::get<2>(peakTuple);
163 
164  if (stopTick > startTick)
165  {
166  std::complex<T> slope = (fftOutputVec[stopTick] - fftOutputVec[startTick]) / T(stopTick - startTick);
167 
168  for(size_t tick = startTick; tick < stopTick; tick++)
169  {
170  std::complex<T> interpVal = fftOutputVec[startTick] + T(tick - startTick) * slope;
171 
172  fftOutputVec[tick] = interpVal;
173  //fftOutputVec[fftDataSize - tick - 1] = interpVal;
174  }
175  }
176  }
177 
178  std::vector<T> tmpVec(corValVec.size());
179 
180  eigenFFT.inv(tmpVec, fftOutputVec);
181 
182  std::transform(corValVec.begin(),corValVec.end(),tmpVec.begin(),corValVec.begin(),std::minus<T>());
183  }
184 
185  return;
186 }
static constexpr Sample_t transform(Sample_t sample)
icarus_signal_processing::WaveformTools< T > fWaveformTool
T abs(T value)
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
template<class T >
void caldata::RawDigitFFTAlg< T >::getFFTCorrection ( std::vector< T > &  corValVec,
size_t  maxBin 
) const

Definition at line 188 of file RawDigitFFTAlg.cxx.

189 {
190  // This version will take FFT of input waveform and then remove bins in the time domain above the
191  // cutoff frequency defined by maxBin passed in above
192  size_t const fftDataSize = corValVec.size();
193 
194  Eigen::FFT<T> eigenFFT;
195 
196  std::vector<std::complex<T>> fftOutputVec(corValVec.size());
197 
198  eigenFFT.fwd(fftOutputVec, corValVec);
199 
200  size_t halfFFTDataSize(fftDataSize/2);
201 
202  std::fill(fftOutputVec.begin() + maxBin, fftOutputVec.begin() + halfFFTDataSize, std::complex<T>(0.,0.));
203  std::fill(fftOutputVec.end() - maxBin - 1, fftOutputVec.end(), std::complex<T>(0.,0.));
204 
205  eigenFFT.inv(corValVec, fftOutputVec);
206 
207  return;
208 }
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
template<class T >
void caldata::RawDigitFFTAlg< T >::initializeHists ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
art::ServiceHandle< art::TFileService > &  tfs 
)

Begin job method.

Definition at line 74 of file RawDigitFFTAlg.cxx.

77 {
78  if (fFillHistograms)
79  {
80  // Define the histograms. Putting semi-colons around the title
81  // causes it to be displayed as the x-axis label if the histogram
82  // is drawn.
83 
84  // hijack hists here
85  double sampleRate = sampling_rate(clockData);
86  double readOutSize = detProp.ReadOutWindowSize();
87  double maxFreq = 1.e6 / (2. * sampleRate);
88  double minFreq = 1.e6 / (2. * sampleRate * readOutSize);
89  int numSamples = readOutSize / 2;
90 
91  // Make a directory for these histograms
92  art::TFileDirectory dir = tfs->mkdir(fHistDirName.c_str());
93 
94  fFFTPowerVec.resize(3);
95  fAveFFTPowerVec.resize(3);
96  fConvFFTPowerVec.resize(3);
97  fFilterFuncVec.resize(3);
98 
99  for(size_t plane = 0; plane < 3; plane++)
100  {
101  size_t numHists = fHiWireByPlane[plane] - fLoWireByPlane[plane];
102 
103  fFFTPowerVec[plane].resize(numHists);
104 
105  for(size_t idx = 0; idx < fHiWireByPlane[plane] - fLoWireByPlane[plane]; idx++)
106  {
107  std::string histName = "FFTPower_" + std::to_string(plane) + "-" + std::to_string(idx);
108 
109  fFFTPowerVec[plane][idx] = dir.make<TProfile>(histName.c_str(), "Power Spectrum;kHz;Power", numSamples, minFreq, maxFreq, 0., 10000.);
110  }
111 
112  std::string histName = "AveFFTPower_" + std::to_string(plane);
113 
114  fAveFFTPowerVec[plane] = dir.make<TProfile>(histName.c_str(), "Power Spectrum;kHz;Power", numSamples, minFreq, maxFreq, 0., 1000.);
115 
116  histName = "ConvFFTPower_" + std::to_string(plane);
117 
118  fConvFFTPowerVec[plane] = dir.make<TProfile>(histName.c_str(), "Power Spectrum;kHz;Power", numSamples, minFreq, maxFreq, 0., 1000.);
119 
120  histName = "FilterFunc_" + std::to_string(plane);
121 
122  fFilterFuncVec[plane] = dir.make<TProfile>(histName.c_str(), "Filter Function;kHz;Power", numSamples, minFreq, maxFreq, 0., 1000.);
123  }
124  }
125 
126  return;
127 }
std::vector< size_t > fLoWireByPlane
Low wire for individual wire histograms.
std::vector< TProfile * > fAveFFTPowerVec
std::vector< TProfile * > fFilterFuncVec
std::vector< size_t > fHiWireByPlane
Hi wire for individual wire histograms.
bool fFillHistograms
if true then will fill diagnostic hists
std::vector< std::vector< TProfile * > > fFFTPowerVec
std::string fHistDirName
If writing histograms, the folder name.
tuple dir
Definition: dropbox.py:28
std::string to_string(WindowPattern const &pattern)
std::vector< TProfile * > fConvFFTPowerVec
art::ServiceHandle< art::TFileService > tfs
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
auto const detProp
template<class T >
void caldata::RawDigitFFTAlg< T >::reconfigure ( fhicl::ParameterSet const &  pset)

Reconfigure method.

Arguments:

pset - Fcl parameter set.

Definition at line 49 of file RawDigitFFTAlg.cxx.

50 {
51  fTransformViewVec = pset.get<std::vector<bool> >("TransformViewVec", std::vector<bool>() = {true,false,false});
52  fFillHistograms = pset.get<bool >("FillHistograms", false);
53  fHistDirName = pset.get<std::string >("HistDirName", "FFT_hists");
54  fLoWireByPlane = pset.get<std::vector<size_t>>("LoWireByPlane", std::vector<size_t>()={0,0,0});
55  fHiWireByPlane = pset.get<std::vector<size_t>>("HiWireByPlane", std::vector<size_t>()={100,100,100});
56 
57  // Implement the tools for handling the responses
58  const fhicl::ParameterSet& filterTools = pset.get<fhicl::ParameterSet>("FilterTools");
59 
60  for(const std::string& filterTool : filterTools.get_pset_names())
61  {
62  const fhicl::ParameterSet& filterToolParamSet = filterTools.get<fhicl::ParameterSet>(filterTool);
63  size_t planeIdx = filterToolParamSet.get<size_t>("Plane");
64 
65  fFilterToolMap.insert(std::pair<size_t,std::unique_ptr<icarus_tool::IFilter>>(planeIdx,art::make_tool<icarus_tool::IFilter>(filterToolParamSet)));
66  fFilterVecMap[planeIdx] = std::vector<std::complex<float>>();
67  }
68 
69  fEigenFFT = std::make_unique<Eigen::FFT<float>>();
70 }
std::map< size_t, std::vector< std::complex< float > > > fFilterVecMap
std::vector< size_t > fLoWireByPlane
Low wire for individual wire histograms.
std::vector< size_t > fHiWireByPlane
Hi wire for individual wire histograms.
bool fFillHistograms
if true then will fill diagnostic hists
std::unique_ptr< Eigen::FFT< float > > fEigenFFT
std::string fHistDirName
If writing histograms, the folder name.
std::vector< bool > fTransformViewVec
apply FFT transform to this view
std::map< size_t, std::unique_ptr< icarus_tool::IFilter > > fFilterToolMap

Member Data Documentation

template<class T >
std::vector<TProfile*> caldata::RawDigitFFTAlg< T >::fAveFFTPowerVec
private

Definition at line 81 of file RawDigitFFTAlg.h.

template<class T >
std::vector<TProfile*> caldata::RawDigitFFTAlg< T >::fConvFFTPowerVec
private

Definition at line 82 of file RawDigitFFTAlg.h.

template<class T >
std::unique_ptr<Eigen::FFT<float> > caldata::RawDigitFFTAlg< T >::fEigenFFT
private

Definition at line 88 of file RawDigitFFTAlg.h.

template<class T >
std::vector<float> caldata::RawDigitFFTAlg< T >::fFFTInputVec
private

Definition at line 76 of file RawDigitFFTAlg.h.

template<class T >
std::vector<std::complex<float> > caldata::RawDigitFFTAlg< T >::fFFTOutputVec
private

Definition at line 77 of file RawDigitFFTAlg.h.

template<class T >
std::vector<std::vector<TProfile*> > caldata::RawDigitFFTAlg< T >::fFFTPowerVec
private

Definition at line 80 of file RawDigitFFTAlg.h.

template<class T >
bool caldata::RawDigitFFTAlg< T >::fFillHistograms
private

if true then will fill diagnostic hists

Definition at line 69 of file RawDigitFFTAlg.h.

template<class T >
std::vector<TProfile*> caldata::RawDigitFFTAlg< T >::fFilterFuncVec
private

Definition at line 83 of file RawDigitFFTAlg.h.

template<class T >
std::map<size_t,std::unique_ptr<icarus_tool::IFilter> > caldata::RawDigitFFTAlg< T >::fFilterToolMap
private

Definition at line 86 of file RawDigitFFTAlg.h.

template<class T >
std::map<size_t,std::vector<std::complex<float> > > caldata::RawDigitFFTAlg< T >::fFilterVecMap
private

Definition at line 75 of file RawDigitFFTAlg.h.

template<class T >
std::string caldata::RawDigitFFTAlg< T >::fHistDirName
private

If writing histograms, the folder name.

Definition at line 70 of file RawDigitFFTAlg.h.

template<class T >
std::vector<size_t> caldata::RawDigitFFTAlg< T >::fHiWireByPlane
private

Hi wire for individual wire histograms.

Definition at line 72 of file RawDigitFFTAlg.h.

template<class T >
std::vector<size_t> caldata::RawDigitFFTAlg< T >::fLoWireByPlane
private

Low wire for individual wire histograms.

Definition at line 71 of file RawDigitFFTAlg.h.

template<class T >
std::vector<float> caldata::RawDigitFFTAlg< T >::fPowerVec
private

Definition at line 78 of file RawDigitFFTAlg.h.

template<class T >
std::vector<bool> caldata::RawDigitFFTAlg< T >::fTransformViewVec
private

apply FFT transform to this view

Definition at line 68 of file RawDigitFFTAlg.h.

template<class T >
icarus_signal_processing::WaveformTools<T> caldata::RawDigitFFTAlg< T >::fWaveformTool
private

Definition at line 85 of file RawDigitFFTAlg.h.


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