#include "Riostream.h"

// 3D quadratic function
Double_t QuadFcn3(Double_t *in, Double_t *p) {
  Double_t x = in[0];
  Double_t y = in[1];
  Double_t z = in[2];
  Double_t results = 
    p[0] +
    p[1]*x + p[2]*y + p[3]*z +
    p[4]*x*y + p[5]*x*z + p[6]*y*z + p[7]*x*x + p[8]*y*y + p[9]*z*z;
  return results;
}

void forRootTalk()
{
  gROOT->Reset();
  gStyle->SetPalette(4,0);
  ifstream in;

  Float_t lg,dg1z,dkg, llh;

  // x=dg1z, y=dkg, z=lambda_gamma
  int nbinsX = 3;
  float xmin = -0.15;
  float xmax =  0.15;
  int nbinsY = 3;
  float ymin = -4.5;
  float ymax = 4.5;
  int nbinsZ = 3;
  float zmin = -0.15;
  float zmax =  0.15;


  TH3F *h3 = new TH3F("h3","",nbinsX, xmin, xmax, nbinsY, ymin, ymax, nbinsZ, zmin, zmax);

  // reading data from input file
  in.open("inputFile.txt");

  while (1) {
    if (!(in >> dg1z)) break; 
    if (!(in >> dkg)) break; 
    if (!(in >> lg)) break; 
    if (!(in >> llh)) break;
    if (!in.good()) break;
    printf("dg1z=%3f, dkg=%3f, lg=%3f, llh=%6f \n",dg1z, dkg,lg,llh);

    h3->Fill(dg1z,dkg,lg,-llh);
  }

  gStyle->SetCanvasPreferGL(kTRUE);
  
  TF3 *fitFcn = new TF3("fitFcn",QuadFcn3,-0.15,0.15,-4.5,4.5,-0.15,0.15,10);
  // make 3D quadratic fit on data
  h3->Fit("fitFcn");
  
  Double_t par[10];
  fitFcn->GetParameters(par);
  
  test =new TCanvas("test","contours",600,600);
  
  TF3 *newFcn = new TF3("newFcn",QuadFcn3,-0.15,0.15,-4.5,4.5,-0.15,0.15,10);
  
  for (int i=0;i<10;++i)
    newFcn->SetParameter(i,par[i]);
  
  newFcn->SetTitle(0);
  newFcn->GetXaxis()->SetTitle("x");
  newFcn->GetYaxis()->SetTitle("y");
  newFcn->GetZaxis()->SetTitle("z");

  
  Double_t acpar[1];
  acpar[0]=par[0] +3.91;
  cout << "acpar[0]: " << acpar[0] << endl;
  
  TPad *box2Pad = new TPad("box2", "box2", 0.02, 0.02, 0.98, 0.82);
  box2Pad->Draw();
  box2Pad->cd();
  // make contour on quadratic fit
  newFcn->SetContour(1,acpar); 
  newFcn->Print();

  newFcn->SetFillColor(kRed);
  newFcn->Draw("glsame");
  
}