/*
 * $Id: usbmax.c,v 1.7 2013/08/13 12:16:02 bnv Exp bnv $
 * $Log: usbmax.c,v $
 * Revision 1.7  2013/08/13 12:16:02  bnv
 * Correcting the 0-bin
 *
 * Revision 1.6  2013/08/08 15:09:08  bnv
 * Correct upper limit
 *
 * Revision 1.5  2012/12/10 13:35:19  bnv
 * Limit printing of z-slices
 *
 * Revision 1.4  2009/10/08 15:51:37  bnv
 * Probably working version with axis selection
 *
 * Revision 1.3  2009/06/19 08:12:31  bnv
 * Added scanning in other axes
 *
 *
 * Revision 1.2  2009/01/14 10:15:35  bnv
 * Working version
 *
 * Revision 1.1  2003/02/18 16:11:39  bnv
 * Initial revision
 *
 */

#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <limits.h>
#include <stdarg.h>

#define X	0
#define Y	1
#define Z	2
#define XYZ	"XYZ"

#define	MAXDET	100
#define RUNTIT	80
#define	RUNTIM	32
#define TITUSB	10
#define	LOPOS(I,A)	(((float)I)*dd[A]+low[A])
#define	HIPOS(I,A)	(((float)I+1.0)*dd[A]+low[A])
#define	MPOS(I,A)	(((float)I+0.5)*dd[A]+low[A])

/* options */
char	*prgname;
float	scanMin[3], scanMax[3];
int	scMin[3], scMax[3];

/* variables */
int	axis=2, axis1=0, axis2=1;// axis direction (default Z)
int	nb;		// detector number;
int	itusbn;		// binning type
int	idusbn;		// scoring distribution

int	nbin[3], rebin[3];	// binning size
float	low[3], high[3], dd[3];	// limits
typedef	char	Str10[TITUSB+1];
Str10	titusb;

int	lntzer;
float	bkusbn;
float	b2usbn;
float	tcusbn;

float	*data  = NULL;
float	*plane = NULL;
float	wctot;
int	nctot;

typedef struct {
	int	pos[3];
	int	Hmin, Hmax;
	int	Vmin, Vmax;
	float	max;
	float	err;
} SliceData;

SliceData	*slice=NULL;

#define FTNGET(V,T,B)	{V=*(T*)B; B+=sizeof(T);}
#define FTNSGET(S,L,B)	{memcpy(S,B,L); S[L]=0; B += L; strip(S);}

/* read a fortran block, return length */
int fortranRead(FILE *f, void *buffer, int maxSize)
{
	int length, length2;
	if (fread(&length,sizeof(length),1,f)<=0)
		return -1;
	if (length>maxSize) {
		fseek(f,sizeof(length),SEEK_CUR);
		return 0;
	}
	if (fread(buffer,length,1,f)<=0)
		return -1;
	if (fread(&length2,sizeof(length2),1,f)<=0)
		return -1;
	if (length!=length2)
		return -1;
	return length;
} /* fortranRead */

/* skip a fortran block, return length */
int fortranSkip(FILE *f)
{
	int length, length2;
	if (fread(&length,sizeof(length),1,f)<=0)
		return -1;
	fseek(f, length, SEEK_CUR);
	if (fread(&length2,sizeof(length2),1,f)<=0)
		return -1;
	if (length!=length2)
		return -1;
	return length;
} /* fortranSkip */

/* read a fortran block an parse it to variables *
int fortranParse(int bufSize, char *buffer, char *pattern, ...)
{
	int length;
	va_list ap;
	char *ptr;

	va_start(ap,pattern);
	for (; *pattern; pattern++) {
		switch (*pattern) {
			case 's':	// string
				break;
			case 'i':	// integer*4
				break;
			case 'f':	// float*4
				break;
			case 'd':	// double*8
				break;
		}
	}
	FTNSGET(runtitle,RUNTIT,pbuffer);
	FTNSGET(runtime,RUNTIM,pbuffer);
	FTNGET(wctot,float,pbuffer);
	FTNGET(nctot,int,pbuffer);
} * fortranParse */

/* strip */
void strip(char *line)
{
	int length = strlen(line)-1;
	while (length>=0 && line[length]==' ')
		length--;
	line[length+1] = 0;
} /* strip */

#if __DEBUG
void hexdump(char *line, int length)
{
	int i,j;
	for (i=j=0; i<length; i++) {
		printf("%02X ",(unsigned char)(line[i]));
		if (++j==16) {
			printf("\n");
			j = 0;
		} else
		if (j==8)
			printf(" ");
	}
	printf("\n");
}
#endif

/* --- scanPlane --- */
void scanPlane(int k)
{
	int	i,j,l;
	int	im=-1, jm=-1;
	float	max = -FLT_MAX;

	// Find maximum and position
	for (j=l=0; j<nbin[axis2]; j++)
		for (i=0; i<nbin[axis1]; i++, l++)
			if (i>=scMin[axis1] && i<scMax[axis1] &&
			    j>=scMin[axis2] && j<scMax[axis2] &&
			    plane[l]>max) {
				max = plane[l];
				im = i;
				jm = j;
			}

	// Find the FHWM/x,y
	float	hmax = max/2.0;
	int	xmin=-1,xmax=-1;
	int	ymin=-1,ymax=-1;
	// Scan on X from 0 to max to find the value where
	//    val[x,y] <= hmax <= val[x+1,y]
	for (i=0; i<nbin[axis1]-1; i++)
		for (j=0; j<nbin[axis2]; j++) {
			if (i<scMin[axis1] || i>scMax[axis1] || j<scMin[axis2] || j>scMax[axis2])
				continue;
			if (plane[j*nbin[axis1]+i]<=hmax &&
				hmax <= plane[j*nbin[axis1]+i+1]) {
					xmin = i;
					goto SEARCH10;
			} else
			if (plane[j*nbin[axis1]+i]>hmax) {
					xmin=i-1;
					goto SEARCH10;
			}
		}
SEARCH10:
	for (i=nbin[axis1]-1; i>=1; i--)
		for (j=0; j<nbin[axis2]; j++) {
			if (i<scMin[axis1] || i>scMax[axis1] || j<scMin[axis2] || j>scMax[axis2])
				continue;
			if (plane[j*nbin[axis1]+i]<=hmax &&
				hmax <= plane[j*nbin[axis1]+i-1]) {
					xmax = i;
					goto SEARCH20;
			} else
			if (plane[j*nbin[axis1]+i]>hmax) {
					xmax = i+1;
					goto SEARCH20;
			}
		}
SEARCH20:
	// Scan on Y
	for (j=0; j<nbin[axis2]-1; j++)
		for (i=0; i<nbin[axis1]-1; i++) {
			if (i<scMin[axis1] || i>scMax[axis1] || j<scMin[axis2] || j>scMax[axis2])
				continue;
			if (plane[j*nbin[axis1]+i]<=hmax &&
				hmax <= plane[j*nbin[axis1]+nbin[axis1]+i]) {
					ymin = j;
					goto SEARCH30;
			} else
			if (plane[j*nbin[axis1]+i]>hmax) {
					ymin = j-1;
					goto SEARCH30;
			}
		}
SEARCH30:
	for (j=nbin[axis2]-1; j>=1; j--)
		for (i=0; i<nbin[axis1]; i++) {
			if (i<scMin[axis1] || i>scMax[axis1] || j<scMin[axis2] || j>scMax[axis2])
				continue;
			if (plane[j*nbin[axis1]+i]<=hmax &&
				hmax <= plane[j*nbin[axis1]-nbin[axis1]+i]) {
					ymax = j;
					goto SEARCH40;
			} else
			if (plane[j*nbin[axis1]+i]>hmax) {
					ymax = j+1;
					goto SEARCH40;
			}
		}
SEARCH40:
	slice[k].pos[axis1] = im;
	slice[k].pos[axis2] = jm;
	slice[k].pos[axis]  = k;
	slice[k].Hmin = xmin;
	slice[k].Hmax = xmax;
	slice[k].Vmin = ymin;
	slice[k].Vmax = ymax;
	slice[k].max = max;
	slice[k].err = 0.0;
} /* scanPlane */

/* --- usage --- */
void usage()
{
	fprintf(stderr,"syntax: %s [options] <infile>\n",prgname);
	fprintf(stderr,"options:\n");
	fprintf(stderr,"\t-h?, --help\tThis help page\n");
	fprintf(stderr,"\t-d, --detector #\tDetector index number starting from 1\n");
	fprintf(stderr,"\t-a, --axis [XYZRP]\tAxis to scan (Default Z)\n");
	fprintf(stderr,"\t-x, --xmin #\tMinimum limit of region in x[cm]\n");
	fprintf(stderr,"\t-X, --xmax #\tMaximum limit of region in x[cm]\n");
	fprintf(stderr,"\t-y, --ymin #\tMinimum limit of region in y[cm]\n");
	fprintf(stderr,"\t-Y, --ymax #\tMaximum limit of region in y[cm]\n");
	fprintf(stderr,"\t-z, --zmin #\tMinimum limit of region in z[cm]\n");
	fprintf(stderr,"\t-Z, --zmax #\tMaximum limit of region in z[cm]\n");
	//fprintf(stderr,"\t    --xrebin #\tRebinx axis by #\n");
	//fprintf(stderr,"\t    --yrebin #\tRebinx axis by #\n");
	//fprintf(stderr,"\t    --zrebin #\tRebinx axis by #\n");
	fprintf(stderr,"\t-O, --old\tWrite old format\n");
	fprintf(stderr,"\t-o, --output file\tFile with the pos,FWHM/hor,vert max, err along z-slices\n");
	fprintf(stderr,"\tinfile\tInput file from the USBSUW output\n");
	fprintf(stderr,"author: Vasilis.Vlachoudis@cern.ch (2002)\n");
	exit(0);
} /* usage */

/* --- main --- */
int main(int ac, char *av[])
{
	FILE	*fin, *fout;
	int	i, j, k, m, n, length, Ndet;
	char	runtitle[RUNTIT+1],
		runtime[RUNTIM+1];
	char	buffer[512], *pbuffer;
	int	dataSize;
	int	detector2Look = -1;
	int	old = 0;

	int	c;

	prgname = av[0];

	fin  = stdin;
	fout = stdout;

	for (j=0; j<3; j++) {
		scanMin[j] = -FLT_MAX;
		scanMax[j] =  FLT_MAX;
		scMin[j]   =  INT_MIN;
		scMax[j]   =  INT_MAX;
		rebin[j]   = 1;
	}

	while (1) {
		int option_index = 0;
		static struct option long_options[] = {
			{"help",   0, 0, 'h'},
			{"detector",1,0, 'd'},
			{"axis",   1, 0, 'a'},
			{"xmin",   1, 0, 'x'},
			{"xmax",   1, 0, 'X'},
			{"ymin",   1, 0, 'y'},
			{"ymax",   1, 0, 'Y'},
			{"zmin",   1, 0, 'z'},
			{"zmax",   1, 0, 'Z'},
			{"old",    0, 0, 'O'},
			{"output", 1, 0, 'o'},
			{"xrebin", 1, 0, 'i'},
			{"yrebin", 1, 0, 'j'},
			{"zrebin", 1, 0, 'k'},
			{0, 0, 0, 0}
		};

		c = getopt_long(ac, av, "h?Oa:d:x:X:y:Y:z:Z:o:i:j:k:",
			long_options, &option_index);
		if (c == -1)
			break;

		switch (c) {
			case 0:
				printf("option %s", long_options[option_index].name);
				if (optarg)
					printf(" with arg %s", optarg);
				printf("\n");
				usage();
				break;

			case '?':
			case 'h':
				usage();
				break;

			case 'a':
				if (optarg[0] == 'x' || optarg[0]=='X' ||
				    optarg[0] == 'r' || optarg[0]=='R')
					axis = X;
				else
				if (optarg[0] == 'y' || optarg[0]=='Y' ||
				    optarg[0] == 'p' || optarg[0]=='P')
					axis = Y;
				else
				if (optarg[0] == 'z' || optarg[0]=='Z')
					axis = Z;
				else {
					fprintf(stderr,"Invalid axis %c",optarg[0]);
					axis = 2;
				}
				axis1 = (axis+1)%3;
				axis2 = (axis+2)%3;
				break;

			case 'd':
				detector2Look = atoi(optarg);
				break;

			case 'x':
				scanMin[X] = atof(optarg);
				break;
			case 'X':
				scanMax[X] = atof(optarg);
				break;

			case 'y':
				scanMin[Y] = atof(optarg);
				break;
			case 'Y':
				scanMax[Y] = atof(optarg);
				break;

			case 'z':
				scanMin[Z] = atof(optarg);
				break;
			case 'Z':
				scanMax[Z] = atof(optarg);
				break;

			case 'i':
				rebin[X] = atoi(optarg);
				break;
			case 'j':
				rebin[Y] = atoi(optarg);
				break;
			case 'k':
				rebin[Z] = atoi(optarg);
				break;

			case 'O':
				old = 1;
				break;

			case 'o':
				if ((fout=fopen(optarg,"w"))==NULL) {
					fprintf(stderr,"Error opening output file %s\n",optarg);
					return 3;
				}
				break;

			default:
				printf("?? getopt returned character code 0%o ??\n", c);
				usage();
		}
	}
	if (optind < ac) {
		if ((fin=fopen(av[optind],"rb"))==NULL) {
			fprintf(stderr,"Error opening input file %s\n",av[2]);
			return 2;
		}
	}

	if (detector2Look<=0 || detector2Look>100) {
		if (ac>1)
			fprintf(stderr,"Invalid detector number %s\n",av[1]);
		usage();
	}
	detector2Look--;

	// First buffer WRITE (1) RUNTIT*80,RUNTIM*32,SNGL(WCTOT),NCTOT
	length = fortranRead(fin,buffer,sizeof(buffer)); pbuffer = buffer;
	FTNSGET(runtitle,RUNTIT,pbuffer);
	FTNSGET(runtime,RUNTIM,pbuffer);
	FTNGET(wctot,float,pbuffer);
	FTNGET(nctot,int,pbuffer);

	printf("Title: %s\n",runtitle);
	printf("Time: %s\n",runtime);
#if __DEBUG
	printf("WCTOT: %g\n",wctot);
	printf("NCTOT: %d\n",nctot);
#endif

	for (i=0; i<MAXDET; i++) {
		/* Read detector information
			WRITE (1) NB,TITUSB(NB),ITUSBN(NB),IDUSBN(NB),
				XLOW(NB),XHIGH(NB),NXBIN(NB),DXUSBN(NB),
				YLOW(NB),YHIGH(NB),NYBIN(NB),DYUSBN(NB),
				ZLOW(NB),ZHIGH(NB),NZBIN(NB),DZUSBN(NB),
				LNTZER(NB),BKUSBN(NB),B2USBN(NB),TCUSBN(NB)
		*/
		length = fortranRead(fin,buffer,sizeof(buffer));
		if (length!=86) {
			Ndet = i;
			break;
		}
		if (length<0) {
			fprintf(stderr,"Error reading input file %s\n",av[2]);
			return -1;
		}

		if (detector2Look != i) {
			fortranSkip(fin);
			continue;
		}

		/* parse buffer */
		pbuffer = buffer;
#if __DEBUG
		printf("L=%d (%x)\n",length,length);
		hexdump(buffer,length);
#endif

		FTNGET(nb,	int,	pbuffer);
		FTNSGET(titusb,TITUSB,pbuffer);
		FTNGET(itusbn,int,	pbuffer);
		FTNGET(idusbn,int,	pbuffer);
		for (j=0; j<3; j++) {
			FTNGET(low[j],	float,	pbuffer);
			FTNGET(high[j],	float,	pbuffer);
			FTNGET(nbin[j],	int,	pbuffer);
			FTNGET(dd[j],	float,	pbuffer);
		}
		FTNGET(lntzer,int,	pbuffer);
		FTNGET(bkusbn,float,	pbuffer);
		FTNGET(b2usbn,float,	pbuffer);
		FTNGET(tcusbn,float,	pbuffer);

		//printf("LNTZER %d\n",lntzer);
		//printf("BKUSBN %g\n",bkusbn);
		//printf("B2USBN %g\n",b2usbn);
		//printf("TCUSBN %g\n",tcusbn);

		printf("---> Processing detector %d\n",i+1);
		for (j=0; j<3; j++) {
			if (scanMin[j] <  low[j])  scanMin[j] = low[j];
			if (scanMax[j] >= high[j]) scanMax[j] = high[j];
			scMin[j] = (int)((scanMin[j]-low[j])/dd[j] + 0.5);
			scMax[j] = (int)((scanMax[j]-low[j])/dd[j] + 0.5);
			if (scMin[j] < 0)	  scMin[j] = 0;
			if (scMax[j] >= nbin[j])  scMax[j] = nbin[j];
			if (scMax[j] <= scMin[j]) scMax[j] = scMin[j] + 1;
		}

		printf("\nDetector: %d %s\n",i+1, titusb);
		printf("Score Type=%d Particle=%d\n", itusbn, idusbn);
		printf("Scanning axis: %c\n",XYZ[axis]);
		printf("    Hor. axis: %c\n",XYZ[axis1]);
		printf("    Ver. axis: %c\n",XYZ[axis2]);
		//printf("Rebinning:\n");
		//for (j=0; j<3; j++)
		//	printf("\t%c: %d",XYZ[j], rebin[j]);
		//
		for (j=0; j<3; j++)
			printf("%c: %g\t%g\t%d\t%g\n", XYZ[j],
				low[j], high[j], nbin[j], dd[j]);

		printf("Looking at region\n");
		for (j=0; j<3; j++)
			printf("%c: %g - %g\t[%d .. %d]\n", XYZ[j],
				scanMin[j], scanMax[j],
				scMin[j]+1, scMax[j]);

		dataSize  = nbin[X]*nbin[Y]*nbin[Z]*sizeof(float);
		int planeSize = nbin[axis1]*nbin[axis2]*sizeof(float);

		data  = (float*)malloc(dataSize);
		plane = (float*)malloc(planeSize);

		slice = (SliceData*)malloc(nbin[axis]*sizeof(SliceData));
		memset(slice, 0, nbin[axis]*sizeof(SliceData));

		// Read the entire array
		if (fortranRead(fin, data, dataSize)<0) return -1;

		if (old) {
			// Write the results
			fprintf(fout,"#Title: %s\n",runtitle);
			fprintf(fout,"#Time: %s\n",runtime);
			fprintf(fout,"#Detector: %d %s\n",nb,titusb);
			fprintf(fout,"#Score Type=%d Particle=%d\n",itusbn,idusbn);
			fprintf(fout,"#Scanning axis: %c\n",XYZ[axis]);
			fprintf(fout,"#    Hor. axis: %c\n",XYZ[axis1]);
			fprintf(fout,"#    Ver. axis: %c\n",XYZ[axis2]);
			//fprintf(fout,"#Rebinning:\n");
			//for (j=0; j<3; j++)
			//	fprintf(fout,"#  %c: %d",XYZ[j], rebin[j]);
			fprintf(fout,"#Limits:\n");
			for (j=0; j<3; j++)
				fprintf(fout,"#%c: %g\t%g\t%d\t%g\n", XYZ[j],
						low[j], high[j], nbin[j], dd[j]);
			fprintf(fout,"#Looking at region\n");
			for (j=0; j<3; j++)
				fprintf(fout,"#%c: %g .. %g\t[%d .. %d]\n", XYZ[j],
					scanMin[j], scanMax[j], scMin[j]+1, scMax[j]);
			fprintf(fout,"# X(1) Y(2) Z(3) HL(4) HH(5) VL(6) VH(7) Max(8) Err%%(9)\n");
		}

		// split it into slices
		for (k=0; k<nbin[axis]; k++) {
			if (k<scMin[axis] || k>=scMax[axis]) continue;
			for (m=0; m<nbin[axis2]; m++)
				for (n=0; n<nbin[axis1]; n++)
					switch (axis) {
						case X:	// n=Y(axis1), m=Z(axis2)
							plane[m*nbin[axis1] + n]
								= data[m*nbin[X]*nbin[Y]
								     + n*nbin[X]
								     + k];
							break;
						case Y:	// n=Z(axis1), m=X(axis2)
							plane[m*nbin[axis1] + n]
								= data[n*nbin[X]*nbin[Y]
								     + k*nbin[X]
								     + m];
							break;
						case Z:	// n=X(axis1), m=Y(axis2)
							plane[m*nbin[axis1] + n]
								= data[k*nbin[X]*nbin[Y]
								     + m*nbin[X]
								     + n];
							break;
					}
			scanPlane(k);
		}
	}

	if (length==14) { // The statistics are following
		printf("\nStatistics are present\n");
		for (i=0; i<Ndet; i++) {
			if (detector2Look != i) {
				fortranSkip(fin);
				continue;
			}
			// Read the entire array
			if (fortranRead(fin, data, dataSize)<0) return -1;

			for (k=0; k<nbin[axis]; k++)
				slice[k].err = data[ slice[k].pos[Z]*nbin[X]*nbin[Y]
						   + slice[k].pos[Y]*nbin[X]
						   + slice[k].pos[X] ];
		}
	}

	// Write tables
	if (old) {
		for (k=0; k<nbin[axis]; k++)
			fprintf(fout,"%lg %lg %lg %lg %lg %g %g %lg %g\n",
				MPOS(slice[k].pos[X],X),
				MPOS(slice[k].pos[Y],Y),
				MPOS(slice[k].pos[Z],Z),
				MPOS(slice[k].Hmin,axis1), MPOS(slice[k].Hmax,axis1),
				MPOS(slice[k].Vmin,axis2), MPOS(slice[k].Vmax,axis2),
				slice[k].max, slice[k].err*100.0);
	} else {
		// 1. Maximum trace
		fprintf(fout,"# Detector n: Max %s\n",titusb);
		fprintf(fout,"# Time: %s\n",runtime);
		fprintf(fout,"# Score Type=%d Particle=%d\n",itusbn,idusbn);
		fprintf(fout,"# Scanning axis: %c\n",XYZ[axis]);
		fprintf(fout,"#     Hor. axis: %c\n",XYZ[axis1]);
		fprintf(fout,"#     Ver. axis: %c\n",XYZ[axis2]);
		//fprintf(fout,"# Rebinning:\n");
		//for (j=0; j<3; j++)
		//	fprintf(fout,"#  %c: %d",XYZ[j], rebin[j]);
		fprintf(fout,"# Limits:\n");
		for (j=0; j<3; j++)
			fprintf(fout,"#%c: %g\t%g\t%d\t%g\n", XYZ[j],
					low[j], high[j], nbin[j], dd[j]);
		fprintf(fout,"# Looking at region\n");
		for (j=0; j<3; j++)
			fprintf(fout,"#%c: %g .. %g\t[%d .. %d]\n", XYZ[j],
				scanMin[j], scanMax[j], scMin[j]+1, scMax[j]);
		for (k=0; k<nbin[axis]; k++) {
			if (k<scMin[axis] || k>=scMax[axis]) continue;
			fprintf(fout,"%lg %lg %lg %lg\n",
				LOPOS(slice[k].pos[axis],axis),
				HIPOS(slice[k].pos[axis],axis),
				slice[k].max, slice[k].err*100.0);
		}

		fprintf(fout,"\n\n# Detector n: Horiz %s\n",titusb);
		for (k=0; k<nbin[axis]; k++) {
			if (k<scMin[axis] || k>=scMax[axis]) continue;
			fprintf(fout,"%lg %lg %lg %lg %lg\n",
				LOPOS(slice[k].pos[axis],axis),
				HIPOS(slice[k].pos[axis],axis),
				MPOS(slice[k].pos[axis1],axis1),
				MPOS(slice[k].Hmin,axis1),
				MPOS(slice[k].Hmax,axis1));
		}

		fprintf(fout,"\n\n# Detector n: Vert %s\n",titusb);
		for (k=0; k<nbin[axis]; k++) {
			if (k<scMin[axis] || k>=scMax[axis]) continue;
			fprintf(fout,"%lg %lg %lg %lg %lg\n",
				LOPOS(slice[k].pos[axis],axis),
				HIPOS(slice[k].pos[axis],axis),
				MPOS(slice[k].pos[axis2],axis2),
				MPOS(slice[k].Vmin,axis2),
				MPOS(slice[k].Vmax,axis2));
		}
	}

	// free data, plane, slice!
	if (data) free(data);
	if (plane) free(plane);
	if (slice) free(slice);

	fclose(fin);
	fclose(fout);
	return 0;
} /* main */