Double_t CalculateNLL(Double_t *par)
{
    ROOT::Math::KahanSum<double> fval(0.0);
    const double mu_min = 1e-12;


    //sum_mu = total number of events predicted model
    double sum_mu = ComputeSumMu(par);

    //sum_N = total number of observed events
    double sum_N  = hist->Integral();

    // rescales model so that the total number of predicted events exactly matches the number of observed events.
    double scale  = sum_N / sum_mu;


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

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

            int N = hist->GetBinContent(ix, iy);

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

            double pdf =
                (1.5 * B_f * M_PI * beta_val * alpha * alpha *
                 Function_Theta(x, y, par)) / s;

             // Luminosity==1e3 
            double mu = scale * pdf * dx * dy * 1e3 * 3.87e8;
            mu = std::max(mu, mu_min);

            fval += mu - N * std::log(mu);
        }
    }
    return fval.Sum();
}




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



void Fit_NLL()
{
    PrepareHistogram();


  ROOT::Math::Minimizer* minimizer =
        ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");

    minimizer->SetMaxFunctionCalls(100000);
    minimizer->SetMaxIterations(100000);
    minimizer->SetTolerance(1e-3);
    minimizer->SetPrintLevel(2);

    ROOT::Math::Functor f(&NLL_wrapper, 3);
    minimizer->SetFunction(f);
    minimizer->SetStrategy(1);


    minimizer->SetLimitedVariable(0, "deltaAy", 1e-6, 1e-4, 1e-1, 1e-4);
    minimizer->SetLimitedVariable(1, "deltaAz", 0.0, 1e-3, -1e-2, 1e-2);
    minimizer->SetLimitedVariable(2, "deltaBz", 0.0, 1e-3, -1e-2, 1e-2);


   
    minimizer->Minimize();
//
    const double *xs = minimizer->X();
    const double *errs = minimizer->Errors();

    std::cout << "Fit results:\n";
    std::cout << "deltaAy = " << xs[0] << " +/- " << errs[0] << std::endl;
    std::cout << "deltaAz = " << xs[1] << " +/- " << errs[1] << std::endl;
    std::cout << "deltaBz = " << xs[2] << "+/- " << errs[2] << std::endl;
}