Hello everyone,
I am having some issues understanding the output of a chi squared minimization that I am performing with Minuit2. I have the stand alone version of Minuit2 6.15.
Below is my code. What I wanted to do was to have the mass(int i, double p0, double p1, double p2, double p3){} provide the relevant eigenvalues to be used in the chi squared minimization. It should be in total seven eigenvalues. But when I echo print the output of mass() to the screen it actually keeps printing numbers.
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FCNBase.h"
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
using namespace ROOT::Minuit2;
using namespace Eigen;
const double PI = 4*atan(1);
// Define minimizer class (inherits from FCNBase)
class Minimizer : public FCNBase {
private:
const std::vector< std::vector<double> > &m_data;
// const std::vector<double> &m_mass;
public:
Minimizer(const std::vector< std::vector<double> > &data) : m_data(data) {}
~Minimizer() {}
// functions to find factorial and double factorial of given number - used in the Vqq matrix element
//Omitted here (but they are in the actual code)
double diagonalization(double n, double L, double p0, double p1, double p2, double p3){
// Declare the number of terms in the expansion and the number of iterations for the numerical integrals
int N1 = 12;
// Declare the variational parameters - values taken from the hybrid decay paper: https://arxiv.org/abs/2005.10850
double alphaInit = 0.005;
// Define the matrices - the dimension can be set to be any of the above indexes since they are generally the same
MatrixXd HMX(N1, N1);
MatrixXd NMX(N1, N1);
MatrixXd KMX(N1, N1);
MatrixXd VMX(N1, N1);
//Declare the variables for the minimum eigenvalue
double minEvalue0 = 0.0; // ground state eigenvalue
double minEvalue1 = 0.0; // first excited state eigenvalue
double minEvalue2 = 0.0; // second excited state eigenvalue
double minEvalue3 = 0.0; // third excited state eigenvalue
std::vector<double> eValues0;
std::vector<double> eValues1;
std::vector<double> eValues2;
std::vector<double> eValues3;
std::vector<double> minEvalues;
// Loop over the variational parameters and the matrix elements (k1=n', l1=m', k2=n, l2=m).
for (double alpha = alphaInit; alpha <= 1.000; alpha += 0.005){
for (int k1 = 1; k1 <= N1; k1++){
for (int k2 = 1; k2 <= N1; k2++){
// Define the full Hamiltonian in configuration space - works for all values of L
// Overlap Matrix
NMX(k1-1,k2-1) = pow(2.0,L + 1)*sqrt(2.0)*pow(k1*k2,0.75 + L/2.0)/(pow(k1+k2,1.5 + L));
// Mass and kinetic terms
KMX(k1-1,k2-1) = (pow(2.0,L + 1)*sqrt(2.0)*pow(k1*k2,0.75 + L/2.0)/(pow(k1+k2,1.5 + L))) * (2.0*p0 + (3.0+2.0*L)*(k1*k2)*alpha*alpha/(p0*(k1+k2)));
// Potential term
VMX(k1-1,k2-1) = (powerfactorial(L)*pow(k1*k2,0.75 + L/2.0)/(doublefactorial(2*L+1)*sqrt(PI)*alpha*pow(k1+k2,2.0+L)))*(2*p1*alpha*alpha*(k1+k2)-3*(L+1)*p2);
}
}
// Diagonalize and print the eigenvalues
HMX = KMX - (4.0/3.0)*VMX;
std::setprecision(8);
GeneralizedSelfAdjointEigenSolver<MatrixXd> es(HMX,NMX);
//pushback eigenvalues into vectors
eValues0.push_back(es.eigenvalues()[0]);
eValues1.push_back(es.eigenvalues()[1]);
eValues2.push_back(es.eigenvalues()[2]);
eValues3.push_back(es.eigenvalues()[3]);
}
// Calculate the minimum eigenvalue and the corresponding parameters
minEvalue0 = eValues0[0];
minEvalue1 = eValues1[0];
minEvalue2 = eValues2[0];
minEvalue3 = eValues3[0];
//cout << "\n";
//cout << "Here is the full list of eigenvalues: " << std::endl;
for(unsigned int i = 0; i < eValues0.size(); i++){
if (eValues0[i] < minEvalue0){
minEvalue0 = eValues0[i];
}
}
for(unsigned int i = 0; i < eValues1.size(); i++){
if (eValues1[i] < minEvalue1){
minEvalue1 = eValues1[i];
}
}
for(unsigned int i = 0; i < eValues2.size(); i++){
if (eValues2[i] < minEvalue2){
minEvalue2 = eValues2[i];
}
}
for(unsigned int i = 0; i < eValues3.size(); i++){
if (eValues3[i] < minEvalue3){
minEvalue3 = eValues3[i];
}
}
//push_back all the minimized eigenvalues into a vector
minEvalues.push_back(minEvalue0);
minEvalues.push_back(minEvalue1);
minEvalues.push_back(minEvalue2);
minEvalues.push_back(minEvalue3);
double eigenvalue = minEvalues[n-1]-p3;//return the desired eigenvalue from minEvalues
return eigenvalue;
}
// Define function that selects the state for diagonalization(n, L, p0, p1, p2, p3)
double mass (int state, double p0, double p1, double p2, double p3){
double n;
double L;
double mesonMass = 0.0;
if (state == 0){
n = 1.0;
L = 0.0;
...
else if (state == 6){
n = 4.0;
L = 2.0;
}
mesonMass = diagonalization(n, L, p0, p1, p2, p3);
return mesonMass;
}
double operator() (const std::vector<double> ¶ms) const {
double X2 = 0.0;
//std::cout << mass(0,params[0], params[1], params[2], params[3]) << std::endl;
//std::vector<double> mass;
for (int i=0; i<m_data[0].size(); i++) {
double n = m_data[0][i];
double L = m_data[1][i];
double f = mass(i, params[0], params[1], params[2], params[3]);
X2 += (m_data[2][i] - f)*(m_data[2][i] - f)/(m_data[2][i]);
std::cout << "data[2][i] = " << m_data[2][i] << std::endl;
std::cout << "mass[i]" << f << std::endl;
}
return X2;
}
//set to 1 for chi-sq and 0.5 for log-likelihood
double Up() const {return 1.0;}
};
int main(){
std::vector< std::vector<double> > data(3);
// (1S)
data[0].push_back(1.0);
data[1].push_back(0.0);
data[2].push_back(3.0687);
// (2S)
data[0].push_back(2.0);
data[1].push_back(0.0);
data[2].push_back(3.6740);
// (1P)
data[0].push_back(2.0);
data[1].push_back(1.0);
data[2].push_back(3.5259);
// (1D)
data[0].push_back(3.0);
data[1].push_back(2.0);
data[2].push_back(3.8282);//data[2].push_back(3.7737);
// (3S)
data[0].push_back(3.0);
data[1].push_back(0.0);
data[2].push_back(4.0390);//4.040);
// (4S)
data[0].push_back(4.0);
data[1].push_back(0.0);
data[2].push_back(4.4214);
// (2D)
data[0].push_back(4.0);
data[1].push_back(2.0);
data[2].push_back(4.1915);
std::vector<double> params(4,1.0);
std::vector<double> err(4,0.1);
params[0] = 1.60; //mc
params[1] = 0.50; //alphaS
params[2] = 0.14; //b
params[3] = 0.25; //mass shift
std::cout << "Initial values of parameters: " << std::endl;
for (int i = 0; i < 4; i++){
std::cout << "params[" << i << "] = " << params[i] << std::endl;
}
std::cout << "Computing the values of the parameters..." << std::endl;
//minimize
Minimizer fcn(data);
Minimizer example(data);
for(int i=0; i<7; i++){
std::cout << example.mass(i,params[0],params[1],params[2],params[3]) << std::endl;
}
//std::cout << "data[0].size " << data[0].size() << std::endl;
MnMigrad migrad(fcn, params, err);
FunctionMinimum minParams = migrad();
std::cout << " chi^2 = " << minParams.Fval() << std::endl;
std::cout << " for the following values:\n " << "mass c: " << minParams.UserState().Value(0) << ", alphaS: " << minParams.UserState().Value(1) << ", b: " \
<< minParams.UserState().Value(2) << ", shift: " << minParams.UserState().Value(3) << std::endl;
std::cout << "Process completed." << std::endl;
return 0;
}
When I print the result of example.mass() I get exactly the seven values that I expected.
3.06935
3.6354
3.9012
4.41176
4.0215
4.34403
4.68653
But std::cout << "mass[i] " << f << std::endl;
makes no sense at all to me, as it keeps printing eigenvalues:
mass[i]2.86003
mass[i]3.41612
mass[i]3.67878
mass[i]4.19538
mass[i]3.80451
mass[i]4.13111
mass[i]4.47528
mass[i]2.86004
mass[i]3.41613
mass[i]3.67879
mass[i]4.19539
mass[i]3.80452
mass[i]4.13112
mass[i]4.47528
...
mass[i]3.13612
mass[i]3.52386
mass[i]3.70813
mass[i]4.0931
mass[i]3.80947
mass[i]4.05334
mass[i]4.30513
etc. What I really don’t understand is why it keeps printing if data[0].size() is 7. Am I missing something about Minuit2 here is is it just my mistake in the way I wrote the code? You can see from my code that I am a Newbie, so any help would be appreciated.
Thank you in advance.
Christian
ROOT Version:
Platform: Linux
Compiler: gcc