#include "stdio.h"
#include "TROOT.h"
//#include "TH1I.h"
#include "TH1F.h"
//#include "TH1D.h"
//#include "TH2I.h"
#include "TH2F.h"
#include "TF1.h"
//#include "TFile.h"
//#include "TBenchmark.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLine.h"

// ------------------------------------------------------------------------------

static int n_chains = 8;
static int n_chips = 6;
static int n_chips_a[5];
static int n_ch = 18;

int ReadIni(char *inifile)
{
	FILE *ifp;
	int params[5], chains = 0, chips = 0;
	for(int i=0; i<5; i++) params[i] = 0;
	
	if( (ifp=fopen(inifile, "r")) == NULL) return 0;
	
	fscanf(ifp, "nchip %d %d %d %d %d", &params[0], &params[1], &params[2], &params[3], &params[4]);
	
	for(int i=0; i<5; i++) {
		n_chips_a[i] = params[i];
		if(params[i] != 0) {
			chains++;
			chips = params[i];
		}	
	}
	
	fclose(ifp);
	
	n_chains = chains;
	n_chips = chips;
	
	//printf("n_chains = %d, n_chips = %d, n_ch = %d, n_chips_a = (%d, %d, %d, %d, %d)\n", 
	//n_chains, n_chips, n_ch, n_chips_a[0], n_chips_a[1], n_chips_a[2], n_chips_a[3], n_chips_a[4]);
	
	return 1;
}

#define HAPDCH n_chains*n_chips*n_ch

// ------------------------------------------------------------------------------

#define READBUFFERLENGTH 50000

// data format
#define MAXDATA 200

#define RUNREC_ID 1
#define ENDREC_ID 2
#define POSREC_ID 3
#define EVTREC_ID 4
#define THRREC_ID 5

typedef struct {
	UInt_t id,len;
	UInt_t time;
	UInt_t fver;
	UInt_t nev,nch,thr,xy;
	Int_t nx,x0,dx,ny,y0,dy;
	Int_t nth,th0,dth;
} RUNREC;
RUNREC *runrec;
RUNREC run;

typedef struct {
	UInt_t id,len;
	UInt_t time;
} ENDREC;
ENDREC *endrec;

typedef struct {
	UInt_t id,len;
	UInt_t time;
	Int_t ix,x,xset,iy,y,yset;
} POSREC;
POSREC *posrec;
POSREC pos;

typedef struct {
	UInt_t id,len;
	UInt_t time;
	UInt_t nev;
	UInt_t data[MAXDATA];
} EVTREC;
EVTREC *evtrec;

typedef struct {
	UInt_t id,len;
	UInt_t time;
	Int_t ith,th;
} THRREC;
THRREC *thrrec;
THRREC thr;

// ------------------------------------------------------------------------------
// Swap 4 bytes of one word
UInt_t Swap4(UInt_t *i)
{
	//printf("IN : i = %u, bytes = %u %u %u %u\n", *i, (*i>>24)&0xff, (*i>>16)&0xff, (*i>>8)&0xff, (*i)&0xff);

	*i = ((*i & 0xff000000) >> 24) +
	     ((*i & 0x00ff0000) >> 8) +
	     ((*i & 0x0000ff00) << 8) +
	     ((*i & 0x000000ff) << 24);

	//printf("OUT: i = %u, bytes = %u %u %u %u\n", *i, (*i>>24)&0xff, (*i>>16)&0xff, (*i>>8)&0xff, (*i)&0xff);
	
	return *i;
}

// ------------------------------------------------------------------------------
// shows treshold scan results and creates param file to align sigma_mult 
// position of noise
int tscan(char* datfile="tscan008", double x4s_avg = 125.0, char *csvfile="param_4s_000", 
			int sigma_mult=4.0, int dbg = 0)
{
	char sbuff[256];
	/*
	char *inifile="tscan008";
	char fullininame[256];
	sprintf(fullininame, "%s.ini", inifile);
	if(!ReadIni(fullininame)) {
		printf("Error opening ini file!\n");
		return 0;
	} else printf("	Opened ini file: %s\n", fullininame);
	*/
	//histograms
	char hname[256];
	double hthr_low = 0.0, hthr_hi = 0.0;
	TH2F *hthr = (TH2F*)gROOT->FindObject("hthr");	
	TH1F *hthr_py[HAPDCH], *hthr_pys = (TH1F*)gROOT->FindObject("hthr_pys");

	//data buffer
	char readbuf[READBUFFERLENGTH];
	runrec = (RUNREC *) readbuf;
	endrec = (ENDREC *) readbuf;
	posrec = (POSREC *) readbuf;
	evtrec = (EVTREC *) readbuf;
	thrrec = (THRREC *) readbuf;
	
	//data file
	FILE *fp;
	char fulldatname[256];
	
	// opens data file
	sprintf(fulldatname, "data/%s.dat", datfile);
	if( (fp=fopen(fulldatname, "rb")) == NULL ) {
		printf("	Cannot open data file: %s !!!\n", fulldatname);
		return -1;
	} else printf("	Opened data file: %s\n", fulldatname);
		
	UInt_t rec_id, rec_len;	
	unsigned ulsize = sizeof(UInt_t);
	int need_swap = 1;
	//read first word and detirmine, if swap is necessary
	fread(&rec_id,  1, ulsize, fp);
	fseek(fp, -1*(int)(ulsize), SEEK_CUR);
	//expecting RUNREC_ID
	if( ((rec_id>>24)&0xff) == RUNREC_ID ) need_swap = 1;
	else if( (rec_id&0xff) == RUNREC_ID ) need_swap = 0;
	else printf("Warrning, unexpected format!\n");	
	if(dbg) printf("	need_swap = %d\n", need_swap);

// -----------------------------------------------
	
	int iw, io, ih, ii;
		
	// loop trough records
	int ceve=0;
	while(1) {
		if( fread(&rec_id,  1, ulsize, fp) < ulsize ) break; if(need_swap) Swap4(&rec_id);
		if( fread(&rec_len, 1, ulsize, fp) < ulsize ) break; if(need_swap) Swap4(&rec_len);
		fseek(fp, -2*(int)(ulsize), SEEK_CUR);
		
		if(dbg) printf("-----------------------------------------------\n");
		if(dbg) printf("[%d] rec_id = %d | rec_len = %d\n", ceve, rec_id, rec_len);
		
		switch(rec_id)
		{		
			case RUNREC_ID:
				if( fread(readbuf, 1, rec_len, fp) < rec_len ) break;
				if(need_swap) for(int i=0; i<(READBUFFERLENGTH - 4); i+=4) Swap4((UInt_t*)&readbuf[i]);
				
				run = *runrec;
				
				if(dbg) {
					printf("RUNREC_ID\n");
					printf("id = %d, len = %d, time = %d\n", run.id, run.len, run.time);
					printf("nev = %d, nch = %d, thr = %d\n", run.nev, run.nch, run.thr);
					printf("nx = %d, x0 = %d, dx = %d\n", run.nx, run.x0, run.dx);
					printf("ny = %d, y0 = %d, dy = %d\n", run.ny, run.y0, run.dy);
					printf("nth = %d, th0 = %d, dth = %d\n", run.nth , run.th0, run.dth);
					
					break;
				}
				
				// create histograms
				hthr_low = run.th0 - 0.5*run.dth;
				hthr_hi = hthr_low + run.dth*run.nth;
				
				sprintf(hname, "hthr");
				hthr = (TH2F*)gROOT->FindObject(hname);
				if(hthr) delete hthr;
				hthr = new TH2F(hname, hname, HAPDCH, -0.5, HAPDCH - 0.5, run.nth, hthr_low, hthr_hi);
				
				for(int i=0; i<HAPDCH; i++) {
					sprintf(hname, "TH CH.%d", i);
					hthr_py[i] = (TH1F*)gROOT->FindObject(hname);
					if(hthr_py[i]) delete hthr_py[i];
					hthr_py[i] = new TH1F(hname, hname, run.nth, hthr_low, hthr_hi);
				}
				
				sprintf(hname, "TH PY");
				hthr_pys = (TH1F*)gROOT->FindObject(hname);
				if(hthr_pys) delete hthr_pys;
				hthr_pys = new TH1F(hname, hname, run.nth, hthr_low, hthr_hi);
				break;
			
			case POSREC_ID:
				if( fread(readbuf, 1, rec_len, fp) < rec_len ) break;
				if(need_swap) for(int i=0; i<(READBUFFERLENGTH - 4); i+=4) Swap4((UInt_t*)&readbuf[i]);
				
				pos = *posrec;
				
				if(dbg) {
					printf("POSREC_ID\n");
					printf("id = %d, len = %d, time = %d\n", posrec->id, posrec->len, posrec->time);
					printf("ix = %d, x = %d, xset = %d\n", posrec->ix, posrec->x, posrec->xset);
					printf("iy = %d, y = %d, yset = %d\n", posrec->iy, posrec->y, posrec->yset);
				} else  printf("	%d, %d\n", posrec->xset, posrec->yset);
				break;
				
			case THRREC_ID:
				if( fread(readbuf, 1, rec_len, fp) < rec_len ) break;
				if(need_swap) for(int i=0; i<(READBUFFERLENGTH - 4); i+=4) Swap4((UInt_t*)&readbuf[i]);
				
				thr = *thrrec;
				
				if(dbg) {
					printf("THRREC_ID\n");	
					printf("id = %d, len = %d, time = %d\n", thrrec->id, thrrec->len, thrrec->time);
					printf("ith = %d, th = %d\n", thrrec->ith, thrrec->th);	
				}					
				break;
				
			case EVTREC_ID:
				if( fread(readbuf, 1, rec_len, fp) < rec_len ) break;
				if(need_swap) for(int i=0; i<(READBUFFERLENGTH - 4); i+=4) Swap4((UInt_t*)&readbuf[i]);
				
				if(dbg) {
					printf("EVTREC_ID\n");
					printf("id = %d, len = %d, nev = %d\n", evtrec->id, evtrec->len, evtrec->nev);
					for(int datai = 0; datai < MAXDATA; datai++) printf("%u ", evtrec->data[datai]);
					printf("\n");
					break;
				}
// events ------------------
				// fill histograms
				for(int j=0; j<n_chains; j++) {
					for(int i=0; i<(n_chips*n_ch); i++) {
						ii= j*(n_chips*n_ch) + i;
						
						iw = j*(int(n_ch*n_chips/8.0) + 1) + i/8;
						io = i%8;
						ih = (evtrec->data[iw]>>(4*io))&0xf;
					
						if(ih != 0) {
							hthr->Fill((double)ii, (double)(thr.th));
							hthr_py[ii]->Fill((double)(thr.th));
						}
					}
				}

// events ------------------				
				break;
			
			case ENDREC_ID:
				if( fread(readbuf, 1, rec_len, fp) < rec_len ) break;
				if(need_swap) for(int i=0; i<(READBUFFERLENGTH - 4); i+=4) Swap4((UInt_t*)&readbuf[i]);
				
				if(dbg) {
					printf("ENDREC_ID\n");
					printf("id = %d, len = %d, time = %d\n", endrec->id, endrec->len, endrec->time);
				} else printf("	ENDREC\n");
				break;
			
			default:
				printf("switch(rec_id): default !!!\n");
				return 0;
				break;
		}
		
		ceve++;
		if( (ceve%50000) == 0) printf("	Current event = %d\n", ceve);
		if(dbg) if( ceve>6 ) break;
		if(feof(fp)) break;
	}
	
	fclose(fp);
// -----------------------------------------------

	if(dbg) return 1;

///*
	// plot data
	
	gStyle->SetStatBorderSize(1);
	gStyle->SetFrameBorderMode(0);
	gStyle->SetCanvasBorderMode(0);
	gStyle->SetPadBorderMode(0);
	gStyle->SetPadColor(0);
	gStyle->SetCanvasColor(0);
	gStyle->SetHistFillColor(18);
	gStyle->SetTitleFillColor(0);
	gStyle->SetStatColor(0);
	gStyle->SetPalette(1, 0);
	
	gStyle->SetPadTopMargin(0.13);
	gStyle->SetPadBottomMargin(0.13);
	gStyle->SetPadLeftMargin(0.13);
	gStyle->SetPadRightMargin(0.15);
	
	gStyle->SetTitleOffset(1.3, "y");
	
	//gStyle->SetOptFit(0);
	//gStyle->SetOptStat(0);
	//gStyle->SetOptTitle(0);
	
	
	TCanvas *c0 = new TCanvas("c0", "Treshold", 0, 0, 1200, 700);
	hthr->Draw("COLZ");
	
	TLine *line_chain = new TLine();
	line_chain->SetLineWidth(2);
	line_chain->SetLineColor(2);
	
	TLine *line_chip = new TLine();
	line_chip->SetLineWidth(1);
	line_chip->SetLineColor(1);
	
	for(int i=0;i<HAPDCH;i++) {
		if( !((i+1)%108) ) line_chain->DrawLine(i+1, 280, i+1, 340);
		if( !((i+1)%18) ) line_chip->DrawLine(i+1, 280, i+1, 320);
	}
	
	line_chip->DrawLine(108, -200, 108, 330);
	
	sprintf(sbuff, "eps/thr2d_%s.gif", datfile);
	c0->SaveAs(sbuff);
	
	// plot treshold scan for each channel 
	// and calculate shifts to align sigma_mult positions
	TCanvas *c1 = new TCanvas("c1", "Treshold Projections", 100, 0, 700, 1000);
	c1->Divide(3, 6);
	int numpads = 18;
	int numeps = 0;
	
	TF1 *thr_g = new TF1("thr_g", "gaus", hthr_low, hthr_hi);
	double thr_par[HAPDCH][3];
	double x4sigma[HAPDCH];
	
	int dead_ch[HAPDCH], hot_ch[HAPDCH];
	int dead_i = 0, hot_i = 0;
	for(int i=0; i<HAPDCH; i++) {	
		dead_ch[i] = -1;
		hot_ch[i] = -1;
	}
	
	// find very noisy channels and save them to file
	FILE *noisyfp;
	sprintf(sbuff, "params/%s.hot", datfile);
	if( (noisyfp = fopen(sbuff, "w")) == NULL) printf("Error creating noise_ch file!\n");
	TH1F *hnoisy = new TH1F("hnoisy", "hnoisy", 250, 0, 250);
	
	for(int i=0; i<HAPDCH; i++) {	
		c1->cd((i%numpads)+1);
		(c1->cd((i%numpads)+1))->SetLogy();
		
		if(hthr_py[i]->GetMaximum() > 10) {
			hthr_py[i]->Fit(thr_g, "Q");
			thr_g->GetParameters(thr_par[i]);
			if(thr_par[i][2] > 60.0)
				hot_ch[hot_i++] = i;
			//else
				//hthr_pys->Add((TH1*)hthr_py[i]);
		} else {
			(c1->cd((i%numpads)+1))->Clear();
			thr_par[i][0] = 0.0;
			thr_par[i][1] = 0.0;
			thr_par[i][2] = 0.0;
			dead_ch[dead_i++] = i;
		}
		x4sigma[i] = thr_par[i][1] + sigma_mult*thr_par[i][2];
		
		//save canvas to eps
		if( !((i+1)%numpads) ) {
			sprintf(sbuff, "eps/thr_%s_%2d.gif", datfile, numeps++);
			c1->SaveAs(sbuff);
		}
		
		double noise = hthr_py[i]->Integral( (hthr_py[i]->GetXaxis())->FindBin(190), (hthr_py[i]->GetXaxis())->FindBin(300) );
		if(noise > 20) fprintf(noisyfp, "%d\n", i+1);
		else hthr_pys->Add((TH1*)hthr_py[i]);
		hnoisy->Fill( noise );
	}
	
	fclose(noisyfp);
	TCanvas *c1c = new TCanvas("c1c", "Over noise", 100, 300, 700, 1000);
	(c1c->cd(0))->SetLogy();
	hnoisy->Draw();
	
	// save hot and dead channel info
	//FILE *hotfp;
	//sprintf(sbuff, "params/%s.hot", datfile);
	//if( (hotfp = fopen(sbuff, "w")) == NULL) printf("Error creating hot_ch file!\n");
	//else for(int i=0; i<hot_i; i++) fprintf(hotfp, "%d\n", hot_ch[i]+1);
	//fclose(hotfp);
	sprintf(sbuff, "params/%s.dead", datfile);
	if( (hotfp = fopen(sbuff, "w")) == NULL) printf("Error creating dead_ch file!\n");
	else for(int i=0; i<dead_i; i++) fprintf(hotfp, "%d\n", dead_ch[i]+1);
	fclose(hotfp);
	/*
	TCanvas *c0a = new TCanvas("c0a", "Treshold Proj", 0, 0, 800, 1000);
	c0a->Divide(1,2);
	c0a->cd(1);
	(hthr->ProjectionY("hthr_py"))->Draw();
	c0a->cd(2); (c0a->cd(2))->SetLogy();
	(hthr->ProjectionY("hthr_py"))->Draw();
	
	sprintf(sbuff, "eps/thr_proj_%s.gif", datfile);
	c0a->SaveAs(sbuff);
	*/
	// plot th proj. only for good channels
	TCanvas *c0b = new TCanvas("c0b", "Treshold Proj", 300, 200, 800, 1000);
	c0b->Divide(1,2);
	c0b->cd(1);
	TH1F *hprojfull = (TH1F*)(hthr->ProjectionY("hthr_py"));
	(hprojfull->GetYaxis())->UnZoom();
	hprojfull->Draw();
	hthr_pys->SetFillColor(40);
	hthr_pys->Draw("SAME");
	c0b->cd(2); (c0b->cd(2))->SetLogy();
	(hprojfull->GetYaxis())->SetRangeUser(1.0, hprojfull->GetMaximum() + hprojfull->GetMaximum()/12.0);
	hprojfull->Draw();
	hthr_pys->SetFillColor(40);
	hthr_pys->Draw("SAME");
	
	sprintf(sbuff, "eps/thr_proj_s_%s.gif", datfile);
	c0b->SaveAs(sbuff);
	
	
	// differences of sigma_mult position from average for each channel
	double x4s_delta[HAPDCH];
	// difference in bits for setting the ASIC
	int bitm[HAPDCH];
	TH1F *hbitm = new TH1F("hbitm", "hbitm", 32, -0.5, 31.5);
	// koeficient to convert mV values to bits
	double x2b_k = 9.0;
	
	for(int i=0; i<HAPDCH; i++) {
		x4s_delta[i] = x4sigma[i] - x4s_avg;
		bitm[i] = int( 8.5 - x4s_delta[i]/x2b_k );
		// make shure very noisy channels get in the right range
		if( bitm[i] > 16 ) bitm[i] = 16;
		if( (-3 < bitm[i])&&(bitm[i] < 0) ) bitm[i] = 0;
		if( bitm[i] < 0 ) bitm[i] += 2;
		if( bitm[i] < -15 ) bitm[i] = -15;
		
		bitm[i] = 32 - bitm[i];
		
		if( bitm[i] > 31 ) bitm[i] -= 32;
		
		hbitm->Fill(bitm[i]);

		//printf("[%3d]%7.1lf,%6.1lf,%6.1lf |%6.1lf,", i, thr_par[i][0], thr_par[i][1], thr_par[i][2], x4sigma[i]);
		//printf("%7.1lf,%3d\n", x4s_delta[i], bitm[i]);
	}
		
	// histogram bitm
	TCanvas *c2 = new TCanvas("c2", "Offset bits", 200, 0, 700, 500);
	hbitm->Draw();
	sprintf(sbuff, "eps/bitm_%s.gif", datfile);
	c2->SaveAs(sbuff);
	
		
	// create treshold map file
	FILE *fpo;
	char fullcsvname[256];
	sprintf(fullcsvname, "params/%s.csv", csvfile);
	if( (fpo=fopen(fullcsvname, "w")) == NULL ) {
		printf("	Cannot create param file !!!\n");
		return 0;
	} else printf("	Created param file %s\n", fullcsvname);
	
	fprintf(fpo, "#\n");       
    fprintf(fpo, "# pts PTS address\n");     
    fprintf(fpo, "pts 0 0x50000000\n");     
    fprintf(fpo, "pts 1 0x51000000\n");       
    fprintf(fpo, "#\n");        
    fprintf(fpo, "# configure PTS filename\n");     
    fprintf(fpo, "configure 0 ../bit/ptshapd.bit\n");    
    fprintf(fpo, "configure 1 ../bit/ptshapd.bit\n");        
    fprintf(fpo, "#\n");        
    fprintf(fpo, "# nchip PTS NCHIP0 NCHIP1 NCHIP2 NCHIP3\n");  
    fprintf(fpo, "# define number of chips per chain\n");  
    fprintf(fpo, "nchip 0 6 6 6 6\n");
    fprintf(fpo, "nchip 1 6 6 6 6\n");
	//fprintf(fpo, "nchip ");
    //for(int i=0;i<5;i++) fprintf(fpo, "%d ", n_chips_a[i]);
	//fprintf(fpo, "\n");

    fprintf(fpo, "#  \n");      
    fprintf(fpo, "# csr1 PTS DATA \n");    
    fprintf(fpo, "csr1 0 0x2f0000\n");    
	fprintf(fpo, "csr1 1 0x2f0000\n");   
    fprintf(fpo, "#\n");        
    fprintf(fpo, "# halfcycle PTS VALUE\n");     
    fprintf(fpo, "# trgdelay PTS VALUE\n");     
    fprintf(fpo, "# optionally, halfcycle and trgdelay can be changed\n"); 
    fprintf(fpo, "halfcycle 0 6\n"); 
    fprintf(fpo, "trgdelay 0 12 # 13 = 0xd\n");      
    fprintf(fpo, "halfcycle 1 6\n");  
    fprintf(fpo, "trgdelay 1 12 # 13 = 0xd\n");  
    fprintf(fpo, "#   \n");     
    fprintf(fpo, "# monitor PTS VALUE\n");     
    fprintf(fpo, "monitor 0 0\n");  
    fprintf(fpo, "monitor 1 0\n");       
    fprintf(fpo, "#\n");        
    fprintf(fpo, "# param PTS CHAIN CHIP CH OFFSET PAR\n"); 
    fprintf(fpo, "#  \n");

	for(int j=0;j<n_chains;j++)
		for(int i=0; i<(n_chips*n_ch); i++) {
			ii= j*(n_chips*n_ch) + i;
			
			if(!(i%18)) fprintf(fpo, "# \n");
			
			fprintf(fpo, "param %d %d %d %d %d 2\n", 
			j/(n_chains/2),
			j%(n_chains/2),
			(i/18)%n_chips,
			i%18,
			bitm[ii]);	
		}
	
	
    fprintf(fpo, "#  \n");      
    fprintf(fpo, "# init PTS \n");     
    fprintf(fpo, "init 0\n"); 
    fprintf(fpo, "init 1\n");   
       
    fprintf(fpo, "# \n");
    fprintf(fpo, "# load_param  \n");     
    fprintf(fpo, "load_param\n");    
    fprintf(fpo, "# EOF \n");
	
	fclose(fpo);
		
// -----------------------------------------------			
	
	return 1;
}