/////////////////////////////////////////////////////
//
// Calcula o grafico de residuos em relacao a uma
// funcao f. Primeiro faz uma copia do grafico 
// original e depois altera os valores Y desse 
// grafico.
//
/////////////////////////////////////////////////////

TGraphErrors* residuo(TGraphErrors* g, TF1* f)
{
  int n = g->GetN();
  TGraphErrors *r = new TGraphErrors(*g);
  double *x = r->GetX();
  double *y = r->GetY();
  double *e = r->GetEY();
  for(int i = 0; i<n;i++)
  {
    y[i] = (y[i]-f(x[i]))/e[i];
    e[i] = 1;
  }  
  r->GetYaxis()->SetTitle("residuos");  
  return r;
}

/////////////////////////////////////////////////////
//
// Calcula o erro em uma funcao em um ponto X qualquer
// considera a matrix de covariancia do ajuste
//
/////////////////////////////////////////////////////

double erro(double x, TF1* f, TMatrixD COV)
{
  int n = f->GetNpar();
  double e = 0;
  for(int i = 0; i<n; i++)
    for(int j = 0; j<n; j++)
    {
      double d1 = f->GradientPar(i,&x);
      double d2 = f->GradientPar(j,&x);      
      e+=d1*d2*COV(i,j);
    }
  return sqrt(e);
}

void ajuste()
{  
  
  //
  // cria um grafico a partir de dados em um arquivo.
  // 
  TGraphErrors *g = new TGraphErrors("dados.txt","%lg %lg %lg %lg");
  g->Draw("AP");
  g->SetLineColor(2);
  g->SetMarkerColor(2);
  g->SetMarkerStyle(20); 
  g->GetXaxis()->SetTitle("T [K]");
  g->GetYaxis()->SetTitle("P [W]");
  g->SetTitle("Potencia em uma lampada de automovel");
  
  //
  // muda a escala para di-log na tela. A variavel
  // gPad corresponde sempre `a tela grafica ativa
  // no root
  //
  gPad->SetLogy();
  gPad->SetLogx();
   
  //
  // brinca com algumas propriedades dos eixos x e y
  // do grafico, tirando expoentes e adicionando mais
  // labels na escala em X
  //
    
  g->GetXaxis()->SetMoreLogLabels();
  g->GetXaxis()->SetNoExponent();
  g->GetYaxis()->SetNoExponent();
      
  // 
  // cria funcao e faz ajuste aos dados
  // na funcao de Fit, N significa para nao atrelar a funcao ao grafico
  // o que gera alguns inconvenientes praticos. R significa para usar
  // o intervalo de definicao da funcao para fazer o ajuste.
  //
  TF1 *f = new TF1("teoria","[0]*(x-300)^[1]+[2]*(x^[3]-300^[3])",400,3000);
  f->SetParameters(3e-4, 1.1, 5e-15, 4); 
  f->Draw("sameL");
  g->Fit(f,"NR");  
    
  cout <<"Chi2 = "<<f->GetChisquare()<<endl;
  cout <<"NDF =  "<<f->GetNDF()<<endl;
  
  // 
  // pega a matriz de covariancia. Por padrao, o root nao
  // guarda as matrizes de covariancia dos ajustes porque
  // elas podem ocupar um espaco de memoria grande. Assim,
  // usando a variavel gMinuit, que corresponde ao objeto
  // responsavem por operacionar os ajustes, podemos ter
  // acesso sempre `a ultima matriz de covariancia do
  // ultimo ajuste realizado e guarda-la, se necessario.
  //

  TMatrixD C(4,4);
  gMinuit->mnemat(C.GetMatrixArray(),C.GetNrows());
  C.Print();
  
  TLegend *L = new TLegend(0.125,0.75,0.43,0.89);
  L->AddEntry(g,"Dados experimentais","p");
  L->AddEntry(f,"Ajuste","l");
  L->Draw();
  
  //
  // desenha grafico de residuos
  //
  new TCanvas();
  TGraphErrors *res = residuo(g,f);
  res->Draw("AP"); 
  
  // 
  // faz um grafico de incertezas relativas no ajuste
  // efetuado, considerando  a matriz de covariancia
  //   
  int n = 50;
  TGraphErrors *t = new TGraphErrors(n);
  for(int i =0;i<n;i++)
  {
    double x = 400 +(3000-400)*i/n;
    double y = f(x);
    double e = erro(x,f,C);
    t->SetPoint(i,x,1);
    t->SetPointError(i,0,e/y);
  }
  
  new TCanvas(); 
  t->Draw("Al");

}