#include "TMinuit.h"

void LoopGenerate_2D (int Ntoys = 10) {

// Datos a Fitear
auto *Hdata  = new TH1D("Hdata" ,"Dos fuentes radioactivas",60,0,900);
Hdata->SetTitle("Dos exponenciales;Tiempo (seg);Cuentas / 15 seg");

// Funcion con la que fitear
TF1* fTrue = new TF1("fTrue","[0]+[1]*exp(-x/[3])+[2]*exp(-x/[4])",0,900);

const bool PrintCanvases = true;   // "false" para no imprimir un archivo pdf

auto c1 = new TCanvas("c1","Dos Fuentes Radioactivas",700,700);
if (PrintCanvases) c1->Print("Out_Generate_2D.pdf["); // open output file

bool stop_at_plot = true;  // Hacer pausa despues de cada plot

int MinuitStatus=0;           // signal of problems
gErrorIgnoreLevel = kWarning; // Dont print "Current canvas added to pdf file"

// Estudio de coberturas: (p3,p4,p34) (CR,LR) (1s,2s,3s)
int ntot=0;
int p3CR1s=0,p3LR1s=0,p3CR2s=0,p3LR2s=0,p3CR3s=0,p3LR3s=0;
int p4CR1s=0,p4LR1s=0,p4CR2s=0,p4LR2s=0,p4CR3s=0,p4LR3s=0;

int InElipse1=0,InElipse2=0,InElipse3=0;
int InMinos1 =0,InMinos2 =0,InMinos3 =0;

/////////// LOOP SOBRE PSEUDO-EXPERIMENTOS ////////////////////

for (int toy = 0; toy < Ntoys; ++toy)  {

fTrue->SetParameters(10,1000,100,40,200); 

TRandom3 r(0); 
int Nevents = r.Poisson(4000);
for (int i = 0; i < Nevents; ++i)  Hdata->Fill(fTrue->GetRandom());               

gStyle->SetOptStat(0);
//Hdata->Draw("E");
c1->SetLogy();
c1->SetGrid();

//Hdata->Fit(fTrue,"N"); 
TFitResultPtr Ptr = Hdata->Fit(fTrue,"SN");
TMatrixDSym cov = Ptr->GetCovarianceMatrix();

// Fitted values of lifetimes
double p3_hat=fTrue->GetParameter(3);
double p4_hat=fTrue->GetParameter(4);

// Donde guarda los errores del fit 
double errorup1,errordn1,errparabolic1,corr1;
double errorup2,errordn2,errparabolic2,corr2;
double errorup3,errordn3,errparabolic3,corr3;

// Ajuste y obtencion de errores para 1,2,3 sigmas

gMinuit->SetErrorDef(2.296); // 1-sigma     2.296
TGraph *gr1b = (TGraph*)gMinuit->Contour(80,3,4); 
gMinuit->mnerrs(4,errorup1,errordn1,errparabolic1,corr1);
/* 
gMinuit->SetErrorDef(6.180);   // 2-sigma    6.180
TGraph *gr2b = (TGraph*)gMinuit->Contour(80,3,4);
gMinuit->mnerrs(4,errorup2,errordn2,errparabolic2,corr2);

gMinuit->SetErrorDef(11.829);   // 3-sigma    11.829
TGraph *gr3b = (TGraph*)gMinuit->Contour(80,3,4); 
gMinuit->mnerrs(4,errorup3,errordn3,errparabolic3,corr3);
*/

/*
  Codigo de arriba para 2 y 3 sigmas remplazado, pues hay que saltear casos 
  que el fit no converge con este mensaje de error:
    Warning in <TM::Contour>: Cannot find more than 4 points, no TGraph returned
    Warning in <TMinuit::Contour>: Returning a TGraph with 4 points only
*/

// Ajuste y obtencion de errores para 2,3 sigmas, salteando si fit no OK

gMinuit->SetErrorDef(6.180); // 2-sigma  6.180
TGraph *gr2b = (TGraph*)gMinuit->Contour(80,3,4); 
MinuitStatus = gMinuit->GetStatus();
if (MinuitStatus!=0) cout<<"MinuitStatus "<<MinuitStatus<<endl<<endl;
if (MinuitStatus!=0) continue;
if (gr2b->GetN() < 40 ) continue;
gMinuit->mnerrs(4,errorup2,errordn2,errparabolic2,corr2);
//printf("Tau4 95.4 CL: %5.1f+/-%5.1f  [%5.1f,%+6.1f]  CI: (%5.1f, %5.1f)\n",
//    p4_hat,errparabolic2,errordn2,errorup2,p4_hat+errordn2,p4_hat+errorup2);

gMinuit->SetErrorDef(11.829); // 3-sigma 11.829
TGraph *gr3b = (TGraph*)gMinuit->Contour(80,3,4); 
MinuitStatus = gMinuit->GetStatus();
if (MinuitStatus!=0) cout<<"MinuitStatus "<<MinuitStatus<<endl<<endl;
if (MinuitStatus!=0) continue;
if (gr3b->GetN() < 40 ) continue;
gMinuit->mnerrs(4,errorup3,errordn3,errparabolic3,corr3);
//printf("Tau4 99.7 CL: %5.1f+/-%5.1f  [%5.1f,%+6.1f]  CI: (%5.1f, %5.1f)\n",
//    p4_hat,errparabolic3,errordn3,errorup3,p4_hat+errordn3,p4_hat+errorup3);

    cout<<endl;

///////// DIBUJA CONTORNOS DE MINOS //////////////
   
   c1->SetLogy(0);
   c1->DrawFrame(30,0,50,800,"68.3% 95.4% 99.7% CL regions for tau3 & tau4;tau3;tau4");

   gr3b->SetFillColorAlpha(590, 0.5);
   //gr3b->SetFillColorAlpha(45, 0.5);
   gr3b->Draw("lf"); // Draw "contour" for 3-sigma
   gr2b->SetFillColorAlpha(406, 0.5);
   //gr2b->SetFillColorAlpha(42, 0.5);
   gr2b->Draw("Lf");  // Draw "contour" for 2-sigma
   gr1b->SetFillColorAlpha(390, 0.5);
   gr1b->Draw("Lf");  // Draw "contour" for 1-sigma

/////////// DIBUJA ELIPSES ////////////////

   double cov11=sqrt(cov(3,3));
   double cov22=sqrt(cov(4,4));
   double cor12=cov(3,4)/sqrt(cov(3,3)*cov(4,4));
   
   RooEllipse* ell1 = new RooEllipse("Parab CI", p3_hat, p4_hat, 1*cov11, 1*cov22, cor12, 100);
   RooEllipse* ell2 = new RooEllipse("Parab CI", p3_hat, p4_hat, 2*cov11, 2*cov22, cor12, 100);
   RooEllipse* ell3 = new RooEllipse("Parab CI", p3_hat, p4_hat, 3*cov11, 3*cov22, cor12, 100);
   
   ell1->SetLineColor(2);
   ell2->SetLineColor(2);
   ell3->SetLineColor(2);
   ell1->SetLineWidth(2);
   ell2->SetLineWidth(2);
   ell3->SetLineWidth(2);
   ell1->SetLineStyle(1);
   ell2->SetLineStyle(1);
   ell3->SetLineStyle(1);
   ell1->SetFillColor(6);
   ell1->Draw("lsame");
   ell2->Draw("lsame");
   ell3->Draw("lsame");

/////////// DIBUJA EL PUNTO FITEADO Y EL VERDADRO ////////////////
   
   TMarker *FittedResult = new TMarker(p3_hat,p4_hat,20);
   FittedResult->SetMarkerColor(kRed);
   FittedResult->SetMarkerSize(1.);
   FittedResult->Draw();

   TMarker *TrueValue = new TMarker(40,200,20);
   TrueValue->SetMarkerColor(kBlue);
   TrueValue->SetMarkerStyle(71);
   TrueValue->SetMarkerSize(2);
   TrueValue->Draw();

   if (PrintCanvases) c1->Print("Out_Generate_2D.pdf"); 
   c1->Update();

   ntot++;

// Check if inside Minos Contours and inside Ellipses
// (uncomment next 6 lines to print which cases are inside)

// if (TMath::IsInside(40., 200.,  80, gr1b->GetX(), gr1b->GetY())) cout<<"Minos1"<<endl;
// if (TMath::IsInside(40., 200.,  80, gr2b->GetX(), gr2b->GetY())) cout<<"Minos2"<<endl;
// if (TMath::IsInside(40., 200.,  80, gr3b->GetX(), gr3b->GetY())) cout<<"Minos3"<<endl;
// if (TMath::IsInside(40., 200., 101, ell1->GetX(), ell1->GetY())) cout<<"Elipse1"<<endl;
// if (TMath::IsInside(40., 200., 101, ell2->GetX(), ell2->GetY())) cout<<"Elipse2"<<endl;
// if (TMath::IsInside(40., 200., 101, ell3->GetX(), ell3->GetY())) cout<<"Elipse3"<<endl;

if (TMath::IsInside(40., 200.,  80, gr1b->GetX(), gr1b->GetY())) InMinos1++;
if (TMath::IsInside(40., 200.,  80, gr2b->GetX(), gr2b->GetY())) InMinos2++;
if (TMath::IsInside(40., 200.,  80, gr3b->GetX(), gr3b->GetY())) InMinos3++;
if (TMath::IsInside(40., 200., 101, ell1->GetX(), ell1->GetY())) InElipse1++;
if (TMath::IsInside(40., 200., 101, ell2->GetX(), ell2->GetY())) InElipse2++;
if (TMath::IsInside(40., 200., 101, ell3->GetX(), ell3->GetY())) InElipse3++;

if( p4_hat-errparabolic1 < 200 && p4_hat+errparabolic1 > 200) p4CR1s++;
if( p4_hat+errordn1      < 200 && p4_hat+errorup1      > 200) p4LR1s++;
if( p4_hat-errparabolic2 < 200 && p4_hat+errparabolic2 > 200) p4CR2s++;
if( p4_hat+errordn2      < 200 && p4_hat+errorup2      > 200) p4LR2s++;
if( p4_hat-errparabolic3 < 200 && p4_hat+errparabolic3 > 200) p4CR3s++;
if( p4_hat+errordn3      < 200 && p4_hat+errorup3      > 200) p4LR3s++;

if(stop_at_plot) {
 int c = getchar();   // pauses: hit <RET> to continue     
 if (c == '.') stop_at_plot = false;  // Hit "." to stop drawing
}

Hdata->Reset();
delete gr1b;
delete gr2b;
delete gr3b;
delete FittedResult;
delete TrueValue;

}  // Fin de los Toys

// Cobertura p4 CR (Crame-Rao) y LR (Likelihood-Ratio) error 1-2-3 sigmas

float CLp4CR1s = 100.*p4CR1s/ntot;
float CLp4LR1s = 100.*p4LR1s/ntot;
float CLp4CR2s = 100.*p4CR2s/ntot;
float CLp4LR2s = 100.*p4LR2s/ntot;
float CLp4CR3s = 100.*p4CR3s/ntot;
float CLp4LR3s = 100.*p4LR3s/ntot;

// Cobertura Minos y Elipse,  1-2-3 sigmas
float CL_InMinos1  = 100.*InMinos1 /ntot;
float CL_InMinos2  = 100.*InMinos2 /ntot;
float CL_InMinos3  = 100.*InMinos3 /ntot;
float CL_InElipse1 = 100.*InElipse1/ntot;
float CL_InElipse2 = 100.*InElipse2/ntot;
float CL_InElipse3 = 100.*InElipse3/ntot;

//////////// PRINT SUMMARY STATISTICS ////////////

cout.precision(2);      // Print using four decimals precision
cout<<fixed;            // ...exactly four decimals

cout<<"Ntot: "<<ntot<<endl<<endl;

cout<<"Parametro 4: Contador CR 1s,2s,3s:  ";
cout<<  p4CR1s<<"  "<<  p4CR2s<<"  "<< p4CR3s<<endl;
cout<<"Parametro 4: Contador LR 1s,2s,3s:  ";
cout<<  p4LR1s<<"  "<<  p4LR2s<<"  "<< p4LR3s<<endl<<endl;

cout<<"Parametro 4: Coverage CR 1s,2s,3s:  ";
cout<<CLp4CR1s<<"  "<<CLp4CR2s<<"  "<<CLp4CR3s<<endl;
cout<<"Parametro 4: Coverage LR 1s,2s,3s:  ";
cout<<CLp4LR1s<<"  "<<CLp4LR2s<<"  "<<CLp4LR3s<<endl<<endl;

cout<<"Contador MINOS 2D para 1s,2s,3s:   ";
cout<<   InMinos1<<" "<<    InMinos2<<" "<<    InMinos3<<endl;
cout<<"Coverage MINOS 2D para 1s,2s,3s:   ";
cout<<CL_InMinos1<<"  "<<CL_InMinos2<<"  "<<CL_InMinos3<<" "<<endl<<endl;

cout<<"Contador Elipse 2D para 1s,2s,3s:  ";
cout<<    InElipse1<<" "<<    InElipse2<<" "<<    InElipse3<<endl;
cout<<"Coverage Elipse 2D para 1s,2s,3s:  ";
cout<< CL_InElipse1<<"  "<<CL_InElipse2<<"  "<<CL_InElipse3<<" "<<endl<<endl;

if (PrintCanvases) c1->Print("Out_Generate_2D.pdf]");  // close file
delete c1;
return;
}