// rms de uma distribuicao e erro na media (Cauchy)
//  ref Cowan, cap 5
//  jtmn, Set 2003

  Double_t cauchy(Double_t r, Double_t a, Double_t gam){
  // r random number 0-1
  Double_t pi = TMath::Pi();
  Double_t y = a + gam*TMath::Tan(pi*(r-0.5))/2.0;
  return y;
}


void err_cauchy()
{
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  TRandom *r3 = new TRandom3();
  r3->SetSeed(0);
  
 // simule uma amostra de N elementos de uma lorentziana
  TCanvas *cv = new TCanvas("cv"," Hist  ", 600,600);
  Int_t Nbins = 100;
  Int_t N  =1000;
  Double_t mu = 0.0;
  Double_t gam = 20.0;
  Double_t W = 1000;
  TH1F *h  = new TH1F("hexp"," ", Nbins,-W,+W);
  for (Int_t i=1;i<=N;i++){
    Double_t r = r3->Rndm();
    Double_t c = cauchy(r,mu,gam); 
   h->Fill(c);
  }
  // mostre o histograma
  h->Draw();
  Double_t hrms = h->GetRMS();
  cout << " RMS da distribuicao = " << hrms << endl;
  cv->Update();

  // gere M amostras de N elementos da mesma lorentziana 
  // faca um histograma da media das M amostras e verifique que a largura
  // e' igual a largura da distribuicao;
  Int_t M=1000;
  Double_t media;
  TH1F *hm  = new TH1F("hm"," ", Nbins,-W,W);
  for (Int_t j=1;j<=M; j++){
    media = 0.0;  
    for (Int_t i=1;i<=N;i++){
     Double_t r = r3->Rndm();
     Double_t x = cauchy(r,mu,gam);
     media = media + x;
    }
    // nao use histograma e getMean por causa dos over e underflows!!!
    media = media/Double_t(N);
    hm->Fill(media);
  }
  Double_t hmrms = hm->GetRMS();
  cout << " RMS das medias = " <<  hmrms << endl; 
  hm->Draw("same");
  cv->Update();
  
}