//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Reclustering. Redo clustering with the help od track information
//
//
// Environment:
//      Software developed for the FOPI-TPC Detector at GSI.
//
// Author List:
//          Martin Berger          TUM     (original author)
//
//-----------------------------------------------------------


// This Class' Header ---------------------------------------
#include "TpcReClusterizerTask.h"

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

//Collaborating Headers -------------------------------------
#include "FairRootManager.h"
#include "FairRuntimeDb.h"
#include "FairRunAna.h"
#include "TClonesArray.h"
#include "TpcPrelimCluster.h"
#include "TpcCluster.h"
#include "TpcSPHit.h"
#include "TpcRiemannTrack.h"
#include "TpcRiemannHit.h"
#include "GFTrack.h"
#include "GFAbsTrackRep.h"
#include "GFException.h"
#include "RKTrackRep.h"
#include "TpcDigiMapper.h"
#include "TpcAlignmentManager.h"
#include "TpcDigiPar.h"
#include "TpcGas.h"
#include "TpcPadPlane.h"
#include "TH2D.h"
#include "TRandom.h"
#include "TStopwatch.h"
#include "McIdCollection.h"
#include "McId.h"

TpcReClusterizerTask::TpcReClusterizerTask()
: fDigiBranchName("TpcDigi"),
  fSPHitBranchName("TpcSPHit"),
  fClusterBranchName("TpcCluster"),
  fReClusterBranchName("TpcReCluster"),
  fTrackBranchName("RiemannTrack"),
  fRecoHitBranchName("TpcSPHit"),
  fOutTrackBranchName("TpcTrackPreFit"),
  fDigiArray(0),
  fSPHitArray(0),
  fClusterArray(0),
  fReClusterArray(0),
  fRecoHitArray(0),
  fTrackArray(0),
  fOutTrackArray(0),
  fHitBranchID(-1),
  fhTrack(),
  fTimeinfo(kFALSE),
  fPersistence(kFALSE),
  fArhpersistence(kFALSE),
  fstepsize(1),
  fCounter(0),
  fOnDigi(kTRUE),
  fIgnoreEdgeDigis(kFALSE),
  fGlobCounter(0),
  fUseFirstDigiPos(kFALSE),
  fUseChamberEdge(kFALSE),
  fUseCosmics(kFALSE),
  fRDigiCutoff(1),
  fMinClusterSize(5),
  fVarStep(kFALSE),
  fVarStepSqrt(kFALSE),
  fStepOff(15),
  fStepSlope(0.13535),
  fStepOff_sqrt(1),
  fStepSlope_sqrt(0.5),
  fStepOff_0_sqrt(0),
  fFac(4),
  fZoff(0),
  fAlign(kTRUE),
  fTrAlign(kTRUE),
  fUseAllDigis(kFALSE),
  fVarMinClusterSize(kFALSE),
  fTrackInit(kFALSE),
  fDevMapName(""),
  fClusterCorr(kFALSE),
  fSmallClusterCorr(kFALSE),
  fSkipUnFitted(kFALSE),
  fuseFullCov(kFALSE),
  fweightedPlaneConstruction(kFALSE),
  fcutCov(kFALSE),
  frecalcDigiPos(kFALSE),
  fDigiSharing(kTRUE),
  fOnlyStep(kFALSE),
  fUsePRF(kTRUE),
  fUseSPHitParam(kFALSE),
  fDoPRF(kTRUE),
  fC(14.),
  fG(600),
  fzJitter(0.1),
  first(true),
  fPlane(),
  fZMax(-1),
  fZMin(-1),
  fuseOnPlaneCov(false),
  fPar(NULL),
  fcorrFactor(1.),
  fMCS1(0),
  fMCS2(0),
  fMCS3(0),
  fPrelimClusterCounter(0)
{
	fcorrector=new TpcClusterCorrector();
}

TpcReClusterizerTask::~TpcReClusterizerTask()
{;}

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

	fDigiArray=(TClonesArray*) ioman->GetObject(fDigiBranchName);
	if(fDigiArray==0)
	{
		Fatal("TpcReClusterizerTask::Init","Digi Array not found");
		return kERROR;
	}

	fSPHitArray=(TClonesArray*) ioman->GetObject(fSPHitBranchName);
	if(fSPHitArray==0)
	{
		Fatal("TpcReClusterizerTask::Init","SPHit Array not found");
		return kERROR;
	}

	// Get Input Arrays --------------------
	fClusterArray=(TClonesArray*) ioman->GetObject(fClusterBranchName);
	if (fClusterArray==0)
	{
		Fatal("TpcReClusterizerTask::Init","Cluster Array not found");
		return kERROR;
	}
	fTrackArray=(TClonesArray*) ioman->GetObject(fTrackBranchName);
	if (fTrackArray==0)
	{
		Fatal("TpcReClusterizerTask::Init","Track Array not found");
	}

	// Create and register output arrays
	fReClusterArray=new TClonesArray("TpcCluster");
	if (fReClusterArray->GetEntries() !=0)
		Fatal("TpcReClusterizerTask::Init","Cluster Array was not empty");
	ioman->Register(fReClusterBranchName,"Tpc",fReClusterArray,fPersistence);

	fRecoHitArray = new TClonesArray("TpcSPHit");
	if (fRecoHitArray->GetEntries() !=0)
		Fatal("TpcReClusterizerTask::Init","SPHit Array was not empty");
	ioman->Register(fRecoHitBranchName,"tpc",fRecoHitArray,fPersistence);

	fOutTrackArray = new TClonesArray("GFTrack");
	if (fOutTrackArray->GetEntries() !=0)
		Fatal("TpcReClusterizerTask::Init","Output track Array was not empty");
	ioman->Register(fOutTrackBranchName,"GenFit",fOutTrackArray,fPersistence);

	fHitBranchID = FairRootManager::Instance()->GetBranchId(fSPHitBranchName);
	fSPHitID = FairRootManager::Instance()->GetBranchId(fRecoHitBranchName);
	fClBrID = FairRootManager::Instance()->GetBranchId(fReClusterBranchName);

	fPlane = fPar->getPadPlane();
	fGas  = fPar->getGas();
	fDriftVel=fPar->getManDriftVel();
	if(fDriftVel<=0){
		fDriftVel = fGas->VDrift();
	}
	fRMin = fPar->getRMin();
	fRMax = fPar->getRMax();
	fZMax = fPar->getZMax();
	//  fZMin = fPar->getZMin();

	if (fClusterCorr)
	{
		if(fDevMapName=="")
		{
			Error("Init","Reading DevMap from parbase\n");
			fDevMapName = fPar->getClustCorrFile();
		}
		std::cout<<"TpcReClusterizerTask::Init() using devmap: "<<fDevMapName.Data()<<"\n";
		fDevMap = new TpcDevmapCyl(fDevMapName,kTRUE);
		if (!fDevMap->loaded())
		{
			Fatal("Init","Deviation Map not found!");
			return kERROR;
		}
		fDevMap->printc();
	}

	if(fSmallClusterCorr)
	{
		fcorrector->SetDigiArray(fDigiArray);
		fcorrector->SetPadPlane(fPlane);
	}

	if(fIgnoreEdgeDigis) {
		fRMin = fRCutoff1;
		fRMax = fRCutoff2;
	}
	return kSUCCESS;
}

void TpcReClusterizerTask::SetParContainers()
{
	if (fVerbose)
		std::cout<<"TpcReClusterizerTask::SetParContainers"<<std::endl;

	// Get run and runtime database
	FairRun* run = FairRun::Instance();
	if ( ! run ) Fatal("SetParContainers", "No analysis run");

	FairRuntimeDb* db = run->GetRuntimeDb();
	if ( ! db ) Fatal("SetParContainers", "No runtime database");

	// Get Tpc digitisation parameter container
	fPar = (TpcDigiPar*) db->getContainer("TpcDigiPar");
	if (! fPar ) Fatal("SetParContainers", "TpcDigiPar not found");
}

void TpcReClusterizerTask::Exec(Option_t* opt)
{
	if(first)
	{
		//this block is to define the z jitter
		McId dummyID(1,1,0);
		McIdCollection dummyColl;
		dummyColl.AddID(dummyID);

		TVector3 zDiff1,zDiff2;
		TpcDigi zDiffDigi1(1,1,1,dummyColl),zDiffDigi2(1,2,1,dummyColl);
		try {
			TpcDigiMapper::getInstance()->map(&zDiffDigi1,zDiff1);
			TpcDigiMapper::getInstance()->map(&zDiffDigi2,zDiff2);
		}
		catch(const std::exception& ex) {
			std::cout<<ex.what()<<std::endl;
		}
		fzJitter = zDiff2.z() - zDiff1.z();
		std::cout<<"zJitter:"<<fzJitter<<"\n";
		//end of z jitter
		first=false;
	}

	if (fVerbose)
		std::cout<<"TpcReClusterizerTask::Exec; Event Number: "<<fCounter++<<std::endl;

	if (fReClusterArray==0) Fatal("TpcReClusterizerTask::Exec","No Cluster out array");
	fReClusterArray->Delete();

	if (fTrackInit)
	{
		if (fOutTrackArray==0) Fatal("TpcReClusterizerTask::Exec","No Track out array");
		fOutTrackArray->Delete();

		if (fRecoHitArray==0) Fatal("TpcReClusterizerTask::Exec","No RecoHit out array");
		fRecoHitArray->Delete();
	}

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

	//Loop over tracks and recluster for every track
	//double gain=fPar->getGain();

	fGlobCounter=0;
	fPrelimClusterCounter=0;

	for (Int_t itr=0;itr<ntracks;itr++)
	{
		// get the track
		if (dynamic_cast<TpcRiemannTrack*>(fTrackArray->At(itr))!=NULL)
		{
			TpcRiemannTrack * trk=(TpcRiemannTrack*)fTrackArray->At(itr);
			RiemannReCluster(trk);
		}
		else
		{
			if(fVerbose)
				std::cout<<"Using genfit tracks for reclustering"<<std::endl;
			GFTrack* trk=(GFTrack*)fTrackArray->At(itr);
			if(fOnlyStep)
				GenfitStepThrough(trk);
			else
				GenfitReCluster(trk);
		}

	} // done loop over tracks
}

Bool_t TpcReClusterizerTask::RiemannReCluster(TpcRiemannTrack* ftrack)
{
	std::cout<<"Warning RiemanTrack based Reclustering is not implemented!"<<std::endl;
	//get number of hits in track
	Int_t nhits=ftrack->getNumHits();
	//check if track is prefitted
	if ( ! ftrack->isFitted())
	{
		if ( fVerbose) std::cout<<" TpcReClusterizerTask::Exec -skipping, track not prefitted"<<std::endl;
		return kFALSE;
	}
	//Loop over the hits
	for ( Int_t ihit=0;ihit<nhits;ihit++)
	{
		TpcRiemannHit* hit = ftrack->getHit(ihit);
	}//done loop over hits
	return kTRUE;
}

Bool_t TpcReClusterizerTask::GenfitReCluster(GFTrack* trk)
{
	TStopwatch timer;
	timer.Start();
	if(fTimeinfo)
		PrintTime(timer,"START");
	FairRootManager * ioman = FairRootManager::Instance();
	if(fVerbose)
		std::cout<<"****************************GENFIT RECLUSTER******************"<<std::endl;
	GFAbsTrackRep* theRep = trk->getCardinalRep();
	if (theRep->getStatusFlag() != 0 && fSkipUnFitted)
	{
		if (fVerbose)
			std::cout<<"ReCl: track not fitted!\n";
		return kFALSE;
	}

	// //check for GEANE trackrep
	// if(dynamic_cast<GeaneTrackRep*>(theRep) != NULL)
	//   {
	//     ((GeaneTrackRep*)theRep)->setPropDir(0); // not needed for RKTrackRep
	//   }

	GFTrackCand cand = trk->getCand();

	bool useSPHits = true;
	std::vector<int> hitIDs = cand.getHitIDs(fHitBranchID);

	// get pos and momentum for start position and direction
	TVector3 pos,mom;
	if (fVerbose) std::cerr<<"get pos and mom"<<std::endl;
	try
	{
		pos = theRep->getPos();
		mom = theRep->getMom();
	}
	catch(GFException& e)
	{
		e.what();
		return kFALSE;
	}

	if (fVerbose)
	{
		std::cout<<"pos: ("<<pos.X()<<","<<pos.Y()<<","<<pos.Z()<<")\n";
		std::cout<<"mom: ("<<mom.X()<<","<<mom.Y()<<","<<mom.Z()<<")\n";
	}
	TVector3 startDir = mom;
	startDir.SetMag(1.);

	TVector3 altFirstPos;
	TVector3 altFirstDir;
	//set reference plane:
	GFDetPlane pl;
	pl.setO(pos);
	pl.setNormal(startDir);
	Int_t cosmics=1;
	//  if (pos*startDir<0)
	//    cosmics=-1;

	if(fTimeinfo)
		PrintTime(timer,"Before Collecting digis");

	if (fVerbose) std::cerr<<"collect digis"<<std::endl;
	// get all digis in the track
	std::map<unsigned int, TpcDigi*> digis;   //<TCA index, pointer>
	unsigned int count = 0;
	TClonesArray* actClArr;
	int foundDigi = -1;
	double maxDist = 0.;  //(cm)

	if (fVerbose) std::cout<<"Cluster in this track: "<<hitIDs.size()<<std::endl;
	//looping over all cluster from the track to collect all digis
	for(int i=0; i<hitIDs.size(); i++)
	{
		TpcCluster* cl;
		TpcSPHit* theHit = (TpcSPHit*) fSPHitArray->At(hitIDs[i]);

		int clBrID = theHit->getClusterBranchID();
		int clID = theHit->getClusterID();
		if (count==0)
		{
			TString clBrName = ioman->GetBranchName(clBrID);
			actClArr = (TClonesArray*) ioman->GetObject(clBrName);
		}

		cl = (TpcCluster*)(actClArr->At(clID));
		count++;

		const std::map<unsigned int,double>* digiMap = cl->getDigiMap();

		for(std::map<unsigned int, double>::const_iterator it=digiMap->begin();it!=digiMap->end();it++)
		{
			// avoid double counting
			if(!digis.count(it->first))
			{
				digis[it->first] = (TpcDigi*) fDigiArray->At(it->first);
			}
			else
			{
				//	      std::cout<<"double counted digi:"<<it->first<<" "<<it->second<<"\n";
				continue;
			}
			TpcDigi* aDigi = digis[it->first];

			//only for first hit (cluster): find first Digi
			if( fUseFirstDigiPos)
			{
				//get Digi coordinates:
				TVector3 dPos;
				try
				{
					TpcDigiMapper::getInstance()->map(aDigi,dPos);
				}
				catch (...)
				{
					std::cerr<<"TpcReclusterizer::Could not get Position of digi! Skipping ..."<<std::endl;
					continue;
				}
				//if (fAlign)
				//dPos = TpcAlignmentManager::getInstance()->localToMaster("tpc",dPos);
				TVector3 pc, dpc;
				try
				{
					theRep->extrapolateToPoint(dPos,pc,dpc);
				}
				catch(GFException& e)
				{
					e.what();
					continue;
				}
				TVector3 dDist = pl.getO() - pc;
				if (fVerbose)
				{
					std::cout<<"poca: "<<pc.X()<<";"<<pc.Y()<<";"<<pc.Z()<<"\n";
					std::cout<<"poca Dist: "<<dDist.Mag()<<"\n";
					std::cout<<"dirdingens: "<<dDist*pl.getNormal()<<"\n";
				}
				if(dDist.Mag() >= maxDist && (dDist*pl.getNormal()) >= 0.)   //we want negative extrap!
				{
					if(fIgnoreEdgeDigis)
						if(dPos.Perp() < fRCutoff1 || dPos.Perp() > fRCutoff2)
							continue;
					foundDigi = (*it).first;
					altFirstPos = dPos;
					altFirstDir = dpc;
					maxDist = dDist.Mag();
					if(fVerbose)
						std::cout<<"Moving start position\n";
				}
			}
		}
	}
	if (fTimeinfo)
		PrintTime(timer,"Collected digis");

	// collect all unused digis
	std::map<unsigned int, TpcDigi*> udigis;   //<TCA index, pointer>
	if(fUseAllDigis)
	{
		for (Int_t idigi=0;idigi<fDigiArray->GetEntries();idigi++)
		{
			TpcDigi* auDigi = (TpcDigi*) fDigiArray->At(idigi);
			// avoid double counting
			if(!digis.count(auDigi->index()) && !udigis.count(auDigi->index()) )
			{
				//std::cout<<"here is another array"<<std::endl;
				udigis[auDigi->index()] = (TpcDigi*) auDigi;
			}
			else
				continue;
		}
		if(fTimeinfo)
			PrintTime(timer,"Collected unused digis");
	}

	if (fVerbose)
	{
		std::cout<<digis.size()<<" digis to process from cluster"<<std::endl;
		TVector3 pos0(((TpcCluster*)(fClusterArray->At(hitIDs[0])))->pos());
		std::cout<<"pos - pos0 mag"<< (pos-pos0).Mag() << std::endl;
	}


	GFDetPlane here,next;
	mom.SetMag(1);
	here.setO(pos);
	here.setNormal(mom);
	if(fVerbose)
	{
		std::cout<<"First Plane:"<<pos.X()<<":"<<pos.Y()<<":"<<pos.Z()<<"\n";
		std::cout<<"First Direction:"<<mom.X()<<":"<<mom.Y()<<":"<<mom.Z()<<"\n";
	}

	//if we successfully found a "first digi" use it as starting pos
	if (foundDigi >= 0 && fUseFirstDigiPos)
	{
		if(fVerbose)
		{
			std::cout<<"Moved plane from cluster position (r="<<here.getO().Perp()<<"): \n    " ;
			here.getO().Print();
			std::cout.flush();
			std::cout<<"\nto first Digi position (r="<<altFirstPos.Perp()<<"): \n    " ;
			altFirstPos.Print();
			std::cout.flush();
			std::cout<<"\n"<<std::endl;
		}
		startDir=-(altFirstPos-here.getO());
		//double pull = 0.99;
		//if(altFirstPos.Perp() > 2*fRMin)
		//	pull = 1.01;
		//      here.setO(altFirstPos*pull);  //pull it inward a little

		//move first plane to fully include first padhit
		altFirstDir.SetMag(1);
		here.setO(altFirstPos-altFirstDir*0.2);
		here.setNormal(altFirstDir);

		//      if (fVerbose)
		//	{
		//	  std::cout<<"\nto first Digi position (r="<<here.getO().Perp()<<"): \n    " ;
		//	  here.getO().Print();
		//	  std::cout.flush();
		//	  std::cout<<"\n"<<std::endl;
		//	}
	}

	TVector3 stopO;
	// if we want to use the chamber edges they are found and used as sart and stop position
	if (fUseChamberEdge)
	{
		double t=CalcIntersect(here.getO(),here.getNormal());
		TVector3 newDir=here.getNormal();
		newDir.SetMag(1);
		if (fVerbose)
			std::cout<<"the shift value is:"<<t<<std::endl;
		TVector3 newO=here.getO()+t*newDir;
		stopO=here.getO()+CalcIntersect(here.getO(),here.getNormal(),true)*newDir;
		if (fVerbose)
		{
			std::cout<<"The new position is:"<<std::endl;
			newO.Print();
			std::cout<<"The radius is: "<<newO.Perp()<<" The max radius: "<<fRMax<<std::endl;
		}
		startDir=-(newO-here.getO());
		here.setO(newO);
		//      here.setNormal(newDir);
	}

	if (fVerbose)
	{
		std::cout<<"New Plane:"<<here.getO().X()<<":"<<here.getO().Y()<<":"<<here.getO().Z()<<"\n";
		std::cout<<"New Direction:"<<here.getNormal().X()<<":"<<here.getNormal().Y()<<":"<<here.getNormal().Z()<<"\n";
	}

	int clcounter=0;
	bool exc(false);
	//std::set<unsigned int> usedDigis;  //bookkeeping for speedup
	std::map<unsigned int,double> usedDigis; //bookkeeping for speedup
	double dir = 1.;


	//walk through the track and collect the clusters
	bool founddigis=false;
	bool shiftPlane=false;
	double shiftMult=1;
	std::vector<TpcPrelimCluster*> prelimClusters;
	if (fTimeinfo)
		PrintTime(timer,"Got Start position");
	if (fVerbose)
		std::cout<<digis.size()<<" digis in cluster, "<<udigis.size()<<" digis not in cluster"<<std::endl;
	next=here;
	TVector3 tmp;
	while(clcounter<100*hitIDs.size())//stepalong until no digi is found between the planes
	{
		bool fabort;
		bool forceabort(false);
		//      if (fTimeinfo)
		//	PrintTime(timer,"Start a step");
		if (fUseChamberEdge)
			fabort=false;
		else
			fabort=true;


		//      std::map<unsigned int, TpcDigi*> digisToCluster;   //<TCA index, pointer>
		if (fVerbose) std::cerr<<"at counter "<<clcounter<<std::endl;

		if(clcounter++>0 && shiftPlane) here=next;//not the first time

		TVector3 normdir=next.getNormal();
		normdir.SetMag(1);
		if(fVarStep)
		{
			fstepsize=fStepOff + fStepSlope * (next.getO().Z()+fZoff);
			if (fVerbose)
				std::cout<<"The Z position to calculate the stepsize: "<<next.getO().Z()<<"\n";
			if (fstepsize<0.1)
				fstepsize=0.1;
		}
		//----------------------------------------------------------------
		if(fVarStepSqrt)
		{
			double fDiff_l = fGas->Dl();
			double fDiff_t = fGas->Dt();

			TVector3 diffVec=TVector3(normdir);
			diffVec[0]*=fDiff_l;
			diffVec[1]*=fDiff_l;
			diffVec[2]*=fDiff_t;

			//fstepsize = fStepOff_sqrt + fFac* ( TMath::Cos(normdir.Theta())* fDiff_t * sqrt( (next.getO().Z() + here.getO().Z() )*0.5 - fStepOff_0_sqrt);
			//std::cout<<"offset: "<<fStepOff_sqrt<<" factor: "<<fFac<<" diffusion: "<<diffVec.Mag()<<" ("<<diffVec.X()<<" , "<<diffVec.Y()<<" , "<<diffVec.Z()<<")"<<"\n";
			//std::cout<<"diff l: "<<fDiff_l<<" diff t: "<<fDiff_t<<"\n";
			fstepsize = fStepOff_sqrt + fFac * diffVec.Mag() * sqrt( (next.getO().Z() + here.getO().Z() )*0.5 - fStepOff_0_sqrt);
			fstepsize*=gRandom->Uniform(0.9,1.1);
			if (fVerbose)
				std::cout<<"The Z position to calculate the stepsize: "<<next.getO().Z()<<"\n";
			if (fstepsize<0.1)
				fstepsize=0.1;
			if (fstepsize != fstepsize)
				fstepsize=1;
		}
		//----------------------------------------------------------------
		TVector3 destination = next.getO()+normdir*fstepsize*shiftMult*dir; //calculate next point
		//std::cout<<"destination.Mag():"<<destination.Mag()<<" shift plane:"<<shiftPlane<<"\n";
		// next.getO() is the old plane origin
		// normdir is the normalized (Mag=1) direction along the track
		// fstepsize is the stepsize (eithe sqrt dependent, fixed or linear)
		// shiftMult is a multiplicator set to 1 if last cluster was formed normally and set to 0.1 if clustersize was too small during last step
		// dir is there to go along the right direction (actually fixed to one)
		if (fVerbose)
		{
			std::cout<<"Start direction is: "<<dir<<" Stepsize: "<<fstepsize<<std::endl;
			std::cout<<"Destination Point is: ";
			destination.Print();
			std::cout<<"Extrap direction os: ";
			normdir.Print();
		}

		TVector3 poca,dirInPoca;
		try //get the next plane
		{
			theRep->extrapolateToPoint(destination,poca,dirInPoca);
			next.setO(poca);
			next.setNormal(dirInPoca);
			//dist=theRep->extrapolate(next);
		}
		catch(GFException& e)
		{
			e.what();
			exc=true;
			break;
		}

		if (fVerbose)
		{
			std::cout<<"Poca: ";
			poca.Print();
		}

		//check if next lies out of the TPC volume
		//      std::cout<<next.getO().Perp()<<" "<<stopO.Perp()<<std::endl;;
		if(fUseChamberEdge)
		{
			TVector3 td=next.getO()-stopO;
			if (next.getO().Perp()>3*fRMax)
				forceabort=true;
			if (fVerbose)
				std::cout<<"distance to last point is: "<<td.Mag()<<std::endl;
			if(next.getO().Perp() > stopO.Perp() && td.Mag()<fRMax)
			{
				if(fVerbose)
					std::cout<<"OUT OF CHAMBER next.Perp(): "<<next.getO().Perp()<<std::endl;

				//TVector3 shift=stopO-next.getO();
				next.setO(stopO);
				if(fVerbose)
					std::cout<<"AFTER SHIFT  next.Perp(): "<<next.getO().Perp()<<std::endl;
				forceabort = true;
			}
		}
		else
		{
			Bool_t move_plane=kFALSE;
			if (fUseCosmics)
			{
				if(next.getO().Perp() > 1.5*fRMax)
				{
					move_plane=kTRUE;
				}
			}
			else
				if(next.getO().Perp() > 1.5*fRMax || next.getO().Perp() < 0.5*fRMin)
					move_plane=kTRUE;

			if (move_plane)
			{
				if(fVerbose)
					std::cout<<"OUT OF CHAMBER  next.Perp(): "<<next.getO().Perp()<<std::endl;
				TVector3 conn = next.getO() - here.getO();
				conn.SetMag(1.);
				double drn, drh, dr;
				dr = next.getO().Perp() - here.getO().Perp();
				double comp = fRMax;
				if(next.getO().Perp()<fRMin)
					comp = fRMin;
				drn = next.getO().Perp() - comp;
				double moveR = 1.1*fabs(drn)*fstepsize/fabs(dr);   //scaled to get safely inside the cage
				TVector3 newpos = next.getO() - (conn*moveR);
				next.setO(newpos);
				if(fVerbose)
					std::cout<<"AFTER SHIFT  next.Perp(): <<"<<next.getO().Perp()<<"("<<fstepsize/fabs(dr)<<")"<<std::endl;
				forceabort = true;
			}
		}

		if (fVerbose)
		{
			std::cout<<"Checking if digis inside step:\n";
			tmp = here.getO();
			std::cout<<"here center:"<<tmp.X()<<":"<<tmp.Y()<<":"<<tmp.Z()<<std::endl;
			tmp = here.getNormal();
			std::cout<<"here normal:"<<tmp.X()<<":"<<tmp.Y()<<":"<<tmp.Z()<<std::endl;
			tmp = next.getO();
			std::cout<<"next center:"<<tmp.X()<<":"<<tmp.Y()<<":"<<tmp.Z()<<std::endl;
			tmp = next.getNormal();
			std::cout<<"next normal:"<<tmp.X()<<":"<<tmp.Y()<<":"<<tmp.Z()<<std::endl;

		}

		double behindHere,behindNext;

		TpcPrelimCluster* pcl = new TpcPrelimCluster(fPlane,fPrelimClusterCounter,fG,fC);
		pcl->setRecalcDigiPos(frecalcDigiPos);
		pcl->setCorrFactor(fcorrFactor);
		pcl->setDoPRF(fDoPRF);
		std::map<unsigned int, TpcDigi*>::const_iterator cit;
		founddigis=false;

		Int_t rounds;
		if(fUseAllDigis)
			rounds=2;
		else
			rounds=1;

		double meanz;

		for(Int_t round=0;round<rounds;round++)
		{
			std::map<unsigned int, TpcDigi*>dmap;
			if(round==0)
				dmap=digis;
			else
				dmap=udigis;

			for(cit=dmap.begin(); cit!=dmap.end();cit++)
			{
				//try to find the digi in the map
				std::map<unsigned int, double>::iterator it_used=usedDigis.find(cit->first);
				//if found, check if already fully shared
				if(it_used!=usedDigis.end())
				{
					//std::cout<<"digi "<<it_used->first<<" is already used:"<<it_used->second<<"\n";
					//check if digi is fully shared
					if(it_used->second >0.99999999)
					{
						//std::cout<<"digi "<<it_used->first<<" is fully shared\n";
						continue;
					}
				}
				//since digi is not fully/at all used, get position
				TVector3 hpos(0.,0.,0.);
				try
				{
					TpcDigiMapper::getInstance()->map(cit->second,hpos);
				}
				catch(const std::exception& ex)
				{
					std::cout<<ex.what()<<std::endl;
					continue;
				}

				if(fIgnoreEdgeDigis)
					if(hpos.Perp() < fRCutoff1 || hpos.Perp() > fRCutoff2)
						continue;

				//if (fAlign)
				//hpos=TpcAlignmentManager::getInstance()->localToMaster("tpc",hpos);

				//check if digi lies in between planes (deprecated)
				behindHere = here.getNormal() * (here.dist(hpos))*dir;
				behindNext = next.getNormal() * (next.dist(hpos))*dir;
				behindHere *= 1./fabs(behindHere);
				behindNext *= 1./fabs(behindNext);

				double ok = behindHere + behindNext;    // 0 if digi between planes

				//get weight of digi
				TpcPad* pad=fPlane->GetPad(((TpcDigi*)fDigiArray->At(cit->first))->padId());
				double weight1= pad->GetAreaFractionBelowPlane(here.getO()-TVector3(0,0,hpos.Z()),here.getNormal());
				double weight2= pad->GetAreaFractionBelowPlane(next.getO()-TVector3(0,0,hpos.Z()),next.getNormal());

				//ok=weight. 0 if outside both planes,1 if fully inside 0<ok<1 if partially inside
				ok=weight1-weight2;


				if (behindNext < 0.)
					fabort = false;  //there are still viable digis up fronts, keep running

				//if (fabs(behindHere+behindNext)>=0.5 and ok>0)
				//	std::cout<<"digi "<<it_used->first<<" :ok from planepos:"<<fabs(behindHere+behindNext)<<" ok from planecut"<<ok<<"\n";

				//if(fabs(ok) < 0.5)   //possibilities are actually -2, 0 ,2 within double precision...
				if(ok>0)
				{
					if(fVerbose)
					{
						//			  std::cout<<"INSIDE!"<<std::endl;
						std::cout<<"INSIDE! "<<cit->first<<":"<<hpos.X()<<":"<<hpos.Y()<<":"<<hpos.Z()<<std::endl;
						std::cout<<"here: ("<<here.getO().X()<<","<<here.getO().Y()<<","<<here.getO().Z()<<")\n";
						std::cout<<"here: ("<<next.getO().X()<<","<<next.getO().Y()<<","<<next.getO().Z()<<")\n";
						std::cout<<"weight1:"<<weight1<<" weight2:"<<weight2<<"\n";
					}
					TVector3 hd=hpos-here.getO();
					if(fUseAllDigis && round==1)
					{
						if(fVerbose)
							std::cout<<"Checking an unclustered Digi"<<std::endl;
						if (fabs(hd.Perp())>fRDigiCutoff)
							continue;
						else if (fVerbose)
							std::cout<<"Unclustered digi is inside radius ("<<fRDigiCutoff<<")"<<std::endl;
					}
					if(fabs(hd.Perp()>fRMax))
						continue;


					meanz+=hpos.Z();

					if((fabs(here.getO().Y()-hpos.Y())<0.15 or fabs(next.getO().Y()-hpos.Y())<0.15) and fVerbose)
					{
						//std::cout<<"behind here:"<<behindHere<<"\n";
						std::cout<<"here O:"<<here.getO().Y()<<" next O:"<<next.getO().Y()<<"\n";
						std::cout<<"hpos Y:"<<hpos.Y()<<"\n";
						std::cout<<"pad Y:"<<pad->y()<<"\n";
						std::cout<<"weight1:"<<weight1<<" weight2:"<<weight2<<"\n\n";
					}

					if(not fDigiSharing)
						ok=1;

					//update share value of digi
					if(it_used!=usedDigis.end())
					{
						it_used->second+=ok; //add up if already partially used
					}
					else
					{
						usedDigis[cit->first]=ok; //new entry in usedDigis with current share
					}

					pcl->addHit((TpcDigi*) fDigiArray->At(cit->first),here,ok);
					founddigis=true;
				}
				else
				{
					if(fVerbose)
					{
						std::cout<<"NOT INSIDE! "<<cit->first<<":"<<hpos.X()<<":"<<hpos.Y()<<":"<<hpos.Z()<<std::endl;
						std::cout<<"here: ("<<here.getO().X()<<","<<here.getO().Y()<<","<<here.getO().Z()<<")\n";
						std::cout<<"here: ("<<next.getO().X()<<","<<next.getO().Y()<<","<<next.getO().Z()<<")\n";
						std::cout<<"weight1:"<<weight1<<" weight2:"<<weight2<<"\n\n";
					}
				}

				/*
				if (fVerbose)
				{
					std::cout<<" BehindHere: "<<behindHere<<":"<<fabort<<"\n";
					std::cout<<" BehindNext: "<<behindNext<<":"<<forceabort<<"\n";
				}*/
			}
			if (fVerbose)
				std::cout<<"1there are "<<usedDigis.size()<<" digis used and "<<udigis.size()+digis.size()<<" digis in total in round: "<<round<<std::endl;
		}

		//if the clustersize is to small but bigger than 0
		// and if there are still digis up ahead

		meanz/=pcl->getNDigis();
		if(fVarMinClusterSize)
		{
			fMinClusterSize=fMCS1+fMCS2*(meanz+fMCS3);
		}

		if(pcl->getNDigis()<fMinClusterSize && founddigis && !(fabort || forceabort))
		{
			if (fVerbose)
				std::cout<<"Cluster was to small:"<<pcl->getNDigis()<<":"<<fMinClusterSize<<std::endl;
			//check if next plane is outside tpc volume
			if (next.getO().Perp() < fRMin or next.getO().Perp() > fRMax or next.getO().Z()>fZMax)
			{
				if( prelimClusters.size()>0 )
				{
					shiftPlane=true;
					founddigis=false;
					if (fVerbose)
						std::cout<<"Merging Cluster backwards\n";
					for (int ud=0;ud<pcl->getNDigis();ud++)
					{
						prelimClusters.back()->addHit((TpcDigi*)pcl->getDigi(ud),here,pcl->getDigiShare(ud));
					}
				}
			}
			else
			{
				founddigis=false;
				shiftPlane=false;
				if (fVerbose)
					std::cout<<"Just move next plane one step further\n";
				shiftMult=0.1;
				for(int ud=0;ud<pcl->getNDigis();ud++)
				{
					usedDigis.erase(((TpcDigi*)pcl->getDigi(ud))->index());
				}
			}
		}

		//if last cluster in track is to small:
		if( pcl->getNDigis()<fMinClusterSize && founddigis && (fabort || forceabort) )
		{
			if( prelimClusters.size()>0 )
			{
				founddigis=false;
				shiftPlane=true;
				if (fVerbose)
					std::cout<<"Merging Cluster backwards\n";
				for (int ud=0;ud<pcl->getNDigis();ud++)
				{
					prelimClusters.back()->addHit((TpcDigi*)pcl->getDigi(ud),here);
				}
			}
		}

		if (fVerbose)
			std::cout<<"2there are "<<usedDigis.size()<<" digis used and "<<udigis.size()+digis.size()<<" digis in total"<<std::endl;

		if (founddigis)
		{
			fPrelimClusterCounter++;
			if(fVerbose)
				std::cout<<"found some digis, pushing back a cluster"<<std::endl;
			double phi = (next.getO()-here.getO()).Phi();
			pcl->setPhi(phi);
			pcl->setDir(next.getO()-here.getO());

			//try to get the detplane exactly between next and here
			TVector3 Midpoca,dirInMidPoca,Midpoint;
			Midpoint=next.getO()+here.getO();
			Midpoint*=0.5;
			bool extrapGood=false;
			try //get the plane at MidPoca
			{
				theRep->extrapolateToPoint(Midpoint,Midpoca,dirInMidPoca);
				extrapGood=true;
			}
			catch(GFException& e)
			{
				e.what();
				//if extrap didnt work set into middle between planes (ignore bending)
			}

			GFDetPlane MidPlane;
			if (extrapGood)
			{
				pcl->setMidPos(Midpoca);
				MidPlane.setO(Midpoca);
				MidPlane.setNormal(dirInMidPoca);
				pcl->setMidPlane(MidPlane);
			}
			else
			{
				pcl->setMidPos(Midpoint);
				MidPlane.setO(Midpoint);
				MidPlane.setNormal(next.getO()-here.getO());
				pcl->setMidPlane(MidPlane);
			}
			prelimClusters.push_back(pcl);
			shiftPlane=true;
			shiftMult=1.;
		}
		else
		{
			delete pcl;
		}

		if (fabort || forceabort)
			break;

	}//done walk through the track

	if (fTimeinfo)
		PrintTime(timer,"Done stepping along the Track");

	if(fVerbose) std::cout<<"now creating real cluster"<<std::endl;

	unsigned int ncl = prelimClusters.size();
	unsigned int size=fReClusterArray->GetEntries()+ncl;
	fReClusterArray->Expand(size);
	if(fVerbose)
		std::cout<<"There are "<<ncl<<" cluster"<<std::endl;

	if (ncl<2 and fVerbose)
	{
		std::cout<<"Number  of recluster: "<<ncl<<" Number of original: "<<hitIDs.size()<<"\n";
		std::cout<<"The fit flag: "<<theRep->getStatusFlag()<<"\n";
	}

	if (fTrackInit and prelimClusters.size()!=0)
		InitTracks(prelimClusters,trk);

	for(unsigned  int icl=0;icl<ncl;icl++)
	{
		if (!fTrackInit)
		{
			TpcCluster* tcl = prelimClusters[icl]->convTpcCluster(0,fzJitter,true,fuseOnPlaneCov); //be aware, sector set to 0 regardless of real sector position regardless of real value AND saveraw set to true !!!!
			new((*fReClusterArray)[fGlobCounter]) TpcCluster(*tcl);
			fGlobCounter++;
			delete tcl;
		}
		delete prelimClusters[icl];
	}
	prelimClusters.clear();
	//std::cout<<"done"<<std::endl;
	return kTRUE;
}

Bool_t TpcReClusterizerTask::GenfitStepThrough(GFTrack* trk)
{
	FairRootManager * ioman = FairRootManager::Instance();
	//std::cout<<"step only\n";
	if(fVerbose)
		std::cout<<"****************************GENFIT STEPTHROUGH******************"<<std::endl;

	GFAbsTrackRep* theRep = trk->getCardinalRep();
	if (theRep->getStatusFlag() != 0 && fSkipUnFitted)
	{
		if (fVerbose)
			std::cout<<"ReCl: track not fitted!\n";
		return kFALSE;
	}

	GFTrackCand cand = trk->getCand();

	bool useSPHits = true;
	std::vector<int> hitIDs = cand.getHitIDs(fHitBranchID);

	// get pos and momentum for start position and direction
	TVector3 pos,mom;
	if (fVerbose) std::cerr<<"get pos and mom"<<std::endl;
	try
	{
		pos = theRep->getPos();
		mom = theRep->getMom();
	}
	catch(GFException& e)
	{
		e.what();
		return kFALSE;
	}

	TVector3 startDir = mom;
	startDir.SetMag(1.);

	TVector3 altFirstPos;

	//set reference plane:
	GFDetPlane pl;
	pl.setO(pos);
	pl.setNormal(startDir);
	Int_t cosmics=1;
	unsigned int count = 0;
	TClonesArray* actClArr;
	std::vector<TpcPrelimCluster*> prelimClusters;

	for(int i=0; i<hitIDs.size(); i++)
	{
		TpcCluster* cl;
		TpcCluster* lastCl;
		TpcCluster* nextCl;
		TpcSPHit* theHit = (TpcSPHit*) fSPHitArray->At(hitIDs[i]);
		TpcSPHit* lastHit;
		TpcSPHit* nextHit;
		lastHit=nextHit=theHit;
		if (i<hitIDs.size()-1)
			nextHit=(TpcSPHit*) fSPHitArray->At(hitIDs[i+1]);
		if(i>0)
			lastHit=(TpcSPHit*) fSPHitArray->At(hitIDs[i-1]);

		int clBrID = theHit->getClusterBranchID();
		int clID = theHit->getClusterID();
		if (count==0)
		{
			TString clBrName = ioman->GetBranchName(clBrID);
			actClArr = (TClonesArray*) ioman->GetObject(clBrName);
		}

		cl = (TpcCluster*)(actClArr->At(clID));
		lastCl=(TpcCluster*)(actClArr->At(lastHit->getClusterID()));
		//std::cout<<"i:"<<i<<" size:"<<hitIDs.size()<<" next:"<<nextHit<<" this:"<<theHit<<" last:"<<lastHit<<"\n";
		nextCl=(TpcCluster*)(actClArr->At(nextHit->getClusterID()));
		TVector3 theDir= (nextCl->pos()-lastCl->pos())*0.5;
		count++;

		TVector3 poca, dirpoca;
		try
		{
			theRep->extrapolateToPoint(cl->pos(),poca,dirpoca);
		}
		catch(GFException& e)
		{
			//if extrapolation is not working, invent poca and dirpoca by useage of clusters
			e.what();
			std::cerr<<"TpcReclusterizer:could not step to cluster\n";
			poca=cl->pos();
			dirpoca=nextCl->pos()-cl->pos();
			if (dirpoca.Mag()==0)
				dirpoca=cl->pos()-lastCl->pos();
			//return kFALSE;
		}
		GFDetPlane clusterPlane;
		clusterPlane.setO(poca);
		clusterPlane.setNormal(dirpoca);
		TpcPrelimCluster* pcl = new TpcPrelimCluster(fPlane,fPrelimClusterCounter,fG,fC);

		pcl->setRecalcDigiPos(frecalcDigiPos);
		pcl->setCorrFactor(fcorrFactor);
		pcl->setMidPos(poca);
		pcl->setMidPlane(clusterPlane);
		pcl->setDir(theDir);
		pcl->setPhi(theDir.Phi());
		const std::map<unsigned int,double>* digiMap = cl->getDigiMap();
		for(std::map<unsigned int, double>::const_iterator it=digiMap->begin();it!=digiMap->end();it++)
		{
			pcl->addHit((TpcDigi*) fDigiArray->At(it->first),clusterPlane,it->second);
			//std::cout<<"digi weight:"<<it->second<<"\n";
		}
		fPrelimClusterCounter++;
		prelimClusters.push_back(pcl);
	}

	if (fTrackInit)
		InitTracks(prelimClusters,trk);

	for(unsigned  int icl=0;icl<prelimClusters.size();icl++)
	{
		delete prelimClusters[icl];
	}
	prelimClusters.clear();
	return kTRUE;
}

TpcCluster* TpcReClusterizerTask::MergeCluster(std::vector<TpcCluster*> ftomerge)
{
	TpcCluster* merged;
	for (int i=0; i<ftomerge.size();i++)
	{
		//      for
	}
	return merged;
}

double TpcReClusterizerTask::CalcIntersect(TVector3 p, TVector3 dir, bool neg)
{
	// set radius to inner radius for NOT cosmics!
	dir.SetMag(1);
	double a = dir.Perp2();
	double b = 2*(dir.X()*p.X()+dir.Y()*p.Y());
	double c = p.Perp2()-(fRMax*fRMax) ;

	double disc = b*b - 4*a*c;
	if (disc>0)
	{
		double t;
		if (!neg)
			t = (-b + TMath::Sqrt(disc)) / (2*a);
		else
			t = (-b - TMath::Sqrt(disc)) / (2*a);
		return t;
	}
	return 0;
}

void TpcReClusterizerTask::PrintTime(TStopwatch watch, TString text)
{
	watch.Stop();
	Double_t rtime = watch.RealTime()*1000;
	Double_t ctime = watch.CpuTime()*1000;
	watch.Start(kFALSE);
	std::cout<<"RealTime (ms): "<<rtime<<" CpuTime (ms): "<<ctime<<";;;"<<text<<"\n";
}

void TpcReClusterizerTask::InitTracks(std::vector<TpcPrelimCluster*> _prelimClusters,GFTrack* _trk)
{
	fRecoHitMap.clear();
	fToWrite.clear();
	// build GFTrackCands
	std::vector<GFTrackCand*> candlist;
	// create GFTrackCands
	GFTrackCand* cand=new GFTrackCand();
	// fill hits into GFTrackCands  and get seed values
	TpcCluster* cluster1;//=new TpcCluster();
	for (int ih = 0; ih<_prelimClusters.size();ih++)
	{
		//std::cout<<"converting cluster number: "<<ih<<"\n";
		TpcCluster* cl = _prelimClusters[ih]->convTpcCluster(0,fzJitter,true,fuseOnPlaneCov); //be aware, sector set to 0 regardless of real sector position, regardless of real value AND saveRaw set to true !!!!
		//std::cout<<"converted cluster\n";
		if(fSmallClusterCorr)
			cl->SetPos(fcorrector->ExecCorrection(cl,true,true));
		if (fClusterCorr)
			ClusterCorr(cl);

		if (ih==0)
		{
			cluster1=new TpcCluster(*cl);
		}
		new((*fReClusterArray)[fGlobCounter]) TpcCluster(*cl);
		TpcSPHit* theHit = new TpcSPHit(cl,fweightedPlaneConstruction,fcutCov);
		theHit->setUseParam(fUseSPHitParam);
		theHit->setClusterBranchID(fClBrID);
		unsigned int nOut = fRecoHitArray->GetEntriesFast();
		new ((*fRecoHitArray)[nOut]) TpcSPHit(*theHit);
		fToWrite.insert(theHit);
		cand->addHit(fSPHitID, nOut);
		fGlobCounter++;
		delete cl;
	}

	// if (fTrAlign)
	//   {
	//     TpcAlignmentManager* alMan=TpcAlignmentManager::getInstance();
	//     pos1=alMan->localToMaster("tpc",pos1);
	//     mom=alMan->localToMasterVect("tpc",mom);

	//     poserr=alMan->localToMasterVect("tpc",poserr);
	//     momerr=alMan->localToMasterVect("tpc",momerr);
	//   }

	cand->setMcTrackId((_trk->getCand()).getMcTrackId());
	candlist.push_back(cand);
	//RK TRACKREP

	// store GFTracks in output array
	TVector3 mom,pos;
	TVector3 momerr;
	//cluster1->sig(0).Print();
	//std::cout<<"err X:"<<cluster1->sig(0).X()<<" err Y:"<<cluster1->sig(0).Y()<<" err Z:"<<cluster1->sig(0).Z()<<"\n";


	//std::cout<<"getting err\n";
	TVector3 poserr=TVector3(cluster1->sig(0));
	TMatrixDSym cov;
	_trk->getPosMomCov(pos,mom,cov);
	//mom=cand->getMomSeed();
	//cov.ResizeTo((_trk->getCand()).getCovSeed());
	//cov=(_trk->getCand()).getCovSeed();
	if(0)
	{
		std::cout<<"TpcReclusterizerTask::InitTracks momentum:\n";
		mom.Print();
		std::cout<<"TpcReclusterizerTask::InitTracks cov:\n";
		cov.Print();
	}
	//poserr=TVector3(TMath::Sqrt(cov(0,0)),TMath::Sqrt(cov(1,1)),TMath::Sqrt(cov(2,2)));
	momerr=TVector3(TMath::Sqrt(cov(3,3)),TMath::Sqrt(cov(4,4)),TMath::Sqrt(cov(5,5)));
	poserr*=0.00001;
	momerr*=0.00001;


	RKTrackRep* rkrep = new RKTrackRep(cluster1->pos(), mom, poserr, momerr,((RKTrackRep*)_trk->getCardinalRep())->getPDG());

	//TMatrixDSym trackcov=rkrep->getCov();
	//trackcov*=25;
	//rkrep->setCov(trackcov);

	//GFTrack* gftrk=new((*fOutTrackArray)[fOutTrackArray->GetEntriesFast()]) GFTrack(*_trk);
	GFTrack* gftrk=new((*fOutTrackArray)[fOutTrackArray->GetEntriesFast()]) GFTrack(rkrep);
	gftrk->setSmoothing(true);
	gftrk->setCandidate(*cand); // here the candidate is copied!
	//((RKTrackRep*)gftrk->getCardinalRep())->setStatusFlag(1);
	//gftrk->blowUpCovs(10000);
	//write out all other recoHits if required, clean up:
	std::map<TpcCluster*, TpcSPHit*>::const_iterator it;
	for(it=fRecoHitMap.begin(); it!=fRecoHitMap.end(); it++)
	{
		TpcSPHit* iHit = it->second;
		if(!fToWrite.count(iHit))
		{
			if(fArhpersistence && fPersistence)
				new ((*fRecoHitArray)[fRecoHitArray->GetEntriesFast()]) TpcSPHit(*iHit);
		}
		delete iHit;  //all of them have been copied during TCA write ...
	}
	delete cluster1;

}

void TpcReClusterizerTask::ClusterCorr(TpcCluster* cl)
{
	TVector3 pos=cl->pos();
	TVector3 corrval=fDevMap->value(pos);
	//corrval.SetZ(0); //dont correct in z position since drift velocity is not properly set in distortion map creation
	cl->SetPos(pos-corrval);
	cl->SetFieldCorr(corrval);
	//std::cout<<"zpos: "<<pos.Z()<<" Z-corrval:"<<corrval.Z()<<"\n";
}