#include <cmath>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TKey.h>
#include <TFile.h>
#include <TF1.h>
#include <TLatex.h>
#include <TPaveLabel.h>
#include <TTree.h>
#include <TApplication.h>
#include <TLine.h>
#include <TGraph2D.h>
#include <TBrowser.h>
// to use NR routines
//#include "nr.h"
//#include "nrutil.h"

using namespace std;

//      Global variables
	Int_t  nbins;				// ROOT parameters for Histos
	unsigned short  data_image[4194304];	// For the data in the EDF file
	int rsize=0;				// side size of image(rsize=2048)  
	unsigned short  data_image_set[40000000];//In Memory data from the EDF file


//  Function Prototypes
void myswap2(unsigned short *);
void edf2text(FILE *);//reading *edf* image into 'text' format
Double_t background(Double_t *, Double_t *);//Routines for fitting data from images
Double_t GaussianPeak(Double_t *, Double_t *);
Double_t fitFunction(Double_t *, Double_t *);
unsigned short ByteSwap4 (unsigned short );//fixing 'endian' story

/**************************************************************************
Read and Display an EDF file
***************************************************************************/

int main(int argc, char **argv) 
	{
	if(argc==1){
    printf("Usage  : %s (option) <input image file>  <setupfile>      <root output file>\n",argv[0]);
    printf("Option : nothing yet \n"); 
    exit(0); 
	}
	
	// Read setup file 
	char setup_file[64];   // setup file
	strcpy(setup_file,argv[2]); //The setup.txt file contains info about the images  
	FILE* setup_file_p=fopen(setup_file,"r");

// prepare ROOT output file 
	char root_file[64];   // setup file
	strcpy(root_file, argv[3]);
	TFile *hfile = new TFile(root_file, "RECREATE");
	
	// Prepares a CANVAS window
	TApplication theApp("App", &argc, argv);
	TCanvas *MyCanvas = new TCanvas("MyCanvas", "The Data Canvas",1050, 730);
	MyCanvas->Connect("Closed()", "TApplication", &theApp, "Terminate()");
	MyCanvas->Update();

//read the info about the images. eg number of images, step angle at which they were taken
	float ftemp;
	Int_t	nLines = 0;
	Int_t ix,iy;
	float input_parameter[20];     // for input parameter file reading
	char description[100];         // for input parameter file reading
		
	cout << "\n\n Reading setup file \n\n";

	while (!feof(setup_file_p) || nLines >> 10) 
		{
		fscanf (setup_file_p,"%f  : %s \n",&ftemp,description);
		input_parameter[nLines] = ftemp;
		cout << input_parameter[nLines] << " : " << description << "\n";
		nLines = nLines+1;
		}
		printf("\n found %d input paramters  \n\n",nLines);


// Unpack data from the setup file
	Double_t ang_step = input_parameter[0];    // angle between steps
	Double_t ang_first = input_parameter[1];    // initial angle(for first image)
	Int_t steps = (Int_t)input_parameter[2];    // number of steps(number of images)
	Double_t dispersion = input_parameter[3];    // for dispersion correction

	
	fclose(setup_file_p);

	
// Read in data, fill histogram
	  
	char input_file[128];   // name of central image file
        char input_file_root[128];//same as above
	strcpy(input_file,argv[1]); 
	sprintf(input_file_root, "%s%d.edf",input_file, steps/2);  // Do setup with the  entrral image
	FILE* input_file_p=fopen(input_file_root,"r");
	cout << " Will read the file  "  <<   input_file_root << "\n\n";
	edf2text(input_file_p);  // Data will be in variable data_image[rsize] of size rsize
	fclose(input_file_p);
	nbins=(Int_t)sqrt((double)rsize/2);   //Determining  number of bins
	cout << "The size of the 2D histogram side is " << nbins << "\n";
	
// Set up other parameters for ROOT objects
	Double_t xlow, xup, ylow, yup;   // Parameters for 2D Histos
	xlow = 0.;
	xup = nbins-1.;
	ylow = 0.;
	yup = nbins-1.;

	cout << "\n\n Making Root Objects \n";
	TH2I *h2_all = new TH2I("h2_all","Histo image",nbins,xlow,xup,nbins,ylow,yup);//2D histo from first image

	for(ix=0;ix<(Int_t)nbins;ix++)	{
		for(iy=0;iy<nbins;iy++) {
				h2_all->SetBinContent(ix+1,iy+1,data_image[ix+nbins*iy]);  // bin 0 is for underflow
			}
		}
       //Set up X and Y projections of the 2D histo h2_all 
	TH1D *h1_x = new TH1D("h1_x","Histo Projextion X",nbins,xlow,xup);
	TH1D *h1_y = new TH1D("h1_y","Histo Projextion Y",nbins,ylow,yup);
	
	MyCanvas->Divide(3,2);
	MyCanvas->cd(1);  //drawing X projection original data on the 1st canvas
	h1_x->SetLineColor(1);
	h1_x->SetStats(kFALSE);
	h1_x = h2_all->ProjectionX("X-Projection", -1, -1, "");
	h1_x->Draw();
	MyCanvas->cd(2);  //drawing Y Projection original data on the 2nd canvas
	h1_y->SetLineColor(2);
	h1_y->SetStats(kFALSE);
	h1_y = h2_all->ProjectionY("Y-Projection", -1, -1, "");
	h1_y->Draw();
	MyCanvas->Draw();
	MyCanvas->Update();

// Automatic sample image boundary detection in the images, to cut out blank edges
	Int_t i;
	Int_t h_Xmax_L, h_Xmax_R;
	Int_t h_Ymax_L, h_Ymax_R;
	Int_t found_edge = 0;
	for(i=2;i<=nbins-1;i++) {
	if ( (h1_x->GetBinContent(i) >= h1_x->GetBinContent(3)*1.3) && found_edge==0 ) {
		h_Xmax_L=i;
		found_edge=1;
		}
	}
	
	found_edge = 0;
	for(i=2;i<=nbins-1;i++) {
	if ((h1_x->GetBinContent(nbins-i) >= h1_x->GetBinContent(3)*1.3) && found_edge==0 ) {
		h_Xmax_R=nbins-i;
		found_edge=1;
		}  
	}

	found_edge = 0;
	for(i=2;i<=nbins-1;i++) {
		if (( h1_y->GetBinContent(i) >= h1_y->GetBinContent(3)*1.3) && found_edge==0 )  {
		h_Ymax_L=i;
		found_edge=1;
		}
	}
	
	found_edge = 0;
	for(i=2;i<=nbins-1;i++) {
	if ((h1_y->GetBinContent(nbins-i) >= h1_y->GetBinContent(3)*1.3) && found_edge==0 )  {
		h_Ymax_R=nbins-i;
		found_edge=1;
		} 
	}

	Int_t nbins1 = (h_Xmax_R-h_Xmax_L);
	if((h_Xmax_R-h_Xmax_L) << (h_Ymax_R-h_Ymax_L)) nbins1 = (h_Ymax_R-h_Ymax_L);
	nbins1 = nbins1 + nbins1/5;
	
	xlow = h_Xmax_L - h_Xmax_L/20;
	xup = xlow + nbins1-1;
	ylow =h_Ymax_L - h_Ymax_L/20;
	yup = ylow + nbins1-1;
	cout << "xlow = " << xlow << "    ylow = " << ylow << " nbins1 = " << nbins1 << "\n";
	
/* //  Ask user for the Crystal limts
	cout << " Enter X-Low, X-High, Y-Low, Y-High     ..... :\n";
	cin >> xlow ;
	cin >> xup ;
	cin >> ylow;
	cin >> yup;
	*/
/*
 //  Enter the Crystal limts
	xlow = 530;
	xup = 1170;
	ylow =1020;
//---------------------------------
	xlow = 750;
	xup = 850;
	ylow =1300;
	nbins1 = xup-xlow+1;
	yup = ylow + nbins1-1;
*/
        //Now we can use a smaller size image of the actual sample after cutting out the blank edges
	TH2I *h2_small = new TH2I("h2_small","Histo",nbins1,xlow,xup,nbins1,ylow,yup);
	for(int ix=(Int_t)xlow;ix<=(Int_t)xup;ix++)	{
		for(int iy=(Int_t)ylow;iy<=(Int_t)yup;iy++) {
				h2_small->SetBinContent(ix-(Int_t)xlow+1,iy-(Int_t)ylow+1,data_image[ix+nbins*iy]);  // bin 0 is for underflow
			}
		}

	MyCanvas->cd(3);  // Draw the sample image on the 3rd canvas to confirm sample edge detection worked

	gPad->SetRightMargin(0.2);
	gPad->SetLeftMargin(0.15);
	gPad->SetTopMargin(0.15);
	gPad->SetBottomMargin(0.15);
	h2_small->Print();
	h2_small->SetMinimum(h2_small->GetMinimum());
	h2_small->SetMaximum(h2_small->GetMaximum());
	h2_small->SetContour(100);
	gStyle->SetPalette(1);
	h2_small->SetStats(kFALSE);
	h2_small->Draw("colz");
	MyCanvas->Draw();
	MyCanvas->Update();
	delete h2_small;
//If everything looks fine so far, then lets read all 82 images and store them in 2D histos which we write into the output TFile
// Now set up the histos you will need.
	TH2F *Histo_2Db = new TH2F("Histo_2Db","X vs Y vs height",nbins1,xlow,xup,nbins1,ylow,yup);
	TH2F *Histo_2Dc = new TH2F("Histo_2Dc","X vs Y vs posn" ,nbins1,xlow,xup,nbins1,ylow,yup);
	TH2F *Histo_2Dd = new TH2F("Histo_2Dd","X vs Y vs width"  ,nbins1,xlow,xup,nbins1,ylow,yup);

	TH2F *Histo_2Dab = new TH2F("Histo_2Db","X vs Y vs height",nbins1,xlow,xup,nbins1,ylow,yup);
	TH2F *Histo_2Dac = new TH2F("Histo_2Dc","X vs Y vs posn" ,nbins1,xlow,xup,nbins1,ylow,yup);
	TH2F *Histo_2Dad = new TH2F("Histo_2Dd","X vs Y vs width"  ,nbins1,xlow,xup,nbins1,ylow,yup);

// create a TF1 with the range from 0 to 3 and 4 parameters
	TF1 *fitFcn = new TF1("fitFcn",fitFunction, 2* ang_step - ang_first, (steps-1)* ang_step - ang_first, 4);
   fitFcn->SetNpx(500);
   fitFcn->SetLineWidth(1);
   fitFcn->SetLineColor(kMagenta);
	

//  Read in the entire data set  .... a set of images for each tilt angle (82 angular positions)
	Int_t  ix1, iy1;
	for (i=1;i<steps+1;i++){//steps=81=umber of images
		sprintf(input_file, "%s%02d.edf",input_file_root,i);//actual hard disk image names
		FILE* input_file_p=fopen(input_file,"r");
		edf2text(input_file_p);  // Data will be in data_array of size rsize
		fclose(input_file_p);
		cout << "read file " << i << "  " << input_file << "\n";//Tell us which image number you are reading
		for(ix=(Int_t)xlow;ix<=(Int_t)xup;ix++)	{
			for(iy=(Int_t)ylow;iy<=(Int_t)yup;iy++) {
				//if( ((ix % 200 ) == 0) && ((iy % 200 ) == 0) ) 
				ix1 = ix - (int) xlow ;
				iy1 = iy - (int) ylow ;
				data_image_set[i+steps*(ix1+nbins1*iy1)] = data_image[ix+nbins*iy];
			}
		}
	}
	
	
	cout <<  " leaving the read routine for the set of images into memory \n";//finished reading all images
        //Now preparing the data for dumping into those 2D histos:
	//loop over each pixel on each image and extract a Rocking Curve and then extract the parameters for that curve
	TH1I *h1_rock_curve = new TH1I("h1_rock_curve","Rocking Curve",steps,1* ang_step - ang_first,steps* ang_step - ang_first);
	TAxis *axis = h1_rock_curve->GetXaxis();
	axis->SetRange((Int_t)(2*ang_step - ang_first), (Int_t)((steps-1)*ang_step - ang_first));
		for(ix=0;ix<nbins1;ix++)	{
			for(iy=0;iy<nbins1;iy++) {
				for(i=1;i<=steps;i++) {
				h1_rock_curve->SetBinContent(i,data_image_set[i+steps*(ix+nbins1*iy)]);
				}
				Double_t p1 = h1_rock_curve->GetMaximum();  // height of Gaussian
				Double_t p2 = h1_rock_curve->GetMean()  ; //  position of Gaussian
				Double_t p3 = h1_rock_curve->GetRMS()/2.4 ; // width of Gaussian
				if( ((ix % 50 ) == 0) && ((iy % 50 ) == 0) ) {
					MyCanvas->cd(4);
					h1_rock_curve->Draw();
					MyCanvas->Update();
	//cout <<"("<<ix<<","<<iy<<")   height= "<<p1<<"     posn= "<<p2<<"     width= "<<p3<<"\n";
				}
				
				Histo_2Dab->SetBinContent(ix+1,iy+1,p1);
				Histo_2Dac->SetBinContent(ix+1,iy+1,p2);
				Histo_2Dad->SetBinContent(ix+1,iy+1,p3);   

	// set start values for some parameters, fit and the re-extract parameters
				fitFcn->SetParameter(0,400); // background
				fitFcn->SetParameter(1,p1); // Gaussian height
				fitFcn->SetParameter(2,p2); // Gaussian peak posn
				fitFcn->SetParameter(3,p3);   //Gaussian width
				fitFcn->SetParLimits(2,p2-5*p3,p2+5*p3);
				fitFcn->SetParLimits(3,p3/2,2*p3);
				h1_rock_curve->Fit("fitFcn","=Q","",2*ang_step-ang_first,(steps-1)*ang_step- ang_first);
				p1 = fitFcn->GetParameter(1);    // height of Gaussian
				p2 = fitFcn->GetParameter(2);    // peak positionwidth of Gaussian
				p2 = p2+((float)iy)/((float)nbins1)*dispersion;	// dispersion correction
				p3 = fitFcn->GetParameter(3);    // width of Gaussian
				/*if( ((ix % 50 ) == 0) && ((iy % 50 ) == 0) ) {
					fitFcn->Draw("same");
		cout <<"("<<ix<<","<<iy<<")   height= "<<p1<<"     posn= "<<p2<<"     width= "<<p3<<"\n";
				}*/
				
				Histo_2Db->SetBinContent(ix+1,iy+1,p1);
				Histo_2Dc->SetBinContent(ix+1,iy+1,p2);
				Histo_2Dd->SetBinContent(ix+1,iy+1,p3);
            }
        }

        Histo_2Db->SetMaximum(1.02*Histo_2Db->GetMaximum());
        Histo_2Db->SetMinimum(0.90*Histo_2Db->GetMinimum());
        Histo_2Dc->SetMaximum(1.02*Histo_2Dc->GetMaximum());
        Histo_2Dc->SetMinimum(0.90*Histo_2Dc->GetMinimum());
        Histo_2Dd->SetMaximum(1.10*Histo_2Dd->GetMaximum());
		Histo_2Dd->SetMinimum(0.90*Histo_2Dd->GetMinimum());

//We want to see the fruit of our labour--- three images
	MyCanvas->cd(4);
	gPad->SetRightMargin(0.2);
	gPad->SetLeftMargin(0.15);
	gPad->SetTopMargin(0.15);
	gPad->SetBottomMargin(0.15);
	Histo_2Db->SetContour(100);
	Histo_2Db->SetStats(kFALSE);   // suppress statistics box
	Histo_2Db->Draw("colz");
	
	MyCanvas->cd(5);
	gPad->SetRightMargin(0.2);
	gPad->SetLeftMargin(0.15);
	gPad->SetTopMargin(0.15);
	gPad->SetBottomMargin(0.15);
	Histo_2Dc->SetContour(100);
	Histo_2Dc->SetStats(kFALSE);   // suppress statistics box
	Histo_2Dc->Draw("colz");
	
	MyCanvas->cd(6);
	gPad->SetRightMargin(0.2);
	gPad->SetLeftMargin(0.15);
	gPad->SetTopMargin(0.15);
	gPad->SetBottomMargin(0.15);
	Histo_2Dd->SetContour(100);
	Histo_2Dd->SetStats(kFALSE);   // suppress statistics box
	Histo_2Dd->Draw("colz");
	MyCanvas->Update();

	// write ROOT output file 
	// write everything to file, overwriting previuos runs
	hfile->Write();
	hfile->Close();
//	TBrowser *T = new TBrowser();
	MyCanvas->Update();

	
//     Clean up at the end of the program
	MyCanvas->WaitPrimitive();   
	// Waits untill you click the middle mouse button in the graphics window 
//	Halt program for debugging
//	char xxx;
//	cout << " Type a letter     ..... ";
//	cin >> xxx;



	cout << "\n\n Bye there \n";

	return 0;//finito!
}

/* Definition of functions */
void myswap2(unsigned short *datum) {
	unsigned char *d, t ;
	d=(unsigned char *)datum	;
	t = d[0];
	d[0] = d[1];
	d[1] = t;
} 
void edf2text(FILE *fp) {
	char myc, mys[256], mys2[256];
	int need_swap = 0, bytes_read, i;
	unsigned char *data_array;			// For the data in the EDF file
	unsigned short data;

	while ((fgets(mys,100,fp))!=NULL && ftell(fp) <1024) {
		//fprintf(stderr,"%s", mys);
		if (sscanf(mys,"ByteOrder = %s",mys2) ) {
			if (strncmp(mys2,"LowByteFirst",12)==0 )
				need_swap = 0;
			else
				need_swap = 1;
		}
		if (sscanf(mys,"Size = %i",&rsize) ) {
			//fprintf(stderr,"got size=%d\n",rsize);
			// not used yet!!!
			break;
		}
	}
	//fprintf(stderr,"Need swap %s (%s)\n",need_swap?"YES":"NO", mys2);

	fseek(fp,0,0);
	while ((myc = fgetc(fp)) !=EOF && myc!='}')
		;
	fgetc(fp); // there is a \n afterwards

	if ((data_array =(unsigned char *) malloc (rsize)) == NULL) {
		printf("Ooops .. No memory !");
		exit(1);
	}

		if (need_swap) {
			bytes_read = fread(data_array, 1, rsize, fp);
			//fprintf(stderr,"bytes read = %d\n", bytes_read);
			for (i=0; i<rsize ; i+=2) {
				data = (unsigned short) data_array[i];
				if (need_swap)
					myswap2(&data);
					data_image[i/2]=data;
				}
			}
		else {
			bytes_read = fread(data_image, 2, rsize/2, fp);
			}

	
	free(data_array);   // Free memory allocated for EDF file
}

// Constant background function with shifted origin
Double_t background(Double_t *x, Double_t *par) {
   return par[0] ;
}


//Gaussian Peak function
Double_t GaussianPeak(Double_t *x, Double_t *par) {
  return par[0]*exp( -((x[0]-par[1])*(x[0]-par[1]))/ (2*par[2]*par[2]));
}

// Sum of background and peak function
Double_t fitFunction(Double_t *x, Double_t *par) {
  return background(x,par) + GaussianPeak(x,&par[1]);
}

// Solve the Endian problem
unsigned short ByteSwap4 (unsigned short nValue)
{
   return (((nValue>> 8)) | (nValue << 8));

}