Segmentation Violation in Minuit2

Newbie
1

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++

2

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

3

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

4

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