#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <signal.h>
#include <unistd.h>
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <fstream>

#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TH1.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TRandom.h>
#include <TSystem.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TCanvasImp.h>
#include <TH1.h>
#include <TFile.h>
#include <TInterpreter.h>
#include <TBenchmark.h>
#include <TNtuple.h>
#include <TMatrixD.h>
#endif

/////////////////////////////////////////////////////////////////////////
struct StatResult {
  Double_t mean;
  Double_t emean;
  Double_t rms;
  Double_t chi2;
  Double_t chi2red;
};
/////////////////////////////////////////////////////////////////////////
Bool_t gPrintResult=1;
void SetPrintResult(Bool_t pr=1) {
  // settare gPrintResult=0 per non stampare i risultati delle funzioni
  gPrintResult=pr;
}
/////////////////////////////////////////////////////////////////////////
StatResult gMeanResult;
StatResult* GetMeanResult() {return &gMeanResult;}
Double_t GetMean(Double_t *x, Int_t N) {
  // calcola il valore medio di N valori x[] e stampa le dev. standard
  // scrive i valori calcolati in gMeanResult 
  Double_t sx=0,xm=0,ss2=0,sqm=0,em=0;
  for (Int_t i=0; i<N; i++)
    sx += x[i];
  xm=sx/N;
  for (Int_t i=0; i<N; i++)
    ss2 += TMath::Power((x[i]-xm),2);
  sqm = TMath::Sqrt(ss2/(N-1));
  em = sqm/TMath::Sqrt(double(N));

  gMeanResult.mean=xm;
  gMeanResult.emean=em;
  gMeanResult.rms=sqm;
  
  if (gPrintResult)
    printf("mean over %d values: %f +- %f\nsqm = %f\n",N,xm,em,sqm);

  return xm;
}
/////////////////////////////////////////////////////////////////////////
Double_t GetMeanW(Double_t *x, Double_t *ex, Int_t N) {
  // ritorna il valore medio pesato di N valori x[] e stampa le dev. standard
  // scrive i valori calcolati in gMeanResult 
  Double_t sx=0,xm=0,ss2=0,em=0,chi2=0;
  for (Int_t i=0; i<N; i++) {
    sx += x[i]/ex[i]/ex[i];
    ss2 += 1/ex[i]/ex[i];
  }
  xm=sx/ss2;
  em = TMath::Sqrt(1/ss2);

  
  for (Int_t i=0; i<N; i++) {
    chi2 += (x[i]-xm)*(x[i]-xm)/ex[i]/ex[i];
  }
  
  gMeanResult.mean=xm;
  gMeanResult.emean=em;
  gMeanResult.chi2=chi2;
  gMeanResult.chi2red=chi2/(N-1);

  if (gPrintResult)
    printf("weighted mean over %d values: %f +- %f\nreduced chi2: %f\n",N,xm,em,chi2/(N-1));

  return xm;
}
/////////////////////////////////////////////////////////////////////////
Double_t Extrapolate(Double_t x,TGraphErrors *g, TF1 *f) {
  // valuta la sigma(y) del valore medio stimato dalla funzione di fit
  // il coefficiente di correlazione era stato messo in f->GetChisquare()
  Double_t y=f->Eval(x);
  Double_t sc=f->GetParError(0);
  Double_t sm=f->GetParError(1);
  Double_t rho=f->GetChisquare();
  Double_t sy2 = x*x*sm*sm + sc*sc + 2*x*rho*sm*sc;
  Double_t sy = TMath::Sqrt(sy2);
  if (gPrintResult)
    printf("Estimated y_mean( %f ): %f +- %f\n",x,y,sy);

  // valuta l'sqm previsto per la singola misura
  Int_t i=0;
  // prende come ey quello piu' vicino
  while (g->GetX()[i]<x && i<(g->GetN()-1)) i++;
  Double_t ey = g->GetErrorY(i);
  Double_t et2 = ey*ey + sy2;
  Double_t et = TMath::Sqrt(et2);
  if (gPrintResult)
    printf("Estimated y_measured( %f ): %f +- %f\n",x,y,et);

  return y;
}
/////////////////////////////////////////////////////////////////////////
TMatrixD *GetLFCovMat(TGraphErrors *g) {
  // calcola la matrice di covarianza per il fit lineare
  Int_t N=g->GetN();
  Double_t el[4]={0,0,0,0};

  for (Int_t i=0; i<N; i++) {
    el[0] += 1./TMath::Power(g->GetErrorY(i),2);
    el[1] += g->GetX()[i]/TMath::Power(g->GetErrorY(i),2);
    el[2] += g->GetX()[i]/TMath::Power(g->GetErrorY(i),2);
    el[3] += g->GetX()[i]*g->GetX()[i]/TMath::Power(g->GetErrorY(i),2);
  }
  TMatrixD *mi = new TMatrixD(2,2,el);
  TMatrixD *m = new TMatrixD((mi->Invert()));
  delete mi;
  return m;
}
/////////////////////////////////////////////////////////////////////////
TF1 *fgp = new TF1("gausprova","gaus",-10,10);

Double_t GetProb(Double_t dist, Double_t sigma) {
  fgp->SetParameters(3.98942280401432814e-01,0,1);
  Double_t prob=2.*fgp->Integral(TMath::Abs(dist/sigma),10);
  return prob;
}
/////////////////////////////////////////////////////////////////////////
Double_t GetChi2(TGraphErrors *g, TF1 *f1) {
  // ritorna il chiquadro del grafico con errori g associato alla
  // funzione di fit f1
  Double_t chi2=0;
  Int_t npar = f1->GetNpar();
  Int_t N=g->GetN();
  if (N <= npar) {
    printf("Error: N points < npar! Cannot proceed\n");
    return -1;
  }
  for (Int_t i=0; i<N; i++) {
    Double_t delta = g->GetY()[i] - f1->Eval(g->GetX()[i]);
    delta /= g->GetErrorY(i);
    chi2 += delta*delta;
  }

  return chi2;
}
/////////////////////////////////////////////////////////////////////////
TGraphErrors* GetResidues(TGraphErrors *g, TF1 *f1) {
  // ritorna il grafico dei residui associato
  // al grafico g e alla funzione di fit f1

  if (!g || !f1) return NULL;

  Int_t N=g->GetN();

  TGraphErrors *r = new TGraphErrors(N,g->GetX(),g->GetY(),g->GetEX(),g->GetEY());
  r->SetMarkerStyle(g->GetMarkerStyle());
  r->SetMarkerColor(g->GetMarkerColor());
  r->SetMarkerSize(g->GetMarkerSize());
  r->SetLineColor(g->GetLineColor());
  r->SetLineWidth(g->GetLineWidth());

  Char_t ng[200];
  strcpy(ng,g->GetName());
  strcat(ng,"_res");
  r->SetName(ng);

  for (Int_t i=0; i<N; i++) {
    Double_t delta = r->GetY()[i] - f1->Eval(r->GetX()[i]);
    r->SetPoint(i,r->GetX()[i],delta);
  }
  
  return r;
}
/////////////////////////////////////////////////////////////////////////
TF1* LinFitW(TGraphErrors *g) {
  // interpolazione lineare pesata
  // ritorna la funzione di fit (retta) o NULL in caso di errore
  // il coeff. di correlazione e' messo in Chisquare

  TMatrixD cm(*(GetLFCovMat(g)));
  Int_t N=g->GetN();
  Double_t el[2]={0,0};
  for (Int_t i=0; i<N; i++) {
    el[0] += g->GetY()[i]/TMath::Power(g->GetErrorY(i),2);
    el[1] += g->GetX()[i]*g->GetY()[i]/TMath::Power(g->GetErrorY(i),2);
  }
  Double_t c = el[0]*cm(0,0) + el[1]*cm(0,1);
  Double_t m = el[0]*cm(1,0) + el[1]*cm(1,1);
  Double_t ec = TMath::Sqrt(cm(0,0));
  Double_t em = TMath::Sqrt(cm(1,1));
  Double_t corr = cm(1,0)/ec/em;

  if (gPrintResult) {
    printf("\nLinear Fit using minimum ChiSquare\n\n");
    printf("   y = mx + c  (%d points):\n   m = %g +- %g\n   c = %g +- %g\n",N,m,em,c,ec );
    printf("   correl. = %f\n",corr); 
  }

  TF1 *p=new TF1("intlinw","pol1",g->GetX()[0], g->GetX()[N-1]);
  p->SetParameter(0,c);
  p->SetParError(0,ec);
  p->SetParameter(1,m);
  p->SetParError(1,em);
  p->SetChisquare(corr);
  
  if (gPrintResult)
    printf("   chi2 ridotto = %f\n\n",GetChi2(g,p)/(N-2)); 

  return p;
}
/////////////////////////////////////////////////////////////////////////
TF1* LinFit(TGraphErrors *g, Bool_t seterrors=1, Int_t ifrom=0, Int_t ito=0) {
  // interpolazione lineare con le formulette classiche
  //    ritorna la funzione di fit (retta) o NULL in caso di errore
  //    il coeff. di correlazione e' messo in Chisquare
  // - se seterrors=1 setta tutti gli errori su y alla sigma_y comune
  // - fitta i punti da ifrom a ito (compresi)

  Int_t N=g->GetN();  
  if (!ito || ito==N) ito=N;
  else ito++; // ito compreso...

  N=ito-ifrom;

  if (N<3) {
    printf("Error: number of points < 3! Cannot proceed\n");
    return NULL;
  }

  Double_t xm=0, ym=0, x2m=0, xym=0, sx=0, sy=0, sxy=0, sx2=0, varx=0, covxy=0;
  for (Int_t i=ifrom; i<ito; i++) {
    sx += g->GetX()[i];
    sy += g->GetY()[i];
    sxy += g->GetX()[i]*g->GetY()[i];
    sx2 += g->GetX()[i]*g->GetX()[i];
  }
  xm = sx/N;
  ym = sy/N;
  xym = sxy/N;
  x2m = sx2/N;
  varx = x2m - xm*xm;
  if (varx==0) {
    printf("Error: Var(x) = 0! Cannot proceed\n");
    return NULL;
  }
  covxy = xym - xm*ym;
  // y = mx + c
  Double_t m,em,c,ec,ey,corr;
  m = covxy/varx;
  c = ym - m*xm;
  ey=0;
  for (Int_t i=ifrom; i<ito; i++)
    ey += TMath::Power( (g->GetY()[i] - m*g->GetX()[i] - c) ,2);
  ey /= (N-2);
  ey = TMath::Sqrt(ey);
  em = ey/TMath::Sqrt(N*varx);
  ec = em*TMath::Sqrt(x2m);
  corr = -xm/TMath::Sqrt(x2m);

  if (seterrors) {
    for (Int_t i=ifrom; i<ito; i++)
      g->SetPointError(i,0,ey);
  }

  if (gPrintResult) {
    printf("\n\nLinear interpolation y = mx + c  (%d points):\n   m = %g +- %g\n   c = %g +- %g\n",N,m,em,c,ec );
    printf("   sigma_y = %f\n   correl. = %f\n",ey,corr); 
  }

  TF1 *p=new TF1("intlin","pol1",g->GetX()[ifrom], g->GetX()[ito-1]);
  p->SetParameter(0,c);
  p->SetParError(0,ec);
  p->SetParameter(1,m);
  p->SetParError(1,em);
  p->SetChisquare(corr);

  return p;
}
/////////////////////////////////////////////////////////////////////////
TF1* LinFitLor(TGraphErrors *g, Int_t ifrom=0, Int_t ito=0) {
  // interpolazione lineare pesata con le formulette classiche di Loreti
  //    ritorna la funzione di fit (retta) o NULL in caso di errore
  //    il coeff. di correlazione e' messo in Chisquare
  // - se seterrors=1 setta tutti gli errori su y alla sigma_y comune
  // - fitta i punti da ifrom a ito (compresi)

  Int_t N=g->GetN();  
  if (!ito || ito==N) ito=N;
  else ito++; // ito compreso...

  N=ito-ifrom;

  if (N<3) {
    printf("Error: number of points < 3! Cannot proceed\n");
    return NULL;
  }

  Double_t xm=0, ym=0, x2m=0, xym=0, sx=0, sy=0, sxy=0, sx2=0, delta=0, covxy=0, sey=0;
  for (Int_t i=ifrom; i<ito; i++) {
    sx += g->GetX()[i]/g->GetEY()[i]/g->GetEY()[i];
    sy += g->GetY()[i]/g->GetEY()[i]/g->GetEY()[i];
    sxy += g->GetX()[i]*g->GetY()[i]/g->GetEY()[i]/g->GetEY()[i];
    sx2 += g->GetX()[i]*g->GetX()[i]/g->GetEY()[i]/g->GetEY()[i];
    sey += 1./g->GetEY()[i]/g->GetEY()[i];
  }
  //xm = sx/N;
  //ym = sy/N;
  //xym = sxy/N;
  //x2m = sx2/N;
  delta = sey*sx2 - sx*sx;
  if (delta==0) {
    printf("Error: Delta = 0! Cannot proceed\n");
    return NULL;
  }
  //covxy = xym - xm*ym;
  // y = mx + c
  Double_t m,em,c,ec,ey,corr;
  c = (sx2*sy-sx*sxy)/delta;
  m = (sey*sxy-sx*sy)/delta;
  //ey=0;

  em = TMath::Sqrt(sey/delta);
  ec = TMath::Sqrt(sx2/delta);
  //corr = -xm/TMath::Sqrt(x2m);


  if (gPrintResult) {
    printf("\n\nWeited linear interpolation (Loreti style) y = mx + c  (%d points):\n   m = %g +- %g\n   c = %g +- %g\n",N,m,em,c,ec );
    //printf("   sigma_y = %f\n   correl. = %f\n",ey,corr); 
  }

  TF1 *p=new TF1("intlin","pol1",g->GetX()[ifrom], g->GetX()[ito-1]);
  p->SetParameter(0,c);
  p->SetParError(0,ec);
  p->SetParameter(1,m);
  p->SetParError(1,em);
  //p->SetChisquare(corr);

  return p;
}
/////////////////////////////////////////////////////////////////////////
void PCheckName(char* s) {
  // check for a good root name
  for (Int_t i=0; i<int(strlen(s)); i++) {
    if (s[i]==95) continue; // '_'
    if (s[i]<48) s[i]='_';
    if (s[i]>57 && s[i]<65) s[i]='_';
    if (s[i]>90 && s[i]<97) s[i]='_';
    if (s[i]>122) s[i]='_';
  }
  if (s[0]<65) s[0]='_'; // start with a number

}
//////////////////////////////////////////////////////////////////
TGraphErrors* LoadGraphErrors(Char_t *fname, Char_t *graphname="") { 
  // legge un grafico con errori da file tipo testo nel formato:
  //
  // # commenti
  // x0   y0   ex0   ey0
  // ..   ..   ..    ..
  // 
  
  TGraphErrors *g = NULL;
  
  ifstream infile;
  infile.open(fname);
  
  if (infile.fail()) {
    printf("file not found!\n");
    return NULL;
  }
  
  Char_t titolo[100];
  strcpy(titolo,"");
  
  Char_t hn[256];
  
  if (!strlen(graphname)) {
    char *pio = strrchr(fname,'/');
    if (pio) strcpy(hn,++pio);
    else strcpy(hn,fname);
    Int_t pospt = strcspn(hn,".");
    if (pospt>90) pospt=10;
    hn[pospt] = 0;
    hn[pospt+1] = 0;
  }  
  else 
    strcpy(hn,graphname);
  
  PCheckName(hn);
  
  Float_t x[8192],y[8192],ex[8192],ey[8192];
  Float_t xx,yy,exx,eyy;
  Char_t st[100];
  Char_t stmp[30];
  Int_t n=0;

  Int_t lcol = 2;
  Int_t lwid = 1;

  Int_t mcol = 4;
  Int_t msty = 20;
  Float_t msiz = 0.8;

  while (infile.getline(st,256)) {
    if (st[0] == '#') {
      if (strstr(st,"Name"))
	sscanf(st,"%s %s %s",stmp, stmp, hn);
      if (strstr(st,"Title")) {
	int n = sscanf(st,"%s %s %s",stmp, stmp,stmp);
	if (n>2) {
	  strcpy(titolo,strstr(st,stmp));
	  char *ret = index(titolo,'\n');
	  if (ret) {*ret=0; *(ret+1)=0;}
	  strcat(titolo," ");
	}
      }
      if (strstr(st,"LineColor"))
	  sscanf(st,"%s %s %d", stmp, stmp, &lcol);
      if (strstr(st,"LineWidth"))
	  sscanf(st,"%s %s %d", stmp, stmp, &lwid);
      if (strstr(st,"MarkerStyle"))
	  sscanf(st,"%s %s %d", stmp, stmp, &msty);
      if (strstr(st,"MarkerColor"))
	  sscanf(st,"%s %s %d", stmp, stmp, &mcol);
      if (strstr(st,"MarkerSize"))
	  sscanf(st,"%s %s %f", stmp, stmp, &msiz);
    }
    else {
      exx=0; eyy=0;
      int nn=sscanf(st,"%f %f %f %f",&xx,&yy,&exx,&eyy);
      if (nn>=2) {
	x[n]=xx;
	y[n]=yy;
	ex[n]=exx;
	ey[n]=eyy;
	n++;
      }
    }
  }


  if (gROOT->FindObject(hn)) gROOT->FindObject(hn)->Delete();
  g = new TGraphErrors(n,x,y,ex,ey);

  g->SetName(hn);
  g->SetTitle(titolo);
  g->SetLineColor(lcol);
  g->SetLineWidth(lwid);
  g->SetMarkerStyle(msty);
  g->SetMarkerColor(mcol);
  g->SetMarkerSize(msiz);

  return g;
}

//////////////////////////////////////////////////////////////////
Bool_t InpYN(Char_t *question) 
{
  // stampa la domanda e procede solo se risponde 'y' o 'Y' (1) o 'n' o 'N' (0)
  printf("%s  (y/n)  ",question);
  fflush(stdout);
  Char_t rsp=' ';
  while (rsp!='y' && rsp!='n' && rsp!='Y' && rsp!='N') {
    rsp=getchar();
    if (rsp!='y' && rsp!='n' && rsp!='Y' && rsp!='N') {
      rsp=getchar();
      printf("please answer y or n!\n");
    }
  } 
  if (rsp=='y' || rsp=='Y') return 1;
  return 0;
}
//////////////////////////////////////////////////////////////////
TH1F* LoadHist(Char_t *fname, Char_t *histname="")
{
  // legge un istogramma da file tipo testo nel formato:
  //
  // # commenti
  // i1    i2   ...  iN  
  // iN+1  ...  ...  ...
  // 
  // con N (numero di colonne = 1,...,10 )
  // -se histname="" usa il nome del file

  FILE *pf = NULL;
  pf = fopen(fname,"r");
  if (!pf) {
    printf("file not found!\n");
    return NULL;
  }
  Char_t hn[256];

  if (!strcmp(histname,"")) {
    char *pio = strrchr(fname,'/');
    if (pio) strcpy(hn,++pio);
    else strcpy(hn,fname);
    Int_t pospt = strcspn(hn,".");
    if (pospt>90) pospt=10;
    hn[pospt] = 0;
    hn[pospt+1] = 0;
  }  
  else 
    strcpy(hn,histname);

  PCheckName(hn);

  TH1F *h = NULL;
  Float_t y[16384];
  Float_t tmp[10];
  Char_t st[256];
  Char_t stmp[30];
  Char_t titolo[100];
  strcpy(titolo,"");
  Int_t n=0;

  Int_t lcol = 1;
  Int_t lwid = 1;
  Float_t xmin = 0;
  Float_t xmax = 0;

  while (fgets(st,255,pf)) {
    if (st[0] == '#') {
      if (strstr(st,"Name"))
	sscanf(st,"%s %s %s",stmp, stmp, hn);
      if (strstr(st,"Title")) {
	int n = sscanf(st,"%s %s %s",stmp, stmp,stmp);
	if (n>2) {
	  strcpy(titolo,strstr(st,stmp));
	  char *ret = index(titolo,'\n');
	  if (ret) {*ret=0; *(ret+1)=0;}
	  strcat(titolo," ");
	}
      }
      if (strstr(st,"LineColor"))
	  sscanf(st,"%s %s  %d",stmp, stmp, &lcol);
      if (strstr(st,"LineWidth"))
	  sscanf(st,"%s %s  %d",stmp, stmp, &lwid);
      if (strstr(st,"Xmin"))
	  sscanf(st,"%s %s  %f",stmp, stmp, &xmin);
      if (strstr(st,"Xmax"))
	  sscanf(st,"%s %s  %f",stmp, stmp, &xmax);
    }
    else {
      Int_t ndata=sscanf(st,"%f %f %f %f %f %f %f %f %f %f",
			 &tmp[0],&tmp[1],&tmp[2],&tmp[3],&tmp[4],
			 &tmp[5],&tmp[6],&tmp[7],&tmp[8],&tmp[9]);
   
      for (Int_t j=0; j<ndata; j++) y[n++]=tmp[j];
    }
  }
  fclose(pf);

  if (gROOT->FindObject(hn)) gROOT->FindObject(hn)->Delete();
  h = new TH1F(hn,titolo,n,0,n);
  h->SetBinContent(0,0);
  for (Int_t i=0; i<n; i++) {
    h->SetBinContent(i+1,y[i]);
  }
  h->SetLineColor(lcol);
  h->SetLineWidth(lwid);
  if (h->GetNbinsX() && (xmin<xmax)) 
    h->GetXaxis()->Set(h->GetNbinsX(),xmin,xmax);
  h->SetEntries(h->Integral());

  sprintf(st,"TH1F *%s = gROOT->FindObject(\"%s\");",hn,hn);
  gROOT->ProcessLine(st);
  return h;
}
//////////////////////////////////////////////////////////////////
Int_t SaveGraphErrors(TGraphErrors *g, Char_t *fname, Int_t attr=1)
{
  // Scrive un TGraphErrors in formato ascii nel file fname
  // - se attr=1 prima di scrivere controlla se c'e' gia' un file
  //   con lo stesso nome
  if (!g) return -1;
  
  FILE *pf = NULL;
  if (attr==1) {
    pf = fopen(fname,"r");
    if (pf) {
      if (!InpYN("file exists: overwrite? "))
	return -1;
    }
  }
  pf = fopen(fname,"w");
  if (!pf) {
    printf("cannot open file!\n");
    return -2;
  }
  
  if (attr==1) {
    fprintf(pf,"# graph options:\n");
    fprintf(pf,"#   Name         %s\n",g->GetName());
    fprintf(pf,"#   Title        %s\n",g->GetTitle());
    fprintf(pf,"#   LineColor    %d\n",g->GetLineColor());
    fprintf(pf,"#   LineWidth    %d\n",g->GetLineWidth());
    fprintf(pf,"#   MarkerStyle  %d\n",g->GetMarkerStyle());
    fprintf(pf,"#   MarkerColor  %d\n",g->GetMarkerColor());
    fprintf(pf,"#   MarkerSize   %f\n",g->GetMarkerSize());
    fprintf(pf,"# graph data:\n");
  }

  for (Int_t i=0; i<g->GetN(); i++)
    fprintf(pf,"%f  %f  %f  %f\n", g->GetX()[i], g->GetY()[i],
	    g->GetEX()[i], g->GetEY()[i]);

  fclose(pf);
  return 0;
}
/////////////////////////////////////////////////////////////////////////
Int_t SaveHist(TH1 *h, Char_t *fname, Int_t attr=1)
{
  // Scrive un TH1 in formato ascii nel file fname
  // - se attr=1 prima di scrivere controlla se c'e' gia' un file
  //   con lo stesso nome

  if (!h) return -1;
  
  FILE *pf = NULL;
  if (attr==1) {
    pf = fopen(fname,"r");
    if (pf) {
      if (!InpYN("file exists: overwrite? "))
	return -1;
    }
  }

  pf = fopen(fname,"w");
  if (!pf) {
    printf("cannot open file!\n");
    return -2;
  }
  if (attr==1) {
    fprintf(pf,"# histo options:\n");
    fprintf(pf,"#   Name       %s\n",h->GetName());
    fprintf(pf,"#   Title      %s\n",h->GetTitle());
    fprintf(pf,"#   LineColor  %d\n",h->GetLineColor());
    fprintf(pf,"#   LineWidth  %d\n",h->GetLineWidth());
    fprintf(pf,"#   Xmin       %f\n",h->GetXaxis()->GetXmin());
    fprintf(pf,"#   Xmax       %f\n",h->GetXaxis()->GetXmax());
    fprintf(pf,"# histo data:\n");
  }

  for (Int_t i=0; i<h->GetNbinsX(); i++)
    fprintf(pf,"%f\n", h->GetBinContent(i+1));
  
  fclose(pf);
  return 0;
}
/////////////////////////////////////////////////////////////////////////