{
/*

This script is only for HAND, for 112 scintillators.


This script is used for determination of the ADC peaks and writing it to the file, witht HV values.
In the code, the ADC is plotted vs position and corrected. Than the peak is determined.
Peak is determined by cutting on top and bottom of the bar of interest. For top and bottom using two bars above and two bars bellow.

The script is using two function read_voltages.C and read_hv_mapping.C

This script is giving two kind of output files: test.dat and ADC_right.dat/ ADC_left.dat.

The first file store the data for curve.C in order to calculate the gain curves for this runs.

This script is applied for one layer each time ( 3 or 4).


The second output file is used for adjust.C script. This output is relevant AFTER the gain curves found, and correct voltages are calculated.

*/


#include <fstream>
#include <iostream>
#include "read_voltages.C"
#include "read_hv_mapping.C"


//void adc(){

gROOT->Reset();
gStyle->SetPalette(1);
gStyle->SetOptFit();

int num,a,b;

cout<<" Enter Layer number "<< endl;

cin>>num;

if(num==1)	{a=10; const b=30;}
if(num==2)	{a=6;  const b=24;} 
if(num==3)	{a=6;  const b=22;}
if(num==4)	{a=4;  const b=12;}
if(num==5)	{a=4;  const b=12;}
if(num==6)	{a=4;  const b=12;}

TCanvas c1;		// for position distribution
c1.Divide(3,a);
TCanvas c2;		// used for corrected ADCs VS position distribution Left
c2.Divide(3,a);
TCanvas c3;		// used for corrected ADCs VS position distribution Right
c3.Divide(3,a);
TCanvas c4;		// used for corrected ADCs distribution left
c4.Divide(3,a);
TCanvas c5;		// used for corrected ADCs distribution right
c5.Divide(3,a);
TCanvas c6;		// Left TDC distribution
c6.Divide(3,a);
TCanvas c7;		// Right TDC distribution
c7.Divide(3,a);
TCanvas c8;		// Left ADC distribution
c8.Divide(3,a);
TCanvas c9;		// Right ADC distribution
c9.Divide(3,a);
TCanvas c10;		// tmp, for profile histograms

TF1 *fit=new TF1("l","landau",800,8000);
TF1 *fit1=new TF1("profL","[0]*exp(x*[1])",-50,50);
TF1 *fit2=new TF1("profR","[0]*exp(x*[1])",-50,50);

double peakL[30]={0},mL[30]={0};
double peakR[30]={0},mR[30]={0};
double pedL[30]={353, 275, 287, 341, 323, 437, 467, 455, 377, 287, 311, 251,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double pedR[30]={329, 383, 305, 299, 467, 473, 683, 401, 557, 539, 539, 341,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double peakRerr[30], peakLerr[30];

char lineL[200],lineR[200],p[200],p_cut[200],left[300],right[300],cut[200],left_1D[200],right_1D[200],lL[200],lR[200],c_lineL[200],c_lineR[200],PROFILE_L[200],PROFILE_R[200];
double pos[30]={0};

double meanL[30]={0},meanR[30]={0},cL[30]={0},cR[30]={0};

for(int i=0;i<b;i++)
	{

	TH1F *tL=new TH1F("tL","tL",100,10,400);
	TH1F *tR=new TH1F("tR","tR",100,10,400);
	TH1F *aL=new TH1F("aL","aL",100,10,20000);
	TH1F *aR=new TH1F("aR","aR",100,10,20000);
	TH1F *tmp=new TH1F("tmp","tmp",100,-50,50);
	TH2F *hL=new TH2F("hL","hL",100,-50,50,100,0,8000);
	TH2F *hR=new TH2F("hR","hR",100,-50,50,100,0,8000);
	TH1F *HL=new TH1F("HL","HL",100,0,7000);
	TH1F *HR=new TH1F("HR","HR",100,0,7000);
	hprofL = new TProfile("hprofL","hprofL",100,-50,50,2500,8000);
	hprofR = new TProfile("hprofR","hprofR",100,-50,50,2500,8000);

	sprintf(lineL,"BB.HP4.lt[%d]>>tL",i);
	sprintf(lineR,"BB.HP4.rt[%d]>>tR",i);

	sprintf(c_lineL,"BB.HP4.lt[%d]>250",i);
	sprintf(c_lineR,"BB.HP4.rt[%d]>250",i);

	sprintf(lL,"BB.HP4.la[%d]>>aL",i);
	sprintf(lR,"BB.HP4.ra[%d]>>aR",i);

	c6.cd(i+1);
	c6.SetLogy(1);
	T->Draw(lineL,c_lineL);
	c7.cd(i+1);
	c7.SetLogy(1);
	T->Draw(lineR,c_lineR);
	c8.cd(i+1);
	T->Draw(lL);
	c9.cd(i+1);
	T->Draw(lR);
	

	sprintf(p,"BB.HP4.lt[%d]-BB.HP4.rt[%d]>>tmp",i,i);
	sprintf(p_cut,"BB.HP4.lt[%d]>0 && BB.HP4.rt[%d]>0",i,i);

	sprintf(PROFILE_L,"BB.HP4.la[%d]-%f:(BB.HP4.lt[%d]-BB.HP4.rt[%d])>>hprofL",i,pedL[i],i,i);
	sprintf(PROFILE_R,"BB.HP4.ra[%d]-%f:(BB.HP4.lt[%d]-BB.HP4.rt[%d])>>hprofR",i,pedR[i],i,i);

	//Draw position distribution
	c1.cd(i+1);
	T->Draw(p,p_cut);
	pos[i]=tmp->GetMean();
	//=======================

	if(i==0)
		{
		sprintf(cut,"BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200",i+1,i+1,i+2,i+2,i,i);
		}
	
	if(i==b-1)
		{
		sprintf(cut,"BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200",i-1,i-1,i-2,i-2,i,i);
		}
	if(i!=0 && i!=11)
		{
		sprintf(cut,"BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200 && BB.HP4.lt[%d]>200 && BB.HP4.rt[%d]>200",i-1,i-1,i+1,i+1,i,i);	
		}


//===================== Draw a profile histogram in order to fit ADCs position dependence ========================
	c10.cd();
	T->Draw(PROFILE_L,cut);
	hprofL->Fit(profL,"QR");
	meanL[i]=hprofL->GetMean();
	cL[i]=profL->GetParameter(1);

	T->Draw(PROFILE_R,cut);
	hprofR->Fit(profR,"QR");
	meanR[i]=hprofR->GetMean();
	cR[i]=profR->GetParameter(1);
//=================================================================================================================


//=============================== Corrected ADCs Vs position  ===============================================
	sprintf(left,"(BB.HP4.la[%d]-%f)/exp((BB.HP4.lt[%d]-BB.HP4.rt[%d]+%f)*%f):(BB.HP4.lt[%d]-BB.HP4.rt[%d])>>hL",i,pedL[i],i,i,meanL[i],cL[i],i,i);
	sprintf(right,"(BB.HP4.ra[%d]-%f)*exp((-BB.HP4.lt[%d]+BB.HP4.rt[%d]-%f)*%f):(BB.HP4.lt[%d]-BB.HP4.rt[%d])>>hR",i,pedR[i],i,i,meanR[i],cR[i],i,i);
	

	// Draw corrected ADCs Vs position
	c2.cd(i+1);
	T->Draw(right,cut);
	c3.cd(i+1);
	T->Draw(left,cut);
	//===========================
//===============================================================================================================




//========================================= Corrected ADCs distribution ===========================================
	sprintf( left_1D,"(BB.HP4.la[%d]-%f)/exp((BB.HP4.lt[%d]-BB.HP4.rt[%d]+%f)*%f)>>HL",i,pedL[i],i,i,meanL[i],cL[i],i,i);
	sprintf(right_1D,"(BB.HP4.ra[%d]-%f)*exp((-BB.HP4.lt[%d]+BB.HP4.rt[%d]-%f)*%f)>>HR",i,pedR[i],i,i,meanR[i],cR[i],i,i);

	c4.cd(i+1);
	T->Draw(left_1D,cut);
	mL[i]=HL->GetMean();
	HL->Fit(l,"QR");
	peakL[i]=l->GetParameter(1);
	peakLerr[i]=l->GetParError(1);


	c5.cd(i+1);
	T->Draw(right_1D,cut);
	mR[i]=HR->GetMean();
	HR->Fit(l,"QR");	
	peakR[i]=l->GetParameter(1);
	peakRerr[i]=l->GetParError(1);

	} // for

//=============================================  Read HV for this Layer  =======================================


double HV1[16][12]={0}, HV2[16][12]={0};	// variables for read_voltages.C
double VL[112]={0},VR[112];

int CrateL[112]={-1}, SlotL[112]={-1}, ChanL[112]={-1}; 
int CrateR[112]={-1}, SlotR[112]={-1}, ChanR[112]={-1};
int Layer[112]={-1},Bar[112]={-1};

read_hv_mapping(CrateL,SlotL,ChanL,CrateR,SlotR,ChanR,Layer,Bar);
read_voltages(HV1,HV2);

int n=0,m=0;
for(int k=0;k<112;k++)
	{
	if(Layer[k]==num)
		{

		if(CrateL[k]==1)
			{
			
			VL[n]=HV1[SlotL[k]][ChanL[k]];
			n++;
			}

		if(CrateL[k]==2)
			{
			VL[n]=HV2[SlotL[k]][ChanL[k]];
			n++;
			}

		if(CrateR[k]==1)
			{

			VR[m]=HV1[SlotR[k]][ChanR[k]];
			m++;
			}

		if(CrateR[k]==2)
			{
			VR[m]=HV2[SlotR[k]][ChanR[k]];
			m++;
			}
		}
	}
	
		

//======================================= Output to file ===================================================

ofstream file,fileL,fileR;
file.open("test.dat", ios::ate);

fileL.open("ADC_left.dat");
fileR.open("ADC_right.dat");

//============== left side ==========
file<<"double peakL[]= {";
for(j=0;j<b;j++)
	{
	file<<peakL[j];
	fileL<<peakL[j]<<endl;

	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;
file<<"double peakLerr[]= {";




for(j=0;j<b;j++)
	{
	file<<peakLerr[j];
	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;

file<<"double VoltL[]= {";
fileL<<"HV"<<endl;

for(j=0;j<b;j++)
	{
	file<<VL[j];
	fileL<<VL[j]<<endl;
	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;


fileL.close();
file<<"\n"<<endl;



//============== right side ===========
file<<"double peakR[]= {";
for(j=0;j<b;j++)
	{
	file<<peakR[j];
	fileR<<peakR[j]<<endl;
	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;
file<<"double peakRerr[]= {";


for(j=0;j<b;j++)
	{
	file<<peakRerr[j];
	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;

file<<"double VoltR[]= {";
fileR<<"HV"<<endl;
for(j=0;j<b;j++)
	{
	file<<VR[j];
	fileR<<VR[j]<<endl;
	if(j!=b-1)	file<<", ";
	}
file<<"};"<<endl;


file.close();
fileR.close();


}