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

#define PI 3.141592
#define TPCRADIUS 2.466

// Reads file with data in the following order:
//
// number of tracks, vx, vy, vz, magnetic field => esd-summary
// mass, charge, px, py, pz => esd-detail
//
// For now, we are considering classical mechanics, so p = mv (without relativistic considerations)

void main() {

	int i, j, numberOfTracks;
	double vx, vy, vz, w, Rx, Ry, Rz, B, radius, t, dt;
	double gamaX[2], gamaY[2];

	FILE *detail, *summary;

	summary = fopen("esd-summary.dat","r"); // filename must be "esd-summary.dat"
	detail = fopen("esd-detail.dat","r"); // filename must be "esd-detail.dat"

	do { // Reads summary file until number of tracks is different than zero (there are tracks in the event)

		fscanf(summary,"%d %lf %lf %lf %lf\n",&numberOfTracks, &vx,&vy,&vz,&B);

	} while (numberOfTracks == 0);

	vx = 0.01 * vx; // Vertexes are now in meters
	vy = 0.01 * vy;
	vz = 0.01 * vz;
	B = 0.1 * B; // Magnetic field is now in Tesla

	double m[numberOfTracks], q[numberOfTracks], px[numberOfTracks], py[numberOfTracks], pz[numberOfTracks];

	for (i=0;i<numberOfTracks;i++) { // Reads only data regarding whichever event was the first to have a certain number of tracks

		fscanf(detail,"%lf %lf %lf %lf %lf\n",&m[i],&q[i],&px[i],&py[i],&pz[i]);

	}

	printf("Number of tracks: %d \nVx = %em \nVy = %em \nVz = %em \nB = %lf T \n", numberOfTracks, vx, vy, vz, B);

	//Verify conservation of momentum:

	double pxSum=0.0, pySum=0.0, pzSum=0.0;


	for(i=0;i<numberOfTracks;i++) {

		pxSum += px[i];
		pySum += py[i];
		pzSum += pz[i];

	}

	printf("Sum of X momentum: %lf \nSum of Y momentum: %lf \nSum of Z momentum: %lf \n\n",pxSum,pySum,pzSum);


	// The following loops fill the position array for each track and for many different instants of time
	// gamaX is for x and y coordinates of a particle with only px momentum
	// gamaY is for x and y coordinates of a particle with only py momentum
	//
	// gama ( particle position vector ):
	//       X coordinate = x coordinate of gamaX + x coordinate of gamaY + primary x vertex
	//       Y coordinate = y coordinate of gamaX + y coordinate of gamaY + primary y vertex
	//       X coordinate = z coordinate of particle (straight line) + primary z vertex

        FILE *positions;

        positions = fopen("positions.txt","w");

	for(i=0;i<numberOfTracks;i++) { // Loops in each track

	
		// Here we consider the momentum to be m * v
		// If momentum is zero, R is zero!!
		Rx = px[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters
		Ry = py[i] / ( q[i] * B ) * 3.335641; // unit conversion to meters

		w = ( q[i] * B ) / m[i]; // Angular frequency

		j = 0;
		dt = fabs(1/w) / 1000.0; // time step is proportional to (1/w) so it takes plenty of steps

		fprintf(positions,"\n"); // skip a line, for organization

		do { // Loops for every instant of time 

		    t = (double)(j * dt);


		    // "particle" with X momentum only:

		    gamaX[0] = Rx * sin(w*t); 			// 0 = x coordinate
		    gamaX[1] = Rx * ( cos(w*t) - 1 ); 		// 1 = y coordinate


		    // "particle with Y momentum only:

		    gamaY[1] = Ry * sin(w*t); 			// 0 = x coordinate
		    gamaY[0] = -Ry * ( cos(w*t) - 1 ); 		// 1 = y coordinate


		    // Main position vector:
		    double gama[3] = { gamaX[0]+gamaY[0] + vx, gamaX[1]+gamaY[1] + vy, (( pz[i] / m[i] ) * t) + vz }; // vertexes have already been converted from cm to m

		    fprintf(positions, "%e %e %e %e\n", t, gama[0], gama[1], gama[2]); // Print using scientific notation (very small numbers)

		    radius = sqrt(gama[0]*gama[0] + gama[1]*gama[1]); // distance of particle from Z axis ("radial distance")

		    j++;


		} while( radius < TPCRADIUS ); // if the instant of time equals half the period, stop writing! (w = 2*PI / period)

	}

	fclose(positions);


} // Program works for one specific event.

// Remember to organize AliAnalysisTaskMyTask.cxx