/*
 * $Id$
 *
 * Author:	Sarah-Jane Lonsdale
 * Date:	02-Dec-2015
 * Version:     5.1 (8-Aug-2017)
 * Update:	Claudia Lederer-Woods
 * Date:	31-May-2018
 * Version:     6.0 (31-May-2018)
 * Update:	Nikolay Sosnin
 * Date:	24-Nov-2022
 * Version:     6.1 (24-Nov-2022)
 */

#include <math.h>
#include <string.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <vector>
#include <TFile.h>
#include <TH1F.h>
#include <TString.h>

using namespace std;

#include "MWDdetector.h"

void FormatHist(TH1D* h, TString name, TString title, TString xtitle, TString ytitle, int color, int width, int marker_color, int marker_style);

bool MWDDetector::parseConfigLine(char* line, const char* settings_file)
{
	cout << "Config: " << line << endl;
	if (!Detector::parseConfigLine(line)) return false;

	// threshold
	char* pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid Threshold" << endl;
		return false;
	}
	threshold = atof(pch);

	// polarity
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid negative polarity" << endl;
		return false;
	}
	polarity = atoi(pch);
	if (polarity>=0)
		polarity = 1;
	else
		polarity = -1;

	// pole zero correction
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid pole zero correction" << endl;
		return false;
	}
	pz = atof(pch);

	// deconvolution window
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid deconvolution window" << endl;
		return false;
	}
	m = atoi(pch);

	// average window
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid average window" << endl;
		return false;
	}
	l = atoi(pch);

	// presample
	pch = strtok(NULL," "); //NS: 27.07.2023
	if (!pch) {
		cerr << "UserInput: Invalid presample" << endl;
		return false;
	}
	presample = atoi(pch);
	
	//averaging windows 2
	pch = strtok(NULL," "); //NS: 27.07.2023
	if (!pch) {
		cerr << "UserInput: Invalid presample" << endl;
		return false;
	}
	window = atoi(pch);

	// presample/averaging window 2
	/*if (!parse2real(&presample, &window)) {
		cerr << "UserInput: Invalid averager 1/2" << endl;
		return false;
	}*/

	// gain/offset
	if (!parse2real(&gain, &offset)) {
		cerr << "UserInput: Invalid gain/offset" << endl;
		return false;
	}

	// gamma flash search start
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid gamma_threshold" << endl;
		return false;
	}
	g_threshold = atoi(pch);

	// minimum amp_threshold
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid amp_threshold" << endl;
		return false;
	}
	amp_threshold = atoi(pch);

	// gamma flash primary (additional) offset
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: Invalid fixed Dead time(ns)" << endl;
		return false;
	}
	gamma_time_primary = atoi(pch);//SL 08/07/17

	// time different for baseline determination
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: time_diff_baselne" << endl;
		return false;
	}
	time_diff_baseline = atoi(pch);//SL 08/07/17

	// time after g flash which has undershoot
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: tailtime" << endl;
		return false;
	}
	tailtime = atoi(pch);//SL 08/07/17


	// Time constant
	pch = strtok(NULL," ");
	if (!pch) {
		cerr << "UserInput: time_const" << endl;
		return false;
	}
	time_const = atoi(pch);  // 28/08/23


	return true;
} // parseConfigLine

int MWDDetector::analysis(
				ntof::lib::ReaderStructEVEH& eveh,	// EVEH event information
				ntof::lib::ReaderStructMODH& modh,	// MODH header information
				ntof::lib::ReaderStructACQC& acqc,	// ACQC pulse record
				PulseVector*   pulsevec,		// vector of pulses
				int movie_number,
				bool html)				// I: movie number
{
	Detector::analysis(eveh, modh, acqc, pulsevec, movie_number, html);

	double rate = modh.getSampleRate();
	int NofPeaks = 0;
	TString name = modh.getDetectorType();
	
	//double g_threshold = 2000; //move to .h file                           //CHANGED
	double tdiffsig = gamma_time_primary; 
	double tdiffbase = time_diff_baseline; // consider moving to h file
	int aver = presample;             //
	//double tailtime = 100000;  // time up to which there is an undershoot after gflash

	double* x          = new double[length_of_movie];
	double* xsmooth    = new double[length_of_movie]; //13.3.18
	double* xsmooth_ma = new double[length_of_movie]; //NS 22.05.2023
	double* y          = new double[length_of_movie];
	double* z          = new double[length_of_movie];
	double* zdiff      = new double[length_of_movie]; //13.3.18
	double* mwd_m      = new double[length_of_movie];
	double* ma_l       = new double[length_of_movie];	// moving average array
	double* mwd_deriv  = new double[length_of_movie];	// mwd derivative

	double timeScale = 1000.0 / rate; // in ns/Sample

	for(int i = 0; i < length_of_movie; i++){
	
		x[i] = polarity * (acqc[i] * gain + offset);
		y[i] = z[i] = zdiff[i] = xsmooth[i] = 0.0;         //13.3.18
	}
                                                //CHANGED
	int startofevent = aver / 2;
	
	//Averaging preamplifier output
	for(int i = startofevent; i < length_of_movie; i++){
		
		xsmooth[i] = 0.;
		
		for(int j = -1 * (aver - 1) / 2; j <= (aver - 1) / 2; j++){
			
			xsmooth[i] += x[i + j];                         
		}
		
		xsmooth[i] /= aver;
	}

	//Extra moving average test
	//const double window = 39.;
	int startofevent2 = window / 2; //NS 22.05.2023
	
	for(int i = startofevent2; i < length_of_movie; i++){
		
		xsmooth_ma[i] = 0.;
		
		for(int j = -1 * (window - 1) / 2; j <= (window - 1) / 2; j++){
			
			xsmooth_ma[i] += xsmooth[i + j];                         
		}
		
		xsmooth_ma[i] /= window;
	}

	// Locating maximum and minimum of derivative
	if(verbose){cout << "Begin amplitude extraction." << endl;}

	// Parameters for semi-gauss discriminator
	//const double time_const = 0; // SL 12/08/16 Emmanuel
	const double pole_zero = 5.e5;
	double a0, a1, b1;

	b1 = exp(-1. / int(time_const));
	a0 = (1. + b1) / 2.0;
	a1 = -1. * (1. + b1) / 2.0;

	// Single pole high pass with pz correction
	for(int i = 1; i < length_of_movie; i++){
		
		//y[i] = b1 * y[i - 1] + a0 * xsmooth[i] + a1 * xsmooth[i - 1] + xsmooth[i - 1] / pole_zero;
		y[i] = b1 * y[i - 1] + a0 * xsmooth_ma[i] + a1 * xsmooth_ma[i - 1] + xsmooth_ma[i - 1] / pz;
	}

	// Single pole low pass filter
	for(int i = 1; i < length_of_movie; i++){
	
		z[i] = b1 * z[i - 1] + a0 * y[i];
	}

	// devirative of filters
	for(int i = 1; i < length_of_movie; i++){
	
		zdiff[i] = -1. * z[i - 1] + z[i];
	}

	//bool beamType = false;
	
	//if(eveh.getBeamType() != 1){beamType = true;}

	double twait = 16000.; // introduce different time window for baseline depending on ded or par CLW 31/05/18

	/*if(eveh.getBeamType() == 2){twait = 13000.;} // dedicated
	if(eveh.getBeamType() == 3){twait = 19000.;} // parasitic
	else{twait = 16000.;}*/     

	////////////////////////////////////////
	//          MWD and filtering         //
	////////////////////////////////////////

	// moving window deconvolution
	for(int i = startofevent + m; i < length_of_movie; i++){ // SL 10/08/16   // 13/03/18 MWD on smoothed 
	
		//double d_m  = xsmooth[i] - xsmooth[i - m];
	  double d_m  = xsmooth_ma[i] - xsmooth_ma[i - m];
		double ma_m = 0.;

		for(int j = (i - m); j < (i - 1); j++) {
		
			//ma_m += xsmooth[j];
			ma_m += xsmooth_ma[j];
		}

		mwd_m[i] = d_m + ma_m / pz;
	}

	// moving average
	//for(int i = l + m + gamma_flash; i < length_of_movie; i++){ // SL 10/08/16
	for(int i = l + m; i < length_of_movie; i++){ // NS 31.07.2023
		
		ma_l[i] = 0.;

		for(int j = (i - l); j < (i - 1); j++){
			
			ma_l[i] += mwd_m[j];
		}
		
		ma_l[i] /= (double)(l);
	}

	////////////////////////////////////////
	// extracting baseline and amplitudes //
	////////////////////////////////////////

	int amptime = 0;
	double zthreshold = 0.;
	double baselineold = 0; //CHANGED
	double timeold = 0;
	double A = 0.;
	double tflash = 0.;
	
	vector<double> Signals;
	Signals.resize(length_of_movie);
	std::fill(Signals.begin(), Signals.end(), 0.);

	for(int i = startofevent + 3 * m; i < length_of_movie; i++){

		double baseline = 0.;

		// z threshold to compensate undershoot at low tofs
		// find number that corresponds to 100 us 

		int t100 = int(tailtime / timeScale);

		if((i + timestamp) >= t100){zthreshold = amp_threshold;} //amp_threshold should be approx 400

		// in undershoot region look for local minimum

		if((i + timestamp) < t100){
		
			zthreshold = 0.;
			
			if(z[i] > amp_threshold){zthreshold = amp_threshold;}

			if(z[i] <= amp_threshold){
			
				for(int h = i - 30; h < i; h++){
					
					if(z[h] < zthreshold){zthreshold = z[h];}
				}
		
				zthreshold = zthreshold + amp_threshold;
			}   
		}

		if(z[i] >= zthreshold && zdiff[i - 1] >= 0. && zdiff[i] <= 0. && timeScale * (i + timestamp) > 10666.){

			//double maxpos = timeScale * (i + timestamp);
			double signaltime = 0.; 
			double baseline1 = 0., baseline2 = 0., baseline3 = 0.;
			
			for(int revk = i - int(3. * time_const); revk <= i - 6; revk++){

				int k = (i - 6) - revk + (i - int(3 * time_const));

				float diff1 = (zdiff[k] - zdiff[k - 3]);
				float diff2 = (zdiff[k + 3] - zdiff[k]);                                       

				// signal time is defined as maximum in z devirative
				if(diff1 >= 0 && diff2 <= 0){

					baseline = 0.;

					// calculate baseline. if signals too close, baseline is baseline of previous signal
					if((timeScale * (k + timestamp) - timeold) > tdiffbase){ /// signal in ns !!
						
						for(int h = k - 40; h < k - 10; h++){
							
							baseline = baseline + ma_l[h];
						}
						
						baseline = baseline / 30.;
					}
					
					if((timeScale * (k + timestamp) - timeold) <= tdiffbase){baseline = baselineold;}
					
					// if signal too close to gamma flash find baseline after signal
					if(timeScale * (k + timestamp) < twait){

						for(int h = k + 2 * m; h < k + 2 * m + 20; h++){
						
							baseline1 = baseline1 + ma_l[h];
						}
						
						baseline1 = baseline1 / 20.;
						baseline = baseline1;

						for(int h = k + m + 5; h < k + m + 5 + 20; h++){
						
							baseline2 = baseline2 + ma_l[h];
						}
						
						baseline2 = baseline2 / 20.;

						if(baseline2 < baseline){baseline = baseline2;}

						for(int h = k - 10; h < k - 5; h++){
						
							baseline3 = baseline3 + ma_l[h];
						}
						
						baseline3 = baseline3 / 5.;  

						if(baseline3 < baseline){baseline = baseline3;}                             

						if(baseline > 0.){baseline = 0.;}
					}

					A = -10000.;

					// find maximum amplitude of moving average
					for(int c = k + m - l; c < k + m + 5; c++){ //CHANGED
					
						if(ma_l[c] > A){A = ma_l[c]; amptime = c;}
					}

					// Calculate area of signal over 2m
					float area = 0.;

					for(int h = k; h <= k + 2 * m; h++){

						area = area + ma_l[h]; 
					}
					
					if(A - baseline > threshold){

						signaltime = timeScale * (k + timestamp);

						if(signaltime - timeold > tdiffsig/* && !(((amptime - k - m + 5) >= -1) && z[i] < A)*/){

							PulseInfo pulse;
							pulse.tof = signaltime; 
							pulse.amp = A - baseline;

							pulse.area = area - 2. * m * baseline; // NEW 
							pulse.area_0 = baseline;      // NEW 
							pulse.fwhm = zthreshold;  // NEW 
							pulse.risetime = z[i];  // NEW
							pulse.best_pulse = amptime - k - m + 5; // NEW 

							pulse.afast = timeScale * (i + timestamp); //NEW //CHANGED
							
							if(modh.getDetectorType() == "DEED"){Signals[i] = A - baseline;}
							else{Signals[i] = A - baseline;}
							
							if(tflash == 0.){
							
								if(signaltime < 12500 && A > g_threshold){pulse.tflash = signaltime; tflash = signaltime;}
							}							
							
							else{pulse.tflash = tflash;}

							pulsevec->push_back(pulse); //CHANGED
							NofPeaks++; //CHANGED
							
							timeold = signaltime;
							baselineold = baseline;
						}
						
						revk = revk + int(4. * time_const); //CHANGED
					}
				}//closer for revk
			} // closer for diff1 if statmement
		} //closer for z[i] loop
	} //closer for int i loop
	
if((eveh.getEventNumber() == 21) && ((modh.getDetectorType() == "DEED" && modh.getDetectorId() == 1))){


/*

if((eveh.getEventNumber() == 21 || eveh.getEventNumber() == 24) && ((modh.getDetectorType() == "DEED" && modh.getDetectorId() == 1) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 7) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 16)|| (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 19) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 23) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 30) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() == 1) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() == 7) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() == 13) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() == 15))){


	if((eveh.getEventNumber() == 21 || eveh.getEventNumber() == 31 || eveh.getEventNumber() == 36 || eveh.getEventNumber() == 37 || eveh.getEventNumber() == 38 || eveh.getEventNumber() == 45) && ((modh.getDetectorType() == "EDET" && modh.getDetectorId() == 1) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 17) || (modh.getDetectorType() == "EDET" && modh.getDetectorId() == 24) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() == 1) || (modh.getDetectorType() == "DEED" && modh.getDetectorId() > 10))){
	*/ 

		TFile* rootout = TFile::Open((TString)"Traces/Trace_" + modh.getDetectorType() + Form("%d_run%d_segment%d.root", modh.getDetectorId(), eveh.getRunNumber(), eveh.getEventNumber()), "recreate");
		
		TH1D* original = new TH1D("original", "original", length_of_movie, 0, length_of_movie);
		TH1D* x_ma1 = new TH1D("x_ma1", "x_ma1", length_of_movie, 0, length_of_movie);
		TH1D* x_ma2 = new TH1D("x_ma2", "x_ma2", length_of_movie, 0, length_of_movie);
		TH1D* hy = new TH1D("y", "y", length_of_movie, 0, length_of_movie);
		TH1D* hz = new TH1D("z", "z", length_of_movie, 0, length_of_movie);
		TH1D* hzdiff = new TH1D("zdiff", "zdiff", length_of_movie, 0, length_of_movie);
		TH1D* hmwd = new TH1D("mwd", "mwd", length_of_movie, 0, length_of_movie);
		TH1D* hma_l = new TH1D("ma_l", "ma_l", length_of_movie, 0, length_of_movie);
		TH1D* signal_id = new TH1D("signal_id", "signal_id", length_of_movie, 0, length_of_movie);
		
		FormatHist(original, "original", "original", "Time (chan.)", "Amplitude (chan.)", 1, 1, 1, 1);
		FormatHist(x_ma1, "x_ma1", "x_ma1", "Time (chan.)", "Amplitude (chan.)", 13, 1, 13, 1);
		FormatHist(x_ma2, "x_ma2", "x_ma2", "Time (chan.)", "Amplitude (chan.)", 2, 1, 2, 1);
		FormatHist(hy, "y", "y", "Time (chan.)", "Amplitude (chan.)", 7, 1, 7, 1);
		FormatHist(hz, "z", "z", "Time (chan.)", "Amplitude (chan.)", 4, 1, 4, 1);
		FormatHist(hzdiff, "zdiff", "zdiff", "Time (chan.)", "Amplitude (chan.)", 46, 1, 46, 1);
		FormatHist(hmwd, "mwd", "mwd", "Time (chan.)", "Amplitude (chan.)", 15, 1, 15, 1);
		FormatHist(hma_l, "ma_l", "ma_l", "Time (chan.)", "Amplitude (chan.)", 6, 1, 6, 1);
		FormatHist(signal_id, "signal_id", "signal_id", "Time (chan.)", "Amplitude (chan.)", 801, 3, 801, 1);
		
		float avg_baseline1 = 0.;
		float avg_baseline2 = 0.;
		float avg_baseline3 = 0.;
		float N = 0.;
		
		for(int i = 200; i < 400; i++){
		
			avg_baseline1 += x[i];
			avg_baseline2 += xsmooth[i];
			avg_baseline3 += xsmooth_ma[i];
			
			N++;
		}
		
		if(N > 0.){
		
			avg_baseline1 /= N;
			avg_baseline2 /= N;
			avg_baseline3 /= N;
		}
		
		else{cout << "WARNING: Average baseline over zero channels!" << endl;}
		
		for(int i = 0; i < length_of_movie; i++){
		
			original->SetBinContent(i, x[i] - avg_baseline1);
			x_ma1->SetBinContent(i, xsmooth[i] - avg_baseline2);
			x_ma2->SetBinContent(i, xsmooth_ma[i] - avg_baseline3);
			hy->SetBinContent(i, y[i]);
			hz->SetBinContent(i, z[i]);
			hzdiff->SetBinContent(i, zdiff[i]);
			hmwd->SetBinContent(i, mwd_m[i]);
			hma_l->SetBinContent(i, ma_l[i]);
			signal_id->SetBinContent(i, Signals[i]);
		}
		
		rootout->Write();
		rootout->Close();
		
		cout << modh.getDetectorType() << modh.getDetectorId() << " avg baseline 1: " << avg_baseline1 << ", avg baseline 2: " << avg_baseline2 << ", avg baseline 3: " << avg_baseline3 << " (MA window 1: " << aver << ", MA window 2: " << window << ")" << endl;
	} 

	
	delete [] x;
	delete [] y;
	delete [] z;
	delete [] zdiff;
	delete [] xsmooth;
	delete [] xsmooth_ma;
	delete [] mwd_m;
	delete [] ma_l;
	delete [] mwd_deriv;

	return NofPeaks;
} /* analysis */

void FormatHist(TH1D* h, TString name, TString title, TString xtitle, TString ytitle, int color, int width, int marker_color, int marker_style){
	
	h->SetName(name);
	h->SetTitle(title);
	
	h->GetXaxis()->SetTitle(xtitle);
	h->GetYaxis()->SetTitle(ytitle);
	
	h->GetYaxis()->SetTitleOffset(1.);
	
	h->SetLineColor(color);
	h->SetLineWidth(width);
	
	h->SetMarkerColor(marker_color);
	h->SetMarkerStyle(marker_style);
	
	h->SetStats(0);
	
	return;
}