//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Track Cleaning. Remove Noise tracks.
//      
//
// Environment:
//      Software developed for the FOPI-TPC Detector at GSI.
//
// Author List:
//          Martin Berger          TUM     (original author)
//
//-----------------------------------------------------------

//This Class' Header
#include "TpcTrackCleaner.h"

//C/C++ Headers ---------------------------------------------
#include <typeinfo>
#include <vector>
#include <exception>

//Collaborating Headers -------------------------------------
#include "FairRootManager.h"
#include "FairRuntimeDb.h"
#include "FairRunAna.h"
#include "TClonesArray.h"
#include "GFTrack.h"
#include "GFException.h"
#include "TpcCosmics.h"

TpcTrackCleaner::TpcTrackCleaner()
:  fVerbose(kFALSE),
   fTrackArray(0),
   fDelTrackArray(0),
   fOutTrackArray(0),
   fInTrackArray(0),
   fTrackBranchName("TpcTrackPostFit"),
   fInBranchName("TpcTrackFitStat_0"),
   fDelTrackBranchName("DelTrack"),
   fOutBranchName("CutCosmics"),
   fThetaCut(kFALSE),
   fClMeanCut(kFALSE),
   fCosmics(kFALSE),
   fPersistence(kTRUE),
   fPhiCut(kFALSE),
   fLengthCut(kFALSE),
   fKeepDel(kFALSE),
   fThetaStart(0),
   fThetaEnd(-1),
   fPhiStart(0),
   fPhiEnd(-1),
   fClMeanStart(0),
   fClMeanStop(-1),
   fLengthStart(0),
   fLengthStop(-1)
{;}

TpcTrackCleaner::~TpcTrackCleaner()
{;}

InitStatus TpcTrackCleaner::Init()
{
	//Get ROOT Manager ---------------------
	FairRootManager * ioman = FairRootManager::Instance();
	if (ioman==0)
	{
		Error("TpcReClsuterizerTask::Init","Root Manager not instantiated!");
		return kERROR;
	}

	if (!fCosmics)
	{
		fTrackArray=(TClonesArray*) ioman->GetObject(fTrackBranchName);
		if (fTrackArray==0)
		{
			Error("TpcReClusterizerTask::Init","Track Array not found");
		}
	}

	//create and register new output arrays
	if (fCosmics)
	{
		fInTrackArray=(TClonesArray*) ioman->GetObject(fInBranchName);
		if (fInTrackArray==0)
		{
			std::cout<<"track array: "<<fInBranchName<<"\n";
			Fatal("Init","Track input Array not found");
		}
		if (fKeepDel)
		{
			fDelTrackArray=new TClonesArray("TpcCosmics");
			ioman->Register(fDelTrackBranchName,"Tpc",fDelTrackArray,fPersistence);
		}

		fOutTrackArray=new TClonesArray("TpcCosmics");
		ioman->Register(fOutBranchName,"Tpc",fOutTrackArray,fPersistence);
	}
	return kSUCCESS;
}

void TpcTrackCleaner::SetParContainers()
{
	;
}

void TpcTrackCleaner::Exec(Option_t * opt)
{
	if (fCosmics)
		ExecCosmics(opt);
}

void TpcTrackCleaner::ExecCosmics(Option_t* opt)
{
	if(fKeepDel)
	{
		if (fDelTrackArray==0) Fatal("Exec","No deleted track array");
		fDelTrackArray->Delete();
	}

	if (fOutTrackArray==0) Fatal("Exec","No Output Array");
	fOutTrackArray->Delete();

	int nTracks=fInTrackArray->GetEntriesFast();
	if(fVerbose)
		std::cout<<"Start loop over "<<nTracks<<" tracks\n";
	for (int nTr=0;nTr<nTracks;nTr++)
	{
		TpcCosmics* cos=(TpcCosmics*)fInTrackArray->At(nTr);

		if (cos->GetStatFlag()!=0)
		{
			if (fVerbose)
				std::cout<<"Bad TrackFlag\n";
			if (fKeepDel)
			{
				unsigned int nBC = fDelTrackArray->GetEntriesFast();
				TpcCosmics* badguy = new ((*fDelTrackArray)[nBC]) TpcCosmics(*cos);
			}
			continue;
		}

		double theta = cos->GetTheta();
		if (fThetaCut)
			if(theta<fThetaStart || (fThetaEnd!=-1 && std::fabs(theta)>fThetaEnd))
			{
				if (fVerbose)
					std::cout<<"Bad Theta angle\n";
				if (fKeepDel)
				{
					unsigned int nBC = fDelTrackArray->GetEntriesFast();
					TpcCosmics* badguy = new ((*fDelTrackArray)[nBC]) TpcCosmics(*cos);
				}
				continue;
			}

		double phi = cos->GetPhi();
		if(fPhiCut)
			if(std::fabs(phi)<fPhiStart || ( fPhiEnd!=-1 && std::fabs(phi)>fPhiEnd))
			{
				if (fVerbose)
					std::cout<<"Bad Phi angle\n";
				if (fKeepDel)
				{
					unsigned int nBC = fDelTrackArray->GetEntriesFast();
					TpcCosmics* badguy = new ((*fDelTrackArray)[nBC]) TpcCosmics(*cos);
				}
				continue;
			}

		double clmeanamp = cos->GetClMeanAmp();
		if(fClMeanCut)
			if(clmeanamp<fClMeanStart || (fClMeanStop!=-1 && clmeanamp>fClMeanStop ) )
			{
				if (fVerbose)
				{
					std::cout<<"Bad Cluster Mean Amp\n";
					std::cout<<"Mean Amp: "<<clmeanamp<<" (min="<<fClMeanStart<<", max="<<fClMeanStop<<")\n";
					std::cout<< (clmeanamp<fClMeanStart) <<" "<<(fClMeanStop!=-1 && clmeanamp>fClMeanStop )<<"\n";
				}
				if(fKeepDel)
				{
					unsigned int nBC = fDelTrackArray->GetEntriesFast();
					TpcCosmics* badguy = new ((*fDelTrackArray)[nBC]) TpcCosmics(*cos);
				}
				continue;
			}

		double tracklen = cos->GetLength();
		if(fLengthCut)
			if(tracklen<fLengthStart ||(fLengthStop!=-1 && tracklen>fLengthStop) )
			{
				if (fVerbose)
				{
					std::cout<<"Bad Tracklengtyh\n";
					std::cout<<"Length="<<tracklen<<" (min="<<fLengthStart<<", max="<<fLengthStop<<")\n";
				}
				if(fKeepDel)
				{
					unsigned int nBC = fDelTrackArray->GetEntriesFast();
					TpcCosmics* badguy = new ((*fDelTrackArray)[nBC]) TpcCosmics(*cos);
				}
				continue;
			}

		unsigned int nGC = fOutTrackArray->GetEntriesFast();
		if (fVerbose)
			std::cout<<"Good track\n";
		TpcCosmics * goodguy = new ((*fOutTrackArray)[nGC]) TpcCosmics(*cos);


	}//end loop over tracks

}

void TpcTrackCleaner::ExecTrack(Option_t* opt)
{
	if (fDelTrackArray==0) Fatal("TpcTrackCleaner::Exec","No deleted track array");
	fDelTrackArray->Delete();

	if (fOutTrackArray==0) Fatal("Exec","mo Output Array");
	fOutTrackArray->Delete();

	Int_t ntracks=fTrackArray->GetEntriesFast();
	if (ntracks>20000)
	{
		std::cout<<"Track Cleaner: ntracks="<<ntracks<<" Evil Event! skipping"<<std::endl;
		return;
	}
	if(fVerbose)
		std::cout<<"Track Cleaning: Looping over "<<ntracks<<" tracks"<<std::endl;

	for (Int_t itr=0;itr<ntracks;itr++)
	{
		GFTrack* trk=(GFTrack*)fTrackArray->At(itr);
		GFAbsTrackRep* theRep = trk->getCardinalRep();
		if (theRep->getStatusFlag() != 0)
		{
			if (fVerbose)
				std::cout<<"Track Cleaner: Track anyway not fitted. going to next one\n";
			continue;
		}

		// get pos and momentum for start position and direction
		TVector3 pos,mom;
		try
		{
			pos = theRep->getPos();
			mom = theRep->getMom();
		}
		catch(GFException& e)
		{
			e.what();
			if (fVerbose) std::cout<<"Track Cleaner: can't get mom and pos of track\n";
			continue;
		}

		//check if track theta angle is inside cuts
		if (fThetaCut)
		{
			double Theta=mom.Theta();
			if(Theta>fThetaStart && Theta<fThetaEnd)
			{
				;
			}
			else
			{
				if (fVerbose) std::cout<<"Track Cleaner: track is inside theta cut\n";
				theRep->setStatusFlag(10);
			}
		}
	}

}