#include "TFile.h"
#include "TH2F.h"
#include "SepPoHandler.h"
#include "SeparationPower.h"
#include <iostream>
#include <set>
#include "g_defines.h"
#include "DEDXFitCreator.h"
#include "DEDXDistribana.h"

//Gets a root file where histograms of the dedx-Distribution is stored
void seppo(string Filename, bool bShowDEDXDistribution=false)
{
	//________________________________________
	//---------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";
	TFile *f = new TFile(Filename.c_str());
	//TFile *f = new TFile("/afs/e18/data/panda/MC/FAIRRoot/RUNMC/rs1/dedxHistosTest.root");
	
	assert(f);
	
	TH2F *h1 = (TH2F *) gROOT->FindObject("p+");
	h1->Draw();
	//TH2F *h2 = (TH2F *) gROOT->FindObject("Summary");
	//h2->Draw();	
	
	//h->SetShowProjectionY(1);
	//f->Close();
	vector <Double_t> vP_set;
	set<Double_t> sP_set;
	//Finde alle Objekte (d.h. die Histogramme ausser die Summary mit der PDGMap
	map<Int_t,string>::const_iterator cit;
	for(cit=PDGMap.begin(); cit!=PDGMap.end(); cit++)	
	{
		TH2F *h = (TH2F *) gROOT->FindObject((cit->second).c_str());
		if(h)	{
			Double_t offset = h->GetBinLowEdge(1);	//to calculate P
			Double_t offset_next = h->GetBinLowEdge(2);
			Double_t max_P = h->GetBinLowEdge(h->GetNbinsX()) + offset_next;
			Double_t binWidth = offset_next - offset;
			cout << "min_P: " << offset << " binWidth: " << binWidth << "max_P: " << max_P << endl;
			Int_t nbin = h->GetNbinsX();	//number of bins
			cout << "nbin: " << nbin << endl;
			
			
			//if(nbin > 100) {nbin=100;}
			Int_t binStep=1;		//immer drei bins zusammenfassen
			assert(binStep%2==1);	//nur ungerade zulassen
			for(Int_t i=0; i<=nbin; i+=binStep)	
			{	
				//das Histogram zu einem Teilchen wird für viele Impulswerte projeziert und 
				Int_t first_bin=i;
				Int_t last_bin=i+binStep;
				cout << "ProjectionY from bin number " << first_bin << " to " << last_bin << endl; 
				TH1F * pHisto=(TH1F*)h->ProjectionY("temp", first_bin,last_bin );  //z.b. bin 1 , 2, 3        //pHisto will get another name in DeDxFromExistingHistogram
				//P has to be calculated
				if(pHisto->GetEntries() < 20)	{
					cout << "Skipped this projection (less than 20 Entries)" << endl;
					continue;
				}
				Int_t dist=(binStep/2) - 0.5;				//distance from first_bin to centerbin TODO: Überprüfen ob int
				Double_t P=(first_bin)*binWidth+(int)offset+dist;	
				
				if (sP_set.insert(P).second )	{
					vP_set.push_back(P);
					CreateDeDxDistributionSummary("Summary", P);
				}
				DeDxFromExistingHistogram(pHisto, cit->second, P);
				
			}
		}
	}
	

	//________________________________________
	//Create Fits and draw Histograms
	{
		vector<TF1*> vFits;
		map<string, TF1*> FitMap;
		CreateFits(vFits, FitMap);
		traceDrawing(FitMap);
		DrawSummary(vFits);
	}
	
	
	//cout << "start fitting" << endl;
	vector<Double_t>::const_iterator v;
	for(v=vP_set.begin(); v!=vP_set.end(); v++)
	{
		cout << "P: " << *v << " - ";
	}
	cout << endl;
	char c;
	cin >> c;
	

	SepPoHandler TheSep;
	//________________________________________
	//Draw the dedx distribution histograms and create gaus fits for them
	{
		map<Double_t, map<string, TF1*> > FitMap;
		CreateDeDxDistributionGausFits(vP_set, FitMap);
		if(bShowDEDXDistribution)	{
			DrawDeDxDistribution(vP_set, FitMap);
		}
		TheSep.FillSeparationPower(FitMap);	//TODO: delete unused objects
		TheSep.Draw();
	}

}