LArFFTW.h
Go to the documentation of this file.
1 #ifndef LARFFTW_H
2 #define LARFFTW_H
3 
4 // C/C++ standard libraries
5 #include <string>
6 #include <vector>
7 #include <complex>
8 #include <algorithm>
9 
10 #include "fftw3.h"
11 
12 #include "messagefacility/MessageLogger/MessageLogger.h"
13 #include "cetlib_except/coded_exception.h"
14 #include "larvecutils/MarqFitAlg/MarqFitAlg.h"
15 
16 namespace util {
17 
18 class LArFFTW {
19 
20  public:
21 
22  using FloatVector = std::vector<float>;
23  using DoubleVector = std::vector<double>;
24  using ComplexVector = std::vector<std::complex<double>>;
25 
26  LArFFTW(int transformSize, const void* fplan, const void* rplan, int fitbins);
27  ~LArFFTW();
28 
29  template <class T> void DoFFT(std::vector<T>& input);
30  template <class T> void DoFFT(std::vector<T>& input, ComplexVector& output);
31  template <class T> void DoInvFFT(std::vector<T>& output);
32  template <class T> void DoInvFFT(ComplexVector& input, std::vector<T>& output);
33 
34  // ... Do convolution calculation (for simulation).
35  template <class T> void Convolute(std::vector<T>& func, const ComplexVector& kern);
36  template <class T> void Convolute(std::vector<T>& func, std::vector<T>& resp);
37 
38  // ... Do deconvolution calculation (for reconstruction).
39  template <class T> void Deconvolute(std::vector<T>& func, const ComplexVector& kern);
40  template <class T> void Deconvolute(std::vector<T>& func, std::vector<T>& resp);
41 
42  // ... Do correlation
43  template <class T> void Correlate(std::vector<T>& func, const ComplexVector& kern);
44  template <class T> void Correlate(std::vector<T>& func, std::vector<T>& resp);
45 
46  void ShiftData(ComplexVector & input, double shift);
47  template <class T> void ShiftData(std::vector<T> & input, double shift);
48 
49  template <class T> void AlignedSum(std::vector<T> & input, std::vector<T> &output,
50  bool add = true);
51  template <class T> T PeakCorrelation(std::vector<T> &shape1,std::vector<T> &shape2);
52 
53  private:
54 
55  ComplexVector fKern; // transformed response function
56  ComplexVector fCompTemp; // temporary complex data
57  std::vector<float> fConvHist; // Fit data histogram
58  int fSize; // size of transform
59  int fFreqSize; // size of frequency space
60  void *fIn;
61  void *fOut;
62  const void *fPlan;
63  void *rIn;
64  void *rOut;
65  const void *rPlan;
66  int fFitBins; // Bins used for peak fit
67 
68  gshf::MarqFitAlg* fMarqFitAlg;
69 };
70 
71 } // end namespace util
72 
73 // -----------------------------------------------------------------------------
74 // ~~~~ Do Forward Fourier Transform - DoFFT( REAL In )
75 // -----------------------------------------------------------------------------
76 template <class T> inline void util::LArFFTW::DoFFT(std::vector<T> & input)
77 {
78  // ..set point
79  for(size_t p = 0; p < input.size(); ++p){
80  ((double *)fIn)[p] = input[p];
81  }
82 
83  // ..transform (using the New-array Execute Functions)
84  fftw_execute_dft_r2c((fftw_plan)fPlan,(double*)fIn,(fftw_complex*)fOut);
85 
86  return;
87 }
88 
89 // -----------------------------------------------------------------------------
90 // ~~~~ Do Forward Fourier Transform - DoFFT( REAL In, COMPLEX Out )
91 // -----------------------------------------------------------------------------
92 template <class T> inline void util::LArFFTW::DoFFT(std::vector<T> & input, ComplexVector& output)
93 {
94  // ..set point
95  for(size_t p = 0; p < input.size(); ++p){
96  ((double *)fIn)[p] = input[p];
97  }
98 
99  // ..transform (using the New-array Execute Functions)
100  fftw_execute_dft_r2c((fftw_plan)fPlan,(double*)fIn,(fftw_complex*)fOut);
101 
102  for(int i = 0; i < fFreqSize; ++i){
103  output[i].real(((fftw_complex*)fOut)[i][0]);
104  output[i].imag(((fftw_complex*)fOut)[i][1]);
105  }
106 
107  return;
108 }
109 
110 // -----------------------------------------------------------------------------
111 // ~~~~ Do Inverse Fourier Transform - DoInvFFT( REAL Out )
112 // -----------------------------------------------------------------------------
113 template <class T> inline void util::LArFFTW::DoInvFFT(std::vector<T> & output)
114 {
115  // ..transform (using the New-array Execute Functions)
116  fftw_execute_dft_c2r((fftw_plan)rPlan,(fftw_complex*)rIn,(double*)rOut);
117 
118  // ..get point real
119  double factor = 1.0/(double) fSize;
120  const double * array = (const double*)(rOut);
121  for(int i = 0; i < fSize; ++i){
122  output[i] = factor*array[i];
123  }
124 
125  return;
126 }
127 
128 // -----------------------------------------------------------------------------
129 // ~~~~ Do Inverse Fourier Transform - DoInvFFT( COMPLEX In, REAL Out )
130 // -----------------------------------------------------------------------------
131 template <class T> inline void util::LArFFTW::DoInvFFT(ComplexVector& input, std::vector<T> & output)
132 {
133  // ..set point complex
134  for(int i = 0; i < fFreqSize; ++i){
135  ((fftw_complex*)rIn)[i][0] = input[i].real();
136  ((fftw_complex*)rIn)[i][1] = input[i].imag();
137  }
138 
139  // ..transform (using the New-array Execute Functions)
140  fftw_execute_dft_c2r((fftw_plan)rPlan,(fftw_complex*)rIn,(double*)rOut);
141 
142  // ..get point real
143  double factor = 1.0/(double) fSize;
144  const double * array = (const double*)(rOut);
145  for(int i = 0; i < fSize; ++i){
146  output[i] = factor*array[i];
147  }
148 
149  return;
150 }
151 
152 // -----------------------------------------------------------------------------
153 // ~~~~ Do Convolution: using transformed response function
154 // -----------------------------------------------------------------------------
155 template <class T>
156 inline void util::LArFFTW::Convolute(std::vector<T>& func,
157  const ComplexVector& kern){
158 
159  // ... Make sure that time series and kernel have the correct size.
160  int n = func.size();
161  if(n != fSize){
162  throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
163  }
164  n = kern.size();
165  if(n != fFreqSize){
166  throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
167  }
168 
169  DoFFT(func);
170 
171  // ..perform the convolution
172  for(int i = 0; i < fFreqSize; ++i){
173  double re = ((fftw_complex*)fOut)[i][0];
174  double im = ((fftw_complex*)fOut)[i][1];
175  ((fftw_complex*)rIn)[i][0] = re*kern[i].real()-im*kern[i].imag();
176  ((fftw_complex*)rIn)[i][1] = re*kern[i].imag()+im*kern[i].real();
177  }
178 
179  DoInvFFT(func);
180 }
181 
182 // -----------------------------------------------------------------------------
183 // ~~~~ Do Convolution: using all time-domain information
184 // -----------------------------------------------------------------------------
185 template <class T>
186 inline void util::LArFFTW::Convolute(std::vector<T>& func1,
187  std::vector<T>& func2){
188 
189  // ... Make sure that time series has the correct size.
190  int n = func1.size();
191  if(n != fSize){
192  throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
193  }
194  n = func2.size();
195  if(n != fSize){
196  throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
197  }
198 
199  DoFFT(func2);
200  for(int i = 0; i < fFreqSize; ++i){
201  fKern[i].real(((fftw_complex*)fOut)[i][0]);
202  fKern[i].imag(((fftw_complex*)fOut)[i][1]);
203  }
204  DoFFT(func1);
205 
206  // ..perform the convolution
207  for(int i = 0; i < fFreqSize; ++i){
208  double re = ((fftw_complex*)fOut)[i][0];
209  double im = ((fftw_complex*)fOut)[i][1];
210  ((fftw_complex*)rIn)[i][0] = re*fKern[i].real()-im*fKern[i].imag();
211  ((fftw_complex*)rIn)[i][1] = re*fKern[i].imag()+im*fKern[i].real();
212  }
213 
214  DoInvFFT(func1);
215 }
216 
217 // -----------------------------------------------------------------------------
218 // ~~~~ Do Deconvolution: using transformed response function
219 // -----------------------------------------------------------------------------
220 template <class T>
221 inline void util::LArFFTW::Deconvolute(std::vector<T>& func,
222  const ComplexVector& kern){
223 
224  // ... Make sure that time series and kernel have the correct size.
225  int n = func.size();
226  if(n != fSize){
227  throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
228  }
229  n = kern.size();
230  if(n != fFreqSize){
231  throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
232  }
233 
234  DoFFT(func);
235 
236  // ..perform the deconvolution
237  double a,b,c,d,e;
238  for(int i = 0; i < fFreqSize; ++i){
239  a = ((fftw_complex*)fOut)[i][0];
240  b = ((fftw_complex*)fOut)[i][1];
241  c = kern[i].real();
242  d = kern[i].imag();
243  e = 1./(c*c+d*d);
244  ((fftw_complex*)rIn)[i][0] = (a*c+b*d)*e;
245  ((fftw_complex*)rIn)[i][1] = (b*c-a*d)*e;
246  }
247 
248  DoInvFFT(func);
249 }
250 
251 // -----------------------------------------------------------------------------
252 // ~~~~ Do Deconvolution: using all time domain information
253 // -----------------------------------------------------------------------------
254 template <class T>
255 inline void util::LArFFTW::Deconvolute(std::vector<T>& func,
256  std::vector<T>& resp){
257 
258  // ... Make sure that time series has the correct size.
259  int n = func.size();
260  if(n != fSize){
261  throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
262  }
263  n = resp.size();
264  if(n != fSize){
265  throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
266  }
267 
268  DoFFT(resp);
269  for(int i = 0; i < fFreqSize; ++i){
270  fKern[i].real(((fftw_complex*)fOut)[i][0]);
271  fKern[i].imag(((fftw_complex*)fOut)[i][1]);
272  }
273  DoFFT(func);
274 
275  // ..perform the deconvolution
276  double a,b,c,d,e;
277  for(int i = 0; i < fFreqSize; ++i){
278  a = ((fftw_complex*)fOut)[i][0];
279  b = ((fftw_complex*)fOut)[i][1];
280  c = fKern[i].real();
281  d = fKern[i].imag();
282  e = 1./(c*c+d*d);
283  ((fftw_complex*)rIn)[i][0] = (a*c+b*d)*e;
284  ((fftw_complex*)rIn)[i][1] = (b*c-a*d)*e;
285  }
286 
287  DoInvFFT(func);
288 
289 }
290 
291 // -----------------------------------------------------------------------------
292 // ~~~~ Do Deconvolution: using transformed response function
293 // -----------------------------------------------------------------------------
294 template <class T>
295 inline void util::LArFFTW::Correlate(std::vector<T>& func,
296  const ComplexVector& kern){
297 
298  // ... Make sure that time series and kernel have the correct size.
299  int n = func.size();
300  if(n != fSize){
301  throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
302  }
303  n = kern.size();
304  if(n != fFreqSize){
305  throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
306  }
307 
308  DoFFT(func);
309 
310  // ..perform the correlation
311  for(int i = 0; i < fFreqSize; ++i){
312  double re = ((fftw_complex*)fOut)[i][0];
313  double im = ((fftw_complex*)fOut)[i][1];
314  ((fftw_complex*)rIn)[i][0] = re*kern[i].real()+im*kern[i].imag();
315  ((fftw_complex*)rIn)[i][1] = -re*kern[i].imag()+im*kern[i].real();
316  }
317 
318  DoInvFFT(func);
319 
320 }
321 
322 // -----------------------------------------------------------------------------
323 // ~~~~ Do Correlation: using all time domain information
324 // -----------------------------------------------------------------------------
325 template <class T>
326 inline void util::LArFFTW::Correlate(std::vector<T>& func1,
327  std::vector<T>& func2){
328 
329  // ... Make sure that time series has the correct size.
330  int n = func1.size();
331  if(n != fSize){
332  throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
333  }
334  n = func2.size();
335  if(n != fSize){
336  throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
337  }
338 
339  DoFFT(func2);
340  for(int i = 0; i < fFreqSize; ++i){
341  fKern[i].real(((fftw_complex*)fOut)[i][0]);
342  fKern[i].imag(((fftw_complex*)fOut)[i][1]);
343  }
344  DoFFT(func1);
345 
346  // ..perform the correlation
347  for(int i = 0; i < fFreqSize; ++i){
348  double re = ((fftw_complex*)fOut)[i][0];
349  double im = ((fftw_complex*)fOut)[i][1];
350  ((fftw_complex*)rIn)[i][0] = re*fKern[i].real()+im*fKern[i].imag();
351  ((fftw_complex*)rIn)[i][1] = -re*fKern[i].imag()+im*fKern[i].real();
352  }
353 
354  DoInvFFT(func1);
355 
356 }
357 
358 // -----------------------------------------------------------------------------
359 // ~~~~ Shifts real vectors using above ShiftData function
360 // -----------------------------------------------------------------------------
361 template <class T>
362 inline void util::LArFFTW::ShiftData(std::vector<T> & input, double shift)
363 {
364  DoFFT(input,fCompTemp);
365  ShiftData(fCompTemp,shift);
366  DoInvFFT(fCompTemp,input);
367 
368  return;
369 }
370 
371 // -----------------------------------------------------------------------------
372 // ~~~~ Scheme for adding two signals which have an arbitrary relative
373 // translation. Shape1 is translated over shape2 and is replaced with the
374 // sum, or the translated result if add = false
375 // -----------------------------------------------------------------------------
376 template <class T> inline void util::LArFFTW::AlignedSum(std::vector<T> & shape1,
377  std::vector<T> & shape2,
378  bool add)
379 {
380  double shift = PeakCorrelation(shape1,shape2);
381 
382  ShiftData(shape1,shift);
383 
384  if(add)for(int i = 0; i < fSize; i++) shape1[i]+=shape2[i];
385 
386  return;
387 }
388 
389 // -----------------------------------------------------------------------------
390 // ~~~~ Returns the length of the translation at which the correlation
391 // of 2 signals is maximal.
392 // -----------------------------------------------------------------------------
393 template <class T> inline T util::LArFFTW::PeakCorrelation(std::vector<T> & shape1,
394  std::vector<T> & shape2)
395 {
396  float chiSqr = std::numeric_limits<float>::max();
397  float dchiSqr = std::numeric_limits<float>::max();
398  const float chiCut = 1e-3;
399  float lambda = 0.001; // Marquardt damping parameter
400  std::vector<float> p;
401 
402  std::vector<T> holder = shape1;
403  Correlate(holder,shape2);
404 
405  int maxT = max_element(holder.begin(), holder.end())-holder.begin();
406  float startT = maxT-fFitBins/2;
407  int offset = 0;
408 
409  for(int i = 0; i < fFitBins; i++) {
410  if(startT+i < 0) offset=fSize;
411  else if(startT+i > fSize) offset=-fSize;
412  else offset = 0;
413  if(holder[i+startT+offset]<=0.) {
414  fConvHist[i]=0.;
415  } else {
416  fConvHist[i]=holder[i+startT+offset];
417  }
418  }
419 
420  p[0] = *max_element(fConvHist.begin(), fConvHist.end());
421  p[1] = fFitBins/2;
422  p[2] = fFitBins/2;
423  float p1 = p[1]; // save initial p[1] guess
424 
425  int fitResult{-1};
426  int trial=0;
427  lambda=-1.; // initialize lambda on first call
428  do{
429  fitResult=fMarqFitAlg->mrqdtfit(lambda, &p[0], &fConvHist[0], 3, fFitBins, chiSqr, dchiSqr);
430  trial++;
431  if(fitResult){
432  mf::LogWarning("LArFFTW") << "Peak Correlation Fitting failed";
433  break;
434  }else if (trial>100){
435  break;
436  }
437  }
438  while (fabs(dchiSqr) >= chiCut);
439  if (!fitResult)p1=p[1]; // if fit succeeded, use fit result
440 
441  return p1 + 0.5 + startT;
442 }
443 #endif
util::LArFFTW::fKern
ComplexVector fKern
Definition: LArFFTW.h:55
util::LArFFTW::DoubleVector
std::vector< double > DoubleVector
Definition: LArFFTW.h:23
util
Namespace for general, non-LArSoft-specific utilities.
Definition: SignalShapingService.h:21
make_norm_csv.output
list output
Definition: make_norm_csv.py:48
util::LArFFTW::fIn
void * fIn
Definition: LArFFTW.h:60
util::LArFFTW::fFreqSize
int fFreqSize
Definition: LArFFTW.h:59
util::LArFFTW::fPlan
const void * fPlan
Definition: LArFFTW.h:62
util::LArFFTW::ComplexVector
std::vector< std::complex< double >> ComplexVector
Definition: LArFFTW.h:24
util::LArFFTW::LArFFTW
LArFFTW(int transformSize, const void *fplan, const void *rplan, int fitbins)
Definition: LArFFTW.cxx:5
util::LArFFTW::rOut
void * rOut
Definition: LArFFTW.h:64
add
Coord add(Coord c1, Coord c2)
Definition: restypedef.cpp:23
bump_copyright.d
string d
Definition: bump_copyright.py:70
util::LArFFTW::PeakCorrelation
T PeakCorrelation(std::vector< T > &shape1, std::vector< T > &shape2)
Definition: LArFFTW.h:393
util::LArFFTW::AlignedSum
void AlignedSum(std::vector< T > &input, std::vector< T > &output, bool add=true)
Definition: LArFFTW.h:376
util::LArFFTW::fConvHist
std::vector< float > fConvHist
Definition: LArFFTW.h:57
util::LArFFTW::Deconvolute
void Deconvolute(std::vector< T > &func, const ComplexVector &kern)
Definition: LArFFTW.h:221
util::LArFFTW::fFitBins
int fFitBins
Definition: LArFFTW.h:66
e
const double e
Definition: gUpMuFluxGen.cxx:165
input
static int input(void)
Definition: code.cpp:15695
util::LArFFTW::rIn
void * rIn
Definition: LArFFTW.h:63
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
a
const double a
Definition: gUpMuFluxGen.cxx:164
util::LArFFTW::ShiftData
void ShiftData(ComplexVector &input, double shift)
Definition: LArFFTW.cxx:46
test.p
p
Definition: test.py:223
util::LArFFTW::Convolute
void Convolute(std::vector< T > &func, const ComplexVector &kern)
Definition: LArFFTW.h:156
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
util::LArFFTW::rPlan
const void * rPlan
Definition: LArFFTW.h:65
util::LArFFTW::fOut
void * fOut
Definition: LArFFTW.h:61
func1
int func1()
Definition: Exception_t.cc:79
util::LArFFTW::Correlate
void Correlate(std::vector< T > &func, const ComplexVector &kern)
Definition: LArFFTW.h:295
gar::dump::array
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
util::LArFFTW::fCompTemp
ComplexVector fCompTemp
Definition: LArFFTW.h:56
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
func2
void func2()
Definition: group.cpp:82
docstring.func
def func()
Definition: docstring.py:7
util::LArFFTW::fMarqFitAlg
gshf::MarqFitAlg * fMarqFitAlg
Definition: LArFFTW.h:68
b
static bool * b
Definition: config.cpp:1043
util::LArFFTW::FloatVector
std::vector< float > FloatVector
Definition: LArFFTW.h:22
util::LArFFTW::DoInvFFT
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
util::LArFFTW::DoFFT
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Generated by   doxygen 1.8.11