void poisson()
{
  //====================================================================
  //   Laboratorio de Fisica Nuclear y de Particulas
  //   
  //   Macro para construir un histograma de frecuencias y ajustar una función de Poisson
  //====================================================================
  
  // Algunas opciones de ROOT
  gROOT->Reset();
  //gROOT->SetStyle("Plain");
  gStyle->SetOptFit(111);
  gStyle->SetOptStat(0);
  
  // Definir entrada  
  // ...Fichero de datos
  string filename = "fondo.dat";
  // ...Contar número de líneas de la cabecera del archivo
  int nlines = 1;
  // ...Definir histogramas
  // ...de anchura de bin constante
  int hist_ini = 0;         // lower edge will be in fact hist_ini-hist_binwidth/2
  int hist_fin = 10;        // upper edge will be in fact hist_ini+hist_binwidth/2
  float hist_binwidth = 1.; // constant bin width
  // ...de anchura de bin user defined (variable)
  const int hist_nbins = 6;
  Double_t edges[hist_nbins+1] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 10.5};
  // Bin 1 corresponds to range [-0.5, 0.5], thus 0 is in the middle
  // Bin 2 corresponds to range [0.5, 1.5], thus 1 is in the middle, etc...
  
  // Crea histograma
  TH1F* hist = new TH1F("hist","Histograma de anchura constante",hist_fin-hist_ini+1,hist_ini-0.5*hist_binwidth,hist_fin+0.5*hist_binwidth) ;  // constant bin width

  // Declarar variables
  char line[200], name[80];
  float x=0,y=0,z=0,u=0;
  int row=0,nmeasu=0;
  float mean=0;
  float zarray[1000]; 
  
  // Abrir fichero ascii y extraer valores
  FILE *f= fopen(filename.c_str(),"r");
  while (fgets(line,200,f)) 
    {
      //printf("%d: %s",row,line);
      if(row>=nlines)
	{
	  // Saltar cabecera
	  sscanf(&line[0],"%f %f %f %f ",&x,&y,&z,&u);
	  // Llenar histograma
	  hist->Fill(z);
	  zarray[nmeasu] = z;
	  mean += z;
	  nmeasu++;
	}
      row++;
    }
  fclose(f);  // Cerrar fichero
  
  // media y su incertidumbre
  mean = mean/nmeasu;
  float var=0;
  for (int i=0;i<nmeasu;i++)
    var += pow(zarray[i]-mean,2);
  var = var/(nmeasu-1)/nmeasu;

  // crea histograma de anchura de bin user defined (variable)
  TH1F* hist_wvar = dynamic_cast<TH1F*>(hist->Rebin(hist_nbins,"hist_wvar",edges));   // variable bin width
  
  // construye correspondientes histogramas esperados
  TH1F* hist_exp = dynamic_cast<TH1F*>(hist->Clone());
  for (int i=1;i<=hist_exp->GetNbinsX();i++) {
    float binCenter = hist_exp->GetBinCenter(i);
    hist_exp->SetBinContent(i,nmeasu*TMath::Poisson(binCenter,mean));
    hist_exp->SetBinError(i,sqrt(nmeasu*TMath::Poisson(binCenter,mean)));
  }
  TH1F* hist_wvar_exp = dynamic_cast<TH1F*>(hist_exp->Rebin(hist_nbins,"hist_wvar_exp",edges));   // variable bin width
  
  // Definir canvas
  TCanvas* c1 = new TCanvas ("c1","c1",640,480);  
  TCanvas* c2 = new TCanvas ("c2","c2",640,480);  

  // Declarar funcion de Poisson
  TF1 *poisson = new TF1("poisson","[0]*TMath::Poisson(x,[1])", hist_ini, hist_fin);
  poisson->SetParameters(1, 1, 1); // MUST set non-zero initial values for parameters
  poisson->FixParameter(0,nmeasu);
  
  // Opciones finales, dibujar histograma y Poisson
  c1->cd();
  c1->SetTicks();
  hist->GetXaxis()->SetTitle("# cuentas");
  hist->GetXaxis()->SetTitleOffset(0.8);  
  hist->GetXaxis()->SetTitleSize(0.05);
  hist->GetYaxis()->SetTitle("frecuencia");
  hist->GetYaxis()->SetTitleOffset(0.8);  
  hist->GetYaxis()->SetTitleSize(0.05);
  hist->SetMarkerStyle(20);
  hist->SetMarkerSize(0.8);
  hist->SetMarkerColor(4);
  hist->SetTitle("");
  hist->Fit("poisson","R"); // "R" = fit between "xmin" and "xmax" of the "f1"
  hist->Draw("E SAME");
  //
  c2->cd();
  c2->SetTicks();
  hist_wvar->GetXaxis()->SetTitle("# cuentas");
  hist_wvar->GetXaxis()->SetTitleOffset(0.8);  
  hist_wvar->GetXaxis()->SetTitleSize(0.05);
  hist_wvar->GetYaxis()->SetTitle("frecuencia");
  hist_wvar->GetYaxis()->SetTitleOffset(0.8);  
  hist_wvar->GetYaxis()->SetTitleSize(0.05);
  hist_wvar->SetMarkerStyle(20);
  hist_wvar->SetMarkerSize(0.8);
  hist_wvar->SetMarkerColor(4);
  hist_wvar->SetTitle("");
  hist_wvar->Draw("*H");
  hist_wvar_exp->SetLineColor(kRed);
  hist_wvar_exp->SetLineWidth(2.);
  hist_wvar_exp->Draw("SAME E");
  TLegend* legend = new TLegend(0.65,0.68,0.90,0.90);  // x1,y1,x2,y2
  legend->AddEntry(hist_wvar,"Observado","l");
  legend->AddEntry(hist_wvar_exp,"Esperado","l");
  legend->Draw();
        
  // Tabla con resultado del ajuste:
  printf("\nMedia de las %d medidas: %4.7f +/- %4.7f\n",nmeasu,mean,sqrt(var));
  printf("Ajuste histograma de anchura de bin constante:\n");
  printf("  Media: %4.7f +/- %4.7f\n",poisson->GetParameter(1),poisson->GetParError(1));
  printf("  Chi2 / ndf: %f / %d\n",hist->GetFunction("poisson")->GetChisquare(),hist->GetFunction("poisson")->GetNDF());
  printf("  CL: %f\n",TMath::Prob(hist->GetFunction("poisson")->GetChisquare(),hist->GetFunction("poisson")->GetNDF()));
  printf("Comparacion histogramas (observado-esperado) de anchura de bin user defined (variable):\n");
  float chi2=0;
  for (int i=1;i<=hist_wvar_exp->GetNbinsX();i++) {
    chi2 += pow((hist_wvar->GetBinContent(i)-hist_wvar_exp->GetBinContent(i))/hist_wvar_exp->GetBinError(i),2);
  }
  int ndof=hist_wvar_exp->GetNbinsX()-1;
  printf("  Chi2 / ndf: %f / %d\n",chi2,ndof);
  printf("  CL: %f\n\n",TMath::Prob(chi2,ndof));
  
  // dump canvas to pdf
  c1->Print("poisson_fit.pdf");  
  c1->Print("poisson_fit.png");
  c2->Print("poisson.pdf");  
  c2->Print("poisson.png");
}