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

#include "TFile.h"
#include "TTree.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TStopwatch.h"
#include "TLegend.h"

#include "TpcPoint.h"	//include Paths in rootlogon.C!
#include "CbmMCTrack.h"
#include "g_defines.h"
#include "CChainGetter.h"


//Reads the Filenames listet in filename and stores them in a vector
//@param FileExtension : an optional string to be added at the end of verry name in vFileNames
//void GetFileNames(string filename, vector<string> &vFileNames, string fileExtension="");

//Get the Root-Tree. If a Root-Tree-File is specified, it will be used. Otherwise the gcstrFileListFile is evaluated    
//TChain * GetTree(TString SingleFilename=""); 

//***************
//Uses the PDGMap to name the histograms, histograms are accesible over their name
void CreateHistograms();

void FillHistograms(string name, Double_t momentum, Double_t dedx);

void FillHistograms(string name, Double_t momentum, vector<Double_t> vDedx);

void  saveHistograms(string strSaveHistoFile);
//*****************

//________________________________________
//the file where the RootTreeFile-List is stored 
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/AutoGenTrees";
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/AltTrees";
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/ShortTrees";
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/AutoGenTreesTest";
//const string gcstrFileListFile="/afs/e18/data/panda/MC/FAIRRoot/RUNMC/rs1/rs1_Gesamt.txt";
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/AutoGenTrees2_short";
//const string gcstrFileListFile="/afs/e18/user/vmichael/pandaroot/pandaroot/AutoSinglePGenTrees";


//this will create a root tree where the dedx histograms for single Particles (one can fit them, make dedx distributions, etc... ) 
//const string g_cstrSaveDeDxParticleHistos="/afs/e18/data/panda/MC/FAIRRoot/RUNMC/rs1/dedxHistosTest.root";

struct TrackInfo
{
	TrackInfo() : nPdgCode(0), nPoint(0), sumDedx(0), sumDe(0), sumDx(0), momentum(0) {}
	~TrackInfo() {}
	
	void Fill(Double_t de, Double_t dx);
	Double_t GetSumDe_SumDx();
	Double_t GetSumDeDx_nPoint();
	
	Int_t nPdgCode;
	Int_t nPoint;
	Double_t sumDedx;
	Double_t sumDe;
	Double_t sumDx;
	Double_t momentum;
	vector<Double_t> vDeDx;
};

Double_t TrackInfo::GetSumDe_SumDx()
{
	if(sumDx)	{
		return sumDe/sumDx;	//[kev/cm]
	}
	else	{
		return -1;
	}	
}

Double_t TrackInfo::GetSumDeDx_nPoint()
{
	if(nPoint)	{
		return sumDedx/nPoint;	//[kev/cm]
	}
	else	{
		return -1;
	}
}

void TrackInfo::Fill(Double_t de, Double_t dx)
{
	de=de*1000*1000;	//keV
	sumDe+=de;
	sumDx+=dx;
	if(dx)	{
		sumDedx+=(de/dx);
		vDeDx.push_back(de/dx);
	}
	nPoint++;
}

void dedx(string strFileListFile, string strSaveHistoFile, TString mwa="" )
{
	TStopwatch timer;
	timer.Start();
	cout << "dedx - START" << endl;
	//________________________________________
	//---------insert Particles in the PDGMap
	

	//PDGMap[111]="Pi_0";	//TODO: kein Fehler bei nicht existierenden Teilchen
	//PDGMap[-13]="Muon+";	
	PDGMap[13]="Muon-";
	PDGMap[-11]="e+";
	//PDGMap[11]="e-";
	//PDGMap[14]="Neutrino";
	//PDGMap[22]="Gamma";
	PDGMap[211]="Pi+";
	PDGMap[2212]="p+";
	PDGMap[321]="K+";
	PDGMap[0]="Summary";	

	CChainGetter cbmsim;
	cbmsim.SetFileListFile(strFileListFile);
	cbmsim.SetFileExtension(".mc.root");
	TChain* pTree=cbmsim.GetChain(mwa);
	CreateHistograms();
	
	
	//Double_t P_set=0.4;
	vector<Double_t> vP_set;
	Double_t tolerance=0.01;
	for(Double_t p=0.00+tolerance; p<3.0; p+=(tolerance))
	{
		vP_set.push_back(p);
	}
	
	//CreateDeDxDistribution(P_set);
	//CreateDeDxDistribution(vP_set);
	
	//________________________________________
	//----------- Read out the TpcPoint Branch of the pTree and fill the data in the histograms
	TClonesArray *pPointArray=NULL;
	pTree->SetBranchAddress("TpcPoint", &pPointArray);
	
	TClonesArray *pTrackArray=NULL;
	pTree->SetBranchAddress("MCTrack", &pTrackArray);
		
	Long64_t nTreeEntries=pTree->GetEntries();
	printf("nTreeEntries: %lld\n", nTreeEntries);
	
	 timer.Stop();
 	 Double_t rtime = timer.RealTime();
 	 Double_t ctime = timer.CpuTime();
	 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
	 timer.Reset();
	 timer.Start();	 
	 
	//in jedem TreeEntry sind die Tracks und die Punkte von einem Event
	//sammle die Punkte, die zu einem Track gehören und bilde Mittelwert de/dx
	//the momentum for a track is known, so i will use a mean value for de/dx for a track
	//eigentlich müsste ich den Energieverlust ja zum Impuls dazurechnen, aber der ist klein
	//Sammle diese Mittelwerte für jedes Teilchen und tu sie in das entsprechende Histogram
	for(int i=0; i<nTreeEntries; i++)
	{
		pTree->GetEntry(i);
		if(i%10000==0)	{
			printf("Reading Entry %d\n", i);
		}
		Int_t nPointArrayEntries=pPointArray->GetEntriesFast();
		//printf("nPointArrayEntries: %d\n", nPointArrayEntries);//for every Track/foundParticle get the mean de/dx and fill result in the histogram
		map<Int_t, TrackInfo> TrackInfoMap;	//für jeden Track alle dedx speichern um damit am Ende den Mittelwert zu bilden
		for(Int_t j=0; j < nPointArrayEntries; j++)				//for every point
		{
			TpcPoint *pPoints=(TpcPoint*)pPointArray->At(j);
			assert(pPoints);
			Int_t nID=pPoints->GetTrackID();				
			//printf(" ---ID: %d--- ", nID);
			CbmMCTrack *pFoundTrack=(CbmMCTrack*)pTrackArray->At(nID);
			assert(pFoundTrack);
			Int_t nPdgCode = pFoundTrack->GetPdgCode();			//see what particle it is

			if(nPdgCode)	{	
				Double_t momentum = ((TVector3)(pFoundTrack->GetMomentum())).Mag();
				Double_t de = pPoints->GetEnergyLoss();
				Double_t dx = pPoints->GetLength();
				
				//according to its TrackID 
				(TrackInfoMap[nID]).Fill(de, dx);
				(TrackInfoMap[nID]).nPdgCode=nPdgCode;
				(TrackInfoMap[nID]).momentum=momentum;
			}
		}
		map<Int_t, TrackInfo >::iterator itTrack;
		for(itTrack=TrackInfoMap.begin(); itTrack!=TrackInfoMap.end(); itTrack++)
		{
			Double_t momentum=(itTrack->second).momentum;
			Int_t nPdgCode=(itTrack->second).nPdgCode;
			Double_t DeDx=(itTrack->second).GetSumDe_SumDx();
			
			if(nPdgCode && PDGMap.find(nPdgCode)!=PDGMap.end())	{
				//________________________________________
				//Fill dE/dx vs. momentum
				FillHistograms(PDGMap[nPdgCode], momentum, DeDx);

				//________________________________________
				//Also fill in the dedx-Distribution Histograms
				//TODO: At this point i should fill  with arbitrary momentum, but there has to be made a decision in which
				//Histogram the data should be filled
				//FillDeDxDistribution(PDGMap[nPdgCode], DeDx, momentum, vP_set, tolerance );
			}
			else if(nPdgCode && PDGMap.find(nPdgCode)==PDGMap.end() )	{
				//printf("Found unmapped Particle: %d!\n", nPdgCode);
			}
		}
	}
	timer.Stop();
 	rtime = timer.RealTime();
 	ctime = timer.CpuTime();
	printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
	timer.Reset();
	timer.Start();
	
	//________________________________________
	//Create Fits and draw Histograms
	{
		//vector<TF1*> vFits;
		//map<string, TF1*> FitMap;
		//CreateFits(vFits, FitMap);
		//traceDrawing(FitMap);
		//DrawSummary(vFits);
	}
	saveHistograms(strSaveHistoFile);
	//SepPoHandler TheSep;
	//________________________________________
	//Draw the dedx distribution histograms and create gaus fits for them
	{
		//map<Double_t, map<string, TF1*> > FitMap;
		//CreateDeDxDistributionGausFits(vP_set, FitMap);
		//DrawDeDxDistribution(vP_set, FitMap);
		//TheSep.FillSeparationPower(FitMap);	//TODO: delete unused objects
		//TheSep.Draw();
	}
	
	timer.Stop();
 	rtime = timer.RealTime();
 	ctime = timer.CpuTime();
	printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
	timer.Reset();
	
}


void CreateHistograms()
{	
	//________________________________________
	//----------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++)
	{
		{
			ostringstream HistogramName;
			HistogramName << cit->second;
			ostringstream HistoTitle;
			HistoTitle << "de/dx";
			//the histograms are accesible over their names, i hope
			TH2F *pHist = new TH2F(HistogramName.str().c_str(), HistoTitle.str().c_str() , 300, 0, 3, 1300, -2, 11);
			pHist->GetXaxis()->SetTitle("|P|");
			pHist->GetYaxis()->SetTitle("dE/dx");
		}
	}
}

void FillHistograms(string name, Double_t momentum, Double_t dedx)
{
	//________________________________________
	//Fill  the Histograms with data
	{
		ostringstream HistoName;
		HistoName << name.c_str();
		TH2F *pCurHisto=(TH2F*)gDirectory->Get(HistoName.str().c_str());
		if(pCurHisto)	{
			pCurHisto->Fill(momentum, dedx );
		}
	}
	
	//________________________________________
	//Also fill in the summary histograms
	{
		ostringstream HistoName;
		HistoName << "Summary";
		TH2F *pSummaryHisto=(TH2F*)gDirectory->Get(HistoName.str().c_str());
		if(pSummaryHisto)	{
			pSummaryHisto->Fill(momentum, dedx );
		}
	}
}

void FillHistograms(string name, Double_t momentum, vector<Double_t> vDedx)
{
	vector<Double_t>::const_iterator cit;
	for(cit=vDedx.begin(); cit!=vDedx.end(); cit++)
	{
		FillHistograms(name, momentum, *cit);
	}
}

void saveHistograms(string strSaveHistoFile)
{
	TObjArray Histos(0);
	map<Int_t,string>::const_iterator cit;
	for(cit=PDGMap.begin(); cit!=PDGMap.end(); cit++)
	{
		ostringstream HistoName;
		HistoName << cit->second;
		printf("Saving HistoName: %s\n", HistoName.str().c_str());
		TH2F *pCurHisto=(TH2F*)gDirectory->Get(HistoName.str().c_str());
		assert(pCurHisto);
		Histos.Add(pCurHisto);
	}
	TFile f(strSaveHistoFile.c_str(),"recreate");
	Histos.Write();
	f.Close();
}