Fitting a function defined by an integral

ROOT
1

I am trying to fit a function of the form f(x) = Integral g(x,y ; a,b )dy with a data set. The limits of the integral are (0,+inf) but instead of inf I can put a value at which the function g(x,y=value) is basically 0.
Of course the integrand cannot be solved analythically otherwise there were not be any problem.
N.B ‘a’ and ‘b’ are free parameters need to be calculated from the fit.
N.B.B: the function g is a function of two variables g(x,y) but after integrating over y i get a 1D function of x.

I tried the following lines but it does not seem to work (just returns the initialization parameters and the fit always converges, it seems like he’s not trying to work out the optimal parameters):

Double_t ef(Double_t *x, Double_t *p) 
{
	return ( sin(2.*2.13*sin(x[0]/2.)*x[1])*x[1] )/( 1. + exp( ( x[1] - p[0] )/p[1] ) );
}

Double_t funct(Double_t *x, Double_t *p) 
{ 
	TF2 *f1=new TF2("f1",ef,0,M_PI,0,10*0.55+4.37,2); 
	f1->SetParameters(p[0],p[1]);
	TF12 *f12=new TF12("f12",f1,x[0],"y"); 
	double B = 5.8e-4; //fm^-4
	double momentum = 2.13; //fm^-1
	Double_t sum=f12->Integral(0,10*0.55+4.37); 
              delete f1;
	return (B/(momentum*momentum))*( pow(sum , 2) )/( pow( sin( x[0]/2. ) , 6 ));;
}

int main (){
.............
	TF1 *myFun=new TF1("myFun",funct,0,M_PI,2); 
	myFun->SetParName(0,"R");						//par names
	myFun->SetParName(1,"a");	
	myFun->SetParameter(0,4.37);
	myFun->SetParLimits(0,0.1,10);						//starting values and par limits
	myFun->SetParameter(1,0.55);
..........
TFitResultPtr r = myGraph -> Fit("myFun", "S");
........
}

Does someone knows how to work this out?
Thank you in advance.
ROOT Version: 6.18
Compiler: c++

2

I guess you need to set “reasonable” initial values for the parameters of your function (and its range).
Your function seems to be almost always almost zero, except close to zero, where it goes to infinity.

//
// root [0] .x trial.cxx
//

const double momentum = 2.13; // fm^-1

Double_t ef(Double_t *x, Double_t *p) {
  return x[1] * sin(2. * momentum * x[1] * sin(x[0] / 2.)) / (1. + exp((x[1] - p[0]) / p[1]));
}

TF2 *f1 = new TF2("f1", ef, 0., M_PI, 0., 10. * 0.55 + 4.37, 2);
TF12 *f12 = new TF12("f12", f1, 0., "y");

Double_t funct(Double_t *x, Double_t *p) {
  const double Bmm = 5.8e-4 / momentum / momentum; // fm^-4
  f1->SetParameters(p[0], p[1]);
  f12->SetXY(x[0]);
  Double_t sum = f12->Integral(0., 10. * 0.55 + 4.37);
  return Bmm * pow(sum, 2) / pow(sin(x[0] / 2.), 6);
}

TF1 *myFun = new TF1("myFun", funct, 0.01, M_PI, 2);

void trial() {
  myFun->SetNpx(10000);
  myFun->SetParNames("R", "a");
  myFun->SetParameters(4.37, 0.55);
  // myFun->SetParLimits(0, 0.1, 10.);
  myFun->Draw();
  gPad->SetLogy(1);
}
3

Thank you very much for your help,
btw I tried your code and drawing the functions with these parameters I find it pretty similar to my datapoints. My x-range goes from 0.523 to 1.221 and the y-coordinates of the data are enclosed in the range between 10^-5 to 1. And also I have double checked my hand calculations from which I extracted the initial values for the parameters and I am oretty sure they are ok…
But unfortunately the behaviour of the fitting algorithm is the same: it doesent seem like he is looking for the minimum of the chi2 function and it just converges to random values…

4

Attach your graph so that one can try to fit it.

5

transformed.txt (464 Bytes)
These are the data points (x, y, errx, erry).
Don’t look at the error bars.
By the way i checked again and perhaps I found better constants for the function (which I paste here in the code). Maybe now the starting values are not ok but from another calculation I expect similar values to the ones that I have inserted here… I hope the code is correct.
Thank you very much for your help :slight_smile:
I’m using the following code:

Double_t ef(Double_t *x, Double_t *p) {
  const double momentum = 405.02/197.; // fm^-1
  return x[1] * sin( abs( 2.*momentum* sin(x[0] / 2.) ) * x[1] ) / (1. + exp((x[1] - p[0]) / p[1]));
}


Double_t funct(Double_t *x, Double_t *p) {
  const double momentum = 405.02/197.; // fm^-1
  TF2 *f1 = new TF2("f1", ef, 0., M_PI, 0., 10. * 0.55 + 4.37, 2);
  TF12 *f12 = new TF12("f12", f1, 0., "y");
  double B =8.34e-5;
  const double Bmm = B / momentum / momentum; // fm^-4
  f1->SetParameters(p[0], p[1]);
  f12->SetXY(x[0]);
  double sum = f12->Integral(0., 10. * 0.55 + 4.37);
  return Bmm * pow(sum, 2) / pow(sin(x[0] / 2.), 6);
}

int main(int argc, char** argv)
{
	if (argc < 2)
	{
		cout << "Not enough parameters: ./eseguibile " << argv[1] << endl;
		return 1;
	}
	std::string fileName (argv[1]);
	fileName.erase( fileName.end()-4, fileName.end() );	

	TApplication* myApp = new TApplication("myApp", NULL, NULL);
	TCanvas* myC = new TCanvas(fileName.c_str(),fileName.c_str(),0,0,700,500);
	TGraphErrors* myGraph = new TGraphErrors(argv[1]);
	
	TF1 *myFun = new TF1("myFun", funct, 0.01, M_PI, 2);
  	myFun->SetParNames("R", "a");
  	myFun->SetParameters(1.12, 0.55);
  	myFun->SetParLimits(0, 0.1, 5);
	myFun->SetParLimits(1, 0.4, 1);
	TFitResultPtr r = myGraph -> Fit("myFun", "S");
	r->Print("V");
  	myC->cd();
	myGraph->Draw("AP");
	myFun->Draw("SAME");
	myGraph->GetXaxis()->SetRangeUser(0.523,1.221);
   	myGraph->GetYaxis()->SetRangeUser(0.00001e-3,1e-3);

	myC->Modified();
	myC->Update();
	gPad->SetLogy(1);
		
	myApp->Run();
	return 0;
}
6

I guess you really need good initial values of parameters (and maybe you need to make sure that the “y_max” used in the integral calculations is big enough).

//
// root [0] .x trial.cxx("/path/to/transformed.txt")
//

const double momentum = 405.02 / 197.; // fm^-1
const double B = 8.34e-5;
const double Bmm = B / momentum / momentum; // fm^-4

Double_t ef(Double_t *x, Double_t *p) {
  return x[1] * sin(2. * momentum * x[1] * sin(x[0] / 2.)) / (1. + exp((x[1] - p[0]) / p[1]));
}

TF2 *f1 = new TF2("f1", ef, 0., TMath::Pi(), 0., 1., 2);
TF12 *f12 = new TF12("f12", f1, 0., "y");

Double_t funct(Double_t *x, Double_t *p) {
  Double_t y_max = 15. * p[1] + p[0];
  // if (y_max < 10.) y_max = 10.;
  f1->SetRange(0., 0., TMath::Pi(), y_max);
  f1->SetParameters(p[0], p[1]);
  f12->SetRange(0., y_max);
  f12->SetXY(x[0]);
  Double_t sum = f12->Integral(0., y_max);
  return Bmm * pow(sum, 2) / pow(sin(x[0] / 2.), 6);
}

// TF1 *myFun = new TF1("myFun", funct, 0.01, TMath::Pi(), 2);
TF1 *myFun = new TF1("myFun", funct, 0.5, 2.5, 2);

void trial(const char *fname = "transformed.txt") {
  myFun->SetNpx(10000);
  myFun->SetParNames("R", "a");
  myFun->SetParameters(1.12, 0.55);
  myFun->SetParLimits(0, 0.1, 10.);
  myFun->SetParLimits(1, 0.1, 10.);
  
  new TCanvas("c1", "myFun + myGraph");
  myFun->Draw();
  gPad->SetLogy(1);
  
  TGraphErrors* myGraph = new TGraphErrors(fname);
  myGraph->Draw("same");
  
  myGraph->Fit(myFun, "EX0"); // "initial pre-fit"
  // myGraph->Fit(myFun); // "final fit"
  
  new TCanvas("c2", "f1");
  f1->Draw("colz");
}
7

Unfortunately the fit did not work, but at this point I suppose that the problem is in the data or in the analytic form of the function.
By the way thank you very much for your help, at least now i know that the code is correct.

8

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.