# Segmentation Fault When Using FunctionMinimum from Minuit2

Hi there,

New to Minuit2. I am posting all of my code here only because I am not sure what’s relevant to my question and what’s not, and I don’t want to leave out some important info. Basically, what I am trying to do is to use Minuit2 to do the chi^2 minimization of for some matrix that I construct with Eigen. Everything compiles until I call the FunctionMinimum, in which case I get Segmentation Fault. Now, as far as I understand it, this kind of error occurs when I am trying to access memory that I am not allowed to access or that does not exist, but I am not sure why I am getting it here; I have the files with all the data in the same folder as the source code.

Any help will be greatly appreciated. (I am coding in Linux)

``````/*
1. Read input file containing the channel files and initial parameters
2. For each channel file read in the data and input parameters for that channel
3. Place the channel files into a vector<string>
4. For each channel file, place the data into a vector<vector<vector<double>>> (indexes are rows, columns and channels)
5. Construct the sigma vector (dimensions 1xn)
6. For each channel calculate the X2_alpha = (data[i][alpha] - sigma_alpha)^2/(data[i][alpha])^2
7. Sum all the X2_alpha; X2 = sum(X2_alpha)
8. Print out the minimized values
*/
#define USE_MATH_DEFINES
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/FCNBase.h"

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/Core>
#include <cmath>
#include <complex>

using namespace std;
using namespace Eigen;

const complex<double> im(0,1);
// Define the minimization class
using namespace ROOT::Minuit2;

class MinFcn : public FCNBase {

private:
const vector< vector< vector < double > > > &_data;
const vector< string > &_channels;
const vector< int > &_numRes;
const vector< double > &_finalMasses;
public:
MinFcn(const vector< vector< vector<double> > > & data, const vector< string > & channels, const vector< int > & numRes, const vector< double > & finalMasses)
: _data(data), _channels(channels), _numRes(numRes), _finalMasses(finalMasses) {}
~MinFcn() {}
double operator() (const std::vector<double> & params) const {

double X2 = 0.0;
// Define all the needed matrices
MatrixXd pf,Sigma;
MatrixXcd P, K, rho, I;
int n = _channels.size();
int R[_numRes.size()]; // number of resonances in channels
double m_f[n];
pf = MatrixXd(1,n);
P = MatrixXcd(1,n);
K = MatrixXcd(n,n);
rho = MatrixXcd(n,n);
I = MatrixXcd::Identity(n,n);
MatrixXcd A = (I + im*rho*K);
MatrixXcd T = P*A.inverse(); // 1xn matrix
MatrixXcd T2;
T2 = MatrixXcd(1,n);
Sigma = MatrixXd(1,n);
for (int j = 0; j < n; ++j){ // j is the j-th file, which is also the j-th column for the matrices
for (int i = 0; i < _data[j].size(); ++i){
double x = _data[j][i];
double y = _data[j][i];
double err = _data[j][i];
double pi = sqrt((x*x)/4- 0.000511*0.000511);

// the number of columns is equal to the number of channels, which is equal to the number of files, channels.size()
pf(0,j) = sqrt((x*x)/4 - m_f[j]*m_f[j]);

//P = MatrixXcd(1,n);
for (int resTerm = 0; resTerm < R[j]; ++resTerm){ //numRes.at(j) = number of resonances for channel j
P(0,j) += params[resTerm]*params[resTerm+R[j]*(j+1)]/(x*x - params[R[j]*(n+1)+resTerm]*params[R[j]*(n+1)+resTerm]+1e-8);
}

//K = MatrixXcd(n,n);
for (int row = 0; row < n; row++){
for (int resTerm = 0; resTerm < R[j]; resTerm++){
K(row,j) += params[resTerm+R[j]*(row+1)]*params[resTerm+R[j]*(row+1+j)]/(x*x - params[R[j]*(n+1)+resTerm]*params[R[j]*(n+1)+resTerm]+1e-8);
}
}

//rho = MatrixXcd(n,n);
for (int row = 0; row < n; row++){
if(row != j){
rho(row,j) = 0;
}
else {
rho(row,j) = sqrt((x*x/4) - m_f[row]*m_f[row]);
}
}

T2(0,j) = norm(T(0,j));
Sigma(0,j) = (3/(64*3.14159*x*x))*(pf(0,j)/pi)*T2(0,j).real(); // 1xn matrix
X2 += (_data[j][i]-Sigma(0,j))*(_data[j][i]-Sigma(0,j))/(_data[j][i]*_data[j][i]);
}
}
return X2;
}
double Up() const {return 1.0;}
};

int main() {
vector<string> channels;
vector<int> numRes;
vector<double> finalMasses;
ifstream inf ("files.dat");
if (inf.is_open()){
string line;
getline(inf,line);
while(inf.good()){
istringstream iss(line);
string channel;
int R;
double fMass;
iss >> channel >> R >> fMass;
// Place channel file names into a vector<string>
channels.push_back(channel);
// Place number of resonances into a vector<int>
numRes.push_back(R);
// Place the final-state masses into a vector<double>
finalMasses.push_back(fMass);
getline(inf,line);
}
}
// Now you have the names of the channels stored in channels = {"Babar_D+_D-.dat", "Babar_D0_D0.dat", etc.}
// and the resonances stored in numRes = {3, 3, etc.}
// Read each file name and open the file to read the data inside of it
// Iterate through the number of files in vector<string> = channels
// Place all of the data into a vector< vector < vector<double> > >
//vector< vector< double > > data2D(3);
vector< vector< vector< double > > > data3D;
for (int i = 0; i < channels.size(); ++i){
vector< vector< double > > data2D(3);
//vector< vector< vector< double > > > data3D;
ifstream ifs(channels.at(i));
if(!ifs){
cerr << "File could not be found." << endl;
exit(1);
}
else if(ifs.is_open()){
string line;
getline(ifs,line);
while(ifs.good()){
istringstream iss (line);
double x;
double y;
double err;
iss >> x >> y >> err;
data2D.push_back(x);
data2D.push_back(y);
data2D.push_back(err);
getline(ifs,line);
}
data3D.push_back(data2D);
//minimize
MinFcn fcn(data3D, channels, numRes, finalMasses);
vector<double> params(24,1.0);
vector<double> err(24,0.1);
//Define parameters for each channel
// Initial e+e- state
params=1.00; //gee1
params=1.00; //gee2
params=1.00; //gee3
// Final DD state
params=1.00; //gDD1
params=1.00; //gDD2
params=1.00; //gDD3
// Final DsDs state
params=1.00; //gDs1
params=1.00; //gDs2
params=1.00; //gDs3
// Final DstarDstar state
params=1.00; //gDstar1
params=1.00; //gDstar2
params=1.00; //gDstar3
// Final DsstarDsstar state
params=1.00; //gDsstar1
params=1.00; //gDsstar2
params=1.00; //gDsstar3
// Final DDstar state
params=1.00; //gDDstar1
params=1.00; //gDDstar2
params=1.00; //gDDstar3
// Final DsDsstar
params=1.00; //gDsDsstar1
params=1.00; //gDsDsstar2
params=1.00; //gDsDsstar3
// Resonance masses
params=3.750; //m1
params=4.140; //m2
params=4.450; //m3
}
for ( int j = 0; j < data3D[i].size(); ++j){
cout << data3D[i][j] << '\t' << data3D[i][j] << '\t' << data3D[i][j] << endl;
//cout << "---------------" << endl;
}
cout << "---------------" << endl;
}
return 0;
}``````

Hi,

can you share a stacktrace?

Cheers,
P

Hi,

I don’t know how to print a stacktrace. Is there a way to get it with g++?

Thanks.

Hi,

ROOT itself should be printing one. The gdb tool can help you if that does not happen.

Cheers

Hi,

Thanks again for your reply. I am not using ROOT. I am using the stand-alone version of Minuit2, in Ubuntu.

Hi,

then the only tool would be gdb.

P

Hi there,

I suppose I ran into the same problem.
And I’ll provide the stacktrace.

_ROOT Version:6.22.08
_Platform:Ubuntu 18.04.06
_Compiler: gcc 7.5.0
_Architecture: x86_64 Linux

``````dsanana.exe: /home/dsdaq/andrea/root_src/math/minuit2/src/VariableMetricBuilder.cxx:272: ROOT::Minuit2::FunctionMinimum ROOT::Minuit2::VariableMetricBuilder::Minimum(const ROOT::Minuit2::MnFcn&, const ROOT::Minuit2::GradientCalculator&, const ROOT::Minuit2::MinimumSeed&, std::vector<ROOT::Minuit2::MinimumState>&, unsigned int, double) const: Assertion `s0.IsValid()' failed.

__GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:51
51	../sysdeps/unix/sysv/linux/raise.c: No such file or directory.
(gdb) where
#0  __GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:51
#1  0x00007ffff4199921 in __GI_abort () at abort.c:79
#2  0x00007ffff418948a in __assert_fail_base (
fmt=0x7ffff4310750 "%s%s%s:%u: %s%sAssertion `%s' failed.\n%n",
assertion=assertion@entry=0x7ffff5575843 "s0.IsValid()",
file=file@entry=0x7ffff5575308 "/home/dsdaq/andrea/root_src/math/minuit2/src/VariableMetricBuilder.cxx", line=line@entry=272,
function=function@entry=0x7ffff5575b20 <ROOT::Minuit2::VariableMetricBuilder::Minimum(ROOT::Minuit2::MnFcn const&, ROOT::Minuit2::GradientCalculator const&, ROOT::Minuit2::MinimumSeed const&, std::vector<ROOT::Minuit2::MinimumState, std::allocator<ROOT::Minuit2::MinimumState> >&, unsigned int, double) const::__PRETTY_FUNCTION__> "ROOT::Minuit2::FunctionMinimum ROOT::Minuit2::VariableMetricBuilder::Minimum(const ROOT::Minuit2::MnFcn&, const ROOT::Minuit2::GradientCalculator&, const ROOT::Minuit2::MinimumSeed&, std::vector<ROOT:"...) at assert.c:92
#3  0x00007ffff4189502 in __GI___assert_fail (assertion=0x7ffff5575843 "s0.IsValid()",
file=0x7ffff5575308 "/home/dsdaq/andrea/root_src/math/minuit2/src/VariableMetricBuilder.cxx",
line=272,
function=0x7ffff5575b20 <ROOT::Minuit2::VariableMetricBuilder::Minimum(ROOT::Minuit2::MnFcn const&, ROOT::Minuit2::GradientCalculator const&, ROOT::Minuit2::MinimumSeed const&, std::vector<ROOT::Minuit2::MinimumState, std::allocator<ROOT::Minuit2::MinimumState> >&, unsigned int, double) const::__PRETTY_FUNCTION__> "ROOT::Minuit2::FunctionMinimum ROOT::Minuit2::VariableMetricBuilder::Minimum(const ROOT::Minuit2::MnFcn&, const ROOT::Minuit2::GradientCalculator&, const ROOT::Minuit2::MinimumSeed&, std::vector<ROOT:"...) at assert.c:101
#4  0x00007ffff555246e in ROOT::Minuit2::VariableMetricBuilder::Minimum(ROOT::Minuit2::MnFcn const&, ROOT::Minuit2::GradientCalculator const&, ROOT::Minuit2::MinimumSeed const&, std::vector<ROOT::Minuit2::MinimumState, std::allocator<ROOT::Minuit2::MinimumState> >&, unsigned int, double) const ()
from /home/dsdaq/andrea/root/lib/libMinuit2.so
#5  0x00007ffff5552ffa in ROOT::Minuit2::VariableMetricBuilder::Minimum(ROOT::Minuit2::MnFcn const&, ROOT::Minuit2::GradientCalculator const&, ROOT::Minuit2::MinimumSeed const&, ROOT::Minuit2::MnStrategy const&, unsigned int, double) const () from /home/dsdaq/andrea/root/lib/libMinuit2.so
#6  0x00007ffff55405ac in ROOT::Minuit2::ModularFunctionMinimizer::Minimize(ROOT::Minuit2::MnFcn const&, ROOT::Minuit2::GradientCalculator const&, ROOT::Minuit2::MinimumSeed const&, ROOT::Minuit2::MnStrategy const&, unsigned int, double) const () from /home/dsdaq/andrea/root/lib/libMinuit2.so
#7  0x00007ffff553f9bd in ROOT::Minuit2::ModularFunctionMinimizer::Minimize(ROOT::Minuit2::FCNBase const&, ROOT::Minuit2::MnUserParameterState const&, ROOT::Minuit2::MnStrategy const&, unsigned int, double) const () from /home/dsdaq/andrea/root/lib/libMinuit2.so
#8  0x00007ffff5504b50 in ROOT::Minuit2::MnApplication::operator()(unsigned int, double) ()
from /home/dsdaq/andrea/root/lib/libMinuit2.so
I can provide a link to git repo with the exact code, but the call to `MnMigrad` that segfaults is trivial.