//Reads data from a input root tree and calibrate each channel histogram and put the calibrated 
//values in another root file.

#include <stdlib.h>
#include <string>
#include <stdio.h>
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include <iostream>
#include <sstream>
#include <fstream>
#include <TSystem.h>
#include <TROOT.h>
#include <math.h>
#include <string.h>
#include <array>
#include <vector>
#include <thread>
#include <algorithm>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TSpectrum.h>
#include <TPolyMarker.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TPad.h>
#include <TFitResult.h>
#include <TMath.h>
#include <time.h>
#include <TThread.h>
#include <TVector.h>
#include "TROOT.h"



    Double_t    fPositionX[10];
    Double_t    fPositionY[10];
    Int_t       nfound,bin,dummy1,dummy2;
    const Int_t nbins = 1000;
    Double_t    a, source[nbins], dest[nbins], calibX[nbins], par[10],I[1024], S[1024];
	int i,j,k;
	Float_t q[1024], qcalib[1024];
	Double_t Sigma[1024], sigerr, SigErr[1024], sigma, ChannelID[1024], resolution[1024]; 
    fstream f;

//peakcutoff Values for both the Detector Arrays
    float peakcutoff1=10;   //Array on Right
    float peakcutoff2=3;    //Array on Left 
    int   sradius1 = 3;     //radius for identified peak fitting mean +- sradius
    int   sradius2 = 1;     //radius for identified peak fitting mean +- sradius
    int   entries;
    
    const int n_channels    = 128;
    const int n_threads     = 16;
    
    array<TH1F*,n_threads> h;
   // <TTree*,n_threads> data;



/******************************************************************************************/
void* calibrate (void* ptr)
{
    
    int* N = (int*) ptr;
    cout<<"Runnig Thread: "<<N[0]<<"\n";
  
    TTree* data = (TTree*) input->Get("data");
	data->SetBranchAddress("q", &q);
	entries = data->GetEntries();
    
    for(j=1;j<=entries;j++)
    {    
        for(auto CHID=N[1]; CHID<=N[2]; CHID++) 
        {   data->GetEntry(j);
            if(q[CHID]>0)
            {
                h[CHID]->Fill(q[CHID]);
            }
            else{continue;}
        }
    }
    
   
    for(auto CHID=N[1]; CHID<=N[2]; CHID++) 
    {
        int nfound = 0;
     
        //Searching for Peaks in the spectrum
        TSpectrum *s = new TSpectrum();
  
        if(CHID <= 64)
        {
        for (i = 0; i < nbins; i++) {source[i]=h[CHID]->GetBinContent(i+1);}
        nfound = s->SearchHighRes(source, dest, nbins, 2, 1, kFALSE, 1, kTRUE, 3);
        }
        
        if(CHID >= 320)
        {
        for (i = 0; i < nbins; i++) {source[i]=h[CHID]->GetBinContent(i+1);}
        nfound = s->SearchHighRes(source, dest, nbins, 2, 3, kFALSE, 1, kTRUE, 3);
        }
    
        //getting bin center for each found peak
        int npeaks = 0;
        Double_t *xpeaks = s->GetPositionX();
        for (i = 0; i < nfound; i++) 
        {
            a=xpeaks[i];
            bin = 1+Int_t(a+0.5);
            fPositionX[i] = h[CHID]->GetBinCenter(bin);
            fPositionY[i] = h[CHID]->GetBinContent(bin);
        }
	
	   
        //Sorting Peak positions for fitting 
        for(i=0;i<nfound;i++)
        {		
            for(j=i+1;j<nfound;j++)
            {
                if(fPositionX[i]>fPositionX[j])
                {
                    dummy2=fPositionX[i];
                    fPositionX[i]=fPositionX[j];
                    fPositionX[j]=dummy2;
                }
            }//cout<<fPositionX[i]<<"\n";
        }
		
        //using only 0 and 511 peak bin centers for fitting	
        Float_t fPosX[nbins], Energy[nbins], A;
        float peakBinNb;	
        A = 0.0;	
		
    	
        for( i=0;i<nfound;i++)
        {
            if(i == nfound-1)
            {
                A = fPositionX[i];
                fPosX[i]=A;
                peakBinNb = A;
            }
            else
            {
			fPosX[i]=0;
			Energy[i]=0;
            }
		//cout<<fPosX[i]<<"\n";
		}
  

        if(CHID >=1 && CHID < 65)
        {
            if(peakBinNb < peakcutoff1)
            {
            //Can not continue with the calibration
            //cout<<"No peak above cutoff, Skipping Channel\n\n";
            Sigma[CHID] = 0.0;
            SigErr[CHID] = 0.0;
            resolution[CHID] = 0.0;
            continue;	
            }
        auto gausini = new TF1("gausini", "gaus", peakBinNb - sradius1, peakBinNb + sradius1);
        auto fitRes = h[CHID]->Fit("gausini","QRS");
        auto gauspol = new TF1("gauspol", "gaus(0)+pol1(3)", peakBinNb - sradius1, peakBinNb + sradius1);
        gauspol->SetParameter(0,fitRes->Parameter(0)); 
        gauspol->SetParameter(1,fitRes->Parameter(1));
        gauspol->SetParameter(2,fitRes->Parameter(2));
	
        auto fitres = h[CHID]->Fit("gauspol","QRS+");
        Double_t mean = (fitres->Parameter(1));
        Double_t sig = (fitres->Parameter(2));
        Double_t sigerr = (fitres->ParError(2));
        Sigma[CHID] = sig;
        SigErr[CHID] = sigerr;
        resolution[CHID] = (sig*2.35*100)/mean;
	
	
        //cout<<"Energy"<<"\t"<<"Position"<<"\n";
		for( i=0;i<nfound;i++)
		{
			if(i == nfound-1)
			{
			fPosX[i] = mean;
			Energy[i]=511;
			}
			else
			{
			fPosX[i]=0;
			Energy[i]=0;
			}
		
		//cout<<Energy[i]<<"\t"<<fPosX[i]<<"\n";
		}
    }
    
	if(CHID >=320 && CHID < 384)
    {
        if(peakBinNb < peakcutoff2)
        {
        //Can not continue with the calibration
        //cout<<"No peak above cutoff, Skipping Channel\n\n";
        Sigma[CHID] = 0.0;
        SigErr[CHID] = 0.0;
        resolution[CHID] = 0.0;
        continue;	
        }
        auto gausini = new TF1("gausini", "gaus", peakBinNb - sradius2, peakBinNb + sradius2);
        auto fitRes = h[CHID]->Fit("gausini","QRS");
        auto gauspol = new TF1("gauspol", "gaus(0)+pol1(3)", peakBinNb - sradius2, peakBinNb + sradius2);
        gauspol->SetParameter(0,fitRes->Parameter(0)); 
        gauspol->SetParameter(1,fitRes->Parameter(1));
        gauspol->SetParameter(2,fitRes->Parameter(2));
        auto fitres = h[CHID]->Fit("gauspol","QRS");
        Double_t mean = (fitres->Parameter(1));
        Double_t sig = (fitres->Parameter(2));
        Double_t sigerr = (fitres->ParError(2));
        Sigma[CHID] = sig;
        SigErr[CHID] = sigerr;
        resolution[CHID] = (sig*2.35*100)/mean;
	
        //cout<<"Energy"<<"\t"<<"Position"<<"\n";
		for( i=0;i<nfound;i++)
		{
			if(i == nfound-1)
			{
			fPosX[i] = mean;
			Energy[i]=511;
			}
			else
			{
			fPosX[i]=0;
			Energy[i]=0;
			}
		
		//cout<<Energy[i]<<"\t"<<fPosX[i]<<"\n";
		}
    }
	
 
//Fitting procedure	
	     
   	TGraph *gr = new TGraph(nfound,fPosX,Energy);
  	
   	gr->SetMarkerColor(4);
   	gr->SetMarkerStyle(21);
  	gr->Draw("ACP");
  	gr->Fit("pol1","Q");
   	
	//Access the fit resuts
  	TF1 *f3 = gr->GetFunction("pol1");
   	Double_t p0=(f3->GetParameter(0));
   	Double_t p1=(f3->GetParameter(1));
	I[CHID] = p0;	
	S[CHID] = p1;
   	
    delete h[CHID];      //deleting histogram
    
	for (i = 0; i < nfound; i++) //Resetting values  
	{
     	 fPositionX[i] = 0;
     	 fPositionY[i] = 0;
	}
    
}
}   
/****************************************************************************************/

void initialize_plots()
{
    for(int i=0; i<384; i++)
    { if(i>=64 && i<320){continue;}
        h[i] = new TH1F("histo", "Spectrum of Gamma Rays;Bins;Counts", nbins, 0, 100);
    }
}


/****************************************************************************************/

void run_threads()
{
  array<TThread*,n_threads> VT;
  array<Long64_t[3],n_threads> par;
  for(int i=0; i<n_threads/2; i++)
  {
    par[i][0] = i;
    par[i][1] = 1+i*(n_channels/n_threads);
    par[i][2] = (i+1)*(n_channels/n_threads);
    VT[i] = new TThread("histograms", calibrate, (void*) par[i]);
  }
  for(int i=n_threads/2; i<n_threads; i++)
  {
    par[i][0] = i;
    par[i][1] = 255+(1+i*(n_channels/n_threads));
    par[i][2] = 255+((i+1)*(n_channels/n_threads));
    VT[i] = new TThread("histograms", calibrate, (void*) par[i]);
  }
  for(int i=0;i<n_threads;i++) VT[i]->Run();
  for(int i=0;i<n_threads;i++) VT[i]->Join();
}

/****************************************************************************************/

void write_fitting_file()
{
   //Saving Fitting parameters in external file
	cout<<"Writing Parameters in external file.\n\n";
	f.open(Form("Fit_para_%s",oname), fstream::out);
	f<<"ChID"<<"\t"<<"Slope"<<"\t"<<"Sigma"<<"\t"<<"Sigma-Error"<<"\t"<<"Resolution (%)"<<"\n";
	for(i=1; i<384; i++)
	{if(i >64 && i <320){continue;}
	f<<i<<"\t"<<S[i]<<"\t"<<Sigma[i]<<"\t"<<SigErr[i]<<"\t"<<resolution[i]<<"\n";
	} 
	f.close();//closing fitting parameter file 
	
 	cout<<"Done.\n\n";
	 
}

/****************************************************************************************/


void calibrate_channels(const char* fname, const char* oname)
{
    clock_t tic = clock();
	gROOT->cd();
    
    
    //Input file 
	TFile *input = new TFile(Form("%s",fname) , "READ");
	
	TTree* data = (TTree*) input->Get("data");
	data->SetBranchAddress("q", &q);
	entries = data->GetEntries();
    cout<<"Number of Entries: "<<entries<<"\n";

    
    initialize_plots();
    
    run_threads();
    
    write_fitting_file();
 
    
 
//Output 
  	TFile *outfile = new TFile(Form("%s.root",oname), "RECREATE");
  	TTree *tree = new TTree("data","data");           // Define TTree
	tree->Branch("q",&qcalib,"qcalib[1024]/F");       // Define Branch
		
	
	for(j=1;j<=entries;j++)
	{
		for(i=1;i<364;i++)
		{
			if(i >64 && i <320){continue;}
			data->GetEntry(j);
			if(q[i]>0)
			{
			qcalib[i]=q[i]*S[i];
			}
			//else{continue;}
        	
		}
    tree->Fill();
	}
	
	input->Close();
	tree->Write();
	outfile->Close();

   
  	clock_t toc = clock();

    printf("Elapsed: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);

    cout<<"bye bye...\n";
	
}