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

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

  // gere M amostras de N elementos da mesma gaussiana 
  // faca um histograma da media das M amostras e verifique que a largura
  // e' sigma/sqrt(N)
  Int_t M=100;
  TH1F *hm  = new TH1F("hm"," ", Nbins,-W,W);
  for (Int_t j=1;j<=M; j++){
    TH1F *ht  = new TH1F("ht"," ", Nbins,-W,W);  
    for (Int_t i=1;i<=N;i++){
     Double_t x = r3->Gaus(mu,sigma);
     ht->Fill(x);
    }
    Double_t m = ht->GetMean();
    hm->Fill(m);
    delete ht;
  }
  Double_t hmrms = hm->GetRMS();
  cout << " RMS das medias = " <<  hmrms << endl; 
  hm->Draw("same");
  cv->Update();
  
}