#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
#include <stddef.h>
#define ABS(n) ((n>0)?n:-n)


void help();

main(argc, argv)
int 	argc;
char 	*argv[];

{


/* program description */

/* The program reads an ASCII file (.dat) and creates
   images from this file. Three images can be created:
	-Total : image of primary and scattered photons.
	-Scatter : image of scattered photons. 
	-Primary : image of primary photons. 

   The program takes 1 argument: the original ASCII filename./*



/* parameter declaration */
//____________________________________________
float pro1=30;
float pro2=2*pro1;
float pro3=3*pro1;
float pro4=4*pro1;
float degree=90/pro1;
int nbima=120;        /* number of images */
int nx = 200;
int ny = 140;
float pixsize = 1;
//____________________________________________

char *s, count[17];
int nbpix;
float xxx ;
char file_name[4];
int *OrderVol, *OrderCrystal , *int1, *NumberRun, *NumberTete;
float *Energy, *float1,*Coordonnees;
char *chars, *chars2;
int *Primary, *Total, *Scatter;
int ncpttotal=0, ncptscatter=0, ncptprimary=0;
char c;
int CrystalOrder, Order;
char caractere1, caractere2, caractere3, caractere4;
int  A, Xx, Yy, neven, Run, Tete, *ima2, ind, e;
short int outex; /*=1 : existing file, = 0 : the output file does not exist */
short int vhelp;
float boundxp, boundxn, boundyp, boundyn, Ener, alpha, beta, gamma, x, y, a, b, cosgamma, absgamma;
long int i, j;
float xip, yip, ip, E, Efin, Einit, Einit2, fov, fovx ,fovy, pasener, Y, X;
FILE *fin,  *fout;


/* parameter initialisation */

outex=1;
pasener=0;
neven=1;


	while (--argc > 0 && (*++argv)[0] == '-')
		for (s = argv[0]+1; *s!= '\0'; s++)
			switch (*s) {
			case 'h':
				help();
				vhelp=1;
				return 1;
				break;
			default:
				printf ("illegal option %c\n", *s);
				argc = 0;
				break;

			}

	if ((argc != 1)&&(vhelp==0)) {
		printf ("usage : benchmark_projections -h filin \n");
		printf ("type benchmark_projections -h for more help \n");
		return 1;
	}






/* open the files */

fin = fopen(*argv,"r");
if(fin==NULL)
	{
		printf("error when opening fin\n");
		return;
	}




/* Memory allocation */

OrderVol = (int *) calloc(neven, sizeof(int));
OrderCrystal = (int *) calloc(neven, sizeof(int));
Energy = (float *) calloc(neven, sizeof(float));
chars = (char *) calloc(neven*4, sizeof(char));
chars2 = (char *) calloc(neven*4, sizeof(char));
int1 = (int *) calloc(1, sizeof(int));
float1 = (float *) calloc(1, sizeof(float));
NumberRun = (int *) calloc(neven, sizeof(int));
NumberTete = (int *) calloc(neven, sizeof(int));
Coordonnees = (float *) calloc(neven*3, sizeof(float));



for(i=0;i<neven;i++)
  {
    OrderVol[i]==0;
    OrderCrystal[i]==0;
    Energy[i]==0;
    //int1[i]==0;
    //float1[i]==0;
    //NumberRun[i]==0;
    //NumberTete[i]==0;
    //for(j=0;i<3;j++)
    //{
    //Coordonnees[j]==0;
    //}
   }
for(i=0;i<neven*4;i++)
  {
    chars[i]=='Z';
  }

for(i=0;i<neven*4;i++)
  {
    chars2[i]=='Z';
  }


/* Threshold Input  */

printf("Enter the lower energy threshold (MeV)\n");
scanf("%f",&Einit);
printf("\n");
printf("Lower energy threshold is %f MeV\n",Einit);
printf("\n");

printf("Enter the upper energy threshold (MeV)\n");
scanf("%f",&Einit2);
printf("\n");
printf("Upper energy threshold is %f MeV\n",Einit2);
printf("\n");



/* computation of the parameters */

pasener=Einit2-Einit;
nbpix=nx*ny;
boundxp=(pixsize*nx)/2;
boundxn=-(pixsize*nx)/2;
boundyp=(pixsize*ny)/2;
boundyn=-(pixsize*ny)/2;
Efin=Einit2;



/* memory allocation */

Primary = (int *) calloc(nbpix*nbima, sizeof(int));
Total = (int *) calloc(nbpix*nbima, sizeof(int));
Scatter = (int *) calloc(nbpix*nbima, sizeof(int));


/*Initialisation  of the output variables */


for(i=0 ; i<nbpix*nbima; i++)
  {
    Primary[i]=0;
  }

for(i=0 ; i<nbpix*nbima; i++)
  {
    Total[i]=0;
  }

for(i=0 ; i<nbpix*nbima; i++)
  {
    Scatter[i]=0;
  }



/*Processing*/

printf("\n\nBeginning of processing\n");


i=0;
fscanf(fin,"%d\n", &NumberRun[i]);

while(!feof(fin))
 {
     for (i=0 ; i<neven ; i++) 
       {
      fscanf(fin,"%d %d\n", &int1[0],&int1[0]);
      fscanf(fin,"%f %f %f\n", &float1[0],&float1[0],&float1[0]);
      fscanf(fin,"%d\n", &NumberTete[i]);
      fscanf(fin,"%d %d\n", &int1[0],&int1[0]);
      fscanf(fin,"%f\n", &float1[0]);
      fscanf(fin,"%f\n", &Energy[i]);
      fscanf(fin,"%f %f %f\n", &Coordonnees[i],&Coordonnees[i+1],&Coordonnees[i+2]);
      fscanf(fin,"%d\n", &OrderVol[i]);
      fscanf(fin,"%d\n", &OrderCrystal[i]);
      fscanf(fin,"%d\n", &int1[0]);
      fscanf(fin,"%d\n", &int1[0]);
      fscanf(fin,"%c\n", &chars[i]);
      fscanf(fin,"%c\n", &chars[i+1]);
      fscanf(fin,"%c\n", &chars[i+2]);
      fscanf(fin,"%c\n", &chars[i+3]);
      if(chars[i+3] != 'L')
	{
	  c = fgetc(fin);
	  while(c != '_')
	    {
	      c = fgetc(fin);
	    }
	  c = fgetc(fin);

	}
if(chars[i+3] == 'L')
  {c='Z';}
if(c == 'P')
  c = fgetc(fin);
      fscanf(fin,"%c\n", &chars2[i]);
      fscanf(fin,"%c\n", &chars2[i+1]);
      fscanf(fin,"%c\n", &chars2[i+2]);
      fscanf(fin,"%c\n", &chars2[i+3]);
      if(chars2[i+3] != 'L')
        {
          c = fgetc(fin);
          while(c != '_')
            {
              c = fgetc(fin);
            }
          c = fgetc(fin);
	  }

	 Run=NumberRun[i];
	 Tete=NumberTete[i];
	 E=Energy[i];
	 ip=Coordonnees[i]; //read x
	 xip=Coordonnees[i+1];//read y
	 yip=Coordonnees[i+2];//read z
	 caractere1=chars[i];
	 caractere2=chars[i+1];
	 caractere3=chars[i+2];
	 caractere4=chars[i+3];
	 Order=OrderVol[i];
	 CrystalOrder=OrderCrystal[i];
				if(Tete==0){A=Run;}
				if(Tete==1){A=pro1+Run;}
				if(Tete==2){A=pro2+Run;}
				if(Tete==3){A=pro3+Run;}

				if(A==0)
				{
				yip=yip;
				xxx=xip;
				}
				if(A==pro1)
				{
				yip=yip;
				xxx=(-1)*ip;
				}
				if(A==pro2)
				{
				xxx=xip*(-1);
				}
				if(A==pro3)
				{
				xxx=ip;
				}

				if(A > 0 && A < pro1 && xip < 0 && ip >0)    // ----------  1.
				{
				  alpha = ((degree*A)/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(x/a);
				  gamma =alpha+beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  xxx=(-1)*b;
				}
				

				if(A > 0 && A < pro1 && ip < 0 && xip >0)     // ----------  2.
				{
				  alpha = ((degree*A)/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(y/a);  
				  gamma =(1.57-alpha)+beta;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  xxx=b;
				}



				if(A > 0 && A < pro1 && xip >= 0 && ip >= 0) // ----------  3.
				{				  
				  alpha = ((degree*A)/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(x/a);
				  gamma =alpha-beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  if(xip != 0 && ip != 0)
				    {
				      if(gamma > 0)
					{
					  xxx=(-1)*b;
					}
				      if(gamma < 0)
					{
					  xxx=b;
					}
				    }
				  if(xip == 0)
				    {

				      xxx=(-1)*b;
				    }
				  if(ip == 0)
				    {
				      xxx=b;
				    } 
				}
				


				if(A > pro1 && A < pro2 && ip > 0 && xip >0)   // ----------  4.
				{
				alpha = ((degree*(A-pro1))/180)*3.1415;
				x=ABS(ip);
				y=ABS(xip);
				a=sqrt(x*x+y*y);
				beta=acos(y/a);
				gamma =alpha+beta ;
				absgamma = ABS(gamma);
				b = a*(sin(absgamma));
				xxx=(-1)*b;
				}



				if(A > pro1 && A < pro2 && ip < 0 && xip < 0)   // ----------  5.
				{
				alpha = ((degree*(A-pro1))/180)*3.1415;
				x=ABS(ip);
				y=ABS(xip);
				a=sqrt(x*x+y*y);
				beta=acos(x/a);
				gamma =(1.57-alpha)+beta ;
				absgamma = ABS(gamma);
				b = a*(sin(absgamma));
				xxx=b;
				}




				if(A > pro1 && A < pro2 && xip >= 0 && ip <= 0)    // ----------  6.
				{				  
				  alpha = ((degree*(A-pro1))/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(y/a);
				  gamma =alpha-beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  if(xip != 0 && ip != 0)
				    {
				      if(gamma > 0)
					{
					  xxx=(-1)*b;
					}
				      if(gamma < 0)
					{
					  xxx=b;
					}
				    }
				  if(xip == 0)
				    {
				      xxx=b;
				    }
				  if(ip == 0)
				    {
				      xxx=(-1)*b;
				    }
				}


				if(A > pro2 && A < pro3 && xip > 0 && ip < 0)   // ----------  7.
				{
				alpha = ((degree*(A-pro2))/180)*3.1415;
				x=ABS(ip);
				y=ABS(xip);				
				a=sqrt(x*x+y*y);
				beta=acos(x/a);
				gamma =alpha+beta ;
				absgamma = ABS(gamma);
				b = a*(sin(absgamma));
				xxx=((-1)*b);				
				}

				if(A > pro2 && A < pro3 && ip > 0 && xip <0)   // ----------  8.
				{
				alpha = ((degree*(A-pro2))/180)*3.1415;
				x=ABS(ip);
				y=ABS(xip);
				a=sqrt(x*x+y*y);
				beta=acos(y/a);
				gamma =(1.57-alpha)+beta ;
				absgamma = ABS(gamma);
				b = a*(sin(absgamma));
				xxx=b;
				}

				if(A > pro2 && A < pro3 && xip <= 0 && ip <= 0)    // ----------  9.
				{				  
				  alpha = ((degree*(A-pro2))/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(x/a);
				  gamma =alpha-beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));

				  if(xip != 0 && ip != 0)
				    {
				      if(gamma > 0)
					{
					  xxx=(-1)*b;
					}
				      if(gamma < 0)
					{
					  xxx=b;
					}
				    }
				  if(xip == 0)
				    {
				      xxx=(-1)*b;
				    }
				  if(ip == 0)
				    {
				      xxx=b;
				    }
				}




				if(A > pro3 && A < pro4 && ip < 0 && xip < 0)   // ----------  10.
				{
				  alpha = ((degree*(A-pro3))/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(y/a);
				  gamma =alpha+beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  xxx=(-1)*b;
				}



				if(A > pro3 && A < pro4 && xip > 0 && ip > 0)   // ----------  11.

				{
				  alpha = ((degree*(A-pro3))/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);
				  beta=acos(x/a); 
				  gamma =(1.57-alpha)+beta ;
				  absgamma = ABS(gamma);
				  b = a*(sin(absgamma));
				  xxx=b;
				}


				if(A > pro3 && A < pro4 && ip >= 0 && xip <= 0)    // ----------  12.

				{
				  alpha = ((degree*(A-pro3))/180)*3.1415;
				  x=ABS(ip);
				  y=ABS(xip);
				  a=sqrt(x*x+y*y);		   
				  beta=acos(y/a);			 
				  gamma =alpha-beta ;			       
				  absgamma = ABS(gamma);				  
				  b = a*(sin(absgamma));				  
				  if(ip != 0 && xip != 0)
				    {
				      if(gamma > 0)
					{				     
					  xxx=(-1)*b;
					}
				      if(gamma < 0)
					{     
					  xxx=b;
					}
				    }
				  if(ip == 0)
				    {
				      xxx=(-1)*b;
				    }
				  if(xip == 0)
				    {
				      xxx=b;
				    }

				}

				if(E>=Einit && E<=Efin)
				  {
				    if ((xxx>=boundxn)&&(xxx<=boundxp)&&(yip>=boundyn)&&(yip<=boundyp))
				      {       
					Y=(yip+boundyp)/pixsize;
					if(Y<0){Y=(-1)*Y;}
					Yy=(floor)(Y);
					X=(xxx+boundxp)/pixsize;
					if(X<0){X=(-1)*X;}
					Xx=(floor)(X);
					ind=(int)(A*nbpix+nx*Yy+Xx);

					if(caractere1=='N' && caractere2=='U' &&  caractere3=='L' && caractere4=='L' && Order==0)
					  {
					    Primary[ind]++;
					    ncptprimary++;
					  }
					if((caractere1!='N' && caractere2!='U' &&  caractere3!='L' && caractere4!='L') || Order>0)
					  {
					    Scatter[ind]++;
					    ncptscatter++;
					  }
					Total[ind]++;
					ncpttotal++;
				      }
				  }
   

fscanf(fin,"%d\n", &NumberRun[i]);


       }
     }

printf("\n\nEnd of processing\n\n");


printf("Number of detected events in the image of primary and scattered photons = %d\n", ncpttotal);
printf("Number of primary events in the image of primary photons= %d\n", ncptprimary);
printf("Number of scattered events  in the image of scattered photons = %d\n", ncptscatter);




/* Writing the results in fout */
sprintf(file_name,"%s","Total");
printf("\nFile Total correponds to the image of primary and scattered photons");
if((fout = fopen(file_name,"wb")) ==NULL)
  {
    printf(" error opening output file%s\n") ;
  }
fwrite(&Total[0],sizeof(int),nbpix*nbima,fout);
fclose(fout);


sprintf(file_name,"%s","Primary");
printf("\nFile Primary correponds to the image of primary photons");
if((fout = fopen(file_name,"wb")) ==NULL)
  {
    printf(" error opening output file%s\n") ;
  }
fwrite(&Primary[0],sizeof(int),nbpix*nbima,fout);
fclose(fout);


sprintf(file_name,"%s","Scatter");
printf("\nFile Scatter correponds to the image of scattered photons\n\n");
if((fout = fopen(file_name,"wb")) ==NULL)
  {
    printf(" error opening output file%s\n") ;
  }
fwrite(&Scatter[0],sizeof(int),nbpix*nbima,fout);
fclose(fout);


	fclose(fin);
   



}



void help()

{
	printf ("usage : benchmark_projections -h filin \n");
	printf ("    -h : help \n");
	printf ("    filin : input file (gateSingles.dat)\n");
	printf ("\n  The program reads an ASCII file (.dat) and creates images from this file. Three images can be created:\n");
	printf ("    -Total : image of primary and scattered photons.\n");
	printf ("    -Scatter : image of scattered photons.\n");
	printf ("    -Primary : image of primary photons. \n");
	printf ("The program takes 1 argument: the original ASCII filename.\n");

}