void ricorsivo_new()
{
//   example of macro to read data from an ascii file and
//   create a root file with an histogram and an ntuple.
   gROOT->Reset();
//#include "Riostream.h"

Int_t i,j,k;

const Int_t NYEARS=57; // numero di anni
const Int_t NSTEPS=100; // numero di passi teorici da calcolare per ogni anno

Int_t year[NYEARS];
Double_t z1[NYEARS],z2[NYEARS],errorz1[NYEARS],errorz2[NYEARS];
Double_t teo1[NYEARS*NSTEPS],teo2[NYEARS*NSTEPS];
Double_t teo1_y[NYEARS],teo2_y[NYEARS];

ifstream in;
// we assume a file basic.dat in the current directory
// this file has 3 columns of float data
in.open("canada.txt");
Int_t nlines = 0;
while (1) {
  in >> year[nlines] >> z1[nlines] >> z2[nlines];
  if (!in.good()) break;
  nlines++;
}
in.close();

// The errors on z values
 for (Int_t i=0;i<NYEARS; i++) {
   errorz1[i]=sqrt(z1[i]);
   errorz2[i]=sqrt(z2[i]);
 }

// visualizzo dati sperimentali:
Double_t maximz1=1.2*maxim(z1,nlines);
Double_t maximz2=1.2*maxim(z2,nlines);
TH2F *h1 = new TH2F("h1","xy",57,1856,1913,100,0,maximz1);
TH2F *h2 = new TH2F("h2","xz",57,1856,1913,100,0,maximz2);
for (i=0;i<NYEARS; i++) {
  h1->Fill(year[i],z1[i]);
  h2->Fill(year[i],z2[i]);
}
c1 = new TCanvas("c1"," ",200,10,600,480);
c1->Divide(1,2);

h1->SetMarkerStyle(21);
h1->SetMarkerSize(0.7);

h2->SetMarkerStyle(22);
h2->SetMarkerSize(0.7);

c1_1->cd();
h1->Draw();
c1_2->cd();
h2->Draw();

// calcolo teorico:

 Double_t alfa=0.7;
 Double_t gamma=0.5;

 Double_t minimo[2];
 Double_t tempmin;
 for(Int_t kk=0;kk<300;kk++){
   for(Int_t ll=0;ll<100;ll++){

     // Double_t par[2]={0.000029, 0.93 };
     Double_t par[2];
     par[0]=0.000001+0.000001*kk;
     par[1]=0.01+0.01*ll;

     teo1[0]=z1[0];
     teo2[0]=z2[0];
 
     for (j=1;j<NYEARS*NSTEPS; j++) { 
       teo1[j]=teo1[j-1]+(alfa*teo1[j-1]-par[0]*teo1[j-1]*teo2[j-1])/NSTEPS;
       teo2[j]=teo2[j-1]+(par[1]*par[0]*teo1[j-1]*teo2[j-1]-gamma*teo2[j-1])/NSTEPS;
     }
     // a noi serve un vettore con la media per ogni anno
     Double_t media1=0;
     Double_t media2=0;
     for (i=0;i<NYEARS; i++) {
       for (k=0;k<NSTEPS; k++) {
	 media1 += teo1[k+NSTEPS*i];
	 media2 += teo2[k+NSTEPS*i];
       }
       media1 = media1/NSTEPS;
       media2 = media2/NSTEPS;
       teo1_y[i]=media1;
       teo2_y[i]=media2;
       //     teo1_y[i]=teo1[NSTEPS*i];
       //     teo2_y[i]=teo2[NSTEPS*i];
     }

     //confronto con i dati:
     Double_t chisq = 0;
     Double_t delta1,delta2;
     for (i=0;i<NYEARS; i++) {
       delta1  = (z1[i]-teo1_y[i]);
       delta2  = (z2[i]-teo2_y[i]);
       chisq += delta1*delta1 + delta2*delta2;
     }
     cout << "-------------------------------------------------------------" << endl;
     cout << "par[0]="<< par[0] << ", par[1]=" << par[1] << ", chisq=" << chisq << endl;
     cout << "-------------------------------------------------------------" << endl;
     //aggiorna il minimo:
     if (kk==1 && ll==1) {
       tempmin=chisq;
       minimo[0]=par[0];
       minimo[1]=par[1];
     } else {
       if (chisq<tempmin) {
	 tempmin=chisq;
	 minimo[0]=par[0];
	 minimo[1]=par[1]; 
       }
     }
   }
 }

 cout << "-------------------------------------------------------------" << endl;
 cout << "MINIMO:" << endl;
 cout << "par[0]="<< minimo[0] << ", par[1]=" << minimo[1] << ", chisq=" << tempmin << endl;
 cout << "-------------------------------------------------------------" << endl;
 
 //ricalcoliamo i vettori solo al minimo:
 par[0]=minimo[0];
 par[1]=minimo[1];

 teo1[0]=z1[0];
 teo2[0]=z2[0];
 
 for (j=1;j<NYEARS*NSTEPS; j++) { 
   teo1[j]=teo1[j-1]+(alfa*teo1[j-1]-par[0]*teo1[j-1]*teo2[j-1])/NSTEPS;
   teo2[j]=teo2[j-1]+(par[1]*par[0]*teo1[j-1]*teo2[j-1]-gamma*teo2[j-1])/NSTEPS;
 }
 // a noi serve un vettore con la media per ogni anno
 Double_t media1=0;
 Double_t media2=0;
 for (i=0;i<NYEARS; i++) {
   for (k=0;k<NSTEPS; k++) {
     media1 += teo1[k+NSTEPS*i];
     media2 += teo2[k+NSTEPS*i];
   }
   media1 = media1/NSTEPS;
   media2 = media2/NSTEPS;
   teo1_y[i]=media1;
   teo2_y[i]=media2;
   //     teo1_y[i]=teo1[NSTEPS*i];
   //     teo2_y[i]=teo2[NSTEPS*i];
 }

     
// visualizzo dati teorici:
Double_t minimt1=minim(teo1_y,NYEARS);
Double_t minimt2=minim(teo2_y,NYEARS);
Double_t maximt1=1.2*maxim(teo1_y,NYEARS);
Double_t maximt2=1.2*maxim(teo2_y,NYEARS);
TH2F *g1 = new TH2F("g1","xy",57,1856,1913,100,minimt1,maximt1);
TH2F *g2 = new TH2F("g2","xz",57,1856,1913,100,minimt2,maximt2);

for (i=0;i<NYEARS; i++) {
  g1->Fill(year[i],teo1_y[i]);
  g2->Fill(year[i],teo2_y[i]);
}

c2 = new TCanvas("c2"," ",800,10,600,480);
c2->Divide(1,2);

g1->SetMarkerStyle(21);
g1->SetMarkerSize(0.7);

g2->SetMarkerStyle(22);
g2->SetMarkerSize(0.7);

c2_1->cd();
g1->Draw();
c2_2->cd();
g2->Draw();


}

Double_t maxim(Double_t v[], Int_t nlines) {
  Double_t temp=0;
  for (Int_t i=0; i<nlines; i++) {
    if (v[i]>temp) {
      temp=v[i];
    }
  }
  return temp;
}
Double_t minim(Double_t v[], Int_t nlines) {
  Double_t temp=0;
  for (Int_t i=0; i<nlines; i++) {
    if (v[i]<temp) {
      temp=v[i];
    }
  }
  return temp;
}
Double_t integ(Double_t v[], Int_t nlines) {
  Double_t sum=0;
  for (Int_t i=0; i<nlines; i++) {
    sum += v[i];
  }
  cout << nlines << "  Integrale calcolato: " << sum << endl;

  return sum;
}
Double_t exp_safe(Double_t x) {
  if (x<-30.) {
    return 0;
  } elseif (x>30) {
    return 200000;
  } else {
    return TMath::Exp(x);
  }
}
Double_t func(Double_t a, Double_t b, Double_t c, Int_t j, Double_t v[]) {
  /* Questa e' la funzione da usare per l'equazione integrale. */
  /* Qui "a" e' il logaritmo del valore INIZIALE della funzione. */
  Double_t temp = a+b*(c*j-integ(v,j));
  return exp_safe(temp); //quando l'argomento e' troppo negativo al posto dell'esponenziale impone zero
}
Double_t func2(Double_t a, Double_t b, Double_t c, Double_t d) {
  /* Questa e' la funzione da usare per l'equazione differenziale. */
  /* Qui "a" e' il logaritmo del valore PRECEDENTE della funzione stessa, "d" il valore precedente dell'altra funzione. */
  Double_t temp=a+b*(c-d);
  return exp_safe(temp); //quando l'argomento e' troppo negativo al posto dell'esponenziale impone zero
}