LArFFT.h

Go to the documentation of this file.

#ifndef LARFFT_H
#define LARFFT_H

#include "TComplex.h"
#include "TFFTRealComplex.h"
#include "TFFTComplexReal.h"
#include "TF1.h"
#include "TH1D.h"
#include <vector>
#include <string>

#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Registry/ActivityRegistry.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Registry/ServiceMacros.h"

namespace util{
    class LArFFT {
    public:
      LArFFT(fhicl::ParameterSet const& pset, art::ActivityRegistry& reg);
      ~LArFFT();
      
      
      template <class T> void         DoFFT(std::vector<T> & input,
                                            std::vector<TComplex> & output);            
      
      template <class T> void         DoInvFFT(std::vector<TComplex> & input,
                                               std::vector<T> & output);
      
      template <class T> void        Deconvolute(std::vector<T> & input,
                                                 std::vector<T> & respFunc);
      
      template <class T> void        Deconvolute(std::vector<T> & input,
                                                 std::vector<TComplex> & kern);
      
      template <class T> void        Convolute(std::vector<T> & input,
                                               std::vector<T> & respFunc);
      
      template <class T> void        Convolute(std::vector<T> & input,
                                               std::vector<TComplex> & kern);
      
      template <class T> void        Correlate(std::vector<T> & input,
                                               std::vector<T> & respFunc);
      
      template <class T> void        Correlate(std::vector<T> & input,
                                               std::vector<TComplex> & kern);
      
      template <class T> void        AlignedSum(std::vector<T> & input,
                                                std::vector<T> &output,
                                                bool add = true);
      
                         void        ShiftData(std::vector<TComplex> & input,
                                               double shift);
  
      template <class T> void        ShiftData(std::vector<T> & input,
                                               double shift);
      
      template <class T> T           PeakCorrelation(std::vector<T> &shape1,
                                                     std::vector<T> &shape2);       
      
      int   FFTSize()          const { return fSize; }
      std::string FFTOptions() const { return fOption; }
      int FFTFitBins()         const { return fFitBins; }

      void ReinitializeFFT(int, std::string, int);

        private:
      
      int                    fSize;       //size of transform
      int                    fFreqSize;   //size of frequency space
      std::string            fOption;     //FFTW setting
      int                    fFitBins;    //Bins used for peak fit
      TF1                   *fPeakFit;    //Gaussian peak function
      TH1D                  *fConvHist;   //Fit data histogram
      std::vector<TComplex>  fCompTemp;   //temporary complex data
      std::vector<TComplex>  fKern;       //transformed response function
       
      TFFTRealComplex       *fFFT;        
      TFFTComplexReal       *fInverseFFT; 

      void InitializeFFT();
      
    }; // class LArFFT

} //namespace util

// "Forward" Fourier Transform
//--------------------------------------------------------
template <class T> inline void util::LArFFT::DoFFT(std::vector<T> & input,                 
                                                   std::vector<TComplex> & output)
{  
  double real      = 0.;  //real value holder  
  double imaginary = 0.;  //imaginary value hold   

  // set the points
  for(size_t p = 0; p < input.size(); ++p)
    fFFT->SetPoint(p, input[p]);    
  
  fFFT->Transform();    
  
  for(int i = 0; i < fFreqSize; ++i){    
    fFFT->GetPointComplex(i, real, imaginary);    
    output[i]=TComplex(real, imaginary);  
  }  
  
  return;
}

//Inverse Fourier Transform
//-------------------------------------------------
template <class T> inline void util::LArFFT::DoInvFFT(std::vector<TComplex> & input,                  
                                                      std::vector<T> & output)
{  
  for(int i = 0; i < fFreqSize; ++i)    
    fInverseFFT->SetPointComplex(i, input[i]);   

  fInverseFFT->Transform();  
  double factor = 1.0/(double) fSize;  
  
  for(int i = 0; i < fSize; ++i)    
    output[i] = factor*fInverseFFT->GetPointReal(i,false);  

  return;
}

//Deconvolution scheme taking all time-domain
//information
//--------------------------------------------------
template <class T> inline void util::LArFFT::Deconvolute(std::vector<T> & input,                         
                                                         std::vector<T> & respFunction)
{  
  DoFFT(respFunction, fKern);  
  DoFFT(input, fCompTemp);  

  for(int i = 0; i < fFreqSize; i++) 
    fCompTemp[i]/=fKern[i];  

  DoInvFFT(fCompTemp, input);  

  return;
}

//Deconvolution scheme using an already transformed
//response function
//saves cpu time if same response function is used
//for many consecutive transforms
//--------------------------------------------------
template <class T> inline void util::LArFFT::Deconvolute(std::vector<T> & input,                         
                                                         std::vector<TComplex> & kern)
{    
  DoFFT(input, fCompTemp);  

  for(int i = 0; i < fFreqSize; i++) 
    fCompTemp[i]/=kern[i];  

  DoInvFFT(fCompTemp, input);  
  
  return;
}

//Convolution scheme taking all time-domain
//information
//--------------------------------------------------
template <class T> inline void util::LArFFT::Convolute(std::vector<T> & shape1,
                                                       std::vector<T> & shape2)
{  
  DoFFT(shape1, fKern);  
  DoFFT(shape2, fCompTemp);  

  for(int i = 0; i < fFreqSize; i++) 
    fCompTemp[i]*=fKern[i];  

  DoInvFFT(fCompTemp, shape1);  

  return;
}

//Convolution scheme using an already transformed
//response function
//saves cpu time if same response function is used
//for many consecutive transforms
//--------------------------------------------------
template <class T> inline void util::LArFFT::Convolute(std::vector<T> & input,
                                                       std::vector<TComplex> & kern)
{  
  DoFFT(input, fCompTemp);  
  
  for(int i = 0; i < fFreqSize; i++) 
    fCompTemp[i]*=kern[i];  
  
  DoInvFFT(fCompTemp, input);    
  
  return;
}

//Correlation taking all time domain data
//--------------------------------------------------
template <class T> inline void util::LArFFT::Correlate(std::vector<T> & shape1,
                                                       std::vector<T> & shape2)
{  
  DoFFT(shape1, fKern);  
  DoFFT(shape2, fCompTemp);  
  
  for(int i = 0; i < fFreqSize; i++)     
    fCompTemp[i]*=TComplex::Conjugate(fKern[i]);  
  
  DoInvFFT(fCompTemp, shape1);  

  return;
}

//Convolution scheme using an already transformed
//response function
//saves cpu time if same response function is used
//for many consecutive transforms
//--------------------------------------------------
template <class T> inline void util::LArFFT::Correlate(std::vector<T> & input,  
                                                       std::vector<TComplex> & kern)
{  
  DoFFT(input, fCompTemp);  

  for(int i = 0; i < fFreqSize; i++)     
    fCompTemp[i]*=TComplex::Conjugate(kern[i]);  

  DoInvFFT(fCompTemp, input);    
  
  return;
}

//Scheme for adding two signals which have an arbitrary
//relative translation.  Shape1 is translated over shape2
//and is replaced with the sum, or the translated result
//if add = false
//--------------------------------------------------
template <class T> inline void util::LArFFT::AlignedSum(std::vector<T> & shape1,
                                                        std::vector<T> & shape2,     
                                                        bool add)
{  
  double shift = PeakCorrelation(shape1,shape2);    
  
  ShiftData(shape1,shift);  
   
  if(add)for(int i = 0; i < fSize; i++) shape1[i]+=shape2[i];   

  return;
}

//Shifts real vectors using above function
//--------------------------------------------------
template <class T> inline void util::LArFFT::ShiftData(std::vector<T> & input,  
                                                       double shift)
{ 
  DoFFT(input,fCompTemp);
  ShiftData(fCompTemp,shift);
  DoInvFFT(fCompTemp,input); 

  return;
}

//Returns the length of the translation at which the correlation
//of 2 signals is maximal.
//--------------------------------------------------
template <class T> inline T util::LArFFT::PeakCorrelation(std::vector<T> & shape1,      
                                                          std::vector<T> & shape2)
{ 
  fConvHist->Reset("ICE");
  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;     
    fConvHist->Fill(i,holder[i+startT+offset]);    
  }  
  
  fPeakFit->SetParameters(fConvHist->GetMaximum(),fFitBins/2,fFitBins/2);  
  fConvHist->Fit(fPeakFit,"QWNR","",0,fFitBins);  
  return fPeakFit->GetParameter(1)+startT;
}

DECLARE_ART_SERVICE(util::LArFFT, LEGACY)
#endif // LARFFT_H