/*
 * PndIsochroneTrackFinder.h
 *
 *  Created on: 10.10.2014
 *      Author: Stockmanns
 */

#ifndef PNDISOCHRONETRACKFINDER_H_
#define PNDISOCHRONETRACKFINDER_H_

#include <set>
#include <map>
#include <vector>
#include <utility>

#include "TObject.h"
#include "TH2D.h"
#include "TGraph.h"

#include "FairHit.h"

#include "PndTrack.h"
#include "PndTrackCand.h"
#include "PndTrackFromCircle.h"
#include "PndSttHit.h"


#include "TClonesArray.h"
class TVector2;

class PndGeneralHoughHit : public FairHit
{
public:
	PndGeneralHoughHit(){};
	~PndGeneralHoughHit(){};
	PndGeneralHoughHit(const FairHit& hit) : FairHit(hit), isochroneRadius(0){};

	void SetIsochroneRadius(Double_t val) {isochroneRadius = val;}
	Double_t GetIsochroneRadius() const {return isochroneRadius;}
	TVector2 GetXYVector(){return TVector2(GetX(), GetY());}
private:
	Double_t isochroneRadius;
};

struct PndCircleTime {
	TVector2 fCircle;
	Double_t fTime;
	Double_t fRRange;
	Double_t fPhiRange;
};

struct PndPeakContent {
	std::set<Int_t> fHitsInPeaks;
	std::set<Int_t> fPeakBins;
	TVector2 fCircleCenter;
	std::pair<Double_t, Double_t> fRPhiRange;
	void AddHits(std::set<Int_t> hits);
	void AddPeak(Int_t peak);
	void AddPeaks(std::set<Int_t> peaks);
	void AddPeakContent(const PndPeakContent& peak);
	TVector2 GetMeanPeakPos() const {return fCircleCenter;}
	std::pair<Double_t, Double_t> GetPeaksRangeRPhi() const {return fRPhiRange;}
	TVector2 FillMeanPeakPos(TH2* h2);
	std::pair<Double_t, Double_t> FillPeaksRangeRPhi(TH2* h2);
	virtual std::ostream& Print(std::ostream& out = std::cout) const {
		out << "Hits: /";
		for (std::set<Int_t>::iterator hitIter = fHitsInPeaks.begin(); hitIter != fHitsInPeaks.end(); hitIter++){
			out << *hitIter << "/";
		}
		out << " Peaks: /";
		for (std::set<Int_t>::iterator peakIter = fPeakBins.begin(); peakIter != fPeakBins.end(); peakIter++){
			std::cout << *peakIter << "/";
		}
		out << " CircleCenter: " << fCircleCenter.X() << "/" << fCircleCenter.Y() << " Mod: " << fCircleCenter.Mod() << " Phi: " << fCircleCenter.Phi()
				<< " RRange: " << fRPhiRange.first << " PhiRange " << fRPhiRange.second << std::endl;
		return out;
	}

	friend std::ostream& operator<< (std::ostream& out, const PndPeakContent& peak){
		peak.Print(out);
		return out;
	}
};

class PndIsochroneTrackFinder : public TObject {
public:
	PndIsochroneTrackFinder();
	virtual ~PndIsochroneTrackFinder();

	void SetBField(Double_t field) {fBField = field;}
	void SetTubeArray(TClonesArray* val){fSttTubeArray = val;}
	void SetAngleRange(std::pair<Double_t, Double_t> range) {fAngleRange = range;}
	void SetAngleSteps(Double_t stepSize) {fAngleSteps = stepSize;}
	void SetPtRange(std::pair<Double_t, Double_t> range) {fPtRange = range;}
	void SetPtStepSize(Double_t stepSize) {fPtStepSize = stepSize;}

	void AddHits(TClonesArray* hits, Int_t branchId);
	std::vector<PndTrack> FindTracks();

	std::vector<PndTrack> FindMvdBasedTracks();

	TH2D* GetHoughHisto(){return fHoughHisto;}
	void SetSkewedBranchName(TString branchName){fSkewedBranchName = branchName;}
	TString GetSkewedBranchName() const {return fSkewedBranchName;}
	void SetSkewedArray (TClonesArray* val){ fSttSkewedArray = val; }
	TClonesArray* GetSkewed() const {return fSttSkewedArray;}
	std::vector<PndTrackCand> GetTrackCands(){return fTrackCands;}
	std::vector<TVector2> GetPeakCoordinates(){return fPeakCoordinates;};
	TGraph GetPeakCoordinateGraph();

	std::vector<TVector2> GetPeakCoordinatesMvd(){return fPeakCoordinatesMvd;};
	TGraph GetPeakCoordinateGraphMvd();
	std::map<Double_t, FairHit>::iterator CheckContinuity(std::map<Double_t, FairHit>& cand, Double_t circleMag);
	Bool_t CheckTrackSimilarity(std::set<Int_t> first, std::set<Int_t> second, Double_t threshold);
	Bool_t CheckPeakDistance(std::set<Int_t> peakBins, Int_t binToCheck, Double_t threshold);

	void AssignSkewedHits(PndTrackCand& cand, std::vector<std::vector<FairHit> > skewedHits, TVector2 circle);

	void Reset(){
		if (fHoughHisto != 0)
			fHoughHisto->Reset();
		fBinsWithContent.clear();
		for (std::map<Int_t, std::set<Int_t> >::iterator iter = fBinToHit.begin(); iter != fBinToHit.end(); ++iter){
			iter->second.clear();
		}
		fBinToHit.clear();
		fHitsUsed.clear();
		fTrackCands.clear();
		fHitsToFit.clear();
		fSkewedHits.clear();
		if (fSttSkewedArray != 0)
			fSttSkewedArray->Delete();
		fPeakCoordinates.clear();
		fPeakCoordinatesMvd.clear();

	}

protected:
	void AddHitsToFit(TClonesArray* hits, Int_t branchId);
	TVector2 IsoCalc(PndGeneralHoughHit& houghHit, Double_t phi);
	std::vector<Double_t> CalcPhiForPt(PndGeneralHoughHit& hit, Double_t pt);
	std::pair<Double_t, Double_t> CalcPhiForT(Double_t t, TVector2 hit);
	std::vector<Double_t> CalcPhiForCircleTimes(TVector2& hit, std::vector<PndCircleTime>& circleTimes);
	std::vector<Double_t> CalcPhiForPeakContents(TVector2& hit, std::vector<PndPeakContent>& peakContents);

	Int_t GetHistoBin(TH2* histo, Double_t xVal, Double_t yVal);
	void InterpolateHisto(TH2* histo, TVector2& in, TVector2& out, Int_t& hitId, Double_t weight = 1.);
	std::multimap<Int_t, Int_t> FindPeaks(std::set<Int_t>& binsWithContent, TH2* histo, Int_t threshold); ///< returns multimap with number of hits in bin and bin number
	std::vector<PndPeakContent> MergePeaks(std::multimap<Int_t, Int_t> peaks, Double_t trackSimilarityCut,  Int_t lowerHitLimit, Bool_t separateTracks = kFALSE);
	PndTrackCand FillTrackCand(Int_t histoBin);
	PndTrackCand FillTrackCand(PndPeakContent& peakContent);
	Double_t CalcArcLength(FairHit hit, TVector2 circleCenter);
	PndTrack CalcPndTrack(PndTrackCand& cand, TVector2 circle);
	FairTrackParP CalcTrackParP(FairHit hit, TVector2 circleCenter);
	TVector3 CalcHitPosInTrack(FairHit hit, TVector2 circleCenter);
	Double_t CalcPt(TVector2 circleCenter);
	TVector2 CalcPtDir(TVector2 hit, TVector2 circleCenter);
	std::vector<TVector2> CalcCircleIntersections(TVector2 circleCenter1, Double_t radius1, TVector2 circleCenter2, Double_t radius2);
	//Bool_t IsEqual(std::set<Int_t> a, std::set<Int_t> b);
	Bool_t ContainsAB(std::set<Int_t> a, std::set<Int_t> b);
	std::vector<PndPeakContent> SeparateTracks(std::set<Int_t> hitsInBin, TVector2 circle, Int_t peakBin);

	std::vector<Double_t> CalcArcLengthRegions(Double_t radius);

	std::vector< std::vector <FairHit> > FindSkewedHits(TVector2 circle);
	std::vector<TVector2> GenerateHoughValues(PndGeneralHoughHit& houghHit);
	std::vector<TVector2> GenerateHoughValuesPhi(PndGeneralHoughHit& houghHit, std::set<Double_t> phiValues);
	std::vector<Double_t> GeneratePtValues();
	std::set<Double_t> GeneratePhiValues(PndGeneralHoughHit& myHit);
	void FillHoughHist(std::vector<TVector2>, Int_t, Double_t);
	void CalcHoughHist(PndGeneralHoughHit& hit, Int_t hitIndex, Bool_t plusIsochrone);
	std::vector<PndTrack> PeaksToTracks(std::multimap<Int_t, Int_t>, Int_t lowerHitLimit);
	std::vector<PndTrack> PeaksToTracks(std::vector<PndPeakContent>& mergedPeaks, Int_t lowerHitLimit);
	PndCircleTime FindCircleTime(Int_t peakBins);
	Double_t CalcTimeStampInBin(Int_t bin);

	std::vector<PndPeakContent> MergePeaksByPos(std::vector<PndPeakContent>& peaks, Double_t peakDistanceCut);

	Double_t RToPt(Double_t r) {
		return 0.3 * fBField * r * 0.01;
	}

	Double_t PtToR(Double_t pt){
		return pt * 100/(0.3 * fBField);
	}


	//void GetStartEndPointOfTube(PndSttHit*, TVector3& startPoint, TVector3& endPoint);

//	void FillHistogram();

private:
	Double_t fBField;
	std::pair<Double_t, Double_t> fAngleRange;
	Double_t fAngleSteps;
	std::pair<Double_t, Double_t> fPtRange;  // range of investigated pt momentum, usually 0 to 15 GeV/c (all possible momenta a track can have)
	Double_t fPtStepSize;  // the step width with which the momentum range is rastered through
	TH2D* fHoughHisto;
	std::set<Int_t> fBinsWithContent; //!
	std::vector<PndGeneralHoughHit> fHitsToFit;
	std::map<Int_t, std::set<Int_t> > fBinToHit; //Map of BinInHisto(key) to index in fHitsToFit
	std::set<Int_t> fHitsUsed;
	Bool_t fClockwise;
	std::vector<PndTrackCand> fTrackCands;
	std::map<TString, Double_t> fWeightMap;
	TClonesArray* fSttTubeArray;
	std::vector<PndSttHit*> fSkewedHits;
	std::map<FairLink, PndGeneralHoughHit> fMapFairLinkHit;
	std::vector<TVector2> fPeakCoordinates;
	std::vector<TVector2> fPeakCoordinatesMvd;
	std::multimap<Int_t, std::pair<PndTrackCand, TVector2> > fSeparatedTracks; ///< NumberOfHitsInTrack, <TrackCand, Circle> >

	PndTrackFromCircle fTrackCalc;


	TString fSkewedBranchName;				//<Output branchName for the calculated hits from the skewed layers
	TClonesArray* fSttSkewedArray;

	TClonesArray* fMvdHitsPixelArray;
	TClonesArray* fMvdHitsStripArray;
	TClonesArray* fSttHitArray;
	TClonesArray* fGemHitArray;

	Int_t fMvdHitsPixelBranchId;
	Int_t fMvdHitsStripBranchId;
	Int_t fSttHitBranchId;
	Int_t fGemHitBranchId;

	Int_t fHitsToFitAfterMvd;	//< first hit after MVD hits at the beginning of the fHitsToFit vector

	Double_t fCutMergePeakDistance;

public:
	ClassDef(PndIsochroneTrackFinder, 1);
};

#endif /* PNDISOCHRONETRACKFINDER_H_ */