LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
SignalShaping.cxx
Go to the documentation of this file.
1 
12 #include "cetlib_except/exception.h"
13 #include <cmath>
14 
15 //----------------------------------------------------------------------
16 // Constructor.
17 //
18 util::SignalShaping::SignalShaping() : fResponseLocked(false), fFilterLocked(false), fNorm(true) {}
19 
20 //----------------------------------------------------------------------
21 // Destructor.
22 //
24 
25 // void util::SignalShaping::ResetDecon()
26 // {
27 // fResponseLocked = false;
28 // fFilterLocked = false;
29 // fResponse.clear();
30 // fConvKernel.clear();
31 // fFilter.clear();
32 // fDeconvKernel.clear();
33 // //Set deconvolution polarity to + as default
34 // fDeconvKernelPolarity = +1;
35 // }
36 
37 //----------------------------------------------------------------------
38 // Reset this class to its default-constructed state.
40 {
41  fResponseLocked = false;
42  fFilterLocked = false;
43  fResponse.clear();
44  fConvKernel.clear();
45  fFilter.clear();
46  fDeconvKernel.clear();
47  //Set deconvolution polarity to + as default
49 }
50 
51 //----------------------------------------------------------------------
52 // Add a time domain response function.
53 void util::SignalShaping::AddResponseFunction(const std::vector<double>& resp, bool ResetResponse)
54 {
55  // Make sure configuration is not locked.
56 
57  if (fResponseLocked) throw cet::exception("SignalShaping") << "Configuration locked.\n";
58 
59  // Get FFT service.
60 
62  int nticks = fft->FFTSize();
63 
64  // Copy new response function into fResponse attribute, and pad or
65  // truncate to correct size.
66 
67  fResponse = resp;
68  fResponse.resize(nticks, 0.);
69 
70  // Is this the first response function?
71 
72  if (fConvKernel.size() == 0 || ResetResponse) {
73 
74  // This is the first response function.
75  // Just calculate the fourier transform.
76 
77  fConvKernel.resize(nticks / 2 + 1);
79  }
80  else {
81 
82  // Not the first response function.
83  // Calculate the fourier transform of new response function.
84 
85  std::vector<TComplex> kern(nticks / 2 + 1);
86  fft->DoFFT(fResponse, kern);
87 
88  // Update overall convolution kernel.
89 
90  if (kern.size() != fConvKernel.size()) {
91  throw cet::exception("SignalShaping") << __func__ << ": inconsistent kernel size, "
92  << kern.size() << " vs. " << fConvKernel.size();
93  }
94  for (unsigned int i = 0; i < kern.size(); ++i)
95  fConvKernel[i] *= kern[i];
96 
97  // Recalculate overall response function.
98 
100  }
101 }
102 
103 //----------------------------------------------------------------------
104 // Shift the response function and convolution kernel by the specified
105 // number of ticks.
107 {
108  // Make sure configuration is not locked.
109 
110  if (fResponseLocked) throw cet::exception("SignalShaping") << "Configuration locked.\n";
111 
112  // Get FFT service.
113 
115 
116  // Update convolution kernel.
117 
118  fft->ShiftData(fConvKernel, ticks);
119 
120  // Recalculate overall response functiion.
121 
123 }
124 
125 //----------------------------------------------------------------------
126 // Set the peak response time to be at the specified tick.
128 {
129  // Make sure configuration is not locked.
130 
131  if (fResponseLocked) throw cet::exception("SignalShaping") << "Configuration locked.\n";
132 
133  // Get FFT service.
134 
136 
137  // Construct a delta-function response centered at tick zero.
138 
139  std::vector<double> delta(fft->FFTSize(), 0.);
140  delta[0] = 1.;
141 
142  // Figure out peak of current overall response.
143 
144  double peak = fft->PeakCorrelation(delta, fResponse);
145 
146  // Shift peak response to desired tick.
147 
148  ShiftResponseTime(tick - peak);
149 }
150 
151 //----------------------------------------------------------------------
152 // Add a frequency domain filter function to cumulative filter function.
153 void util::SignalShaping::AddFilterFunction(const std::vector<TComplex>& filt)
154 {
155  // Make sure configuration is not locked.
156 
157  if (fFilterLocked) throw cet::exception("SignalShaping") << "Configuration locked.\n";
158 
159  // Get FFT service.
160 
162 
163  // If this is the first filter function, just copy the filter function.
164  // Otherwise, update the overall filter function.
165 
166  if (fFilter.size() == 0) {
167  fFilter = filt;
168  fFilter.resize(fft->FFTSize() / 2 + 1);
169  }
170  else {
171  unsigned int n = std::min(fFilter.size(), filt.size());
172  for (unsigned int i = 0; i < n; ++i)
173  fFilter[i] *= filt[i];
174  for (unsigned int i = n; i < fFilter.size(); ++i)
175  fFilter[i] = 0.;
176  }
177 }
178 
179 //----------------------------------------------------------------------
180 // Add a DeconvKernel Polarity Flag to decide how to normalize
182 {
183 
184  if ((pol != 1) and (pol != -1)) {
185  throw cet::exception("SignalShaping")
186  << __func__ << ": DeconvKernelPolarity should be +1 or -1 (got " << pol
187  << "). Setting to +1\n";
189  }
190 
191  else
192  fDeconvKernelPolarity = pol;
193 }
194 
195 //----------------------------------------------------------------------
196 // Test and lock the response and convolution kernel.
198 {
199  // Do nothing if the response is already locked.
200 
201  if (!fResponseLocked) {
202 
203  // Get FFT service.
204 
206 
207  // Make sure response has been configured.
208 
209  if (fResponse.size() == 0)
210  throw cet::exception("SignalShaping") << "Response has not been configured.\n";
211 
212  // Make sure response and convolution kernel have the correct
213  // size (should always be the case if we get here).
214 
215  unsigned int n = fft->FFTSize();
216  if (fResponse.size() != n)
217  throw cet::exception("SignalShaping")
218  << __func__ << ": inconsistent kernel size, " << fResponse.size() << " vs. " << n << "\n";
219  if (2 * (fConvKernel.size() - 1) != n)
220  throw cet::exception("SignalShaping")
221  << __func__ << ": unexpected FFT size, " << n << " vs. expected "
222  << (2 * (fConvKernel.size() - 1)) << "\n";
223 
224  // Set the lock flag.
225 
226  fResponseLocked = true;
227  }
228 }
229 
230 //----------------------------------------------------------------------
231 // Calculate the deconvolution kernel as the ratio
232 // of the filter function and convolution kernel.
234 {
235  // Make sure configuration is not locked.
236 
237  if (fFilterLocked) throw cet::exception("SignalShaping") << "Configuration locked.\n";
238 
239  // Lock response configuration.
240 
241  LockResponse();
242 
243  // Get FFT service.
244 
246 
247  // Make sure filter function has been configured.
248 
249  if (fFilter.size() == 0)
250  throw cet::exception("SignalShaping") << "Filter function has not been configured.\n";
251 
252  // Make sure filter function has the correct size.
253  // (Should always be the case if we get here.)
254 
255  unsigned int n = fft->FFTSize();
256  if (2 * (fFilter.size() - 1) != n)
257  if (fFilter.size() != fConvKernel.size()) {
258  throw cet::exception("SignalShaping") << __func__ << ": inconsistent size, " << fFilter.size()
259  << " vs. " << fConvKernel.size() << "\n";
260  }
261 
262  // Calculate deconvolution kernel as the ratio of the
263  // filter function and the convolution kernel.
264 
266  for (unsigned int i = 0; i < fDeconvKernel.size(); ++i) {
267  if (std::abs(fConvKernel[i].Re()) <= 0.0001 && std::abs(fConvKernel[i].Im()) <= 0.0001) {
268  fDeconvKernel[i] = 0.;
269  }
270  else {
271  fDeconvKernel[i] /= fConvKernel[i];
272  }
273  }
274 
275  // Normalize the deconvolution kernel.
276 
277  // Calculate the unnormalized deconvoluted response
278  // (inverse FFT of filter function).
279 
280  std::vector<double> deconv(n, 0.);
281  fft->DoInvFFT(const_cast<std::vector<TComplex>&>(fFilter), deconv);
282 
283  if (fNorm) {
284  // Find the peak value of the response
285  // Should normally be at zero, but don't assume that.
286  // Use DeconvKernelPolairty to find what to normalize to
287  double peak_response = 0;
288  if (fDeconvKernelPolarity == -1) peak_response = 4096;
289  for (unsigned int i = 0; i < fResponse.size(); ++i) {
290  if ((fResponse[i] > peak_response) and (fDeconvKernelPolarity == 1))
291  peak_response = fResponse[i];
292  else if ((fResponse[i] < peak_response) and (fDeconvKernelPolarity == -1))
293  peak_response = fResponse[i];
294  }
295  if (fDeconvKernelPolarity == -1) peak_response *= -1;
296  if (peak_response <= 0.) {
297  throw cet::exception("SignalShaping")
298  << __func__ << ": peak should always be positive (got " << peak_response << ")\n";
299  }
300 
301  // Find the peak value of the deconvoluted response
302  // Should normally be at zero, but don't assume that.
303 
304  double peak_deconv = 0.;
305  for (unsigned int i = 0; i < deconv.size(); ++i) {
306  if (deconv[i] > peak_deconv) peak_deconv = deconv[i];
307  }
308  if (peak_deconv <= 0.) {
309  throw cet::exception("SignalShaping")
310  << __func__ << ": deconvolution peak should always be positive (got " << peak_deconv
311  << ")\n";
312  }
313 
314  // Multiply the deconvolution kernel by a factor such that
315  // (Peak of response) = (Peak of deconvoluted response).
316 
317  double ratio = peak_response / peak_deconv;
318  for (unsigned int i = 0; i < fDeconvKernel.size(); ++i)
319  fDeconvKernel[i] *= ratio;
320  }
321  // Set the lock flag.
322 
323  fFilterLocked = true;
324 }
void ShiftData(std::vector< TComplex > &input, double shift)
Definition: LArFFT.cc:116
std::vector< TComplex > fConvKernel
constexpr auto abs(T v)
Returns the absolute value of the argument.
void AddResponseFunction(const std::vector< double > &resp, bool ResetResponse=false)
void LockResponse() const
std::vector< TComplex > fFilter
void DoFFT(std::vector< T > &input, std::vector< TComplex > &output)
Definition: LArFFT.h:95
T PeakCorrelation(std::vector< T > &shape1, std::vector< T > &shape2)
Definition: LArFFT.h:270
tick ticks
Alias for common language habits.
Definition: electronics.h:76
void DoInvFFT(std::vector< TComplex > &input, std::vector< T > &output)
Definition: LArFFT.h:117
Generic class for shaping signals on wires.
int FFTSize() const
Definition: LArFFT.h:66
void CalculateDeconvKernel() const
void SetDeconvKernelPolarity(int pol)
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
void ShiftResponseTime(double ticks)
std::vector< TComplex > fDeconvKernel
void SetPeakResponseTime(double tick)
std::vector< double > fResponse
Char_t n[5]
TString filt[ntarg]
Definition: Style.C:28
void AddFilterFunction(const std::vector< TComplex > &filt)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33