{
  //Define Histogram binning and limits
  const Int_t nDims = 2;
  Int_t nBins[nDims] = {100,100};
  Double_t min[nDims] = {1,1};
  Double_t max[nDims] = {nBins[0],nBins[1]};
  
  //Define 2D Function to Plot
  TF2 *f2 = new TF2("f2","sin([0]*x)*sin([1]*y)",min[0],max[0],min[1],max[1]);
  f2->SetParameter(0,0.25);  //frequency factor in x
  f2->SetParameter(1,0.5);   //frequency factor in y


  //Declare 2D Histogram
  TH1 *hist = new TH2F("hist","hist",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);

  Double_t x,y;
  Int_t gbin; 
  for (Int_t i=0; i<=nBins[0]+1; i++) {
    for (Int_t j=0; j<=nBins[1]; j++) {
      x = (Double_t(i));
      y = (Double_t(j));

      gbin = hist->GetBin(i,j);

      hist->SetBinContent(gbin, f2->Eval(x,y));
    }
  }
  TCanvas *Tcanvas0 = new TCanvas("Tcanvas0", "image",0,0,800,600);
  hist->Draw("COLZ");
  

  //Perform FFT and Draw Magnitude
  //This histogram output is as expected and understood
  TH1 *hm =0;
  TVirtualFFT::SetTransform(0);
  hm = hist->FFT(hm, "MAG");
  
  TCanvas *Tcanvas1 = new TCanvas("Tcanvas1", "ROOT mag_coef",0,0,800,600);
  hm->SetTitle("Magnitude of the 1st transform");
  hm->Draw("COLZ");
  

  //Dump coefficients into arrays
  //These arrays and the following histograms are not understood...
  //...Perhaps histograms not filled properly
  Double_t *re_full = new Double_t[nBins[0]][nBins[1]];
  Double_t *im_full = new Double_t[nBins[0]][nBins[1]];
  TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
  fft->GetPointsComplex(re_full,im_full);
  
  //Fill 2D histogram with real coefficients array
  TH1 *re_coef = new TH2F("re_coef","re_coef",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);
  for (int i=0; i<nBins[0]; i++) {
    for (int j=0; j<nBins[1]; j++) {
      re_coef->SetBinContent(i+1,j+1, re_full[i+1][j+1]);
    }
  }
  TCanvas *Tcanvas2 = new TCanvas("Tcanvas2", "re_coef",0,0,800,600);
  re_coef->SetTitle("Re_Coefs"); 
  re_coef->Draw("COLZ");
  
  //Fill 2D histogram with imaginary coefficients array
  TH1 *im_coef = new TH2F("im_coef","im_coef",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);
  for (int i=0; i<nBins[0]; i++) {
    for (int j=0; j<nBins[1]; j++) {
      im_coef->SetBinContent(i+1,j+1, im_full[i+1][j+1]);
    }
  }
  TCanvas *Tcanvas3 = new TCanvas("Tcanvas3", "im_coef",0,0,800,600);
  im_coef->SetTitle("Im_Coefs");
  im_coef->Draw("COLZ");
  
  //Calculate magnitude using above histograms
  TH1 *mag_coef = new TH2F("mag_coef","mag_coef",nBins[0],min[0],max[0],nBins[1],min[1],max[1]);
  float re,im;
  for (int i=1; i<=nBins[0]; i++) {
    for (int j=1; j<=nBins[1]; j++) {
      re = re_coef->GetBinContent(i,j);
      im = im_coef->GetBinContent(i,j);

      mag_coef->SetBinContent(i+1,j+1, sqrt(re*re+im*im));
    }
  }
  TCanvas *Tcanvas4 = new TCanvas("Tcanvas4", "mag_coef",0,0,800,600);
  mag_coef->SetTitle("Mag_Coefs");
  mag_coef->Draw("COLZ");


  
  //Now let's make a backward transform:
  TVirtualFFT *fft_back = TVirtualFFT::FFT(nDims,nBins, "C2R M K");
  fft_back->SetPointsComplex(re_full,im_full);
  fft_back->Transform();
  
  TH1 *hb = 0;
  //Let's look at the output
  hb = TH1::TransformHisto(fft_back,hb,"Re");

  TCanvas *Tcanvas5 = new TCanvas("Tcanvas5", "iFFT image",0,0,800,600);
  hb->SetTitle("The backward transform result");
  hb->SetBins(nBins[0],min[0],max[0],nBins[1],min[1],max[1]);
  hb->Scale(1./(nBins[0]*nBins[1]));
  hb->Draw("COLZ");
  
}