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 >
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.

  : fSize    (transformSize)
  , fPlan    (fplan)
  , rPlan    (rplan)
  , fFitBins (fitbins)
{

  fFreqSize = fSize/2+1;

  // ... Real-Complex
  fIn = fftw_malloc(sizeof(double)*fSize);
  fOut= fftw_malloc(sizeof(fftw_complex)*fFreqSize);

  // ... Complex-Real
  rIn = fftw_malloc(sizeof(fftw_complex)*fFreqSize);
  rOut= fftw_malloc(sizeof(double)*fSize);

  // ... allocate other data vectors
  fCompTemp.resize(fFreqSize);
  fKern.resize(fFreqSize);
  fConvHist.resize(fFitBins);
}

util::LArFFTW::~LArFFTW

(

)

Definition at line 28 of file LArFFTW.cxx.

{
  fPlan = 0;
  fftw_free(fIn);
  fIn = 0;
  fftw_free((fftw_complex*)fOut);
  fOut = 0;

  rPlan = 0;
  fftw_free((fftw_complex*)rIn);
  rIn = 0;
  fftw_free(rOut);
  rOut = 0;
}

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.

{
  double shift = PeakCorrelation(shape1,shape2);

  ShiftData(shape1,shift);

  if(add)for(int i = 0; i < fSize; i++) shape1[i]+=shape2[i];

  return;
}

template<class T >

void util::LArFFTW::Convolute ( std::vector< T > &  func, const ComplexVector &  kern  )

inline

Definition at line 156 of file LArFFTW.h.

                                                                    {

  // ... Make sure that time series and kernel have the correct size.
  int n = func.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
  }
  n = kern.size();
  if(n != fFreqSize){
    throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
  }

  DoFFT(func);

  // ..perform the convolution
  for(int i = 0; i < fFreqSize; ++i){
    double re = ((fftw_complex*)fOut)[i][0];
    double im = ((fftw_complex*)fOut)[i][1];
    ((fftw_complex*)rIn)[i][0] = re*kern[i].real()-im*kern[i].imag();
    ((fftw_complex*)rIn)[i][1] = re*kern[i].imag()+im*kern[i].real();
  }

  DoInvFFT(func);
}

template<class T >

void util::LArFFTW::Convolute ( std::vector< T > &  func, std::vector< T > &  resp  )

inline

Definition at line 186 of file LArFFTW.h.

                                                              {

  // ... Make sure that time series has the correct size.
  int n = func1.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
  }
  n = func2.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
  }

  DoFFT(func2);
  for(int i = 0; i < fFreqSize; ++i){
    fKern[i].real(((fftw_complex*)fOut)[i][0]);
    fKern[i].imag(((fftw_complex*)fOut)[i][1]);
  }
  DoFFT(func1);

  // ..perform the convolution
  for(int i = 0; i < fFreqSize; ++i){
    double re = ((fftw_complex*)fOut)[i][0];
    double im = ((fftw_complex*)fOut)[i][1];
    ((fftw_complex*)rIn)[i][0] = re*fKern[i].real()-im*fKern[i].imag();
    ((fftw_complex*)rIn)[i][1] = re*fKern[i].imag()+im*fKern[i].real();
  }

  DoInvFFT(func1);
}

template<class T >

void util::LArFFTW::Correlate ( std::vector< T > &  func, const ComplexVector &  kern  )

inline

Definition at line 295 of file LArFFTW.h.

                                                                    {

  // ... Make sure that time series and kernel have the correct size.
  int n = func.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
  }
  n = kern.size();
  if(n != fFreqSize){
    throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
  }

  DoFFT(func);

  // ..perform the correlation
  for(int i = 0; i < fFreqSize; ++i){
    double re = ((fftw_complex*)fOut)[i][0];
    double im = ((fftw_complex*)fOut)[i][1];
    ((fftw_complex*)rIn)[i][0] =  re*kern[i].real()+im*kern[i].imag();
    ((fftw_complex*)rIn)[i][1] = -re*kern[i].imag()+im*kern[i].real();
  }

  DoInvFFT(func);

}

template<class T >

void util::LArFFTW::Correlate ( std::vector< T > &  func, std::vector< T > &  resp  )

inline

Definition at line 326 of file LArFFTW.h.

                                                              {

  // ... Make sure that time series has the correct size.
  int n = func1.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
  }
  n = func2.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
  }

  DoFFT(func2);
  for(int i = 0; i < fFreqSize; ++i){
    fKern[i].real(((fftw_complex*)fOut)[i][0]);
    fKern[i].imag(((fftw_complex*)fOut)[i][1]);
  }
  DoFFT(func1);

  // ..perform the correlation
  for(int i = 0; i < fFreqSize; ++i){
    double re = ((fftw_complex*)fOut)[i][0];
    double im = ((fftw_complex*)fOut)[i][1];
    ((fftw_complex*)rIn)[i][0] =  re*fKern[i].real()+im*fKern[i].imag();
    ((fftw_complex*)rIn)[i][1] = -re*fKern[i].imag()+im*fKern[i].real();
  }

  DoInvFFT(func1);

}

template<class T >

void util::LArFFTW::Deconvolute ( std::vector< T > &  func, const ComplexVector &  kern  )

inline

Definition at line 221 of file LArFFTW.h.

                                                                      {

  // ... Make sure that time series and kernel have the correct size.
  int n = func.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad time series size = " << n << "\n";
  }
  n = kern.size();
  if(n != fFreqSize){
    throw cet::exception("LArFFTW") << "Bad kernel size = " << n << "\n";
  }

  DoFFT(func);

  // ..perform the deconvolution
  double a,b,c,d,e;
  for(int i = 0; i < fFreqSize; ++i){
    a = ((fftw_complex*)fOut)[i][0];
    b = ((fftw_complex*)fOut)[i][1];
    c = kern[i].real();
    d = kern[i].imag();
    e = 1./(c*c+d*d);
    ((fftw_complex*)rIn)[i][0] = (a*c+b*d)*e;
    ((fftw_complex*)rIn)[i][1] = (b*c-a*d)*e;
  }

  DoInvFFT(func);
}

template<class T >

void util::LArFFTW::Deconvolute ( std::vector< T > &  func, std::vector< T > &  resp  )

inline

Definition at line 255 of file LArFFTW.h.

                                                               {

  // ... Make sure that time series has the correct size.
  int n = func.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 1st time series size = " << n << "\n";
  }
  n = resp.size();
  if(n != fSize){
    throw cet::exception("LArFFTW") << "Bad 2nd time series size = " << n << "\n";
  }

  DoFFT(resp);
  for(int i = 0; i < fFreqSize; ++i){
    fKern[i].real(((fftw_complex*)fOut)[i][0]);
    fKern[i].imag(((fftw_complex*)fOut)[i][1]);
  }
  DoFFT(func);

  // ..perform the deconvolution
  double a,b,c,d,e;
  for(int i = 0; i < fFreqSize; ++i){
    a = ((fftw_complex*)fOut)[i][0];
    b = ((fftw_complex*)fOut)[i][1];
    c = fKern[i].real();
    d = fKern[i].imag();
    e = 1./(c*c+d*d);
    ((fftw_complex*)rIn)[i][0] = (a*c+b*d)*e;
    ((fftw_complex*)rIn)[i][1] = (b*c-a*d)*e;
  }

  DoInvFFT(func);

}

template<class T >

void util::LArFFTW::DoFFT ( std::vector< T > &  input )

inline

Definition at line 76 of file LArFFTW.h.

{
  // ..set point
  for(size_t p = 0; p < input.size(); ++p){
    ((double *)fIn)[p] = input[p];
  }

  // ..transform (using the New-array Execute Functions)
  fftw_execute_dft_r2c((fftw_plan)fPlan,(double*)fIn,(fftw_complex*)fOut);

  return;
}

template<class T >

void util::LArFFTW::DoFFT ( std::vector< T > &  input, ComplexVector &  output  )

inline

Definition at line 92 of file LArFFTW.h.

{
  // ..set point
  for(size_t p = 0; p < input.size(); ++p){
    ((double *)fIn)[p] = input[p];
  }

  // ..transform (using the New-array Execute Functions)
  fftw_execute_dft_r2c((fftw_plan)fPlan,(double*)fIn,(fftw_complex*)fOut);

  for(int i = 0; i < fFreqSize; ++i){
    output[i].real(((fftw_complex*)fOut)[i][0]);
    output[i].imag(((fftw_complex*)fOut)[i][1]);
  }

  return;
}

template<class T >

void util::LArFFTW::DoInvFFT ( std::vector< T > &  output )

inline

Definition at line 113 of file LArFFTW.h.

{
  // ..transform (using the New-array Execute Functions)
  fftw_execute_dft_c2r((fftw_plan)rPlan,(fftw_complex*)rIn,(double*)rOut);

  // ..get point real
  double factor = 1.0/(double) fSize;
  const double * array =  (const double*)(rOut);
  for(int i = 0; i < fSize; ++i){
    output[i] = factor*array[i];
  }

  return;
}

template<class T >

void util::LArFFTW::DoInvFFT ( ComplexVector &  input, std::vector< T > &  output  )

inline

Definition at line 131 of file LArFFTW.h.

{
  // ..set point complex
  for(int i = 0; i < fFreqSize; ++i){
    ((fftw_complex*)rIn)[i][0] = input[i].real();
    ((fftw_complex*)rIn)[i][1] = input[i].imag();
  }

  // ..transform (using the New-array Execute Functions)
  fftw_execute_dft_c2r((fftw_plan)rPlan,(fftw_complex*)rIn,(double*)rOut);

  // ..get point real
  double factor = 1.0/(double) fSize;
  const double * array =  (const double*)(rOut);
  for(int i = 0; i < fSize; ++i){
    output[i] = factor*array[i];
  }

  return;
}

template<class T >

T util::LArFFTW::PeakCorrelation ( std::vector< T > &  shape1, std::vector< T > &  shape2  )

inline

Definition at line 393 of file LArFFTW.h.

{
  float chiSqr = std::numeric_limits<float>::max();
  float dchiSqr = std::numeric_limits<float>::max();
  const float chiCut   = 1e-3;
  float lambda  = 0.001;        // Marquardt damping parameter
  std::vector<float> p;

  std::vector<T> holder = shape1;
  Correlate(holder,shape2);

  int   maxT   = max_element(holder.begin(), holder.end())-holder.begin();
  float startT = maxT-fFitBins/2;
  int   offset = 0;

  for(int i = 0; i < fFitBins; i++) {
    if(startT+i < 0) offset=fSize;
    else if(startT+i > fSize) offset=-fSize;
    else offset = 0;
    if(holder[i+startT+offset]<=0.) {
      fConvHist[i]=0.;
    } else {
      fConvHist[i]=holder[i+startT+offset];
    }
  }

  p[0] = *max_element(fConvHist.begin(), fConvHist.end());
  p[1] = fFitBins/2;
  p[2] = fFitBins/2;
  float p1 = p[1];      // save initial p[1] guess

  int fitResult{-1};
  int trial=0;
  lambda=-1.;           // initialize lambda on first call
  do{
    fitResult=fMarqFitAlg->mrqdtfit(lambda, &p[0], &fConvHist[0], 3, fFitBins, chiSqr, dchiSqr);
    trial++;
    if(fitResult){
        mf::LogWarning("LArFFTW") << "Peak Correlation Fitting failed";
        break;
    }else if (trial>100){
        break;
    }
  }
  while (fabs(dchiSqr) >= chiCut);
  if (!fitResult)p1=p[1]; // if fit succeeded, use fit result

  return p1 + 0.5 + startT;
}

void util::LArFFTW::ShiftData	(	ComplexVector &	input,
		double	shift
	)

Definition at line 46 of file LArFFTW.cxx.

{
  double factor = -2.0*std::acos(-1)*shift/(double)fSize;

  for(int i = 0; i < fFreqSize; i++){
    input[i] *= std::exp(std::complex<double>(0,factor*(double)i));
  }

  return;
}

template<class T >

void util::LArFFTW::ShiftData ( std::vector< T > &  input, double  shift  )

inline

Definition at line 362 of file LArFFTW.h.

{
  DoFFT(input,fCompTemp);
  ShiftData(fCompTemp,shift);
  DoInvFFT(fCompTemp,input);

  return;
}

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:

• lardata/lardata/Utilities/LArFFTW.h
• lardata/lardata/Utilities/LArFFTW.cxx