// -------------------------------------------------------------------------
// -----                CbmSttMinuitTrackFitter source file            -----
// -----                  Created 29/03/06  by R. Castelijns           -----
// -------------------------------------------------------------------------
#include <string>
#include <sstream>

#include "CbmSttMinuitTrackFitter.h"
#include "CbmRootManager.h"
#include "CbmTask.h"
#include "CbmSttGeomPoint.h"

#include "TArc.h"
#include "TF1.h"
#include "TH2.h"
#include "TH3.h"
#include "TLine.h"
#include "TMarker.h"
#include "TGraphErrors.h"
using namespace std;    

CbmSttMinuitTrackFitter::CbmSttMinuitTrackFitter()
{
    rootoutput = kFALSE;
    fEventCounter = 0;
    fVerbose = 0;
}

CbmSttMinuitTrackFitter::CbmSttMinuitTrackFitter(Int_t verbose)
{
    fVerbose = verbose;

    if (verbose < 3)
	rootoutput = kFALSE;
    else
	rootoutput = kTRUE;

    fEventCounter = 0;
}

CbmSttMinuitTrackFitter::~CbmSttMinuitTrackFitter()
{
}

void CbmSttMinuitTrackFitter::Init()
{
    fEventCounter = 0;

  // Get and check CbmRootManager
  CbmRootManager
    *ioman = CbmRootManager::Instance();

  if (! ioman) 
    {
      cout << "-E- CbmSttMinuitTrackFitter::Init: "
	   << "RootManager not instantised!" << endl;
      //      return kFATAL;
    }
}

Double_t CbmSttMinuitTrackFitter::GetHitRadius(Int_t hitNo)
{
    CbmSttHit
	*pMhit = NULL;
    
    pMhit = GetHitFromCollections(hitNo);
	
    return sqrt(pMhit->GetX() * pMhit->GetX() + pMhit->GetY() * pMhit->GetY());
}

Double_t CbmSttMinuitTrackFitter::GetHitAngle(Int_t hitNo, Double_t dCenter, 
					      Double_t phiCenter, Double_t radius)
{
    CbmSttHit
	*pMhit = NULL;

    pMhit = GetHitFromCollections(hitNo);
	
    Double_t
	xCenter = (dCenter + radius) * cos(phiCenter),
	yCenter = (dCenter + radius) * sin(phiCenter);
	
    // get XY of wire
    Double_t
	deltaX = pMhit->GetX() - xCenter,
	deltaY = pMhit->GetY() - yCenter;
	
    Double_t
	angle = atan(deltaY / deltaX);
	
    // bring angle into the usual frame of reference
    if (deltaX < 0.)
    {
	if (deltaY < 0.)
	    angle -= dPi;
	else
	    angle += dPi;
    }

    return angle;
}

void CbmSttMinuitTrackFitter::OrderHitsByR(map<Double_t, Int_t> &hitMap)
{
    // sort the hits by r position
    CbmSttHit
	*pMhit = NULL;
    
    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	pMhit = GetHitFromCollections(iHit);
	if ( ! pMhit ) 
	    continue;
	
	hitMap[sqrt(pMhit->GetX() * pMhit->GetX() + pMhit->GetY() * pMhit->GetY())] = iHit;
    }
}

Bool_t CbmSttMinuitTrackFitter::GetCharge(Double_t dCenter, Double_t phiCenter, Double_t radius, Int_t &charge, Int_t type)
{
    map<Double_t, Int_t>
	hitMap;
    
    if (type > 0)
	OrderHitsByR(hitMap);
    else
	OrderHitsByZ(hitMap);
    
    map<Double_t, Int_t>::iterator
	hitMapStart = hitMap.begin(),
	hitMapEnd = hitMap.end();

    hitMapEnd--;

    CbmSttHit
	*startHit = GetHitFromCollections(hitMapStart->second),
	*endHit = GetHitFromCollections(hitMapEnd->second);

    CbmSttGeomPoint
	vertex(0., 0., 0.),
	firstPoint(startHit->GetX(), startHit->GetY(), startHit->GetZ()), 
	lastPoint(endHit->GetX(), endHit->GetY(), endHit->GetZ());

    Bool_t
	retval = kTRUE;

    Double_t
	angle,
	oldAngle,
	difAngle;

    Int_t 
	votesForPositive = 0,
	votesForNegative = 0,
	hitNo = 0;

    map<Double_t, Int_t>::iterator
	hitMapIter = hitMap.begin();

    while (hitMapIter != hitMap.end())
    {
	// get hit's angle
	angle = GetHitAngle(hitMapIter->second, dCenter, phiCenter, radius);
	
	// store the first hit's angle as a reference (hit with smallest z coordinate)
	if (hitNo == 0)
	{
	    refAngle = angle;
	}

	// calculate the angle difference between the reference hit and this hit
	difAngle = angle - refAngle;

	while (difAngle < 0)
	{
	    difAngle += 2 * dPi;
	}

	// cast votes for criterium 1
	if (hitNo > 0)
	{
	    if (difAngle > oldAngle)
		votesForPositive++;
	    else
		votesForNegative++;
	}

	oldAngle = difAngle;

	hitNo++;
	hitMapIter++;
    }

    // 2/3 majority is enough to decide on track
    if (votesForPositive > 2 * votesForNegative)
    {
	retval = kTRUE;

	// when swimming the track from -z to z we get a positive rotation
	if (firstPoint.DistanceTo(vertex) < lastPoint.DistanceTo(vertex))
	{
	    charge = -1;
	}
	else
	{
	    charge = 1;
	}
    }
    else if (votesForNegative > 2 * votesForPositive)
    {
	retval = kTRUE;

	// when swimming the track from -z to z we get a positive rotation
	if (firstPoint.DistanceTo(vertex) < lastPoint.DistanceTo(vertex))
	{
	    charge = 1;
	}
	else
	{
	    charge = -1;
	}
    }
    else
    {
	charge = 0;
	retval = kFALSE;
    }

    return retval;
}

void CbmSttMinuitTrackFitter::OrderHitsByZ(map<Double_t, Int_t> &hitMap)
{
    // sort the hits by z position
    CbmSttHit
	*pMhit = NULL;

    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	pMhit = GetHitFromCollections(iHit);
	if ( ! pMhit ) 
	    continue;

	hitMap[pMhit->GetZ()] = iHit;
    }
}

Double_t CbmSttMinuitTrackFitter::CalculateRadialHitCoordinatesCircular(Double_t dCenter, Double_t phiCenter, Double_t radius)
{
    Int_t
	charge = 0;


    if (GetCharge(dCenter, phiCenter, radius, charge, 0))
    {
	map<Double_t, Int_t>
	    hitMap;

	OrderHitsByZ(hitMap);

	map<Double_t, Int_t>::iterator
	    hitMapIter = hitMap.begin();

	Double_t
	    angle,
	    difAngle,
	    oldAngle,
	    accumulatedAngle = 0.;

	while (hitMapIter != hitMap.end())
	{
	    angle = GetHitAngle(hitMapIter->second, dCenter, phiCenter, radius);

	    // calculate the angle difference between this hit and the last on 0,2pi
	    if (hitMapIter == hitMap.begin())
		difAngle = 0.;
	    else
		difAngle = angle - oldAngle;

	    while (difAngle < -1 * dPi)
		difAngle += 2 * dPi;
	    while (difAngle >= dPi)
		difAngle -= 2 * dPi;

	    accumulatedAngle += difAngle;
	    oldAngle = angle;
	    
	    CbmSttHit
		*pMhit = GetHitFromCollections(hitMapIter->second);

	    pMhit->SetRadial(radius * accumulatedAngle);

	    hitMapIter++;
	}
    }
    else
    {
	if (GetCharge(dCenter, phiCenter, radius, charge, 1))
	{
	    cout << "falling back to method 2" << endl;
	}
	else
	{
	    cout << "cannot determine charge" << endl;
	}
    }
    return charge;
}

void CbmSttMinuitTrackFitter::plotAllIsochrones(CbmSttTrack* pTrack)
{
    for (Int_t hitCounter = 0; hitCounter < pTrack->GetNofHits(); hitCounter++)
    {
	Int_t
	    hitNumber = pTrack->GetHitIndex(hitCounter);

	CbmSttHit
	    *pHit = GetHitFromCollections(hitNumber);
	
	if (!pHit)
	    continue;

	if (rootoutput)
	{
	    TArc  
		myArc;
	    
	    Double_t
		radius = pHit->GetIsochrone();
	    
	    if (radius == 0)
		radius = 1.0;
	    
	    myArc.SetFillStyle(0);
	    myArc.SetLineColor(4);
	    myArc.DrawArc(pHit->GetX(), pHit->GetY(), radius);
	}
    }
}

Int_t CbmSttMinuitTrackFitter::DoFit(CbmSttTrack* pTrack, Int_t pidHypo)
{    
  cout << "hello " << pTrack->GetNofHits() << endl;
    fFirstHit = kTRUE;
    fRefAngle = 0.;

    fEventCounter++;
    fTrack = pTrack;
    
    if (rootoutput)
    {
	eventCanvas = new TCanvas("eventcanvas", "eventcanvas", 800, 600); 
	eventCanvas->Divide(2, 2);
	cout << "-W- CbmSttMinuitTrackFitter::Init: dividing canvas!"
	     << endl;

	eventCanvas->cd(1);
	TH2F myRange("range", "range", 100, -50., 50., 100, -50., 50);
	myRange.DrawCopy();
	plotAllIsochrones(pTrack);
	plotAllStraws();

	eventCanvas->Update();
    }
    
    Double_t
	rstart = pTrack->GetParamLast()->GetTx(),
	dstart = pTrack->GetParamLast()->GetX(),
	phistart = pTrack->GetParamLast()->GetY(),
	dstart2 = pTrack->GetParamLast()->GetZ(),
	thetastart = pTrack->GetParamLast()->GetTy();

    cout << "[[[[[[[start[[[[[[[[[[[[[[[[[ " << rstart << endl;

    TMinuit 
      minimizer(3);
    
    minimizer.SetFCN(fcnRadialCircular);
    minimizer.SetErrorDef(1);
    
    minimizer.DefineParameter(0, "d0", dstart, 0.1, 0., 0.);
    minimizer.DefineParameter(1, "phi0", phistart, 0.1, 0., 0.);
    minimizer.DefineParameter(2, "r", rstart, 0.1, 0., 0.);
    
    minimizer.SetObjectFit(this);
    
    minimizer.SetPrintLevel(-1);
    minimizer.SetMaxIterations(500);

    minimizer.Migrad();
    
    Double_t
	chisquare,
	resultsRadial[3], 
	errorsRadial[3]; 
    
    minimizer.GetParameter(0, resultsRadial[0], errorsRadial[0]);
    minimizer.GetParameter(1, resultsRadial[1], errorsRadial[1]);
    minimizer.GetParameter(2, resultsRadial[2], errorsRadial[2]);
    minimizer.Eval(3, NULL, chisquare, resultsRadial, 0); 

    cout << "[[[[[[[end[[[[[[[[[[[[[[[[[ " << resultsRadial[2] << endl;

    if (rootoutput)
    {
	eventCanvas->cd(1);
	
	double
	  xstart = (dstart + rstart) * cos(phistart),
	  ystart = (dstart + rstart) * sin(phistart),
	  xend = (resultsRadial[0] + resultsRadial[2]) * cos(resultsRadial[1]),
	  yend = (resultsRadial[0] + resultsRadial[2]) * sin(resultsRadial[1]);

	TArc
	    myArc;
	
	myArc.SetFillStyle(0);
	myArc.SetLineColor(2);
	myArc.SetLineWidth(1);
	myArc.DrawArc(xend, yend, resultsRadial[2]);

	myArc.SetLineColor(3);
	myArc.SetLineWidth(1);
	myArc.DrawArc(xstart, ystart, rstart);

	cout << "fitted result: dpr  " << resultsRadial[0] << "  " << resultsRadial[1] << "  " << resultsRadial[2] << endl;
	eventCanvas->Show();
	eventCanvas->Update();
    }
  
    // recalculate hit positions along radial vector using y-axis crossing as origin
    Double_t
	charge = CalculateRadialHitCoordinatesCircular(resultsRadial[0], resultsRadial[1], resultsRadial[2]);

    GetLinearSeed(dstart2, thetastart);

    if (rootoutput)
    {
	eventCanvas->cd(3);
	plotAllIsochronesZR();

	Double_t
	    z1 = dstart2 * cos(thetastart) - 100. * sin(thetastart),
	    r1 = dstart2 * sin(thetastart) + 100. * cos(thetastart),
	    z2 = dstart2 * cos(thetastart) + 100. * sin(thetastart),
	    r2 = dstart2 * sin(thetastart) - 100. * cos(thetastart);

	TLine
	    myLine;

	myLine.SetLineColor(2);
	myLine.DrawLine(z1, r1, z2, r2);
    }

    TMinuit 
      minimizer2(2);
    
    minimizer2.SetFCN(fcnLongitudinal);
    minimizer2.SetErrorDef(1);
    
    minimizer2.DefineParameter(0, "d1", dstart2, 0.1, 0., 0.);
    minimizer2.DefineParameter(1, "theta", thetastart, 0.1, 0., 0.);

    minimizer2.SetObjectFit(this);
    
    minimizer2.SetPrintLevel(-1);
    minimizer2.SetMaxIterations(500);

    minimizer2.Migrad();
    
    Double_t
	chisquareLongitudinal = 0.,
	resultsLongitudinal[2], 
	errorsLongitudinal[2]; 
    
    minimizer2.GetParameter(0, resultsLongitudinal[0], errorsLongitudinal[0]);
    minimizer2.GetParameter(1, resultsLongitudinal[1], errorsLongitudinal[1]);
    minimizer2.Eval(2, NULL, chisquareLongitudinal, resultsLongitudinal, 0); 

    Double_t
	z1 = resultsLongitudinal[0] * cos(resultsLongitudinal[1]) - 10000. * sin(resultsLongitudinal[1]),
	r1 = resultsLongitudinal[0] * sin(resultsLongitudinal[1]) + 10000. * cos(resultsLongitudinal[1]),
	z2 = resultsLongitudinal[0] * cos(resultsLongitudinal[1]) + 10000. * sin(resultsLongitudinal[1]),
	r2 = resultsLongitudinal[0] * sin(resultsLongitudinal[1]) - 10000. * cos(resultsLongitudinal[1]);

    if (rootoutput)
    {
	TLine
	    myLine;

	myLine.SetLineColor(3);
	myLine.DrawLine(z1, r1, z2, r2);
	
	char
	    goOnChar;
	
	eventCanvas->Update();
	eventCanvas->Show();
	
	stringstream
	    conv;
	
	string
	    eventCounterStr;
	
	conv << fEventCounter;
	conv >> eventCounterStr;
	
	string
	    fileName = "fitterView" + eventCounterStr + ".C";
	
	eventCanvas->SaveSource(fileName.c_str()); 
	
	cout << "press any key to continue: " << endl;
	cin >> goOnChar;
	
	eventCanvas->Clear();
	delete eventCanvas;
    }

    Double_t
	zslope = (r2 - r1) / (z2 - z1),
	zoffset = r1 - zslope * z1;

    Double_t
	alpha = resultsRadial[2] / zslope,
	z0 = (-1 * zoffset / zslope) - alpha * refAngle;
        
    // fill the CbmTrackInfo object and do some checks on fit convergence
    pTrack->SetChi2Long(chisquareLongitudinal);
    pTrack->SetChi2Rad(chisquare);
    pTrack->SetNDF(fTrack->GetNofHits());
  
    pTrack->GetParamLast()->SetX(resultsRadial[0]);                
    pTrack->GetParamLast()->SetY(resultsRadial[1]);
    pTrack->GetParamLast()->SetZ(z0);                
    pTrack->GetParamLast()->SetTx(resultsRadial[2]);
    pTrack->GetParamLast()->SetTy(alpha);                
    pTrack->GetParamLast()->SetQp(charge);
    
    return 0;
}

Bool_t CbmSttMinuitTrackFitter::putStraw(Double_t xpos, Double_t ypos, Double_t radius)
{
    Double_t
      pipeDiam = 0.42,
      tubeOuterDiam = 1.032,
      CoverThickness = 0.1,
      outerRadius = 42.,
      innerRadius = 15.;

    if ((sqrt(xpos * xpos + ypos * ypos) < outerRadius - CoverThickness - (tubeOuterDiam / 2.)) &&
	(sqrt(xpos * xpos + ypos * ypos) > innerRadius + CoverThickness + (tubeOuterDiam / 2.)) &&
	(sqrt(ypos * ypos) > (pipeDiam + tubeOuterDiam / 2.)))
    {
	eventCanvas->cd(1);
	
	TArc  
	    myArc;
	
	myArc.SetFillStyle(0);
	myArc.SetLineColor(5);
	myArc.DrawArc(xpos, ypos, radius);
    }
    else
    {
	return false;
    }
    
    return true;
}


void CbmSttMinuitTrackFitter::plotAllStraws()
{
    Double_t
	tubeOuterDiam = 1.01,
	sttCenterX = 0.,
	sttCenterY = 0.;
    
    Int_t
	ringmax = -1;
    
    Bool_t
	started = kFALSE,
	goOn = kTRUE;
    
    Int_t
	ringteller = 18;
    
    while (goOn) 
    {
	goOn = kFALSE;
	
	Double_t
	    sqrt3 = sqrt(3.),
	    radius = tubeOuterDiam / 2.;
	
	Double_t
	    zpos = radius;
	
	// place first
	Double_t
	    xpos = sttCenterX + ((ringteller) * 2 * radius),  
	    ypos = sttCenterY;
	
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos -= radius;
	    ypos += sqrt3 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;
	}
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos -= 2 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;
	}
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos -= radius;
	    ypos -= sqrt3 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;	   
	}
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos += radius;
	    ypos -= sqrt3 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;
	}
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos += 2 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;
	}
	for(Int_t i = 0; i < ringteller; i++) 
	{
	    xpos += radius;
	    ypos += sqrt3 * radius;
	    if (putStraw(xpos, ypos, zpos))
		goOn = kTRUE;
	}
	
	if (goOn)
	    started = kTRUE;
	
	if (!started)
	    goOn = kTRUE;
	
	ringteller++;
	if ((ringmax > 0) && (ringteller == ringmax))
	    goOn = kFALSE;
    }
}

void fcnRadialCircular(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  Int_t
    elementCounter = 0;

  const CbmSttMinuitTrackFitter 
    *mama = (CbmSttMinuitTrackFitter *)gMinuit->GetObjectFit();

  CbmSttTrack
    *fTrack = mama->GetTrack();

  Double_t 
    chisq = 0.,
    measuredRadius = 0.,
    normalizedDistance = 0.,
    radialResolution = 0.,
    distanceOfClosestApproach = 0.,
    x_0 = 0.,
    y_0 = 0.,
    xCircleCenter = 0.,
    yCircleCenter = 0.,
    rCircle = 0.; 

  CbmSttHit
    *pMhit = NULL;
  
  for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
      // get index of hit
      Int_t
	iHit = fTrack->GetHitIndex(i);
      
      // get hit
      pMhit = mama->GetHitFromCollections(iHit);

      if ( ! pMhit ) 
	continue;

      elementCounter++;

      // get isochrone
      measuredRadius = pMhit->GetIsochrone();
      radialResolution = pMhit->GetIsochroneError();

      // get XY of wire
      x_0 = pMhit->GetX();
      y_0 = pMhit->GetY();
 
      xCircleCenter = (par[0] + par[2]) * cos(par[1]);
      yCircleCenter = (par[0] + par[2]) * sin(par[1]);
      rCircle = par[2];

      distanceOfClosestApproach = sqrt((x_0 - xCircleCenter) * (x_0 - xCircleCenter) + 
				       (y_0 - yCircleCenter) * (y_0 - yCircleCenter)) - rCircle;

      distanceOfClosestApproach = sqrt(distanceOfClosestApproach * distanceOfClosestApproach);
      
      if (radialResolution > 0.)
	{
	    normalizedDistance = ((distanceOfClosestApproach - measuredRadius) * (distanceOfClosestApproach - measuredRadius)) / 
		(radialResolution * radialResolution);
	}
      else
	{
	    normalizedDistance = (distanceOfClosestApproach * distanceOfClosestApproach) / 
	      (pMhit->GetDx() * pMhit->GetDx() + pMhit->GetDy() * pMhit->GetDy());
	}
      chisq += normalizedDistance;
    }

  chisq /= elementCounter;
  f = chisq;
}

void fcnLongitudinal(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
    const CbmSttMinuitTrackFitter 
	*mama = (CbmSttMinuitTrackFitter *)gMinuit->GetObjectFit();
    
    CbmSttTrack
	*fTrack = mama->GetTrack();
    
    Int_t
	elementCounter = 0;
    
    Double_t 
	chisq = 0.,
	normalizedDistance = 0.,
	distanceOfClosestApproach = 0.,
	z_0, 
	r_0,
	z_1,
	r_1,
	z_2,
	r_2; 
    
    CbmSttHit
	*pMhit = NULL;
    
    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	pMhit = mama->GetHitFromCollections(iHit);
	
	if (!pMhit)
	    continue;
	
	elementCounter++;
	
	// get ZR of wire
	z_0 = pMhit->GetZ();
	r_0 = pMhit->GetRadial();

	z_1 = par[0] * cos(par[1]) - 100. * sin(par[1]);
	r_1 = par[0] * sin(par[1]) + 100. * cos(par[1]);
	
	z_2 = par[0] * cos(par[1]) + 100. * sin(par[1]); 
	r_2 = par[0] * sin(par[1]) - 100. * cos(par[1]); 

	Double_t
	    slope = ((r_2 - r_1) / (z_2 - z_1)),
	    offset = r_1 - slope * z_1;

	distanceOfClosestApproach = z_0 - ((r_0 - offset) / slope);
	distanceOfClosestApproach = sqrt(distanceOfClosestApproach * distanceOfClosestApproach);

	if (pMhit->GetDz() > 0.)
	{
	    normalizedDistance = (distanceOfClosestApproach * distanceOfClosestApproach) / 
		(pMhit->GetDz() * pMhit->GetDz());
	}
	else
	{
	    cout << "-W- CbmSttMinuitTrackFitter::fcnLongitudinal: Position resolution is 0, not possible!"
		 << endl;
	}

	chisq += normalizedDistance;
    }
    
    chisq /= elementCounter;
    f = chisq;
}

void CbmSttMinuitTrackFitter::Extrapolate(CbmSttTrack* track, Double_t r, CbmTrackParam *param )
{
  cout << "-W- CbmSttMinuitTrackFitter::Extrapolate: Not yet implemented, sorry!"
       << endl;
}

CbmSttHit* CbmSttMinuitTrackFitter::GetHitFromCollections(Int_t hitCounter) const
{
    CbmSttHit
	*retval = NULL;
 
    Int_t
	relativeCounter = hitCounter;

    for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++)
    {
	Int_t
	    size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast();
	if (relativeCounter < size)
	{
	    retval = (CbmSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter);
	    break;
	}
	else
	{
	    relativeCounter -= size;
	}
    }
    return retval;
}

void CbmSttMinuitTrackFitter::GetLinearSeed(Double_t &dist, Double_t &theta)
{
    Double_t
	maxD = 0.,
	maxR = 0.,
	minR = 0.,
	maxZ = 0.,
	minZ = 0.;
	
    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	CbmSttHit
	    *pMhit = GetHitFromCollections(iHit);

	if (!pMhit)
	    continue;
	
	for (Int_t ii = 0; ii < fTrack->GetNofHits(); ii++)
	{
	    // get index of hit
	    Int_t
		iiHit = fTrack->GetHitIndex(ii);
	
	    // get hit
	    CbmSttHit
		*pMhit2 = GetHitFromCollections(iiHit);

	    if (!pMhit2)
		continue;
	    
	    Double_t
		dist = sqrt((pMhit->GetZ() - pMhit2->GetZ()) * (pMhit->GetZ() - pMhit2->GetZ()) +
			    (pMhit->GetRadial() - pMhit2->GetRadial()) * (pMhit->GetRadial() - pMhit2->GetRadial()));

	    if (dist > maxD)
	    {
		maxD = dist;
		maxZ = pMhit->GetZ();
		maxR = pMhit->GetRadial();
		minZ = pMhit2->GetZ();
		minR = pMhit2->GetRadial();
	    }
	}
    }

    // get normal form
    theta = atan((maxZ - minZ) / (minR - maxR));

    if ((minR - maxR) < 0.)
    {
	if ((maxZ - minZ) < 0.)
	    theta -= dPi;
	else
	    theta += dPi;
    }

    dist = (minR - (((maxR - minR) / (maxZ  - minZ)) * minZ)) * sin(theta);
}

void CbmSttMinuitTrackFitter::plotAllIsochronesZR()
{
    eventCanvas->cd(3);

    Double_t
	maxR = -100000.,
	minR = 100000.,
	maxZ = -100000.,
	minZ = 100000.;

    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	CbmSttHit
	    *pMhit = GetHitFromCollections(iHit);

	if (!pMhit)
	    continue;

	for (Int_t ii = 0; ii < fTrack->GetNofHits(); ii++)
	{
	    // get index of hit
	    Int_t
		iiHit = fTrack->GetHitIndex(ii);
	
	    // get hit
	    CbmSttHit
		*pMhit2 = GetHitFromCollections(iiHit);

	    if (!pMhit2)
		continue;
	    
	    if (maxZ < pMhit->GetZ())
	    {
		maxZ = pMhit->GetZ();
	    }
	    if (maxR < pMhit->GetRadial())
	    {
		maxR = pMhit->GetRadial();
	    }
	    if (minZ > pMhit->GetZ())
	    {
		minZ = pMhit->GetZ();
	    }
	    if (minR > pMhit->GetRadial())
	    {
		minR = pMhit->GetRadial();
	    }
	}
    }
    
    TH2F myRange2("range2", "range2", 100, minZ - 10., maxZ + 10., 100, minR - 10., maxR + 10.);
    myRange2.DrawCopy();
    
    for (Int_t i = 0; i < fTrack->GetNofHits(); i++)
    {
	// get index of hit
	Int_t
	    iHit = fTrack->GetHitIndex(i);
	
	// get hit
	CbmSttHit
	    *pMhit = GetHitFromCollections(iHit);

	cout << "z coordinate of hit: " << pMhit->GetZ() << endl;
	
	if (!pMhit)
	    continue;
	
	TLine
	    myLine;

	myLine.DrawLine(pMhit->GetZ() - pMhit->GetDz(), pMhit->GetRadial(), pMhit->GetZ() + pMhit->GetDz(), pMhit->GetRadial());
	myLine.DrawLine(pMhit->GetZ(), pMhit->GetRadial() - 0.1, pMhit->GetZ(), pMhit->GetRadial() + 0.1);
    }
}

ClassImp(CbmSttTrackFitter)