// The function I want to minimize
double myFunc(double x, double y) {
    return x*x*x*x + x*x + y*y + x*y;
}    

// A wrapper for Minuit around my function.
//    nDim -- The number of dimensions
//    gout -- The calculated gradient (DON'T CALCULATE!)
//    result -- The function value
//    point -- The point to calculate the function at
//    iflag -- A flag: 1) initialize, 2) calculate 3) end
// The function value should always be calculated, so you can usually ignore
// the value of iflag.  
void fcn(int& nDim, double* gout, double& result, double* point, int iflag) {
    result = myFunc(point[0],point[1]);
}

void minimizerExample() {
    // Create a new minimizer.  The argument is the number dimensions.
    TFitter* minimizer = new TFitter(2);
    { 
        double p1 = -1;
        minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);
    }
    minimizer->SetFCN(fcn);
    minimizer->SetParameter(0,"X",10,1,0,0);
    minimizer->SetParameter(1,"Y",10,1,0,0);

    // Run the simplex minimizer to get close to the minimum
    minimizer->ExecuteCommand("SIMPLEX",0,0);

    // Run the migrad minimizer (an extended Powell's method) to improve the
    // fit.
    minimizer->ExecuteCommand("MIGRAD",0,0);

    double bestX = minimizer->GetParameter(0);
    double bestY = minimizer->GetParameter(1);
    double minimum = myFunc(bestX, bestY);
    
    // Scan the X parameter to find it's uncertainty
    double errX;
    for (errX=0; errX<10; errX+=0.001) {
        minimizer->SetParameter(0,"X",errX,0,0,0);
        minimizer->SetParameter(1,"Y",1,1,0,0);
        minimizer->ExecuteCommand("MIGRAD",0,0);
        double m = myFunc(minimizer->GetParameter(0),
                          minimizer->GetParameter(1));
        if (m-minimum > 1.0) break;
    }

    // Scan the Y parameter to find it's uncertainty
    double errY;
    for (errY=0; errY<10; errY+=0.001) {
        minimizer->SetParameter(0,"X",1,1,0,0);
        minimizer->SetParameter(1,"Y",errY,0,0,0);
        minimizer->ExecuteCommand("MIGRAD",0,0);
        double m = myFunc(minimizer->GetParameter(0),
                          minimizer->GetParameter(1));
        if (m-minimum > 1.0) break;
    }

    std::cout << bestX << "+-" << errX << std::endl;
    std::cout << bestY << "+-" << errY << std::endl;
    std::cout << minimum << std::endl;

}