/*
 * PndAnalysisForwardTask.cpp
 *
 *  Created on: Jul 15, 2016
 *      Author: kibellus
 */

#include <PndAnalysisForwardTask.h>

PndAnalysisForwardTask::~PndAnalysisForwardTask() {
	// TODO Auto-generated destructor stub
}

InitStatus PndAnalysisForwardTask::Init(){
	fIoman = FairRootManager::Instance();
	fHits = (TClonesArray*) fIoman->GetObject("FTSHit");
	fSolution = (TClonesArray*) fIoman->GetObject("FinalSolution");
	fMCTracks = (TClonesArray*) fIoman->GetObject("MCTrack");
	fIdealTrack = (TClonesArray*) fIoman->GetObject("IdealTrack");
	return kSUCCESS;
}

void PndAnalysisForwardTask::Exec(Option_t*){
	cout << "####################  Event nr " << eventNum << "  ####################" << endl;

	//get the data
	map<Int_t,vector<PndFtsHit*>> mcToHits = getMcTracks();
	map<Int_t,vector<PndFtsHit*>*> hitsToTrack = getReconstructedTracks();

	//analyze the tracks
	cout << "Found " <<mcToHits.size() << " MC-Tracks"<< endl;
	for(map<Int_t,vector<PndFtsHit*>>::iterator i=mcToHits.begin();i!=mcToHits.end();i++){
		Int_t mcTrackID = i->first;
		vector<PndFtsHit*> hits = i->second;
		cout << "-----------  Analyze MC-Track " << i->first <<"  -----------" << endl;
		analyzeMCTracks(hits,hitsToTrack,mcTrackID, mcToHits.size());
	}

	anaSpecialCases(mcToHits);
	eventNum++;
}

void PndAnalysisForwardTask::FinishEvent(){
	fHits->Delete();
}

void PndAnalysisForwardTask::Finish(){
	cout << "--------------------- Quality Results: ----------------------" << endl;
	Int_t trackCount = completeClean+completeUnclean+incompleteClean+incompleteUnclean+notFound+notEnoughFound;
	Double_t completeCleanRel = ((Double_t) completeClean)/trackCount*100;
	Double_t completeUncleanRel = ((Double_t) completeUnclean)/trackCount*100;
	Double_t incompleteCleanRel = ((Double_t) incompleteClean)/trackCount*100;
	Double_t incompleteUncleanRel = ((Double_t) incompleteUnclean)/trackCount*100;
	Double_t notFoundRel = ((Double_t) notFound)/trackCount*100;
	Double_t notEnoughFoundRel = ((Double_t)notEnoughFound)/trackCount*100;
	Double_t foundTracksRel = 100-notFoundRel-notEnoughFoundRel;
	cout << "Tracks: " << trackCount << endl;
	cout << "Found " << trackCount-notFound-notEnoughFound << " Tracks (" << foundTracksRel << "%)" << endl;
	cout << "	Found complete and clean     : " << completeClean << " (" << completeCleanRel << "%)" << endl;
	cout << "	Found complete and unclean   : " << completeUnclean << " (" << completeUncleanRel << "%)" << endl;
	cout << "	Found incomplete and clean   : " << incompleteClean << " (" << incompleteCleanRel << "%)" << endl;
	cout << "	Found incomplete and unclean : " << incompleteUnclean << " (" << incompleteUncleanRel << "%)" << endl;
	cout << "Not found(under 70%): " << notEnoughFound << " (" << notEnoughFoundRel << "%)" << endl;
	cout << "Not found : " << notFound << " (" << notFoundRel << "%)" << endl;

	cout << "correlation:" << endl;
	cout << "momentum:" << endl;
	for(int i=0;i<10;i++){
		cout << momentumArray[i] << "  ";
	}
	cout << endl;
	cout << "hits:" << endl;
	for(int i=0;i<10;i++){
		cout << hitArray[i] << "  ";
	}
	cout << endl;
	cout << "TrackCount:" << endl;
	for(int i=0;i<10;i++){
		cout << trackCountArray[i] << "  ";
	}
	cout << endl;


	cout << "sum of all tracks:" << endl;
	cout << "momentum:" << endl;
	for(int i=0;i<10;i++){
		cout << momentumArray2[i] << "  ";
	}
	cout << endl;
	cout << "hits:" << endl;
	for(int i=0;i<10;i++){
		cout << hitArray2[i] << "  ";
	}
	cout << endl;
	cout << "TrackCount:" << endl;
	for(int i=0;i<10;i++){
		cout << trackCountArray2[i] << "  ";
	}
	cout << endl;

	//write to file
	ofstream ofs;
	ofs.open("results.csv",ios::app);

	ofs << trackCount << ";";
	ofs << completeClean<< ";";
	ofs << completeUnclean<< ";";
	ofs << incompleteClean<< ";";
	ofs << incompleteUnclean<< ";";
	ofs << notEnoughFound<< ";";
	ofs << notFound << ";" << endl;

	Double_t ohRel =  100*((Double_t)oneHitPerLayer)/trackCount;
	Double_t mhRel =  100*((Double_t)moreThan3Hits)/trackCount;
	Double_t mtRel =  100*((Double_t)twoTracksOneTube)/trackCount;
	cout << "Only one Hit per layer: " << oneHitPerLayer << ": (" << ohRel << "%)" << endl;
	cout << "More than 2 Hits per layer: " << moreThan3Hits << ": (" << mhRel << "%)" << endl;
	cout << "TwoTracksOneTube " << twoTracksOneTube << ": (" << mtRel << "%)" << endl;
}

map<Int_t,vector<PndFtsHit*>> PndAnalysisForwardTask::getMcTracks(){
	Int_t branchID = fIoman->GetBranchId("FTSHit"); //right hits
	//Int_t mcTrackBranchId = fIoman->GetBranchId("MCTrack"); //[R.K.03/2017] unused variable
	map<Int_t,vector<PndFtsHit*>> mcTracksToHits;
	mcToHitCount.clear();
	momentum.clear();
	//cout << "Found " << fHits->GetEntries() << " FTSHits " << endl;
	/*for(int i=0;i<fHits->GetEntries();i++){
		PndFtsHit* hit = (PndFtsHit*) fHits->At(i);
		hit->SetEntryNr(FairLink(-1, fIoman->GetEntryNr(), fIoman->GetBranchId("FTSHit"), i));
		FairLink link = hit->GetEntryNr();
		PndFtsPoint* p = (PndFtsPoint*) fIoman->GetCloneOfLinkData(link);
		FairMultiLinkedData mld = p->GetLinksWithType(mcTrackBranchId);
		FairLink linksToMc = mld.GetLink(0);
		int trackID = linksToMc.GetIndex();
		mcTracksToHits[trackID].push_back(hit);
	}*/


	for(int i=0;i<fIdealTrack->GetEntries();i++){
		PndTrack* t = (PndTrack*) fIdealTrack->At(i);
		Double_t mom = t->GetParamFirst().GetMomentum().Mag();
		momentum[i]=mom;
		//if(mom < 0.3) continue;
		PndTrackCand cand = t->GetTrackCand();
		set<FairLink> linksToHits = cand.GetLinksWithType(branchID).GetLinks();
		mcToHitCount[i]=linksToHits.size();
		for (std::set<FairLink>::iterator it = linksToHits.begin();it!=linksToHits.end();it++) {
			FairLink link = *it;
			PndFtsHit* hit = (PndFtsHit*) fIoman->GetCloneOfLinkData(link);
			mcTracksToHits[i].push_back(hit);
		}
	}

	//get the tracks with more than 6 hits
	/*map<Int_t,vector<PndFtsHit*>> result;
	for(map<Int_t,vector<PndFtsHit*>>::iterator i=mcTracksToHits.begin();i!=mcTracksToHits.end();i++){
		if(i->second.size()>=6){ //min 1 FTS-Station
			result[i->first] = i->second;
		}
		mcToHitCount[i->first]=i->second.size();
	}*/
	//if(result.size()>5) result.clear();

	return mcTracksToHits;
}

map<Int_t,vector<PndFtsHit*>*> PndAnalysisForwardTask::getReconstructedTracks(){
	map<Int_t,vector<PndFtsHit*>*> hitToTrack;
	cout << "Found " << fSolution->GetEntries() << " reconstructed Tracks" << endl;
	for(int i=0;i<fSolution->GetEntries();i++){
		//get the track
		PndTrack* track = (PndTrack*)fSolution->At(i);
		PndTrackCand* trackCand = track->GetTrackCandPtr();
		//get the hits
		vector<PndFtsHit*> *hits = new vector<PndFtsHit*>();
		std::set<FairLink> linksToHits = trackCand->GetLinksWithType(fIoman->GetBranchId("FTSHit")).GetLinks();
		for (std::set<FairLink>::iterator iter = linksToHits.begin();
				iter != linksToHits.end(); iter++) {
			FairLink link = *iter;
			PndFtsHit* hit = (PndFtsHit*) fIoman->GetCloneOfLinkData(link);
			hit->SetEntryNr(*iter);
			hits->push_back(hit);
		}
		//save the hits
		for(size_t j=0;j<hits->size();j++){
			PndFtsHit* hit = (*hits)[j];
			hitToTrack[hit->GetTubeID()]=hits;
		}
	}
	return hitToTrack;
}

void PndAnalysisForwardTask::analyzeMCTracks(vector<PndFtsHit*> hits,map<Int_t,vector<PndFtsHit*>*> tracks,Int_t trackID,Int_t trackCount){
	//get all reconstructed tracks with hits from this mc track
	vector<vector<PndFtsHit*>*> tracksToMc;
	cout << "Track has " << hits.size() << " Hits" << endl;
	for(size_t i=0;i<hits.size();i++){
		vector<PndFtsHit*> *track = tracks[hits[i]->GetTubeID()];
		if(track==0) continue;

		Bool_t exists = kFALSE;
		for(size_t j=0;j<tracksToMc.size();j++){
			if(tracksToMc[j]==track){
				exists=kTRUE;
				break;
			}
		}

		if(!exists)tracksToMc.push_back(track);
	}
	//all tracks
	//hits
	Int_t h2 = mcToHitCount[trackID];
	if(h2>48) hitArray2[8]++;
	else if(h2==48) hitArray2[7]++;
	else if(h2>40) hitArray2[6]++;
	else if(h2==40) hitArray2[5]++;
	else hitArray2[h2/8]++;
	//momentum
	Double_t m2 = momentum[trackID];
	if(m2>0.9) momentumArray2[9]++;
	else momentumArray2[(Int_t)(m2*10)]++;
	//check track count
	if(trackCount>9) trackCountArray2[9]++;
	else trackCountArray2[trackCount-1]++;

	cout << "Found " << tracksToMc.size() << " reconstructed tracks to MC-track" << endl;
	if(tracksToMc.size()==0){
		cout << ">>Track not found                            --! (0 Hits)" << endl;
		notFound++;
		//check hit corralation
		Int_t h = mcToHitCount[trackID];
		if(h>48) hitArray[8]++;
		else if(h==48) hitArray[7]++;
		else if(h>40) hitArray[6]++;
		else if(h==40) hitArray[5]++;
		else hitArray[h/8]++;
		//check momentum corralation
		Double_t m = momentum[trackID];
		if(m>0.9) momentumArray[9]++;
		else momentumArray[(Int_t)(m*10)]++;
		//check track count
		if(trackCount>9) trackCountArray[9]++;
		else trackCountArray[trackCount-1]++;
	} else {
		vector<Int_t> bestResult = analyzeMCTrack(hits,tracksToMc[0]);
		for(size_t i=1;i<tracksToMc.size();i++){
			vector<Int_t> nextResult = analyzeMCTrack(hits,tracksToMc[i]);
			if(bestResult[1]<nextResult[1]) bestResult = nextResult;
		}

		cout << "Mc-Hits: " << bestResult[0] << endl;
		cout << "Correct found: " << bestResult[1] << endl;
		cout << "Hits not found: " << bestResult[2] << endl;
		cout << "Wrong hits: " << bestResult[3] << endl;
		if(bestResult[0]==bestResult[1]){ //complete found
			if(bestResult[3]==0) {
				completeClean++;
				cout << ">>Found complete track, track is clean(++)"<< endl;
			} else {
				completeUnclean++;
				cout << ">>Found complete track, track is not clean(+-)"<< endl;
			}
		} else if(bestResult[0]*0.7 < bestResult[1]){ //fond more than 70%
			if(bestResult[3]==0) {
				incompleteClean++;
				cout << ">>Found incomplete track, track is clean(-+)"<< endl;
			} else {
				incompleteUnclean++;
				cout << ">>Found incomplete track, track is not clean(--)"<< endl;
			}
			cout << "Found Hits: (" << bestResult[1] << "/" << bestResult[0] <<")" << endl;
		} else {
			notEnoughFound++;
			cout << ">>Track not found(Under 70%)                        --!" << endl;
			cout << "Found Hits: (" << bestResult[1] << "/" << bestResult[0] <<")" << endl;
			//check hit corralation
			Int_t h = mcToHitCount[trackID];
			if(h>48) hitArray[8]++;
			else if(h==48) hitArray[7]++;
			else if(h>40) hitArray[6]++;
			else if(h==40) hitArray[5]++;
			else hitArray[h/8]++;
			//check momentum corralation
			Double_t m = momentum[trackID];
			if(m>0.9) momentumArray[9]++;
			else momentumArray[(Int_t)(m*10)]++;
			//check track count
			if(trackCount>9) trackCountArray[9]++;
			else trackCountArray[trackCount-1]++;
		}
	}
}

vector<Int_t> PndAnalysisForwardTask::analyzeMCTrack(vector<PndFtsHit*> mcHits,vector<PndFtsHit*>* reconstructedHits){
	Int_t mcHitCount = mcHits.size();
	Int_t correctHits = 0;
	Int_t hitsNotFound = 0;
	Int_t wrongHits = 0;

	//search for found mc hits
	for(size_t i=0;i<mcHits.size();i++){
		PndFtsHit *mcHit = mcHits[i];

		Int_t hitFound = kFALSE;
		for(size_t j=0;j<reconstructedHits->size();j++){
			PndFtsHit *recoHit = (*reconstructedHits)[j];
			if(mcHit->GetTubeID()==recoHit->GetTubeID()){
				hitFound = kTRUE;
				break;
			}
		}

		if(hitFound)correctHits++;
		else hitsNotFound++;
	}

	for(size_t i=0;i<reconstructedHits->size();i++){
		PndFtsHit* recoHit = (*reconstructedHits)[i];
		Int_t hitFound = kFALSE;
		for(size_t j=0;j<mcHits.size();j++){
			PndFtsHit *mcHit = mcHits[j];
			if(mcHit->GetTubeID()==recoHit->GetTubeID()){
				hitFound = kTRUE;
				break;
			}
		}
		if(!hitFound)wrongHits++;
	}

	vector<Int_t> result;
	result.push_back(mcHitCount);
	result.push_back(correctHits);
	result.push_back(hitsNotFound);
	result.push_back(wrongHits);
	return result;
}

void PndAnalysisForwardTask::anaSpecialCases(map<Int_t,vector<PndFtsHit*>> mcTracks){

	for(map<Int_t,vector<PndFtsHit*>>::iterator i=mcTracks.begin();i!=mcTracks.end();i++){
		Int_t layer[24];
		Int_t tubeID[10000];
		for(int j=0;j<24;j++) layer[j]=0;
		for(int j=0;j<10000;j++) tubeID[j]=0;
		for(size_t j=0;j<i->second.size();j++){
			PndFtsHit* hit = i->second[j];
			Int_t layerID = (hit->GetLayerID() - 1) / 2;
			layer[layerID]++;
			tubeID[hit->GetTubeID()]++;
		}

		for(int j=0;j<24;j++){
			if(layer[j]==1){
				oneHitPerLayer++;
				break;
			}
		}

		for(int j=0;j<24;j++){
			if(layer[j]>2){
				moreThan3Hits++;
				break;
			}
		}

		for(int j=0;j<10000;j++){
			if(tubeID[j]>1){
				twoTracksOneTube++;
			}
		}
	}
}