Eval() on TGraph interpolates in a curve

[atha.stam]

---

ROOT Version: 6.18/04
Platform: Ubuntu 18.04
Compiler: gcc 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04)

---

Dear co-rooters,

I’m facing an issue with TGraph::Eval(). I’m creating a TGraph from an ascii file and my goal is to change the x-axis grid (aka binning) of that TGraph. The new x-values are also read from another ascii file.

I thought the easiest way is to use TGraph::Eval and the TGraph TGraph(const char* filename, const char* format = "%lg %lg", Option_t* option = "") constructor, however the interpolation isn’t linear but it curves, as shown below in red line/points

My code is the following(I’m including the data files)
Any idea on what might be wrong?

Thanks in advance!

#include "TGraph.h"
#include <Math/Interpolator.h>

#include <iostream>
#include <fstream>
#include <vector>

void interpolate(){

	TGraph *g_ENDF = new TGraph("Hf_ENDF.dat", "%lg %lg");
	TGraph *g_new  = new TGraph("Hf_new.dat" , "%lg %lg");
	
	g_ENDF->Sort();
	g_new ->Sort();
	
//	std::ifstream input("Hf_ENDF.dat");
//	std::vector<double> EN, XS;
//	int   lines=0;
//	double en, xs;
//	while ( !input.eof() ){
//		input >> en >> xs;
//		EN.push_back(en);
//		XS.push_back(xs);
//		lines++;
//	}
//		
//	TGraph *g_ENDF = new TGraph(lines, &EN[0], &XS[0]);
	g_ENDF->Draw("APL");
	g_ENDF->SetLineColor(kBlack);
	g_ENDF->SetLineWidth(4);
	
//		
//	ROOT::Math::Interpolator inter (lines, ROOT::Math::Interpolation::kLINEAR);
//	inter.SetData(lines, &EN[0], &XS[0]);


//	std::ifstream input_new("Hf_new.dat");
//	std::vector<double> EN_new, XS_new;
//	int   lines_new=0;
//	double en_new, xs_new;
//	while ( !input_new.eof() ){
//		input_new >> en_new >> xs_new;
//		EN_new.push_back(en_new);
//		XS_new.push_back( inter.Eval( en_new ) );
//		lines_new++;
//	}
////	
//	TGraph *g_new = new TGraph(lines_new, &EN_new[0], &XS_new[0]);
	g_new->SetLineColor(kRed);
	g_new->SetLineWidth(2);
	
	int     n = g_new->GetN();
	double *x = g_new->GetX();
	double *y = g_new->GetY();
	
	int     n_endf = g_ENDF->GetN();
	double *x_endf = g_ENDF->GetX();
	double *y_endf = g_ENDF->GetY();
	
	for (int i=0; i<n; ++i){
		y[i] = g_ENDF->Eval( x[i] );
		g_new->SetPoint(i, x[i], y[i]);
	}
	
	g_new->Draw("PLsame");
	g_new->SetName("g_new");
	//g_new->Print();

}	

Interpolation.zip (149.1 KB)

[couet]

When I run the macro you posted I get an empty canvas.

[dastudillo]

See how Eval works. If it doesn’t do what you want, you need another way.

[atha.stam]

That’s bizarre… Maybe it doesn’t read the files properly?
I recopied the macro just in case, something went wrong with pasting.

Also, I’m first compiling it

.L interpolate.C++

and then executing it

interpolate()

[atha.stam]

TGraph::Eval is supposed to do a liner interpolation between points, right?
At least if you run it like Eval( x ) or Eval(x, 0, "").
Am I missing something?

It seems to be working nicely in flat areas, but when the y values rapidly change, I get this curvy behaviour…

[atha.stam]

I also tried with ROOT::Math::Interpolator and although I sorted my file using sort -k 1g Hf_ENDF.dat and checked for duplicates using awk 'a[$1]++{print $1}' Hf_ENDF.dat, the interpolator complains

Error in <GSLError>: Error 4 in interp.c at 83 : x values must be strictly increasing
Warning in <ROOT::Math::GSLInterpolator::Eval>: input domain error

Any idea on what might be wrong?

Thanks in advance!

The version of the code with the interpolator is:

#include "TGraph.h"
#include <Math/Interpolator.h>

#include <iostream>
#include <fstream>
#include <vector>

void interpolate(){


	std::ifstream input("Hf_ENDF.dat");
	std::vector<double> EN, XS;
	int   lines=0;
	double en, xs;
        while ( 1 ){
		input >> en >> xs;
		if( input.eof() ) break;
		cout << en << ", " << xs << endl;
		EN.push_back(en);
		XS.push_back(xs);
		lines++;
	}
		
	TGraph *g_ENDF = new TGraph(lines, &EN[0], &XS[0]);
	g_ENDF->Draw("APL");
	g_ENDF->SetLineColor(kBlack);
	g_ENDF->SetLineWidth(4);
	
		
	ROOT::Math::Interpolator inter (lines, ROOT::Math::Interpolation::kLINEAR);
	inter.SetData(lines, &EN[0], &XS[0]);
	
	cout << inter.Eval( 1.666e-1 ) << endl;

	std::ifstream input_new("Hf_new.dat");
	std::vector<double> EN_new, XS_new;
	int   lines_new=0;
	double en_new, xs_new;
	while ( !input_new.eof() ){
		input_new >> en_new >> xs_new;
		EN_new.push_back(en_new);
		XS_new.push_back( inter.Eval( en_new ) );
		lines_new++;
	}
	
	TGraph *g_new = new TGraph(lines_new, &EN_new[0], &XS_new[0]);
	g_new->SetLineColor(kRed);
	g_new->SetLineWidth(2);

	g_new->Draw("PLsame");
	g_new->SetName("g_new");

}

[dastudillo]

1 - Read the second file (HF_new.dat) in the same way you read the first one.
2 - The first entry in the second file is outside the range (in the x-axis) of the first one, so there cannot be an “interpolation”; remove that point or increase its x (or add a lower-x point in the first file).

[atha.stam]

Thanks for your reply.
I did what you recommended and I still don’t get a reasonable linear interpolation.

I’m re-attaching the files and the updated code.
Can you replicate the issue?
I also tried the same code on lxplus (root version 6.22 I believe) and I get the same behaviour.

#include "TH1D.h"
#include "TGraph.h"
#include <Math/Interpolator.h>

#include <iostream>
#include <fstream>
#include <vector>

void interpolate(){
	
	double bin_content = 0.;
	std::ifstream input_ENDF("Hf_ENDF.dat");
	std::vector<double> EN_ENDF, XS_ENDF;
	int   lines_ENDF=0;
	double en_ENDF, xs_ENDF;
	while (1){
		input_ENDF >> en_ENDF >> xs_ENDF;
		if( input_ENDF.eof() ) break;
		EN_ENDF.push_back(en_ENDF);
		XS_ENDF.push_back(xs_ENDF);
		lines_ENDF++;
	}	
		
	TGraph *g_ENDF = new TGraph(lines_ENDF, &EN_ENDF[0], &XS_ENDF[0]);
	g_ENDF->Draw("APL");
	g_ENDF->SetLineColor(kBlack);
	g_ENDF->SetLineWidth(4);
    g_ENDF->GetXaxis()->SetLimits(1e-6, 1e-5);|
    g_ENDF->SetMarkerStyle(20);|
    g_ENDF->SetMarkerSize(3);|


	
		
	ROOT::Math::Interpolator inter (lines_ENDF, ROOT::Math::Interpolation::kLINEAR);
	inter.SetData(lines_ENDF, &EN_ENDF[0], &XS_ENDF[0]);


	std::ifstream input_new("Hf_new_binning.dat");
	std::vector<double> EN_new, XS_new;
	int   lines_new=0;
	double en_new, xs_new;
	while ( !input_new.eof() ){
		input_new >> en_new >> xs_new;
		if( input_new.eof() ) break;
		EN_new.push_back(en_new);
		bin_content = inter.Eval(en_new);
		XS_new.push_back( bin_content );
		lines_new++;
	}
	
	TGraph *g_new = new TGraph(lines_new, &EN_new[0], &XS_new[0]);
	g_new->SetLineColor(kRed);
	g_new->SetLineWidth(2);
	
	g_new->Draw("Psame");
	g_new->SetName("g_new");
	g_new->SetMarkerColor(kRed);
	g_new->SetMarkerStyle(20);

}

The thing is that it looks fine when the y-axis isn’t log. Is it maybe related?

Interpolation.zip (130.8 KB)

[atha.stam]

Thanks for the reply.
I used your way of “filling” the TGraphs and setting the logy axis but still I see the"curvy" interpolation.

#include "TH1D.h"
#include "TGraph.h"
#include "TPad.h"
#include <Math/Interpolator.h>

#include <iostream>
#include <fstream>
#include <vector>

void interpolate(){

	double bin_content = 0.;
	std::ifstream input_ENDF("Hf_ENDF.dat");
	std::vector<double> EN_ENDF, XS_ENDF;
	int   lines_ENDF=0;
	double en_ENDF, xs_ENDF;
	while (1){
		input_ENDF >> en_ENDF >> xs_ENDF;
		if( input_ENDF.eof() ) break;
		EN_ENDF.push_back(en_ENDF);
		XS_ENDF.push_back(xs_ENDF);
		lines_ENDF++;
	}
	
	
		
	TGraph *g_ENDF = new TGraph(lines_ENDF);//, &EN_ENDF[0], &XS_ENDF[0]);
	for(int i=0; i<lines_ENDF;++i){
		g_ENDF->SetPoint(i, EN_ENDF[i], XS_ENDF[i]);
	}
	g_ENDF->Draw("AL*");
	g_ENDF->SetLineColor(kBlack);
	g_ENDF->SetLineWidth(4);

	
		
	ROOT::Math::Interpolator inter (lines_ENDF, ROOT::Math::Interpolation::kLINEAR);
	inter.SetData(lines_ENDF, &EN_ENDF[0], &XS_ENDF[0]);

	std::ifstream input_new("Hf_new_binning.dat");
	std::vector<double> EN_new, XS_new;
	int   lines_new=0;
	double en_new, xs_new;
	while ( !input_new.eof() ){
		input_new >> en_new >> xs_new;
		if( input_new.eof() ) break;
		EN_new.push_back(en_new);
		bin_content = inter.Eval(en_new);
		XS_new.push_back( bin_content );
		lines_new++;
	}
	
	TGraph *g_new = new TGraph(lines_new);//, &EN_new[0], &XS_new[0]);
	for(int i=0; i<lines_new;++i){
		g_new->SetPoint(i, EN_new[i], XS_new[i]);
	}
	g_new->SetLineColor(kRed);
	g_new->SetLineWidth(2);
	g_new->Draw("L*");
	g_new->SetName("g_new");
	g_ENDF->GetXaxis()->SetLimits(1e-6, 1e-5);
	g_ENDF->SetMarkerStyle(20);
	g_ENDF->SetMarkerSize(3);
	g_new->SetMarkerColor(kRed);
	g_new->SetMarkerStyle(20);
	gPad->SetLogy(1);

}

[Wile_E_Coyote]

{
  TGraph *g10 = new TGraph(10);
  for (int i = 1; i <= 10; i++) g10->SetPoint(i - 1, i, i);
  g10->SetMarkerColor(kBlue); g10->SetLineColor(kBlue);
  TGraph *g2 = new TGraph(2);
  g2->SetPoint(0, 1., 1.); g2->SetPoint(1, 10., 10.);
  g2->SetMarkerColor(kRed); g2->SetLineColor(kRed);
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(2, 2);
  c->cd(1);
  g10->Draw("AL*"); g2->Draw("L*");
  c->cd(2);
  g10->Draw("AL*"); g2->Draw("L*"); gPad->SetLogx(1);
  c->cd(3);
  g10->Draw("AL*"); g2->Draw("L*"); gPad->SetLogy(1);
  c->cd(4);
  g10->Draw("AL*"); g2->Draw("L*"); gPad->SetLogx(1); gPad->SetLogy(1);
  c->cd(0);
}