#include "rootlogon.h"

#include <TRandom.h>
#include <TMath.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TLine.h>

TF1* moyal = new TF1("moyal", "1/sqrt(2*pi)*exp(-1.0/2*((1.27+(x-7.5))+exp(-(1.27+(x-7.5)))))", 0, 20);

void once(double* mean, double* meanmedian)
{
  double vals[20];
  for(unsigned int i = 0; i < 20; ++i)
  {
    vals[i] = moyal->GetRandom();
  }

  std::sort(vals, vals+20);

  *mean = TMath::Mean(vals, vals+20);
  *meanmedian = TMath::Mean(vals, vals+11);
}

int main(int argc, char* argv[])
{
  rootlogon();

  TH1F* mean = new TH1F("mean", "", 125, 0, 20);
  TH1F* meanmedian = new TH1F("median", "", 125, 0, 20);

  for(unsigned int i = 0; i < 100000; ++i)
  {
    double meanval, meanmedianval;
    once(&meanval, &meanmedianval);

    mean->Fill(meanval);
    meanmedian->Fill(meanmedianval);
  }

  mean->Scale(1/mean->Integral());
  meanmedian->Scale(1/meanmedian->Integral());

  TCanvas* cnv = new TCanvas;

  moyal->SetTitle("");
  moyal->Draw();
  TLine line(7.5, 0.0, 7.5, 0.253);
  line.SetLineColor(kRed);
  line.Draw("same");
  cnv->SaveAs("out/moyal.eps");

  mean->SetLineColor(kRed);
  meanmedian->SetLineColor(kBlue);
  mean->GetXaxis()->SetRangeUser(3,10);
  meanmedian->GetXaxis()->SetRangeUser(3,10);
  meanmedian->Draw("");
  mean->Draw("same");
  cnv->SetMargin(0.15, 0.9, 0.10, 0.9);
  cnv->SaveAs("out/bias.eps");
}