double lnL(double* x, int n, double lambda, double xi)
{
	double sum = 0;
	for (int i = 0; i < n; i++)
		sum += log(pow(lambda,2) + pow(x[i]-xi, 2));
	return n*log(lambda) - n*log(TMath::Pi()) - sum;
}

void SetTextOpts(TH1* h)
{
	h->SetTitleSize(0.05,"xyz");
	h->SetLabelSize(0.05,"xyz");
	h->SetTitleFont(132,"xyz");
	h->SetLabelFont(132,"xyz");
}

void upf_cv9_BW()
{
	double lambda = 5;
	double xi     = 2;
	const int N = 50;
	double x[N];

	int m = 100;
	double lambdaMin = 0.1;
	double lambdaMax = 10;
	double dLambda   = (lambdaMax-lambdaMin) / m;
	double xiMin = -40;
	double xiMax =  40;
	double dXi   = (xiMax-xiMin) / m;

	printf("\nUsed for generation of %i x: lambda = %.2f, estimated xi = %.2f\n",
			N,lambda,xi);

	TH1D* h = new TH1D("h",";#it{x};entries",10,xiMin,xiMax);

	TRandom3* r = new TRandom3();
	for (int i = 0; i < N; i++)
	{
		x[i] = r->BreitWigner(xi,lambda*2);
		h->Fill(x[i]);
	}

	TGraph2D* gr = new TGraph2D(m*m);
	TGraph2D* grMax = new TGraph2D(1);

	double lnLMax, lambdaEst, xiEst;

	for (int j = 0; j < m; j++)
	{
		double lambda_j = lambdaMin + dLambda*j;
		for(int k = 0; k < m; k++)
		{
			double xi_k = xiMin + dXi*k;
			double lnL_jk = lnL(x,N,lambda_j,xi_k);
			gr->SetPoint(k+m*j, lambda_j, xi_k, lnL_jk);
			if (k+m*j == 0 || lnL_jk > lnLMax)
			{
				lnLMax    = lnL_jk;
				lambdaEst = lambda_j;
				xiEst     = xi_k;
			}
		}
	}
	printf("\nEstimated: lambda = %.2f, xi = %.2f\n\n",lambdaEst,xiEst);

	gStyle->SetOptStat(0);
	gStyle->SetPalette(55);
	TCanvas* c = new TCanvas("c","",1400,500);
	c->Divide(3);

	c->cd(1);
	gPad->SetLeftMargin(0.12);
	SetTextOpts(h);
	h->Draw();

	c->cd(2);
	gr->SetTitle("");
	gr->SetMarkerStyle(20);
	gr->SetMarkerSize(0.5);
	gr->Draw("pcol");
	gr->GetHistogram()->GetXaxis()->SetTitle("#lambda");
	gr->GetHistogram()->GetYaxis()->SetTitle("#xi");
	gr->GetHistogram()->GetZaxis()->SetTitle("ln(#it{L})");
	gr->GetHistogram()->GetZaxis()->SetTitleOffset(1.5);
	gr->GetHistogram()->GetZaxis()->SetRangeUser(gr->GetMinimum(),lnLMax+1);
	SetTextOpts(gr->GetHistogram());

	c->cd(3);
	TGraph2D* grZoom = (TGraph2D*)gr->Clone();
	grZoom->Draw("pcol");
	grZoom->GetHistogram()->GetXaxis()->SetTitle("#lambda");
	grZoom->GetHistogram()->GetYaxis()->SetTitle("#xi");
	grZoom->GetHistogram()->GetZaxis()->SetTitle("ln(#it{L})");
	grZoom->GetHistogram()->GetZaxis()->SetTitleOffset(1.5);
	SetTextOpts(grZoom->GetHistogram());
	grZoom->GetHistogram()->GetXaxis()->SetRangeUser(3.0,7.0);
	grZoom->GetHistogram()->GetYaxis()->SetRangeUser(-4.0,5.0);
	grZoom->GetHistogram()->GetZaxis()->SetRangeUser(gr->GetMinimum(), lnLMax+1);
}