Double_t Gaussian(Double_t *mES,Double_t* par)
	{
	Double_t y, g;     //Gaussian parameters
	
		if(par[2] = 0)
			{
			y = (mES[0] - par[1])/(par[2]);
			g = par[0] *TMath::Exp(-0.5*y*y);
			}
			else
			{
			g = 0.0;
			}
	return g;
	}
Double_t Argus(Double_t*mES, Double_t* par)
	
	{
	Double_t x,u,u1,v,arg1,v1; // Argus parameters
	Double_t Ebeam = 5.29;	//same
	
	x = mES[0] / Ebeam;
	
	if (x < 1.0 && x > -1.0)
			{
			u = (1.0 -(1.0*(x*x)) );
			v = 1.0*(u-(1.0*(u*u)) );
			v1 = 1.0*sqrt(v);
			arg1 = par[1]*u;
			}
	Double_t fitty = par[0] *v1*exp(-1.0*arg1);
	
	return fitty;
	
	}

Double_t Total(Double_t* mES, Double_t* par) 

	{
	Double_t y, f;     //Gaussian parameters
	Double_t x,u,u1,v,arg1,v1; // Argus parameters
	Double_t Ebeam = 5.29;	//same
	
	x = mES[0] / Ebeam;
	
	if (x < 1.0 && x > -1.0)
			{
			u = (1.0 -(1.0*(x*x)) );
			v = 1.0*(u-(1.0*(u*u)) );
			v1 = 1.0*sqrt(v);
			arg1 = par[1]*u;
			}
	Double_t fitval = par[0] *v1*exp(-1.0*arg1);
	
	
	
		y = (mES[0] - par[3])/(par[4]);
		f = par[2] *TMath::Exp(-0.5*y*y);
	
		
			
	return fitval + f;	
	
	}

	

#include <iostream>
#include <fstream>
using namespace std;

void Fit()
{


char name[40],title[40];
char BestCutFiles[40];
sprintf(BestCutFiles,"dataFiles1.txt");

Int_t D_nEntries_MC,D_nEntries_Off,D_nEntries_On,D_nBins,D_Size,D_CutRes;
Int_t
D_nBins_mES_MC,D_nBins_deltaE_MC,D_nBins_mES_Off,D_nBins_deltaE_Off,D_nBins_mES_On,D_nBins_deltaE_On;
Double_t D_minfisher_MC, D_minfisher_On, D_minfisher_Off,D_maxfisher_MC,D_maxfisher_On,D_maxfisher_Off,P_NBB;
Double_t D_MIN,D_MAX;
/////////////////////////////////////////////
//Edges of graphs
//////////////////////////////////////////////

Double_t D_Low_edge_mES_MC,D_High_edge_mES_MC,D_Low_edge_deltaE_MC,D_High_edge_deltaE_MC;
Double_t D_Low_edge_mES_Off,D_High_edge_mES_Off,D_Low_edge_deltaE_Off,D_High_edge_deltaE_Off;
Double_t D_Low_edge_mES_On,D_High_edge_mES_On,D_Low_edge_deltaE_On,D_High_edge_deltaE_On;

//////////////////////////////////////////
//Box range
/////////////////////////////////////////

Double_t D_LEdge_mES_SR,D_HEdge_mES_SR,D_LEdge_deltaE_SR,D_HEdge_deltaE_SR;
Double_t D_LEdge_mES_GSB,D_HEdge_mES_GSB,D_LEdge_deltaE_GSB,D_HEdge_deltaE_GSB;

////////////////////////////////////////////
//Lengths of boxes
//////////////////////////////////////////
Double_t D_lengthSR_mES, D_lengthGSB_mES, D_lengthSR_deltaE,D_lengthGSB_deltaE;
Double_t D_areaSR,D_areaGSB;

//////////////////////////////////////////////////
//Best values
//////////////////////////////////////////////////
Double_t D_Ratio_Counting,D_Ratio_Integral;
Int_t D_imax;
Double_t D_BranchFrac,D_NBb1, D_BestFisher,D_SIGFracArrayBest,D_BGFracArrayBest,D_OnResSIGFracArrayBest;
Int_t D_OnResSRcntArrayBest,D_SRcntArrayBest,D_MCSRcntArrayBest,D_NBgBest,D_GSBcntArrayBest,D_MCGSBcntArrayBest;
Double_t D_Xi,D_Amp;
char trash[50];

ifstream read;
read.open(BestCutFiles,ios::in);

read >>trash>>D_nEntries_MC;
read >>trash>>D_nEntries_Off;
read >>trash>>D_nEntries_On;

read >>trash>>D_minfisher_MC;
read >>trash>>D_minfisher_Off;
read >>trash>>D_minfisher_On;
read >>trash>>D_maxfisher_MC;
read >>trash>>D_maxfisher_Off;
read >>trash>>D_maxfisher_On;
read >>trash>>D_MIN;
read >>trash>>D_MAX;

read >>trash>>D_CutRes;

read >>trash>>D_Size;

read >>trash>>D_nBins_mES_MC;
read >>trash>>D_nBins_deltaE_MC;
read >>trash>>D_Low_edge_mES_MC;
read >>trash>>D_High_edge_mES_MC;
read >>trash>>D_Low_edge_deltaE_MC;
read >>trash>>D_High_edge_deltaE_MC;

read >>trash>>D_nBins_mES_Off;
read >>trash>>D_nBins_deltaE_Off;
read >>trash>>D_Low_edge_mES_Off;
read >>trash>>D_High_edge_mES_Off;
read >>trash>>D_Low_edge_deltaE_Off;
read >>trash>>D_High_edge_deltaE_Off;

read >>trash>>D_nBins_mES_On;
read >>trash>>D_nBins_deltaE_On;
read >>trash>>D_Low_edge_mES_On;
read >>trash>>D_High_edge_mES_On;
read >>trash>>D_Low_edge_deltaE_On;
read >>trash>>D_High_edge_deltaE_On;

read >>trash>>D_LEdge_mES_SR;
read >>trash>>D_HEdge_mES_SR;
read >>trash>>D_LEdge_deltaE_SR;
read >>trash>>D_HEdge_deltaE_SR;


read >>trash>>D_LEdge_mES_GSB;
read >>trash>>D_HEdge_mES_GSB;
read >>trash>>D_LEdge_deltaE_GSB;
read >>trash>>D_HEdge_deltaE_GSB;

read >>trash>>D_lengthSR_mES;
read >>trash>>D_lengthGSB_mES;
read >>trash>>D_lengthSR_deltaE;
read >>trash>>D_lengthGSB_deltaE;

read >>trash>>D_areaSR;
read >>trash>>D_areaGSB;

read >>trash>>D_Ratio_Counting;
read >>trash>>D_Ratio_Integral;
read >>trash>>D_NBb1;
read >>trash>>D_BranchFrac;

read >>trash>>D_imax;
read >>trash>>D_BestFisher;
read >>trash>>D_Xi;
read >>trash>>D_Amp;
read >>trash>>D_SIGFracArrayBest;
read >>trash>>D_BGFracArrayBest;
read >>trash>>D_OnResSIGFracArrayBest;

read >>trash>>D_OnResSRcntArrayBest;
read >>trash>>D_SRcntArrayBest;
read >>trash>>D_MCSRcntArrayBest;
read >>trash>>D_NBgBest;
read >>trash>>D_GSBcntArrayBest;
read >>trash>>D_MCGSBcntArrayBest;

read.close();



	
	
	///////////////////////////////////////////////////////////////////////////////////
	//home_Run1_data
	//////////////////////////////////////////////////////////////////////////////////////
	
	//TFile *file  = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-SP-1218-Run1-Sig-noMVACut.root");
	//TFile *file1 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-Offpeak-data-bg-noMVACut.root");
	//TFile *file2 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-Run1-OnPeak-all-noMVACut.root");
	
	//home_Run2_data
	//////////////////////////////////////////////////////////////////////////////////////
	//TFile *file  = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-SP-1218-Run1-Sig-noMVACut.root");
	//TFile *file1 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-Run2-OffPeak-all-noMVACut.root");
	//TFile *file2 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-Run2-OnPeak-all-noMVACut.root");
	
	TFile *file  = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-SP-1218-RunAll-Sig-noMVACut.root");
	TFile *file1 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-RunAll-OffPeakNew-noMVACut.root");
	TFile *file2 = new TFile("/home/robert/uni/root/work/rootFiles/Kpp-RunAll-OnPeak-noMVACut.root");
	
	
	file  ->ls();
	file1 ->ls();
	file2 ->ls();
	
	TTree *tree   = (TTree*) file   -> Get("fitTree");
	TTree *tree1  = (TTree*) file1  -> Get("fitTree");
	TTree *tree2 =  (TTree*) file2  -> Get("fitTree");
	
Int_t nEvents_MC  = (Int_t) tree  -> GetEntries();
Int_t nEvents_Off = (Int_t) tree1 -> GetEntries();
Int_t nEvents_On = (Int_t) tree2 -> GetEntries();

Double_t fisher_MC  ,  mES_MC  , deltaE_MC;
Double_t fisher_Off ,  mES_Off , deltaE_Off;	
Double_t fisher_On  ,  mES_On  , deltaE_On;


tree   ->SetBranchAddress("fisher", &fisher_MC);
tree   ->SetBranchAddress("mES", &mES_MC);
tree   ->SetBranchAddress("deltaE", &deltaE_MC);

tree1  ->SetBranchAddress("fisher", &fisher_Off);
tree1  ->SetBranchAddress("mES", &mES_Off);
tree1  ->SetBranchAddress("deltaE", &deltaE_Off);

tree2  ->SetBranchAddress("fisher", &fisher_On);
tree2  ->SetBranchAddress("mES", &mES_On);
tree2  ->SetBranchAddress("deltaE", &deltaE_On);

cout<<"nEvents_MC:   "<<nEvents_MC<<endl;
cout<<"nEvents_ON:   "<<nEvents_On<<endl;
cout<<"nEvents_Off:   "<<nEvents_Off<<endl;
Int_t i,j,counter;
Int_t nBins = 50;
Int_t nBins_Kst = 50;//500;
Int_t nBins_Rho = 50;//500;
Int_t nbins_mES    = 50;//D_nBins_mES_MC;
Int_t nbins_deltaE = 50;//D_nBins_deltaE_MC;






Double_t minfisher_MC       = D_minfisher_MC;
Double_t maxfisher_MC       = D_maxfisher_MC;
Double_t minmES_MC          = tree  ->GetMinimum("mES");
Double_t maxmES_MC          = tree  ->GetMaximum("mES");
Double_t mindeltaE_MC       = tree  ->GetMinimum("deltaE");
Double_t maxdeltaE_MC       = tree  ->GetMaximum("deltaE");



Double_t minfisher_Off     = D_minfisher_Off;
Double_t maxfisher_Off     = D_maxfisher_Off;
Double_t minmES_Off        = tree1 ->GetMinimum("mES");
Double_t maxmES_Off        = tree1 ->GetMaximum("mES");
Double_t mindeltaE_Off     = tree1 ->GetMinimum("deltaE");
Double_t maxdeltaE_Off     = tree1 ->GetMaximum("deltaE");



Double_t minfisher_On     = D_minfisher_On;
Double_t maxfisher_On     = D_maxfisher_On;
Double_t minmES_On        = tree2 ->GetMinimum("mES");
Double_t maxmES_On        = tree2 ->GetMaximum("mES");
Double_t mindeltaE_On     = tree2 ->GetMinimum("deltaE");
Double_t maxdeltaE_On     = tree2 ->GetMaximum("deltaE");


//gSystem->Exit(0);
////////////////////////////////////////////////////////////////////////////////////
//Edges_of_SR&GSB_Regions
//////////////////////////////////////////////////////////////////////////////////


Double_t LmES_SR = D_LEdge_mES_SR;
Double_t HmES_SR = D_HEdge_mES_SR;
Double_t LmES_GSB = D_LEdge_mES_GSB;
Double_t HmES_GSB = D_HEdge_mES_GSB;

Double_t LdeltaE_SR = D_LEdge_deltaE_SR;
Double_t HdeltaE_SR = D_HEdge_deltaE_SR;
Double_t LdeltaE_GSB = D_LEdge_deltaE_GSB;
Double_t HdeltaE_GSB = D_HEdge_deltaE_GSB;


Double_t RBest = D_Ratio_Counting;


Int_t m;


Double_t BestCut = D_BestFisher;	
	

cout<<"bestfisher:           "<<BestCut<<endl;
TH2D *hhist_MC = new TH2D("hhist_MC" , "" , nbins_mES , 5.225 , maxmES_MC , nbins_deltaE , mindeltaE_MC , maxdeltaE_MC);
TH2D *hhist_Off = new TH2D("hhist_Off" , "" , nbins_mES , 5.225, maxmES_Off , nbins_deltaE , mindeltaE_Off , maxdeltaE_Off);
TH2D *hhist_On = new TH2D("hhist_On" , "" , nbins_mES , 5.225 , maxmES_On , nbins_deltaE , mindeltaE_On , maxdeltaE_On);
	

///////////////////////////////////////////////////////////////////////////

cout<<"Appling_mES_deltaE_cuts_for_SR_and_GSB_regions:           "<<endl;

////////////////////////////////////////////////////////////////////////////

for(j = 0; j<nEvents_MC; j++)
		{
		tree->GetEntry(j);
		
		if(fisher_MC < BestCut)
			{	
			
			
			if(mES_MC > 5.225)
				{
				if (deltaE_MC < HdeltaE_SR && deltaE_MC > LdeltaE_SR )		
					
					{
					
					hhist_MC->Fill(mES_MC , deltaE_MC);
							
					}
					
				
				}
			}
		}
//hhist_MC->Draw();


for(j = 0; j<nEvents_Off; j++)
		{
		tree1->GetEntry(j);
		
		if(fisher_Off < BestCut)
			{	
			
			if(mES_Off > 5.225)
				{
			
				if ((deltaE_Off < HdeltaE_SR) && (deltaE_Off > LdeltaE_SR) )		
					
					{
					
					hhist_Off->Fill(mES_Off , deltaE_Off);
							
					}
					
				}
		
			}
		}


for(j = 0; j<nEvents_On; j++)
		{
		tree2->GetEntry(j);
		
		if(fisher_On < BestCut)
			{
			
			if(mES_On > 5.225)
				{
				if (deltaE_On < HdeltaE_SR   && deltaE_On > LdeltaE_SR )		
					
					{
					
					hhist_On->Fill(mES_On , deltaE_On);
											
					}
					
		
			
				}
	
			}
	
	
	

	
		}


TH1D *PmES_On;
TH1D *PmES_Off;
TH1D *PmES_MC;

PmES_On =(TH1D*) hhist_On -> ProjectionX("a" , 0 , nbins_mES);
PmES_Off =(TH1D*) hhist_Off -> ProjectionX("a" , 0 , nbins_mES);
PmES_MC = (TH1D*) hhist_MC -> ProjectionX("a" , 0 , nbins_mES);

Double_t scale = 1.0*(1.0/(Double_t) 1000);
PmES_MC->Scale(scale);

Double_t pA0 = 21.4884;
Double_t pA1 = 43.7694;

Double_t p0 = 21.4884;
Double_t p1 = 43.7694;
Double_t p2 = 811;
Double_t p3 = 5.27826;
Double_t p4 = 0.00255814;

Double_t pg0 = 5.27826;
Double_t pg1 = 0.00255814;
Double_t pg2 = 800;


TCanvas *f1 = new TCanvas("f1","f1");
TF1 *gauss = new TF1("gausspeak",Gaussian,5.227,5.29,3);
gausspeak->SetParameters(pg0,pg1,pg2);
gausspeak->SetParNames("centre1","sigma1","peakmax1");
gausspeak->SetLineColor(kRed);

PmES_MC->Fit("gausspeak","r");

TCanvas *f2 = new TCanvas("f2","f2");
TF1 *Argus = new TF1("fit",Argus,5.227,5.29,2);
fit->SetParameters(pA0,pA1);
fit->SetParNames("Aamp","Axi");
fit->SetLineColor(kBlue);

PmES_Off->Fit("fit","r");

TCanvas *f3 = new TCanvas("f3","f3");
TF1 *Bothfit = new TF1("Fitname",Total,5.227,5.29,5);
Fitname->SetParameters(p0,p1,p2,p3,p4);
Fitname->SetParNames("amp","xi","peakmax","centre","sigma");
Fitname->SetLineColor(kRed);

PmES_On->Fit("Fitname","r");





}