Minuit fit values of the parameters are way far from the initialized values

ROOT
1

_ROOT Version: 6.14/06
_Platform: MacOS X Mojave
Compiler: Not Provided

Hi,
I am trying to use Minuit to fit three parameters from a chi-square distribution. But the results from Minuit don’t make sense at all and are really off from their initialized values.
This is what my code looks like:

#include <iostream>
#include <fstream>
#include <vector>
#include <iomanip>
#include <cmath>
#include "beta_decay.h"

#include <TMinuit.h>
#include <TFile.h>
#include <TMinuit.h>
#include <TApplication.h>
#include <TROOT.h>

#define N_BINS 500

using namespace std;
static Params p = Params();
static DecayFunction pol = DecayFunction(p);
std::vector<double> a;
std::vector<double> da;
std::vector<double> A;
std::vector<double> dA;
std::vector<double> x_avg;

double A_e(double x, double *par)
{
    double A;
    double beta, eps, R, E0, Ee;
    double kappa_p = 2.792847356 - 1.;
    double kappa_n = -1.91304272;
    double f2, g2, f3;
    double bBSM;
    double me = 0.510998928 * 1000.; // [keV]
    double mn = 939.565379 * 1000.;  // [keV]
    double mp = 938.272046 * 1000.;  // [keV]
    double costhAvg = 0.5;
    double g_s = 1.;
    double g_t = 1.;

    f2 = (kappa_p - kappa_n) / 2.;
    f3 = 0.;
    g2 = 0;
    bBSM = 0.;

    E0 = (mn - mp) - ((mn - mp) * (mn - mp) - me * me) / (2. * mn);
    // beta = sqrt(T * T + 2. * me * T) / (T + me);
    eps = (me * me) / (mn * mn);
     R = E0 / mn;<a class="attachment" href="/uploads/default/original/3X/9/b/9b4cd401d9e0e80751b22d0b9d964c2b19bb4251.txt">output.txt</a> (5.9 KB)
 

    A = (2. * par[0] * (1 - par[0]) + 2 * 16 * g_t * g_t * par[2] * par[2] + 8 * g_s * par[1] * g_t * par[2]) /
        (1. + 3. * par[0] * par[0] + g_s * g_s * par[1] * par[1] + 3 * 16 * g_t * g_t * par[2] * par[2]);

    A += (eps / (R * x) * (4. * par[0] * par[0] * (1. - par[0]) * (1. + par[0] + 2. * f2) + 4. * par[0] * (1. - par[0]) * (par[0] * g2 - f3)) +
          R * ((2. / 3.) * (1. + par[0] + 2. * (f2 + g2)) *
               (3. * par[0] * par[0] + 2. * par[0] - 1.)) +
          (R * x) * ((2. / 3.) * (1. + par[0] + 2. * f2) * (1. - 5. * par[0] - 9. * par[0] * par[0] - 3. * par[0] * par[0] * par[0]) +
                     (4. / 3.) * g2 * (1. + par[0] + 3. * par[0] * par[0] + 3. * par[0] * par[0] * par[0]))) /
         ((1. + 3. * par[0] * par[0]) * (1. + 3. * par[0] * par[0]));

    return A;
}

// function for chi2
void chi2Function(int &npar, double *gin, double &f, double *par, int iflag)
{
    double delta = 0.;
    double chi_square = 0.;

    for (auto i = 0; i < (N_BINS - 5); ++i)
    {
        delta = (A[i] - A_e(x_avg[i], par)) / dA[i];
        chi_square += delta * delta;
    }
    f = chi_square;
}

// main function
int main()
{
    double dalpha_e_hold;
    double alpha_e_hold;
    double x_avg_hold;
    double A_hold;
    double dA_hold;
    double beta_hold;

    ifstream fin1("electron_asymm.dat");
    
    while (fin1 >> alpha_e_hold >> dalpha_e_hold >> beta_hold >> x_avg_hold >> A_hold >> dA_hold)
    {
        x_avg.push_back(x_avg_hold);
        A.push_back(A_hold);
        dA.push_back(dA_hold);
    }

    TMinuit minuit(3);
    minuit.SetFCN(chi2Function);

    double arglist[10];
    int ierflg = 0;
    arglist[0] = 3.67; // Value of UP
    minuit.mnexcm("SET ERR", arglist, 1, ierflg);

    arglist[0] = 5000;
    arglist[1] = 1.;

    minuit.mnparm(0, "lambda", 1.2701, 0.00001, 0.0, 0.0, ierflg);
    minuit.mnparm(1, "epsilon_s", 0.0, 0.0001, 0.0, 0.0, ierflg);
    minuit.mnparm(2, "epsilon_t", 0.0, 0.0001, 0.0, 0.0, ierflg);

    minuit.mnexcm("MIGRAD", arglist, 2, ierflg);

    minuit.mnexcm("MINOS", arglist, 2, ierflg);

 }

I am attaching the output I got when I run this code in the terminal.output.txt (5.9 KB)

2

I think @moneta can help you with this question.

3

HI,
You should carefully check your function to minimize. I see the derivatives of epsilon_s and epsilon_t are zero. Is this expected ? Maybe also the initial parameter point is not the best one.

Lorenzo

4

@moneta Thanks for your reply,
I tried running again my code and like you suggested I changed the initial parameter values for \epsilon_s and \epsilon_t but this time I removed minuit.mnexcm("MINOS",...) command from the code and it converged with MIGRAD.
I don’t think the derivatives of epsilon_s and epsilon_t should be zero but I think my minimizer is working fine because for the first parameter I got the fitted value within one \sigma standard deviation. But the errors on \epsilon_s and \epsilon_t are way too big.
One last thing is that I still have problem with the error matrix calculations with HESSE.
This is my improved code looks like:

#include <iostream>
#include <fstream>
#include <vector>
#include <iomanip>
#include <cmath>
#include "beta_decay.h"

#include <TMinuit.h>
#include <TFile.h>
#include <TMinuit.h>
#include <TApplication.h>
#include <TROOT.h>

#define N_BINS 500

using namespace std;
// static Params p = Params();
// static DecayFunction pol = DecayFunction(p);
std::vector<double> a;
std::vector<double> da;
std::vector<double> A;
std::vector<double> dA;
std::vector<double> x_avg;

double A_e(double x, double *par)
{
    double A;
    double beta, eps, R, E0, Ee;
    double kappa_p = 2.792847356 - 1.;
    double kappa_n = -1.91304272;
    double f2, g2, f3;
    double bBSM;
    double me = 0.510998928 * 1000.; // [keV]
    double mn = 939.565379 * 1000.;  // [keV]
    double mp = 938.272046 * 1000.;  // [keV]
    double costhAvg = 0.5;
    double g_s = 1.;
    double g_t = 1.;

    f2 = (kappa_p - kappa_n) / 2.;
    f3 = 0.;
    g2 = 0;
    bBSM = 0.;

    E0 = (mn - mp) - ((mn - mp) * (mn - mp) - me * me) / (2. * mn);
    // beta = sqrt(T * T + 2. * me * T) / (T + me);
    eps = (me * me) / (mn * mn);
    // x = (T + me) / E0;
    R = E0 / mn;
    // Ee = T + me;

    A = (2. * par[0] * (1 - par[0]) + 2 * 16 * g_t * g_t * par[2] * par[2] + 8 * g_s * par[1] * g_t * par[2]) /
        (1. + 3. * par[0] * par[0] + g_s * g_s * par[1] * par[1] + 3 * 16 * g_t * g_t * par[2] * par[2]);

    A += (eps / (R * x) * (4. * par[0] * par[0] * (1. - par[0]) * (1. + par[0] + 2. * f2) + 4. * par[0] * (1. - par[0]) * (par[0] * g2 - f3)) +
          R * ((2. / 3.) * (1. + par[0] + 2. * (f2 + g2)) *
               (3. * par[0] * par[0] + 2. * par[0] - 1.)) +
          (R * x) * ((2. / 3.) * (1. + par[0] + 2. * f2) * (1. - 5. * par[0] - 9. * par[0] * par[0] - 3. * par[0] * par[0] * par[0]) +
                     (4. / 3.) * g2 * (1. + par[0] + 3. * par[0] * par[0] + 3. * par[0] * par[0] * par[0]))) /
         ((1. + 3. * par[0] * par[0]) * (1. + 3. * par[0] * par[0]));

    // A = A / (1. + bBSM * me / Ee);

    return A;
}

// function for chi2
void chi2Function(int &npar, double *gin, double &f, double *par, int iflag)
{
    double delta = 0.;
    double chi_square = 0.;

    for (auto i = 5; i < (x_avg.size() - 5); ++i)
    {
        delta = (A[i] - A_e(x_avg[i], par)) / dA[i];
        chi_square += delta * delta;
    }
    std::cout << chi_square << "    " << *par << "    " << *(par + 1)
              << "    " << *(par + 2) << '\n';
    f = chi_square;
}

// main function
int main()
{
    double dalpha_e_hold;
    double alpha_e_hold;
    double x_avg_hold;
    double A_hold;
    double dA_hold;
    double beta_hold;

    ifstream fin1("electron_asymm.dat");

    while (fin1 >> alpha_e_hold >> dalpha_e_hold >> beta_hold >> x_avg_hold >> A_hold >> dA_hold)
    {
        x_avg.push_back(x_avg_hold);
        A.push_back(A_hold);
        dA.push_back(dA_hold);
    }

    TMinuit minuit(3);
    minuit.SetFCN(chi2Function);

    double arglist[10];
    int ierflg = 0;
    arglist[0] = 3.67; // Value of UP
    minuit.mnexcm("SET ERR", arglist, 1, ierflg);

    arglist[0] = 1000;
    arglist[1] = 1;

    minuit.mnparm(0, "lambda", 1.2723, 0.00001, 0.0, 0.0, ierflg);
    minuit.mnparm(1, "epsilon_s", 0.000001, 0.0000001, 0.0, 0.0, ierflg);
    minuit.mnparm(2, "epsilon_t", 0.000001, 0.0000001, 0.0, 0.0, ierflg);

    minuit.mnexcm("MIGRAD", arglist, 2, ierflg);

    // minuit.mnexcm("MINOS", arglist, 2, ierflg);
  
}

Here is my output from the terminal:output2.txt (8.1 KB)

5

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.