/**
* @file TConvolution.C
* @brief .
*
* @author 
*
* @date 
*/
#include <iostream>
#include <sstream>

#include <TError.h>
#include <TVirtualFFT.h>
#include <TComplex.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TCanvas.h>
#include <TString.h>

#define DEBUG 0
#define INFO_OUT_TAG "TConvolution> "
#include <utilities/debugging.hh>
#include "TConvolution.hh"

using std::endl;

ClassImp(TConvolution)

TConvolution::TConvolution(Int_t Ndim, Int_t NbinsX, Int_t NbinsY, 
			   TConvolution::BufStrat bs,
			   Double_t bufFracX, Double_t bufFracY,
			   Option_t* opt)
  : fDim(Ndim),
    fFFTOpt(opt),
    fFFTr2cA(0),
    fFFTr2cB(0),
    fFFTc2r(0),
    fIndex(0),
    fBufStrat(bs),
    fBufFrac(0),
    fNbin(0),
    fNbuf(0),
    fNtot(0),
    fShiftBinA(0),
    fShiftBinB(0),
    fShiftValA(0),
    fShiftValB(0)
{
  // check dimensionality
  if((fDim != 1) && (fDim != 2))
    Error("TConvolution", "Number of dimensions = %d, must be 1 or 2!", fDim);
  
  fIndex = new Int_t[2];

  fNbin = new Int_t[2];
  fNbin[0] = NbinsX;
  fNbin[1] = NbinsY;

  fBufFrac = new Double_t[2];
  fBufFrac[0] = bufFracX;
  fBufFrac[1] = bufFracY;
  DEBUG_OUT << "BufFrac = " << fBufFrac[0] << ", " << fBufFrac[1] << endl;

  fNbuf = new Int_t[2];
  fNtot = new Int_t[2];
  CalcNumberOfBins();
  InitFFT();

  fShiftValA = new Double_t[2];
  fShiftValA[0] = 0;
  fShiftValA[1] = 0;
  fShiftBinA = new Int_t[2];
  fShiftBinA[0] = 0;
  fShiftBinA[1] = 0;

  fShiftValB = new Double_t[2];
  fShiftValB[0] = 0;
  fShiftValB[1] = 0;
  fShiftBinB = new Int_t[2];
  fShiftBinB[0] = 0;
  fShiftBinB[1] = 0;
}

/** 
    \parameter option: choosing how much time should be spent in planning the transform:
    Possible options:
    - "ES" (from "estimate") : no time in preparing the transform, 
    but probably sub-optimal performance
    - "M"  (from "measure") : some time spend in finding the optimal way to do the transform
    - "P" (from "patient") : more time spend in finding the optimal way to do the transform
    - "EX" (from "exhaustive") : the most optimal way is found.
    This option should be chosen depending on how many transforms of the
    same size and type are going to be done. Planning is only done once.
*/
TConvolution::TConvolution(const TH1* histoA, const TH1* histoB, 
			   TConvolution::BufStrat bs,
			   Double_t bufFracX, Double_t bufFracY,
			   Option_t* FFT_option)
  : fDim(histoA->GetDimension()),
    fFFTOpt(FFT_option),
    fFFTr2cA(0),
    fFFTr2cB(0),
    fFFTc2r(0),
    fIndex(0),
    fBufStrat(bs),
    fBufFrac(0),
    fNbin(0),
    fNbuf(0),
    fNtot(0),
    fShiftBinA(0),
    fShiftBinB(0),
    fShiftValA(0),
    fShiftValB(0)
{
  DEBUG_OUT << endl;
  // check dimensionality
  if((fDim != 1) && (fDim != 2))
    Error("TConvolution", "Number of dimensions = %d, must be 1 or 2!", fDim);

  CheckConsistency(histoA, histoB);

  fIndex = new Int_t[2];

  fNbin = new Int_t[2];
  fNbin[0] = histoA->GetNbinsX();
  fNbin[1] = histoA->GetNbinsY();
  DEBUG_OUT << "Nbin = " << fNbin[0] << ", " << fNbin[1] << endl;

  fBufFrac = new Double_t[2];
  fBufFrac[0] = bufFracX;
  fBufFrac[1] = bufFracY;
  DEBUG_OUT << "BufFrac = " << fBufFrac[0] << ", " << fBufFrac[1] << endl;

  fNbuf = new Int_t[2];
  fNtot = new Int_t[2];
  CalcNumberOfBins();
  InitFFT();

  fShiftValA = new Double_t[2];
  fShiftValA[0] = 0;
  fShiftValA[1] = 0;
  fShiftBinA = new Int_t[2];
  CalcShiftBins(histoA, fShiftValA, fShiftBinA);

  fShiftValB = new Double_t[2];
  fShiftValB[0] = 0;
  fShiftValB[1] = 0;
  fShiftBinB = new Int_t[2];
  CalcShiftBins(histoB, fShiftValB, fShiftBinB);

  SetHistoA(histoA);
  SetHistoB(histoB);
}

TConvolution::TConvolution(const TConvolution& other)
  : fDim(other.fDim),
    fFFTOpt(other.fFFTOpt),
    fFFTr2cA(other.fFFTr2cA),
    fFFTr2cB(other.fFFTr2cB),
    fFFTc2r(other.fFFTc2r),
    fIndex(other.fIndex),
    fBufStrat(other.fBufStrat),
    fBufFrac(other.fBufFrac),
    fNbin(other.fNbin),
    fNbuf(other.fNbuf),
    fNtot(other.fNtot),
    fShiftBinA(other.fShiftBinA),
    fShiftBinB(other.fShiftBinB),
    fShiftValA(other.fShiftValA),
    fShiftValB(other.fShiftValB)
{}

TConvolution::~TConvolution()
{
  if(fFFTOpt){ delete fFFTOpt; fFFTOpt = 0; }
  if(fFFTr2cA){ delete fFFTr2cA; fFFTr2cA = 0; }
  if(fFFTr2cB){ delete fFFTr2cB; fFFTr2cB = 0; }
  if(fFFTc2r){ delete fFFTc2r; fFFTc2r = 0; }
  if(fIndex){ delete [] fIndex; fIndex = 0; }
  if(fBufFrac){ delete [] fBufFrac; fBufFrac = 0; }
  if(fNbin){ delete [] fNbin; fNbin = 0; }
  if(fNbuf){ delete [] fNbuf; fNbuf = 0; }
  if(fNtot){ delete [] fNtot; fNtot = 0; }
  if(fShiftBinA){ delete [] fShiftBinA; fShiftBinA = 0; }
  if(fShiftBinB){ delete [] fShiftBinB; fShiftBinB = 0; }
  if(fShiftValA){ delete [] fShiftValA; fShiftValA = 0; }
  if(fShiftValB){ delete [] fShiftValB; fShiftValB = 0; }
}

bool TConvolution::is2D() const
{
  if(!fFFTr2cA) return false;
  return fFFTr2cA->GetNdim() == 2;
}

std::string TConvolution::BufStratToString(BufStrat bs)
{
  std::string sbs;
  switch(bs){
  case None : return "None";
  case Extend : return "Extend";
  case Mirror : return "Mirror";
  case Flat : return "Flat";
  default: return "None";
  }
}

TConvolution::BufStrat TConvolution::StringToBufStrat(std::string sbs)
{
  TString ssbs(sbs);
  ssbs.ToUpper();
  if(ssbs == "NONE") return None;
  else if(ssbs == "EXTEND") return Extend;
  else if(ssbs == "MIRROR") return Mirror;
  else if(ssbs == "FLAT") return Flat;
  else return None;
}

std::string TConvolution::CommaSeparatedListOfBufStrats()
{
  std::stringstream s;
  s << "None,Extend,Mirror,Flat";
  return s.str();
}

void TConvolution::SetBufferStrategy(TConvolution::BufStrat bs) 
{ 
  fBufStrat = bs; 
}

void TConvolution::SetBufferFraction(Double_t bufFracX,
				     Double_t bufFracY)
{
  if(!fBufFrac) fBufFrac = new Double_t[2];

  fBufFrac[0] = bufFracX;
  fBufFrac[1] = bufFracY;
}

void TConvolution::SetShiftA(TH1* histo, 
			     Double_t shiftX, Double_t shiftY)
{
  if(!fShiftValA) fShiftValA = new Double_t[2];
  fShiftValA[0] = shiftX;
  fShiftValA[1] = shiftY;

  // update
  CalcShiftBins(histo, fShiftValA, fShiftBinA);
}

void TConvolution::SetShiftB(TH1* histo, 
			     Double_t shiftX, Double_t shiftY)
{
  if(!fShiftValB) fShiftValB = new Double_t[2];
  fShiftValB[0] = shiftX;
  fShiftValB[1] = shiftY;

  // update
  CalcShiftBins(histo, fShiftValB, fShiftBinB);
}

void TConvolution::SetHistoA(const TH1* histo)
{
  ScanHisto(histo, fShiftBinA, fFFTr2cA);
}

void TConvolution::SetHistoB(const TH1* histo)
{
  ScanHisto(histo, fShiftBinB, fFFTr2cB);
}

void TConvolution::Convolve(TH1* hOut)
{
  MultiplyAndTransform();
  ScanResult(hOut);
}

bool TConvolution::CheckConsistency(const TH1* hA, const TH1* hB)
{
  Int_t Ndim = hA->GetDimension();
  if((Ndim>2) || (hB->GetDimension()>2)){
    Error("CheckConsistency", "Histos must have dimension either 1 or 2.");
    return false;
  }
  if(Ndim != hB->GetDimension()){
    Error("CheckConsistency", "Histos must have same dimension.");
    return false;
  }
    
  if(hA->GetNbinsX() != hB->GetNbinsX()){
    Error("CheckConsistency", "Histos must have same number of X-bins.");
    return false;
  }
  if(hA->GetNbinsY() != hB->GetNbinsY()){
    Error("CheckConsistency", "Histos must have same number of Y-bins.");
    return false;
  }

  return true;
}

/**
   \param transform The Transofrmation to diagnose 
   - "A" is transform A
   - "B" is transform B
   - "P" is the Product transform
 */
void TConvolution::DiagnosticPlot(std::string transform, 
				 std::string name)
{
  if(transform == "A") debugFFT(name.c_str(), fFFTr2cA, true);
  if(transform == "B") debugFFT(name.c_str(), fFFTr2cB, true);
  if(transform == "P") debugFFT(name.c_str(), fFFTc2r, false);
}

void TConvolution::CalcNumberOfBins()
{

  //Nbuf = new Int_t[2];
  fNbuf[0] = Int_t((fNbin[0]*fBufFrac[0])/2 + 0.5) ;
  fNbuf[1] = Int_t((fNbin[1]*fBufFrac[1])/2 + 0.5) ;
  if(fBufStrat == None || fBufStrat == Extend){
    fNbuf[0] = 0 ;
    fNbuf[1] = 0 ;
  }
  DEBUG_OUT << "Nbuf = " << fNbuf[0] << ", " << fNbuf[1] << endl;

  //NbinsTotal = new Int_t[2];
  fNtot[0] = fNbin[0] + 2*fNbuf[0] ;
  fNtot[1] = fNbin[1] + 2*fNbuf[1] ;
  DEBUG_OUT << "Ntot = " << fNtot[0] << ", " << fNtot[1] << endl;
}

void TConvolution::InitFFT()
{
  TString rsopt = TString("R2CK") + TString(fFFTOpt);
  fFFTr2cA = TVirtualFFT::FFT(fDim, fNtot, rsopt.Data());
  fFFTr2cB = TVirtualFFT::FFT(fDim, fNtot, rsopt.Data());
  TString csopt = TString("C2RK") + TString(fFFTOpt);
  fFFTc2r  = TVirtualFFT::FFT(fDim, fNtot, csopt.Data());
}

void TConvolution::CalcShiftBins(const TH1* histo, 
				 const Double_t* ShiftVal,
				 Int_t* zeroBin)
{
  // Find bin ID that contains zero value
  //zeroBin = new Int_t[2]; 
  zeroBin[0] = 0;
  zeroBin[1] = 0;
  Int_t* binShift = new Int_t[2]; 
  binShift[0] = 0;
  binShift[1] = 0;

//   // find center bin
//   if(!%2) zeroBin[0] = fNbin[0]/2;
//   else zeroBin[0] = 1+fNbin[0]/2;
//   if(!NbinsY%2) centerBinNumY = NbinsY/2;
//   else centerBinNumY = 1+NbinsY/2;

  // X
  Double_t Xmin = histo->GetXaxis()->GetXmin();
  Double_t Xmax = histo->GetXaxis()->GetXmax();
  Double_t dX = Xmax - Xmin;
  Double_t bw = (Xmax - Xmin) / fNtot[0] ;
  if (Xmax>=0 && Xmin<=0) {
    zeroBin[0] = histo->GetXaxis()->FindFixBin(0) ;
  } else if (Xmin>0) {
    zeroBin[0] = Int_t(-Xmin/bw) ;
  } else {
    zeroBin[0] = Int_t(-1*Xmax/bw) ;
  }

  // bin shift
  binShift[0] = Int_t((fNtot[0] * ShiftVal[0]) / dX) ;
  zeroBin[0] += binShift[0] ;
  while(zeroBin[0]>=fNtot[0]) zeroBin[0] -= fNtot[0] ;
  while(zeroBin[0]<0) zeroBin[0] += fNtot[0] ;

  // Y
  if(is2D()){
    Double_t Ymin = histo->GetYaxis()->GetXmin();
    Double_t Ymax = histo->GetYaxis()->GetXmax();
    Double_t dY = Ymax - Ymin;
    Double_t bw = (Ymax - Ymin) / fNtot[1] ;
    if (Ymax>=0 && Ymin<=0) {
      zeroBin[1] = histo->GetYaxis()->FindFixBin(0) ;
    } else if (Ymin>0) {
      zeroBin[1] = Int_t(-Ymin/bw) ;
    } else {
      zeroBin[1] = Int_t(-1*Ymax/bw) ;
    }
    // bin shift
    binShift[1] = Int_t((fNtot[1] * ShiftVal[1]) / dY) ;
    zeroBin[1] += binShift[1] ;
    while(zeroBin[1]>=fNtot[1]) zeroBin[1] -= fNtot[1] ;
    while(zeroBin[1]<0) zeroBin[1] += fNtot[1] ;
  }

  DEBUG_OUT << "zeroBin = " << zeroBin[0] 
	    << ", " << zeroBin[1] << endl;
  DEBUG_OUT << "binShift = " << binShift[0] 
	    << ", " << binShift[1] << endl;

  delete [] binShift;
}


void TConvolution::ScanHisto(const TH1* histo,
			     Int_t* zeroBin, 
			     TVirtualFFT* fft)
{
  // check consistancy
  if(!fft){
    Error("ScanHisto", "FFT object not found.");
    return;
  }
  if(histo->GetDimension() != fft->GetNdim()){
    Error("ScanHisto", "Dimensionally inconsistant objects.");
    return;
  }
  if(histo->GetNbinsX()+2*fNbuf[0] != fft->GetN()[0]){
    Error("ScanHisto", "Inconsistant number of X bins.");
    return;
  }
  if(histo->GetDimension() == 2){
    if(histo->GetNbinsY()+2*fNbuf[0] != fft->GetN()[1]){
      Error("ScanHisto", "Inconsistant number of Y bins.");
      return;
    }
  }


  // Allocate array of sampling size plus optional buffer zones
  Int_t NallBins = fNtot[0]*fNtot[1];
  // First scan hist into temp array 
  Double_t * tmp = new Double_t[NallBins] ;
  Int_t bin = 0;
  Double_t val = 0;

  // center
  for (Int_t i=fNbuf[0]; i<fNtot[0]-fNbuf[0]; i++) {
    for (Int_t j=fNbuf[1]; j<fNtot[1]-fNbuf[1]; j++) {
      bin = i + fNtot[0]*j;
      val = histo->GetBinContent(i+1-fNbuf[0], j+1-fNbuf[1]);
      tmp[bin] = val;
    }
  }

  //
  // FLAT
  //
  if(fBufStrat == Flat){
 
    // left 
    for (Int_t i=0; i<fNbuf[0]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(1, 1);
	tmp[bin] = val;
      }
      // middle 
      for (Int_t j=fNbuf[1]; j<fNtot[1]-fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(1, j+1-fNbuf[1]);
	tmp[bin] = val;
      }
      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(1, fNbin[1]);
	tmp[bin] = val;
      }
    }

    // right 
    for (Int_t i=fNbin[0]+fNbuf[0]; i<fNtot[1]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0], 1);
	tmp[bin] = val;
      }
      // middle 
      for (Int_t j=fNbuf[1]; j<fNtot[1]-fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0], j+1-fNbuf[1]);
	tmp[bin] = val;
      }
      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0], fNbin[1]);
	tmp[bin] = val;
      }
    }

    // middle 
    for (Int_t i=fNbuf[0]; i<fNbin[0]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(i+1-fNbuf[0], 1);
	tmp[bin] = val;
      }
      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(i+1-fNbuf[0], fNbin[1]);
	tmp[bin] = val;
      }
    }
  } // end flat


  //
  // Mirror
  //
  if(fBufStrat == Mirror){
 
    // left 
    DEBUG_OUT_L(2) << "left: " << std::endl;
    for (Int_t i=0; i<fNbuf[0]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbuf[0]-i, fNbuf[1]-j);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << fNbuf[0]-i
		       << ", biny = " << fNbuf[1]-j << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }

      // middle 
      for (Int_t j=fNbuf[1]; j<fNtot[1]-fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbuf[0]-i, j+1-fNbuf[1]);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << fNbuf[0]-i
		       << ", biny = " << j+1-fNbuf[1] << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }

      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbuf[0]-i, 
				   fNbin[1]-(j-fNtot[1]+fNbuf[1])); 
	tmp[bin] = val;
 	DEBUG_OUT_L(2) << "binx = " << fNbuf[0]-i
		       << ", biny = " << fNbin[1]-(j-fNtot[1]+fNbuf[1]) << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }
    }

    // right 
    DEBUG_OUT_L(2) << "right: " << std::endl;
    for (Int_t i=fNtot[0]-fNbuf[0]; i<fNtot[0]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0]-(i-fNtot[0]+fNbuf[0]), 
				   1);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << fNbin[0]-(i-fNtot[0]+fNbuf[0])
		       << ", biny = " << 1 << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }

      // middle 
      for (Int_t j=fNbuf[1]; j<fNtot[1]-fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0]-(i-fNtot[0]+fNbuf[0]), 
				   j+1-fNbuf[1]);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << fNbin[0]-(i-fNtot[0]+fNbuf[0])
		       << ", biny = " << j+1-fNbuf[1] << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }

      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(fNbin[0]-(i-fNtot[0]+fNbuf[0]), 
				   fNbin[1]);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << fNbin[0]-(i-fNtot[0]+fNbuf[0])
		       << ", biny = " << fNbin[1] << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }
    }

    // middle 
    DEBUG_OUT_L(2) << "middle: " << std::endl;
    for (Int_t i=fNbuf[0]; i<fNbin[0]; i++) {
      // bottom 
      for (Int_t j=0; j<fNbuf[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(i+1-fNbuf[0], j+1-fNbuf[1]);
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << i+1-fNbuf[0]
		       << ", biny = " << j+1-fNbuf[1] << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }

      // top 
      for (Int_t j=fNtot[1]-fNbuf[1]; j<fNtot[1]; j++) {
	bin = i + fNtot[0]*j;
	val = histo->GetBinContent(i+1-fNbuf[0], 
				   fNbin[1]-(j-fNtot[1]+fNbuf[1]));
	tmp[bin] = val;
	DEBUG_OUT_L(2) << "binx = " << i+1-fNbuf[0]
		       << ", biny = " << fNbin[1]-(j-fNtot[1]+fNbuf[1]) << std::endl;
	DEBUG_OUT_L(2) << "tmp[" << bin 
		       << "] = " << val << std::endl;
      }
    }
  } // end mirror



  //
  // Scan function and store values in FFT
  //
  //   Double_t* array = new Double_t[NallBins] ;
  for (Int_t binx=0 ; binx<fNtot[0] ; binx++) {
    // Cyclically shift writing location by zero bin position    
    Int_t j = binx - (zeroBin[0]) ;
    if (j<0) j += fNtot[0] ;
    if (j>=fNtot[0]) j -= fNtot[0] ;
    
    for (Int_t biny=0 ; biny<fNtot[1] ; biny++) {
      // Cyclically shift writing location by zero bin position    
      Int_t k = biny - (zeroBin[1]) ;
      if (k<0) k+= fNtot[1] ;
      if (k>=fNtot[1]) k-= fNtot[1] ;

      bin = binx + fNtot[0]*biny;
      //       array[bin] = tmp[j+fNtot[0]*k] ;
      fft->SetPoint(bin, tmp[j+fNtot[0]*k]);
    }
  }

  delete [] tmp;

  //
  // Transform
  //
  fft->Transform();
}

void TConvolution::complexProduct(Double_t re1, Double_t im1, 
				  Double_t re2, Double_t im2,
				  Double_t& re, Double_t& im)
{
  re = re1*re2 - im1*im2 ;
  im = re1*im2 + re2*im1 ;
}

void TConvolution::MultiplyAndTransform()
{
  if(!fFFTc2r || !fFFTr2cA || !fFFTr2cB) 
    return;

  // multiply
  Double_t re1,re2,im1,im2,re,im;

  if(is2D()){
    for (Int_t i=0 ; i<fFFTc2r->GetN()[0] ; i++) {
      fIndex[0] = i; 
      for (Int_t j=0 ; j<fFFTc2r->GetN()[1]/2+1 ; j++) {
	fIndex[1] = j;
	fFFTr2cA->GetPointComplex(fIndex,re1,im1) ;
	fFFTr2cB->GetPointComplex(fIndex,re2,im2) ;
	complexProduct(re1, im1, re2, im2, re, im);
	fFFTc2r->SetPoint(fIndex,re,im) ;
      }
    }
  }

  else{
    fIndex[1] = 0;
    for (Int_t i=0 ; i<fFFTc2r->GetN()[0]/2+1 ; i++) {
      fIndex[0] = i; 
      fFFTr2cA->GetPointComplex(fIndex,re1,im1) ;
      fFFTr2cB->GetPointComplex(fIndex,re2,im2) ;
      complexProduct(re1, im1, re2, im2, re, im);

      TComplex t(re,im) ;
      fFFTc2r->SetPointComplex(fIndex[0],t) ;
      //fFFTc2r->SetPoint(fIndex,re,im) ;
    }
  }

  // transform
  fFFTc2r->Transform();
}

void TConvolution::ScanResult(TH1* histo) const
{
  Int_t Ndim = fFFTc2r->GetNdim();

  if(!histo || histo->GetDimension()!=Ndim){
    Int_t * n = fFFTc2r->GetN();
    TString name = "TranformedHisto";
    if (Ndim==1)
      histo = new TH1D(name, name, n[0], 0, n[0]);
    else if (Ndim==2)
      histo = new TH2D(name, name, n[0], 0, n[0], n[1], 0, n[1]);
  }

  Int_t* totShift = new Int_t[2];
//   totShift[0] = fShiftBinA[0] + (fNtot[0]-fNbin[0])/2;
//   totShift[1] = fShiftBinA[1] + (fNtot[1]-fNbin[1])/2;
  totShift[0] = fShiftBinA[0]; // + fNbuf[0];
  totShift[1] = fShiftBinA[1]; // + fNbuf[1];
  DEBUG_OUT << "totShift = " << totShift[0] << ", " << totShift[1] << endl;

  Int_t Nall = fNtot[0]*fNtot[1];
  //if(!fIsRenorm) Nall = 1;
  
  Double_t val = 0;
  for (Int_t binx = 0 ; binx<fNbin[0] ; binx++) {
    // Cyclically shift array back so that bin containing zero is back in fShiftBinA
    Int_t j = binx + totShift[0] ;
    while (j<0) j += fNtot[0] ;
    while (j>=fNtot[0]) j -= fNtot[0] ;
    //DEBUG_OUT << binx << ", " << zeroBin[0] << ", " << j << endl;

    for (Int_t biny = 0 ; biny<fNbin[1] ; biny++) {
      // Cyclically shift array back so that bin containing zero is back in fShiftBinA
      Int_t k = biny + totShift[1] ;
      while (k<0) k += fNtot[1] ;
      while (k>=fNtot[1]) k -= fNtot[1] ;

      //bin = binx + fNtot[0]*biny;
      fIndex[0] = j; fIndex[1] = k; //j+fNtot[0]*k
      val = fFFTc2r->GetPointReal(fIndex) / Nall;

      if(is2D()) histo->SetBinContent(binx+1, biny+1, val);
      else histo->SetBinContent(binx+1, val);
    }
  }

  delete [] totShift;
}

TH1 * TConvolution::GetTransformedHisto(const char* name, const char* title, 
					TVirtualFFT* fft, Option_t* option)
{
  TH1 * histo = 0;
  if(!fft) return histo;

  if(is2D()) histo = new TH2D(name, title, 
			      fft->GetN()[0], 0, fft->GetN()[0]-1, 
			      fft->GetN()[1], 0, fft->GetN()[1]-1);
  else histo  = new TH1D(name, title, 
			 fft->GetN()[0], 0, fft->GetN()[0]-1);

  TH1::TransformHisto(fft, histo, option);
  histo->SetStats(kFALSE);
  histo->GetXaxis()->SetTitle("Bin Number");
  if(is2D()) histo->GetYaxis()->SetTitle("Bin Number");

  return histo;
}

void TConvolution::debugFFT(const char* name, TVirtualFFT * fft, 
			   bool isComplex)
{
  if(!fft) return;

  char hName[512];
  char hTitle[512];
  sprintf(hName, "c%sFFT", name);
  sprintf(hTitle, "Fourier Space of %s", name);
  TCanvas * cFFT = 0;

  char drawopt[128] = "";
  if(is2D()) sprintf(drawopt, "COLZ");

  if(isComplex){
    cFFT = new TCanvas(hName, hTitle);
    cFFT->Divide(2,2);
    cFFT->cd(1);
      sprintf(hName, "hRE%s", name);
      sprintf(hTitle, "FFT of %s: Real Part", name);
      TH1 * hRE = GetTransformedHisto(hName, hTitle, fft, "RE");
      hRE->Draw(drawopt);
    cFFT->cd(2);
      sprintf(hName, "hIM%s", name);
      sprintf(hTitle, "FFT of %s: Imaginary Part", name);
      TH1 * hIM = GetTransformedHisto(hName, hTitle, fft, "IM");
      hIM->Draw(drawopt);
    cFFT->cd(3);
      sprintf(hName, "hMAG%s", name);
      sprintf(hTitle, "FFT of %s: Magnitude", name);
      TH1 * hMAG = GetTransformedHisto(hName, hTitle, fft, "MAG");
      hMAG->Draw(drawopt);
    cFFT->cd(4);
      sprintf(hName, "hPH%s", name);
      sprintf(hTitle, "FFT of %s: Phase", name);
      TH1 * hPH = GetTransformedHisto(hName, hTitle, fft, "PH");
      hPH->Draw(drawopt);
  }
  else{
    cFFT = new TCanvas(hName, hTitle, 610, 600);
    sprintf(hName, "hRE%s", name);
    sprintf(hTitle, "FFT of %s: Real Part", name);
    TH1 * hRE = GetTransformedHisto(hName, hTitle, fft, "RE");
    hRE->Draw(drawopt);
  }

  sprintf(hName, "cDebugFFT%s.png", name);
  cFFT->SaveAs(hName);
}

std::ostream& operator<<(std::ostream& os, const TConvolution::BufStrat& bs)
{
  os << TConvolution::BufStratToString(bs).c_str();
  return os;
}