/*
  Some utility functions used by FF_RAA_brick.C

  Author: M. van Leeuwen, Utrecht University
*/

TGraph *divide_graph(TGraph *gr, TF1* f1) {
  Int_t n=gr->GetN();
  TGraph *rat_gr=new TGraph(n);
  for (Int_t i=0; i<n; i++) {
    Double_t x,y;
    gr->GetPoint(i,x,y);
    Double_t y_fun=f1->Eval(x);
    if ( y_fun != 0) {
      rat_gr->SetPoint(i,x,y/y_fun);
    }
  }
  return rat_gr;
}

TGraph *divide_graphs_int(TGraph *num, TGraph *denom) {
  Int_t j = 0;
  Double_t x1,y1,x2,y2, y_int;
  denom->GetPoint(j,x1,y1);
  j++;
  denom->GetPoint(j,x2,y2);
  j++;
  TGraph *gr = new TGraph();
  for (Int_t i = 0; i < num->GetN(); i++) {
    Double_t x,y;
    num->GetPoint(i,x,y);
    while (x2 < x && j < denom->GetN()) {
      x1 = x2; y1 = y2;
      denom->GetPoint(j,x2,y2);
      j++;
    }
    if ((x2-x1)==0)
      cout << "Warning in divide_graphs_int: x2-x1=0, i=" << i << " j=" << j << endl;
    if (x <= x1+(x2-x1)*1.05 && x > x1) {  // allow small tolerance for last point
      y_int = (x2-x)/(x2-x1)*y1 + (x-x1)/(x2-x1)*y2;
      //cout << "x " << x << " x1 " << x1 << " x2 " << x2 << " y " << y << " y_int " << y_int << endl;
      if (y_int == 0)
	cout << "Warning in divide_graphs_int: y_int = 0" << endl;
      else 
	gr->SetPoint(gr->GetN(),x,y/y_int);
    }
  }
  return gr;
}