void Gauss2D() {

  // Output file
  TString filename;
  filename = Form("Gauss2D.root");
  TFile  *ifile = (TFile*) gROOT->GetListOfFiles()->FindObject(filename.Data());
  if (!ifile) ifile = new TFile(filename,"rewrite");

  ifile->cd();

  // Definition of the charge density
  Int_t NBinsZ = 200;
  Float_t zmin = -20;
  Float_t zmax = 20;
  Float_t zrange = zmax-zmin;

  Int_t NBinsX = 200;
  Float_t xmin = -20;
  Float_t xmax = 20;
  Float_t xrange = xmax-xmin;

  Float_t kzmin = 0.;
  Float_t kzmax = (TMath::TwoPi() / zrange);
  // Float_t kzmax = (1.0 / zrange);
  
  Float_t kxmin = 0.;
  Float_t kxmax = (TMath::TwoPi() / xrange);
  // Float_t kxmax = (1.0 / xrange);

  // Charge density function:
  TF2 *fRho2D = new TF2("fRho2D","TMath::Gaus(x)*TMath::Gaus(y)",zmin,zmax,xmin,xmax);
  fRho2D->SetNpx(NBinsZ);
  fRho2D->SetNpy(NBinsX);
  // --------------------------------------

  //  Fast Fourier transform
  cout << Form(" Filling data for charge density ..." ) << endl;
  Int_t dims[2] = {NBinsZ,NBinsX};
  //  Double_t *data = new Double_t[NBinsZ*NBinsX];
  Double_t *data = new Double_t[NBinsZ*2*(NBinsX/2+1)];
  // Extra padding! (see http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data) 
  
  for(Int_t i=0; i<NBinsZ; i++) {
    
    for(Int_t j=0; j<NBinsX; j++) {

      Int_t index =  i * NBinsX + j;
      
      data[index]  = fRho2D->Eval((zrange/NBinsZ) * (i-(NBinsZ/2)+0.5) ,
				  (xrange/NBinsX) * (j-(NBinsX/2)+0.5));
      
      // Add background charge
      // data[index] -= 1;
    }
  }
  
  cout << Form("   done!" ) << endl;
  
  cout << Form(" Fourier transform ...") << endl;
  TVirtualFFT *fft = TVirtualFFT::FFT(2, dims, "R2C ES K");
  fft->SetPoints(data);
  fft->Transform();  
  cout << Form("   done!" ) << endl;

  TH2F *hFFTRE = 0;
  hFFTRE = (TH2F*) TH1::TransformHisto(fft, hFFTRE, "RE");
  hFFTRE->SetName("hFFTRE");
  hFFTRE->Scale(1.0/TMath::Sqrt(NBinsZ*NBinsX));
  hFFTRE->SetBins(NBinsZ,kzmin,kzmax,NBinsX,kxmin,kxmax);
  hFFTRE->GetZaxis()->SetRangeUser(-1,1);
  
  TH2F *hFFTIM = 0;
  hFFTIM = (TH2F*) TH1::TransformHisto(fft, hFFTIM, "IM");
  hFFTIM->SetName("hFFTIM");
  hFFTIM->Scale(1.0/TMath::Sqrt(NBinsZ*NBinsX));    
  hFFTIM->SetBins(NBinsZ,kzmin,kzmax,NBinsX,kxmin,kxmax);
  hFFTIM->GetZaxis()->SetRangeUser(-1,1);

  // -----------------------------------------------------
 
  cout << Form(" Solving Gauss equation in Fourier space" )  << endl;
  // Equation: \phi(kz,kx) = \rho(kz,kx)/(kz^2 + kx^2) 
  
  fft->GetPoints(data);    
  for(Int_t i=0; i<dims[0]; i++) {
    for(Int_t j=0; j<(dims[1]/2+1); j++) {

      Int_t index =  2 * (i * (dims[1]/2+1)  + j);
      //  cout << Form(" i = %4i  j = %4i  index = %10i",i,j,index) << endl;
      
      Float_t kz;
      if(i<dims[0]/2)
	kz =  kzmax * (i);
      else
	kz =  kzmax * (i-dims[0]);     

      Float_t kx = kxmax * (j) ;
      
      Float_t k2 = kx*kx + kz*kz;

      if(k2==0) {
	data[index] = 0.0;
	data[index+1] = 0.0;      
      } else {
	data[index] /= k2;
	data[index+1] /= k2;      
      }
    }
  }
  cout << Form("   done!" ) << endl;
  
  cout << Form(" Inverse Fourier transform ..." ) << endl;
  TVirtualFFT *fft_back = TVirtualFFT::FFT(2, dims, "C2R ES K");
  fft_back->SetPoints(data);
  fft_back->Transform();  
  cout << Form("   done!" ) << endl;

  // Get Phi data directly from the fft_back object
  TH2F *hPhi2D = 0;
  hPhi2D = (TH2F*) TH1::TransformHisto(fft_back, hPhi2D, "RE");
  hPhi2D->SetName("hPhi2D");
  hPhi2D->Scale(1.0/(NBinsZ*NBinsX));    
  hPhi2D->SetBins(NBinsZ,zmin,zmax,NBinsX,xmin,xmax);

  // Calculate electric fields
    
  // Alternative way to get Phi
  // TH2F *hPhi2D = new TH2F("hPhi2D","",NBinsZ,zmin,zmax,NBinsX,xmin,xmax);

  cout << Form(" Take the gradient of Phi ..." ) << endl;
  // Directly from the data array

  // Back fourier transformed data (real)
  Double_t *re_back = fft_back->GetPointsReal();
  
  // Ez: Derivative in z-dimension
  TH2F *hEz2D = new TH2F("hEz2D","",NBinsZ,zmin,zmax,NBinsX,xmin,xmax);
  Double_t dz = hPhi2D->GetXaxis()->GetBinWidth(1);
  for(Int_t j=0; j<NBinsX; j++) {
    for(Int_t i=0; i<NBinsZ; i++) {
      Double_t der = 0.;
      if(i>2 && i<NBinsZ-2) {
	Int_t indip1 =  (i+1) * NBinsX + j;
	Int_t indip2 =  (i+2) * NBinsX + j;
	Int_t indim1 =  (i-1) * NBinsX + j;
	Int_t indim2 =  (i-2) * NBinsX + j;
	
	der =  ( 4.0 / 3.0 * (re_back[indip1] - re_back[indim1]) / (2.0 * dz)
	       - 1.0 / 3.0 * (re_back[indip2] - re_back[indim2]) / (4.0 * dz) );
	
      }
      
      hEz2D->SetBinContent(i+1,j+1,-der);
    }
  }
  hEz2D->Scale(1.0/(NBinsZ*NBinsX));

  // Ex: Derivative in x-dimension
  TH2F *hEx2D = new TH2F("hEx2D","",NBinsZ,zmin,zmax,NBinsX,xmin,xmax);
  Double_t dx = hPhi2D->GetYaxis()->GetBinWidth(1);
  for(Int_t i=0; i<NBinsZ; i++) {
    for(Int_t j=0; j<NBinsX; j++) {
      
      // Int_t index =  i * NBinsX + j; 
      // hPhi2D->SetBinContent(i+1,j+1,re_back[index]);
      
      Double_t der = 0.;
      if(j>2 && j<NBinsX-2) {
	Int_t indjp1 =  i * NBinsX + (j+1);
	Int_t indjp2 =  i * NBinsX + (j+2);
	Int_t indjm1 =  i * NBinsX + (j-1);
	Int_t indjm2 =  i * NBinsX + (j-2);
	
	der =  ( 4.0 / 3.0 * (re_back[indjp1] - re_back[indjm1]) / (2.0 * dx)
	       - 1.0 / 3.0 * (re_back[indjp2] - re_back[indjm2]) / (4.0 * dx) );
	
      }
      
      hEx2D->SetBinContent(i+1,j+1,-der);
    }
  }
  hEx2D->Scale(1.0/(NBinsZ*NBinsX));


  cout << Form("   done!" ) << endl;

  gStyle->SetPalette(kBird);

  gStyle->SetCanvasPreferGL(kTRUE);
  hEz2D->Draw("glsurf2");

  // Analytical formula of E field on-axis
  TF1 *f = new TF1("f","(1.0-TMath::Exp(-x*x/2.0))/x",xmin,xmax);
  f->SetNpx(NBinsZ);
  
  TH1F *hEz1D_0 = (TH1F*) hEz2D->ProjectionX("hEz1D_0",NBinsX/2,NBinsX/2);
  
  TH1F *hEx1D_0 = (TH1F*) hEx2D->ProjectionY("hEx1D_0",NBinsZ/2,NBinsZ/2);

  // Check with analytical solution
  // f->Draw();
  // hEz1D_0->Draw("same L");
  // hEx1D_0->Draw("same L");

}