#include <iostream>
#include "TFile.h"
#include "TROOT.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "PndSdsMCPoint.h"
#include "PndMCTrack.h"
#include "TH1D.h"
#include "TApplication.h"
#include "TF1.h"
#include "TStyle.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TFile.h"
using namespace std;
//void ReadTree();

/* void ReadTree(){
 double renergyloss;
 int rparticletype;
 bool risprimary;
 double rposx;
 double rposy;
 double rposz;
 double rposr;
 double renergy;
 double rvertexx;
 double rvertexy;
 double rvertexz;

 DataTree->SetBranchAddress("energyloss",&renergyloss);

 DataTree->SetBranchAddress("particletype",&rparticletype);

 DataTree->SetBranchAddress("isprimary",&risprimary);

 DataTree->SetBranchAddress("posx",&rposx);
 DataTree->SetBranchAddress("posy",&rposy);
 DataTree->SetBranchAddress("posz",&rposz);
 DataTree->SetBranchAddress("posr",&rposr);

 DataTree->SetBranchAddress("vertexx",&rvertexx);
 DataTree->SetBranchAddress("vertexy",&rvertexy);
 DataTree->SetBranchAddress("vertexz",&rvertexz);
 } */

int main() {
	TApplication* HistoApp = new TApplication("HistoApp", 0, 0);
	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
	TFile* MC = new TFile("Lumi_MC.root", "read");
	TFile Dataout("datenanalyse.root", "recreate");
	gStyle->SetOptFit(111);

	if (MC->IsZombie()) {
		cout << "File not here!" << endl;
		return 0;
	}

	TTree* DataTree;

	double Lostevents = 0;
	double vertexx;
	double vertexy;
	double vertexz;
	double energyloss;
	double MaxEnergyloss;
	double SumEnergyloss = 0;
	double energy;
	int particletype;
	bool isprimary;
	double posx;
	double posy;
	double posz;
	double posr;
	int detectorid;

	TTree* MCDat = (TTree*) MC->Get("cbmsim");
	TH1D* HEnergyloss = new TH1D("Energyloss", "Energyloss per Point", 10000, 0.5,
			MaxEnergyloss);
	Dataout.cd();
	DataTree = new TTree("DataTree", "Particledatatree");
	DataTree->Branch("energyloss", &energyloss);
	DataTree->Branch("particletype", &particletype);
	DataTree->Branch("isprimary", &isprimary);
	DataTree->Branch("posx", &posx);
	DataTree->Branch("posy", &posy);
	DataTree->Branch("posz", &posz);
	DataTree->Branch("posr", &posr);
	DataTree->Branch("vertexx", &vertexx);
	DataTree->Branch("vertexy", &vertexy);
	DataTree->Branch("vertexz", &vertexz);
	DataTree->Branch("energy", &energy);
	DataTree->Branch("detectorid", &detectorid);

	if (MCDat) {
		TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
		MCDat->SetBranchAddress("MCTrack", &true_tracks); //True Track to compare

		TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
		MCDat->SetBranchAddress("LMDPoint", &true_points); //True Points to compare

		cout << "Tree found!" << MCDat->GetEntries() << endl;

		int nEntries = MCDat->GetEntries();
		for (int iEvent = 0; iEvent < nEntries; iEvent++) {
			MCDat->GetEntry(iEvent);
			cout << iEvent << " Found " << true_tracks->GetEntries()
					<< " Tracks" << endl;
			cout << iEvent << " Found " << true_points->GetEntries()
					<< " Points" << endl;

			double Eventeloss = 0;
			int MaxPoints = true_points->GetEntries();
			for (int iPoint = 0; iPoint < MaxPoints; iPoint++) {

				PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(
						iPoint);

				int TrackId = mcpoint->GetTrackID();
				if (TrackId < 0) {
					cout << "warning no Track ID!" << endl;
					continue;
				}
				PndMCTrack* mctrack = (PndMCTrack*) true_tracks->At(TrackId);
				if (!mctrack->IsGeneratorCreated()) {
					Lostevents = Lostevents + 1;
				//	continue;
				}

				energyloss = mcpoint->GetEnergyLoss();
				if (MaxEnergyloss < energyloss)
					MaxEnergyloss = energyloss;
				Eventeloss = Eventeloss + energyloss;
				SumEnergyloss = SumEnergyloss + energyloss;
				cout << SumEnergyloss << " ";
				cout << Lostevents << " " << "Lost events ";
				//	cout<<energyloss;
				particletype = mctrack->GetPdgCode();

				isprimary = mctrack->IsGeneratorCreated();

				posx = mcpoint->GetX();
				posy = mcpoint->GetY();
				posz = mcpoint->GetZ();
				posr = sqrt((posx) * (posx) + (posy) * (posy));

				TLorentzVector p4mom(mctrack->Get4Momentum());
				energy = p4mom.Energy();

				TVector3 vertex(mctrack->GetStartVertex());
				vertexx = vertex.X();
				vertexy = vertex.Y();
				vertexz = vertex.Z();

				detectorid = mcpoint->GetSensorID();

				DataTree->Fill();

			//	HEnergyloss->Fill(mcpoint->GetEnergyLoss() * 1e6);
				HEnergyloss->Fill(Eventeloss);

// cout<<"Energieverlust"<<mcpoint->GetEnergyLoss()<<endl;
			}
		}
		HEnergyloss->Draw();
		HEnergyloss->GetXaxis()->SetTitle("Energyloss [keV]");
		HEnergyloss->GetYaxis()->SetTitle("Events");
		//	TF1* Fitlandau = new TF1("Fitlandau","landau",0.5,40);

		//	TF1* Fit1 = new TF1("Fit1","([1]/(sqrt(2*TMath::Pi())*[0])*exp(-0.5*((x-[2])**2)/([0]**2)))",0.5,40);
		//	Fit1->SetParNames("#sigma","Events","#mu");
		//	Fit1->SetParameters(3,540,8);
		//	HEnergyloss->Fit(Fitlandau);
	}

	Dataout.Write();
	HistoApp->Run();
	return 0;
}