//Limits using the asymptotic formulae from the corresponding paper

TH1F *h_data = new TH1F();
TH1F *h_signal = new TH1F();
TH1F *h_bkg = new TH1F();

double log_likelihood(TH1F *h_data, TH1F *h_signal, TH1F *h_bkg, int i,double x,int flag)   
{
   double val=0;                                                                     
   TF1 *f= new TF1("f","TMath::Poisson(x,[0])",0,1e5);

   double s=h_signal->GetBinContent(i);
   double b=h_bkg->GetBinContent(i);

   f->SetParameter(0,((x*s)+(b))); 

   if(flag==0) //asimov_data
   {
      val=TMath::Log(f->Eval(h_bkg->GetBinContent(i)));
   }

   if(flag==1) //observed_data
   {
      val=TMath::Log(f->Eval(h_data->GetBinContent(i)));
   }
   

   return val;

}

double likelihood(const double *par)
{
   double lnL=0;
   double value=0;

   TF1 *f1= new TF1("f1","TMath::Poisson(x,[0])",0,1e5); 
   
    for(int i=1;i<=h_data->GetNbinsX();i++)
    {
       double s=h_signal->GetBinContent(i);
       double b=h_bkg->GetBinContent(i);

       f1->SetParameter(0,(par[0]*s)+b);

       value=TMath::Log((f1->Eval(h_data->GetBinContent(i))));

       lnL=lnL+value;
    }

    lnL=-2*lnL;

    return lnL;

 }

 double likelihood_asimov(const double *par)
{
   double lnL=0;
   double value=0;

   TF1 *f1= new TF1("f1","TMath::Poisson(x,[0])",0,1e5); 
   
    for(int i=1;i<=h_data->GetNbinsX();i++)
    {
       double s=h_signal->GetBinContent(i);
       double b=h_bkg->GetBinContent(i);

       f1->SetParameter(0,(par[0]*s)+b);

       value=TMath::Log((f1->Eval(h_bkg->GetBinContent(i))));

       lnL=lnL+value;
    }

    lnL=-2*lnL;

    return lnL;

 }


double test_statistic(double x,double mu_hat, int flag)
{
   double val=0;
   
   if(mu_hat<=x)
   {
      for(int i=1;i<=h_data->GetNbinsX();i++)
      {
         val=val+(log_likelihood(h_data,h_signal,h_bkg,i,x,flag)-(log_likelihood(h_data,h_signal,h_bkg,i,mu_hat,flag)));
      }

         
      val=-2*val;
   }

   if(mu_hat>x)
   {
      val=0;
   }


   return val;
 }

 double cumulative(double test_mu,double mu_hat,double mu_prime,double sigma,int flag)
 {
    double val=0;

    double arg1 = sqrt(test_statistic(test_mu,mu_hat,flag));

    double arg2 = (test_mu-mu_prime)/(sigma);

    val = ROOT::Math::normal_cdf(arg1-arg2);

    return val;
 }


void asymptotics()
{
   TFile *file = new TFile("drell_yan_sr_histos.root","read");
   
   h_data=(TH1F*)file->Get("data_DY_SR");
   h_signal=(TH1F*)file->Get("hvt_m500_302266_DY_SR");
   h_bkg=(TH1F*)file->Get("Asimov_DY_SR");

   double test_mu = 0.00684735;

   ROOT::Math::Minimizer *min_asimov = new ROOT::Minuit2::Minuit2Minimizer();

   min_asimov->SetMaxFunctionCalls(10000);
   min_asimov->SetTolerance(0.001);
   min_asimov->SetPrintLevel(1);

   ROOT::Math::Functor f_asimov(&likelihood_asimov,1);

   double variable_asimov[1] = {0};
   double step_asimov[1] = {0.01};

   min_asimov->SetFunction(f_asimov);

   min_asimov->SetVariable(0,"Mu",variable_asimov[0], step_asimov[0]);

   min_asimov->Minimize();

   const double *xs_asimov= min_asimov->X();
   const double *x_asimov = min_asimov->Errors();

   double mu_hat_asimov=xs_asimov[0];

   ROOT::Math::Minimizer *min = new ROOT::Minuit2::Minuit2Minimizer();

   min->SetMaxFunctionCalls(10000);
   min->SetTolerance(0.001);
   min->SetPrintLevel(1);

   ROOT::Math::Functor f(&likelihood,1);

   double variable[1] = {0};
   double step[1] = {0.01};

   min->SetFunction(f);

   min->SetVariable(0,"Mu",variable[0], step[0]);

   min->Minimize();

   const double *xs= min->X();
   const double *x = min->Errors();

   double mu_hat=xs[0];

   double sigma = (test_mu-mu_hat_asimov)/(sqrt(test_statistic(test_mu,mu_hat_asimov,0)));

   double clsb = 1-(cumulative(test_mu,mu_hat,test_mu,sigma,1));
   double clb = 1-(cumulative(test_mu,mu_hat,0,sigma,1));

   cout << endl << "CLs+b: " << clsb << endl;
   cout << "CLb: " << clb << endl;
   cout << "CLs: " << clsb/clb << endl;

   double clsb_exp = 1-(cumulative(test_mu,mu_hat_asimov,test_mu,sigma,0));
   double clb_exp = 1-(cumulative(test_mu,mu_hat_asimov,0,sigma,0));

   cout << endl << "CLs+b_exp: " << clsb_exp << endl;
   cout << "CLb_exp: " << clb_exp << endl;
   cout << "CLs_exp: " << clsb_exp/clb_exp << endl;

   cout << endl <<"Sigma: " << sigma << endl;

   cout << "Expected Limit-2σ : " << test_mu - (2*sigma) << endl;

   
}