#include "TMath.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"


double RdiDt(double tau, double t1, double t2)
{
  //std::cout<<" son qui "<<std::endl;
  return exp(-t1/tau)-exp(-t2/tau);
}

double invMyfunction(double xmeas, TF1* f)
{
  double step = (350.-150.)/10000;
  double tau;
  double x;
  for (int i=0; i<10000; ++i)
    {
      tau = 150.+float(i)*step;
      x = f->Eval(tau);
      if (x<xmeas)
	{
	  return tau;
	}
    }
  return -9999;
}
// Macro myfunc.C
Double_t myfunction(Double_t *x, Double_t *par)
{
   double tau =x[0];
   double t1=par[0];
   double t2=par[1];
   //std::cout<<" in myfunction "<<tau<<" t1, t2 = "<<t1<<" "<<t2<<std::endl;
   return par[2]*RdiDt(tau,t1,t2);
}
Double_t myfunctionMin(Double_t *x, Double_t *par)
{
   double tau =x[0];
   double t1=par[0];
   double t2=par[1];
   //std::cout<<" in myfunction "<<tau<<" t1, t2 = "<<t1<<" "<<t2<<std::endl;
   double xx = par[2]*RdiDt(tau,t1,t2);
   return xx-par[3]*sqrt(xx);
}
Double_t myfunctionMax(Double_t *x, Double_t *par)
{
   double tau =x[0];
   double t1=par[0];
   double t2=par[1];
   //std::cout<<" in myfunction "<<tau<<" t1, t2 = "<<t1<<" "<<t2<<std::endl;
   double xx = par[2]*RdiDt(tau,t1,t2);
   return xx+par[3]*sqrt(xx);//par[2]*RdiDt(tau,t1,t2)+par[3]*sqrt(par[2]*RdiDt(tau,t1,t2));
}
void Neyman_esperimento()
{
  TRandom* gen = new TRandom(4564567);
  TH2F* hbox = new TH2F("hbox","hbox",100,120.,320.,100,1500.,3000.); 

  double tauvero = 220.;
  double t1 = 10.;
  double t2 = 60.;

  double RdiDtvero = RdiDt(tauvero, t1, t2);// exp(-t1/tauvero)-exp(-t2/tauvero);
  std::cout<<"tau = "<<tauvero<<" Frazioni di eventi tra t1="<<t1<<" e t2="<<t2<<" -->>> "<<RdiDtvero<<std::endl;

  TF1* f1    = new TF1("NevDt",   myfunction,   50.,1000.,3);
  TF1* f1min = new TF1("NevDtmin",myfunctionMin,50.,1000.,4);
  TF1* f1max = new TF1("NevDtmax",myfunctionMax,50.,1000.,4);
  int N=10000;

  hbox->Draw();
  f1->SetParameter(0,t1);
  f1->SetParameter(1,t2);
  f1->SetParameter(2,N);
  f1->SetParNames("t1", "t2");
  f1->Draw("SAME");
  f1min->SetParameter(0,t1);
  f1min->SetParameter(1,t2);
  f1min->SetParameter(2,N);
  f1min->SetParameter(3,1.645); //90% CL 
  f1max->SetParameter(0,t1);
  f1max->SetParameter(1,t2);
  f1max->SetParameter(2,N);
  f1max->SetParameter(3,1.645); //90% CL
  f1min->SetLineColor(kBlue);
  f1max->SetLineColor(kBlue);
  f1min->Draw("same");
  f1max->Draw("same");
  double yy    = f1->Eval(tauvero);
  double yymin = f1min->Eval(tauvero);
  double yymax = f1max->Eval(tauvero);
  std::cout<<" Eventi attesi tra 10 e 60 s (per tau vero = 220 s) = "<<yy<<" min, max (90% CL) -> "<<yymin<<" "<<yymax<<std::endl;
  
  double NdiDtvero = RdiDtvero * N;

  //return;


  std::cout<<" Starting the simulation"<<std::endl;
  std::cout<<" Cumulativa a -1sigma = ";
  double x1=NdiDtvero-sqrt(NdiDtvero);
  double xx = ROOT::Math::normal_cdf(x1, sqrt(NdiDtvero), NdiDtvero);
  std::cout<<xx<<std::endl;

  //return;
  double xmeas;
  double taumin;
  double taumax;

  xmeas = 1850.;
  taumin = invMyfunction(xmeas, f1min);      
  taumax = invMyfunction(xmeas, f1max);
  std::cout<<" x_misurato = "<<xmeas<<" tau_min, max -> "<<taumin<<" "<<taumax<<std::endl;

  
  int inCLint = 0;
  int NEXP = 10000;
  for (int i=0; i<NEXP; ++i)
    {
      xmeas = gen->Gaus(NdiDtvero, sqrt(NdiDtvero));

      taumin = invMyfunction(xmeas, f1min);      
      taumax = invMyfunction(xmeas, f1max);

      // htaumin -> Fill(taumin);
      // htaumax -> Fill(taumax);

      if (tauvero<taumax && tauvero>taumin) inCLint++;
      if (i%1000==0) std::cout<<" experiment n. "<<i<<" measurement = "<<xmeas<<" taumin/max "<<taumin<<"/"<<taumax<<" tauveroInClInterval "<<  (tauvero<taumax && tauvero>taumin)<<std::endl;
	
    }
  double fracInCLint = double (inCLint)/double(NEXP);
  std::cout<<" NEXP = "<<NEXP<<" fracInCLint = "<<fracInCLint<<std::endl; 
  
}