Double_t CalculateNLL(Double_t *par)
{
    ROOT::Math::KahanSum<double> fval(0.0);

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


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

            Int_t MC_N = hist->GetBinContent(ix, iy);

            Double_t x_f_2      = 0.04*ix - 0.02;
            Double_t costheta_2 = 0.04*iy - 1.02;

            // PDF is in GeV^-2 we multiply by 3.87e8 to make it in pb 
            Double_t pdf =
                (1.5 * B_f * M_PI * beta_val * alpha * alpha *
                 Function_Theta(x_f_2, costheta_2, par)) / s;

            // Luminosity = 1000 ab^-1 ---> converted to pb-1 : 1e9
            Double_t mu = pdf * 0.04 * 0.04 * 1e9 * 3.87e8;

            if (mu > 0 && MC_N > 0) {
                fval += mu - MC_N * std::log(mu);
                 // can be used or not
                fval +=logFactorial(MC_N);
            }




        }
    }
    return fval.Sum();
}


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


void Fit_chis_WithoutScan_Delta_AY_Reco()
{
    PrepareHistogram();

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

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

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


    minimizer->SetVariable(0, "deltaAy", 0.0, 1e-5);
    minimizer->SetVariable(1, "deltaAz", 0.0, 1e-5);
    minimizer->SetVariable(2, "deltaBz", 0.0, 1e-5);


    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;
}