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

This script calculate the position offset for HAND scintillators.
The variables in this function matching the THaScintillator.C module.

The output is writen to file: shifts.dat
The structure of the output matching the db_NA.nd.p_.dat data base file.

Pay attention that the position offset in the data base file is defined from the starting edge of the fit and not from the center
this is due to the hand libraries for the analyzer


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



Double_t fitfun(Double_t *x, Double_t *par)
	{
        if ((x[0]>par[0])  && (x[0]<par[0]+par[1])) 
		{
	
        	 Double_t fitval = par[2];

		} 
	else Double_t fitval=0;

      	return fitval;
   	}
	

void hand_position_offsets()
	{
gStyle->SetOptFit(1);

//	TFile *f=TFile::Open("~/E01015/ROOTfiles/e01015_ND_1704.root");
	TFile *f=TFile::Open("~/E01015/Rootfiles/1704_all_0.root");
/*	TFile *f=TFile::Open("~/e05102/ROOTfiles/e05102_phys_matrix3488.root");
	TChain *c=new TChain("T");
	c.Add("~/e05102/ROOTfiles/e05102_phys_3488_a.root");
	c.Add("~/e05102/ROOTfiles/e05102_phys_3488_1.root");
	c.Add("~/e05102/ROOTfiles/e05102_phys_3488_2.root");
	c.Add("~/e05102/ROOTfiles/e05102_phys_3488_3.root");
	c.Add("~/e05102/ROOTfiles/e05102_phys_3488_4.root");
*/
	TCanvas c1;
	Float_t Shifts[32], Scales[32],m;
	int channels;
	char *dummy,*dummy1,*str;
	dummy=(char *)malloc(1000);
	dummy1=(char *)malloc(1000);
	str=(char *)malloc(1000);
	TH1F H("H","H",100,-50,50);
	TF1 *func = new TF1("fit",fitfun,-30,30,3);
	func->SetParNames("Start","Width","Height");

int pl,evtype;
cout<<" ENTER PLANE NUMBER "<<endl;
cin>>pl;
cout<<" ENTER EVENTTYPE  (singels)  "<<endl;
cin>>evtype;

if(pl==1) channels=30;
if(pl==2) channels=24;
if(pl==3) channels=22;
if(pl==4) channels=12;
if(pl==5) channels=12;
if(pl==6) channels=12;


	for (Int_t i=0;i<channels;i++) 
		{
	       
		sprintf(dummy, "NA.nd.p%d.lt_c[%d]-NA.nd.p%d.rt_c[%d]>>H",pl,i,pl,i);
		sprintf(dummy1,"NA.nd.p%d.rt_c[%d]>0 && NA.nd.p%d.lt_c[%d]>0 && fEvtHdr.fEvtType==%d",pl,i,pl,i,evtype);
        		
		T->Draw(dummy,dummy1);
		m=H.GetMean();
		TH1F h("h","h",60,m-30,m+30);
		sprintf(dummy, "NA.nd.p%d.lt_c[%d]-NA.nd.p%d.rt_c[%d]>>h",pl,i,pl,i);
		T->Draw(dummy,dummy1);


gets(dummy);
		func->SetParameters(h->GetMean()-1.5*h->GetRMS(),h->GetRMS()*3,h->GetMaximum()/2.);
		h->Fit("fit");
		Scales[i]=func->GetParameter(1);
        	Shifts[i]=-(func->GetParameter(0)+Scales[i]/2.);
		c1.Paint();
		c1.Update();
	        printf("%i\n",i);
		gets(dummy);  
      		}
      
	ofstream file;
	file.open("shifts.dat");
	for (i=0;i<channels;i++) 
		{
		sprintf(str,"%.3f   %.3f                                 - Paddle %i   (Left/Right)",Shifts[i]+Scales[2]/2,Scales[i],i+1);
		file<<str<<endl;
		cout<<"  bar --- "<<i+1<<"  Shift === "<<Shifts[i]<<"\n"<<endl;
		cout<<"            bar --- "<<i+1<<"  Scales === "<<Scales[i]<<"\n"<<endl;

      		}
	file.close();

	printf("\n");
    
}