/* compiles by `gcc grid.c -I. -L. -lrdmc -lm -o grid' */

#include <rdmc.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#define VERSION "0.01"
#define MAX_IN 256
#define SPEEDOFLIGHT 0.299792458

int iom=339;
int minhit=10;
int fn=2;

main(int arg_c, char *arg_a[]){
  char geomname[MAX_IN]="geom";
  FILE *fg;
  mcfile *fi;
  array ar;
  mevt ev;
  mtrack mt;
  int i, j, flag=0, geomfix=0, length, error=0, num;
  double x, y, z, t, xom, yom, zom, sum, amp, amx, tin;


  if(arg_c==1){
    printf("This creates datafiles for muon geometry calibration. Use `%s -h' for help\n", arg_a[0]);
    return 1;
  }

  for(i=1;i<arg_c;i++){
    if(*arg_a[i]=='-'){
      if(!(strcmp(arg_a[i], "-h") && strcmp(arg_a[i], "-help") && strcmp(arg_a[i], "--help"))){
	printf("\n\
This program computes the grid values of the OM contribution to the likelihood\n\
Available options are: -h  = -help  print this message and quit\n\
                       -om=[OM number] selects OM to calibrate\n\
                       -geom=[file name of the geometry file]\n\
                       -fn=[number of the fitted track to use]\n\
\n\
"); return 1;
      }

      else if(!strncmp(arg_a[i], "-om=", 4)){
	iom=atoi(arg_a[i]+4)-1;
      }

      else if(!strncmp(arg_a[i], "-geom=", 6)){
	length=strlen(arg_a[i]+6);
	if(length>=MAX_IN || length<=0){
	  fprintf(stderr, "Incorrect usage of `-geom=[geom file]' option\n");
	  return 1;
	}
	strcpy(geomname, arg_a[i]+6);
	geomfix=1;
      }

      else if(!strncmp(arg_a[i], "-fn=", 4)){
	fn=atoi(arg_a[i]+4)-1;
      }

      else{
	fprintf(stderr, "Unrecognized option `%s', exiting\n", arg_a[i]);
	return 1;
      }
    }
    else flag=1;
  }

  if(flag){
    fprintf(stderr, "attempting to do OM %i with fit result #%i\n", iom+1, fn+1);

    for(i=1; i<arg_c; i++) if(*arg_a[i]!='-'){

      if((fi=rdmc_mcopen(arg_a[i], "r", F2000_F))==NULL){
	fprintf(stderr, "cannot open input file `%s', exiting\n", arg_a[i]);
	return 1;
      }

      rdmc_rarr(fi, &ar);
      if(ar.n_fit<=fn || fn<0){
	if(ar.n_fit>0) fprintf(stderr, "fn must be one of 1 ... %i\n", ar.n_fit);
	else fprintf(stderr, "No fits defined, cannot continue\n");
	return 1;
      }

      rdmc_init_mevt(&ev);

      if(geomfix){
	if((fg=fopen(geomname,"r"))==NULL){
	  fprintf(stderr, "cannot open output file `%s', exiting\n", geomname);
	  fclose(fg);
	  return 1;
	}

	while(1){
	  error=fscanf(fg, "%i %lg %lg %lg", &i, &x, &y, &z);
	  if(error==4){
	    ar.x[i-1]=x;
	    ar.y[i-1]=y;
	    ar.z[i-1]=z;
	  }
	  else if(error==EOF){
	    fprintf(stderr, "EOF in geom file `%s' reached   ", geomname);
	    fclose(fg);
	    break;
	  }
	  else{
	    fprintf(stderr, "Format error %i in geom file `%s'\n", error, geomname);
	    fclose(fg);
	    return 1;
	  }
	}
      }

      xom=ar.x[iom];
      yom=ar.y[iom];
      zom=ar.z[iom];
      fprintf(stderr, "default coordinates of OM %i are x=%g y=%g z=%g\n", iom+1, xom, yom, zom);

      while(1){
	error=rdmc_revt(fi, &ar, &ev);
	if(error==EOF){
	  fprintf(stderr, "EOF reached\n");
	  break;
	}
	if(error==RDMC_ILF){
	  fprintf(stderr, "File format error\n");
	  break;
	}
	if(error){
	  fprintf(stderr, "Unknown error %i\n", error);
	  break;
	}

	if(ev.nhits<minhit) continue;

	mt=ev.rec[fn];
	mt=ev.gen[1];
	num=0;
	amx=0;
	sum=0;
	tin=0;

	for(j=0;j<ev.nhits;j++) if(ev.h[j].ch==iom){
	  num++; amp=ev.h[j].amp>0?ev.h[j].amp:1; sum+=amp;
	  if(amx==0){ amx=amp; tin=ev.h[j].t; }
	}
	t=(xom-mt.x)*mt.px+(yom-mt.y)*mt.py+(zom-mt.z)*mt.pz;
	x=mt.x+mt.px*t-xom;
	y=mt.y+mt.py*t-yom;
	z=mt.z+mt.pz*t-zom;
	tin-=mt.t+t/SPEEDOFLIGHT;

	int Nall, Ndir;
	float cogX, cogY, cogZ, Rwei, Rave, Rmax, Ldir, Shnh;

	rdmc_get_direct_hits(&ar, &mt, &ev,
     "-100.:1000.;c -15.:25.;c -15.:75.;x -15.:75.;y -15.:75.;z -15.:75.;R -15.:75.;D -15.:75.;M -15.:75.;V",
			     &Nall, &Ndir, &cogX, &cogY, &cogZ, &Rave, &Rmax, &Ldir, &Shnh);

	t=(cogX-mt.x)*mt.px+(cogY-mt.y)*mt.py+(cogZ-mt.z)*mt.pz;
	cogX-=mt.x+mt.px*t;
	cogY-=mt.y+mt.py*t;
	cogZ-=mt.z+mt.pz*t;
	Rwei=sqrt(cogX*cogX+cogY*cogY+cogZ*cogZ);

	printf("%i %g %g %g %g %g %g %g %g", num, sum, x, y, z, mt.px, mt.py, mt.pz, tin);
	printf("  %i %i %g %g %g %g %g %g\n", Nall, Ndir, Rwei, Rave, Rmax, Ldir, Shnh, ev.fresult[fn].val[1]);
      }

      rdmc_mcclose(fi);
    }
  }
}