Hi,
I hope to quantify the area/integral under a TGraph and was just wondering if there is any function already implemented in ROOT somewhere?
or if anyone has done this could they offer some advice on the best way to go about this?
Thanks in advance,
Ben
Given that you don’t know anything about spacing between different x values, I’d recommend using the trapezoid rule.
en.wikipedia.org/wiki/Trapezoidal_rule
The only thing I’d check is that you need to make sure you get your points in order (e.g., either increasing or decreasing X value).
Charles
Hi,
there is no direct function to do it.
You can compute it from a function (like a TF1) obtained for example by fitting the TGraph or alternatively by computing directly the area by summing the trapezoids obtained from the (x,y) points of the TGraph.
You can also smooth your original TGraph in a new TGraph with more points, by using TGraphSmooth.
Best Regards
Lorenzo
One could easily add a TGraph::Integral(Int_t first=0, Int_t last=0,Option_t *option="") along the lines suggested by Charles.
However, what do you expect in the most frequent case illustrated in the picture in attachment? Should integral return the colored area?
Rene
[quote=“brun”]One could easily add a TGraph::Integral(Int_t first=0, Int_t last=0,Option_t *option="") along the lines suggested by Charles.
However, what do you expect in the most frequent case illustrated in the picture in attachment? Should integral return the colored area?
[/quote]
I would (personally) expect Integral() to do the same as the TH1::Integral (i.e., not what your picture illustrates).
Computing TGraph::Integral like in TH1 will work for graphs having points in increasing order only. What do you do for closed graphs like in the picture in attachment? Most of the request for TGraph::Integral were so far for this type of graph.
Rene
[quote=“brun”]Computing TGraph::Integral like in TH1 will work for graphs having points in increasing order only. What do you do for closed graphs like in the picture in attachment? Most of the request for TGraph::Integral were so far for this type of graph.
[/quote]
How easy is it to tell that the TGraph is “closed” or not? Can you give it two different options (i.e. a “TH1-like” integration option and a “assume-it-is-a-closed-tgraph” option". Of course, calculating the latter is not as easy.
Yes, one has to check
-if the graph has points in increasing order in all cases (if not -> sort
into a temp graph)
-if the request is a a closed graph (main use is to compute the integral
inside a TCutG). This is a bit more difficult.
The user could instruct integral about the above cases or specify the
type of integral that he/she requires.
If somebody is willing to contribute, let me know.
Rene
The following example shows a possible way to compute the surface of a closed TGraph:
Double_t graph_surface(TGraph *g, Int_t nb)
{
Double_t s = 0;
Double_t *x = g->GetX();
Double_t *y = g->GetY();
Int_t n = g->GetN();
Double_t xmin = x[0], xmax = x[0];
Double_t ymin = y[0], ymax = y[0];
Int_t i;
for (i=1; i<n; i++) {
if (xmin>x[i]) xmin = x[i];
if (xmax<x[i]) xmax = x[i];
if (ymin>y[i]) ymin = y[i];
if (ymax<y[i]) ymax = y[i];
}
Double_t dx = (xmax-xmin)/nb;
Double_t dy = (ymax-ymin)/nb;
Double_t dx2 = dx/2;
Double_t dy2 = dy/2;
Double_t ds = dx*dy;
Double_t xi, yi, xc, yc;
for (xi=xmin; xi<xmax; xi=xi+dx) {
for (yi=ymin; yi<ymax; yi=yi+dy) {
xc = xi+dx2;
yc = yi+dy2;
if (TMath::IsInside(xc, yc, n, x, y)) {
s = s + ds;
}
}
}
return s;
}
void integral(Int_t nb=50)
{
TGraph *g = new TGraph(10);
g->SetPoint(0,0.1278736,0.537155);
g->SetPoint(1,0.5502874,0.7707006);
g->SetPoint(2,0.4683908,0.5647558);
g->SetPoint(3,0.2571839,0.7346072);
g->SetPoint(4,0.2988506,0.3821656);
g->SetPoint(5,0.4784483,0.507431);
g->SetPoint(6,0.2183908,0.4946921);
g->SetPoint(7,0.4022989,0.552017);
g->SetPoint(8,0.2600575,0.5583864);
g->SetPoint(9,0.1278736,0.537155);
printf("surface = %g\n",graph_surface(g,nb));
}The parameter of the main macro is the presicion used to compute the surface. Example:
root [0] .x integral.C(100)
surface = 0.0397341Hi! I have to compare the integral values of the curves of which I attached a screenshot
Al curves have the same number of points. So can I simply call TGraph::Integral() in order to get the integral between 284 and 900 nm?
Many thanks
C
Hi! Any comment on my previous post? Thanks
C
If you want the integral under the graph curve I would create first a TF1 from the Graph and integrate the function, as in this example
double x[] = { 1,2,3,4,5,6,7,8,9,10};
double y[] = { 1,2,3,4,5,6,7,8,9,10};
TGraph g(10,x,y);
TF1 f1("f",[&](double *x, double *){ return g.Eval(x[0]); },1,10,0);
double integral = f1.Integral(1,10);Lorenzo
[quote=“moneta”]If you want the integral under the graph curve I would create first a TF1 from the Graph and integrate the function, as in this example
double x[] = { 1,2,3,4,5,6,7,8,9,10};
double y[] = { 1,2,3,4,5,6,7,8,9,10};
TGraph g(10,x,y);
TF1 f1("f",[&](double *x, double *){ return g.Eval(x[0]); },1,10,0);
double integral = f1.Integral(1,10);Lorenzo[/quote]
Dear Lorenzo,
I have difficulties understandin the syntax “TF1 f1(“f”,[&](double *x, double *){ return g.Eval(x[0]); },1,10,0);”
do basically I define a function with name “f”, a ROOT::Math::ParamFunctor f and then the min and max value in x and the npar?
Thanks in advance for clarifying
C
It is using a lambda function, a new features of C+±11.
See root.cern.ch/doc/master/classTF1.html#F5 and
en.cppreference.com/w/cpp/language/lambda for understanding what a lambda is
Best Regards
Lorenzo
[quote=“moneta”]It is using a lambda function, a new features of C+±11.
See root.cern.ch/doc/master/classTF1.html#F5 and
en.cppreference.com/w/cpp/language/lambda for understanding what a lambda is
Best Regards
Lorenzo[/quote]
Thanks for answering. I will check. However I already encountered a problem. Using the same syntaz as you did I am getting an error:
TF1::TF1(“f”,&{return g60_5->Eval(lambda[0]);},290,900,0) ;
any idea why?
Which version of ROOT are you using ? If it is version 5 this works only if you have compiled it with C+±11 support.
Lorenzo
[quote=“moneta”]Which version of ROOT are you using ? If it is version 5 this works only if you have compiled it with C+±11 support.
Lorenzo[/quote]
I use root_v5.34.30 whih has been installed at our machines. I am not sure whether it got compiled with that support. May I ask you how to compile with the suppert mentioned by you orhow to find out whether ot is or not?
OK, if you are using this old ROOT version then you don;t have support for these new features.
In that case you can create a free function or a C++ functor as described in
root.cern.ch/doc/master/classTF1.html#F5
struct MyFunction {
MyFunction(TGraph & g) : fGraph(&g) {}
double operator() (double *x, double *) { return fGraph->Eval(x[0]); }
TGraph * fGraph;
};Then at the prompt you can do
.L MyFunction.h
double x[] = { 1,2,3,4,5,6,7,8,9,10};
double y[] = { 1,2,3,4,5,6,7,8,9,10};
TGraph g(10,x,y);
MyFunction f(g);
TF1 f1("f1",f,1,10,0,"MyFunction");Lorenzo
Great! Many thanks! I will tell you whether it works!
Dear Lorenzo
I tried your code and I was still getting an error:
Error: Can’t call MyFunction::MyFunction(g60_5) in current scope new_lauyout.C:7110:
Possible candidates are…
(in MyFunction)
MyFunction.h 2:1 0 public: MyFunction MyFunction::MyFunction(TGraph& g);
*** Interpreter error recovered ***
As for the .h file, I just opened a file named it MyFunction.h and pasted the code from you below without changing anything
so since I define the TGraph like this
TGraph *g60_5=new TGraph();,
I changed the code to
MyFunction f(*g60_5) (is this actually correct?)
That solved, I got another error:
Error: Can’t call TF1::TF1(“f1”,f,1,10,0,“MyFunction”) in current scope (tmpfile):1:
So I changed it to (“f1”,&f,1,10,0,“MyFunction”)
and in this way I did not get any error messages and could compute the integral. But I am not sure whether I was am doing the right thing.
[quote=“moneta”]OK, if you are using this old ROOT version then you don;t have support for these new features.
In that case you can create a free function or a C++ functor as described in
root.cern.ch/doc/master/classTF1.html#F5
struct MyFunction {
MyFunction(TGraph & g) : fGraph(&g) {}
double operator() (double *x, double *) { return fGraph->Eval(x[0]); }
TGraph * fGraph;
};Then at the prompt you can do
.L MyFunction.h
double x[] = { 1,2,3,4,5,6,7,8,9,10};
double y[] = { 1,2,3,4,5,6,7,8,9,10};
TGraph g(10,x,y);
MyFunction f(g);
TF1 f1("f1",f,1,10,0,"MyFunction");Lorenzo[/quote]