/*
uint64_t PCt;
  uint64_t Inter;
  UShort_t beam;
  UShort_t MWPC;
  uint32_t NCh;

  int32_t  chn;
  uint32_t amp;
*/

#include <cstdint>
#include <iostream>
#include "Math/GSLSimAnMinimizer.h"
#include "Math/Functor.h"
#include <vector>
#include <cmath>
#include <string>
#include <TSpectrum.h>
#include <ctime>
//#include "TObject.h

using namespace std;



struct simple_event {
  int32_t  chn;
  uint32_t amp;
};

struct event {
  uint64_t PCt;
  uint32_t amp;
  uint64_t time;
  vector <simple_event> simult;
  vector <simple_event> MWPC;
  vector <simple_event> side;
		
};
struct coef {
	Double_t A;
	Double_t B;
};

struct peak {
	Double_t xpeak;
	Double_t sigma;
};



vector <simple_event> find_front (UInt_t ampl[], Int_t ev_chan[], Int_t NCh) {
	vector<simple_event> front;
    for (auto i=0;i<NCh;i++) {
		if (ev_chan[i] >=1 && ev_chan[i] <= 48) {
		simple_event ev;
		ev.amp = ampl[i];
		ev.chn = ev_chan[i];

		front.push_back(ev);
		}
    }
    return front;
  };

 
vector <simple_event> find_back (UInt_t ampl[], Int_t ev_chan[], Int_t NCh) {
	vector<simple_event> back;
    for (auto i=0;i<NCh;i++) {
		if (ev_chan[i] >=-128 && ev_chan[i] <= -1) {
		simple_event ev;
		ev.amp = ampl[i];
		ev.chn = -ev_chan[i];
		if ((-ev_chan[i] >= 1)&&(-ev_chan[i] <= 16)) ev.chn = ev.chn+16;
		if ((-ev_chan[i] >= 17)&&(-ev_chan[i] <= 32)) ev.chn = ev.chn-16;
		back.push_back(ev);
		}
    }
    return back;
  };




vector <simple_event> find_side (UInt_t ampl[], Int_t ev_chan[], Int_t NCh) {
	vector<simple_event> side;
    for (auto i=0;i<NCh;i++) {
		if (ev_chan[i] >=103 && ev_chan[i] <= 108) {
		simple_event ev;
		ev.amp = ampl[i];
		ev.chn = ev_chan[i]-102;
		side.push_back(ev);
		}
    }
    return side;
  };



vector <simple_event> find_MWPC (UInt_t ampl[], Int_t ev_chan[], Int_t NCh) {
	vector<simple_event> MWPC;
    for (auto i=0;i<NCh;i++) {
		if (ev_chan[i] == 121 || ev_chan[i] == 122) {
		simple_event ev;
		ev.amp = ampl[i];
		ev.chn = ev_chan[i]-120;
		MWPC.push_back(ev);
		}
    }
    return MWPC;
  };

vector <simple_event> find_veto (UInt_t ampl[], Int_t ev_chan[], Int_t NCh) {
	vector<simple_event> veto;
    for (auto i=0;i<NCh;i++) {
		if (ev_chan[i] == 112) {
		simple_event ev;
		ev.amp = ampl[i];
		ev.chn = ev_chan[i]-111;
		veto.push_back(ev);
		}
    }
    return veto;
  };





bool check_MWPC (vector <simple_event> MWPC) {
	bool check = false;
	for (auto i = 0; i<MWPC.size(); i++) if ((MWPC[i].chn == 1 && MWPC[i].amp > 300) || (MWPC[i].chn == 2 && MWPC[i].amp > 1000)) check = true;  
	return check;
  };

bool check_veto (vector <simple_event> veto) {
	bool veto_check = false;
	for (auto i = 0; i<veto.size(); i++) if (veto[i].chn == 1 && veto[i].amp < 500) veto_check = true; 
	return veto_check;
  };
	
vector<Double_t> out_time(ULong64_t time) {
	Double_t day, hour, min, sec;
	Double_t start = (time - 1451595600)/(60.0*60*24);
	hour = 24*modf(start, &day);
	min = 60*modf(hour, &hour);
	sec = 60*modf(min, &min);
	//cout <<"Day " << day << ", time " << hour << ":" << min << ":" << sec << endl;
	vector<Double_t> date = {day,hour,min,sec};
	return date;
};



TH1F Smoothing(TH1F* hist, int smooth,  Int_t xmin, Int_t xmax) {
	
	Int_t i;
	Int_t binmin = hist->GetXaxis()->FindBin(xmin);
	Int_t binmax = hist->GetXaxis()->FindBin(xmax);
	Int_t nbins = binmax - binmin;

	string name = hist->GetName() ;
	name = name + "_smooth";
	
	Double_t source [nbins];
   	for (i = 0; i < nbins; i++) source[i]=hist->GetBinContent(binmin + i );

   	TSpectrum *sp = new TSpectrum();

   	sp->SmoothMarkov(source,nbins,smooth);

	TH1F d (name.c_str(),"Smoothed hist",nbins,xmin,xmax);
		
   	for (i = 0; i < nbins; i++) d.SetBinContent(i+1,source[i]);
		
	return d;
		
	
}


void  Background(TH1F* hist,  Int_t window, Int_t smooth, Int_t xmin, Int_t xmax) { 
   Int_t i;  
   Int_t nbins = hist->GetNbinsX();
   Double_t source [nbins];
   TH1F d ("d","",nbins,xmin,xmax);
   TSpectrum *s = new TSpectrum();

   for (i = 0; i < nbins; i++) source[i]=hist->GetBinContent(i + 1);

   s->Background(source,nbins,window,1,2,kFALSE, smooth,kFALSE);

   for (i = 0; i < nbins; i++) d.SetBinContent(i + 1,source[i]);

    hist->Add(&d,-1); 
}


void Deconvolution(TH1F* hist, Double_t sigma, Int_t iter_numb, Int_t xmin, Int_t xmax) { 
    Int_t i;
    Int_t nbins = hist->GetNbinsX();
    Double_t source[nbins];
    Double_t response[nbins];
  
    string func = "TMath::Gaus(x," + to_string(xmin+3*sigma) + "," + to_string(sigma) + ")"; 
    TF1 *gaus = new TF1("gaus",func.c_str(),xmin,xmax); 
    TH1F h1 ("h1","Impulse response function",nbins,xmin,xmax); 
    
    h1.FillRandom("gaus",100000); 
    for (i = 0; i < nbins; i++) source[i]=hist->GetBinContent(i + 1);
    for (i = 0; i < nbins; i++) response[i]=h1.GetBinContent(i + 1);
 
    TSpectrum *s = new TSpectrum();
    s->Deconvolution(source,response,nbins,iter_numb,1,1);

    for (i = 0; i < nbins; i++) hist->SetBinContent(i + 1,source[i]);
 }





vector<peak> find_peaks (TH1F* hist, Int_t N_peak, Double_t sigmaf, Double_t tresh,  Int_t xmin, Int_t xmax) { 
	
	Int_t nfound=1;
	Int_t np=100;
	TSpectrum *spec = new TSpectrum(np);
		while (nfound!=N_peak) {		
			nfound = spec->Search(hist,sigmaf,"noMarkov nodraw",tresh);
			tresh = tresh+0.001;
			if (tresh>0.4) {
				//cout << hist->GetName() << "  !!!!!!!!!!!!!!!" << endl;
			 	break;
				}					
		}
	
     	Double_t * xpeaks = spec->GetPositionX();
	Double_t xsort[nfound];
	Int_t index[nfound];
     	TMath::Sort(nfound,xpeaks,index, kFALSE); 
     	for (int i=0; i <nfound; i++)  xsort[i] = xpeaks[index[i]];		


     	vector<peak> peaks(nfound);
			
     	for(int i=0;i<nfound;i++) peaks[i].xpeak = xsort[i];
	
	return peaks;	 	
}

vector<peak> fit_spectrum (vector<peak> xpeaks, TH1F*hist, Int_t xmin, Int_t xmax, Double_t par) { 
	hist->GetXaxis()->SetRangeUser(xmin, xmax);
	Int_t nfound = xpeaks.size();
     	TF1 **fp = new TF1*[nfound];
     	vector<peak> peaks(nfound);
	
			
     	for(int i=0;i<nfound;i++) {
      		 string  fname = hist->GetName();
      		 fname = fname + "fp"+ to_string(i);
		 fp[i] = new TF1(fname.c_str(),"gaus",xpeaks[i].xpeak-par, xpeaks[i].xpeak+par);
        	 hist->Fit(fname.c_str(),"QR+"); 
         	 peaks[i].xpeak= fp[i]->GetParameter(1);
        	 peaks[i].sigma = fp[i]->GetParameter(2);
     		}

	
	return peaks;
	
	
}

class peak_min {	
public:
	TH1F *hist = new TH1F();
	//Int_t N_peaks;
	TH1F smooth;
	vector<Double_t> diff;
	Int_t xmin, xmax;
	vector <peak> f_peaks;
 	double *pars;
	double operator()(const double *par) { 
	/*
	par0	Int_t smooth = 3;
	par1	Int_t window = 12;
	par2	Int_t smooth_back = 9;
	par3	Int_t IRF_sig = 2;
	par4	Int_t dec = 10;
	par5	Int_t sigma = 1;
	par6	Double_t tresh = 0.001;
	*/	
		TH1F *hist_smooth = new TH1F();
		vector <peak> f_peaks;
		
		*hist_smooth =  Smoothing(hist, par[0], xmin, xmax);
		Background(hist_smooth, par[1], par[2], xmin, xmax);
		Deconvolution(hist_smooth,  par[3], par[4]*100, xmin, xmax);
		f_peaks =  find_peaks (hist_smooth, diff.size() + 1,par[5], par[6]/100, xmin, xmax); 
	
		double square = 0;
		vector <Double_t> diff_est;
	
		if (f_peaks.size() < diff.size()+1) {
			square = 200;
			delete hist_smooth;
			
			return square;
		}
		for (Int_t j =0; j<f_peaks.size()-1; j++) {
			Double_t diff_j = f_peaks[j+1].xpeak - f_peaks[j].xpeak;
			diff_est.push_back(diff_j); 
		}
		
		for (Int_t j = 0; j<f_peaks.size()-1; j++) {
			square = square + (diff[j] - diff_est[j])*(diff[j] - diff_est[j]);
		}

		delete hist_smooth;
		//cout << square << endl;
		return square;
	}
	void analysis() {
		smooth =  Smoothing(hist, pars[0], xmin, xmax);
		Background(&smooth, pars[1], pars[2], xmin, xmax);
		Deconvolution(&smooth,  pars[3], pars[4], xmin, xmax);
		f_peaks =  find_peaks (&smooth, diff.size() + 1,pars[5], pars[6], xmin, xmax); 
	}

};



uint64_t time_start = 1481359608;
//uint64_t time_start = 0;
uint64_t time1 = 1481382780;
uint64_t time2 = 1481418960;
uint64_t time3 = 1481432340;
uint64_t time_finish = 1481436708;


void read(vector<string> name) { 
  ofstream out;
  out.open("log_read.txt");
  
  const Int_t NCh_max = 1000;
  ULong64_t PCt, Inter, empty;
  UShort_t beam, MWPC; 
  UInt_t NCh, ampl_1;
  UInt_t ampl[NCh_max] = {0};
  Int_t ev_chan[NCh_max] = {0};
  Int_t  ev_chan_1;
  Int_t Nch_c= 0;
  
  TTree all("all", "All file");
  all.Branch("PCt",&PCt,"PCt/l");
  all.Branch("Inter",&Inter,"Inter/l");
  all.Branch("beam",&beam,"beam/s");
  all.Branch("MWPC",&MWPC,"MWPC/s");
  all.Branch("NCh",&NCh,"NCh/i");
  //all.Branch("veto",&check_veto,"veto/O");
  all.Branch("ev_chan", ev_chan, "ev_chan[NCh]/I");
  all.Branch("ampl", ampl, "ampl[NCh]/i");
 
  Int_t c = 0;
 for (Int_t i = 0; i< name.size(); i++) {
   ifstream is1(name[i], ios::in | ios::binary);
      while(!is1.eof()) { 
		
        	is1.read((char*)&PCt,sizeof(PCt));
        	is1.read((char*)&empty,sizeof(empty));
        	is1.read((char*)&Inter,sizeof(Inter));
      	   	is1.read((char*)&beam,sizeof(beam));
      		is1.read((char*)&MWPC,sizeof(MWPC));
      		is1.read((char*)&NCh,sizeof(NCh));
		if(NCh>Nch_c) Nch_c = NCh;	
      		for(Int_t j = 0;j<NCh;j++) {
	  		is1.read((char*)&ev_chan_1,sizeof(ev_chan_1));
	  		is1.read((char*)&ampl_1,sizeof(ampl_1));
	  		ev_chan[j] = ev_chan_1;
	  		ampl[j] = ampl_1;
    		  }
      	      
     if (PCt > time_start) all.Fill();
      	
      
       
      if ( PCt>time_start && NCh > 100) { 
		vector<Double_t> data_time = out_time(PCt);
		out << name[i] << endl;
		out << "Day: " << data_time[0] << " Time: " << data_time[1] << ":" << data_time[2] << ":" << data_time[3] << endl;
		out << "Inter " << Inter << " MWPC " << MWPC << " NCh " << NCh << endl;
		for (Int_t j =0; j<NCh; j++) out << ev_chan[j] << "  " << ampl[j] << endl;
		out << endl;
		c++;
		
      }
	 for (Int_t j = 0; j<NCh_max; j++) {
			ev_chan[j] = 0;
	  		ampl[j] = 0;
		}
}
  is1.close();
  out << c << endl;
  Int_t cc = all.GetEntries();
  cout << c << " " << cc << " " << Nch_c << endl;
}
cout << "NCH MAX " << Nch_c << endl;  

TFile f("output.root","recreate");
all.Write();
f.Close();
out.close();
}

void dist() {
	ofstream out;
	out.open("log_dist.txt");
	TFile *f = new TFile("output.root","update");
	TTree *all =0;
	f->GetObject("all",all);

	const Int_t NCh_max = 1000;
	ULong64_t PCt, Inter, empty;
	UShort_t beam, MWPC; 
  	UInt_t NCh, ampl_1;
  	UInt_t ampl[NCh_max] = {0}, ampl_back[NCh_max] = {0}, ampl_MWPC[NCh_max] = {0}, ampl_side[NCh_max] = {0}, ampl_front[NCh_max] = {0};
 	Int_t ev_chan[NCh_max] = {0}, ev_chan_back[NCh_max] = {0}, ev_chan_MWPC[NCh_max] = {0}, ev_chan_side[NCh_max] = {0}, ev_chan_front[NCh_max] = {0};
	Bool_t MWPC_check, veto_check;
	//Int_t back_ampl[NCh_max] = {0};
 	//Int_t back_chan[NCh_max] = {0};

 	Int_t ev_chan_1, N_back, N_MWPC, N_side, N_front, N_veto, cc;
	Int_t c = 0;
	all->SetBranchAddress("PCt",&PCt);
	all->SetBranchAddress("Inter",&Inter);
	all->SetBranchAddress("beam",&beam);
	all->SetBranchAddress("MWPC",&MWPC);
	all->SetBranchAddress("NCh",&NCh);
	all->SetBranchAddress("ev_chan",ev_chan);
	all->SetBranchAddress("ampl",ampl);

	ULong64_t nn = all->GetEntries();
	
	vector <simple_event> back_f, side_f, front_f, MWPC_f, veto_f;


	TTree *front_all = new TTree("front_all", "Front events");
	TTree *back_all = new TTree("back_all", "Back events");
	TTree *side_all = new TTree("side_all", "Side events");
	TTree *MWPC_all = new TTree("MWPC_all", "MWPC events");
	TTree *veto_all = new TTree("veto_all", "Veto events");

	front_all->SetAutoFlush(0);
	back_all->SetAutoFlush(0);
	side_all->SetAutoFlush(0);
	MWPC_all->SetAutoFlush(0);
	veto_all->SetAutoFlush(0);

	//all.Branch("PCt",&PCt,"PCt/l");
  	front_all->Branch("Inter",&Inter,"Inter/l");
  	front_all->Branch("beam",&beam,"beam/s");
  	front_all->Branch("MWPC",&MWPC_check,"MWPC/O");
  	front_all->Branch("NCh",&NCh,"NCh/i");

 	front_all->Branch("ev_chan", &ev_chan_1, "ev_chan/I");
 	front_all->Branch("ampl", &ampl_1, "ampl/i");

	front_all->Branch("N_back",&N_back, "N_back/I");
	front_all->Branch("ev_chan_back", ev_chan_back, "ev_chan_back[N_back]/I");
 	front_all->Branch("ampl_back", ampl_back, "ampl_back[N_back]/i");

	front_all->Branch("N_side",&N_side, "N_side/I");
	front_all->Branch("ev_chan_side", ev_chan_side, "ev_chan_side[N_side]/I");
 	front_all->Branch("ampl_side", ampl_side, "ampl_side[N_side]/i");

	front_all->Branch("N_MWPC",&N_MWPC, "N_MWPC/I");
	front_all->Branch("ev_chan_MWPC", ev_chan_MWPC, "ev_chan_MWPC[N_MWPC]/I");
 	front_all->Branch("ampl_MWPC", ampl_MWPC, "ampl_MWPC[N_MWPC]/i");



	back_all->Branch("Inter",&Inter,"Inter/l");
  	back_all->Branch("beam",&beam,"beam/s");
  	back_all->Branch("MWPC",&MWPC_check,"MWPC/O");
  	back_all->Branch("NCh",&NCh,"NCh/i");

 	back_all->Branch("ev_chan", &ev_chan_1, "ev_chan/I");
 	back_all->Branch("ampl", &ampl_1, "ampl/i");

	back_all->Branch("N_front",&N_front, "N_front/I");
	back_all->Branch("ev_chan_front", ev_chan_front, "ev_chan_front[N_front]/I");
 	back_all->Branch("ampl_front", ampl_front, "ampl_front[N_front]/i");

	back_all->Branch("N_side",&N_side, "N_side/I");
	back_all->Branch("ev_chan_side", ev_chan_side, "ev_chan_side[N_side]/I");
 	back_all->Branch("ampl_side", ampl_side, "ampl_side[N_side]/i");

	back_all->Branch("N_MWPC",&N_MWPC, "N_MWPC/I");
	back_all->Branch("ev_chan_MWPC", ev_chan_MWPC, "ev_chan_MWPC[N_MWPC]/I");
 	back_all->Branch("ampl_MWPC", ampl_MWPC, "ampl_MWPC[N_MWPC]/i");


		
	side_all->Branch("Inter",&Inter,"Inter/l");
  	side_all->Branch("beam",&beam,"beam/s");
  	side_all->Branch("MWPC",&MWPC_check,"MWPC/O");
 	side_all->Branch("ev_chan", &ev_chan_1, "ev_chan/I");
 	side_all->Branch("ampl", &ampl_1, "ampl/i");


	
	MWPC_all->Branch("Inter",&Inter,"Inter/l");
  	MWPC_all->Branch("beam",&beam,"beam/s");
 	MWPC_all->Branch("ev_chan", &ev_chan_1, "ev_chan/I");
 	MWPC_all->Branch("ampl", &ampl_1, "ampl/i");



	veto_all->Branch("Inter",&Inter,"Inter/l");
  	veto_all->Branch("beam",&beam,"beam/s");
 	veto_all->Branch("ev_chan", &ev_chan_1, "ev_chan/I");
 	veto_all->Branch("ampl", &ampl_1, "ampl/i");
	
	
	cout << nn << endl;
	for( ULong64_t i = 0; i<nn; i++) {
		all->GetEntry(i);
		front_f = find_front (ampl, ev_chan, NCh);
		back_f = find_back (ampl, ev_chan, NCh);
		side_f = find_side (ampl, ev_chan, NCh);
		MWPC_f = find_MWPC (ampl, ev_chan, NCh);
		veto_f = find_veto (ampl, ev_chan, NCh);
		N_back = back_f.size();
		N_side = side_f.size();
		N_front = front_f.size();
		N_MWPC = MWPC_f.size();
		N_veto = veto_f.size();
		if (N_MWPC != MWPC) {
			vector<Double_t> data_time = out_time(PCt);
			out << "Day: " << data_time[0] << " Time: " << data_time[1] << ":" << data_time[2] << ":" << data_time[3] << endl;
			out << "Inter " << Inter << " MWPC " << MWPC << " NCh " << NCh << endl;
			for (Int_t j =0; j<NCh; j++) out << ev_chan[j] << "  " << ampl[j] << endl;
			out << endl;
			c++;
     	 		}
			
		MWPC_check = check_MWPC(MWPC_f);
		if (N_MWPC == 0) MWPC_check = false;
	
		for(Int_t ii = 0; ii < N_back; ii++) {
			ampl_back[ii] = back_f[ii].amp;
			ev_chan_back[ii] = back_f[ii].chn; 
			}
	
		for(Int_t ii = 0; ii < N_side; ii++) {
			ampl_side[ii] = side_f[ii].amp;
			ev_chan_side[ii] = side_f[ii].chn; 
			}
		for(auto ii = 0; ii < N_MWPC; ii++) {
			ampl_MWPC[ii] = MWPC_f[ii].amp;
			ev_chan_MWPC[ii] = MWPC_f[ii].chn; 
			}
		for(Int_t ii = 0; ii < N_front; ii++) {	
			ampl_front[ii] = front_f[ii].amp;
			ev_chan_front[ii] = front_f[ii].chn; 
			}
		
		MWPC_check = check_MWPC(MWPC_f);
		veto_check = check_veto(veto_f);

		
		if (!veto_check) {
			for(Int_t ii = 0; ii < N_front; ii++) {	
				ev_chan_1 = ev_chan_front[ii];
				ampl_1 = ampl_front[ii];	
				front_all->Fill();
			}	
			for(Int_t ii = 0; ii < N_back; ii++) {
				ev_chan_1 = ev_chan_back[ii];
				ampl_1 = ampl_back[ii];
				back_all->Fill();
			}
			for(Int_t ii = 0; ii < N_MWPC; ii++) {
				ev_chan_1 = ev_chan_MWPC[ii];
				ampl_1 = ampl_MWPC[ii];
				MWPC_all->Fill();
			}
			for(Int_t ii = 0; ii < N_side; ii++) {
				ev_chan_1 = ev_chan_side[ii];
				ampl_1 = ampl_side[ii];
				side_all->Fill();
			}
		
		}

		for(Int_t ii = 0; ii < NCh; ii++) {
			if (ev_chan[ii] == 112) {
				ev_chan_1 = ev_chan[ii];
				ampl_1 = ampl[ii];
				veto_all->Fill();
			}
		}
	}
/*
		for (Int_t j = 0; j<NCh;j++) {
			if (ev_chan[j] >=1 && ev_chan[j] <= 48 && !veto_check) {
				ev_chan_1 = ev_chan[j];
				ampl_1 = ampl[j];
				front_all->Fill();
			}
			if (ev_chan[j] >=-128 && ev_chan[j] <= -1 && !veto_check) {
				ev_chan_1 = -ev_chan[j];
				ampl_1 = ampl[j];
				if ((-ev_chan_1 >= 1)&&(-ev_chan_1 <= 16)) ev_chan_1 = ev_chan_1+16;
				if ((-ev_chan_1 >= 17)&&(-ev_chan_1 <= 32)) ev_chan_1 = ev_chan_1-16;
				back_all->Fill();
	   	 	}	
		}
*/
					

	out << cc << "   " << c << endl;
	
	front_all->Write();
	back_all->Write();
	side_all->Write();
	MWPC_all->Write();
	veto_all->Write();
	f->ls();
	f->Close();	
}

void hist() {
	TFile *f = new TFile("output.root","update");

	TTree *front_all = 0;
	TTree *back_all = 0;
	TTree *side_all = 0;
	TTree *MWPC_all = 0;
	TTree *veto_all = 0;
	

	f->GetObject("front_all",front_all);
	f->GetObject("back_all",back_all);
	f->GetObject("side_all", side_all);
	f->GetObject("MWPC_all", MWPC_all);
	f->GetObject("veto_all", veto_all);

	Int_t xmin = 0;
	Int_t xmax = 32767;
	Int_t nbins = xmax;
	f->mkdir("Front hist");
	f->cd("Front hist");
	TH1F **front_hist = new TH1F*[48];
	TH1F **front_hist_MWPC = new TH1F*[48];
	TH1F **front_hist_no_MWPC = new TH1F*[48];
	for (Int_t i = 0; i<48;i++) {
		string name = "front_" + to_string(i);
		front_hist[i] = new  TH1F(name.c_str(),"Front spectrum",nbins,xmin,xmax); 
		name = "front_MWPC_" + to_string(i);
		front_hist_MWPC[i] = new TH1F(name.c_str(),"Front spectrum MWPC",nbins,xmin,xmax);
		name = "front_no_MWPC_"+to_string(i);
		front_hist_no_MWPC[i] = new TH1F(name.c_str(),"Front spectrum no MWPC",nbins,xmin,xmax);
	}
	
	f->cd();
	f->mkdir("Back hist");
	f->cd("Back hist");
	TH1F **back_hist = new TH1F*[128];
	TH1F **back_hist_MWPC = new TH1F*[128];
	TH1F **back_hist_no_MWPC = new TH1F*[128];
	for (Int_t i = 0; i<128;i++) {
		string name = "back_" + to_string(i);
		back_hist[i] = new  TH1F(name.c_str(),"Back spectrum",(nbins-3)/4,xmin,xmax-3); 
		name = "back_MWPC_" + to_string(i);
		back_hist_MWPC[i] = new TH1F(name.c_str(),"Back spectrum MWPC",(nbins-3)/4,xmin,xmax-3);
		name = "back_no_MWPC_"+to_string(i);
		back_hist_no_MWPC[i] = new TH1F(name.c_str(),"Back spectrum no MWPC",(nbins-3)/4,xmin,xmax-3);
	}

	
	f->cd();
	f->mkdir("Side hist");
	f->cd("Side hist");
	TH1F **side_hist = new TH1F*[6];
	TH1F **side_hist_MWPC = new TH1F*[6];
	TH1F **side_hist_no_MWPC = new TH1F*[6];
	for (Int_t i = 0; i<6; i++) {
		string name = "side_" + to_string(i);
		side_hist[i] = new TH1F(name.c_str(),"Side spectrum",(nbins-1)/2,xmin,xmax-1);
		name = "side_MWPC_" + to_string(i);
		side_hist_MWPC[i] = new TH1F(name.c_str(),"Side spectrum MWPC",(nbins-1)/2,xmin,xmax-1);
		name = "side_no_MWPC_"+to_string(i);
		side_hist_no_MWPC[i] = new TH1F(name.c_str(),"Side spectrum no MWPC",(nbins-1)/2,xmin,xmax-1);
	}
	
	f->cd();
	TH1F **MWPC_hist = new TH1F*[2];
	for (Int_t i = 0; i<2; i++) {
		string name = "MWPC_" + to_string(i);
		MWPC_hist[i] = new TH1F(name.c_str(),"MWPC spectrum",nbins,xmin,xmax);
	}		


	TH1F *veto_hist = new TH1F("veto_hist","Veto spectrum",nbins,xmax,xmax);

	const Int_t NCh_max = 1000;
	ULong64_t PCt, Inter, empty;
	UShort_t beam, MWPC; 
  	UInt_t NCh, ampl_1;
  	
	Bool_t MWPC_check, veto_check;
	Int_t ev_chan_1, N_back, N_MWPC, N_side, N_front, N_veto;
	
	front_all->SetBranchAddress("MWPC",&MWPC_check);
 	front_all->SetBranchAddress("ev_chan", &ev_chan_1);
 	front_all->SetBranchAddress("ampl", &ampl_1);

	back_all->SetBranchAddress("MWPC",&MWPC_check);
 	back_all->SetBranchAddress("ev_chan", &ev_chan_1);
 	back_all->SetBranchAddress("ampl", &ampl_1);

	side_all->SetBranchAddress("MWPC",&MWPC_check);
 	side_all->SetBranchAddress("ev_chan", &ev_chan_1);
 	side_all->SetBranchAddress("ampl", &ampl_1);

 	MWPC_all->SetBranchAddress("ev_chan", &ev_chan_1);
 	MWPC_all->SetBranchAddress("ampl", &ampl_1);

 	veto_all->SetBranchAddress("ev_chan", &ev_chan_1);
 	veto_all->SetBranchAddress("ampl", &ampl_1);
	
	N_front = front_all->GetEntries();
	N_back = back_all->GetEntries();
	N_side = side_all->GetEntries();
	N_MWPC = MWPC_all->GetEntries();
	N_veto = veto_all->GetEntries();

	for(Int_t i=0; i<N_front;i++) {
		front_all->GetEntry(i);
		front_hist[ev_chan_1-1]->Fill(ampl_1);
		if (MWPC_check) { 
			front_hist_MWPC[ev_chan_1-1]->Fill(ampl_1); 
		}  else { 
			front_hist_no_MWPC[ev_chan_1-1]->Fill(ampl_1); 
		}
	}
	front_hist_no_MWPC[5]->Draw();
	for (Int_t i = 0;i<N_back;i++) {
		back_all->GetEntry(i);
		if (ampl_1 > xmax-3) ampl_1 = xmax-3;
		back_hist[ev_chan_1-1]->Fill(ampl_1);
		if (MWPC_check) {
			back_hist_MWPC[ev_chan_1-1]->Fill(ampl_1);
		} else {
			back_hist_no_MWPC[ev_chan_1-1]->Fill(ampl_1);
		}
	}

	for (Int_t i=0; i<N_side;i++) {
		side_all->GetEntry(i);
		if (ampl_1 > xmax-1) ampl_1 = xmax-1;
		side_hist[ev_chan_1-1]->Fill(ampl_1);
		if (MWPC_check) {
			side_hist_MWPC[ev_chan_1-1]->Fill(ampl_1);
		} else { 
			side_hist_no_MWPC[ev_chan_1-1]->Fill(ampl_1);
		}
	}

	for (Int_t i=0; i<N_MWPC;i++) {
		MWPC_all->GetEntry(i);
		MWPC_hist[ev_chan_1-1]->Fill(ampl_1);
	}

	for (Int_t i=0; i<N_veto;i++) {
		veto_all->GetEntry(i);
		veto_hist->Fill(ampl_1);
	}

	
	f->Write();
	f->Close();				
}

void NumMin(peak_min peak) {

	ROOT::Math::GSLSimAnMinimizer min;
	ROOT::Math::Functor funct(peak,7);

	double step[7] = {1, 3, 3,.5, 1,.5, 0.2};
	double variable[7] = {5, 3, 6, 2, 1, 2, 0.1};
	//cout << "test0" << endl;
	min.SetFunction(funct);
	//min->SetMaxFunctionCalls(100000);
	min.SetMaxIterations(10);
	min.SetTolerance(0.1);

	min.SetVariable(0,"Smooth",variable[0], step[0]);
	min.SetVariable(1,"Window",variable[1], step[1]);
	min.SetVariable(2,"Smooth back",variable[2], step[2]);
	min.SetVariable(3,"IRF_sig",variable[3], step[3]);
	min.SetVariable(4,"Decon iter",variable[4], step[4]);
	min.SetVariable(5,"Sigma",variable[5], step[5]);
	min.SetVariable(6,"Tresh",variable[6], step[6]);

	min.SetVariableLimits(0,0,12);
	min.SetVariableLimits(1,0,18);
	min.SetVariableLimits(2,0,18);
	min.SetVariableLimits(3,1,3);
	min.SetVariableLimits(4,0,20);
	min.SetVariableLimits(5,1,5);
	min.SetVariableLimits(6,0.1,20);

	min.SetPrintLevel(0);
	cout << "test00" << endl;
	min.Minimize(); 
	cout << "test01" << endl;
	const double *xs  = min.X();
  	cout << "test02" << endl;
	for(int i=0; i<7; i++){
		cout << "test03" << endl;
		peak.pars[i] = xs[i];
		if(i==4) peak.pars[i]=xs[i]*100;
		if(i==6) peak.pars[i] = xs[i]/100;
		cout << peak.pars[i] << endl;
	}
	cout << "test04" << endl;
}

void find_peak () {
	TFile f ("output.root");
	
	TH1F **front_hist_no_MWPC = new TH1F*[48];
	
	
	f.cd("Front hist");

	for (Int_t i = 0; i<48;i++) {
		string name = "front_no_MWPC_"+to_string(i) + ";1";
		gDirectory->GetObject(name.c_str(),front_hist_no_MWPC[i]);
		front_hist_no_MWPC[i]->SetDirectory(0);
	}

	f.cd();
	f.cd("Back hist");
	TH1F **back_hist_no_MWPC = new TH1F*[128];
	for (Int_t i = 0; i<128;i++) {
		string name = "back_no_MWPC_"+to_string(i) + ";1";
		gDirectory->GetObject(name.c_str(),back_hist_no_MWPC[i]);
		back_hist_no_MWPC[i]->SetDirectory(0);
	}
	f.Close();
	TFile ff ("smooth.root","update");
	TH1F **front_smooth = new TH1F*[48];

	
	/*
	TCanvas *c = new TCanvas("c","",200,10,700,1000);
  	TPad *pad1 = new TPad("pad1","",0.05,0.52,0.95,0.97);
 	TPad *pad2 = new TPad("pad2","",0.05,0.02,0.95,0.47);
  
 	pad1->Draw();
 	pad2->Draw();
  	pad1->cd();

	front_hist_no_MWPC[22]->GetXaxis()->SetRangeUser(300,550);
	front_hist_no_MWPC[22]->Draw();
	*/
	
	

	
	vector< vector<peak> > front_peaks(48,vector<peak>());

	Int_t xmin = 280;
	Int_t xmax = 550;
	Int_t smooth = 3;
	Int_t smooth_back = 9;
	Int_t window = 12;
	Int_t dec = 10;
	Int_t IRF_sig = 2;
	Double_t tresh = 0.001;
	Int_t N_peaks = 11;
	Int_t sigma_f = 1;
	Double_t par = 2.5;
		
	front_smooth[0] = new TH1F();
	*front_smooth[0] =  Smoothing(front_hist_no_MWPC[0], smooth, xmin, xmax);
	Background(front_smooth[0],  window, smooth_back, xmin, xmax);
	Deconvolution(front_smooth[0],  IRF_sig, dec, xmin, xmax);
	front_peaks[0] =  find_peaks (front_smooth[0], N_peaks,sigma_f, tresh, xmin, xmax); 
	string str = to_string(tresh);
	str.erase ( str.find_last_not_of('0') + 1, std::string::npos );
	string name = "front_" + to_string(0) + "_sm_" + to_string(smooth) + "_b_" + to_string(window) + "_" + to_string(smooth_back) + "_dec_" + to_string(IRF_sig) + "_" + to_string(dec) + "_pf_" + to_string(N_peaks) + "_" + str + "_" + to_string(sigma_f);
	front_smooth[0]->SetName(name.c_str());
	front_peaks[0] = fit_spectrum(front_peaks[0],front_smooth[0],xmin,xmax, par);
	
	front_hist_no_MWPC[0]->GetXaxis()->SetRangeUser(xmin,xmax);
	front_hist_no_MWPC[0]->Write();
	front_smooth[0]->Write();
	//vector<peak> fit_spectrum (vector<peak> xpeaks, TH1F*hist, Int_t xmin, Int_t xmax, Double_t par) 
	

	vector <Double_t> diffrn;
	for (Int_t j =0; j<front_peaks[0].size() - 1; j++) {
		Double_t diff_j = front_peaks[0][j+1].xpeak - front_peaks[0][j].xpeak;
		diffrn.push_back(diff_j); 
	}

	peak_min peak;
		
	peak.diff = diffrn;
	peak.xmin = xmin;
	peak.xmax = xmax;
		peak.hist = front_hist_no_MWPC[1];
		NumMin (peak);
		cout << "test1" << endl;
		peak.analysis();
		cout << "test2" << endl;
		*front_smooth[1] = peak.smooth;
		cout << "test3" << endl;
		front_smooth[1]->Write();
		cout << "test4" << endl;
	
		
}
		
		/*
class peak_min {	
public:
	TH1F *hist = new TH1F();
	//Int_t N_peaks;
	TH1F *smooth = new TH1F();
	vector<Double_t> diff;
	Int_t xmin, xmax;
 	double pars[7];
	vector <peak> f_peaks;
	double operator()(const double *par) { 


		TH1F *hist_smooth = new TH1F();
		vector <peak> f_peaks;
		
		*hist_smooth =  Smoothing(hist, par[0], xmin, xmax);
		Background(hist_smooth, par[1], par[2], xmin, xmax);
		Deconvolution(hist_smooth,  par[3], par[4], xmin, xmax);
		f_peaks =  find_peaks (hist_smooth, diff.size() + 1,par[5], par[6], xmin, xmax); 
	
		double square = 0;
		vector <Double_t> diff_est;
	
		if (f_peaks.size() < diff.size()+1) {
			square = 200;
			delete hist_smooth;
			cout << square << endl;
			return square;
		}
		for (Int_t j =0; j<f_peaks.size()-1; j++) {
			Double_t diff_j = f_peaks[j+1].xpeak - f_peaks[j].xpeak;
			diff_est.push_back(diff_j); 
		}
		
		for (Int_t j = 0; j<f_peaks.size()-1; j++) {
			square = square + (diff[j] - diff_est[j])*(diff[j] - diff_est[j]);
		}
		pars = par;
		cout << square << endl;
		delete hist_smooth;
		return square;
	}
	TH1F anal() {
		TH1F hist_smooth;
	
		
		hist_smooth =  Smoothing(hist, pars[0], xmin, xmax);
		Background(hist_smooth, pars[1], pars[2], xmin, xmax);
		Deconvolution(hist_smooth,  pars[3], pars[4], xmin, xmax);
		f_peaks =  find_peaks (hist_smooth, diff.size() + 1,pars[5], pars[6], xmin, xmax); 
		return hist_smooth
	}
};	
};
 
*/


	
	/*
	
	for (Int_t i = 1; i<48; i++) {
		
		front_smooth[i] = new TH1F();
		*front_smooth[i] =  Smoothing(front_hist_no_MWPC[i], smooth, xmin, xmax);
		Background(front_smooth[i],  window, smooth_back, xmin, xmax);
		Deconvolution(front_smooth[i],  IRF_sig, dec, xmin, xmax);
		front_peaks[i] =  find_peaks (front_smooth[i], N_peaks,sigma_f, tresh, xmin, xmax); 

		str = to_string(tresh);
		str.erase ( str.find_last_not_of('0') + 1, std::string::npos );
		name = "front_" + to_string(i) + "_sm_" + to_string(smooth) + "_b_" + to_string(window) + "_" + to_string(smooth_back) + "_dec_" + to_string(IRF_sig) + "_" + to_string(dec) + "_pf_" + to_string(N_peaks) + "_" + str + "_" + to_string(sigma_f);
		front_smooth[i]->SetName(name.c_str());
		

		//Double_t diffr  = front_peaks[0][N_peak-1].xpeak - front_peaks[i][N_peak-1].xpeak;
		for (Int_t j =1; j<front_peaks[i].size(); j++) front_peaks[i][N_peaks-1-j].xpeak = front_peaks[i][N_peaks-j].xpeak - diffrn[N_peaks -1 - j]; 
		front_peaks[i] = fit_spectrum(front_peaks[i],front_smooth[i],xmin,xmax, par);
		
		front_hist_no_MWPC[i]->GetXaxis()->SetRangeUser(xmin,xmax);
		front_hist_no_MWPC[i]->Write();	
		front_smooth[i]->Write();
	
	}
	*/
	/*
	vector < vector <Double_t>> diff(48,vector<Double_t>());
	for (Int_t i = 0; i<48; i++) {
		for (Int_t j =0; j<front_peaks[i].size() - 1; j++) {
			Double_t diff_j = front_peaks[i][j+1].xpeak - front_peaks[i][j].xpeak;
			diff[i].push_back(diff_j); 
		}
	}
	vector <Double_t> square;
	square.push_back(0);
	for (Int_t i = 0; i<47; i++) {
		Double_t diff_j = 0;
		for (Int_t j =0; j<front_peaks[i].size() - 1; j++) {
			diff_j = diff_j + (diff[0][j] - diff[i+1][j])*(diff[0][j] - diff[i+1][j]);
		}
		square.push_back(diff_j);
	}
	
	Double_t mean[diff[1].size()];
	Double_t sigma[diff[1].size()];
	Double_t sigma_0[diff[1].size()];
	
	for (Int_t j = 0; j<diff[1].size(); j++) {
		mean[j] =0;
		sigma[j] = 0;
		sigma_0[j] =0;
		for (Int_t i =0; i<48; i++) {
			mean[j] = mean[j] + diff[i][j];
		}
	}
	
	ofstream out;
	out.open("log_find.txt");
	for (Int_t j =0; j<diff[1].size(); j++) mean[j] = mean[j]/48;


	for (Int_t j = 0; j<diff[1].size(); j++) {
		for (Int_t i =0; i<48; i++) {
			sigma[j] = sigma[j] + (diff[i][j] - mean[j])*(diff[i][j] - mean[j]);
		}
	}

	for (Int_t j = 0; j<diff[1].size(); j++) {
		for (Int_t i =1; i<48; i++) {
			sigma_0[j] = sigma_0[j] + (diff[i][j] - diff[0][j])*(diff[i][j] - diff[0][j]);
			
		}
		
	}

	for (Int_t j =0; j<diff[1].size(); j++) {
		sigma[j] = sigma[j]/48;
		sigma_0[j] = sigma_0[j]/47;
		sigma[j] = sqrt(sigma[j]);
		sigma_0[j] = sqrt(sigma_0[j]);
	}


	

	for (Int_t j =0; j<diff[1].size(); j++) out << mean[j] << "  ";	
	out << endl;	
	for (Int_t j =0; j<diff[1].size(); j++) out << sigma[j] << "  ";
	out << endl;
	for (Int_t j =0; j<diff[1].size(); j++) out << sigma_0[j] << "  ";
	out << endl;	
	for (Int_t i = 0; i<48; i++) {
		out << "Strip Number " << i << endl;
		for (Int_t j =0; j<diff[i].size(); j++) out << diff[i][j] << "  ";
		out << "Diff sq " << square[i];
		out << endl;
		out << endl;
	}
	double *v = new double;
	*v = 10;
	*/
	//testf fs(4.5);
	//fs.x = 4.5;
	//double a = fs(v);
	//cout << a << " "<< fs.xx() << fs.yy() << endl;

	/*
	TH1F *diffs = new TH1F("name.c_str()","Diffs",30,0,30);;
	for (Int_t j =0; j<47; j++) diffs->Fill(square[j+1]);
	diffs->Write();
	*/
	
	
	

//}