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[1][i][alpha] - sigma_alpha)^2/(data[2][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/MnMigrad.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][0].size(); ++i){
double x = _data[j][0][i];
double y = _data[j][1][i];
double err = _data[j][2][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][1][i]-Sigma(0,j))*(_data[j][1][i]-Sigma(0,j))/(_data[j][2][i]*_data[j][2][i]);
}
}
return X2;
}
double Up() const {return 1.0;}
};
int main() {
// Read the input file
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[0].push_back(x);
data2D[1].push_back(y);
data2D[2].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[0]=1.00; //gee1
params[1]=1.00; //gee2
params[2]=1.00; //gee3
// Final DD state
params[3]=1.00; //gDD1
params[4]=1.00; //gDD2
params[5]=1.00; //gDD3
// Final DsDs state
params[6]=1.00; //gDs1
params[7]=1.00; //gDs2
params[8]=1.00; //gDs3
// Final DstarDstar state
params[9]=1.00; //gDstar1
params[10]=1.00; //gDstar2
params[11]=1.00; //gDstar3
// Final DsstarDsstar state
params[12]=1.00; //gDsstar1
params[13]=1.00; //gDsstar2
params[14]=1.00; //gDsstar3
// Final DDstar state
params[15]=1.00; //gDDstar1
params[16]=1.00; //gDDstar2
params[17]=1.00; //gDDstar3
// Final DsDsstar
params[18]=1.00; //gDsDsstar1
params[19]=1.00; //gDsDsstar2
params[20]=1.00; //gDsDsstar3
// Resonance masses
params[21]=3.750; //m1
params[22]=4.140; //m2
params[23]=4.450; //m3
MnMigrad migrad(fcn, params, err);
FunctionMinimum minParams = migrad();
}
for ( int j = 0; j < data3D[i][0].size(); ++j){
cout << data3D[i][0][j] << '\t' << data3D[i][1][j] << '\t' << data3D[i][2][j] << endl;
//cout << "---------------" << endl;
}
cout << "---------------" << endl;
}
return 0;
}