#include <iostream.h>

int xft_res(Char_t* filename="119632.hbk")
{
	//Clear out variables on the stack (avoids much grief)
	gROOT->Reset();
	gStyle->SetHistLineWidth(1);	
	gStyle->SetPadGridX(true);
	gStyle->SetPadGridY(true);
	gStyle->SetLabelSize(0.05,"XY");
	gStyle->SetTitleSize(0.08, "XY");
	gStyle->SetTitleOffset(0.4, "XY");
	gStyle->SetFuncColor(6);
	gStyle->SetFuncWidth(3.);
	
	gStyle->SetOptStat(0); //Integral, Overflow, Underflow, RMS, Mean, Nent, Name
	gStyle->SetOptFit(11);     //probability, Chi2, errors, name/values of parameters
	gStyle->SetStatX(0.995);		//Stat box x position (top right hand corner)	
	gStyle->SetStatY(0.995); 		//Stat box y position 		
	gStyle->SetStatW(0.20);	 		//Stat box width as fraction of pad size	
	gStyle->SetStatH(0.17);	 		//Size of each line in stat box	

   gStyle->SetTextFont(62);
   gStyle->SetTextColor(4);
   gStyle->SetTextSize(0.06);
   gStyle->SetTextAlign(12);

	//Open histogram file	
	cout << "filename is "<< filename << endl;	
 	TFile* f= new TFile(filename);

	//Open file in which to store pictures
	TPostScript *ps = new TPostScript("XFT_resolution.ps",111);

	//Open Canvas
	TCanvas *myc = new TCanvas("MyC","Test Canvas",700,700);
	
//-------------------------------------------------------------------
//Pt and phi resolution
	TH1F *h2001 = (TH1F*)f->Get("dPtHist");
	TH1F *h2002 = (TH1F*)f->Get("dPhiHist");
	TH2F *h2003 = (TH2F*)f->Get("dPtBinPhiHist");

	TH1F *h2011 = (TH1F*)f->Get("MdPtHist");
	TH1F *h2012 = (TH1F*)f->Get("MdPhiHist");
	TH2F *h2013 = (TH2F*)f->Get("dPtBinPhiNeg");

	TH1F *h2021 = (TH1F*)f->Get("PdPtHist");
	TH1F *h2022 = (TH1F*)f->Get("PdPhiHist");
	TH2F *h2023 = (TH2F*)f->Get("dPtBinPhiPos");

	TH2F *h3000 = (TH2F*)f->Get("offphid0");

	h2001->SetXTitle("pT resolution (1/GeV/c)");
 	h2001->SetYTitle("Events");	
	h2011->SetXTitle("pT resolution (1/GeV/c)");
 	h2011->SetYTitle("Events");	
	h2021->SetXTitle("pT resolution (1/GeV/c)");
 	h2021->SetYTitle("Events");	

	h2002->SetXTitle("phi resolution (radians)");
 	h2002->SetYTitle("Events");	
	h2012->SetXTitle("phi resolution (radians)");
 	h2012->SetYTitle("Events");	
	h2022->SetXTitle("phi resolution (radians)");
 	h2022->SetYTitle("Events");	

	h3000->SetXTitle("offline phi (radians)");
 	h3000->SetYTitle("offline d0 (cm)");	


	ps->NewPage();	
	myc->Clear();
 	myc->Divide(1,3);
  
	myc->cd(1);
	h2001->Fit("gaus");
	
	char unit[20], quantity[50];
	sprintf(unit,"%%/GeV/c"); sprintf(quantity,"Transverse momentum");
	xft_showstats(h2001,quantity,unit,100);

	myc->cd(2);
	h2002->Fit("gaus");
	sprintf(unit,"milliradians"); sprintf(quantity,"Phi");
	xft_showstats(h2002,quantity,unit,1000);


	myc->cd(3);
	h2003->ProfileX("px2003");
	TF1 *f1 = new TF1("f1","[0]*sin(x+[1])+[2]",0.,6.29);
	px2003->SetMaximum(.3);
	px2003->SetMinimum(-.3);
	px2003->SetXTitle("Phi (radians)");
 	px2003->SetYTitle("L3-XFT PtBin");	
	px2003->Fit("f1");


	myc->Update();
//---------------------------------------------------------------------------------
	char *ch = new char[1]; printf("Hit Enter to continue"); gets(ch);
	ps->NewPage();	
	myc->Clear();
 	myc->Divide(1,3);
  
	myc->cd(1);
	h2011->Fit("gaus");
	sprintf(unit,"%%/GeV/c"); sprintf(quantity,"Transverse momentum");
	xft_showstats(h2011,quantity,unit,100);

	myc->cd(2);
	h2012->Fit("gaus");
	sprintf(unit,"milliradians"); sprintf(quantity,"Phi");
	xft_showstats(h2012,quantity,unit,1000);


	myc->cd(3);
	h2013->ProfileX("px2013");
	px2013->SetMaximum(-0.0);
	px2013->SetMinimum(-1.5);
	px2013->SetXTitle("Phi (radians)");
 	px2013->SetYTitle("L3-XFT PtBin");	
	px2013->Fit("f1");

	myc->Update();
//---------------------------------------------------------------------------------
	char *ch = new char[1]; printf("Hit Enter to continue"); gets(ch);
	ps->NewPage();	
	myc->Clear();
 	myc->Divide(1,3);
  
	myc->cd(1);
	h2021->Fit("gaus");
	sprintf(unit,"%%/GeV/c"); sprintf(quantity,"Transverse momentum");
	xft_showstats(h2021,quantity,unit,100);

	myc->cd(2);
	h2022->Fit("gaus");
	sprintf(unit,"milliradians"); sprintf(quantity,"Phi");
	xft_showstats(h2022,quantity,unit,1000);


	myc->cd(3);
	h2023->ProfileX("px2023");
	px2023->SetMaximum(1.5);
	px2023->SetMinimum(0.0);
	px2023->SetXTitle("Phi (radians)");
 	px2023->SetYTitle("L3-XFT PtBin");	
	px2023->Fit("f1");


	myc->Update();
//---------------------------------------------------------------------------------
	char *ch = new char[1]; printf("Hit Enter to continue"); gets(ch);
	ps->NewPage();	
	myc->Clear();
 	myc->Divide(1,2);
  
	myc->cd(1);
	h3000->Draw();
	
	myc->cd(2);
	h3000->ProfileX("px3000");
	px3000->SetMaximum(.6);
	px3000->SetMinimum(-.6);
	px3000->SetXTitle("Phi (radians)");
 	px3000->SetYTitle("Offline d0 (cm)");		
	px3000->Fit("f1");
	xft_showbeamposition(px3000);

	myc->Update();
	ps->Close();
	char *ch = new char[1]; printf(" Hit Enter to continue \n");	gets(ch);	

	gStyle->SetStatX(0.995);		//Stat box x position (top right hand corner)	
	gStyle->SetStatY(0.995); 		//Stat box y position 		
	gStyle->SetStatW(0.09);	 		//Stat box width as fraction of pad size	
	gStyle->SetStatH(0.09);	 		//Size of each line in stat box	
	return 0;
}
int xft_showstats(TH1* h1, char* quantity, char* units, Float_t mult=1)
{	
	Double_t mean  = mult * h1->GetFunction("gaus")->GetParameter(1);
	Double_t emean = mult * h1->GetFunction("gaus")->GetParError(1);
	Double_t res   = mult * h1->GetFunction("gaus")->GetParameter(2);
	Double_t eres  = mult * h1->GetFunction("gaus")->GetParError(2);
	
	TText *t = new TText();
 	t->SetNDC();
	char toffset[100], tres[100];
	sprintf(toffset,"%5.2f +- %6.3f %s",mean,emean,units);
	sprintf(tres,"%5.2f +- %6.3f %s",res,eres,units);
	t->DrawText(0.12,.85,quantity);
	t->DrawText(0.12,.75,"Offset");
   t->DrawText(0.22,.75,toffset);
	t->DrawText(0.12,0.65,"Resolution");
  	t->DrawText(0.22,.65,tres);


 return 0;
}
int xft_showbeamposition(TH1* h1)
{	
	Double_t ampl  = h1->GetFunction("f1")->GetParameter(0);
	Double_t eampl = h1->GetFunction("f1")->GetParError(0);
	Double_t phase  = h1->GetFunction("f1")->GetParameter(1);
	Double_t ephase = h1->GetFunction("f1")->GetParError(1);
	Double_t off   = h1->GetFunction("f1")->GetParameter(2);
	Double_t eoff  = h1->GetFunction("f1")->GetParError(2);
	
	phase *= (180./3.1416);
	//Get round root's crazy fitting (amplitude can be negative!)
	if(ampl < 0.){
		ampl = -ampl;
		phase = phase + 180.;
		if(phase > 360.) phase = phase - 360.;
	}
		
	Double_t angle = -phase + 360.;
	angle = angle + 180.;
	if(angle>360.) angle = angle -360.;
	
	Double_t bx = fabs(ampl)*cos(angle);
	Double_t by = fabs(ampl)*sin(angle);
	Double_t ebx = sqrt( sin(angle)*sin(angle)*eampl*eampl +
	 			  ampl*ampl*cos(angle)*cos(angle)*ephase*ephase);
	Double_t eby = sqrt( cos(angle)*cos(angle)*eampl*eampl +
	 			  ampl*ampl*sin(angle)*sin(angle)*ephase*ephase);

	char temp1[100],temp2[100],temp3[100];
	sprintf(temp1,"radius = %6.3f +- %6.3f cm",ampl,eampl);
	sprintf(temp2,"angle = %6.3f +- %6.3f degrees",angle,ephase);	
//	sprintf(temp3,"Position (%6.3f+-%6.3f, %6.3f+-%6.3f) cm",bx,ebx,by,eby);
	
	printf("%s, %s, %s",temp1,temp2,temp3);
	TText *t = new TText();
   t->SetTextSize(0.04);
	t->SetNDC();
	t->DrawText(0.35,.85,temp1);
	t->DrawText(0.35,.75,temp2);
	t->DrawText(0.35,.65,temp3);



 return 0;
}