// random sample = N "measurements of Poisson variable"
const Int_t N = 20;
const Int_t NPE = 1000;

void ml_bootstrap_gauss()
{
  // true parameters
  double mean  = 4.5;
  double sigma = 1.5;

  // first experiment
  double firstmean, firstsigma;
  doPE( firstmean, firstsigma, mean, sigma );
  cout << "First mean:  " << firstmean << endl;
  cout << "First sigma: " << firstsigma << endl;

  TH1* h = new TH1F( "h", "Parametrisches Bootstrapping", 100, 0, 10 );

  double newmean  = 0;
  double newsigma = 0;
  double PEs[NPE];
  for( Int_t i = 0; i < NPE; ++i ) {
    doPE( newmean, newsigma, firstmean, firstsigma );
    PEs[i] = newmean;
    h->Fill( newmean );
  }
  cout << "Standard deviation of mean: " << TMath::RMS( NPE, PEs )
       << "\t (expected: " << sigma/sqrt(N) << ")" << endl;
  h->Draw();

}

void doPE( Double_t& meanout, Double_t& sigmaout,
	   const Double_t& meanin, const Double_t& sigmain )
{
  Double_t x[N];
  for( Int_t i = 0; i < N; ++i ) {
    x[i] = gRandom->Gaus( meanin, sigmain );
  }
  meanout  = TMath::Mean( N, x );
  sigmaout = TMath::RMS( N, x );
}