///////////////////////////////////////////////////////////////////////////////
// File:	cqw_kin.c
// Summary:	2 Body relativistic constant qw kinematics
//		(w,q) + b -> 2 + x
// Author(s):	Konrad Aniol ( kaniol@calstatela.edu )
//		Juan Carlos Cornejo ( cornejo@jlab.org ) -- rewrite with
//		comments
// Notes:	To compile this file using the GNU C Compiler type
//		gcc -lm -o cqw_kin cqw_kin.c

///////////////////////////////////////////////////////////////////////////////
// Includes
#include<stdio.h>	// Standard C I/O Library
#include<stdlib.h>	// Standard C Library
#include<math.h>	// Standard C Math library

///////////////////////////////////////////////////////////////////////////////
// Definitions for pre-processor ( these are usually constants )
#define _PI 		3.141592654	// The constant PI
#define _ELECTRON_MASS 	0.511e-3 	// In GeV/c^2
#define _PROTON_MASS 	0.938783	// In GeV/c^2
#define _3HE_MASS	2.809413	// In GeV/c^2
#define _4HE_MASS	3.72840		// In GeV/c^2
#define _208PB_MASS	193.728988	// In GeV/c^2
#define _2H_MASS	1.876124	// In GeV/c^2
#define _207TI_MASS	192.798214	// In GeV/c^2

#define _MAX_PX_VALUES	20		// Maximum number of Px values

///////////////////////////////////////////////////////////////////////////////
// Function Pre-headers

/////
/// Converts a value from radians to degrees
/////
double toDegrees( double rad );

/////
/// Converts a value from degrees to radians
/////
double toRadians( double deg );

/////
/// Calculates kinematics based on constant (q,w)
/////
void kin_const_qw( double* values );

/////
/// Begins a new set of kinematic calculations based on constant (q,w)
/////
//void processNewKin( FILE *fout, bool firstPass );
void processNewKin(FILE *fout, short firstPass );

///////////////////////////////////////////////////////////////////////////////
/// MAIN FUNCTION

/////
/// main function
/////
int main( int argc, char** argv )
{
	// Create an array to hold whether the user wants to create a new
	// kinematics or not ( default is no )
	char cont[2];
	cont[0] = 'n';

	// Create an array to hold whether the user wants to store the
	// results in a file or not ( default is no )
	char storeInFile[1];
	storeInFile[2] = 'n';
	FILE *fout;

	// Boolean which holds whether it is the first pass through the
	// kinematics calculations ( default is true )
	short firstPass = 1;

	// Print welcome message and instructions
        printf( "This program calculates the kinematics for a constant qw\n" );
        printf( "The energies are in GeV and momenta in GeV/c\n" );
        printf( "To properly display the results, a terminal window of ");
        printf( "at least 100 characters accross is required.\n\n");

	// Ask if the user wants the results stored in a file so they can be
	// read by the geometry and beam file creation script
	printf( "Would you like to store the results in a file? ( y/n )\n" );
	printf( "This is usefull if you want to create beam and geometry" );
	printf( " files for use with GEANT.\n" );
	scanf( "%s", storeInFile );

	// Check the user's answer
	if ( storeInFile[0] == 'y' ) { // Create the file then
		// Ask for a File Name from user
		printf( "Please enter a filename for the output:\n" );
		char fileName[100];
		scanf( "%s", fileName ); // Read in Filename

		fout = fopen( fileName, "w" );	// Open the file in write mode

		// Check that the file was opened sucessfully
		if ( fout == NULL ) { // File not opened sucessfully
			printf( "Error opening file for writting.  ");
			printf( "Exiting now.\n\n");
			return -1; // Exit with errors.
		}
	}
	
	// Process the kinematics until the user says no ( or just 'n' )
	// or honestly, anything but 'y'
	do {
		cont[0] = 'n';
		processNewKin( fout, firstPass ); // Kinematics calculation
		
		// Ask if the user wants to continue
		printf( "Enter other kinematic settings? ( y/n )\n" );
		scanf( "%s", cont );

		// Update firstPass boolean
		firstPass = 0;
	} while( cont[0] == 'y' );

	return 0;  // Exit sucessfully
}


void processNewKin( FILE *fout, short firstPass )
{
	// Let us first keep track of all the variables we will be using
	
	// Mass variables ( in GeV/c^2 )
	// mb = target mass ( 3,4He or 208Pb for now )
	// ma = incoming mass ( in the case of JLab it's always the electron )
	// m2 = ejected mass ( more than likely a proton )
	// mx0 = Mass of the resulting nuclei ( at ground state )
	// mx = Mass of the resulting nuclei ( with excitation energy )
	double mb,ma,m2,mx0,mx;
	ma = _ELECTRON_MASS; // A standard at JLab

	double x; // this will be a holding value for many things
	double pxvalue[_MAX_PX_VALUES];	// Array to hold the Px values
	double qpass[13]; // Array to store all the results
	
	// Other variables, to be described later... I hope!
	double px,p2,q,w,wmn,wmx,dw,v[4],Q2,xB;
        double thqx,csthqx,thq2,csthq2,th1,thq,sign,eps,thp1,thp2;
        double y,thqmin = 12.5,pa0,ea0;
        double e2,p2_x,p2_x2,t2_x,e2_x,e2_x2,s,sv,csth2x;
        double p20,p2value[20],vpass[13], csthq,csth1,p1;
        char framp2[2],mk_file[2];
        thqmin = toRadians( thqmin );

	int i = 0, j = 0; // Random indices useful for later on
	int ireact;	// Target Z

	// Print welcome message and instructions
	printf( "This program calculates the kinematics for a constant qw\n" );
	printf( "The energies are in GeV and momenta in GeV/c\n" );
	printf( "To properly display the results, a terminal window of ");
	printf( "at least 100 characters accross is required.\n\n");

	// Begin taking input

	// Select target
	printf( "Please select a target:\n3\t - 3He\n4\t - 4He\n208\t - 208Pb\n" );
	scanf( "%d", &ireact ); // Read integer from STDIN
	
	// Setup kinematics for selected target
	switch( ireact ) {
		case 3:	// 3He Selected
			mb = _3HE_MASS;
			m2 = _PROTON_MASS;
			mx0 = _2H_MASS;
			break;
		case 4: // 4He Selected
			mb = _4HE_MASS;
			m2 = _PROTON_MASS;
			mx0 = _3HE_MASS;
			break;
		case 208: // 208Pb Selected
			mb = _208PB_MASS;
			m2 = _PROTON_MASS;
			mx0 = _207TI_MASS;
			break;
		default: // Exit program, because no valid mass selected
			printf( "\nAn incorrect selection was made.  Exiting now.\n" );
			exit( 1 ); // Exit with errors.
	}

	// Read excitation engery from STDIN
	printf( "Enter excitation energy x in GeV\n" );
	scanf( "%lf", &x );
	mx = mx0 + x; // Add excitation energy to ground state mass
	x = (mb*mb - mx*mx )/2./mb; // I honestly don't know why this is so
	// I have to check in on it

	// Read in q and w on the same line
	printf( "For csthqx=-1. the absolute maximum is px= %3lf\n", x );
	printf( "enter q and w ( spaced out )\n" );
	scanf( "%lf %lf", &q, &w ); // Read two integers from STDIN

	// Now read in desired momenta & loop until the user enters a negative
	// number.  Currently we are limited to _MAX_PX_VALUES
	i = 0;	// initialize our index
	px = 0.; // initialize the read-in momentum
	printf( "Enter all the desired momentum (less than %d entries)\n", _MAX_PX_VALUES);
	printf( "Hit enter after each entry. " );
	printf( "To stop entering values simply enter a negative number\n" );
	while( ( px >= 0. ) && ( i < _MAX_PX_VALUES ) ) {
		scanf( "%lf" , &px );
		pxvalue[i] = px;
		i++;
	}
	i--;	// Decrease the count, so that it's still within our array
	
	// Read ea0 ( Energy Beam )
	printf( "Enter the energy beam (ea0 ):\n" );
	scanf( "%lf", &ea0 );

	// Now it's time to calculate the results, and display them

	// Display ea0 & table header
	printf( "\n\nResults:\n" );
	if ( firstPass )
		printf( "Ea0 = %06.4lf GeV\n", ea0 );
	printf("  px      q       w      Q2     thqx      p1     th1     thq");
	printf("     thp1     thp2     eps      p2      xB\n");

	// Check if we want to write the results to the file, if so, write header
	// on the file as well.
	if ( fout != NULL && firstPass ) {
		fprintf( fout, "Ea0 = %06.4lf\n", ea0 );
		fprintf( fout, "  px      q       w      Q2     thqx      p1     th1     thq");
		fprintf( fout, "     thp1     thp2     eps      p2      xB\n");
	}

	// Calculate initial pa0
	pa0 = sqrt( ea0*ea0 - ma*ma );

	// Begin filling storage array
	qpass[0] = pa0;
	qpass[1] = q;
	qpass[2] = w;
	qpass[3] = m2;
	qpass[4] = mb;
	qpass[5] = mx;

	// Now loop through the entered momenta and continue filling the array
	// Remember, i holds the number of entered momenta
	for( j = 0; j < i; j++ ) {
		px = qpass[6] = pxvalue[j]; // initialize with current momentum
		kin_const_qw( qpass ); // Calculate kinematics for constant q,w

		// Retrieve settings from kin_const_qw() stored in an array
		csthq2 = qpass[7];
		csthqx = qpass[8];
		p2 = qpass[9];
		csthq = qpass[10];
		p1 = qpass[11];
		csth1 = qpass[12];

		// Now calculate the new angles
		//printf( "\n\ncsth1:%lf\n\n",csth1 );
		th1 = toDegrees( acos( csth1 ) );
		thqx = toDegrees( acos( csthqx ) );
		thq = toDegrees( acos( csthq ) );
		thq2 = toDegrees( acos( csthq2 ) );
		thp1 = thq + thq2;	// Forward px angle
		thp2 = thq - thq2;	// Backwards px angle

		// Special case for px = 0  The angles are that of q
		if ( px == 0 ) {
			thqx = 0.;
			thp1 = thq;
			thp2 = thq;
		}

		// Calculate Q^2
		Q2 = q*q - w*w;

		x = tan( toRadians( th1/2.) ); // temporary holding value
		eps = 1./( 1. + (2.*q*q*x*x)/Q2 );
		xB = Q2/( 2. * w * m2 );

		// Now output the results
		printf("%.3lf | %.3lf | %.3lf | %.3lf | %05.2lf | %.3lf | %.2lf | %.2lf | %.3lf | %.3lf | %.3lf | %.3lf | %.3lf\n",
		        px,q,w,Q2,thqx,p1,th1,thq,thp1,thp2,eps,p2,xB);
		// Again, if we are to output it to a file, do so again
		if( fout != NULL ) {
			fprintf( fout, "%.3lf | %.3lf | %.3lf | %.3lf | %05.2lf | %.3lf | %.2lf | %.2lf | %.3lf | %.3lf | %.3lf | %.3lf | %.3lf\n",
				px,q,w,Q2,thqx,p1,th1,thq,thp1,thp2,eps,p2,xB);
		}
	

	}

}

///////////////////////////////////////////////////////////////////////////////
double toDegrees( double rad )
{
	return rad * (180./_PI);
}
///////////////////////////////////////////////////////////////////////////////
double toRadians( double deg )
{
	return deg * (_PI/180.);
}

///////////////////////////////////////////////////////////////////////////////
void kin_const_qw( double* values )
{
	// Redefine some local variables again, look at the beginning of the
	// main() function for the definitions of the variables.
	double pa0, q, w, m2, mb, mx, px, csthq2, csthqx, p2, csthq, csth1, p1;
	double e0, ex, e2, th1, x, e1, ea0, Q2, ma;
	ma = _ELECTRON_MASS;  // Once again, ma is the electron mass at JLab

	// Now initialize some variables from the passed values
	pa0 = values[0];
	q = values[1];
	w = values[2];
	m2 = values[3];
	mb = values[4];
	mx = values[5];
	px = values[6];

	// Calculate initial energy
	e0 = w + mb;

	// Calculate the missing energy = sqrt(px^2 + mx^2) in GeV
	ex = sqrt( px*px + mx*mx );

	// Calculate the energy of the ejectile
	e2 = e0 - ex;

	// Calculate the ejectiles momentum
	p2 = sqrt( e2*e2 - m2*m2 );

	// Calculate the angle ( theta ) of the ejectile and missing mass
	csthq2 = ( q*q + p2*p2 - px*px )/( 2.*q*p2 );
	csthqx = ( q*q + px*px - p2*p2 )/( 2.*q*px );
	
	// Store the results so far
	values[7] = csthq2;
	values[8] = csthqx;
	values[9] = p2;

	// Now let us calculate the electron kinematics
	
	// initial energy
	ea0 = sqrt( ma*ma + pa0*pa0 );

	e1 = ea0 - w; // scattered electron energy
	p1 = sqrt( e1*e1 - ma*ma ); // scattered electron momentum

	Q2= q*q - w*w; // Q^2 calculation

	// Electron anagle calculation
	th1 = 2.*asin( sqrt( Q2 / ( 4.*pa0*p1 ) ) );
	csth1 = cos( th1 );

	csthq = (pa0*pa0 + q*q - p1*p1 ) / ( 2.*pa0*q );
	
	// Store results in the array
	values[10] = csthq;
	values[11] = p1;
	values[12] = csth1;
	
	// Ok, and it ends here, just return
	return;
}

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