DETECTOR SPECIFIC PARAMETERS (Lines may be commented with '#' sign!)
===================================================================================================================================================================================================================================
DETECTOR   DETECTOR   DETECTOR STEP   TIMING    MIXED     EXPAND   SMOOTHING     TIME          G-FLASH    G-FLASH     G-FLASH    G-FLASH   BASELINE   BASELINE   AMPLITUDE   AMPLITUDE   AREA/AMP.   AREA/AMP.   SIGNAL WIDTH   SIGNAL WIDTH    NUMBER OF     PULSE SHAPE
  NAME      NUMBER     CLASS   SIZE   FILTER   POLARITY   PULSES    FILTER       LIMIT         OPTION    THRESHOLD   MIN_WIDTH   WINDOW     OPTION     FILTER     OPTION     THRESHOLD   LOW THR.    HIGH THR.     LOW THR.       HIGH THR.    PULSE SHAPES     ADDRESS
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PKUP          0         PSA   350/6      0        0          3       100     100000              0         100.        1.         0          -1       300        0              100         0.0        2000              1            4000            0

#################################
#    MWD PARAMETERS V6.2    #
#    Date: 10-Aug-2023      #
#################################
#################################

#                             Threshold   Polarity  PoleZeroCorr    DeconWindow      AvergWin    MAWindows  gain/offset   Gamma_Threshold    Amp_Threshold    fixed Dead time(ns)  time_diff_baselne    tailtime      time_const

EDET         1         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         2         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         3         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         4         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         5         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         6         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         7         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         8         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         9         MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         10        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         11        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         12        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         13        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         14        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         15        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20
EDET         16        MWD      1600            -1         5E5           150              40      1 60         16.0/0.0        100                6500             550                    1500            2      20

EDET         17        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         18        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         19        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         20        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         21        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         22        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         23        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         24        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         25        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         26        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         27        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         28        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         29        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         30        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         31        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20
EDET         32        MWD      800             1         5E5           150             40      1 60         16.0/0.0        100                6800             550                    1500            2      20

DEED         1         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200            550                    1500            2      20
DEED         2         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         3         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         4         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         5         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         6         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         7         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         8         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         9         MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         10        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         11        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         12        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         13        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         14        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         15        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20
DEED         16        MWD       30             -1         5E5            150             40      1 51         1.0/0.0         100                200             550                    1500            2      20

/*
 * $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.;