All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Types | Public Member Functions | Private Attributes | List of all members
util::LArFFTW Class Reference

#include <LArFFTW.h>

Public Types

using FloatVector = std::vector< float >
 
using DoubleVector = std::vector< double >
 
using ComplexVector = std::vector< std::complex< double >>
 

Public Member Functions

 LArFFTW (int transformSize, const void *fplan, const void *rplan, int fitbins)
 
 ~LArFFTW ()
 
template<class T >
void DoFFT (std::vector< T > &input)
 
template<class T >
void DoFFT (std::vector< T > &input, ComplexVector &output)
 
template<class T >
void DoInvFFT (std::vector< T > &output)
 
template<class T >
void DoInvFFT (ComplexVector &input, std::vector< T > &output)
 
template<class T >
void Convolute (std::vector< T > &func, const ComplexVector &kern)
 
template<class T >
void Convolute (std::vector< T > &func, std::vector< T > &resp)
 
template<class T >
void Deconvolute (std::vector< T > &func, const ComplexVector &kern)
 
template<class T >
void Deconvolute (std::vector< T > &func, std::vector< T > &resp)
 
template<class T >
void Correlate (std::vector< T > &func, const ComplexVector &kern)
 
template<class T >
void Correlate (std::vector< T > &func, std::vector< T > &resp)
 
void ShiftData (ComplexVector &input, double shift)
 
template<class T >
void ShiftData (std::vector< T > &input, double shift)
 
template<class T >
void AlignedSum (std::vector< T > &input, std::vector< T > &output, bool add=true)
 
template<class T >
PeakCorrelation (std::vector< T > &shape1, std::vector< T > &shape2)
 

Private Attributes

ComplexVector fKern
 
ComplexVector fCompTemp
 
std::vector< float > fConvHist
 
int fSize
 
int fFreqSize
 
void * fIn
 
void * fOut
 
const void * fPlan
 
void * rIn
 
void * rOut
 
const void * rPlan
 
int fFitBins
 
gshf::MarqFitAlg * fMarqFitAlg
 

Detailed Description

Definition at line 18 of file LArFFTW.h.

Member Typedef Documentation

using util::LArFFTW::ComplexVector = std::vector<std::complex<double>>

Definition at line 24 of file LArFFTW.h.

using util::LArFFTW::DoubleVector = std::vector<double>

Definition at line 23 of file LArFFTW.h.

using util::LArFFTW::FloatVector = std::vector<float>

Definition at line 22 of file LArFFTW.h.

Constructor & Destructor Documentation

util::LArFFTW::LArFFTW ( int  transformSize,
const void *  fplan,
const void *  rplan,
int  fitbins 
)

Definition at line 5 of file LArFFTW.cxx.

6  : fSize (transformSize)
7  , fPlan (fplan)
8  , rPlan (rplan)
9  , fFitBins (fitbins)
10 {
11 
12  fFreqSize = fSize/2+1;
13 
14  // ... Real-Complex
15  fIn = fftw_malloc(sizeof(double)*fSize);
16  fOut= fftw_malloc(sizeof(fftw_complex)*fFreqSize);
17 
18  // ... Complex-Real
19  rIn = fftw_malloc(sizeof(fftw_complex)*fFreqSize);
20  rOut= fftw_malloc(sizeof(double)*fSize);
21 
22  // ... allocate other data vectors
23  fCompTemp.resize(fFreqSize);
24  fKern.resize(fFreqSize);
25  fConvHist.resize(fFitBins);
26 }
ComplexVector fKern
Definition: LArFFTW.h:55
void * fIn
Definition: LArFFTW.h:60
int fFreqSize
Definition: LArFFTW.h:59
const void * fPlan
Definition: LArFFTW.h:62
void * rOut
Definition: LArFFTW.h:64
std::vector< float > fConvHist
Definition: LArFFTW.h:57
int fFitBins
Definition: LArFFTW.h:66
void * rIn
Definition: LArFFTW.h:63
const void * rPlan
Definition: LArFFTW.h:65
void * fOut
Definition: LArFFTW.h:61
ComplexVector fCompTemp
Definition: LArFFTW.h:56
util::LArFFTW::~LArFFTW ( )

Definition at line 28 of file LArFFTW.cxx.

29 {
30  fPlan = 0;
31  fftw_free(fIn);
32  fIn = 0;
33  fftw_free((fftw_complex*)fOut);
34  fOut = 0;
35 
36  rPlan = 0;
37  fftw_free((fftw_complex*)rIn);
38  rIn = 0;
39  fftw_free(rOut);
40  rOut = 0;
41 }
void * fIn
Definition: LArFFTW.h:60
const void * fPlan
Definition: LArFFTW.h:62
void * rOut
Definition: LArFFTW.h:64
void * rIn
Definition: LArFFTW.h:63
const void * rPlan
Definition: LArFFTW.h:65
void * fOut
Definition: LArFFTW.h:61

Member Function Documentation

template<class T >
void util::LArFFTW::AlignedSum ( std::vector< T > &  input,
std::vector< T > &  output,
bool  add = true 
)
inline

Definition at line 376 of file LArFFTW.h.

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 }
T PeakCorrelation(std::vector< T > &shape1, std::vector< T > &shape2)
Definition: LArFFTW.h:393
shift
Definition: fcl_checks.sh:26
void ShiftData(ComplexVector &input, double shift)
Definition: LArFFTW.cxx:46
template<class T >
void util::LArFFTW::Convolute ( std::vector< T > &  func,
const ComplexVector kern 
)
inline

Definition at line 156 of file LArFFTW.h.

157  {
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 }
int fFreqSize
Definition: LArFFTW.h:59
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::Convolute ( std::vector< T > &  func,
std::vector< T > &  resp 
)
inline

Definition at line 186 of file LArFFTW.h.

187  {
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 }
ComplexVector fKern
Definition: LArFFTW.h:55
int fFreqSize
Definition: LArFFTW.h:59
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::Correlate ( std::vector< T > &  func,
const ComplexVector kern 
)
inline

Definition at line 295 of file LArFFTW.h.

296  {
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 }
int fFreqSize
Definition: LArFFTW.h:59
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::Correlate ( std::vector< T > &  func,
std::vector< T > &  resp 
)
inline

Definition at line 326 of file LArFFTW.h.

327  {
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 }
ComplexVector fKern
Definition: LArFFTW.h:55
int fFreqSize
Definition: LArFFTW.h:59
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::Deconvolute ( std::vector< T > &  func,
const ComplexVector kern 
)
inline

Definition at line 221 of file LArFFTW.h.

222  {
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 }
int fFreqSize
Definition: LArFFTW.h:59
process_name gaushit a
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
do i e
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::Deconvolute ( std::vector< T > &  func,
std::vector< T > &  resp 
)
inline

Definition at line 255 of file LArFFTW.h.

256  {
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 }
ComplexVector fKern
Definition: LArFFTW.h:55
int fFreqSize
Definition: LArFFTW.h:59
process_name gaushit a
void * rIn
Definition: LArFFTW.h:63
void * fOut
Definition: LArFFTW.h:61
do i e
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76
template<class T >
void util::LArFFTW::DoFFT ( std::vector< T > &  input)
inline

Definition at line 76 of file LArFFTW.h.

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 }
void * fIn
Definition: LArFFTW.h:60
const void * fPlan
Definition: LArFFTW.h:62
pdgs p
Definition: selectors.fcl:22
void * fOut
Definition: LArFFTW.h:61
template<class T >
void util::LArFFTW::DoFFT ( std::vector< T > &  input,
ComplexVector output 
)
inline

Definition at line 92 of file LArFFTW.h.

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 }
void * fIn
Definition: LArFFTW.h:60
int fFreqSize
Definition: LArFFTW.h:59
const void * fPlan
Definition: LArFFTW.h:62
pdgs p
Definition: selectors.fcl:22
void * fOut
Definition: LArFFTW.h:61
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
template<class T >
void util::LArFFTW::DoInvFFT ( std::vector< T > &  output)
inline

Definition at line 113 of file LArFFTW.h.

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 }
void * rOut
Definition: LArFFTW.h:64
then echo ***************************************echo array
Definition: find_fhicl.sh:28
void * rIn
Definition: LArFFTW.h:63
const void * rPlan
Definition: LArFFTW.h:65
template<class T >
void util::LArFFTW::DoInvFFT ( ComplexVector input,
std::vector< T > &  output 
)
inline

Definition at line 131 of file LArFFTW.h.

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 }
int fFreqSize
Definition: LArFFTW.h:59
void * rOut
Definition: LArFFTW.h:64
then echo ***************************************echo array
Definition: find_fhicl.sh:28
void * rIn
Definition: LArFFTW.h:63
const void * rPlan
Definition: LArFFTW.h:65
template<class T >
T util::LArFFTW::PeakCorrelation ( std::vector< T > &  shape1,
std::vector< T > &  shape2 
)
inline

Definition at line 393 of file LArFFTW.h.

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 }
BEGIN_PROLOG TPC Trig offset(g4 rise time) ProjectToHeight
Definition: CORSIKAGen.fcl:7
double lambda(double a, double b, double c)
pdgs p
Definition: selectors.fcl:22
std::vector< float > fConvHist
Definition: LArFFTW.h:57
int fFitBins
Definition: LArFFTW.h:66
void Correlate(std::vector< T > &func, const ComplexVector &kern)
Definition: LArFFTW.h:295
gshf::MarqFitAlg * fMarqFitAlg
Definition: LArFFTW.h:68
do i e
physics associatedGroupsWithLeft p1
void util::LArFFTW::ShiftData ( ComplexVector input,
double  shift 
)

Definition at line 46 of file LArFFTW.cxx.

47 {
48  double factor = -2.0*std::acos(-1)*shift/(double)fSize;
49 
50  for(int i = 0; i < fFreqSize; i++){
51  input[i] *= std::exp(std::complex<double>(0,factor*(double)i));
52  }
53 
54  return;
55 }
int fFreqSize
Definition: LArFFTW.h:59
shift
Definition: fcl_checks.sh:26
template<class T >
void util::LArFFTW::ShiftData ( std::vector< T > &  input,
double  shift 
)
inline

Definition at line 362 of file LArFFTW.h.

363 {
364  DoFFT(input,fCompTemp);
366  DoInvFFT(fCompTemp,input);
367 
368  return;
369 }
shift
Definition: fcl_checks.sh:26
void ShiftData(ComplexVector &input, double shift)
Definition: LArFFTW.cxx:46
ComplexVector fCompTemp
Definition: LArFFTW.h:56
void DoInvFFT(std::vector< T > &output)
Definition: LArFFTW.h:113
void DoFFT(std::vector< T > &input)
Definition: LArFFTW.h:76

Member Data Documentation

ComplexVector util::LArFFTW::fCompTemp
private

Definition at line 56 of file LArFFTW.h.

std::vector<float> util::LArFFTW::fConvHist
private

Definition at line 57 of file LArFFTW.h.

int util::LArFFTW::fFitBins
private

Definition at line 66 of file LArFFTW.h.

int util::LArFFTW::fFreqSize
private

Definition at line 59 of file LArFFTW.h.

void* util::LArFFTW::fIn
private

Definition at line 60 of file LArFFTW.h.

ComplexVector util::LArFFTW::fKern
private

Definition at line 55 of file LArFFTW.h.

gshf::MarqFitAlg* util::LArFFTW::fMarqFitAlg
private

Definition at line 68 of file LArFFTW.h.

void* util::LArFFTW::fOut
private

Definition at line 61 of file LArFFTW.h.

const void* util::LArFFTW::fPlan
private

Definition at line 62 of file LArFFTW.h.

int util::LArFFTW::fSize
private

Definition at line 58 of file LArFFTW.h.

void* util::LArFFTW::rIn
private

Definition at line 63 of file LArFFTW.h.

void* util::LArFFTW::rOut
private

Definition at line 64 of file LArFFTW.h.

const void* util::LArFFTW::rPlan
private

Definition at line 65 of file LArFFTW.h.


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