#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <cmath>
#include "TColor.h"
#include "TStyle.h"
#include <TFile.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TSystem.h>
#include "TPaletteAxis.h"
#include <TRandom3.h>
#include <TF1.h>
#include <TF2.h>
#include <TMinuit.h>
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TFitter.h"
#include "TVirtualFitter.h"
#include "TRegexp.h"
#include "TMath.h"
#include "TLegend.h"
#include "TROOT.h"
//#include "process_results.h"
#include "TCanvas.h"
using namespace std;


Double_t GausErfBIF1(Double_t *x, Double_t *par)
{
if(par[2] > 0 && par[4] > 0&& par[6]>=0 && TMath::Sqrt((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[3])*(x[1]-par[3]))<=0.5 )
//if(((TMath::Sqrt((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[3])*(x[1]-par[3])))<3)&&par[6]>=0)

   {
  double rx=0.;
  double ry=0.;
  double rz=0.;
  double den=0;
  double g=0;
  rx=(x[0]-par[1])/par[2];
  ry=(x[1]-par[3])/par[4];
  rz=2.*par[5]*rx*ry;
  den= 1.-par[5]*par[5];
  g=TMath::Sqrt( (rx*rx +ry*ry - rz)/(2.*den) );


  double fitval = (par[0]/ (2.*par[6])) * (TMath::Erf(par[6]+g) - TMath::Erf(g-par[6]));

  return fitval;
   }
else
  return 0;
}



Double_t GausErfBIF2(Double_t *x, Double_t *par)
{
if(par[2] > 0 && par[4] > 0&& par[6]>=0 && TMath::Sqrt((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[3])*(x[1]-par[3]))<=3 )
//if(((TMath::Sqrt((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[3])*(x[1]-par[3])))<3)&&par[6]>=0)

   {
  double rx=0.;
  double ry=0.;
  double rz=0.;
  double den=0;
  double g=0;
  rx=(x[0]-par[1])/par[2];
  ry=(x[1]-par[3])/par[4];
  rz=2.*par[5]*rx*ry;
  den= 1.-par[5]*par[5];
  g=TMath::Sqrt( (rx*rx +ry*ry - rz)/(2.*den) );


  double fitval = (par[0]/ (2.*par[6])) * (TMath::Erf(par[6]+g) - TMath::Erf(g-par[6]));

  return fitval;
   }
else
  return 0;
}






void run(){
  

   TF2 *g2D1=new TF2("g2dBIF1",GausErfBIF1,-6.,6.,-6.,6.,7);
   g2D1->FixParameter(0,1.43020e+04);
   g2D1->FixParameter(1,-2.83076e-01);
   g2D1->FixParameter(2,4.23765e-01);
   g2D1->FixParameter(3,-2.15459e-01);
   g2D1->FixParameter(4,4.22232e-01);
   g2D1->FixParameter(5,6.84517e-03);
   g2D1->FixParameter(6,1.58316e+00);
  
   g2D1->Draw("lego2");
   TF2 *g2D2=new TF2("g2dBIF2",GausErfBIF2,-6.,6.,-6.,6.,7);
   
   //while(h2->GetBinContent(106,106)<14000) h2->Fill(r.Gaus(0.,SX),r.Gaus(0.,SY));
   g2D2->FixParameter(0,1.43020e+04);
   g2D2->FixParameter(1,-2.83076e-01);
   g2D2->FixParameter(2,4.23765e-01);
   g2D2->FixParameter(3,-2.15459e-01);
   g2D2->FixParameter(4,4.22232e-01);
   g2D2->FixParameter(5,6.84517e-03);
   g2D2->FixParameter(6,1.58316e+00);
   cout << g2D1->Integral(-3,3,-3,3) <<endl;
   cout << g2D2->Integral(-3,3,-3,3)<<endl;
  cout <<  g2D1->Integral(-3,3,-3,3)/g2D2->Integral(-3,3,-3,3)<<endl;

  
}