#include "TMinuit.h"

void Exp_2D() {

// Set ode pares de datos con sus errores ////////////////////////////////////

const Int_t n = 59;

double x[n] = {15,30,45,60,75,90,105,120,135,150,165,180,195,210,225,
      240,255,270,285,300,315,330,345,360,375,390,405,420,435,450,
      465,480,495,510,525,540,555,570,585,600,615,630,645,660,675,
      690,705,720,735,750,765,780,795,810,825,840,855,870,885};
double y[n] = {775,479,380,302,185,157,137,119,110,89,74,61,66,68,48,
	          54,51,46,55,29,28,37,49,26,35,29,31,24,25,35,
              24,30,26,28,21,18,20,27,17,17,14,17,24,11,22,
	          17,12,10,13,16, 9, 9,14,21,17,13,12,18,10};
double ex[n] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
					0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
					0,0,0,0,0,0,0,0,0};
double ey[n] = {27.8,21.9,19.5,17.4,13.6,12.5,11.7,10.9,10.5,9.4,
		8.6,7.8,8.1,8.2,6.9,7.3,7.1,6.8,7.4,5.4,5.3,6.1,7.0,5.1,5.9,
	    5.4,5.6,4.9,5.0,5.9,4.9,5.5,5.1,5.3,4.6,4.2,4.5,5.2,4.1,4.1,
	    3.7,4.1,4.9,3.3,4.7,4.1,3.5,3.2,3.6,4,3,3.0,3.7,4.6,4.1,3.6,
	    3.5,4.2,3.2};
	    	    
////////////////////////////////////////////////////////////////////////
// Set default Minuit in order to use gMinuit
TVirtualFitter::SetDefaultFitter("Minuit");

TCanvas *c1 = new TCanvas("c1","Dos decaimientos exponenciales",300,10,1200,1000);
c1->Divide(1,2);
c1->cd(1);

TGraphErrors *gr = new TGraphErrors(n,x,y,ex,ey);
   
// Cosmetic stuff //////////////////////////////////////////////////////
gr->SetLineColor(4);
gr->SetLineWidth(2);
gr->SetMarkerColor(4);
gr->SetMarkerSize(0.5);
gr->SetMarkerStyle(20);
gr->SetTitle("Dos decaimientos exponenciales;Tiempo (segundos);Cuentas");
gr->Draw("AP");

gPad->SetLogy();
gPad->SetGridx();
gPad->SetGridy();

////////////////////////////////////////////////////////////////////////

gPad->Print("Out_Exp_2D.pdf["); // abre un archivo
gPad->Print("Out_Exp_2D.pdf"); 

////////////////////////////////////////////////////////////////////////
// Defino la función que voy auar para el ajuste ///////////////////////
auto double_exp = [](double *x, double *p) {
   double xx = x[0];
   double f = p[0] + p[1] * exp(-xx/p[3])+ p[2] * exp(-xx/p[4]);
   return f;
};

TF1* f1 = new TF1("double_exp",double_exp,0,900,5);
// Valores iniciales de los parámetros
f1->SetParameters(10,1000,100,40,200); 

// Calculate the best set of parameters
gPad->WaitPrimitive();
TFitResultPtr r = gr->Fit(f1,"S E");
TMatrixDSym sigma = r->GetCovarianceMatrix();
r->Print("V");
gPad->WaitPrimitive();
gPad->Print("Out_Exp_2D.pdf"); 

// Optimizo los valores de los parámetros
double p0_hat=f1->GetParameter(0); 
double p1_hat=f1->GetParameter(1);
double p2_hat=f1->GetParameter(2);
double p3_hat=f1->GetParameter(3);
double p4_hat=f1->GetParameter(4);

// Fijo los parámetros que no juegan en las elipses en sus valores óptimos
f1->FixParameter(0, p0_hat); // parametro [0] ahora esta fijo en p0_hat
f1->FixParameter(1, p1_hat); // parametro [1] ahora esta fijo en p1_hat
f1->FixParameter(2, p2_hat); // parametro [2] ahora esta fijo en p2_hat

gr->Fit(f1); // Ajusta de nuevo

c1->cd(2);

////////////////////////////////////////////////////////////////////////

gMinuit->SetErrorDef(1); // 1-sigma
// Hace uel contorno de 3 y 4 usando 80 puntos
TGraph *gr1b = (TGraph*)gMinuit->Contour(80,3,4); 

gMinuit->SetErrorDef(4); // 2-sigma  /
//¿Porque ponemos 4 para dos sigmas aca?
TGraph *gr2b = (TGraph*)gMinuit->Contour(80,3,4);

gr2b->SetTitle("Regiones de 68.3% y 95.4% CL para los parametros tau1 & tau2;tau1;tau2");

gr2b->SetFillColor(42);
gr2b->Draw("Alf"); // Draw "ellipse" for 2-sigma

gr1b->SetFillColor(38);
gr1b->Draw("Lf"); // Draw "ellipse" for 1-sigma

gPad->SetLogy(0);
gPad->SetGridx();
gPad->SetGridy();

TAxis *axis = gr2b->GetXaxis();
axis->SetLimits(28.,40.);               
gr2b->GetHistogram()->SetMinimum(120.); 
gr2b->GetHistogram()->SetMaximum(320.); 

gPad->Update();
gPad->Print("Out_Exp_2D.pdf"); 

cout<<endl<<"Presionar enter para ver las elipses"<<endl;
getchar();

/////////// Dibuja las elipses ////////////////

double sigma11=sqrt(sigma(3,3)); 
double sigma22=sqrt(sigma(4,4));
double cor12=sigma(3,4)/sqrt(sigma(3,3)*sigma(4,4));

RooEllipse* ell1 = new RooEllipse("Parab CI", p3_hat, p4_hat, 1*sigma11, 1*sigma22, cor12, 100);
RooEllipse* ell2 = new RooEllipse("Parab CI", p3_hat, p4_hat, 2*sigma11, 2*sigma22, cor12, 100);

ell1->SetLineColor(2);
ell2->SetLineColor(2);
ell1->SetLineWidth(2);
ell2->SetLineWidth(2);
ell1->SetLineStyle(1);
ell2->SetLineStyle(1);
ell1->SetFillColor(6);
ell1->Draw("lsame");
ell2->Draw("lsame");

cout<<endl<<"No se parecen ... algo va mal, no? modificá mínimamente el código para arreglarlo!"<<endl<<endl<<endl;

////////////////////////////////////////////

gPad->Print("Out_Exp_2D.pdf"); 
gPad->Print("Out_Exp_2D.pdf]");  // close file

}