void Fill1(){
  // Gaussian model
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);
  gStyle->SetLabelSize(0.05,"XY");
  gStyle->SetTitleSize(0.05,"XYZ");
  gStyle->SetTitleOffset(1,"XYZ");

  double mux = 0;
  double sigmax = 1;
  TF1 *fg = new TF1("fg", "[0]*TMath::Gaus(x,[1],[2])", -5, 5);
  //fg->SetParameters(1. / (sqrt(2 * TMath::Pi())*sigmax), mux, sigmax);
  fg->SetParameters(100. / (sqrt(2 * TMath::Pi())*sigmax), mux, sigmax);
  fg->GetXaxis()->SetTitle("x"); fg->GetYaxis()->SetTitle("f_{x}(x)");
  fg->SetTitle("");
  fg->SetLineColor(2);
  //fg->SetFillStyle(1001);
  //fg->SetFillColor(4);
  fg->SetNpx(1000);

  TH1F *h = new TH1F("h","",100,-5,5);
  h->FillRandom("gaus",1000);
  h->SetMarkerStyle(20);
  h->GetXaxis()->SetTitle("x"); h->GetYaxis()->SetTitle("Events");
  h->Draw("e");
  fg->Draw("same");

  TLegend *l = new TLegend(0.7, 0.6, 0.9, 0.8); l->SetBorderSize(0);
  l->AddEntry(h, "Data", "p");
  l->AddEntry(fg, "Model", "l");
  l->Draw();
}

void Fill2(){
  // Exp. model

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);
  gStyle->SetLabelSize(0.05,"XY");
  gStyle->SetTitleSize(0.05,"XYZ");
  gStyle->SetTitleOffset(1,"XYZ");

  // Hist
  int nbin = 100;
  double xmin = 0;
  double xmax = 18;
  int ntot = 500; // npexp
  TH1F *h = new TH1F("h","",nbin,xmin,xmax);
  h->SetMarkerStyle(20);
  h->GetXaxis()->SetTitle("x"); h->GetYaxis()->SetTitle("Events");

  // Model
  TF1 *fg = new TF1("fg", "[0]*TMath::Exp(-[1]*x)",xmin,xmax);
  double lambda=0.40;
  fg->SetParameter(0,lambda*ntot*h->GetBinWidth(1));
  fg->SetParameter(1,lambda);
  fg->SetNpx(1000);
  h->FillRandom("fg",ntot);

  // Draw
  h->Draw("e");
  lambda = 0.35; // different model
  fg->SetParameter(0,lambda*ntot*h->GetBinWidth(1));
  fg->Draw("same");
  TLegend *l = new TLegend(0.7, 0.6, 0.9, 0.8); l->SetBorderSize(0);
  l->AddEntry(h, "Data", "p");
  l->AddEntry(fg, "Model", "l");
  l->Draw();
}

void Fill3(){
	// Calc. variance of exp. model with pseudo-exp

	gStyle->SetOptStat(0);
	gStyle->SetOptFit(1111);
	gStyle->SetLabelSize(0.05, "XY");
	gStyle->SetTitleSize(0.05, "XYZ");
	gStyle->SetTitleOffset(1, "XYZ");

	int Npexp = 500;
	int Nsample = 200;

	// Model
	double fmin = 0;
	double fmax = 20;
	TF1 *fg = new TF1("fg", "([1]/[0])*TMath::Exp(-x/[0])", fmin, fmax);
	double tau = 2;
	fg->SetParameter(0, tau);
	fg->SetParameter(1, 1);
	fg->SetNpx(1000);
	TH1F *hg = new TH1F("hg", "", 100, fmin, fmax);
	hg->SetMarkerStyle(20); hg->GetXaxis()->SetTitle("x"); hg->GetYaxis()->SetTitle("Events");

	// Hist
	int nbin = 100;
	double xmin = 1;
	double xmax = 3;
	TH1F *h = new TH1F("h", "", nbin, xmin, xmax);
	h->SetMarkerStyle(20);
	h->GetXaxis()->SetTitle("#hat{#tau}"); h->GetYaxis()->SetTitle("N experiment");

	// calc. likelihood estimator on tau (=mean)
	double tau_mean = 0;
	double tau_var = 0;
	for (int i = 0; i < Npexp; i++){
		double tau_hat = 0;
		for (int j = 0; j < Nsample; j++){
			double t = fg->GetRandom();
			tau_hat += t / Nsample;
			if ( i==0 ) hg->Fill(t);
		}
		h->Fill(tau_hat);
		tau_mean += tau_hat / Npexp;
		tau_var += pow(tau_hat, 2) / Npexp;
	}
	tau_var = tau_var - pow(tau_mean, 2);

	cout << "Mean = " << tau_mean << ", variance = " << tau_var << ", std = " << sqrt(tau_var) << endl;
	cout << "Expected std = " << tau / sqrt(Nsample);

	new TCanvas(); fg->Draw();
	new TCanvas(); hg->Draw("e1"); fg->SetParameter(1, Nsample*hg->GetBinWidth(1)); fg->Draw("same");
	new TCanvas(); h->SetLineColor(50);  h->SetFillColor(5); h->Draw();

}