/*
This script is calculating the New Voltages needed for BB scintillators base on the Gain Curves measured in the TestLab.
The procedure as following:

1. Calculate the desired voltages for each scintillator based on the Gain Curves.
2. Take a measurement with the calculated voltages.
3. find the punch through points / elastic protons.
4. Run this script with updated information for the included scripts.
5. Apply the new voltages, and check the ADC channel.

6. Make additional iterations if needed


read_hv_mapping.C - read the mapping DB.
read_voltages.C   - read the voltages saved by HV GUI.
read_gains.C      - read the gain curves parameters from DB file
read_adc.C      - read the punch through points for E and dE planes

*/

#include <iostream>
#include "read_hv_mapping.C"
#include "read_voltages.C"
#include "read_gains.C"
#include "read_adc.C"



void main()	{


	
	//=================  read_hv_mapping.C  ======================
	int EL_crate[24]={0}, EL_slot[24]={0}, EL_ch[24]={0};
	int ER_crate[24]={0}, ER_slot[24]={0}, ER_ch[24]={0};
	int dEL_crate[24]={0}, dEL_slot[24]={0}, dEL_ch[24]={0};
	int dER_crate[24]={0}, dER_slot[24]={0}, dER_ch[24]={0};

	//================== read_voltages.C    ======================
	const int V_rows=16,V_cols=12;				// Number of rows and col (number of channels in card). in the *.set file
	float Voltages[V_rows][V_cols]={0};			// Matrix of voltages as it appears in *.set file

	int crate[20]={0},sl[20]={0};				// crate - is a HV crate. This info is read from the *.set file


	//=================== read_gains.C      ======================

	double c_EL[24]={0},c_ER[24]={0},c_dEL[24]={0},c_dER[24]={0};		//  K*V**n+const  :   K = c_ , n = p_ , const = cnst_
	double cnst_EL[24]={0},cnst_ER[24]={0},cnst_dEL[24]={0},cnst_dER[24]={0};
	int p_EL[24],p_ER[24],p_dEL[24],p_dER[24];

	//=================== read_adc.C      ======================

	float EL[24]={0}, ER[24]={0}, dEL[24]={0}, dER[24]={0};		// ADC values.

	//=====================================================


read_hv_mapping(&EL_crate[], &EL_slot[], &EL_ch[], &ER_crate[], &ER_slot[], &ER_ch[], &dEL_crate[], &dEL_slot[], &dEL_ch[], &dER_crate[], &dER_slot[], &dER_ch[]);
read_voltages( &Voltages[][],V_rows,V_cols,&crate[],&sl[]);
read_gains(&c_EL[],&c_ER[],&c_dEL[],&c_dER[],&cnst_EL[],&cnst_ER[],&cnst_dEL[],&cnst_dER[],&p_EL[],&p_ER[],&p_dEL[],&p_dER[]);
read_adc( &EL[], &ER[], &dEL[], &dER[]);

//===============================================================================================

int i,j,k;
float EL_V_new[24]={0},ER_V_new[24]={0},dEL_V_new[24]={0},dER_V_new[24]={0};
int	el_crate,el_slt,el_ch,
	er_crate,er_slt,er_ch,
	del_crate,del_slt,del_ch,
	der_crate,der_slt,der_ch;

float	V_el, ADC_el_expected, ADC_el_actual, D_el, V_el_th, D_ADC_el, D_V_el,
	V_er, ADC_er_expected, ADC_er_actual, D_er, V_er_th, D_ADC_er, D_V_er,
	V_del, ADC_del_expected, ADC_del_actual, D_del, V_del_th, D_ADC_del, D_V_del,
	V_der, ADC_der_expected, ADC_der_actual, D_der, V_der_th, D_ADC_der, D_V_der;

float slope;

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

	int el_crate=EL_crate[i];
	int el_slt=EL_slot[i];
	int el_ch=EL_ch[i];

	if(crate[EL_slot[i]] == el_crate)	//check if the mapping match the actual crate
		{
		V_el=Voltages[el_slt][el_ch];
		}
	else cout<<"PROBLEM 1"<<endl;

	ADC_el_expected	=c_EL[i]*V_el**p_EL[i]+cnst_EL[i];
	ADC_el_actual	=EL[i];
	D_V_el=(ADC_el_actual-ADC_el_expected)/(c_EL[i]*p_EL[i]*V_el**(p_EL[i]-1));

	EL_V_new[i]=V_el+D_V_el;



	//==== For E right	

	er_crate=ER_crate[i];
	er_slt=ER_slot[i];
	er_ch=ER_ch[i];

	if(crate[er_slt]==er_crate)	//check if the mapping match the actual crate
		{
		V_er=Voltages[er_slt][er_ch];
		}
	else cout<<"PROBLEM 2"<<endl;

	ADC_er_expected	=c_ER[i]*V_er**p_ER[i]+cnst_ER[i];
	ADC_er_actual	=ER[i];
	D_V_er=(ADC_er_actual-ADC_er_expected)/(c_ER[i]*p_ER[i]*V_er**(p_ER[i]-1));

	ER_V_new[i]=V_er+D_V_er;




	//===== For dE left

	del_crate=dEL_crate[i];
	del_slt=dEL_slot[i];
	del_ch=dEL_ch[i];

	if(crate[del_slt]==del_crate)	//check if the mapping match the actual crate
		{
		V_del=Voltages[del_slt][del_ch];
		}
	else cout<<"PROBLEM 3"<<endl;

	ADC_del_expected	=c_dEL[i]*V_del**p_dEL[i]+cnst_dEL[i];
	ADC_del_actual	=dEL[i];
	D_V_del=(ADC_del_actual-ADC_del_expected)/(c_dEL[i]*p_dEL[i]*V_el**(p_dEL[i]-1));
	
	dEL_V_new[i]=V_del+D_V_del;




	//===== For dE right

	der_crate=dER_crate[i];
	der_slt=dER_slot[i];
	der_ch=dER_ch[i];

	if(crate[der_slt]==der_crate)	//check if the mapping match the actual crate
		{
		V_der=Voltages[der_slt][der_ch];
		}
	else cout<<"PROBLEM 4"<<endl;

	ADC_der_expected	=c_dER[i]*V_der**p_dER[i]+cnst_dER[i];
	ADC_der_actual	=ER[i];

	D_V_der=(ADC_der_actual-ADC_der_expected)/(c_dER[i]*p_dER[i]*V_der**(p_dER[i]-1));	

	dER_V_new[i]=V_der+D_V_der;
	}

//------------------------------------ WRITE TO FILE------------------------------------------------------
int s1[2]={0}, s2[2]={0}, s3[2]={0}, s4[2]={0};

//==== this for is needed only to define the slot number to the scintillator type: left, right, E or dE.
for(i=0;i<24;i++)
	{
	if(s1[0]!=EL_slot[i])
		{
		s1[1]=EL_slot[i];
		}
	else s1[0]=EL_slot[i];

	if(s2[0]!=ER_slot[i])
		{
		s2[1]=ER_slot[i];
		}
	else s2[0]=ER_slot[i];

	if(s3[0]!=dEL_slot[i])
		{
		s3[1]=dEL_slot[i];
		}
	else s3[0]=dEL_slot[i];


	if(s4[0]!=dER_slot[i])
		{
		s4[1]=dER_slot[i];
		}
	else s4[0]=dER_slot[i];
	}

ofstream file;
char *line,*val;
line=(char*)malloc(1000);
val= (char*)malloc(1000);
remove("new_V.dat");
k=0;
for(i=0;i<V_rows;i++)
	{
	file.open("new_V.dat", ios::app);
	sprintf(line,"bbps%d:2002 S%d DV ",crate[i],sl[i]);
	file<<line;
	file.close();
	file.open("new_V.dat", ios::app);
	
	for(j=0;j<V_cols;j++)
		{
		if(i==s1[0]||i==s1[1])	{	sprintf(val,"%.1f ",EL_V_new[k]); }//cout<<" line 1"<<endl;}
		if(i==s2[0]||i==s2[1])	{	sprintf(val,"%.1f ",ER_V_new[k]); }//cout<<" line 2"<<endl;}
		if(i==s3[0]||i==s3[1])  {	sprintf(val,"%.1f ",dEL_V_new[k]);} //cout<<" line 3"<<endl;}
		if(i==s4[0]||i==s4[1])  {	sprintf(val,"%.1f ",dER_V_new[k]);} //cout<<" line 4"<<endl;}
		if(i!=s1[0] && i!=s1[1] && i!=s2[0] && i!=s2[1] && i!=s3[0] && i!=s3[1] && i!=s4[0] && i!=s4[1]) sprintf(val,"%.1f ",Voltages[i][j]);
		file<<val;
		}
		file<<endl;
		file.close();
	if(k<23) k++;
	else k=0;
	}
		
}