/*

This script calculate the desired voltages for HAND scintillators in two cases:

1. Cosmic run=1;
2. Beam (run=0)

The difference is the desired ADC channel at each layer. For cosmic, the energy deposit is different between layer because of different heights of bars. The ADC channel location is added based on the minimum ionizing particlas with 2.2 MeV/cm Energy deposit.

In the beam case, all bars have same width, however the kinetic energy of protons is decreasing as we moving from layer 1 to 6. Correct values for the ADC channels must be added to the code.


This script must be run with "read_gain.C"
and gain curve parameter must be available in the corresponding .dat file.


The output of this script are two files, one for left side and second for the right side, with voltage values.


*/






#include "read_gains.C"
#include <iostream>

void desired_voltage()	{

const int nbars[7]={32,30,24,22,12,12,12};



//======================== read_gains.C ====================
double  c_L1[30]={0},c_R1[30]={0},
	c_L2[24]={0},c_R2[24]={0},
	c_L3[22]={0},c_R3[22]={0},
	c_L4[12]={0},c_R4[12]={0},
	c_L5[12]={0},c_R5[12]={0},
	c_L6[12]={0},c_R6[12]={0};		//  K*V**n+const  :   K = c_ , n = p_ , const = cnst_

double  cnst_L1[30]={0},cnst_R1[30]={0},
	cnst_L2[24]={0},cnst_R2[24]={0},
	cnst_L3[22]={0},cnst_R3[22]={0},
	cnst_L4[12]={0},cnst_R4[12]={0},
	cnst_L5[12]={0},cnst_R5[12]={0},
	cnst_L6[12]={0},cnst_R6[12]={0};

int	p_L1[30],p_R1[30],
	p_L2[24],p_R2[24],
	p_L3[22],p_R3[22],
	p_L4[12],p_R4[12],
	p_L5[12],p_R5[12],
	p_L6[12],p_R6[12];
//============================================================

double V1R[30]={0}, V2R[24]={0}, V3R[22]={0}, V4R[12]={0}, V5R[12]={0}, V6R[12]={0};
double V1L[30]={0}, V2L[24]={0}, V3L[22]={0}, V4L[12]={0}, V5L[12]={0}, V6L[12]={0};

int run=-1;
cout<<" It's cosmic run ? (yes =1 no =0) "<<endl;
cin>>run;

double DL,DR, A;

if(run==1) {double ADC[7]={0,1200,1500,1800,3000,3000,3000};}
if(run==0) {double ADC[7]={0,1200,1200,1200,1200,1200,1200};}


if(run!=1 && run!=0) {cout<<"Bad Choise"<<endl; return;}

read_gains(	&c_L1[], &cnst_L1[], &p_L1[],
		&c_L2[], &cnst_L2[], &p_L2[],
		&c_L3[], &cnst_L3[], &p_L3[],
		&c_L4[], &cnst_L4[], &p_L4[],
		&c_L5[], &cnst_L5[], &p_L5[],
		&c_L6[], &cnst_L6[], &p_L6[],
		&c_R1[], &cnst_R1[], &p_R1[],
		&c_R2[], &cnst_R2[], &p_R2[],
		&c_R3[], &cnst_R3[], &p_R3[],
		&c_R4[], &cnst_R4[], &p_R4[],
		&c_R5[], &cnst_R5[], &p_R5[],
		&c_R6[], &cnst_R6[], &p_R6[]);

//cout<<"cnst = "<< cnst_R4[6]<<endl;

for(int i=1;i<7;i++)
	{
	for(int j=0;j<nbars[i];j++)
		{
		if(i==1)
			{
			DL=ADC[i]-cnst_L1[j];
			DR=ADC[i]-cnst_R1[j];
			if(DL>=0 && c_L1[j]>0)	{V1L[j]=(DL/c_L1[j])**(1./p_L1[j]);}
			else	cout<<" Problem with layer 1, LEFT bar "<<j<<endl;
			if(DR>=0 && c_R1[j]>0)	{V1R[j]=(DR/c_R1[j])**(1/p_R1[j]);}
			else	cout<<" Problem with layer 1, RIGHT bar "<<j<<endl;	
			}

		if(i==2)
			{
			DL=ADC[i]-cnst_L2[j];
			DR=ADC[i]-cnst_R2[j];
			if(DL>=0 && c_L2[j]>0)	{V2L[j]=(DL/c_L2[j])**(1./p_L2[j]);}
			else	cout<<" Problem with layer 2, LEFT bar "<<j<<endl;
			if(DR>=0 && c_R2[j]>0)	{V2R[j]=(DR/c_R2[j])**(1./p_R2[j]);}
			else	cout<<" Problem with layer 2, RIGHT bar "<<j<<endl;	
			}

		if(i==3)
			{

			if(run==1)
				{
				if(j<6||j>15)
					{
					DL=ADC[i]-cnst_L3[j];
					DR=ADC[i]-cnst_R3[j];
					if(DL>=0 && c_L3[j]>0)	{V3L[j]=(DL/c_L3[j])**(1./p_L3[j]);}
					else	cout<<" Problem with layer 3, LEFT bar "<<j<<endl;
					if(DR>=0 && c_R3[j]>0)	{V3R[j]=(DR/c_R3[j])**(1./p_R3[j]);}
					else	cout<<" Problem with layer 3, RIGHT bar "<<j<<endl;
					}
				if((j>=6 && j<=9)||(j>=12 && j<=15))
					{
					DL=ADC[i-1]-cnst_L3[j];
					DR=ADC[i-1]-cnst_R3[j];
					if(DL>=0 && c_L3[j]>0)	{V3L[j]=(DL/c_L3[j])**(1./p_L3[j]);}
					else	cout<<" Problem with layer 3, LEFT bar "<<j<<endl;
					if(DR>=0 && c_R3[j]>0)	{V3R[j]=(DR/c_R3[j])**(1./p_R3[j]);}
					else	cout<<" Problem with layer 3, RIGHT bar "<<j<<endl;
					}
				if(j==10 || j==11)
					{
					DL=ADC[i-2]-cnst_L3[j];
					DR=ADC[i-2]-cnst_R3[j];
					if(DL>=0 && c_L3[j]>0)	{V3L[j]=(DL/c_L3[j])**(1./p_L3[j]);}
					else	cout<<" Problem with layer 3, LEFT bar "<<j<<endl;
					if(DR>=0 && c_R3[j]>0)	{V3R[j]=(DR/c_R3[j])**(1./p_R3[j]);}
					else	cout<<" Problem with layer 3, RIGHT bar "<<j<<endl;
					}
				}
			else
				{
				DL=ADC[i]-cnst_L3[j];
				DR=ADC[i]-cnst_R3[j];
				if(DL>=0 && c_L3[j]>0)	{V3L[j]=(DL/c_L3[j])**(1./p_L3[j]);}
				else	cout<<" Problem with layer 3, LEFT bar "<<j<<endl;
				if(DR>=0 && c_R3[j]>0)	{V3R[j]=(DR/c_R3[j])**(1./p_R3[j]);}
				else	cout<<" Problem with layer 3, RIGHT bar "<<j<<endl;
				}
			
			}

		if(i==4)
			{
			DL=ADC[i]-cnst_L4[j];
			DR=ADC[i]-cnst_R4[j];
			if(DL>=0 && c_L4[j]>0)	{V4L[j]=(DL/c_L4[j])**(1./p_L4[j]);}
			else	cout<<" Problem with layer 4, LEFT bar "<<j<<endl;
			if(DR>=0 && c_R4[j]>0)	{V4R[j]=(DR/c_R4[j])**(1./p_R4[j]);}
			else	cout<<" Problem with layer 4, RIGHTbar "<<j<<endl;
			}

		if(i==5)
			{
			DL=ADC[i]-cnst_L5[j];
			DR=ADC[i]-cnst_R5[j];
			if(DL>=0 && c_L5[j]>0)	{V5L[j]=(DL/c_L5[j])**(1./p_L5[j]);}
			else	cout<<" Problem with layer 5, LEFT bar "<<j<<endl;
			if(DR>=0 && c_R5[j]>0)	{V5R[j]=(DR/c_R5[j])**(1./p_R5[j]);}
			else	cout<<" Problem with layer 5, RIGHTbar "<<j<<endl;
			}

		if(i==6)
			{
			DL=ADC[i]-cnst_L6[j];
			DR=ADC[i]-cnst_R6[j];
			if(DL>=0 && c_L6[j]>0)	{V6L[j]=(DL/c_L6[j])**(1./p_L6[j]);}
			else	cout<<" Problem with layer 6, LEFT bar "<<j<<endl;
			if(DR>=0 && c_R6[j]>0)	{V6R[j]=(DR/c_R6[j])**(1./p_R6[j]);}
			else	cout<<" Problem with layer 6, RIGHTbar "<<j<<endl;
			}

		} // for j
	}//for i


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

char *lineL, *lineR, *vL, *vR;
lineL=(char*)malloc(1000);
lineR=(char*)malloc(1000);
vL=(char*)malloc(1000);
vR=(char*)malloc(1000);

ofstream fileL, fileR;
remove("Voltages_left.dat");
remove("Voltages_right.dat");
for(i=1;i<7;i++)
	{
	fileL.open("Voltages_left.dat", ios::app);
	fileR.open("Voltages_right.dat",ios::app);
	sprintf(lineL," Left side, Layer %d : ",i);
	sprintf(lineR," Right side, Layer %d : ",i);
	fileL<<lineL;
	fileR<<lineR;
	for(j=0;j<nbars[i];j++)
		{
		if(i==1)
			{
			sprintf(vL," %.1f ",V1L[j]);
			sprintf(vR," %.1f ",V1R[j]);
			fileL<<vL;
			fileR<<vR;
			}

		if(i==2)
			{
			sprintf(vL," %.1f ",V2L[j]);
			sprintf(vR," %.1f ",V2R[j]);

			fileL<<vL;
			fileR<<vR;
			}

		if(i==3)
			{
			sprintf(vL," %.1f ",V3L[j]);
			sprintf(vR," %.1f ",V3R[j]);

			fileL<<vL;
			fileR<<vR;
			}

		if(i==4)
			{
			sprintf(vL," %.1f ",V4L[j]);
			sprintf(vR," %.1f ",V4R[j]);

			fileL<<vL;
			fileR<<vR;
			}

		if(i==5)
			{
			sprintf(vL," %.1f ",V5L[j]);
			sprintf(vR," %.1f ",V5R[j]);

			fileL<<vL;
			fileR<<vR;
			}

		if(i==6)
			{
			sprintf(vL," %.1f ",V6L[j]);
			sprintf(vR," %.1f ",V6R[j]);

			fileL<<vL;
			fileR<<vR;
			}
		
		}
	fileL<<endl;
	fileR<<endl;
	fileL.close();
	fileR.close();
	}
			
}//void