void PrepareHistogram()

{

  hist = new TH2D("hist", "MC Distribution", 20, 0.12, 1.0, 20, -0.9, 0.9);

  double N_gen = tree->GetEntries();
  double sigma_MG =0.087752;
  double luminosity = 1e6;    // pb
  double weight = luminosity * sigma_MG / N_gen;


for (int i = 0; i < N_gen; i++) {
         tree->GetEntry(i);
         hist->Fill(x_f, -costheta,weight);
}
 hist->SetDirectory(0);
}


Double_t CalculateNLL(const Double_t *par)
{

double deltaAy = par[0];
double deltaAz = par[1];
double deltaBz = par[2];
double deltaBy = par[3];
double kappa   = par[4];


    ROOT::Math::KahanSum<double> fval(0.0);
    const double mu_min = 1e-12;

    int nx = hist->GetNbinsX();
    int ny = hist->GetNbinsY();

    for (int ix = 1; ix <= nx; ix++) {
        for (int iy = 1; iy <= ny; iy++) {

            double x  = hist->GetXaxis()->GetBinCenter(ix);
            double cosTheta = hist->GetYaxis()->GetBinCenter(iy);

            double dx = hist->GetXaxis()->GetBinWidth(ix);
            double dy = hist->GetYaxis()->GetBinWidth(iy);

            double S = Function_Theta(x, cosTheta, par);

            double prefactor =
                (3.0 * M_PI* beta_val * alpha * alpha * B_f) / (2.0 * s);


            double mu = kappa*prefactor*  S * dx * dy * 1e6 *0.387e9;

            mu = std::max(mu, mu_min);
            int N = hist->GetBinContent(ix, iy);

 
            fval += mu - N * std::log(mu);

        }
    }

    return fval.Sum();
}




#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"

double NLL_wrapper(const double *par)
{
     return CalculateNLL(const_cast<double*>(par));

}


void Likelihood_From_Old() {

        PrepareHistogram();

        double par_SM[4] = {0,0,0,1.0};

        std::cout << "Histogram integral = " << hist->Integral() << std::endl;
        std::cout << "SM analytic total (kappa=1) = " << ComputeSumMu(par_SM) * 1e6 * 0.387e9 << std::endl;



        auto minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Minimize");


     ROOT::Math::Functor f(&NLL_wrapper, 5);
     minimizer->SetFunction(f);


 
  
    minimizer->SetMaxFunctionCalls(50000);
    minimizer->SetMaxIterations(50000);
    minimizer->SetTolerance(1e-6);
    minimizer->SetPrintLevel(2);
    minimizer->SetStrategy(2);


  
   minimizer->SetVariable(0, "deltaAy", 0.000000, 0.000001);
   minimizer->SetVariable(1, "deltaAz", 0.000000, 0.000001);
   minimizer->SetVariable(2, "deltaBz", 0.000000, 0.000001);
   minimizer->SetFixedVariable(3, "deltaBy", 0.00000);
   //minimizer->SetFixedVariable(4, "kappa", 1.0);
   minimizer->SetVariable(4, "kappa",1.0, 0.01);

   minimizer->Minimize();


   double NLL_min = minimizer->MinValue();

   const double* best = minimizer->X();
   const double* err  = minimizer->Errors();


}