void landau_gaus_test(void)
{
	gROOT->Reset();

	TH1D *sum_gaus = new TH1D("sum_gaus","", 1000, 0.0, 100.0);
	TH1D *cl_gaus = new TH1D("cl_gaus","",100, 0.0, 1.0);
	
	TH1D* sum_landau = (TH1D*)sum_gaus->Clone("sum_landau");
	TH1D* cl_landau = (TH1D*)cl_gaus->Clone("cl_landau");
	
	TRandom *rand = new TRandom;
	
	Int_t Nmax = 100000;
	for(Int_t i=0; i<Nmax; i++){
	
		// dE/dx gives Landau
		double mean1 = 3.0;
		double sigma1 = 2.0;
		double dE_landau = (rand->Landau(mean1, sigma1) - mean1)/sigma1;
		double dE_gaus = (rand->Gaus(mean1, sigma1) - mean1)/sigma1;
		
		// Tracking chisq is sum of ~ 24 independent gaussian
		int Nhits = 24;
		int Nparam = 5;
		int Ndof = Nhits - Nparam;
		double trk_chisq = 0.0;
		for(int j=0; j<Ndof; j++){
			trk_chisq += pow(rand->Gaus(), 2.0);
		}
		
		// TOF is simple gaussian(?)
		double tof = rand->Gaus();
		
		double fom_gaus = dE_gaus*dE_gaus + trk_chisq + tof*tof;
		double fom_landau = dE_landau*dE_landau + trk_chisq + tof*tof;
		sum_gaus->Fill(fom_gaus);
		sum_landau->Fill(fom_landau);

		// Use TMath::Prob to get probability for this
		cl_gaus->Fill(TMath::Prob(fom_gaus, Ndof+2));
		cl_landau->Fill(TMath::Prob(fom_landau, Ndof+2));
	}

	// Integrate FOM distribution
	TH1D *fom_gaus_int = (TH1D*)sum_gaus->Clone("fom_gaus_int");
	TH1D *fom_landau_int = (TH1D*)sum_landau->Clone("fom_landau_int");
	for(int bin=2; bin<=fom_gaus_int->GetNbinsX(); bin++){
		fom_gaus_int->SetBinContent(bin, fom_gaus_int->GetBinContent(bin-1)+fom_gaus_int->GetBinContent(bin));
		fom_landau_int->SetBinContent(bin, fom_landau_int->GetBinContent(bin-1)+fom_landau_int->GetBinContent(bin));
	}
	fom_gaus_int->Scale(1.0/fom_gaus_int->GetBinContent(fom_gaus_int->GetNbinsX()));
	fom_landau_int->Scale(1.0/fom_landau_int->GetBinContent(fom_landau_int->GetNbinsX()));
	
	TCanvas *c1 = new TCanvas("c1","", 600, 1000);
	c1->Divide(1,3);
	
	TLegend *leg = new TLegend(0.6, 0.55, 0.88, 0.75);
	leg->SetFillColor(kWhite);
	
	c1->cd(1);
	sum_gaus->SetLineColor(kRed);
	sum_landau->SetLineColor(kBlue);
	leg->AddEntry(sum_gaus, "21 Gauss 0 Landau");
	leg->AddEntry(sum_landau, "20 Gauss 1 Landau");
	sum_gaus->SetTitle("FOM (figure of merit)");
	sum_gaus->SetStats(0);
	sum_gaus->Draw();
	sum_landau->Draw("same");
	leg->Draw();
		
	c1->cd(2);
	fom_gaus_int->SetLineColor(kRed);
	fom_landau_int->SetLineColor(kBlue);
	fom_gaus_int->SetTitle("FOM Integral Fraction");
	fom_gaus_int->SetStats(0);
	fom_gaus_int->Draw();
	fom_landau_int->Draw("same");
	leg->Draw();

	c1->cd(3);
	cl_gaus->SetLineColor(kRed);
	cl_landau->SetLineColor(kBlue);
	cl_gaus->GetYaxis()->SetRangeUser(0.0, cl_gaus->GetMaximum()*2.0);
	cl_gaus->SetTitle("Confidence Level (assuming normal errors)");
	cl_gaus->SetStats(0);
	cl_gaus->Draw();
	cl_landau->Draw("same");
	leg->Draw();

	c1->SaveAs("landau_gaus_test.pdf");
	c1->SaveAs("landau_gaus_test.png");
}