#include <stdio.h>
#include <string.h>
//#include <conio.h>
#include <math.h>
#include <iostream.h>
#include "distrib.h"// contains the uniform, Poisson, Exponential, Normal distributions realizations
double const pi=3.14159265;
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//Parameters taken from configuration file config.txt
 double Z_Position;
 double OrientationAngle;
 int    FiberNumber;
 double FiberDiameter;
 double Pitch;
 double PhotonYield;
 int    FiberBundling;
 double QuantumEfficiency;
 double PMT_Gain;
 int GainDistribution;
 double GainSpread;
 int RandomGen;
 double Delta_Z;
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
const int a=6,b=11;//define rows to be read from G4BL file: read rows only from 6 to 11
double outputdata[4096];// output file data array
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
double GetY(double y0,double z, double *pt)
{  double y;
	if(*(pt+4)>=0)
		{
		   y=y0+z**(pt+25);
		   return y;
		}
		else
		{
		   y=y0-z**(pt+25);
		   return y;
		}
}
double GetX(double x0, double z, double *pt)
{
	double x;
	double tanfi=*(pt+32)/*tangent fi*/;
	if(*(pt+3)>=0) x=x0+z*tanfi;
	else           x=x0-z*tanfi;
	return x;
}
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////

void SearchPosition(FILE *fp)//analizes G4BL file and searches data corresponding to concrete detector
{                            //which is determined by Z_Position value

int i,k=0,t=1;
char str[200],str1[20];
double d;
      fgets(str,200,fp); //missing first three header non-data strings
      fgets(str,200,fp);
      fgets(str,200,fp);
	 while(t)
	{
		for(i=0;i<a+3;i++)
		{
			fscanf(fp,"%s",str1);
		}
			d=atof(str1);
			if(d!=Z_Position)
			{
				fgets(str,200,fp);//finds the necessary line, but now pointer is shifted on the whole line ahead
				k++;                  // we need to rewind file and place pointer in the right position
			}
			else break;
	}
	rewind(fp);// places the pointer in the beginning of file
	fgets(str,200,fp);// missing first three header strings again :)
	fgets(str,200,fp);
	fgets(str,200,fp);
	for(i=0;i<k;i++)
	{
		fgets(str,200,fp);// places the pointer in the beginning of necessary string
	}

}
/////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////

//function  reads initial parameters from input.txt according its filling rules
//returns a pointer on first element of initial data array
double* InitialData(char *s2)
{
FILE *fp;
int i=0;
double p=0, inputdata[20];
char s[50], string[100];
double *ptr=inputdata;
      if((fp=fopen(s2,"r"))==NULL)
	       {
		printf("Cannot open file\n");
		exit(1);
	       }
      while(!feof(fp))
  {
	fgets(string,100,fp);//here and below this command is used only for move the pointer on the next line in the file
	fscanf(fp,"%s",s);//read data between two spaces and put it into small string s; pointer is being replaced inside the big string
	p=atof(s);// convert char string to double number
	inputdata[i]=p;
//	printf("%f\n",inputdata[i]);
	fgets(string,100,fp);
	fgets(string,100,fp);
	i++;
  }
		fclose(fp);

	 Z_Position=inputdata[0];
	 OrientationAngle=inputdata[1];
	 FiberNumber=inputdata[2];
	 FiberDiameter=inputdata[3];
	 Pitch=inputdata[4];
	 PhotonYield=inputdata[5];
	 FiberBundling=inputdata[6];
	 QuantumEfficiency=inputdata[7];
	 PMT_Gain=inputdata[8];
	 GainDistribution=inputdata[9];
	 GainSpread=inputdata[10];
	 RandomGen=inputdata[11];
	 Delta_Z=inputdata[12];

return  ptr;


}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//function gives a number of current histogram bin (the fiber, current muon has passed through)
int HistNumber(double x, int flag)
{    //zero-fiber is in the center of beamline crossection; if we have 64 fibers there are: 31 top, 1 zero, 32 bottom
     //coordinates ZERO is in the center of zero-fiber
int number, z;

      if(Pitch>=FiberDiameter)
      {
	 if(x>=0) z=int(floor(fabs(x)))+1;
	 else z=-int((floor(fabs(x))+1));

	 if(z%2==0) number=z/2;
	 else
	 {
		if(z>=0) number=(z-1)/2;
		else number=(z+1)/2;
	 }
       }
       if((Pitch>=FiberDiameter/2)&&(Pitch<FiberDiameter))
       {
	  if(x>=0) z=int(floor(fabs(x)))+1;
	 else z=-int((floor(fabs(x))+1));

	 if(z%2==0) number=z;
	 else
	 {
		if(z>=0) number=(z-1);
		else number=(z+1);
	 }
	 if(flag==1) number+=1;
       }
     return number;
}

//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
double Degrees_to_Rad(double x)//converts degrees to radians
{
   double y=(pi*x)*1.0/180.0;
   return y;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
double* Correction(double *pt)//makes corrections before muon passes through detector's second layer
{
    double R=FiberDiameter/2;
    double x=*pt,y=*(pt+1)/*y coordinate*/;
    *pt=GetX(x,sqrt(3)*R,pt);
    *(pt+1)=GetY(y,sqrt(3)*R,pt);
    return pt;
}
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//this function rotates the fibers on the given angle relatively to initial position
// and puts new coordinates in data array
double* Orientation(double *data)
{
    double x=*data, y=*(data+1),px=*(data+3),py=*(data+4);
    *data=x*cos(Degrees_to_Rad(OrientationAngle))+y*sin(Degrees_to_Rad(OrientationAngle));
    *(data+1)=y*cos(Degrees_to_Rad(OrientationAngle))-x*sin(Degrees_to_Rad(OrientationAngle));
    *(data+3)=px*cos(Degrees_to_Rad(OrientationAngle))+py*sin(Degrees_to_Rad(OrientationAngle));
    *(data+4)=py*cos(Degrees_to_Rad(OrientationAngle))-px*sin(Degrees_to_Rad(OrientationAngle));

    return data;
}
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//if detectors consists of two layers Bundling() function can form the bundle
//to unite some fibers wich will correspond to the same histogram bin
int Bundling(int s)
{
int number=s+FiberNumber/2;
number=int(floor(number/FiberBundling)-ceil(0.5*FiberNumber/FiberBundling));
return number;
}
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////

//function returns absolute meaning of initial Y-coordinate of a particle
//relatively to the fiber this particle is going to pass through
double FractionalPart(double x, double R)
{
double result;
int z=floor(fabs(x));
	if(z%2==0) result=fabs(x)-z;//here we don't forget that we measure all distances from the center of zero-fiber
	else result=R-(fabs(x)-z);
return result;
}
////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//function operates with space and momentum coordinates from G4BeamLine output file
//counts muon pathlenght in fibers, simulates photon generation
//simulates electron production in PMT
double* GetData(double *data, int flag)
{
int ii, histnum/*histogram bin number*/, PhotonNumber,PhotonCounter=0,s;
double x,y,z,px,py,pz,yy;//coordinates of muons
double costeta, sinteta, tgteta, cosfi,sinfi;//angles shown in documentation
double L; //muon pathlength thru fiber
double ElectronNumber;//number of electrons produced by produced in PMT
double Lenght=100.0;// effective fiber length, mm
double dist,sum=0,rand;
double const R=FiberDiameter/2;//fiber radius
double *ptr=data;//pointer on data array
	x=data[0]*1000; y=data[1]*1000;/*from meters to millimeters*/ z=data[2];//space coordinates determine muon position
	px=data[3]; py=data[4]; pz=data[5];// momentum coordinates determine moving directions
	//all the angles are described in Documentations
	costeta=pz/sqrt((py*py+pz*pz));
	data[23]=costeta;
	cosfi=pz/sqrt((px*px+pz*pz));
	data[24]=cosfi;
	sinfi=sqrt(1-cosfi*cosfi);
	data[25]=sinfi;
	sinteta=sqrt(1-costeta*costeta);
	data[26]=sinteta;
	tgteta=sinteta/costeta;
	data[27]=tgteta;
	data[32]=sinfi/cosfi;
		 y=GetY(y,Delta_Z,ptr);
		 x=GetX(x,Delta_Z,ptr);
		 s=HistNumber(GetY(y,R,ptr)/R,flag);
		if(FiberBundling==1)//means no bundling
			histnum=s;
		else
			histnum=Bundling(s);

		   //	if((x>=0)&&(px>=0)||((x<0)&&(px<0)))
		   //		dist=fabs(Lenght/2)-fabs(y);// distance beetween muon enter point and the ending of fiber wire
		   //	if((x>=0)&&(px<0)||((x<0)&&(px>=0)))
		   //		dist=fabs(Lenght/2)+fabs(y);


			if(fabs(y)>R)
				yy=FractionalPart(y/R,R)*R;
			if(fabs(y)<=R)
				yy=fabs(y);
		L=(2*R/cosfi)*sqrt((1-pow((yy-R*tgteta),2.0)/(R*R*(1+tgteta*tgteta))));//main formula, just geometry :)
			    //	if(dist<L*sinfi) L=dist/sinfi;
					data[31]=L;
     if(RandomGen==1)//here and below: random generator is switched on
     PhotonNumber=Poisson(PhotonYield*L);//photon number is distributed with Poisson distribution
     else PhotonNumber=10*L;//random generator is switched off; this option is used only for program testing

//we have to take into account Quantum Efficiency of electron production process
	  for(ii=0;ii<PhotonNumber;ii++)
	  {
		rand=genrand_real1();//generates value distributed with Uniform distribution
		if(rand<=QuantumEfficiency/100) PhotonCounter++;//counts photons which can produce electrons by striking PMT window
	  }

	   PhotonNumber=PhotonCounter;
	   // printf("\t%i\t%f\t%f\t%i\n\n",histnum,y,L,PhotonNumber);

			 if (GainDistribution==0)// if the electrons distributed with the Exponential one
			{
				for(ii=0;ii<PhotonNumber;ii++)
				  {          if(RandomGen==1)
						 ElectronNumber=ExpDist(1.0/PMT_Gain);//electron number is distributed with exponential dist.
						else ElectronNumber=39345;//for testing
						sum+=ElectronNumber;//summing all the electrons
						outputdata[histnum+FiberNumber/2]+=ElectronNumber;// distributes electrons between the bins and stores this data
					      //	printf("%8.6f\t%f  ", ElectronNumber,outputdata[histnum+FiberNumber/2] );
				  }
			}
			 if (GainDistribution==1)//if the electrons distributed with the Gaussian one
			{
				 for(ii=0;ii<PhotonNumber;ii++)
				{
						 ElectronNumber=Normal(PMT_Gain,GainSpread);//electron number is distributed with exponential dist. with mean A
						 sum+=ElectronNumber;
						 outputdata[histnum+FiberNumber/2]+=ElectronNumber;
						// printf("%8.6f  ", ElectronNumber);
				}
			}
			 if((GainDistribution!=1)&&(GainDistribution!=0))
			{printf("GainDistribution is not defined correctly");}
	  data[28]+=sum;
	  data[29]+=sum*y;
	  data[30]+=sum*y*y;//data for counting MEAN and SIGMA
	  return ptr;

}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[]){
char  s1[200], string[200];//s2[50]/*="config.txt"*/,s3[50]/*="output.txt"*/;
double *pt/*data array pointer*/,aa,bb,cc,sigma,mean,meansqr;

FILE *input, *output;
int i=0,ii,counter=0,flag;
double p1;
double data[64];//values 0-22 reserved for data which can be collected from G4BL file
input=fopen(argv[1],"r+");//opens G4BL file
//clrscr();
//printf("Enter the name of configuration file:   ");
//scanf("%s",s2);
//printf("Enter the name of output file:   ");
//scanf("%s",s3);
  if(argc!=4)
  {
	printf("You forgot to type the arguments\n");
	exit(1);
  }
	InitialData(argv[2]);
		for(ii=0;ii<64;ii++)
			data[ii]=0;//cleans the data array from trash

			SearchPosition(input);//searchs the detector in G4BL file and places the pointer in the beginning of the string
				while(!feof(input))
				{
					fscanf(input,"%s",s1);//reads data in char format
					p1=atof(s1);//converts to double

						if((i>=a)&&(i<=b))//stores necessary data in array
						{
							 data[i-a]=p1;

						  //	printf("%f\n",data[i-a]);
						}
						if((i==b)&&(data[2]!=Z_Position))
						{
							break;
						}
					i++;
						if(i==b+1)//operates with data collected at the previous step
						{
							flag=0;
							pt=GetData(Orientation(data),flag);//rotating and simulating
							if((Pitch>=FiberDiameter/2)&&(Pitch<FiberDiameter))
							{
								flag=1;
								pt=GetData(Correction(pt),flag);
							}
							i=0;
							fgets(string, 200,input);//replaces the pointer in the beginning of the next string
						      counter++;
						}
				}
			aa=*(pt+28); bb=*(pt+29); cc=*(pt+30);
			mean=bb/aa;
			meansqr=cc/aa;
			sigma=sqrt(meansqr-pow(mean,2));
			printf("\n   Well done! \n\n mean x=%f\n sigma=%f\n",mean,sigma);
			*(pt+28)=0; *(pt+29)=0; *(pt+30)=0;

		fclose(input);

	 if((output=fopen(argv[3],"w+"))==NULL)
	 printf("Cannot open file.\n");
	 for(int kk=0; kk<FiberNumber; kk++)
	 fprintf(output,"%i %f\n", kk-FiberNumber/2, floor(outputdata[kk]));
	 fprintf(output,"#MEAN = %f\n#SIGMA = %f\n",mean, sigma);

   fclose(output);

//getchar();
return 0;
}