/*
 * Per risolvere l'esercizio, bisogna dare le seguenit istruzioni in root:
 * .L fitRisonanza.C (carica le funzioni)
 * riordina("datiRisonanza.txt","datiOrdinati.txt")
 * fitRisonanza("datiOrdinati.txt")
 */


double BWfun(double *x, double*p) {
  // la vriabile x[0] e' la pulazione omega
  // i parametri sono p[0]=F, p[1]=omega0, p[2]=alpha
  double d1 = x[0]*x[0]-p[1]*p[1];
  double d2 = x[0]*p[2];
  return p[0]/sqrt(d1*d1+d2*d2);
}

TF1 *risonanza=new TF1("risonanza",BWfun,0.,3.,3);

void riordina(char* filein, char* fileout) {
  double *x=new double[100]; // creiamo un vettore abbastanza lungo
  double *y=new double[100];
  // lettura dei dati
  int N=0;
  ifstream in(filein);
  while (1) {
    in >> x[N] >> y[N];
    if (in.eof()||in.bad()) break;
    N++;
  }
  // ordinamento (simple)
  for (int i=0; i<N-1; i++) {
    for (int j=i+1; j<N; j++) {
      if ( x[i]>x[j] ) {
	double dum=x[i]; x[i]=x[j]; x[j]=dum; // si scambiano le x
	dum=y[i]; y[i]=y[j]; y[j]=dum; // ricordarsi di scambiare anche le y
      }
    }
  }
  // scrittura dei dati
  int i=0;
  ofstream out(fileout);
  for ( int i=0; i<N; i++) out <<  x[i] << " " <<  y[i] << endl;
  out.close();
  return;
}


double fitRisonanza(char* filename) {
  TGraphErrors *g = new TGraphErrors;
  ifstream f(filename);
  int i=0;
  double frequenza, ampiezza;
  while (1) {
    f >> frequenza >> ampiezza;
    if ( f.eof() || f.bad() ) break;
    g->SetPoint(i,2*TMath::Pi()*frequenza,ampiezza);
    g->SetPointError(i,0.,0.01);
    i++;
  }
  risonanza->SetParameter(0,0.2);
  risonanza->SetParameter(1,11.);
  risonanza->SetParameter(2,0.01);
  g->Fit(risonanza,"Q");
  g->Draw("AL");
  double F = risonanza->GetParameter(0);
  double eF= risonanza->GetParError(0);
  cout << "F=      " << F << "+-" << eF << endl;
  double w = risonanza->GetParameter(1);
  double ew= risonanza->GetParError(1);
  cout << "omega0= " << w << "+-" << ew << endl;
  double a = risonanza->GetParameter(2);
  double ea= risonanza->GetParError(2);
  cout << "alpha=  " << a << "+-" << ea << endl;
}