Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 2807)
+++ trunk/MagicSoft/Mars/Changelog	(revision 2808)
@@ -11,8 +11,7 @@
 
 
-
  2004/01/14: Markus Gaug
   
-   * manalysis/MFFT.[h,cc]
+   * mtools/MFFT.[h,cc]
      - new class to perform Fast Fourier Transforms
 
Index: trunk/MagicSoft/Mars/manalysis/MFFT.cc
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MFFT.cc	(revision 2807)
+++ 	(revision )
@@ -1,752 +1,0 @@
-/* ======================================================================== *\
-!
-! *
-! * This file is part of MARS, the MAGIC Analysis and Reconstruction
-! * Software. It is distributed to you in the hope that it can be a useful
-! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
-! * It is distributed WITHOUT ANY WARRANTY.
-! *
-! * Permission to use, copy, modify and distribute this software and its
-! * documentation for any purpose is hereby granted without fee,
-! * provided that the above copyright notice appear in all copies and
-! * that both that copyright notice and this permission notice appear
-! * in supporting documentation. It is provided "as is" without express
-! * or implied warranty.
-! *
-!
-!
-!   Author(s): Markus Gaug 01/2004  <mailto:markus@ifae.es>
-!
-!   Copyright: MAGIC Software Development, 2001-2004
-!
-!
-\* ======================================================================== */
-
-//////////////////////////////////////////////////////////////////////////////
-//                                                                          //
-//  Fast Fourier Transforms                                                 //
-//                                                                          //
-//  (Most of the code is adapted from Numerical Recipies in C++, 2nd ed.,   //
-//  pp. 509-563)                                                            //
-//                                                                          //
-//  Usage:                                                                  //
-//                                                                          //
-//  1) Functions RealFunctionFFT:  (FOURIER TRANSFORM)                      //
-//     * Take as argument arrays of real numbers,                           // 
-//       in some cases the dimension of the array has to be given separately//
-//     * Return a COMPLEX array with the following meaning:                 //
-//       array[0]: The value of F(0) (has only real component)              //
-//       array[1]: The value of F(N/2) (has only real component)            //
-//       array[2i]: The real part of F(i)                                   //
-//       array[2i+1]: The imaginary part of F(i)                            //
-//     * Note that F(N-i)* = F(i), therefore only the positive frequency    //
-//       half is stored.                                                    //
-//     * The dimension MUST be an integer power of 2,                       //
-//       otherwise, the array will be shortened!!                           //
-//                                                                          //
-//  2) Functions RealFunctionIFFT:  (INVERSER FOURIER TRANSFORM)            // 
-//     * Take as argument a COMPLEX array                                   //
-//       of Fourier-transformed REAL numbers                                //
-//       with the following meaning:                                        //
-//       array[0]: The value of F(0) (has only real component)              //
-//       array[1]: The value of F(N/2) (has only real component)            //
-//       array[2i]: The real part of F(i)                                   //
-//       array[2i+1]: The imaginary part of F(i)                            //
-//     * Returns the original complex array of dimension 2N-1               //                                
-//                                                                          //
-//  3) Functions PowerSpectrumDensity:                                      //
-//     * Return a histogram with the spectral density, i.e.                 //
-//       P(k) = 1/(N*N) * |F(k)|*|F(k)|                                     //
-//     * The histogram is ranged between 0 and 1./(2*binwidth)              //
-//     * The number of bins equals N/2+1                                    //
-//     * Note that histograms with unequal binwidth can not yet be treated! //
-//     * If the PSD does NOT CONVERGE to 0 at the maximum bin,              //
-//       you HAVE TO sample your data finer!                                //
-//                                                                          //
-//////////////////////////////////////////////////////////////////////////////
-
-#include "MFFT.h"
-
-#include "TMath.h"
-
-#include "MLog.h"
-#include "MLogManip.h"
-
-ClassImp(MFFT);
-
-using namespace std;
-
-// ---------------------------------------------------------------------------
-//
-//  Default Constructor
-//  Initializes random number generator and default variables
-//
-MFFT::MFFT() : fDim(0)
-{
-}
-
-// --------------------------------------------------------------------------
-//
-//  Destructor. 
-//
-MFFT::~MFFT()
-{
-}
-
-void MFFT::TransformF(const Int_t isign)
-{
-  
-  UInt_t   n,mmax,m,j,istep,i;
-  Float_t wtemp,wr,wpr,wpi,wi,theta;
-  Float_t tempr,tempi;
-  
-  Int_t nn = fDim/2;
-  n = nn << 1;
-
-  //
-  // The bit-reversal section of the routine:
-  // Exchange the two complex numbers
-  //
-  j=1;
-  for (i=1;i<n;i+=2) {
-    if (j > i) {
-      Swap(fDataF[j-1],fDataF[i-1]);
-      Swap(fDataF[j],fDataF[i]);
-    }
-    m=nn;
-    while (m >= 2 && j > m) {
-      j -= m;
-      m >>= 1;
-    }
-    j += m;
-  }
-  // 
-  // Here begins the Danielson-Lanczos section of the routine
-  //
-  mmax=2;
-  while (n > mmax) {         // Outer loop executed log_2(nn) times
-
-    istep = mmax << 1;
-    //
-    // Initialize the trigonometric recurrence:
-    //
-    theta = isign*(6.28318530717959/mmax);
-
-    wtemp = TMath::Sin(0.5*theta);
-    wpr   = -2.0*wtemp*wtemp;
-    wpi   = TMath::Sin(theta);
-
-    wr=1.0;
-    wi=0.0;
-
-    for (m=1; m<mmax; m+=2) {
-      for (i=m; i<=n; i+=istep) {
-        // 
-        // The Danielson-Lanczos formula:
-        //
-        j          = i+mmax;
-        tempr      = wr*fDataF[j-1] - wi*fDataF[j];
-        tempi      = wr*fDataF[j]   + wi*fDataF[j-1];
-        fDataF[j-1] = fDataF[i-1]   - tempr;
-        fDataF[j]   = fDataF[i]     - tempi;
-        fDataF[i-1] += tempr;
-        fDataF[i]   += tempi;
-      }
-
-      //
-      // Trigonometric recurrence
-      //
-      wr = (wtemp=wr)*wpr - wi*wpi+wr;
-      wi = wi*wpr         + wtemp*wpi+wi;
-
-    }
-    mmax=istep;
-  }
-}
-
-
-void MFFT::TransformD(const Int_t isign)
-{
-  
-  UInt_t   n,mmax,m,j,istep,i;
-  Double_t wtemp,wr,wpr,wpi,wi,theta;
-  Double_t tempr,tempi;
-  
-  Int_t nn = fDim/2;
-  n = nn << 1;
-
-  //
-  // The bit-reversal section of the routine:
-  // Exchange the two complex numbers
-  //
-  j=1;
-  for (i=1;i<n;i+=2) {
-    if (j > i) {
-      Swap(fDataD[j-1],fDataD[i-1]);
-      Swap(fDataD[j],fDataD[i]);
-    }
-    m=nn;
-    while (m >= 2 && j > m) {
-      j -= m;
-      m >>= 1;
-    }
-    j += m;
-  }
-  // 
-  // Here begins the Danielson-Lanczos section of the routine
-  //
-  mmax=2;
-  while (n > mmax) {         // Outer loop executed log_2(nn) times
-
-    istep = mmax << 1;
-    //
-    // Initialize the trigonometric recurrence:
-    //
-    theta = isign*(6.28318530717959/mmax);
-
-    wtemp = TMath::Sin(0.5*theta);
-    wpr   = -2.0*wtemp*wtemp;
-    wpi   = TMath::Sin(theta);
-
-    wr=1.0;
-    wi=0.0;
-
-    for (m=1; m<mmax; m+=2) {
-      for (i=m; i<=n; i+=istep) {
-        // 
-        // The Danielson-Lanczos formula:
-        //
-        j          = i+mmax;
-        tempr      = wr*fDataD[j-1] - wi*fDataD[j];
-        tempi      = wr*fDataD[j]   + wi*fDataD[j-1];
-        fDataD[j-1] = fDataD[i-1]   - tempr;
-        fDataD[j]   = fDataD[i]     - tempi;
-        fDataD[i-1] += tempr;
-        fDataD[i]   += tempi;
-      }
-
-      //
-      // Trigonometric recurrence
-      //
-      wr = (wtemp=wr)*wpr - wi*wpi+wr;
-      wi = wi*wpr         + wtemp*wpi+wi;
-
-    }
-    mmax=istep;
-  }
-}
-
-//
-// Calculates the Fourier transform of a set of n real-valued data points. 
-// Replaces this data (which is stored in array data[1..n]) by the positive
-// frequency half of its complex Fourier transform. The real-valued first 
-// and last components of the complex transform are returned as elements 
-// data[1] and data[2], respectively. n must be a power of 2. This routine
-// also calculates the inverse transform of a complex data array if it is 
-// the transform of real data. (Result in this case mus be multiplied by 
-// 2/n.). From NUMERICAL RECIPES IN C++.
-//
-void MFFT::RealFTF(const Int_t isign)
-{
-  
-  Int_t    i,i1,i2,i3,i4;
-  Float_t  c1=0.5,c2,h1r,h1i,h2r,h2i;
-  Float_t wr,wi,wpr,wpi,wtemp,theta;
-
-  //
-  // Initialize the recurrence
-  //
-  theta=3.141592653589793/(Double_t) (fDim>>1);
-
-  if(isign==1) // forward transform
-    {
-      c2    = -0.5;
-      TransformF(1);
-    }
-  else         // set up backward transform
-    {
-      c2    = 0.5;
-      theta = -theta;
-    }
-
-  wtemp = TMath::Sin(0.5*theta);
-  wpr   = -2.0*wtemp*wtemp;
-  wpi   = TMath::Sin(theta);
-
-  wr    = 1.0 + wpr;
-  wi    = wpi;
-
-  for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
-    {
-
-      i2 = 1 + (i1 = i+i);
-      i4 = 1 + (i3 = fDim-i1);
-
-      //
-      // The two separate transforms are separated out of the data
-      //
-      h1r  =  c1*(fDataF[i1]+fDataF[i3]);
-      h1i  =  c1*(fDataF[i2]-fDataF[i4]);
-      h2r  = -c2*(fDataF[i2]+fDataF[i4]);
-      h2i  =  c2*(fDataF[i1]-fDataF[i3]);
-
-      //
-      // Here, they are recombined to from the true transform 
-      // of the orginal real data
-      //
-      fDataF[i1] =  h1r + wr*h2r - wi*h2i;
-      fDataF[i2] =  h1i + wr*h2i + wi*h2r;
-      fDataF[i3] =  h1r - wr*h2r + wi*h2i;
-      fDataF[i4] = -h1i + wr*h2i + wi*h2r;
-      
-      //
-      // The recurrence
-      //
-      wr = (wtemp=wr)*wpr - wi*wpi + wr;
-      wi =    wi*wpr   + wtemp*wpi + wi;
-    }
-
-  //
-  // Squeeze the first and last data together to get them all 
-  // within the original array
-  //
-  if(isign==1)
-    {
-      fDataF[0] = (h1r=fDataF[0]) + fDataF[1];
-      fDataF[1] =     h1r  -      fDataF[1];
-    }
-  else
-    {
-
-      fDataF[0] = c1*((h1r=fDataF[0]) + fDataF[1]);
-      fDataF[1] = c1*(h1r-fDataF[1]);
-
-      //
-      // The inverse transform for the case isign = -1
-      //
-      TransformF(-1);  
-
-      //
-      // normalize correctly (not done in original NR's)
-      //
-      for(i=1;i<=fDim;i++)
-        fDataF[i] *= (2./fDim);
-    }
-}
-void MFFT::RealFTD(const Int_t isign)
-{
-  
-  Int_t    i,i1,i2,i3,i4;
-  Float_t  c1=0.5,c2,h1r,h1i,h2r,h2i;
-  Double_t wr,wi,wpr,wpi,wtemp,theta;
-
-  //
-  // Initialize the recurrence
-  //
-  theta=3.141592653589793/(Double_t) (fDim>>1);
-
-  if(isign==1) // forward transform
-    {
-      c2    = -0.5;
-      TransformD(1);
-    }
-  else         // set up backward transform
-    {
-      c2    = 0.5;
-      theta = -theta;
-    }
-
-  wtemp = TMath::Sin(0.5*theta);
-  wpr   = -2.0*wtemp*wtemp;
-  wpi   = TMath::Sin(theta);
-
-  wr    = 1.0 + wpr;
-  wi    = wpi;
-
-  for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
-    {
-
-      i2 = 1 + (i1 = i+i);
-      i4 = 1 + (i3 = fDim-i1);
-
-      //
-      // The two separate transforms are separated out of the data
-      //
-      h1r  =  c1*(fDataD[i1]+fDataD[i3]);
-      h1i  =  c1*(fDataD[i2]-fDataD[i4]);
-      h2r  = -c2*(fDataD[i2]+fDataD[i4]);
-      h2i  =  c2*(fDataD[i1]-fDataD[i3]);
-
-      //
-      // Here, they are recombined to from the true transform 
-      // of the orginal real data
-      //
-      fDataD[i1] =  h1r + wr*h2r - wi*h2i;
-      fDataD[i2] =  h1i + wr*h2i + wi*h2r;
-      fDataD[i3] =  h1r - wr*h2r + wi*h2i;
-      fDataD[i4] = -h1i + wr*h2i + wi*h2r;
-      
-      //
-      // The recurrence
-      //
-      wr = (wtemp=wr)*wpr - wi*wpi + wr;
-      wi =    wi*wpr   + wtemp*wpi + wi;
-    }
-
-  //
-  // Squeeze the first and last data together to get them all 
-  // within the original array
-  //
-  if(isign==1)
-    {
-      fDataD[0] = (h1r=fDataD[0]) + fDataD[1];
-      fDataD[1] =     h1r  -      fDataD[1];
-    }
-  else
-    {
-
-      fDataD[0] = c1*((h1r=fDataD[0]) + fDataD[1]);
-      fDataD[1] = c1*(h1r-fDataD[1]);
-
-      //
-      // The inverse transform for the case isign = -1
-      //
-      TransformD(-1);  
-      
-      //
-      // normalize correctly (not done in original NR's)
-      //
-      for(i=1;i<=fDim;i++)
-        fDataD[i] *= (2./fDim);
-    }
-}
-
-
-//
-// Fast Fourier Transform for float arrays
-//
-Float_t* MFFT::RealFunctionFFT(const Int_t n, const Float_t *data)
-{
-
-  fDim = n;
-  CheckDim(n);
-
-  fDataF.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataF[i] = data[i];
-
-  RealFTF(1);
-  
-  return fDataF.GetArray();
-
-}
-
-//
-// Fast Inverse Fourier Transform for float arrays
-//
-Float_t* MFFT::RealFunctionIFFT(const Int_t n, const Float_t *data)
-{
-
-  fDim = n;
-  CheckDim(fDim);
-  
-  fDataF.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataF[i] = data[i];
-
-  RealFTF(-1);
-  
-  return fDataF.GetArray();
-
-}
-
-//
-// Fast Fourier Transform for double arrays
-//
-Double_t* MFFT::RealFunctionFFT(const Int_t n, const Double_t *data)
-{
-
-  fDim = n;
-  CheckDim(n);
-
-  fDataD.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataD[i] = data[i];
-
-  RealFTD(1);
-  
-  return fDataD.GetArray();
-
-}
-
-//
-// Fast Inverse Fourier Transform for double arrays
-//
-Double_t* MFFT::RealFunctionIFFT(const Int_t n, const Double_t *data)
-{
-
-  fDim = n;
-  CheckDim(fDim);
-  
-  fDataD.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataD[i] = data[i];
-
-  RealFTD(-1);
-  
-  return fDataD.GetArray();
-
-}
-
-//
-// Fast Fourier Transform for TArrayF's
-//
-TArrayF* MFFT::RealFunctionFFT(const TArrayF *data)
-{
-
-  fDim = data->GetSize();
-  CheckDim(fDim);
-
-  fDataF.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataF[i] = data->At(i);
-
-  RealFTF(1);
-  
-  return new TArrayF(fDim,fDataF.GetArray());
-
-}
-
-//
-// Inverse Fast Fourier Transform for TArrayF's
-//
-TArrayF* MFFT::RealFunctionIFFT(const TArrayF *data)
-{
-
-  fDim = data->GetSize();
-  CheckDim(fDim);
-  
-  fDataF.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataF[i] = data->At(i);
-
-  RealFTF(-1);
-
-  return new TArrayF(fDim,fDataF.GetArray());
-}
-
-
-//
-// Fast Fourier Transform for TArrayD's
-//
-TArrayD* MFFT::RealFunctionFFT(const TArrayD *data)
-{
-
-  fDim = data->GetSize();
-  CheckDim(fDim);
-
-  fDataD.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataD[i] = data->At(i);
-
-  RealFTD(1);
-  
-  return new TArrayD(fDim,fDataD.GetArray());
-
-}
-
-//
-// Inverse Fast Fourier Transform for TArrayD's
-//
-TArrayD* MFFT::RealFunctionIFFT(const TArrayD *data)
-{
-
-  fDim = data->GetSize();
-  CheckDim(fDim);
-  
-  fDataD.Set(fDim);
-  //
-  // Clone the array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataD[i] = data->At(i);
-
-  RealFTD(-1);
-
-  return new TArrayD(fDim,fDataD.GetArray());
-}
-
-
-//
-// Power Spectrum Density Calculation
-//
-TH1D* MFFT::PowerSpectrumDensity(const TH1D *hist)
-{
-
-  TH1D *newhist = (TH1D*)CheckHist(hist,1);
-
-  fDataD.Set(fDim);
-  //
-  // Copy the hist into an array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataD[i] = hist->GetBinContent(i);
-
-  RealFTD(1);
-
-  Int_t dim2 = fDim*fDim;
-  Double_t c02;
-  Double_t ck2;
-  Double_t cn2;
-  //
-  // Fill the new histogram: 
-  //
-  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
-  //    (stored in fData{0])
-  //
-  c02 = fDataD[0]*fDataD[0];
-  newhist->Fill(c02/dim2);
-  //
-  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)| + |C(N-k)|*|C(N-k)|)
-  //
-  for (Int_t k=2;k<fDim-2;k+=2)
-    {
-
-      Int_t ki  = k+1;
-      ck2 = (fDataD[k]*fDataD[k] + fDataD[ki]*fDataD[ki]);
-      newhist->Fill(ck2/dim2);
-    }
-  //
-  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
-  //    (stored in fData[1])
-  //
-  cn2 = (fDataD[1]*fDataD[1]);
-  newhist->Fill(cn2/dim2);
-
-  return newhist;
-}
-
-
-//
-// Power Spectrum Density calculation
-//
-TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
-{
-
-  TH1F *newhist = (TH1F*)CheckHist(hist,0);
-
-  fDataF.Set(fDim);
-  //
-  // Copy the hist into an array
-  //
-  for (Int_t i=0;i<fDim;i++)
-    fDataF[i] = hist->GetBinContent(i);
-
-  RealFTF(1);
-
-  Int_t dim2 = fDim*fDim;
-  Float_t c02;
-  Float_t ck2;
-  Float_t cn2;
-  //
-  // Fill the new histogram: 
-  //
-  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
-  //
-  c02 = (fDataF[0]*fDataF[0]);
-  newhist->Fill(c02/dim2);
-  //
-  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
-  //
-  for (Int_t k=2;k<fDim;k+=2)
-    {
-      Int_t ki  = k+1;
-      ck2 = (fDataF[k]*fDataF[k] + fDataF[ki]*fDataF[ki]);
-      newhist->Fill(ck2/dim2);
-    }
-  //
-  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
-  //
-  cn2 = (fDataF[1]*fDataF[1]);
-  newhist->Fill(cn2/dim2);
-  
-  return newhist;
-}
-
-
-void MFFT::CheckDim(Int_t a)
-{
-
-  // If even number, return 0
-  if (a==2)  return;
-
-  // If odd number, return the closest power of 2
-  if (a & 1) 
-    { 
-      Int_t b = 1; 
-      while (b < fDim/2+1)
-        b <<= 1; 
-
-      fDim = b;
-      *fLog << warn << "Dimension of Data is not a multiple of 2, will take only first " 
-            << fDim << " entries! " << endl;
-      return; 
-    }
-
-  CheckDim(a>>1);
-}
-
-TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag)
-{
-  
-  TString name = hist->GetName();
-  name += " PSD";
-  TString title = hist->GetTitle();
-  title += " - Power Spectrum Density";
-
-  // number of entries
-  fDim = hist->GetNbinsX();
-  CheckDim(fDim);
-
-  // Step width
-  Double_t delta = hist->GetBinWidth(1);
-  
-  // Nyquist frequency
-  Axis_t fcrit = 1./(2.*delta);
-  Axis_t low = 0.;
-  Axis_t up  = fcrit;
-
-  switch (flag)
-    {
-    case 0: 
-      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
-      break;
-    case 1:
-      return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up);
-      break;
-    default:
-      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
-      break;
-    }
-}
Index: trunk/MagicSoft/Mars/manalysis/MFFT.h
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MFFT.h	(revision 2807)
+++ 	(revision )
@@ -1,66 +1,0 @@
-#ifndef MARS_MFFT
-#define MARS_MFFT
-
-#ifndef MARS_MParContainer
-#include "MParContainer.h"
-#endif
-
-#ifndef ROOT_TArrayF
-#include "TArrayF.h"
-#endif
-
-#ifndef ROOT_TArrayD
-#include "TArrayD.h"
-#endif
-
-#ifndef ROOT_TH1F
-#include "TH1F.h"
-#endif
-
-#ifndef ROOT_TH1D
-#include "TH1D.h"
-#endif
-
-class MFFT : public MParContainer
-{
-private:
-
-  void Swap(Float_t &a,  Float_t &b)    { Float_t  c = a;  a = b;  b = c;  }
-  void Swap(Double_t &a, Double_t &b)   { Double_t c = a;  a = b;  b = c;  }
-
-  void TransformF(const Int_t isign);
-  void TransformD(const Int_t isign);  
-  void RealFTF(const Int_t isign);
-  void RealFTD(const Int_t isign);
-
-  void CheckDim(Int_t a);
-  TH1 *CheckHist(const TH1 *hist, const Int_t flag);
-  
-  Int_t   fDim;
-  TArrayF fDataF;
-  TArrayD fDataD;  
-
-public:
-
-  MFFT();
-  ~MFFT();
-
-  TArrayF*  RealFunctionFFT( const TArrayF *data);
-  TArrayF*  RealFunctionIFFT(const TArrayF *data);  
-  
-  TArrayD*  RealFunctionFFT( const TArrayD *data);
-  TArrayD*  RealFunctionIFFT(const TArrayD *data);  
-  
-  Float_t*  RealFunctionFFT( const Int_t n, const Float_t *data);
-  Float_t*  RealFunctionIFFT(const Int_t n, const Float_t *data);  
-  
-  Double_t* RealFunctionFFT( const Int_t n, const Double_t *data);
-  Double_t* RealFunctionIFFT(const Int_t n, const Double_t *data);  
-  
-  TH1F* PowerSpectrumDensity(const TH1F *hist);
-  TH1D* PowerSpectrumDensity(const TH1D *hist);
-  
-  ClassDef(MFFT,0)  // Class to perform a Fast Fourier Transform
-};
-    
-#endif
Index: trunk/MagicSoft/Mars/mtools/MFFT.cc
===================================================================
--- trunk/MagicSoft/Mars/mtools/MFFT.cc	(revision 2808)
+++ trunk/MagicSoft/Mars/mtools/MFFT.cc	(revision 2808)
@@ -0,0 +1,752 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Markus Gaug 01/2004  <mailto:markus@ifae.es>
+!
+!   Copyright: MAGIC Software Development, 2001-2004
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//                                                                          //
+//  Fast Fourier Transforms                                                 //
+//                                                                          //
+//  (Most of the code is adapted from Numerical Recipies in C++, 2nd ed.,   //
+//  pp. 509-563)                                                            //
+//                                                                          //
+//  Usage:                                                                  //
+//                                                                          //
+//  1) Functions RealFunctionFFT:  (FOURIER TRANSFORM)                      //
+//     * Take as argument arrays of real numbers,                           // 
+//       in some cases the dimension of the array has to be given separately//
+//     * Return a COMPLEX array with the following meaning:                 //
+//       array[0]: The value of F(0) (has only real component)              //
+//       array[1]: The value of F(N/2) (has only real component)            //
+//       array[2i]: The real part of F(i)                                   //
+//       array[2i+1]: The imaginary part of F(i)                            //
+//     * Note that F(N-i)* = F(i), therefore only the positive frequency    //
+//       half is stored.                                                    //
+//     * The dimension MUST be an integer power of 2,                       //
+//       otherwise, the array will be shortened!!                           //
+//                                                                          //
+//  2) Functions RealFunctionIFFT:  (INVERSER FOURIER TRANSFORM)            // 
+//     * Take as argument a COMPLEX array                                   //
+//       of Fourier-transformed REAL numbers                                //
+//       with the following meaning:                                        //
+//       array[0]: The value of F(0) (has only real component)              //
+//       array[1]: The value of F(N/2) (has only real component)            //
+//       array[2i]: The real part of F(i)                                   //
+//       array[2i+1]: The imaginary part of F(i)                            //
+//     * Returns the original complex array of dimension 2N-1               //                                
+//                                                                          //
+//  3) Functions PowerSpectrumDensity:                                      //
+//     * Return a histogram with the spectral density, i.e.                 //
+//       P(k) = 1/(N*N) * |F(k)|*|F(k)|                                     //
+//     * The histogram is ranged between 0 and 1./(2*binwidth)              //
+//     * The number of bins equals N/2+1                                    //
+//     * Note that histograms with unequal binwidth can not yet be treated! //
+//     * If the PSD does NOT CONVERGE to 0 at the maximum bin,              //
+//       you HAVE TO sample your data finer!                                //
+//                                                                          //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "MFFT.h"
+
+#include "TMath.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MFFT);
+
+using namespace std;
+
+// ---------------------------------------------------------------------------
+//
+//  Default Constructor
+//  Initializes random number generator and default variables
+//
+MFFT::MFFT() : fDim(0)
+{
+}
+
+// --------------------------------------------------------------------------
+//
+//  Destructor. 
+//
+MFFT::~MFFT()
+{
+}
+
+void MFFT::TransformF(const Int_t isign)
+{
+  
+  UInt_t   n,mmax,m,j,istep,i;
+  Float_t wtemp,wr,wpr,wpi,wi,theta;
+  Float_t tempr,tempi;
+  
+  Int_t nn = fDim/2;
+  n = nn << 1;
+
+  //
+  // The bit-reversal section of the routine:
+  // Exchange the two complex numbers
+  //
+  j=1;
+  for (i=1;i<n;i+=2) {
+    if (j > i) {
+      Swap(fDataF[j-1],fDataF[i-1]);
+      Swap(fDataF[j],fDataF[i]);
+    }
+    m=nn;
+    while (m >= 2 && j > m) {
+      j -= m;
+      m >>= 1;
+    }
+    j += m;
+  }
+  // 
+  // Here begins the Danielson-Lanczos section of the routine
+  //
+  mmax=2;
+  while (n > mmax) {         // Outer loop executed log_2(nn) times
+
+    istep = mmax << 1;
+    //
+    // Initialize the trigonometric recurrence:
+    //
+    theta = isign*(6.28318530717959/mmax);
+
+    wtemp = TMath::Sin(0.5*theta);
+    wpr   = -2.0*wtemp*wtemp;
+    wpi   = TMath::Sin(theta);
+
+    wr=1.0;
+    wi=0.0;
+
+    for (m=1; m<mmax; m+=2) {
+      for (i=m; i<=n; i+=istep) {
+        // 
+        // The Danielson-Lanczos formula:
+        //
+        j          = i+mmax;
+        tempr      = wr*fDataF[j-1] - wi*fDataF[j];
+        tempi      = wr*fDataF[j]   + wi*fDataF[j-1];
+        fDataF[j-1] = fDataF[i-1]   - tempr;
+        fDataF[j]   = fDataF[i]     - tempi;
+        fDataF[i-1] += tempr;
+        fDataF[i]   += tempi;
+      }
+
+      //
+      // Trigonometric recurrence
+      //
+      wr = (wtemp=wr)*wpr - wi*wpi+wr;
+      wi = wi*wpr         + wtemp*wpi+wi;
+
+    }
+    mmax=istep;
+  }
+}
+
+
+void MFFT::TransformD(const Int_t isign)
+{
+  
+  UInt_t   n,mmax,m,j,istep,i;
+  Double_t wtemp,wr,wpr,wpi,wi,theta;
+  Double_t tempr,tempi;
+  
+  Int_t nn = fDim/2;
+  n = nn << 1;
+
+  //
+  // The bit-reversal section of the routine:
+  // Exchange the two complex numbers
+  //
+  j=1;
+  for (i=1;i<n;i+=2) {
+    if (j > i) {
+      Swap(fDataD[j-1],fDataD[i-1]);
+      Swap(fDataD[j],fDataD[i]);
+    }
+    m=nn;
+    while (m >= 2 && j > m) {
+      j -= m;
+      m >>= 1;
+    }
+    j += m;
+  }
+  // 
+  // Here begins the Danielson-Lanczos section of the routine
+  //
+  mmax=2;
+  while (n > mmax) {         // Outer loop executed log_2(nn) times
+
+    istep = mmax << 1;
+    //
+    // Initialize the trigonometric recurrence:
+    //
+    theta = isign*(6.28318530717959/mmax);
+
+    wtemp = TMath::Sin(0.5*theta);
+    wpr   = -2.0*wtemp*wtemp;
+    wpi   = TMath::Sin(theta);
+
+    wr=1.0;
+    wi=0.0;
+
+    for (m=1; m<mmax; m+=2) {
+      for (i=m; i<=n; i+=istep) {
+        // 
+        // The Danielson-Lanczos formula:
+        //
+        j          = i+mmax;
+        tempr      = wr*fDataD[j-1] - wi*fDataD[j];
+        tempi      = wr*fDataD[j]   + wi*fDataD[j-1];
+        fDataD[j-1] = fDataD[i-1]   - tempr;
+        fDataD[j]   = fDataD[i]     - tempi;
+        fDataD[i-1] += tempr;
+        fDataD[i]   += tempi;
+      }
+
+      //
+      // Trigonometric recurrence
+      //
+      wr = (wtemp=wr)*wpr - wi*wpi+wr;
+      wi = wi*wpr         + wtemp*wpi+wi;
+
+    }
+    mmax=istep;
+  }
+}
+
+//
+// Calculates the Fourier transform of a set of n real-valued data points. 
+// Replaces this data (which is stored in array data[1..n]) by the positive
+// frequency half of its complex Fourier transform. The real-valued first 
+// and last components of the complex transform are returned as elements 
+// data[1] and data[2], respectively. n must be a power of 2. This routine
+// also calculates the inverse transform of a complex data array if it is 
+// the transform of real data. (Result in this case mus be multiplied by 
+// 2/n.). From NUMERICAL RECIPES IN C++.
+//
+void MFFT::RealFTF(const Int_t isign)
+{
+  
+  Int_t    i,i1,i2,i3,i4;
+  Float_t  c1=0.5,c2,h1r,h1i,h2r,h2i;
+  Float_t wr,wi,wpr,wpi,wtemp,theta;
+
+  //
+  // Initialize the recurrence
+  //
+  theta = TMath::Pi() / (Double_t)(fDim>>1);
+
+  if(isign==1) // forward transform
+    {
+      c2    = -0.5;
+      TransformF(1);
+    }
+  else         // set up backward transform
+    {
+      c2    = 0.5;
+      theta = -theta;
+    }
+
+  wtemp = TMath::Sin(0.5*theta);
+  wpr   = -2.0*wtemp*wtemp;
+  wpi   = TMath::Sin(theta);
+
+  wr    = 1.0 + wpr;
+  wi    = wpi;
+
+  for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
+    {
+
+      i2 = 1 + (i1 = i+i);
+      i4 = 1 + (i3 = fDim-i1);
+
+      //
+      // The two separate transforms are separated out of the data
+      //
+      h1r  =  c1*(fDataF[i1]+fDataF[i3]);
+      h1i  =  c1*(fDataF[i2]-fDataF[i4]);
+      h2r  = -c2*(fDataF[i2]+fDataF[i4]);
+      h2i  =  c2*(fDataF[i1]-fDataF[i3]);
+
+      //
+      // Here, they are recombined to from the true transform 
+      // of the orginal real data
+      //
+      fDataF[i1] =  h1r + wr*h2r - wi*h2i;
+      fDataF[i2] =  h1i + wr*h2i + wi*h2r;
+      fDataF[i3] =  h1r - wr*h2r + wi*h2i;
+      fDataF[i4] = -h1i + wr*h2i + wi*h2r;
+      
+      //
+      // The recurrence
+      //
+      wr = (wtemp=wr)*wpr - wi*wpi + wr;
+      wi =    wi*wpr   + wtemp*wpi + wi;
+    }
+
+  //
+  // Squeeze the first and last data together to get them all 
+  // within the original array
+  //
+  if(isign==1)
+    {
+      fDataF[0] = (h1r=fDataF[0]) + fDataF[1];
+      fDataF[1] =     h1r  -      fDataF[1];
+    }
+  else
+    {
+
+      fDataF[0] = c1*((h1r=fDataF[0]) + fDataF[1]);
+      fDataF[1] = c1*(h1r-fDataF[1]);
+
+      //
+      // The inverse transform for the case isign = -1
+      //
+      TransformF(-1);  
+
+      //
+      // normalize correctly (not done in original NR's)
+      //
+      for(i=1;i<=fDim;i++)
+        fDataF[i] *= (2./fDim);
+    }
+}
+void MFFT::RealFTD(const Int_t isign)
+{
+  
+  Int_t    i,i1,i2,i3,i4;
+  Float_t  c1=0.5,c2,h1r,h1i,h2r,h2i;
+  Double_t wr,wi,wpr,wpi,wtemp,theta;
+
+  //
+  // Initialize the recurrence
+  //
+  theta=3.141592653589793/(Double_t) (fDim>>1);
+
+  if(isign==1) // forward transform
+    {
+      c2    = -0.5;
+      TransformD(1);
+    }
+  else         // set up backward transform
+    {
+      c2    = 0.5;
+      theta = -theta;
+    }
+
+  wtemp = TMath::Sin(0.5*theta);
+  wpr   = -2.0*wtemp*wtemp;
+  wpi   = TMath::Sin(theta);
+
+  wr    = 1.0 + wpr;
+  wi    = wpi;
+
+  for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
+    {
+
+      i2 = 1 + (i1 = i+i);
+      i4 = 1 + (i3 = fDim-i1);
+
+      //
+      // The two separate transforms are separated out of the data
+      //
+      h1r  =  c1*(fDataD[i1]+fDataD[i3]);
+      h1i  =  c1*(fDataD[i2]-fDataD[i4]);
+      h2r  = -c2*(fDataD[i2]+fDataD[i4]);
+      h2i  =  c2*(fDataD[i1]-fDataD[i3]);
+
+      //
+      // Here, they are recombined to from the true transform 
+      // of the orginal real data
+      //
+      fDataD[i1] =  h1r + wr*h2r - wi*h2i;
+      fDataD[i2] =  h1i + wr*h2i + wi*h2r;
+      fDataD[i3] =  h1r - wr*h2r + wi*h2i;
+      fDataD[i4] = -h1i + wr*h2i + wi*h2r;
+      
+      //
+      // The recurrence
+      //
+      wr = (wtemp=wr)*wpr - wi*wpi + wr;
+      wi =    wi*wpr   + wtemp*wpi + wi;
+    }
+
+  //
+  // Squeeze the first and last data together to get them all 
+  // within the original array
+  //
+  if(isign==1)
+    {
+      fDataD[0] = (h1r=fDataD[0]) + fDataD[1];
+      fDataD[1] =     h1r  -      fDataD[1];
+    }
+  else
+    {
+
+      fDataD[0] = c1*((h1r=fDataD[0]) + fDataD[1]);
+      fDataD[1] = c1*(h1r-fDataD[1]);
+
+      //
+      // The inverse transform for the case isign = -1
+      //
+      TransformD(-1);  
+      
+      //
+      // normalize correctly (not done in original NR's)
+      //
+      for(i=1;i<=fDim;i++)
+        fDataD[i] *= (2./fDim);
+    }
+}
+
+
+//
+// Fast Fourier Transform for float arrays
+//
+Float_t* MFFT::RealFunctionFFT(const Int_t n, const Float_t *data)
+{
+
+  fDim = n;
+  CheckDim(n);
+
+  fDataF.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = data[i];
+
+  RealFTF(1);
+  
+  return fDataF.GetArray();
+
+}
+
+//
+// Fast Inverse Fourier Transform for float arrays
+//
+Float_t* MFFT::RealFunctionIFFT(const Int_t n, const Float_t *data)
+{
+
+  fDim = n;
+  CheckDim(fDim);
+  
+  fDataF.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = data[i];
+
+  RealFTF(-1);
+  
+  return fDataF.GetArray();
+
+}
+
+//
+// Fast Fourier Transform for double arrays
+//
+Double_t* MFFT::RealFunctionFFT(const Int_t n, const Double_t *data)
+{
+
+  fDim = n;
+  CheckDim(n);
+
+  fDataD.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = data[i];
+
+  RealFTD(1);
+  
+  return fDataD.GetArray();
+
+}
+
+//
+// Fast Inverse Fourier Transform for double arrays
+//
+Double_t* MFFT::RealFunctionIFFT(const Int_t n, const Double_t *data)
+{
+
+  fDim = n;
+  CheckDim(fDim);
+  
+  fDataD.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = data[i];
+
+  RealFTD(-1);
+  
+  return fDataD.GetArray();
+
+}
+
+//
+// Fast Fourier Transform for TArrayF's
+//
+TArrayF* MFFT::RealFunctionFFT(const TArrayF *data)
+{
+
+  fDim = data->GetSize();
+  CheckDim(fDim);
+
+  fDataF.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = data->At(i);
+
+  RealFTF(1);
+  
+  return new TArrayF(fDim,fDataF.GetArray());
+
+}
+
+//
+// Inverse Fast Fourier Transform for TArrayF's
+//
+TArrayF* MFFT::RealFunctionIFFT(const TArrayF *data)
+{
+
+  fDim = data->GetSize();
+  CheckDim(fDim);
+  
+  fDataF.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = data->At(i);
+
+  RealFTF(-1);
+
+  return new TArrayF(fDim,fDataF.GetArray());
+}
+
+
+//
+// Fast Fourier Transform for TArrayD's
+//
+TArrayD* MFFT::RealFunctionFFT(const TArrayD *data)
+{
+
+  fDim = data->GetSize();
+  CheckDim(fDim);
+
+  fDataD.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = data->At(i);
+
+  RealFTD(1);
+  
+  return new TArrayD(fDim,fDataD.GetArray());
+
+}
+
+//
+// Inverse Fast Fourier Transform for TArrayD's
+//
+TArrayD* MFFT::RealFunctionIFFT(const TArrayD *data)
+{
+
+  fDim = data->GetSize();
+  CheckDim(fDim);
+  
+  fDataD.Set(fDim);
+  //
+  // Clone the array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = data->At(i);
+
+  RealFTD(-1);
+
+  return new TArrayD(fDim,fDataD.GetArray());
+}
+
+
+//
+// Power Spectrum Density Calculation
+//
+TH1D* MFFT::PowerSpectrumDensity(const TH1D *hist)
+{
+
+  TH1D *newhist = (TH1D*)CheckHist(hist,1);
+
+  fDataD.Set(fDim);
+  //
+  // Copy the hist into an array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = hist->GetBinContent(i);
+
+  RealFTD(1);
+
+  Int_t dim2 = fDim*fDim;
+  Double_t c02;
+  Double_t ck2;
+  Double_t cn2;
+  //
+  // Fill the new histogram: 
+  //
+  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
+  //    (stored in fData{0])
+  //
+  c02 = fDataD[0]*fDataD[0];
+  newhist->Fill(c02/dim2);
+  //
+  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)| + |C(N-k)|*|C(N-k)|)
+  //
+  for (Int_t k=2;k<fDim-2;k+=2)
+    {
+
+      Int_t ki  = k+1;
+      ck2 = (fDataD[k]*fDataD[k] + fDataD[ki]*fDataD[ki]);
+      newhist->Fill(ck2/dim2);
+    }
+  //
+  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
+  //    (stored in fData[1])
+  //
+  cn2 = (fDataD[1]*fDataD[1]);
+  newhist->Fill(cn2/dim2);
+
+  return newhist;
+}
+
+
+//
+// Power Spectrum Density calculation
+//
+TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
+{
+
+  TH1F *newhist = (TH1F*)CheckHist(hist,0);
+
+  fDataF.Set(fDim);
+  //
+  // Copy the hist into an array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = hist->GetBinContent(i);
+
+  RealFTF(1);
+
+  Int_t dim2 = fDim*fDim;
+  Float_t c02;
+  Float_t ck2;
+  Float_t cn2;
+  //
+  // Fill the new histogram: 
+  //
+  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
+  //
+  c02 = (fDataF[0]*fDataF[0]);
+  newhist->Fill(c02/dim2);
+  //
+  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
+  //
+  for (Int_t k=2;k<fDim;k+=2)
+    {
+      Int_t ki  = k+1;
+      ck2 = (fDataF[k]*fDataF[k] + fDataF[ki]*fDataF[ki]);
+      newhist->Fill(ck2/dim2);
+    }
+  //
+  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
+  //
+  cn2 = (fDataF[1]*fDataF[1]);
+  newhist->Fill(cn2/dim2);
+  
+  return newhist;
+}
+
+
+void MFFT::CheckDim(Int_t a)
+{
+
+  // If even number, return 0
+  if (a==2)  return;
+
+  // If odd number, return the closest power of 2
+  if (a & 1) 
+    { 
+      Int_t b = 1; 
+      while (b < fDim/2+1)
+        b <<= 1; 
+
+      fDim = b;
+      *fLog << warn << "Dimension of Data is not a multiple of 2, will take only first " 
+            << fDim << " entries! " << endl;
+      return; 
+    }
+
+  CheckDim(a>>1);
+}
+
+TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag)
+{
+  
+  TString name = hist->GetName();
+  name += " PSD";
+  TString title = hist->GetTitle();
+  title += " - Power Spectrum Density";
+
+  // number of entries
+  fDim = hist->GetNbinsX();
+  CheckDim(fDim);
+
+  // Step width
+  Double_t delta = hist->GetBinWidth(1);
+  
+  // Nyquist frequency
+  Axis_t fcrit = 1./(2.*delta);
+  Axis_t low = 0.;
+  Axis_t up  = fcrit;
+
+  switch (flag)
+    {
+    case 0: 
+      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
+      break;
+    case 1:
+      return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up);
+      break;
+    default:
+      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
+      break;
+    }
+}
Index: trunk/MagicSoft/Mars/mtools/MFFT.h
===================================================================
--- trunk/MagicSoft/Mars/mtools/MFFT.h	(revision 2808)
+++ trunk/MagicSoft/Mars/mtools/MFFT.h	(revision 2808)
@@ -0,0 +1,66 @@
+#ifndef MARS_MFFT
+#define MARS_MFFT
+
+#ifndef MARS_MParContainer
+#include "MParContainer.h"
+#endif
+
+#ifndef ROOT_TArrayF
+#include "TArrayF.h"
+#endif
+
+#ifndef ROOT_TArrayD
+#include "TArrayD.h"
+#endif
+
+#ifndef ROOT_TH1F
+#include "TH1F.h"
+#endif
+
+#ifndef ROOT_TH1D
+#include "TH1D.h"
+#endif
+
+class MFFT : public MParContainer
+{
+private:
+
+  void Swap(Float_t &a,  Float_t &b)    { Float_t  c = a;  a = b;  b = c;  }
+  void Swap(Double_t &a, Double_t &b)   { Double_t c = a;  a = b;  b = c;  }
+
+  void TransformF(const Int_t isign);
+  void TransformD(const Int_t isign);  
+  void RealFTF(const Int_t isign);
+  void RealFTD(const Int_t isign);
+
+  void CheckDim(Int_t a);
+  TH1 *CheckHist(const TH1 *hist, const Int_t flag);
+  
+  Int_t   fDim;
+  TArrayF fDataF;
+  TArrayD fDataD;  
+
+public:
+
+  MFFT();
+  ~MFFT();
+
+  TArrayF*  RealFunctionFFT( const TArrayF *data);
+  TArrayF*  RealFunctionIFFT(const TArrayF *data);  
+  
+  TArrayD*  RealFunctionFFT( const TArrayD *data);
+  TArrayD*  RealFunctionIFFT(const TArrayD *data);  
+  
+  Float_t*  RealFunctionFFT( const Int_t n, const Float_t *data);
+  Float_t*  RealFunctionIFFT(const Int_t n, const Float_t *data);  
+  
+  Double_t* RealFunctionFFT( const Int_t n, const Double_t *data);
+  Double_t* RealFunctionIFFT(const Int_t n, const Double_t *data);  
+  
+  TH1F* PowerSpectrumDensity(const TH1F *hist);
+  TH1D* PowerSpectrumDensity(const TH1D *hist);
+  
+  ClassDef(MFFT,0)  // Class to perform a Fast Fourier Transform
+};
+    
+#endif