// Code for 1-cut
#include <oxstd.h>
#include <oxprob.h>
#include <oxdraw.h>

#import "../code/sim_1cut"

PlotBiasCorrection(const cM, const cT)
{
	decl mbeta_sel, mbeta_bc2step_sel, i, psi = <4;2;1;0>, cbars = 100;
	
	[mbeta_sel, mbeta_bc2step_sel] = Run_1cut(cM, <4;2;1> / sqrt(cT), 4, cT, 0.05);
	
	SetDrawWindow("1CutSmallUncond");
	for (i = 0; i < 4; ++i)
	{
		DrawDensity(i, mbeta_sel[i][], "", TRUE, TRUE, FALSE, FALSE, FALSE, cbars, 2);
		DrawTitle(i, sprint("$\psi = ", psi[i], "$"));
	}
	ShowDrawWindow();

	SetDrawWindow("1CutSmallCond");
	for (i = 0; i < 4; ++i)
	{
		DrawDensity(i, mbeta_sel[i][vecindex(mbeta_sel[i][] .!= 0)], "", TRUE, TRUE, FALSE, FALSE, FALSE, cbars, 2);
		DrawDensity(i, mbeta_bc2step_sel[i][vecindex(mbeta_bc2step_sel[i][] .!= 0)], "", TRUE, TRUE, FALSE, FALSE, FALSE, cbars, 3);
		DrawTitle(i, sprint("$\psi = ", psi[i], "$"));
	}
	ShowDrawWindow();
}
main()
{
	// impact of bias correction on small 1-cut experiment
	PlotBiasCorrection(1000000, 100);
}