/*********************
 *                        
 *  kinematics.c            
 *                    
 **********************
 *
 * */

#include <stdio.h>
#include <string.h>
#include <signal.h>
#include <errno.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <math.h>

#include <genkin.h>

/*
 #include <itypes.h>
 #include <hdr.h>
 #include <param.h>
 #include <disData.h>
 #include <disIO.h>
 #include <dataIO.h>
 #include <utility.h>
 #include <esr.h>
 #include <omega_system.h>
 #include <particleType.h>
 #include <eventType.h>

 */

#define RESTFRAME -1
#define PARENTFRAME +1

double SQ(double x) {
	double z;

	z = (x) * (x);
	return (z);
}

/**************
 * DotProduct.c
 **************
 */
double DotProduct3(const vector3_t *p1, const vector3_t *p2) {
	return (p1->x * p2->x + p1->y * p2->y + p1->z * p2->z);
}

/****************
 * CrossProduct.c
 ****************
 */
vector3_t CrossProduct3(const vector3_t *p1, const vector3_t *p2) {
	vector3_t c;

	c.x = p1->y * p2->z - p1->z * p2->y;
	c.y = -(p1->x * p2->z - p1->z * p2->x);
	c.z = p1->x * p2->y - p1->y * p2->x;

	return c;
}

/*******************
 *   get_beta()
 *
 *******************/
vector4_t get_beta(vector4_t *boost, int sign) {

	/* find beta 4vector where beta->t = gamma
	 *
	 */
	vector4_t beta;

	beta.space.x = sign * (boost->space.x / boost->t);
	beta.space.y = sign * (boost->space.y / boost->t);
	beta.space.z = sign * (boost->space.z / boost->t);
	/* gamma = E/m */

	beta.t = (boost->t)
			/ sqrt(
					(boost->t) * (boost->t)
							- ((boost->space.x) * (boost->space.x)
									+ (boost->space.y) * (boost->space.y)
									+ (boost->space.z) * (boost->space.z)));
	return beta;
}

/***********
 * lorentz.c
 ***********
 */

vector4_t lorentz(const vector4_t *beta, const vector4_t *pin) {
	vector4_t ret;
	double d, c, c2;
	d = DotProduct3(&(beta->space), &(pin->space));
	c = beta->t / (beta->t + 1.0);
	c2 = c * d + pin->t;
	ret.space.x = pin->space.x + beta->space.x * beta->t * c2;
	ret.space.y = pin->space.y + beta->space.y * beta->t * c2;
	ret.space.z = pin->space.z + beta->space.z * beta->t * c2;
	ret.t = beta->t * (pin->t + d);
	return (ret);
}
/**************
 * CMmomentum.c
 **************/

double CMmomentum(double cm_engy, double m1, double m2) {
	double A, B, C, D, E;

	A = cm_engy * cm_engy;
	B = (m1 + m2) * (m1 + m2);
	C = (m1 - m2) * (m1 - m2);
	D = (A - B) * (A - C);
	E = sqrt(D) / (2.0 * cm_engy);
	return (E);

}
/***********
 * energy.c
 **********/

double energy(double mass, const vector3_t *p) {
	double E, pmagsq;
	pmagsq = SQ(p->x) + SQ(p->y) + SQ(p->z);
	E = sqrt(mass * mass + pmagsq);
	return (E);
}

/********
 * v3mag.c
 ********/

double v3mag(const vector3_t *p) {
	double mag;
	mag = sqrt((p->x) * (p->x) + (p->y) * (p->y) + (p->z) * (p->z));
	return (mag);
}

/*
 ***********************
 *                     *
 *  Sum4vec()   *
 *                     *
 ***********************
 */

vector4_t Sum4vec(vector4_t *vec4, int nvec4) {
	int i;
	vector4_t temp4;

	temp4.space.x = 0;
	temp4.space.y = 0;
	temp4.space.z = 0;
	temp4.t = 0;

	for (i = 0; i < nvec4; i++) {
		temp4.space.x += (vec4 + i)->space.x;
		temp4.space.y += (vec4 + i)->space.y;
		temp4.space.z += (vec4 + i)->space.z;
		temp4.t += (vec4 + i)->t;
	}
	return temp4;
}

/************************
 *
 *  eff_mass()
 *
 ************************/
double eff_mass(vector4_t *v, int nparticles) {
	int i;
	double mass = -1; /* set to -1 for debugging */
	vector4_t vsum;

	/*
	 * initilize vsum to zero
	 */

	vsum.t = 0;
	vsum.space.x = 0;
	vsum.space.y = 0;
	vsum.space.z = 0;

	/*
	 *Sum the nparticles four vectors
	 */
	for (i = 0; i < nparticles; i++) {
		vsum.t += v[i].t;
		vsum.space.x += v[i].space.x;
		vsum.space.y += v[i].space.y;
		vsum.space.z += v[i].space.z;
	}
	/*
	 *Calculate the effective mass
	 */
	mass = sqrt(
			vsum.t * vsum.t
					- (vsum.space.x * vsum.space.x + vsum.space.y * vsum.space.y
							+ vsum.space.z * vsum.space.z));
	return mass;
}

/*
 ***********************
 * lambda3pi()         *
 ***********************
 */
int lambda3pi(vector4_t *vec, int nvec, double *lambda) {
	int i;
	vector4_t parent, beta, vecp[3];
	vector3_t analyzer;

	if (nvec != 3) {
		fprintf(stderr, "Error:: lambda3pi(): you do not have 3 particles!\n");
		*lambda = -1.0;
	} else {
		/* define parent
		 */
		parent = Sum4vec(vec, nvec);

		/*
		 *  Boost to the parent restframe
		 */

		beta = get_beta(&parent, RESTFRAME);
		for (i = 0; i < nvec; i++)
			vecp[i] = lorentz(&beta, &vec[i]);

		/* get the normal to the decay plane
		 */
		analyzer = CrossProduct3(&(vecp->space), &((vecp + 1)->space));

		*lambda =
				DotProduct3(&analyzer, &analyzer)
						/ ((3.0 / 4.0)
								* SQ(SQ(eff_mass(vec, nvec)) / 9.0 - SQ(PIMASS)));
	} /* end of else */
	return 1;
}
/*
 ***********************
 *                     *
 *  helicityAngles()   *
 *                     *
 ***********************
 */

int helicityAngles(vector4_t *vec, int nvec, double *theta, double *phi) {

	int i;
	vector3_t x, y, z, xhel, yhel, zhel, analyzer;
	vector4_t vecp[3], beta, parent;

	/* define parent
	 */
	parent = Sum4vec(vec, nvec);

	/* define lab frame */
	x.x = 1;
	x.y = 0;
	x.z = 0;
	y.x = 0;
	y.y = 1;
	y.z = 0;
	z.x = 0;
	z.y = 0;
	z.z = 1;

	/* define helicity frame */
	zhel = parent.space;
	yhel = CrossProduct3(&z, &(parent.space));
	xhel = CrossProduct3(&yhel, &zhel);/* right handed */

	/* Note that the helicity frame is invariant to the boost
	 * since zhel is along and yhel is normal to the boost.
	 */

	/*
	 *  Boost to the parent restframe
	 */

	beta = get_beta(&parent, RESTFRAME);
	for (i = 0; i < nvec; i++)
		vecp[i] = lorentz(&beta, &vec[i]);

	/*
	 * get the helicity angles of the analyzer
	 *  (the particle that defines the angles)
	 */

	switch (nvec) {/* define the analyzer */
	case 2:
		analyzer = vecp->space; /* 3momentum of the 1st particle */
		break;
	case 3:
		/* use the normal to the decay plane */
		analyzer = CrossProduct3(&(vecp->space), &((vecp + 1)->space));
		break;
	default:
		fprintf(stderr, "Error(helicityAngles()):: to many decay particles\n");
		exit(-1);
	}

	*theta = acos(
			DotProduct3(&analyzer, &zhel) / (v3mag(&analyzer) * v3mag(&zhel)));
	*phi = atan2(DotProduct3(&analyzer, &yhel), DotProduct3(&analyzer, &xhel));

	return 1;

}

/*
 ***********************
 *                     *
 *  END OF FILE        *
 *                     *
 ***********************
 */