#include <sstream>
#include <fstream>
#include <map>
#include <iostream>

#include "TCanvas.h"
#include "TF1.h"
#include "TGraph.h"
#include "TLegend.h"

#include "g_defines.h"
#include "DEDXDistribana.h"

Int_t g_cnMinEntries=15;	// this is the minimum number of entrys a histogram has to have to be fitted with a gaus function

void CreateDeDxDistribution(vector<Double_t> vP_set)
{
	vector<Double_t>::const_iterator cit;
	for(cit=vP_set.begin(); cit!=vP_set.end(); cit++)
	{
		CreateDeDxDistribution(*cit);
	}
}

void CreateDeDxDistribution(const Double_t P_set)
{
	//________________________________________
	//----------Create a Histogram for every Particle and a Summary histograms if wanted
	map<Int_t,string>::const_iterator cit;
	for(cit=PDGMap.begin(); cit!=PDGMap.end(); cit++)
	{
		if(cit->first)	{
			ostringstream HistogramName;
			HistogramName << cit->second << "_de_dx_" << P_set;
			ostringstream HistoTitle;
			HistoTitle << "de/dx for set |P|: " << P_set << "GeV/c";
			
			//the histograms are accesible over their names, i hope
			TH1F *pHist = new TH1F(HistogramName.str().c_str(), HistoTitle.str().c_str(), 700, 0,13 );
		//	dedxlog << "Creating Histogram: " << HistogramName.str().c_str() << endl;
			assert(pHist);
			pHist->GetXaxis()->SetTitle("dE/dx");
		}
		else	{
			CreateDeDxDistributionSummary(cit->second, P_set);
		}
	}
}


void DeDxFromExistingHistogram(TH1F *pHisto, string strName, Double_t P)
{
	assert(pHisto);
	{
		ostringstream HistoName;
		HistoName << strName.c_str() << "_de_dx_" << P;
		pHisto->SetNameTitle(HistoName.str().c_str(), HistoName.str().c_str());
		pHisto->GetXaxis()->SetTitle("dE/dx");
		cout << "Got existing Histogram " << HistoName.str().c_str() << endl;
	}
}


void CreateDeDxDistributionSummary(string strName, Double_t P)
{
	if(strName!="Summary")	{ return; }
	ostringstream HistogramName;
	HistogramName << "_de_dx_distSummary_" << P;
	ostringstream HistoTitle;
	HistoTitle << "de/dx Distribution Summary for |P| = " << P << " GeV/c";
	
	//the histograms are accesible over their names, i hope
	TH1F *pHist = new TH1F(HistogramName.str().c_str(), HistoTitle.str().c_str(), 700, 0,13 );
	assert(pHist);
	pHist->GetXaxis()->SetTitle("dE/dx");
	cout << "Created DeDx Distribution Summary" << endl;
}

void FillDeDxDistribution(string strName, Double_t input, Double_t P, const vector<Double_t> vP_set, Double_t tolerance)
{
	vector<Double_t>::const_iterator cit;
	bool bFilledIn=0;
	for(cit=vP_set.begin(); cit!=vP_set.end(); cit++)
	{
		bFilledIn=FillDeDxDistribution(strName, input, P,*cit , tolerance);
		if(bFilledIn)	{
			break;
		}
	}
	if(!bFilledIn)	{
	//	dedxlog << "Momentum " << P << " is not filled in!" << endl;
	}
}

bool FillDeDxDistribution(string strName, Double_t input, Double_t P,const Double_t P_set, Double_t tolerance)
{
	if( !(P >= P_set-tolerance && P <= P_set+tolerance) )	{ 
		return 0;	
	} 
	
	//Fill the distribution for one Particle (strName)
	{
		ostringstream HistoName;
		HistoName << strName.c_str() << "_de_dx_" << P_set;
		TH1F *pHisto=(TH1F*)gDirectory->Get(HistoName.str().c_str());
		if(pHisto)	{
			pHisto->Fill( input );
			return 1;
		}
	}
	
	//Fill the summary
	{
		ostringstream HistoName;
		HistoName << "_de_dx_distSummary_" << P_set;
		TH1F *pHisto=(TH1F*)gDirectory->Get(HistoName.str().c_str());
		if(pHisto)	{
			pHisto->Fill( input );
			return 1;
		}
	}
	
	return 0;
}

void CreateDeDxDistributionGausFits(const vector<Double_t> vP_set, map<Double_t, map<string, TF1*> > &FitMap)
{
	vector<Double_t>::const_iterator cit;
	for(cit=vP_set.begin(); cit!=vP_set.end(); cit++)
	{
		map<string, TF1*> FitMap_setP;
		CreateDeDxDistributionGausFits(*cit, FitMap_setP);
		FitMap[*cit]=FitMap_setP;
	}
}

void CreateDeDxDistributionGausFits(const Double_t P_set, map<string, TF1*> &FitMap)
{
	short int Colors[]={kRed, kGreen, kBlue, kMagenta, kYellow, kCyan, kRed, kGreen, kBlue, kMagenta, kYellow, kCyan};
	assert( sizeof(Colors)/sizeof(short int) >= PDGMap.size() );
	
	Int_t nColor=0;
	map<Int_t,string>::const_iterator cit;
	for(cit=PDGMap.begin(); cit!=PDGMap.end(); cit++)
	{
		if(cit->first)	{
			ostringstream HistoName;
			HistoName << cit->second << "_de_dx_" << P_set;
			TH1F *pCurHisto=(TH1F*)gDirectory->Get(HistoName.str().c_str());
			if(!pCurHisto)	{
				cout << "CreateDeDxDistributionGausFits - Skipping: " << HistoName.str().c_str() << endl;
				continue;
			}
			if( pCurHisto->GetEntries() >= g_cnMinEntries )	{
				cout << "CreateDeDxDistributionGausFits - Creating Gaus Fits for " << HistoName.str().c_str() << endl;
				Double_t upperFitLimit=pCurHisto->GetBinLowEdge(pCurHisto->GetNbinsX());
				cout << "CreateDeDxDistributionGausFits - Using upper fit limit: " << upperFitLimit << endl;
				TF1 *fgaus = new TF1("fgaus","gaus",0 , upperFitLimit+100);
				assert(fgaus);
				pCurHisto->Fit("fgaus","RNO");
				fgaus->SetLineColor(Colors[nColor++]);
				ostringstream FunctionName;
				FunctionName << cit->second << "_" << P_set;
				FitMap[FunctionName.str()]=fgaus;
			}
			else	{
				cout << "CreateDeDxDistributionGausFits - Skipping: " << HistoName.str().c_str() << endl;
				cout << "Reason: Less than " << g_cnMinEntries << " Entries" << endl;
			}
		}
	}
}

void DrawDeDxDistribution(const vector <Double_t> vP_set,map<Double_t, map<string, TF1*> > &FitMap)
{
	vector<Double_t>::const_iterator cit;
	for(cit=vP_set.begin(); cit!=vP_set.end(); cit++)
	{ 
		DrawDeDxDistribution(*cit, FitMap[*cit]);
	}
}

void DrawDeDxDistribution(const Double_t P_set, map<string, TF1*> &FitMap)
{
	//cout << "DrawDeDxDistribution -Start" << endl;
	//________________________________________
	//Draw DeDx-Distributions for the momentum P_set in new Canvas
	ostringstream CanvasName;
	CanvasName << "dedxDistribution_Canvas" << P_set;
	ostringstream CanvasTitle;
	CanvasTitle << "dedx for set P value " << P_set << " GeV/c";
	TCanvas *pCanvas = new TCanvas(CanvasName.str().c_str(), CanvasTitle.str().c_str() , 30,30,500,350);
	assert(pCanvas);
	
	if( (PDGMap.size())%2 == 0  )	{	
		pCanvas->Divide(2,(PDGMap.size())/2);
	}
	else	{
		pCanvas->Divide(2,(PDGMap.size()+1)/2);
	}	
	
	Int_t nCanvas=1;
	map<Int_t,string>::const_iterator cit;
	for(cit=PDGMap.begin(); cit!=PDGMap.end(); cit++)	//draw in a subpad for every particle + summary
	{
		if(cit->first)	{
			pCanvas->cd(nCanvas++);
			ostringstream HistoName;
			HistoName << cit->second << "_de_dx_" << P_set;
			TH1F *pCurHisto=(TH1F*)gDirectory->Get(HistoName.str().c_str());
			if(pCurHisto) {
				pCurHisto->Draw();
				ostringstream FunctionName;
				FunctionName << (cit->second).c_str() << "_" << P_set;
				if(FitMap.find(FunctionName.str())!=FitMap.end())	{
					FitMap[FunctionName.str()]->Draw("SAME");
				}
			}
		}
		else	{	//the summary
			pCanvas->cd(nCanvas++);
			ostringstream HistoName;
			HistoName << "_de_dx_distSummary_" << P_set;
			TH1F *pCurHisto=(TH1F*)gDirectory->Get(HistoName.str().c_str());
			if(pCurHisto)	{
				pCurHisto->Draw();
				
				TLegend *legend=new TLegend(0.39,0.6,1.1,0.8);
				legend->SetTextFont(72);
				legend->SetTextSize(0.04);
	
				map<string, TF1*>::const_iterator cit_FM;
				for(cit_FM=FitMap.begin(); cit_FM!=FitMap.end(); cit_FM++)
				{				
					ostringstream FunctionName;
					FunctionName << (cit->second).c_str() << "_" << P_set;
					(cit_FM->second)->Draw("SAME");
					ostringstream legendText;
					legendText << (cit_FM->first).c_str() << " Mean: " << (cit_FM->second)->GetParameter("Mean") << ", Sigma: " << (cit_FM->second)->GetParameter("Sigma");
					legend->AddEntry( cit_FM->second, legendText.str().c_str() );
					
				}
				legend->Draw();
			}	
		}
	}
	//cout << "DrawDeDxDistribution -Ende" << endl;
}