{
/*// =================================================================================================================

	Checking the gain of the scintillators. In this program we plot dE vs E in order to see punch through point.
	The program runs over dE ADC's channels and ploting the distribution of the E plane ADC. The program takes geometric sum of the dE ADCs and ploting left and right side of E bar. In addition position cut on E bar is used.
	The program calculate the mean position on the E plane ADC distribution or it can find gausian fit for the peak location. this is choosen by "type".
	The punch through point is defined (in this script) as maximum peak position of E bar.

*/// =================================================================================================================

gROOT->Reset();

double nchannels=4000;
double step=50;		// Steps over dE plane
double first_channel=150;

char *lineL,*lineR,*cut,*d,*strL,*strR,*offset,*c_offset;

lineL=(char*)malloc(1000);	// ADC distribution of left side ofE bar PMT. 
lineR=(char*)malloc(1000);	// ADC distribution of right side of E bar 
cut=(char*)malloc(1000);	// cut on the dE strip, position cut of the E bar and event type cut
d=(char*)malloc(1000);		 
strL=(char*)malloc(1000);	// used for drawing dE/E(left)
strR=(char*)malloc(1000);	//  used for drawing dE/E(right)
offset=(char*)malloc(1000);	// used to find position offset LT - RT
c_offset=(char*)malloc(1000);	// cut for position offset

double numL[40],adc_chL[40],maxL=-1,heightL;	// maxL = ADC position of the punch through for E left side. heightL - same for dE 
double numR[40],adc_chR[40],maxR=-1,heightR;	// maxR = ADC position of the punch through for E right side. heightR - same for dE

double pos=first_channel;	//first ADC channel on the dE plane, used to cut off the minimum ionizing
double off;			// position offset of the bar
int bar;			// bar to make analysis
int evtype=1;			// event type, singles
int type=2;			// method for punch through point determination. 1- gauss fit, 2 - mean value
int k=0;

TH1 *hl=new TH1F("hl","hl",100,100,2500);
TH1 *hr=new TH1F("hr","hr",100,100,2500);
TH1 *Hpos=new TH1F("Hpos","Hpos",50,-10,10);
TF1 *g=new TF1("g","gaus",150,2500);



cout<<" Enter bar to analyze "<<endl;

cin>>bar;


TCanvas c1;
sprintf(offset,"BB.tp.e.LT[%d]-BB.tp.e.RT[%d]>>Hpos",bar,bar);
sprintf(c_offset,"BB.tp.e.LT[%d]>-10000 && BB.tp.e.RT[%d]>-10000",bar,bar);
T->Draw(offset,c_offset);
off=Hpos->GetMean();


TCanvas c2;
c2.Divide(2,1);




// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = // 
for(int i=0;i<15;i++)
	{
	cout<<"Position === "<<pos<<endl;
	sprintf(lineL,"BB.tp.e.LApedc[%d]>>hl",bar);
	sprintf(lineR,"BB.tp.e.RApedc[%d]>>hr",bar);

	sprintf(cut,"TMath::Abs(sqrt(BB.tp.de.RApedc[%d]*BB.tp.de.LApedc[%d])-%f)<25 && fEvtHdr.fEvtType==%d && TMath::Abs(BB.tp.e.LT[%d]-BB.tp.e.RT[%d]-%f)<1",bar,bar,pos,evtype,bar,bar,off);

	c2.cd(1);
	T->Draw(lineL,cut);
	if(hl->GetEntries()>10)
		{
		

		if(type==1)
			{
			hl->Fit("g","RQ");
			adc_chL[k]=g->GetParameter(1);
			TLine *ll=new TLine(adc_chL[k],0,adc_chL[k],g->GetParameter(0));
			}
		if(type==2)
			{
			adc_chL[k]=hl->GetMean();
			TLine *ll=new TLine(adc_chL[k],0,adc_chL[k],hl->GetMaximum());
			}
		ll->SetLineColor(2);
		ll->SetLineWidth(3);
		ll->Draw();
		numL[k]=pos;
		c2.Paint();
		c2.Update();
		
		} //if

	c2.cd(2);
	T->Draw(lineR,cut);
	if(hr->GetEntries()>10)
		{
		
		if(type==1)
			{
			hr->Fit("g","RQ");
			adc_chR[k]=g->GetParameter(1);
			TLine *lr=new TLine(adc_chR[k],0,adc_chR[k],g->GetParameter(0));
			}
		if(type==2)
			{
			adc_chR[k]=hr->GetMean();
			TLine *lr=new TLine(adc_chR[k],0,adc_chR[k],hr->GetMaximum());
			}

	
		lr->SetLineColor(2);
		lr->SetLineWidth(3);
		lr->Draw();
		numR[k]=pos;
		c2.Paint();
		c2.Update();
		//gets(d);
		} //if

	

	k++;
	pos=pos+step;
		
	} //for

// ================================ calculate the maxR and maxL ===================//
for(int j=0;j<=k;j++)
	{
	if(adc_chL[j]>maxL) 
		{
		maxL=adc_chL[j];
		heightL=numL[j];
		}
	if(adc_chR[j]>maxR) 
		{
		maxR=adc_chR[j];
		heightR=numR[j];
		}
	}



//================= Draw the dE/E plot for left and right side with lines pointing on the punch through ppoints =========//
TCanvas c3;

c3.Divide(2,1);

TH2 *hL=new TH2F("hL","hL",100,0,3000,100,0,1500);
TH2 *hR=new TH2F("hR","hR",100,0,3000,100,0,1500);

sprintf(strL,"TMath::Sqrt(BB.tp.de.RApedc[%d]*BB.tp.de.LApedc[%d]):BB.tp.e.LApedc[%d]>>hL",bar,bar,bar,bar);
sprintf(strR,"TMath::Sqrt(BB.tp.de.RApedc[%d]*BB.tp.de.LApedc[%d]):BB.tp.e.RApedc[%d]>>hR",bar,bar,bar,bar);

c3.cd(1);
T->Draw(strL);
TLine *lL=new TLine(maxL,0,maxL,heightL);
TLine *lL1=new TLine(0,heightL,maxL,heightL);
lL->SetLineColor(2);
lL1->SetLineColor(2);
lL->SetLineWidth(3);
lL1->SetLineWidth(3);
lL->Draw();
lL1->Draw();

c3.cd(2);
T->Draw(strR);
TLine *lR=new TLine(maxR,0,maxR,heightR);
TLine *lR1=new TLine(0,heightR,maxR,heightR);
lR->SetLineColor(2);
lR1->SetLineColor(2);
lR->SetLineWidth(3);
lR1->SetLineWidth(3);
lR->Draw();
lR1->Draw();

}