/*
 *  MCP_sim.cpp
 *  
 *
 *  Created by David Yu on 7/30/07.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */	

#include <cstdlib>
#include <math.h>
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include "quartic.cpp"
using namespace std;

#define TUBE_RADIUS 1.0e-5
#define ELECTRIC_FIELD_X (2500000*sin(CHEVRON_ANGLE*2*PI/360))
#define ELECTRIC_FIELD_Z (-2500000*cos(CHEVRON_ANGLE*2*PI/360))
#define TUBE_LENGTH .0008
#define ENERGY_MAX (400*1.602e-19)
#define MEAN_SECONDARY_ENERGY (8*1.602e-19)
#define ELECTRON_CHARGE 1.602e-19
#define ELECTRON_MASS 9.11e-31
#define N 8
#define ALPHA 0.62
#define BETA 0.65
#define DELTA_MAX 2.1
#define PI 3.1415926
#define CHEVRON_ANGLE 8

struct Electron 
{
	double x1, y1, z1, vx1, vy1, vz1, t1;
	double x2, y2, z2, vx2, vy2, vz2, t2;
	int numberOfSecondaries;
};

int numberOfSecondaries ( Electron electronInitial);
void trajectory ( Electron & trajectoryElectron );
int factorial (int n);

//trajectory: argument is an electron with initial variables defined, final variables undefined. trajectory() fills the final variables. Also, fills the number of secondaries. 
void trajectory ( Electron & trajectoryElectron )
{
	//cout << "Starting trajectory calculations...\n";
	double x1 = trajectoryElectron.x1;
	double y1 = trajectoryElectron.y1;
	double z1 = trajectoryElectron.z1;
	double vx1 = trajectoryElectron.vx1;
	double vy1 = trajectoryElectron.vy1;
	double vz1 = trajectoryElectron.vz1;
	double tInitial = trajectoryElectron.t1;
	//cout << "Initial data (x1,y1,z1,vx1,vy1,vz1) = (" << x1 << "," << y1 << "," << z1 << "," << vx1 << "," << vy1 << "," << vz1 << ")\n";
	
	//cout << "Solving quartic...\n";
	double a = -4*ELECTRON_MASS/(ELECTRON_CHARGE*ELECTRIC_FIELD_X)*vx1;
	double b = 4*pow(ELECTRON_MASS/(ELECTRON_CHARGE*ELECTRIC_FIELD_X),2)*(-1*ELECTRON_CHARGE*ELECTRIC_FIELD_X/ELECTRON_MASS*x1+pow(vx1,2)+pow(vy1,2));
	double c = 8*pow(ELECTRON_MASS/(ELECTRON_CHARGE*ELECTRIC_FIELD_X),2)*(vx1*x1+vy1*y1);
	double d = (pow(x1,2)+pow(y1,2)-pow(TUBE_RADIUS,2)) * 4.0*pow(ELECTRON_MASS/(ELECTRON_CHARGE*ELECTRIC_FIELD_X),2);
	//cout << "Quartic = x^4 + " << a << "x^3 + " << b << "x^2 + " << c << "x + " << d << "\n";
	double t;
	//cout << "d=" << d << " and distance = " << sqrt(pow(x1,2)+pow(y1,2))-TUBE_RADIUS << "\n";
	if ( abs(sqrt(pow(x1,2)+pow(y1,2))-TUBE_RADIUS) <= 1e-8 )    //This is annoying. Since the equations are solved numerically, d is never actually 0; a cutoff is needed to determine if the quartic actually reduces to a cubic.
	{
		//cout << "Reduces to cubic...\n";
		double discriminant = 4*pow(a,3)*c-pow(a,2)*pow(b,2)+4*pow(b,3)-18*a*b*c+27*pow(c,2);
		//cout << "Discriminant is " << discriminant << "\n";
		if (discriminant > 0)
		{
			t = (cbrt(4.5*a*b-pow(a,3)-13.5*c+1.5*pow(3*discriminant,0.5))+cbrt(4.5*a*b-pow(a,3)-13.5*c-1.5*sqrt(3*discriminant))-a)/3;
		} else if (discriminant < 0)
		{
			double r = sqrt( pow((4.5*a*b-pow(a,3)-13.5*c),2)-6.75*discriminant );
			double theta = atan2(1.5*sqrt(-3*discriminant) , (4.5*a*b-pow(a,3)-13.5*c));
			double t0 = (-a + 2*cbrt(r)*cos(theta/3))/3;
			double t1 = (-a + 2*cbrt(r)*cos(theta/3-2*PI/3))/3;
			double t2 = (-a + 2*cbrt(r)*cos(theta/3+2*PI/3))/3;
			
			//t = least positive root
			if ( (t0 > 0) && (t0 <= t1 || t1 < 0) && (t0 <= t2 || t2 < 0) ) { t = t0; }
			else if ( (t1 > 0) && (t1 <= t0 || t0 < 0) && (t1 <= t2 || t2 < 0) ) { t = t1; }
			else if ( (t2 > 0) && (t2 <= t0 || t0 < 0) && (t2 <= t1 || t1 < 0) ) { t = t2; }
			else { cout << "86: Your logic is fuzzy!"; exit(1); }
		} else 
		{
			cout << "Discriminant = 0!\n";
			exit(1);
		}
	} else 
	{
		complex<double> t0c, t1c, t2c, t3c;
		quartic(a,b,c,d,t0c,t1c,t2c,t3c);
		double t0, t1, t2, t3;
		if ( abs(imag(t0c))<1e-13 ) { t0 = real(t0c); } else { t0 = -1; }
		if ( abs(imag(t1c))<1e-13 ) { t1 = real(t1c); } else { t1 = -1; }
		if ( abs(imag(t2c))<1e-13 ) { t2 = real(t2c); } else { t2 = -1; }
		if ( abs(imag(t3c))<1e-13 ) { t3 = real(t3c); } else { t3 = -1; }
		//cout << "Complex solutions are " << t0c << ", " << t1c << ", " << t2c << ", " << t3c << "\n";
		//cout << "Real solutions are " << t0 << ", " << t1 << ", " << t2 << ", " << t3 << "\n";
		
		if ( (t0>0) && (t0<=t1 || t1<0) && (t0<=t2 || t2<0) && (t0<=t3 || t3<0) ) { t=t0; }
		else if ( (t1>0) && (t1<=t0 || t0<0) && (t1<=t2 || t2<0) && (t1<=t3 || t3<0) ) { t=t1; }
		else if ( (t2>0) && (t2<=t0 || t0<0) && (t2<=t1 || t1<0) && (t2<=t3 || t3<0) ) { t=t2; }
		else if ( (t3>0) && (t3<=t0 || t0<0) && (t3<=t1 || t1<0) && (t3<=t2 || t2<0) ) { t=t3; }
		else { cout << "107: Your logic is fuzzy!\n"; exit(1); }
	}
		
	
	
	//cout << "Time of trajectory = " << t << "\n";
	trajectoryElectron.x2 = x1 + vx1*t + -0.5*ELECTRON_CHARGE/ELECTRON_MASS*ELECTRIC_FIELD_X*pow(t,2);
	trajectoryElectron.y2 = y1 + vy1*t;
	trajectoryElectron.z2 = z1 + vz1*t + -0.5*ELECTRON_CHARGE/ELECTRON_MASS*ELECTRIC_FIELD_Z*pow(t,2);
	trajectoryElectron.vx2 = vx1 + -ELECTRON_CHARGE/ELECTRON_MASS*ELECTRIC_FIELD_X*t;
	trajectoryElectron.vy2 = vy1;
	trajectoryElectron.vz2 = vz1 + -ELECTRON_CHARGE/ELECTRON_MASS*ELECTRIC_FIELD_Z*t;
	trajectoryElectron.t2 = tInitial + t;
	trajectoryElectron.numberOfSecondaries = numberOfSecondaries(trajectoryElectron);
	if ( trajectoryElectron.z2 < 0 )
	{
		trajectoryElectron.numberOfSecondaries = 0;
	}
	//cout << "Final data (x2,y2,z2,vx2,vy2,vz2) = (" << trajectoryElectron.x2 << "," << trajectoryElectron.y2 << "," << trajectoryElectron.z2 << "," << trajectoryElectron.vx2 << "," << trajectoryElectron.vy2 << "," << trajectoryElectron.vz2 << ")\n";
	//cout << "Number of secondaries = " << trajectoryElectron.numberOfSecondaries << "\n";
	//cout << "**********Check: x^2+y^2 = " << pow(trajectoryElectron.x2,2)+pow(trajectoryElectron.y2,2) << "\n";
	
}

/**************************************************************************************************/
/**************************************************************************************************/

//The bend in the chevron
void reverseTrajectory ( Electron & trajectoryElectron )
{
	double a = -0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS;
	double b = trajectoryElectron.vz1;
	double c = trajectoryElectron.z1 - TUBE_LENGTH/2;
	double t = (-b+sqrt(pow(b,2)-4*a*c))/(2*a);
	//cout << "t = " << t << "\n";
	trajectoryElectron.x1 = -1*(-0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_X/ELECTRON_MASS*pow(t,2)+trajectoryElectron.vx1*t+trajectoryElectron.x1);
	trajectoryElectron.y1 = trajectoryElectron.vy1*t+trajectoryElectron.y1;
	trajectoryElectron.z1 = -0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS*pow(t,2)+trajectoryElectron.vz1*t+trajectoryElectron.z1;
	trajectoryElectron.vx1 = -1*(-ELECTRON_CHARGE*ELECTRIC_FIELD_X/ELECTRON_MASS*t + trajectoryElectron.vx1);
	trajectoryElectron.vz1 = -ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS*t + trajectoryElectron.vz1;
	trajectoryElectron.t1 = trajectoryElectron.t1 + t;
	trajectory(trajectoryElectron);
}
		

/**************************************************************************************************/
/**************************************************************************************************/

//The end of the tube
void endOfTube ( Electron & trajectoryElectron)
{
	double a = -0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS;
	double b = trajectoryElectron.vz1;
	double c = trajectoryElectron.z1 - TUBE_LENGTH;
	double t = (-b+sqrt(pow(b,2)-4*a*c))/(2*a);
	trajectoryElectron.x2 = -1*(-0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_X/ELECTRON_MASS*pow(t,2)+trajectoryElectron.vx1*t+trajectoryElectron.x1);
	trajectoryElectron.y2 = trajectoryElectron.vy1*t+trajectoryElectron.y1;
	trajectoryElectron.z2 = -0.5*ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS*pow(t,2)+trajectoryElectron.vz1*t+trajectoryElectron.z1;
	trajectoryElectron.vx2 = -1*(-ELECTRON_CHARGE*ELECTRIC_FIELD_X/ELECTRON_MASS*t + trajectoryElectron.vx1);
	trajectoryElectron.vz2 = -ELECTRON_CHARGE*ELECTRIC_FIELD_Z/ELECTRON_MASS*t + trajectoryElectron.vz1;
	trajectoryElectron.t2 = trajectoryElectron.t1 + t;
	trajectoryElectron.numberOfSecondaries = 0;
}


/**************************************************************************************************/
/**************************************************************************************************/


//takes an initial electron, returns the number of electrons yielded in the collision
int numberOfSecondaries ( Electron electronPrimary )
{
	//cout << "Computing number of secondaries...\n";
	double x1 = electronPrimary.x2;
	double y1 = electronPrimary.y2;
	double z1 = electronPrimary.z2;
	double vx1 = electronPrimary.vx2;
	double vy1 = electronPrimary.vy2;
	double vz1 = electronPrimary.vz2;
	double energyPrimary = ELECTRON_MASS/2*(pow(vx1,2)+pow(vy1,2)+pow(vz1,2));
	double cosThetaPrimary = -(vx1*x1+vy1*y1)/(sqrt(pow(x1,2)+pow(y1,2))*sqrt(pow(vx1,2)+pow(vy1,2)+pow(vz1,2)));

	double deltaMean = DELTA_MAX * pow(energyPrimary/ENERGY_MAX*sqrt(fabs(cosThetaPrimary)),BETA)*exp(ALPHA*(1-cosThetaPrimary)+BETA*(1-energyPrimary/DELTA_MAX*sqrt(fabs(cosThetaPrimary))));
	//cout << "***energyPrimary = " << energyPrimary/ELECTRON_CHARGE << "eV, deltaMean = " << deltaMean << "\n";
	//Poisson distribution: given mean deltaMean, feed random number into Poisson distribution, yielding eta electrons.
	double randomnumber = rand();
	for (int i=1; i<10; i++)
	{
		randomnumber = rand();
	}//Get well into the random number chain
	randomnumber = randomnumber/(RAND_MAX + 1.0);
	int eta = 0;
	double integral = 0;
	while ( integral < randomnumber )
	{
		integral += pow(deltaMean,eta)*exp(-1*deltaMean)/factorial(eta);
		eta += 1;
	}
	//cout << eta-1 << " secondary electrons produced\n";
	return eta-1;
}

/**************************************************************************************************/
/**************************************************************************************************/

int factorial (int n)
{
	if (n == 0)
	{
		return 1;
	} else
	{
		return n*factorial(n-1);
	}
}

/**************************************************************************************************/
/**************************************************************************************************/

int main (int argc, char *argv[])
{
	//Variables
	Electron currentElectron;
	int m = 0; //index of nextElectronArray
	double energy, velocity, theta, phi, vxprime, vyprime, vzprime; //variables used in computing the properties of secondary electrons
	int nextStageNumber=1, lastStageNumber = 1; //number of electrons in the next stage
	double randomnumber1, randomnumber2, randomnumber3;


	//initialize random number generator
	srand((unsigned)time(0));

	//read in command line arguments
	if (argc != 8) { cout << "Wrong number of arguments\n"; return(0); }
	
	Electron firstElectron;
	firstElectron.x1 = strtod(argv[1],NULL);
	firstElectron.y1 = strtod(argv[2],NULL);
	firstElectron.z1 = strtod(argv[3],NULL);
	firstElectron.vx1 = strtod(argv[4],NULL);
	firstElectron.vy1 = strtod(argv[5],NULL);
	firstElectron.vz1 = strtod(argv[6],NULL);
	firstElectron.t1 = strtod(argv[7],NULL);
	trajectory(firstElectron); //fill in the final coordinates of the first electron
	//cout << "firstElectron: vx2=" << firstElectron.vx2 << ", vy2=" << firstElectron.vy2 << "\n";
	//electronArray will contain the current stage of electrons
	//cout << "First collision: " << firstElectron.x2 << ", " << firstElectron.y2 << ", " << firstElectron.z2 << ", " << firstElectron.t2 << "\n";

	Electron * electronArray;
	electronArray = new Electron [1];
	electronArray[0] = firstElectron;
	
	for (int i=1; i<=N; i++)
	{
		//calculate total number of secondaries in next stage. create a temporary array of this size.
		lastStageNumber = nextStageNumber;
		nextStageNumber = 0;
		//cout << "*****************************************************************************************************************\n";
		cout << "On stage " << i << "\n";
		//cout << "*****************************************************************************************************************\n";

		for (int j=0; j < lastStageNumber ; j++)
		{
			nextStageNumber += electronArray[j].numberOfSecondaries;
		}
		cout << "nextStageNumber = " << nextStageNumber << " and lastStageNumber = " << lastStageNumber << "\n";
		Electron * nextElectronArray;
		nextElectronArray = new Electron [nextStageNumber];
		
		//Fill nextElectronArray
		m=0;
		for (int k=0; k<lastStageNumber; k++)
		{
			currentElectron = electronArray[k];
			for (int l=1; l<=currentElectron.numberOfSecondaries; l++)
			{
				randomnumber1 = rand()/(RAND_MAX+1.0);
				randomnumber2 = rand()/(RAND_MAX+1.0);
				randomnumber3 = rand()/(RAND_MAX+1.0);
				
				energy = pow(-2*pow(MEAN_SECONDARY_ENERGY,2)*log(1-randomnumber1),0.5);
				velocity = pow(2*energy/ELECTRON_MASS,0.5);
				theta = asin(2*randomnumber2-1);
				phi = PI*(randomnumber3-0.5);
				//cout << "Energy = " << energy/ELECTRON_CHARGE << "eV, velocity = " << velocity << ", theta = " << theta << ", phi = " << phi << "\n";
				
				nextElectronArray[m].x1 = currentElectron.x2;
				nextElectronArray[m].y1 = currentElectron.y2;
				nextElectronArray[m].z1 = currentElectron.z2;
				nextElectronArray[m].t1 = currentElectron.t2;
				
				vxprime = velocity*cos(theta);
				vyprime = velocity*sin(theta)*cos(phi);
				vzprime = velocity*sin(theta)*sin(phi);
				nextElectronArray[m].vx1 = -1*vxprime*(currentElectron.x2/sqrt(pow(currentElectron.x2,2)+pow(currentElectron.y2,2)))-vyprime*(currentElectron.y2/sqrt(pow(currentElectron.x2,2)+pow(currentElectron.y2,2)));
				nextElectronArray[m].vy1 = -1*vxprime*(currentElectron.y2/sqrt(pow(currentElectron.x2,2)+pow(currentElectron.y2,2))) + vyprime*(currentElectron.x2/sqrt(pow(currentElectron.x2,2)+pow(currentElectron.y2,2)));
				nextElectronArray[m].vz1 = vzprime;
				
				trajectory(nextElectronArray[m]);
				if (nextElectronArray[m].z1 <= TUBE_LENGTH/2 && nextElectronArray[m].z2 >= TUBE_LENGTH/2)
				{
					reverseTrajectory(nextElectronArray[m]);
				}
				if (nextElectronArray[m].z2 >= TUBE_LENGTH)
				{
					cout << "Electron left tube\n";
					endOfTube(nextElectronArray[m]);
				}
				//cout << "Generated electron, stage " << i << "number " << m << " with data: \n x1 = " << nextElectronArray[m].x1 << "\n y1 = " << nextElectronArray[m].y1 << "\n z1 = " << nextElectronArray[m].z1 << "\n vx1 = " << nextElectronArray[m].vx1 << "\n vy1 = " << nextElectronArray[m].vy1 << "\n vz1 = " << nextElectronArray[m].vz1 << "\n t1 = " << nextElectronArray[m].t1 << "\n x2 = " << nextElectronArray[m].x2 << "\n y2 = " << nextElectronArray[m].y2 << "\n z2 = " << nextElectronArray[m].z2 << "\n vx2 = " << nextElectronArray[m].vx2 << "\n vy2 = " << nextElectronArray[m].vy2 << "\n vz2 = " << nextElectronArray[m].vz2 << "\n t2 = " << nextElectronArray[m].t2 << "\n numberOfSecondaries = " << nextElectronArray[m].numberOfSecondaries << "\n";
				m++;
			}
		}
		//delete [] electronArray;
		//Electron * electronArray;
		electronArray = new Electron [nextStageNumber];
		for (int n=0; n<nextStageNumber; n++)
		{
			electronArray[n] = nextElectronArray[n];
		}
		//delete [] nextElectronArray;
	}
	
	//smear electronArray by Gaussians
	double histo [10000];
	for (int i=0; i<10000; i++)
	{
		histo[i]=0;
	}
	double time_avg;
	double gain;
	double tt;
	double totalgain = 0;
	double sigma_tt;
	for (int i=0; i<nextStageNumber; i++)
	{
		if (electronArray[i].z2 < TUBE_LENGTH)
		{
			double vz_avg=3281349.20821114;
			double sigma_vz_avg=977034.491248129;
			double zz = (electronArray[i].z1+electronArray[i].z2)/2.0;
			gain = exp(17269.4*(TUBE_LENGTH-zz));//fitted to gain curve
			totalgain += gain;
			sigma_tt = sigma_vz_avg*(TUBE_LENGTH-zz)/pow(vz_avg,2);//=sigma_vz*length/vz^2
			tt = (TUBE_LENGTH-zz)/vz_avg;
			//cout << "sigma_tt = " << sigma_tt << ", tt = " << tt << "\n";
			for (int j=0; j<10000; j++)
			{
				histo[j] += gain/(sigma_tt*sqrt(2*PI))*exp(-pow((j*4e-14-tt),2)/(2*pow(sigma_tt,2)))*4e-14;
				//cout << (int)1.0/(tts*sqrt(2*PI))*exp(-1*pow((j*4e-12-time_avg),2)/(2*pow(tts,2))) << "\n";
			}
		}
	}
	cout << "total gain is " << totalgain << "\n";
	for (int i=0; i<10000; i++)
	{
		//cout << histo[i] << "\n";
	}

	//Files to write out to
	ofstream myfile;
	string file = "single_event.txt";
	const char *filename = file.c_str();
	myfile.open(filename,ios::out);
	for (int k=0; k<10000; k++)
	{
		myfile << k*1e-13 << " " << histo[k] << "\n";
	}
	for (int k=0; k<nextStageNumber; k++)
	{
		myfile << electronArray[k].z2/electronArray[k].t2 << "\n";
	}
	myfile.close();
	
/*	
	for (int i=1; i<sizeof(electronArray)/sizeof(*electronArray); i++)
	{
		cout << "(" << electronArray[i].x1 << ", " << electronArray[i].y1 << ", " << electronArray[i].z1 << ", " << electronArray[i].vx1 << ", " << electronArray[i].vy1 << ", " << electronArray[i].vz1 << ")\n";
	}
*/
/*	file = "single_event_full_remainder.txt";
	const char *filename2 = file.c_str();
	myfile.open(filename2,ios::out);
	for (int i=1; i<nextStageNumber; i++)
	{
		myfile << electronArray[i].x2 << " " << electronArray[i].y2 << " " << electronArray[i].z2 << " " << electronArray[i].vx2 << " " << electronArray[i].vy2 << " " << electronArray[i].vz2 << " " << electronArray[i].t2 << "\n";
	}*/

}