#!/bin/bash
G4BL_DIR=`g4bl-config --dir`
if ! test -r "$G4BL_DIR/.config"
then
	# some necessary code is not present
	echo test55 -- Compiling BLMinimize omitted
	exit 0
fi

source $G4BL_DIR/.config

if test $SYSTEM = WIN32 -o \
	! -r $G4BL_DIR/include/gsl/gsl_multimin.h -o \
	! -r $G4BL_DIR/lib/libgsl.a -o \
        ! -r $SOURCE/source/include/BLMinimize.hh || \
	! type g++ >/dev/null 2>&1
then
	# some necessary code is not present
	echo test55 -- Compiling BLMinimize omitted
	exit 0
fi
echo "test55: compile and run BLMinimize (1 sec)"

cat >test55.cc <<\!
#include <assert.h>
#include <math.h>
#include <vector>
#include "BLMinimize.hh"
double gaussran(double sigma=1.0)
{
	static double A[35] = { // Gaussian, mean=0 sigma=1
		-0.3836,1.3049,0.1900,-0.2952,0.3521,-1.5032,-0.8877,
		0.0026,0.5015,-0.0375,1.2435,-0.9606,-0.7007,0.3411,
		-0.9010,1.5105,1.3040,0.0778,0.8959,-2.0137,0.2735,
		0.6965,-1.8167,0.4734,-0.9461,0.0591,0.7711,2.5389,
		1.0801,-0.2345,-0.3373,-0.1851,-0.2302,-0.2688,-0.0623
	};
	static int i=0;
	return A[++i%35]*sigma;
}
class Track : public BLMinimizeFunction {
public:
	int nHit;
	double x[20];
	double y[20];
	double z[20];
	virtual double operator()(const std::vector<double> &arg) {
		double x0 = arg[0];
		double y0 = arg[1];
		double unused = arg[2];
		assert(fabs(unused-9.876)<0.01);
		double dxdz = arg[3];
		double dydz = arg[4];
		double chisq = 0.0;
		for(int i=0; i<nHit; ++i) {
			double dx = x0 + dxdz*z[i] - x[i];
			double dy = y0 + dydz*z[i] - y[i];
			chisq += dx*dx/4.0 + dy*dy/4.0;
		}
		return chisq/(2*nHit-4);
	}
	void setup(double x0, double y0, double dxdz, double dydz)
	{
		nHit = 7;
		for(int i=0; i<nHit; ++i) {
			z[i] = i*1000;
			x[i] = x0 + dxdz*z[i] + gaussran(2.0);
			y[i] = y0 + dydz*z[i] + gaussran(2.0);
		}
		printf("%d\t%.4f\t%.4f\t%d\t%.2f\t%.2f\t%.4f\t%.4f\n",
			0,0.0,0.0,0,x0,y0,dxdz,dydz);
	}
};
void G4Exception(char const* origin, char const *code, 
	G4ExceptionSeverity severity, char const *desc)
{
	printf("G4EXCEPTION: %s %s %s\n",origin, code, desc);
	exit(99);
}
int main(int argc, char *argv[])
{
	std::vector<double> x;
	x.push_back(1.0);
	x.push_back(1.0);
	x.push_back(0.0); // unused - fixed in fit
	x.push_back(0.002);
	x.push_back(0.002);

	Track t;
	BLMinimize m(0);  // use x[] for scales
	m.setFunc(t,x);

	x[2] = 9.876; // fixed value for unused entry

	int status=0;
	printf("#ev\tChisqDf\tSize\tIter\tx0\ty0\tdxdz\tdydz\n");

#define GO(N,X0,Y0,DXDZ,DYDZ) \
	x[0]=X0; x[1]=Y0; x[3]=DXDZ; x[4]=DYDZ; \
	t.setup(X0,Y0,DXDZ,DYDZ); \
	assert(m.minimize(x,0.001,200)==0);  \
	printf("%d\t%.4f\t%.4f\t%d\t%.2f\t%.2f\t%.4f\t%.4f\n", \
		N,m.value(),m.getSize(),m.getIterations(), \
		x[0],x[1],x[3],x[4]);

	GO(1,0.0,0.0,0.10,0.10);
	GO(2,10.0,-10.0,0.10,-0.10);
	GO(3,300.0,500.0,0.010,0.050);
	GO(4,-70.0,440.0,0.010,-0.080);
	GO(5,19.0,-135.0,0.066,-0.144);
	GO(6,312.0,58.0,0.010,0.050);

	return status;
}
!

cat >test55.good <<-\!
#ev	ChisqDf	Size	Iter	x0	y0	dxdz	dydz
0	0.0000	0.0000	0	0.00	0.00	0.1000	0.1000
1	0.8297	0.0009	65	0.27	0.64	0.0998	0.0998
0	0.0000	0.0000	0	10.00	-10.00	0.1000	-0.1000
2	1.8834	0.0009	64	9.77	-9.10	0.1004	-0.1002
0	0.0000	0.0000	0	300.00	500.00	0.0100	0.0500
3	0.4893	0.0010	94	299.69	500.06	0.0100	0.0499
0	0.0000	0.0000	0	-70.00	440.00	0.0100	-0.0800
4	1.2532	0.0009	69	-69.39	439.93	0.0100	-0.0800
0	0.0000	0.0000	0	19.00	-135.00	0.0660	-0.1440
5	0.9820	0.0009	77	17.36	-133.15	0.0664	-0.1444
0	0.0000	0.0000	0	312.00	58.00	0.0100	0.0500
6	0.8297	0.0009	65	312.27	58.64	0.0098	0.0498
!

if g++ -I$G4BL_DIR/include -I$G4BL_DIR/include/Geant4 -I$SOURCE/source/include \
   -o test55.exe test55.cc $SOURCE/source/src/BLMinimize.cc -L$G4BL_DIR/lib\
   -lgsl -lgslcblas -lm >test55.out1 2>&1
then
	./test55.exe >test55.out
else
	echo "****** test55: compilation FAILED -- test55.* files left intact"
	exit 1
fi

if ! ./comparecols test55.out test55.good 1 0.001 0.001 10 0.1 0.1 0.001 0.001
then
	echo "****** test55: Output Files DIFFER -- test55.* files left intact"
	exit 1
fi

# clean up and exit
rm test55.*
exit 0