#include "TSpline.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TFrame.h"
   
Double_t Integral(TGraph *gr, Double_t a, Double_t b, Option_t *option="",Double_t epsilon=1e-6);
void gint() {
   TCanvas *c1 = new TCanvas("c1","gerrors2",200,10,700,500);
   c1->SetFillColor(42);
   c1->SetGrid();

      // draw a frame to define the range
   TH1F *hr = c1->DrawFrame(-0.4,0,1.2,12);
   hr->SetXTitle("X title");
   hr->SetYTitle("Y title");
   c1->GetFrame()->SetFillColor(21);
   c1->GetFrame()->SetBorderSize(12);
   Int_t n1 = 10;
   Double_t x1[]  = {-0.22, 0.05, 0.25, 0.35, 0.5, 0.61,0.7,0.85,0.89,0.95};
   Double_t y1[]  = {1,2.9,5.6,7.4,9,9.6,8.7,6.3,4.5,1};
   Double_t ex1[] = {.05,.1,.07,.07,.04,.05,.06,.07,.08,.05};
   Double_t ey1[] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8};
   TGraphErrors *gr1 = new TGraphErrors(n1,x1,y1,ex1,ey1);
   gr1->SetMarkerColor(kBlue);
   gr1->SetMarkerStyle(21);
   gr1->Draw("LP");

      // create second graph
   Int_t n2 = 10;
   Float_t x2[]  = {-0.28, 0.005, 0.19, 0.29, 0.45, 0.56,0.65,0.80,0.90,1.01};
   Float_t y2[]  = {0.82,3.86,7,9,10,10.55,9.64,7.26,5.42,2};
   Float_t ex2[] = {.04,.12,.08,.06,.05,.04,.07,.06,.08,.04};
   Float_t ey2[] = {.6,.8,.7,.4,.3,.3,.4,.5,.6,.7};
   TGraphErrors *gr2 = new TGraphErrors(n2,x2,y2,ex2,ey2);
   gr2->SetMarkerColor(kRed);
   gr2->SetMarkerStyle(20);
   gr2->Draw("LP");
   printf("gr1 integral = %g\n",Integral(gr1,0,0.6,"3"));
   printf("gr2 integral = %g\n",Integral(gr2,0,0.6,"5"));
}
//______________________________________________________________________________
Double_t Integral(TGraph *gr, Double_t a, Double_t b, Option_t *option, Double_t epsilon)
{
   //Compute the integral of TGraph *gr from a to b.
   //if option "5 is specified a quintic spline is used to approximate
   //the graph, otherwise a cubic spline is used.
   //The method used is described in TF1::Integral
   //This function assumes that 
   //  -the points in the graph are sorted along X
   //  -there is only one value of Y for a given value of X in [a,b].
   //if ga = gr(a) and gb=gr(b), the value returned excludes the area
   //outside the line connecting (a,ga) and (b,gb) by default.
   //To also include this area, specify the option "out"

   const Double_t kHF = 0.5;
   const Double_t kCST = 5./1000;

   Double_t x[12] = { 0.96028985649753623,  0.79666647741362674,
                      0.52553240991632899,  0.18343464249564980,
                      0.98940093499164993,  0.94457502307323258,
                      0.86563120238783174,  0.75540440835500303,
                      0.61787624440264375,  0.45801677765722739,
                      0.28160355077925891,  0.09501250983763744};

   Double_t w[12] = { 0.10122853629037626,  0.22238103445337447,
                      0.31370664587788729,  0.36268378337836198,
                      0.02715245941175409,  0.06225352393864789,
                      0.09515851168249278,  0.12462897125553387,
                      0.14959598881657673,  0.16915651939500254,
                      0.18260341504492359,  0.18945061045506850};

   Double_t h, aconst, bb, aa, c1, c2, u, s8, s16, f1, f2;
   Double_t xx;
   Int_t i;
   h = 0;
   if (b == a) return h;
   //Evaluate a spline smoothing this graph
   TSpline *spline = 0;
   TString opt = option;
   opt.ToLower();
   if (opt.Contains("5")) spline= new TSpline3("integral",gr);
   else                   spline= new TSpline5("integral",gr);
   Bool_t absValue = kFALSE;
   if (opt.Contains("abs")) absValue = kTRUE;
   
   //Evaluate the integral
   Double_t ga = gr->Eval(a,spline);
   Double_t gb = gr->Eval(b,spline);
   aconst = kCST/TMath::Abs(b-a);
   bb = a;
CASE1:
   aa = bb;
   bb = b;
CASE2:
   c1 = kHF*(bb+aa);
   c2 = kHF*(bb-aa);
   s8 = 0;
   for (i=0;i<4;i++) {
      u     = c2*x[i];
      xx = c1+u;
      f1    = gr->Eval(xx,spline);
      if (absValue) f1 = TMath::Abs(f1);
      xx = c1-u;
      f2    = gr->Eval(xx,spline);
      if (absValue) f2 = TMath::Abs(f2);
      s8   += w[i]*(f1 + f2);
   }
   s16 = 0;
   for (i=4;i<12;i++) {
      u     = c2*x[i];
      xx = c1+u;
      f1    = gr->Eval(xx,spline);
      if (absValue) f1 = TMath::Abs(f1);
      xx = c1-u;
      f2    = gr->Eval(xx,spline);
      if (absValue) f2 = TMath::Abs(f2);
      s16  += w[i]*(f1 + f2);
   }
   s16 = c2*s16;
   if (TMath::Abs(s16-c2*s8) <= epsilon*(1. + TMath::Abs(s16))) {
      h += s16;
      if(bb != b) goto CASE1;
   } else {
      bb = c1;
      if(1. + aconst*TMath::Abs(c2) != 1) goto CASE2;
      h = s8;
   }
   delete spline;
   if (!opt.Contains("out")) return h;
   else                      return h - 0.5*(ga+gb)/(b-a);
}