#include "ScintillatorArray.hpp"

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <boost/cstdint.hpp>


ScintillatorArray::ScintillatorArray(const std::string &config_dir,
				         const std::string &name,
				         bool no_files)
	: Processor(name, config_dir, no_files, std::vector<int>())
{
	NAME_PARAMETER(n_detectors);
	NAME_PARAMETER(ztarget);
	NAME_PARAMETER(ztarget_offset);
	NAME_PARAMETER(ytarget_offset);
	NAME_PARAMETER(xtarget_offset);
	NAME_PARAMETER(zin);
	NAME_PARAMETER(zout);

	// read the parameter file before setting the input/output channels
	// because they depend on some of the parameters
	init();
	std::vector<std::string> ignore;
//	ignore.push_back("position");
	ignore.push_back("Encartx");
	ignore.push_back("Encarty");
	ignore.push_back("Encartz");
	read_parameters(ignore);

//	NAME_PARAMETER_ARRAY(position, parameter(n_detectors));
	NAME_PARAMETER_ARRAY(Encartx, 10);
	NAME_PARAMETER_ARRAY(Encarty, 10);
	NAME_PARAMETER_ARRAY(Encartz, 10);

	// specify in and output
	NAME_OUTPUT_CHANNEL(e1_multiplicity);
	NAME_OUTPUT_CHANNEL(e2_multiplicity);
	NAME_OUTPUT_CHANNEL(t_multiplicity);

	NAME_INPUT_ARRAY(amplitude1 , parameter(n_detectors));
	NAME_INPUT_ARRAY(amplitude2 , parameter(n_detectors)); 
	NAME_INPUT_ARRAY(tdc_channel, parameter(n_detectors));  
      
	NAME_INPUT_CHANNEL(beta_measured);
	NAME_INPUT_CHANNEL(xin);
	NAME_INPUT_CHANNEL(yin);
	NAME_INPUT_CHANNEL(xtarget);
	NAME_INPUT_CHANNEL(ytarget);
	NAME_INPUT_CHANNEL(xout);
	NAME_INPUT_CHANNEL(yout);
	
	NAME_OUTPUT_ARRAY(energy1    , parameter(n_detectors));
	NAME_OUTPUT_ARRAY(energy2    , parameter(n_detectors));
	NAME_OUTPUT_ARRAY(time       , parameter(n_detectors));
	NAME_OUTPUT_ARRAY(energy1_DC , parameter(n_detectors));
	NAME_OUTPUT_ARRAY(energy2_DC , parameter(n_detectors));
	
	NAME_COEFFICIENTS_ARRAY(cal_energy1 , parameter(n_detectors));
	NAME_COEFFICIENTS_ARRAY(cal_energy2 , parameter(n_detectors));
	NAME_COEFFICIENTS_ARRAY(cal_time    , parameter(n_detectors));
	
	init();
	read_parameters();
	read_coefficients();
}

ScintillatorArray::~ScintillatorArray()
{
}

void ScintillatorArray::process(prespec::viscon::Interface &viscon_interface, int trigger)
{
	if ( input_valid(xtarget) && 
	     input_valid(ytarget) && 
		 input_valid(xin)    && 
		 input_valid(yin)    && 
		 input_valid(xout)    && 
		 input_valid(yout)    && 
		 input_valid(beta_measured) )
	{

		double r_xy = sqrt( (( input_value(xtarget) - input_value(xin))*( input_value(xtarget) - input_value(xin)))+(( input_value(ytarget) - input_value(yin))*( input_value(ytarget) - input_value(yin))) );
		double rx_in = input_value(xtarget) - input_value(xin);
		double ry_in = input_value(ytarget) - input_value(yin);
		double rz_in = input_value(ztarget) - input_value(zin);
		double r_in = sqrt( (rx_in*rx_in)+(ry_in*ry_in)+(rz_in*rz_in) );
		double sin_psi = r_xy / r_in;
		double psi =	asin(sin_psi);
		double cos_Theta = ( input_value(xtarget) - input_value(xin)) / r_xy;
		double sin_Theta = ( input_value(ytarget) - input_value(yin)) / r_xy;



		double dx_ta = (cos(psi)*sin(psi)*cos_Theta) /	parameter(ztarget_offset);
		double dy_ta = (cos(psi)*sin(psi)*sin_Theta) /	parameter(ztarget_offset);

		// the outgoing particle vector
		double x_p = input_value(xout) - ( input_value(xtarget) + dx_ta			    );
		double y_p = input_value(yout) - ( input_value(ytarget) + dy_ta			    );
		double z_p = parameter(zout)   - ( parameter(ztarget)   + parameter(ztarget_offset) );

		double xtarg = input_value(xtarget);
		double xtargoff = parameter(xtarget_offset);
		double ytarg = input_value(ytarget);
		double ytargoff = parameter(ytarget_offset);
		double ztarg = parameter(ztarget);
		double ztargoff = parameter(ztarget_offset);

	
	int Ne1 = input_array_size(energy1);
	set_output(e1_multiplicity, Ne1);
	for (int i = 0; i < Ne1; ++i)
	{
		int index = input_array_index(amplitude1, i);
		double calib_energy1 = calibrate(input_array_value(amplitude1, i), cal_energy1, index);
		double calib_energy1_DC = 0;
		double cos_theta_DC = 0;

	
		// the gamma vector en1
		double x_g_1 = parameter_array(Encartx, index) - ( xtarg + xtargoff );
		double y_g_1 = parameter_array(Encarty, index) - ( ytarg + ytargoff );
		double z_g_1 = parameter_array(Encartz, index) - ( ztarg + ztargoff );



		double ct1 = costheta(x_p, y_p, z_p,      x_g_1, y_g_1, z_g_1);

		correctDoppler(calib_energy1, input_value(beta_measured), ct1, calib_energy1_DC, cos_theta_DC);
		fill_output_array(energy1,    index, calib_energy1);
		fill_output_array(energy1_DC, index, calib_energy1_DC);
	}

	int Ne2 = input_array_size(energy2);
	set_output(e2_multiplicity, Ne2);
	for (int ii = 0; ii < Ne2; ++ii)
	{
		int index2 = input_array_index(amplitude2, ii);
		double calib_energy2 = calibrate(input_array_value(amplitude2, ii), cal_energy2, index2);
		double calib_energy2_DC = 0;
		double cos_theta_DC = 0;

		// the gamma vector en1
		double x_g_2 = parameter_array(Encartx, index2) - ( input_value(xtarget) + parameter(xtarget_offset) );
		double y_g_2 = parameter_array(Encarty, index2) - ( input_value(ytarget) + parameter(ytarget_offset) );
		double z_g_2 = parameter_array(Encartz, index2) - ( parameter(ztarget)   + parameter(ztarget_offset) );

		double ct2 = costheta(x_p, y_p, z_p,      x_g_2, y_g_2, z_g_2);

		correctDoppler(calib_energy2, input_value(beta_measured), ct2, 
					   calib_energy2_DC, cos_theta_DC);
		fill_output_array(energy2,    index2, calib_energy2);
		fill_output_array(energy2_DC, index2, calib_energy2_DC);
	}
	
	
	int Nt = input_array_size(tdc_channel);
	set_output(t_multiplicity, Nt);
	for (int i = 0; i < Nt; ++i)
	{
		int index = input_array_index(tdc_channel, i);
		fill_output_array(time, index, calibrate(input_array_value(tdc_channel, i), cal_time, index));
	}

    }
}

double ScintillatorArray::costheta(double x_p, double y_p,double z_p, 
								   double x_g, double y_g,double z_g)
{
	// calculate the angle between particle trajectory and gamma trajectory
	double rg_rp =  x_g*x_p + y_g*y_p + z_g*z_p ;
	double rg = sqrt(x_g*x_g + y_g*y_g + z_g*z_g);
	double rp = sqrt(x_p*x_p + y_p*y_p + z_p*z_p);
	return rg_rp / rg / rp;
}

void ScintillatorArray::correctDoppler(double E, double beta, double ct,
				   					   double &Erest, double &ct_rest)
{
	// do the Doppler correction
	double bct   = beta*ct;
	double gamma = 1./sqrt(1-beta*beta);
	Erest        = E * gamma * (1-bct);
	ct_rest      = (ct - beta)/(1-bct);		
}
