#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 <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 "TROOT.h"

  

   

void calibrate1(const char* fname, const char* oname)//, const float peakcutoff=10)
{	
    
   //ROOT::EnableImplicitMT(8);  
	
    clock_t tic = clock();
	gROOT->cd();
	Double_t    fPositionX[10];
    Double_t    fPositionY[10];
    //Double_t    Energy[10], counts[1024];
    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,CHID;
	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 peakcutoff=10;   //Array on Right
    int   sradius = 2;      //radius for identified peak fitting mean +- sradius
    
	
	
	TFile *input = new TFile(Form("%s",fname), "READ");
	
	TTree *data = (TTree*) input->Get("data");
	data->SetBranchAddress("q", &q);
	int entries = data->GetEntries();
   //	entries = entries/2;
	cout<<entries<<"\n";

	
	TH1F* h1 = new TH1F("h1", "Spectrum of Gamma Rays", nbins, 0, 100);


	for(CHID = 1;CHID<384;CHID++)
	{
        cout<<"#####################################################################"<<"\n";
		
        ChannelID[CHID]=CHID;
		cout<<"Performing fit for Channel ID ="<<CHID<<"\n"; 		
		for(j=1;j<=entries;j++)
		{
            data->GetEntry(j);
			
			if(q[CHID]>0)
			{
                h1->Fill(q[CHID]);
			}
			else{continue;}
        }	
	
		
	//h1->Draw("L");

  int nfound = 0;
     
 	//Searching for Peaks in the spectrum
	TSpectrum *s = new TSpectrum();
  
  	for (i = 0; i < nbins; i++) {source[i]=h1->GetBinContent(i+1);}
   	nfound = s->SearchHighRes(source, dest, nbins, 8, 1, kFALSE, 3, kTRUE, 3);
   //nfound = s->Search(h1,1,"",0.2);

	printf("Found %d total Peaks\n",nfound);
	
	//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);
        fPositionX[i] = h1->GetBinCenter(bin);
     	fPositionY[i] = h1->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(peakBinNb < peakcutoff)
        {
        //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 -sradius, peakBinNb + sradius);
	auto fitRes = h1->Fit("gausini","RS");
	auto gauspol = new TF1("gauspol", "gaus(0)+pol1(3)", peakBinNb - sradius, peakBinNb + sradius);
	gauspol->SetParameter(0,fitRes->Parameter(0)); 
	gauspol->SetParameter(1,fitRes->Parameter(1));
	gauspol->SetParameter(2,fitRes->Parameter(2));
	auto fitres = h1->Fit("gauspol","RS");
    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");
   	
	//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;
   	
cout<<"#####################################################################"<<"\n";
		
	for (i = 0; i < nfound; i++) //Resetting values  
	{
     	 fPositionX[i] = 0;
     	 fPositionY[i] = 0;
	}
    
	} //closing CHID main loop over all the channels
	


//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(CHID=1;CHID<65;CHID++)
	{if(CHID >64 && CHID <320){continue;}
	f<<CHID<<"\t"<<S[CHID]<<"\t"<<Sigma[CHID]<<"\t"<<SigErr[CHID]<<"\t"<<resolution[CHID]<<"\n";
	} 
	f.close();//closing fitting parameter file 
	
 	cout<<"Done.\n\n";
	

	cout<<"Writing root file.\n\n";
    
//Output 
  	TFile *outfile = new TFile(Form("%s.root",oname), "RECREATE");
  
  	  
  	// Define TTree
  	TTree *tree = new TTree("data","data");
	
	tree->Branch("q",&qcalib,"qcalib[1024]/F");
		
	
	for(j=1;j<=entries;j++)
	{
		for(CHID=1;CHID<65;CHID++)
		{
			data->GetEntry(j);
			if(q[CHID]>0)
			{
			qcalib[CHID]=q[CHID]*S[CHID];
			}
			//else{continue;}
        	
		}
    tree->Fill();
	}
	
	

	//fclose(ptr_infile);
	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";
	
}