#include <string>
#include <vector>
#include <cmath>
#include <glob.h>
#include <stdio.h>
#include <stdlib.h>
#include <libgen.h>
#include <iostream>
#include <fstream>
#include <regex>
#include <TApplication.h>
#include <TChain.h>
#include <TChainElement.h>
#include <TTree.h>
#include <TBranch.h>
#include <TLeaf.h>
#include <TCut.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TF1.h>
#include <TFitResult.h>
#include <TGraphErrors.h>
#include <TSpectrum.h>
#include <TPolyMarker.h>
#include <TText.h>
#include <TROOT.h>
#include <TComplex.h>
#include <TObjArray.h>
#include <TMinuit.h>
#include <TMath.h>
#include <TVirtualPad.h>
#include <TModel.h>


#define	DEBUG			1
#define	DPRINTF			if (DEBUG & 1) printf
#define	dprintf			if (DEBUG & 2) printf

#define	TDC_LLD			128
#define	TDC_ULD			4095
#define	TDC_PEAK		3500
#define	TDC_CALIB_WIDTH		500
#define	TDC_CALIB_LEFT		(TDC_PEAK-TDC_CALIB_WIDTH)
#define	TDC_CALIB_RIGHT		(TDC_PEAK+TDC_CALIB_WIDTH)

#define	SHOW_SPEC		1
#define	SUBTRACT_BGND		0

#define	SIGMA			2.0

#define	ProjSpecResol		0.01
#define	ProjSpecSigma		SIGMA

#define	BASE_LT			3
#define	BASE_RT			3

#define	MAX(a, b)		((a > b) ? a : b)
#define	MIN(a, b)		((a < b) ? a : b)

#define	NumBins			8192
#define	HistXmin		1
#define	HistXmax		8191
#define	ProjSpecName		"SPEC"
#define	ProjSpecLineWidth	2
#define	ProjSpecLineColor	kMagenta

#define	RECALIBRATE		0.5

#define	SkewL			0
#define	SkewR			0
#define	SkewedGauss		(SkewL+SkewR)
#define	QuadBack		0
#define	MIN_AREA		50

#define	NumPolyMax		5
#define	NumPoly			3

//a global variable to store the current histogram	
TH1F* gCurrentHistogram = nullptr;
using namespace std;
Double_t X[100]; // the fitted peaks
Double_t selected_Peaks[6]; // peaks used for calibration in calib();
 
void errorExit (const char*fmt, ...)
{
va_list ap;
	va_start (ap, fmt);
	vprintf (fmt, ap);
	va_end (ap);
	exit (0);
}


// globFiles: a library to get "(file.*)"
std::vector<std::string> globFiles(const std::string& pattern)
{
    std::vector<std::string> files;
    glob_t globResult;

    if (glob(pattern.c_str(), GLOB_TILDE, nullptr, &globResult) == 0) {
        for (size_t i = 0; i < globResult.gl_pathc; ++i) {
            files.push_back(globResult.gl_pathv[i]);
        }
        globfree(&globResult);
    }

    return files;
}
TChain *t = new TChain ("RoseNIAS");

std::vector<std::string> ListBranches ()
{		
		
		TObjArray *branches = t -> GetListOfBranches();
		//create the search pattern (step1)
		std::cout<<"Enter the Clover number (e.g. 22 04 would mean CL_22_E04)"<<std::endl;
		std::string searchNumber;
		std::string searchCrystal;
		std::cin>>searchNumber>>searchCrystal;
		
		std::string searchPattern = "_" + searchNumber + "_E" + searchCrystal; 
		// all matching clovers for the search (step1)
		
		std::vector<std::string> matchingBranches; 
		for (int i = 0; i < branches->GetEntries(); ++i) 
		{
       		TBranch *branch = dynamic_cast<TBranch *>(branches->At(i));
        		if (branch) {
            			std::string branchName = branch->GetName();
            			// Check if branchName contains the searchPattern
            			if (branchName.find(searchPattern) != std::string::npos) {
                		matchingBranches.push_back(branchName);
                		std::cout << "Branch :  " << branchName << std::endl; }
							 }
		}
		if (!matchingBranches.empty())
		{
			return matchingBranches;
		}
		else
		{
			std::cout << "No matching branches found." << std::endl;
			return std::vector<std::string>(); // Return an empty string if no matching branches
		}
}		



//	PROGRAM TO FIT A GAUSSIAN
Int_t xx1, xx2;

enum Parameters {
Peak, Centroid, Sigma, SkewLpar, SkewRpar, BackA, BackB, BackC, Area, BackArea, PeakY, NumPar };

Double_t background (Double_t *x, Double_t *par)
{
Double_t back;
	back = par[0]+ par[1]*x[0];
	if (QuadBack) back += par[2]*x[0]*x[0];
	return back;
}

Double_t gauss (Double_t*x, Double_t *par)
{
Double_t peak, arg;
	
	arg = (x[0] - par [1]) / par [2];
	if (SkewL && (arg <= -par[3]))
		peak = exp (0.5*par[3]*par[3] + par[3]*arg);
	else if (SkewR && (arg >= par[4]))
		peak = exp (0.5*par[4]*par[4] + par[4]*arg);
	else
		peak = exp (-0.5*arg*arg);
	return par[0]*peak;
}

Double_t skewedGauss (Double_t *x, Double_t *par)
{
	return (gauss(x, par) + background (x, &par [3+SkewedGauss]));
}

Double_t fitPeak (TH1F*hist, double u, Double_t*par)
{

TF1		*f1, *f2;
Double_t	chi=1e10;
Double_t	x1, y1, x2, y2;
TFitResultPtr	fitres;
int		result, ndf, pol1Used = 0;

    x1 = u - 1.5 * SIGMA;
    x2 = u + 1.5 * SIGMA;

    f1 = new TF1 ("f1", "gaus", x1, x2);
    f1->SetLineColor (kGreen);
    f1->SetLineWidth (3);
    fitres = hist->Fit (f1, "QSR0");
    result = int (fitres);
    f1->GetParameters (par);
    cout<<"LOOK: the parameters par[n] are: "<<par[0]<< " "<<par[1]<<" "<<par[2]<<endl;
    chi = f1->GetChisquare ();
    ndf = f1->GetNDF ();

    x1 = u - BASE_LT * par [2];
    x2 = u + BASE_RT * par [2];

    {
	y1 = hist->GetBinContent (x1);
	y2 = hist->GetBinContent (x2);
	par [3] = par [4] = 1.;
	par [4+SkewedGauss] = (y2 - y1) / (x2 - x1);
	par [3+SkewedGauss] = y1 - par [4+SkewedGauss] * x1;
	par [5+SkewedGauss] = par [4+SkewedGauss] * 1e-4;
//	par [3] = par [4] = par [5] = par [6] = par [7] = par [8] = 1.;
	f2 = new TF1 ("f2", skewedGauss, x1, x2, 5+SkewedGauss+QuadBack);
	f2->SetLineColor (kRed);
	f2->SetLineWidth (3);

	f2->SetParameters (par);
	fitres = hist->Fit (f2, "QR0");
	fitres = hist->Fit (f2, "QSR0");
	result = int (fitres);
	f2->GetParameters (par);
	if (! result) {
	    chi = f2->GetChisquare ();
	    ndf = f2->GetNDF ();
	    pol1Used = 1;
		cout<<"the value of ChiSquare is: "<<chi<<endl;
	}
    }


    dprintf ("%8.2lf %8.2lf %8.2lf CHI=%.1lf\n",
		par [0], par [1], par [2], chi);
    if ((par [2] < 1) || (par [0] < MIN_AREA) || result
		|| (fabs (u-par [1]) > (BASE_RT*SIGMA)))
	return -1;

    for (int p=0; p<5; ++p) dprintf ("\tpar[%d] = %.1lf\n", p, par [p]);
    dprintf ("RES=%d CHI=%.1lf NDF=%d RedChi=%.1lf\n",
		result, chi, ndf, chi / ndf);

    if (pol1Used) f2->Draw ("same"); else f1->Draw ("same");

    TF1*sig = new TF1 ("sig", gauss, x1, x2, SkewedGauss+3);
    sig->SetParameters (par);
    sig->SetLineColor (kTeal);
    sig->SetLineWidth (3);
    sig->Draw ("same");

    TF1*back = new TF1 ("back", background, x1, x2, 2+QuadBack);
    back->SetParameters (&par [3+SkewedGauss]);
    back->SetLineColor (kBlue);
    back->SetLineWidth (3);
    back->Draw ("same");

    par [BackArea] = 0.;
    if (pol1Used)
	par [BackArea] = back->Integral (x1, x2);
    xx1 = (par [1] - par [Sigma] * 2.35 - 0.5);
    xx2 = (par [1] + par [Sigma] * 2.35 + 0.5);
    par [Peak] = sig->Integral (x1, x2);

    dprintf ("***u=%8.1lf (%5.0lf-%5.0lf) A=%8.0lf Area=%8.0lf\n",
		u, x1, x2, par [Peak], par [Area]);

    return chi / ndf;
}
// end of fitting funciton

Double_t IdentifyPeaks(TH1F *H, Double_t *X)

{
/*TSpectrum *spec = new TSpectrum();
Int_t numP = spec -> Search(H, 3, "g", 0.05);
H->Draw("same");
*/
//alternate option with ShowPeaks()	
	Int_t	numP = H->ShowPeaks (ProjSpecSigma, "", ProjSpecResol);
	if (! numP) {cout<<"the code didnt work";}
    cout<<"the number of peaks is: "<<numP<<endl;
    
    TList	*funcP = H->GetListOfFunctions ();
    if (! funcP) return 0;
    TPolyMarker	*mX = (TPolyMarker*)funcP->FindObject ("TPolyMarker");
    if (! mX) return 0;
    Double_t	tX, *peaks = mX->GetX ();
    if (! peaks) return 0;
	cout<< "Peaks have been found and now are being stored in X[]"<<endl;
    for (int p=0; p<numP; ++p) {
	for (int q=p+1; q<numP; ++q) {
	    if (peaks [p] > peaks [q]) {
		tX = peaks [p];
		peaks [p] = peaks [q];
		peaks [q] = tX;
	    }
	}
    }
	TText *markerPrint = new TText();
	
    Int_t numX = 0;
	
	Double_t par[10];
    for (int i=0; i<numP; i++) {
    cout<<"peak "<< i <<" is at "<<peaks[i]<<"\n"<<endl;
    Double_t rChi = fitPeak (H, peaks [i], par);
    cout<<"the value of rchi is "<<rChi<<endl;
    if (rChi <0) {dprintf ("\t\tthe loop for rChi::%3d. %8.2lf %8.2lf %8.2lf %8.2lf\n",
			i+1, peaks [i], rChi, par [1], par [SkewLpar]);
	    continue;
	}

		X [numX] = par[1];
	
    if (X [numX] > 0) 
    { printf ("the loop for X[numX]::%3d. %8.2lf %8.2lf %8.2lf %8.2lf\n", i+1, peaks [i], X [numX], par [Sigma], par [SkewLpar]);
	
	cout<<"the value of numX is: "<<numX<<endl;
	numX++;
	
    }
	
	}
	// Writing on the canvas, the number of peaks fitted to be used for calibration
	for (int j = 0; j < numX; j++)
	{
		int jj = j+1;
	Double_t y_cord_marker = H -> GetBinContent(X[j]) + 150;
	Double_t x_cord_marker = X[j] - 2;
	TText *markerPrint = new TText(x_cord_marker, y_cord_marker, Form("%d", jj));
	markerPrint -> SetTextSize(0.03);
	markerPrint -> Draw("SAME");
	}
   return numX;

}


//void BinClick(Int_t event, Int_t x, Int_t y, Int_t /*px*/, Int_t /*py*/) {
//    if (event == kButton1Down) { // Check if left mouse button is clicked
//        TCanvas *c1 = (TCanvas*)gPad->GetCanvas();
//        TH1F *h = (TH1F*)c1->GetPrimitive("h"); 
//        if (h) {
//            Int_t bin = h->GetXaxis()->FindBin(c1->AbsPixeltoX(x));
//            std::cout << "Clicked on bin number: " << bin << std::endl;
//        }
//    }
//}
/*int BinClick ()
{
	Double_t x1x = 0;
	Double_t y1y = 0;
	Int_t event = gPad -> GetEvent();
	int px = gPad -> GetEventX();
	int py = gPad -> GetEventY();
	TObject *select = gPad -> GetSelected();
	if (!select) return -1;
	Double_t x1x_to_pixel = gPad -> AbsPixeltoX(px);
	Double_t y1y_to_pixel = gPad -> AbsPixeltoY(py);
	x1x = gPad -> PadtoX(x1x_to_pixel);
	y1y = gPad -> PadtoY(y1y_to_pixel);
	return x1x;
}
void SetBinClickHandler() {
    gPad->AddExec("mouse_click", "BinClick");
}
*/
TH1F* ShowSpectra (const std::string& branchName)
{	TCanvas *c = new TCanvas ("Canvas 1", branchName.c_str());
	TH1F *h = new TH1F ("h", branchName.c_str(), 8192,0,8191);
	t -> Draw((branchName+">>h(8192,0,8191)").c_str());
	
	//SetBinClickHandler();
	//gPad -> WaitPrimitive();
	//gPad->AddExec("mouse_click", "BinClick");
	return h;
}


void ShowSpectraGrid(const std::vector<std::vector<std::string>> &branchNames, int rows, int columns)
{
    TCanvas *c1 = new TCanvas("Canvas 1", "Grid");
    c1->Divide(rows, columns);

    // Create an array of histograms
    TH1F *histograms[rows * columns];

    for (int i = 0; i < rows * columns; ++i)
    {
        c1->cd(i + 1);
        // Create a new histogram for each branch
        histograms[i] = new TH1F(("h" + std::to_string(i)).c_str(), branchNames[i][0].c_str(), 8192, 0, 8191);
        t->Draw((branchNames[i][0] + ">>h" + std::to_string(i)).c_str());
        histograms[i]->Draw();
        c1->Update();
    }

    c1->Update();
    c1->Draw();
}

void ShowSpectraOverlay(const std::vector<std::vector<std::string>> &branchNames, int numSpectra)
{
    TCanvas *c2 = new TCanvas("CanvasOverlay", "Overlay");
    c2->cd();

    // Create a legend to identify each spectrum
    TLegend *legend = new TLegend(0.7, 0.7, 0.9, 0.9);

    // Create an array of colors for each spectrum
    int colors[] = {kRed, kBlue, kGreen, kOrange, kMagenta, kCyan};
	
	TH1F *hist[numSpectra];
    // Loop over the selected spectra
    for (int i = 0; i < numSpectra; ++i)
    {
        // Draw the spectrum with a different color
        t->SetLineColor(colors[i]);
        hist[i] = new TH1F (("hist" + std::to_string(i)).c_str(), branchNames[i][0].c_str(), 8192, 0, 8191);
        if (i == 0)
        {
		t->Draw((branchNames[i][0] + ">>hist" + std::to_string(i)).c_str());
		hist[i] -> SetLineColor(colors[i]);
	}
         else
            { t->Draw((branchNames[i][0] + ">>hist" + std::to_string(i)).c_str(), "", "same");
            hist[i] -> SetLineColor(colors[i]);
		}
        
        //t->Draw((branchNames[i][0] + ">>h(8192,0,8191)").c_str(), "", i > 0 ? "same" : "");

        // Add the spectrum to the legend
	legend->AddEntry(hist[i], branchNames[i][0].c_str(), "l");
    }

    // Draw the legend
    legend->Draw();
    c2->Update();
    c2->Draw();
}

void view(){
		std::cout<<"Do you wish to view OR overlay spectra? (v/o)"<<std::endl;
		char xx;
		std::cin>>xx;
		
		if(xx=='v' || xx == 'V')
		{	
			// view spectrum code
			std::cout<<"Do you want to see a single spectrum or different spectra in a grid? (g/s)"<<std::endl;
			char yy;
			std::cin>>yy;
			
			if(yy == 'g' || yy == 'G')
			{
				// call the function to make a grid and show the spectra
				std::cout<<"Split the window into a grid of? (e.g. for 6, enter 2 3)";
				int rows, columns;
				std::cin>> rows >> columns;
				
				
				std::vector<std::vector<std::string>> branches;	
					for (int i = 0; i<(rows*columns); ++i)
					{
						std::vector<std::string> branchNames_for_grid = ListBranches();
						branches.push_back(branchNames_for_grid);
					}
					std::cout<<"The signals you chose for the grid are: "<<std::endl;
					for (int i = 0; i<(rows* columns); ++i)
					{
						std::cout<<branches[i][0]<<std::endl;
					}
					ShowSpectraGrid(branches, rows, columns);
			}		
			else if (yy == 's' || yy == 'S')
			{
				//call the function to show a single spectrum
				std::vector<std::string> selectedBranch = ListBranches ();
				
				TH1F* h = ShowSpectra(selectedBranch[0]);
				cout << "Do you wish to fit the peaks? Type   ---> F();" << endl;
				
				gCurrentHistogram = h;
				
			}
			else 
			{
				std::cout<<"Please choose a valid response"<<std::endl;
			}
			
		}
		
		else if (xx=='o' || xx =='O')
		{
			// overlay spectrum code
			std::cout << "Enter the number of spectra to overlay: ";
			int numSpectra;
			std::cin >> numSpectra;

			// Get the branches for overlay
			std::vector<std::vector<std::string>> overlayBranches;
			for (int i = 0; i < numSpectra; ++i)
			{
				std::vector<std::string> branchNames_for_overlay = ListBranches();
				overlayBranches.push_back(branchNames_for_overlay);
			}

			// Call the function to show overlay spectra
			ShowSpectraOverlay(overlayBranches, numSpectra);
			
		}
		
		else
		{
			std::cout<<"Please choose a valid response"<<std::endl;
			// some unknown functionality (to be added later)
		}
		
		std::cout<<"If you wish to view other spctra, type -->  view();"<<std::endl;
	
}

void cal ()
{
	Double_t energiesEu152[15];
	energiesEu152[0] = 121.78;
	energiesEu152[1] = 244.70;
	energiesEu152[2] = 344.28;
	energiesEu152[3] = 778.90;
	energiesEu152[4] = 964.08;
	energiesEu152[5] = 1408.01;

	cout<<"Calibrate by manually feeding the energies or pre-set values? (m/a)";
	char c2;
	cin>>c2;
	if(c2 == 'm'||c2 == 'M')
	{
	//option 1: select the peaks manually
	}
	else if (c2 == 'a'||c2 == 'A')
	{}
	//option 2: select from the set of pre-fitted values
	cout<<"the pre-set energies of Eu152 are the following: (for other energies, change in script) \n";
	for (int i=0; i<6; i++)
	{
		cout<<"energy #"<<i+1<<"\t"<<energiesEu152[i]<<"\n"<<endl;
	}


	//check if X has values
	cout<<"LOOK: the array X is still accessible"<<endl;
	for (int l = 0; l<20; l++)
	{
	cout<<X[l]<<"\n";
	}
	cout<<endl;
	Double_t select_Peak_number[6];
	cout<<"Select the peaks for calib by peak number: (e.g. 1 4 6 7...)"<<endl;
	for (int j =0; j<6; ++j)
	{
	cin>>select_Peak_number[j];
	}
	
	for (int i=0; i<6; ++i)
	{
	int jj = select_Peak_number[i]; //converted array to an integer
	selected_Peaks[i] = X[jj-1];
	cout<<"THE SELECTED ENERGY FOR PEAK #"<<i+1<<" IS: "<<selected_Peaks[i]<<"\n";
	}
	cout<<"do you wish to continue with these peaks or repeat? (y/n)";
	char c3;
	cin>>c3;
	if(c3 == 'y'||c3 == 'Y')
	{
		//code for performing calibration
		
    
    const int num_points = 6;
 
    
    // Define the model function
    auto model = [&](double x, const double* params) {
        return params[0] * x * x + params[1] * x + params[2] + params[3] * sqrt(x);
    };

    // Define the initial guess for parameters
    double params[4] = {1.0, 1.0, 1.0, 1.0};

    // Perform Chi-square minimization
    double chiSquare = 0.0;
    for (int i = 0; i < num_points; ++i) {
        double residual = energiesEu152[i] - model(selected_Peaks[i], params);
        chiSquare += residual * residual;
    }

    // Degrees of freedom
    int dof = num_points - 4; // Assuming 4 parameters

    // Reduced Chi-square
    double reducedChiSquare = chiSquare / dof;

    // Output the results
    cout << "Chi-square: " << chiSquare << endl;
    cout << "Degrees of Freedom: " << dof << endl;
    cout << "Reduced Chi-square: " << reducedChiSquare << endl;

    // You can further analyze the results or use params array for the optimized parameters

	}
	else if (c3 == 'n'||c3 == 'N')
	{
		cal();
	}

}



void F()
{
if (gCurrentHistogram)
{
	Int_t numX;
	
	numX = IdentifyPeaks(gCurrentHistogram, X);
	
	for (int i= 0; i<numX; i++)
	{
	cout<<"peak # "<<i+1<<": \t"<<X[i]<<"\n";
	}
	cout<<"if you wish to calibrate with the marked peaks, type  ----> cal();"<<endl;
	
}
else {cout<<"No histogram available."<<endl;}
}



void Prabhat_Script_mod9()
{
std:: string fileNames;
std::cout<<"Add Files from the Present Working Directory: ";
std::cin>>fileNames;

// use globFiles to get a vector of matching files

std::vector<std::string> matchingFiles=globFiles(fileNames);

for (const auto& file: matchingFiles)
{
t -> Add (file.c_str());
}

t -> ls();
//t -> Print();

		char zz;
		std::cout<<"Do you want to list all the branches individually?"<<std::endl;
		std::cin>>zz;

			if (zz == 'y' || zz == 'Y')
			{
 				TObjArray *branches = t -> GetListOfBranches();
				for (int i = 0; i<branches-> GetEntries(); ++i)
				{TBranch *branch = dynamic_cast<TBranch*>(branches->At(i)); 
				if (branch){
					std::cout<<"Branch "<<i<<": "<<branch->GetName()<<std::endl;
				   	   }
				}
			}
view();
}