////////////////////////////////////////////////////////////////////////////
//	AZMI		13/05/05		Berlin
// 	t0 calibration 
//	
//	
//	o Fit lower end of each tdc spectra (exclude noisy channels)
//		o with a half/half gaussian func
//		  (function in file twogausfit.cxx)
//		o common mean but varying left and right sigma
//		o define tzero as steepest point on slope 
//		  (time when exponential in gaussian == 1/2)
//	o Plot tdc spectra with fit
//	o collect array of tzero values for each histogram
//	o plot histogram of tzero values
//	o plot histogram of deviation of mean tzeros from avrg tzero
//
//	to execute: Drkam_monitor("filename",int bins,double xlow,double xmax)
//	drift times from 0 - 1100 so similar range should be chosen for histogarm
//	should work for any binning
//
//	til line <226> declarations and functions
//	line <230> main program starts
//	lines <230> - <361> reading file and filling histogram
//	lines <363> - <471> FITTING histograms and extracting t0
//	lines <473> - <528> summary plots for t0
//
////////////////////////////////////////////////////////////////////////////

//new things learnt:
// o pass by reference (passing array into fit_tdc to extract fit parameters)
//	To pass the array; myarray[3]
// 		in main prog: func(otherarguments, myarray, 3);
//		the function declaration: func(types arguments, int array[], int size)
// 	In C though, u just pass the array as myarray[3]
//
// o TSpectrum to find position of peaks on a spectrum 
//

#include "TH1F.h"
#include "TCanvas.h"
#include "TF1.h"
#include <TSpectrum.h>
#include <TMath.h>
#include <TStyle.h>
#include <TPaveStats.h>
#include <TPaveLabel.h>
#include <TPaveText.h>
#include <iostream.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "twogausfit.cxx"

// prototypes
//
int chtoi16(char ch);
int read_event(FILE *inp, struct time_hit *timehit);
double mean(double matrix[], int size);
int rounding(double n);
int firstentry(TH1F *h1);

// define structure
//
struct time_hit 
{
    int numtdc;
    char month[4],year[5],day[3],time[10];
    int count[16];
    unsigned int hit[16][16];
};

// global variable
//
struct time_hit timehit0[4];

// convertor: 4 TDC (4x16 channels) --> tube numbers (72 tubes)
// 8 tubes are nor read out (see the cable scheme)
// first tube number is 0
int tdc2tube[4][16]={62,55,71,63,48,64,56,65,49,57,50,66,58,51,67,59,
                     25,33,26,42,34,27,43,35,52,68,60,69,53,61,54,70,
                     28,44,36,45,29,37,30,46,38,31,47,39,24,40,32,41,
                      0, 8, 1, 9, 2,10, 3,11, 4,12, 5,13, 6,14, 7,15};

// mean() calculates average (the mean) of all non zero elements in an array
// Purpose: calculate average of tzeroes (excluding non zero entries which are for noisy channels)
//		o used to find deviation of all tzeroes
//
double mean(double matrix[], int size)
{
    double avrg=0;
    int n=0;
    double sum=0;
	    
    for (int i =0; i<size; i++)
    {
	    if (matrix[i] != 0.0)
	    {
		    sum=sum+matrix[i];
		    n++;
	    }
	    
    }
    avrg=sum/n;
    cout<<endl<<"func _ avrg tzero : "<<avrg<<endl;
    return avrg;
}


// rounding() rounds up or down a double to an integer
// Purpose: round up tzeroes
//	o they have an integer value error (so higher precision makes no sense)
//
int rounding(double n)
{
	double diff = n - int(n);
	
	if (diff < 0.5) { return abs(floor(n)); }
	else { return abs(ceil(n));}
	
}

// firstentry() returns the bin number before the first non-zero entry of a histogram
// purpose: to effectively estimate a lower limit for fitting
//	o some tdc spectra start well after x=0
//	o so setting a number (e.g. x=0) as lower limit gives funny fit
//
int firstentry(TH1F *h1)
{
	int entries;			//no of oentries in bin
	int nbins=h1->GetNbinsX();	//number of bins
	int x; 				//x co-ordinate of bin
	
	for (int i=1; i<=nbins; i++)	//binning starts at nbin=1		
	{
		entries = h1->GetBinContent(i+1);	// entries in current bin
		if (entries !=0)			// for first non-zero bin
		{
			x= h1 -> GetBinCenter(i);	//returns x co-ord of last empty bin
			return (x);			
		}
	}
}


// thepeak() compares an array of real numbers to a given real and returns the closest element of the array 
// Purpose: return the x position of the peak closest to the ...
//		o modal drift time ...? 
//		o observed peak position of 200 ..?
//
float thepeak(float array[], int size, float n)
// array[] is array holding numbers, size = size of array, n = ideal element (for comparison)
{
	int min=0;
	float diff[size];
	
	for (int i =0; i<size; i++)
	{
		diff[i] = fabs(n - array[i]);	// find difference between element and ideal number
		if ( diff[i] < diff[min] )	// sort for smallest elemnt (ie closest number)
			min = i;
	}
	
	return array[min] ;			// return value closest to ideal number
	
}

// chtoi16() converts hexadecimal notation into decimal numbers 1 to 16
//
int chtoi16(char ch)
{  
	if(ch>='0' && ch<='9') return(ch-'0');
	if(ch>='A' && ch<='F') return(ch-'A'+10);
	if(ch>='a' && ch<='f') return(ch-'a'+10);
	return(-1);
}


// read_event() reads the hits of one event of one TDC
// 
int read_event(FILE *inp, struct time_hit *timehit)
{   
    int i,ii,nch;
    unsigned int hit;
    if(inp==0){ printf("Something wrong with inp file\n");  return(3); }
    char tim1[10],tim2[10],tim3[10],tim4[10],tim5[10];
    char chtmp,ch[11];

    for(i=0;i<16;i++) 
    { 
	    timehit->count[i]=0;	// initialisation
    }
    // feof: Check if End Of File has been reached
    if(feof(inp)) { printf("EOF achieved\n"); return(1); }
    // fgetc: Get next character from a stream
    chtmp=fgetc(inp);
    if(feof(inp)) { printf("EOF achieved\n"); return(1); }
    // fscanf: 	Read formatted data from a stream
    fscanf(inp,"%s%s%s%s%s%c",tim1,tim2,tim3,tim4,tim5,&chtmp);// Read date
    // strcpy: 	Copy string
    strcpy(timehit->day,tim3);
    strcpy(timehit->month,tim2);
    //cout<<"month="<<timehit->month<<endl;
    strcpy(timehit->time,tim4);
    strcpy(timehit->year,tim5);
    for(i=0;i<11;i++) 			// why 11 here and 10 on line 82?
    {					// FIRST LINE
    	ch[i]=fgetc(inp);		// -> HEADER
    }					//
    ii=chtoi16(ch[2]);	
    timehit->numtdc=ii;			// Set TDC number
    //cout<<ii<<"  ";
    do 		// Read all data 
    {
        for(i=0;i<10;i++) 		// Read one line of data 
	{
		ch[i]=fgetc(inp); if(ch[i]=='\n') return(0);
        }
        chtmp=fgetc(inp);
        unsigned long tmp=1;
        for(i=3,hit=0;i>=0;i--)		// Read time beginning at the end 
	{				// and convert it into decimal time
            hit+=chtoi16(ch[i+6])*tmp;	// hit is "time of hit"
            tmp*=16;
        }
        nch=chtoi16(ch[3]);		// channel number
        i=(timehit->count[nch])++;	// counts in channel +1
        timehit->hit[nch][i]=hit;	
    }
    while(1);
    return(0);
}


// main program
//
int Drkam_monitor(char *name, int bins, double xlow, double xmax)
{
    char naam[4], title[100];
    
    /***************** FOR A 4 CHANNEL PER CANVAS OUTPUT*******************/               	            		
    
    TCanvas *can[16];
    int j;
    {
	    for (int i=0; i<16; i++)
	    //int i = 15; //for one canvas (last 4 histo)
	    {
		    sprintf(naam, "can %d", i);
		    j=i/4;
		    sprintf(title, "Canvas for TDC %d", j);
		    can[i]= new TCanvas(naam,title,1200,900);
		    //can[i]= new TCanvas(naam,title,1);
		    can[i]->Divide(1,4);
	    }
    }
    /*------------------------------------------------------------------------*/
    
    
    /****************** HISTOGRAM FOR EACH CHANNEL OF EACH TDC ****************/
    
    //for small times
    TH1F *TDCsmalltime[64];
    int hist =0;
    for (int i=0; i<4; i++)
    {
	    for (int j=0; j<16; j++)
	    {
		    sprintf(naam, "%d", hist);
		    sprintf(title, "TDC spectra for TDC %d - channel %d ... tube # %d;Time /ns;Total count rate",i,j, tdc2tube[i][j]);
		    TDCsmalltime[hist] = new TH1F(naam,title,bins,xlow,xmax);
		    hist = hist+1;
	    }
    }
    
    // number of nano seconds per bin
    // the range divided by number of bins
    double binsize = (xmax-xlow)/bins;
    /***************************************************************************/
    
    /************Include Overflow and Underflow on canvas **************/
    TStyle *mystyle = new TStyle();
    mystyle->Copy(*gStyle);
    mystyle->SetOptStat(111111);
    //mystyle->SetOptFit(1111);
    mystyle->cd();
    
    
    //main program
    // Opening file
    FILE *inp;
    printf("%s\n",name);
    // fopen: 	Open a file	r: read, t: text mode
    inp=fopen(name,"rt");
    if(!inp)			// Check File pointer
    {
        printf("Error opening file!\n");
        return(1);
    }
    
    int tube = 0; 	//tube number
    int err;
    int event=0;
    int hittime=0;    	//time from tdc (taken in commonstop mode)
    int delay=1080;  	//delay btwn triggersignal and common stop. what is true delay? 1080 gives better histogram than 1000
    int signaltime=0; 	//corrected drift time
    int allentries=0; 	//all smalltime
    int tmp0; 		//tmp variable to sort for smallest time
    int entry[64];	//number of entries in each tdc spectra
    
    for (int i=0; i<64; i++) {entry[i]=0;}	//initialise entry array
    
    do 	// run through all events
    {  	
	//event++;
	//cout<<"event: "<<event<<endl;
        err=read_event(inp,&(timehit0[0]));		// read TDC 0
	if(!err) err=read_event(inp,&(timehit0[1]));  	// read TDC 1
	if(!err) err=read_event(inp,&(timehit0[2]));  	// read TDC 2
	if(!err) err=read_event(inp,&(timehit0[3]));  	// read TDC 3
	//cout<<"err: "<< err<<endl;
	if(!err)
	{ 
		event++;
		//cout<<"event: "<<event<<endl;
		hist=0;
		for (int i=0; i<4; i++) //each TDC
		{
			for (int j=0; j<16; j++) //each channel
			{
				//convert channel to tube
				tube = tdc2tube[i][j];
				//initialise tmp.
				tmp0=0;
				if (timehit0[i].count[j]!=0)
				{tmp0=delay;}
				for (int k=0; k<timehit0[i].count[j]; k++) //till there are no more entries for that channel
				{
					hittime=timehit0[i].hit[j][k];	//time from tdc (common stop mode)
					//if hit within time window
					if (hittime<delay && hittime>0)
					{
						signaltime=delay-timehit0[i].hit[j][k]; //drift time
						//sort for smallest times
						if (tmp0>signaltime)
							{tmp0=signaltime;}
					}
				}
				
				//fill TDCsmalltime for smallest times
				if(tmp0<delay && tmp0>0)  //exclude times outside window
				{ 
					allentries++;
					TDCsmalltime[hist]->Fill(tmp0);
					entry[hist]++;
					//exclude results for noisytubes
					//if (tube!=71 && tube!=31 && tube!=30  && tube!=8 && tube!=7 && tube!=0)
						//{alltubes-> Fill(tmp0); entries++;}
				}
				hist = hist+1;
			}
			
		}
	}
    }
    while(!err);
    
    /******************* FOR A 4 CHANNEL PER CANVAS OUTPUT ********************/
    
    hist=0;
    
    //Fitting
    
    //the calculated tzero from parameters of each tdc spectra fit
    double tzero[64];	for (int i=0; i<64; i++) {tzero[i]=0;}
    //the error on tzero for each tdc spectra fit
    double errt0[64];	for (int i=0; i<64; i++) {errt0[i]=0;}
    
    double para[2]={0,0}; 	//to extract mean and sigmaleft parameters
    double error[3]={0,0,0};	//to extract mean and sigmaleft errors and their correlation error
    
    TPaveText *pt[64];		//Text labe for customised stats box
    
    for (int i=0; i<16; i++)  		//each canvas
    //int i = 15;  hist = 60;		//test program for just last few histograms
    	{
		for (int j=0; j<4; j++)  //each pad of the canvas
		{
		
		
		can[i]->cd(j+1);
		TDCsmalltime[hist] -> SetLineColor(kBlue);
		TDCsmalltime[hist] -> Draw();
		
		if (hist != 2 && hist !=41 && hist !=48 && hist !=49 && hist !=62 && hist!=1 && hist!=42)	//don't bother with fit for noisy tubes ... and 'funny' tubes 
		{	
			TDCsmalltime[hist] -> SetStats(kFALSE);
			
			//TSpectrum to find peaks
			TSpectrum *spectrum = new TSpectrum();				//Tspectrum with 2 peaks
			int peakwidth = 10;
			if (abs(binsize)==1) { peakwidth=50; }				// decrease resolution for peak finding
			
			int npeaks = spectrum -> Search(TDCsmalltime[hist],peakwidth,"new");		//#peaks on histogram with sigma peakwidth
			float *xpeaks = spectrum -> GetPositionX();			//get x positions of all peaks
			
			
			//estimate parameters
			//
			int maxbin = TDCsmalltime[hist]->GetMaximumBin(); 	//bin with most entries
			int maxentry = TDCsmalltime[hist] -> GetBinContent(maxbin); //number of entries at modal drift time
			
			float peakposition = thepeak(xpeaks, npeaks, 200);	//return peak posn nearest 200 (the 'by sight' peak)
			
			int fmean = peakposition;			//estimated mean = position of peak
			int fconstant = maxentry;		//estimated constant = maximum count
			int fsigmaleft = fmean / 3;		//estimated sigma left ... 1/3 of mean ~
			int fsigmaright = fmean * 1/2;		//estimated sigma right ... 1/2 mean ~
			int fuplim = fmean + 100;		//large uplim => bad fit
			int flowlim = firstentry(TDCsmalltime[hist]);		//estimated lowlim = bin before first entry
			
			
			
			//fit histogram with above estimated parameters
			//
			fit_tdc(TDCsmalltime[hist],flowlim,fuplim,fsigmaleft,fsigmaright,fmean,fconstant,para,error);
			
			
			//calc tzero from returned parameters
			//
			tzero[hist] = para[0] - para[1]*sqrt(log(4)); //log is the natural logarithm
			errt0[hist] = sqrt((pow(error[0],2))+(pow(error[1],2)*log(4))+(2*error[2]*sqrt(log(4)))); //error on tzero
			
			cout<<"****************************"<<endl;
			cout<<endl<<" tzero["<<hist<<"] = "<<rounding(tzero[hist])<<endl;
			cout<<"****************************"<<endl;
			
			//a nicety?
			double fit_mean = para[0];		double fit_errmean = error[0];
			double fit_sigma = para[1];		double fit_errsigma = error[1];
			double fit_tzero = tzero[hist];		double fit_errtzero = errt0[hist];
			
			/******************* Some labels for the canvas ********************/
			
			//for each channel: events, entries, mean, sigma, tzero
			char ptevent[20], ptentries[20], ptmean[40], ptsigma[40], pttzero[40];
			sprintf(ptevent, "events: %d ", event);
			sprintf(ptentries, "entries: %d ",entry[hist]);
			sprintf(ptmean, "Mean: %.0f #pm %.0f", fit_mean, fit_errmean);
			sprintf(ptsigma, "Sigma: %.0f #pm %.0f ", fit_sigma, fit_errsigma);
			sprintf(pttzero, "t_zero: %.0f #pm %.0f ", fit_tzero, fit_errtzero);
			
			pt[hist] = new TPaveText(0.80,0.40,0.90,0.90,"blNDC");
			pt[hist]->SetName("more stats");
			pt[hist]->SetBorderSize(2);
			pt[hist]->SetFillColor(19);
			pt[hist]->SetTextAlign(12);
			TText *text = pt[hist]->AddText("Some stats");
			pt[hist] -> AddLine(0,0.84,0,0.84);
			text = pt[hist]-> AddText(ptevent);
			text = pt[hist]-> AddText(ptentries);
			text = pt[hist]-> AddText(ptmean);
			text = pt[hist]-> AddText(ptsigma);
			text = pt[hist]-> AddText(pttzero);
			pt[hist]->Draw();
			
			//update labels additions to canvas
			gPad->Modified();
			gPad->cd();
			
			spectrum->~TSpectrum();
		}
		
		hist=hist+1;
		}
	}
	
	// integer values of tzeroes
	int inttzero[64]; 
	for (int i=0; i<64; i++)
	{
	    inttzero[i]=0;
	    inttzero[i]=rounding(tzero[i]);
	}
	
	
	
    
    /********************** tzero plots ***********************/
    
    TCanvas *Ctzero = new TCanvas("tzeroes", "Deviation of tzeroes", 1000,500);
    Ctzero->Divide(2,1);
    
    //********* histogram of tzeroes (the value of each tzero)
    sprintf(title, "Plot of calibrated tzeroes when tdc spectra have %.1fns perbin; (1ns per bin) .......... tzero / ns; count", binsize);
    TH1F *tzeroplot = new TH1F("tzero",title,61, 100.5, 160.5);
    
    for (int i=0; i<64 ;i++)
    {
	    if (tzero[i] != 0.0)	//noisy channels have tzero==0
	    {
		    tzeroplot->Fill(tzero[i]);
	    }
    }
    
    Ctzero->cd(1);
    tzeroplot->Draw();
    
    
    //********** histogram of deviation from mean (average - mean)
    TH1F *deviation = new TH1F("deviation","..and the deviation from the average tzero; (1ns per bin) ....... average - tzero /ns; count",61, -30.5,30.5);
    
    double average = mean(tzero,64);	//function mean, calculates mean of tzeroes
    double dev=0;
    for (int i=0; i<64; i++)
    {
	    if (tzero[i] != 0.0)	//noisy channels have tzero==0
	    {
	    dev=(int)average-inttzero[i];	//difference between average and t_zero
	    deviation -> Fill(dev);
	    }
    }
    Ctzero->cd(2);
    deviation->Draw();
    
    
    for (int i=0; i<64; i++)
    {
	    printf(" %4d   ",inttzero[i]);
	    if (i==15  || i==31 || i==47 || i==63)
			    printf("\n");
	    
    }
    
    //cout<<"tzero average: "<<average<<endl;
    //cout<<"binsize = "<<binsize<<endl;

   /*------------------------------------------------------------------------*/
    
    
    
    /************** Check prints *********************
    
    for (int n=0; n<npeaks; n++)
    {
	cout<<endl<<" peak "<<hist<<":"<<n<<" at x posn "<<xpeaks[n]<<endl;
    }
			
    cout<<"THE PEAK POSITION: "<<peakposition<<endl;
			
    
    //printf("\n Fit statistics returned from function to main via array: \n main_mean: %f \t main_errmean: %f \n main_sigma: %f \t main_errsigma: %f \t main_correrr: %f \n", para[0],error[0],para[1],error[1],error[2]);
    //printf("\n\n");
			
    //printf("\n Parameters to input into fit function \n mean= %d \t constant= %.2f \t sl= %d \t sr= %d \t uplim= %d \t lowlim=%d \n\n",mean,constant,sl, sr,uplim, lowlim);
    
    for (int i=0; i<64; i++)
    {
	    printf(" %4.1f   ",tzero[i]);
	    if (i==15  || i==31 || i==47 || i==63)
			    printf("\n");
	    
    }
    
    cout<<"Rounded up values of tzero:"<<endl;
	for (int i=0; i<64; i++)
	{
	    printf(" %3d   ",inttzero[i]);
	    if (i==15  || i==31 || i==47 || i==63)
			    printf("\n");
	    
	}
    
    /*
    cout<<"calculated tzero,(from mean&sigma extracted from fit_tdc): "<<tzero<<endl;
    printf("\n");
    cout<<"total number of small times: "<<allentries<<endl;
    cout<<"events: "<<event<<endl;
    cout<<"entries in histogram (smalltimes in window excluding noisy channels): "<<entries<<endl;
    //printf("Entries: %d \n",tmp);
    /**/
    
    fclose(inp);
    printf("Driftkammer Monitor: Normal stop\n");
    return(0);
}