Home | News | Documentation | Download

Segmentation Violation in Minuit2

Hello everyone,

I am running Minuit2 to do a minimization. The code compiles fine but when I run it I get a Segmentation Violation error.

Below is the entire code:

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <iomanip>
#include <utility>
#include <Eigen/Dense>
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnScan.h"
#include "Minuit2/FCNBase.h"

using namespace ROOT::Minuit2;
using namespace Eigen;

const double PI = 4*atan(1);
// min and max values for variational parameters
const double alphaMin = 0.005;
const double alphaMax = 2.000;
const double betaMin = 0.005;
const double betaMax = 2.000;

// struct that defines the state's properties
struct state{//name, hplet, mass, width, n, ell, s, L, J
  std::string name;
  int hplet;
  double mass;
  double width;//maybe(?)
  int n;
  int ell;
  int s;
  int ELL;
  int jay;
  state(std::string Name, int Hplet, double Mass, double Width, int N, int l, int S, int L, int J)
    : name(Name), hplet(Hplet), mass(Mass), width(Width), n(N), ell(l), s(S), ELL(L), jay(J) {}
};

// vector containing all the particles
std::vector<state> lgt; //lgt = values from lattice gauge theory

void fillHBottom(std::vector<state> & lgt){
//lgt.push_back( state("H-name", H-plet, mass, width, n, ell, s, ELL, jay) )
}

void fillHCharm(std::vector<state> & lgt){
//lgt.push_back(state("H-name",H-pl,mass,width, n, ell,s,ELL,jay) )
  lgt.push_back(state("H1 1--", 1, 4.411, 0.0, 1, 0, 0, 1, 1));
  lgt.push_back(state("H1 0-+", 1, 4.279, 0.0, 1, 0, 1, 1, 0));
  lgt.push_back(state("H1 1-+", 1, 4.310, 0.0, 1, 0, 1, 1, 1));
  lgt.push_back(state("H1 2-+", 1, 4.456, 0.0, 1, 0, 1, 1, 2));
  lgt.push_back(state("H2 1++", 2, 4.470, 0.0, 1, 1, 0, 1, 1));
  lgt.push_back(state("H2 0+-", 2, 4.437, 0.0, 1, 1, 1, 1, 0));
  lgt.push_back(state("H2 1+-", 2, 4.438, 0.0, 1, 1, 1, 1, 1));
  lgt.push_back(state("H2 2+-", 2, 4.502, 0.0, 1, 1, 1, 1, 2));
  lgt.push_back(state("H3 0++", 3, 4.591, 0.0, 1, 1, 0, 0, 0));
  lgt.push_back(state("H3 1+-", 3, 4.571, 0.0, 1, 1, 1, 0, 1));
  lgt.push_back(state("H4 1++", 4, 4.623, 0.0, 1, 1, 0, 2, 2));
  lgt.push_back(state("H4 0+-", 4, 4.665, 0.0, 1, 1, 1, 2, 1));
  lgt.push_back(state("H4 1+-", 4, 4.631, 0.0, 1, 1, 1, 2, 2));
  lgt.push_back(state("H4 2+-", 4, 4.644, 0.0, 1, 1, 1, 2, 3));
}

// Define the omega(k) function - gluon dispersion
double omega (const std::vector<double> &params, double k){
// mg = params[3], bg = params[4]
	return sqrt(k*k + params[3]*params[3]*exp(-k/params[4]));
}

// Define the 1D iterated integral for the gluon kinetic energy - use midpoint
double iteratedKg (const std::vector<double> &params, double ki, double kf, double beta, int n, int M1, int M2){
	double widthk = (kf-ki)/n;
	double sum =  0.0;
	for (int i = 1; i < n; i++){
        // SHO basis
		sum += (4.0/(3.0*sqrt(PI)*pow(M1*M2,1.25)*pow(beta,5.0)))* ( pow((ki+(i+0.5)*widthk),4.0) * (omega(params, ki+(i+0.5)*widthk) + (ki+(i+0.5)*widthk)*(ki+(i+0.5)*widthk)/omega(params, ki+(i+0.5)*widthk)) * exp((-1)*(ki+(i+0.5)*widthk)*(ki+(i+0.5)*widthk)*(M1+M2)/(2*(M1*M2)*beta*beta)) );
	}
	return widthk*sum;
}

double diagonalization(const state &hybrid, const std::vector<double> &params) {
  //double mq = params[0]; // quark mass
  //double aS = params[1]; // Coulomb parameter
  //double b = params[2]; // confining parameter
  //double mg = params[3]; // gluon mass
  //double bg = params[4]; // gluon parameter

  // Declare the number of terms in the expansion and the number of iterations for the numerical integrals
  int N1 = 1;
  int iterations = 100;

  // Define the matrices and the eigenvector - the expressions work for any H-plet
    MatrixXd H(N1*N1, N1*N1);
    VectorXd eVect;
    MatrixXd Kq(N1*N1, N1*N1);
    MatrixXd Vqq(N1*N1, N1*N1);
    MatrixXd Kg(N1*N1, N1*N1);
    MatrixXd Vqg0(N1*N1, N1*N1);
    MatrixXd Vqg2(N1*N1, N1*N1);
    MatrixXd N(N1*N1, N1*N1);

  //Declare the variables for the minimum eigenvalue
  double minEvalue = 0.0;
  double minAlpha;
  double minBeta;
  std::vector<double> eValues;
  std::vector<double> alphas;
  std::vector<double> betas;
  std::vector<double> minEvalues;
  std::vector<double> eVector;
  std::vector<double> minEVector;
  std::vector<std::vector<double> > eVectors;

  // Loop over the variational parameters and the matrix elements (k1=n', l1=m', k2=n, l2=m).
	for	(double beta = betaMin; beta <= betaMax; beta += 0.01){
	  for (double alpha = alphaMin; alpha <= alphaMax; alpha += 0.01){
	    for (int k1 = 1; k1 <= N1; k1++){
	      for (int k2 = 1; k2 <= N1; k2++){
					for (int l1 = 1; l1 <= N1; l1++){
		  			for(int l2 = 1; l2 <= N1; l2++){
		    			// Quark kinetic energy
		    			Kq(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = (4.0*pow(2.0,hybrid.ell)*pow(k1*k2,0.75 + hybrid.ell/2.0)*pow(l1*l2,1.25)/(params[0]*pow((k1+k2),2.5 + hybrid.ell)*pow((l1+l2),3.5))) * (4.0*(l1+l2)*(2.0*params[0]*params[0]*(k1+k2)+(pow(2.0,hybrid.ell+1) + 1.0)*k1*k2*alpha*alpha)+5.0*l1*l2*(k1+k2)*beta*beta);
		    			// Quark interaction term
		    			Vqq(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = (pow(2.0, 2*hybrid.ell+2.0)*sqrt(2.0/PI)*(pow(k1*k2,0.75 + hybrid.ell/2.0)*pow(l1*l2,1.25))/(pow(3.0,hybrid.ell + 1)*alpha*pow((k1+k2),2.0 + hybrid.ell)*pow((l1+l2),2.5))) * (2.0*alpha*alpha*params[1]*(k1+k2) - pow(2.0, hybrid.ell)*3.0*params[2]);
              // Quark-gluon term for H1
              Vqg0(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = (pow(2.0,3*hybrid.ell + 1)*sqrt(2.0/PI)*pow(k1*k2,0.75 + hybrid.ell/2.0)*pow(l1*l2,1.25)/(pow(k1+k2,2.0 + hybrid.ell)*pow(l1+l2,3.0)*alpha*beta*pow((4.0*(k1+k2)*alpha*alpha+(l1+l2)*beta*beta),1.5 + hybrid.ell)))*( 192.0*params[2]*pow(alpha,4.0 + 2.0*hybrid.ell)*pow((k1+k2),2.0 + hybrid.ell) - 4.0*pow((k1+k2),hybrid.ell + 1.0)*(l1+l2)*pow(alpha,2.0 + 2.0*hybrid.ell)*beta*beta*(16.0*params[1]*(k1+k2)*alpha*alpha - (pow(2.0,hybrid.ell+1.0)+1.0)*7.0*params[2]) + pow(2.0, hybrid.ell + 1.0)*pow(beta,4.0 + 2.0*hybrid.ell)*pow((l1+l2),2.0 + hybrid.ell)*((pow(2.0,2.0 + hybrid.ell) -1.0)*params[2]-8.0*params[1]*(k1+k2)*alpha*alpha) - hybrid.ell*pow((l1+l2),3.0)*(4.0*params[1]*(k1+k2)*alpha*alpha - 3.0*params[2])*pow(beta,6.0) );
              // Numerical gluon kinetic energy
              Kg(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = (pow(2.0, hybrid.ell + 1)*sqrt(2.0)*pow(k1*k2,0.75 + hybrid.ell/2.0)/pow(k1+k2,1.5 + hybrid.ell)) * iteratedKg(params,0.0001,10.0,beta,iterations,l1,l2);
              // Overlap matrix
              N(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = 16.0*pow(2.0,hybrid.ell)*pow(k1*k2,0.75 + hybrid.ell/2)*pow(l1*l2,1.25)/(pow((k1+k2),1.5 + hybrid.ell/2.0)*pow((l1+l2),2.5));
		    			// Quark-gluon second term for H2-3-4
              Vqg2(k1-1+N1*(l1-1),k2-1+N1*(l2-1)) = (-80.0)*sqrt(2.0/PI)*alpha*beta*pow(k1*k2*l1*l2,1.25)*( 4.0*params[2]*(k1+k2)*alpha*alpha + (l1+l2)*beta*beta*(params[2]+8.0*params[1]*(k1+k2)*alpha*alpha) )/(pow((k1+k2)*(l1+l2),2.0)*pow((4.0*(k1+k2)*alpha*alpha+(l1+l2)*beta*beta),2.5));
						}
					}
				}
			}

		// Diagonalize
		// Print the matrix and eigenvalue
		// Color factors are already included in the expressions and the Vqg is multiplied by 2 to account for the qbar g term
    	if (hybrid.hplet == 1){
      	H = Kq + Kg + Vqq + 2.0*Vqg0;
      }
      else if (hybrid.hplet == 2){
      	H = Kq + Kg + Vqq + 2.0*(Vqg0 + 0.1*Vqg2);
      }// The expressions for H3 and H4 are the same as for H2, only the value of total L changes (hence the coefficient of Vqg).
      else if (hybrid.hplet == 3){
      	H = Kq + Kg + Vqq + 2.0*(Vqg0 - 0.2*Vqg2);
      }
      else if (hybrid.hplet == 4){
      	H = Kq + Kg + Vqq + 2.0*(Vqg0 - 0.02*Vqg2);
      }
			
			std::setprecision(8);
      GeneralizedSelfAdjointEigenSolver<MatrixXd> es(H,N);
      eVect = es.eigenvectors().col(hybrid.n - 1);
			// This maps an Eigen::VectorXd to an std::vector<double>
      eVector.resize(eVect.size());
      Map<VectorXd>(eVector.data(), eVector.size()) = eVect;
      // Save all values in vectors
      eValues.push_back(es.eigenvalues()[hybrid.n - 1]);
      alphas.push_back(alpha);
      betas.push_back(beta);
		}
	}

	// Calculate the minimum eigenvalue and the corresponding parameters
	minEvalue = eValues[0];
  minAlpha = alphas[0];
  minBeta = betas[0];
  minEVector = eVectors[0];

  for(int i = 0; i < eValues.size(); i++){
  	if (eValues[i] < minEvalue){
    	minEvalue = eValues[i];
			minAlpha = alphas[i];
      minBeta = betas[i];
      minEVector = eVectors[i];
    }
  }

  return minEvalue;

}

// Define average error routine
double avgError(const std::vector<state> &lgt, const std::vector<double> &params){
  double avgErr = 0.0;
  double sum = 0.0;
  int ct = 0;
  for (int i=0; i< lgt.size(); i++) {
    double mass_exp = lgt[i].mass;
    double mass_thy = diagonalization(lgt[i], params);
    sum += std::abs(mass_exp - mass_thy);
    ct++;
  }
  avgErr = sum/ct;
  return 1000*avgErr; //MeV
}

// Define minimizer class (inherits from FCNBase)
class Minimizer : public FCNBase {

private:

  const std::vector<state> &m_lgt;

public:

  Minimizer(const std::vector<state> &lgt) : m_lgt(lgt) {}

  ~Minimizer() {}
  double operator() (const std::vector<double> &params) const {
    double X2 = 0.0;
    for (int i=0; i< m_lgt.size(); i++) {
      double mass_exp = m_lgt[i].mass;
      double mass_thy = diagonalization(m_lgt[i], params);
      X2 += (mass_exp - mass_thy)*(mass_exp - mass_thy);
      std::cout << params[0] << " || " << params[1] << " || " << params[2] << " || " << params[3] << " || " << params[4] << std::endl;
      //std::cout << "exp: " << mass_exp << " || thy: " << mass_thy << std::endl;
    }
    std::cout << "chi^2 = " << X2 << std::endl;
    std::cout << "\n";
    return X2;
  }
  //set to 1 for chi-sq and 0.5 for log-likelihood
  double Up() const {return 1.0;}
};

int main(){
  std::vector<state> lgt;
  fillHCharm(lgt);
  //fillHBottom(lgt);

  //Create FCN function
  Minimizer fcn(lgt);
  //mass q: 4.78514, alphaS: 0.382944, b: 0.170503, shift: 1.67881
  //mass q: 4.7851, alphaS: 0.382327, b: 0.170471, shift: 1.66752
  //mass q: 1.55998, alphaS: 0.661734, b: 0.121585
  //mass q: 1.55806, alphaS: 0.633959, b: 0.120329, sigma: 0.636216
  //mass q: 4.79024, alphaS: 0.39042, b: 0.169599, sigma: 1.51675
  //std::vector<double> inPars = {4.80, 0.50, 0.17, 0.6, 6.0};
  std::vector<double> inPars = {1.60, 0.594, 0.16, 0.6, 6.0};
  //Create parameters
  double mq = inPars[0];
  double alphaS = inPars[1];
  double b = inPars[2];
  double mg = inPars[3];
  double bg = inPars[4];

  double err = 0.2; // initial error of params

  MnUserParameters uParams;//use params[i]/5 relative error
  uParams.Add("mq", mq, err*mq);
  uParams.Add("alphaS", alphaS, err*alphaS);
  uParams.Add("b", b, err*b);
  uParams.Add("mg", mg, err*mg);
  uParams.Add("bg", bg, err*bg);

  //MnStrategy strategy;
  //strategy.SetHighStrategy();
std::cout << "Hello." << std::endl;
  std::cout << diagonalization(lgt[0], inPars) << std::endl;

/*
  MnMigrad migrad(fcn, uParams);

  FunctionMinimum min = migrad();
  std::cout << "Ciao" << std::endl;

 std::cout << "\n";
  std::cout << "Variational parameter range: [" << alphaMin << ", " << alphaMax << "]" << std::endl;
  std::cout << "CHISQUARED MINIMIZATION RESULTS:" << std::endl;
  std::cout << "\n";
  std::cout << min << std::endl;
  std::cout << "\n";
  std::cout << "****************************************" << std::endl;
  std::cout << "Initial values of parameters:" << "\n";
  std::cout << "mass q: " << inPars[0] << ", alphaS: " << inPars[1] << ", b: " << inPars[2] << ", sigma: " << inPars[3] << std::endl;
  std::cout << "Initial error step-size: " << err << std::endl;
  //std::cout << "Initial average error: " << avgError(pdg, inPars) << " MeV" << std::endl;
  std::cout << " chi^2 = " << min.Fval() << std::endl;
  std::cout << " for the following optimized values:\n " << "mass q: " << min.UserState().Value(0) << ", alphaS: " << min.UserState().Value(1) << ", b: " \
            << min.UserState().Value(2) << ", mg " << min.UserState().Value(3) << ", bg " << min.UserState().Value(4) << std::endl;

  std::vector<double> outPars;
  outPars.push_back(min.UserState().Value(0));
  outPars.push_back(min.UserState().Value(1));
  outPars.push_back(min.UserState().Value(2));
  outPars.push_back(min.UserState().Value(3));
  outPars.push_back(min.UserState().Value(4));
  //std::cout << "Final average error: " << avgError(pdg, outPars) << " MeV" << std::endl;
  std::cout << "****************************************" << std::endl;
  std::cout << "\n";
  std::cout << std::setw(10) << "Name " << std::setw(15) << "Exp mass " << std::setw(18) << "Init. thy mass " << std::setw(18) << "Fin. thy mass" << std::endl;
  for (int i = 0; i < lgt.size(); i++){
    std::cout << std::setw(10) << lgt[i].name << std::setw(15) << lgt[i].mass << std::setw(15) << std::setprecision(8) << diagonalization(lgt[i], inPars)
              << std::setw(15) << diagonalization(lgt[i], outPars) << std::endl;
  }
  std::cout <<"****************************************" << std::endl;
  std::cout << "\n";
 std::cout << "Ciao" << std::endl;
*/
 return 0;
}

and here is the entire stack trace

Hello.

 *** Break *** segmentation violation



===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0  0x00007fb8db5e8687 in __GI___waitpid (pid=22522, stat_loc=stat_loc
entry=0x7ffcaa909028, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:30
#1  0x00007fb8db553067 in do_system (line=<optimized out>) at ../sysdeps/posix/system.c:149
#2  0x00007fb8da690893 in TUnixSystem::StackTrace() () from /home/christian/root/lib/libCore.so.6.24
#3  0x00007fb8da6933a5 in TUnixSystem::DispatchSignals(ESignals) () from /home/christian/root/lib/libCore.so.6.24
#4  <signal handler called>
#5  0x000055fbf8fdacf4 in std::vector<double, std::allocator<double> >::size() const ()
#6  0x000055fbf8fdaadd in std::vector<double, std::allocator<double> >::operator=(std::vector<double, std::allocator<double> > const&) ()
#7  0x000055fbf8fd7c5c in diagonalization(state const&, std::vector<double, std::allocator<double> > const&) ()
#8  0x000055fbf8fd854c in main ()
===========================================================


The lines below might hint at the cause of the crash.
You may get help by asking at the ROOT forum https://root.cern.ch/forum
Only if you are really convinced it is a bug in ROOT then please submit a
report at https://root.cern.ch/bugs Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5  0x000055fbf8fdacf4 in std::vector<double, std::allocator<double> >::size() const ()
#6  0x000055fbf8fdaadd in std::vector<double, std::allocator<double> >::operator=(std::vector<double, std::allocator<double> > const&) ()
#7  0x000055fbf8fd7c5c in diagonalization(state const&, std::vector<double, std::allocator<double> > const&) ()
#8  0x000055fbf8fd854c in main ()
===========================================================

I have never gotten an error like this, so I don’t know what do make of it. If I understand it correctly, it is telling me that there is some sort of problem with my diagonalization() routine, though I checked and cannot find any.

I ran the diagonalization() routine independently with another program and it runs just fine though. Any suggestions or advice?

Thank you.
Christian

_ROOT Version: 6.24.02
_Platform: Linux Mint
_Compiler: g++


Hi,
It is maybe a problem with the vector size. Maybe Minuit2 things the vector has less dimension that what is expected in your diagonalisation routine. I would check if all the parameters have been defined correctly and the size of the parameter vector that Minuit passes to the FCN function.

Cheers

Lorenzo

Hi Lorenzo,

Thank you for your reply. So you are saying that maybe Minuit2 is not recognizing the size of the prams vector that I pass to it? That is strange, bcause I have another very similar program in which I pass the same type of vector to the diagonalization() routing and I don’t get that error.

double diagonalization(const state &meson, const std::vector<double> &params) {
  double mq = params[0];
  double aS = params[1];
  double b = params[2];
  double sigma = params[3];
  double hypConst = 32.0*PI*aS/(9.0*mq*mq);
  // Declare the number of terms in the expansion and the number of iterations for the numerical integrals
  int N1 = 12;

  // 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);
  MatrixXd Vhyp(N1, N1);
  VectorXd eVect;

  //Declare the variables for the minimum eigenvalue
  //double eigenMass = 0.0;
  double minEvalue = 0.0;
  double minAlpha;
  std::vector<double> eValues;
  std::vector<double> alphas;
  std::vector<double> minEvalues;
  std::vector<double> eVector;
    std::vector<double> minEVector;
  std::vector<std::vector<double> > eVectors;

  // Loop over the variational parameters and the matrix elements (k1=n', l1=m', k2=n, l2=m).
  for (double alpha = alphaMin; alpha <= alphaMax; 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,meson.ell+1)*sqrt(2.0)*pow(k1*k2,0.75 + meson.ell/2.0)/(pow(k1+k2,1.5 + meson.ell));
        // Mass and kinetic terms
        KMX(k1-1,k2-1) = NMX(k1-1,k2-1) * (2.0*mq + (3.0+2.0*meson.ell)*(k1*k2)*alpha*alpha/(mq*(k1+k2)));
        // Potential term
        VMX(k1-1,k2-1) = (powerfactorial(meson.ell)*pow(k1*k2,0.75 + static_cast<double>(meson.ell)/2.0)/(doublefactorial(2*meson.ell+1)*sqrt(PI)*alpha*pow(k1+k2,2.0+m\
eson.ell)))*(2*aS*alpha*alpha*(k1+k2)-3*(meson.ell+1)*b);
        Vhyp(k1-1,k2-1) = hypConst*pow(2.0, meson.ell+1)*pow(k1*k2, 0.75 + static_cast<double>(meson.ell)/2)*sqrt(2.0/PI)*pow(sigma, 3.0)*pow(alpha, 3+2*meson.ell)/(po\
w(pow(k1+k2, 1.5+meson.ell)*alpha*alpha + 2.0*sigma*sigma, 1.5 + meson.ell));
        HMX(k1-1,k2-1) = KMX(k1-1,k2-1) - (4.0/3.0)*VMX(k1-1,k2-1) + (meson.s*(meson.s + 1)/2.0 - 0.75)*Vhyp(k1-1,k2-1);
      }
    }

    // Diagonalize and print the eigenvalues
    //HMX = KMX - (4.0/3.0)*VMX;// + (meson.s*(meson.s + 1)/2.0 - 0.75)*Vhyp;

      std::setprecision(8);
      GeneralizedSelfAdjointEigenSolver<MatrixXd> es(HMX,NMX);
      //if (es.eigenvalues()[meson.n - 1] < minEvalue){
      //minEvalue = es.eigenvalues()[meson.n - 1];
      //}
      // choose the eigenvector associated with the n-1 eigenvalue
      eVect = es.eigenvectors().col(meson.n - 1);
      // This maps an Eigen::VectorXd to an std::vector<double>
      eVector.resize(eVect.size());
      Map<VectorXd>(eVector.data(), eVector.size()) = eVect;
      // Save all values in vectors
      eValues.push_back(es.eigenvalues()[meson.n -1]);
      alphas.push_back(alpha);
      eVectors.push_back(eVector);
  }
  // Calculate the minimum eigenvalue, eigenvector and the parameter
  minEvalue = eValues[0];
  minAlpha = alphas[0];
  minEVector = eVectors[0];
  for(unsigned int i = 0; i < eValues.size(); i++){
    if (eValues[i] < minEvalue){
      minEvalue = eValues[i];
      minAlpha = alphas[i];
      minEVector = eVectors[i];
    }
  }

  double SdotL = (meson.jay*(meson.jay+1) - meson.ell*(meson.ell+1) - meson.s*(meson.s+1))*VspinOrbit(eVector, 0.01, 20.0, minAlpha, 1000, meson.ell, params)/2.0;
  double T = 0.0;
  if (meson.s == 1 && meson.ell > 0){
    if (meson.jay == meson.ell - 1){
      T = -((meson.ell + 1)/(6.0*(2*meson.ell - 1)))*Vtensor(eVector, 0.01, 20.0, minAlpha, 1000, meson.ell, params);
    }
    else if (meson.jay == meson.ell){
      T = (1.0/6.0)*Vtensor(eVector, 0.01, 20.0, minAlpha, 1000, meson.ell, params);
    }
    else if (meson.jay == meson.ell + 1){
      -((meson.ell)/(6.0*(2*meson.ell + 3)))*Vtensor(eVector, 0.01, 20.0, minAlpha, 1000, meson.ell, params);
    }
  }

  //std::cout << T << std::endl;
  return minEvalue + SdotL + T;
}

The function above is the similar diagonalization() function I use in the other program without getting an error. The params vector is defined and passed the same way in both cases.

Thanks,
Christian

Hi,
Are you still having a problem here ? If this is the case I would suggest using a debugger or valgrind to understand better where the problem originates.
Cheers
Lorenzo