#include <iostream>
#include <fstream>
#include <sstream>
#include <stdio.h>
#include <string>
#include <vector>
#include <map>
#include <iomanip>
#include <sstream>
#include <time.h>
#include "TF1.h"
#include "TF2.h"
#include "TH1.h"
#include "TH2.h"
#include "TAxis.h"
#include "TMath.h"
#include "TEllipse.h"
#include "TStyle.h"
#include "TFile.h"
#include "TColor.h"
#include "TSpectrum.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TMinuit.h"

// This code does the following :
//  - creates a 2-D function of an annulus.
//  - creates a histogram and fills it randomly via the function.  
//  - draws the function and the histogram as a surf2 plot. 
// This example can be executed via the interpreter or ACLIC
//   root > .x rotateAnnulusFunc.C
//   root > .x rotateAnnulusFunc.C++


 Double_t g3(Double_t *x, Double_t *par) {                                             
  Double_t x1 = Double_t((x[0]-par[7])*cos(par[3])+(x[1]-par[8])*sin(par[3])); 
  Double_t x2 = Double_t(-(x[0]-par[7])*sin(par[3])+(x[1]-par[8])*cos(par[3])); 
  Double_t  a = Double_t(par[1]); 
  Double_t  b = Double_t(par[2]); 
  Double_t r1 = Double_t((x1)*(x1)/(a*a));
  Double_t r2 = Double_t((x2)*(x2)/(b*b));     
  Double_t A1 = Double_t((r1+r2-par[5])/par[4]);
  Double_t A2 = Double_t((r1+r2-par[6])/par[4]);
  //return par[0]*(+1/(exp(A1)+1));
  //return par[0]*(-(1/(exp(A2)+1))); 
  return par[0]*(-(1/(exp(A2)+1))+1/(exp(A1)+1));  
}

 Double_t fun3(Double_t *x, Double_t *par) {                     // Both x and par are vectors. 
    Double_t *p1 = &par[0];                                      // p1 is a pointer to the address of parameter 0. 
    Double_t result = g3(x,p1); 
    return result;
 }

void rotateAnnulusFunc() {
    gStyle->SetPalette(1);
    Double_t w = 1400;
    Double_t h = 700;
    TCanvas * can_fun = new TCanvas("Puck1", "Puck2", w, h );         // name, title, width, height. 
    can_fun  ->Divide(2,1); 
    Double_t pi = 3.14159265359; 
    const Int_t npar = 9;
    Double_t clip = 1.; 
    Double_t x0 = 2.8; Double_t x1 = clip*13.2; 
    Double_t y0 = 4.2; Double_t y1 = clip*13.6; 
   TF2 *f3 = new TF2("f3",fun3,x0,x1,y0,y1, npar);    
   Double_t fitPass[npar] =
   {clip*4., clip*1.5, clip*1.0, pi/8., clip*0.3,clip*1.5,clip*1.3, clip*6.3, clip*7.}; //2*pi/180*60.
   //amp[0]   , a[1] ,b[2]     , theta[3]    , kbT[4]   , r_out[5], r_in[6], x0[6]  , y0[7]     
   f3->SetParameters(fitPass);
   //Create an histogram and fill it randomly with f3
   TH2F *h3 = new TH2F("h3","from f3",100,x0,x1,100,y0, y1);
   Int_t nentries = 10000000;
   h3->FillRandom("f3",nentries);
   // Write histogram into root file. 
   char nameFile[60];
   sprintf(nameFile, "%d_TransformAnnulus.root", 22222);
   TFile* file1 = new TFile( nameFile, "RECREATE");     
   h3 -> Write(); 
   file1 -> Close(); 
   delete file1; 
   // Read histogram from root file. 
   char nameFile1[30];
   sprintf(nameFile1, "%d_TransformAnnulus.root", 22222);
   TFile* file = new TFile(nameFile1, "READ"); 
   TH2F * histo = (TH2F*)file->Get("h3"); 
   can_fun -> cd(1); 
   //histo -> Draw("surf2"); 
   f3->SetNpx(10000); 
   f3->Draw("surf2"); 
   can_fun -> cd(2); 
   h3 -> Draw("surf2"); 
}