// -------------------------------------------------------------------------
// -----                CbmSttTrackFinderHough source file             -----
// -----                  Created 28/03/06  by R. Castelijns           -----
// -------------------------------------------------------------------------


// C++ includes
#include "iostream.h"
#include <map>
#include <string>
#include <sstream>
#include <algorithm>
using namespace std;

// ROOT includes
#include "TClonesArray.h"
#include "TArc.h"
#include "TLine.h"
#include "TH2.h"
#include "TH3.h"
#include "TStopwatch.h"
#include "TGraph.h"
#include "TF1.h"

// CBM includes
#include "CbmMCPoint.h"
#include "CbmMCTrack.h"
#include "CbmRootManager.h"
#include "CbmSttHit.h"
#include "CbmSttTrack.h"
#include "CbmSttTrackFinderHough.h"

#define CbmSttHoughThetaBins 100
#define CbmSttHoughPhiBins 500

#define CbmSttHoughThreshold 1
#define CbmSttHoughCDistBins 150
#define CbmSttHoughCPhiBins 150
#define CbmSttHoughCRadBins 50
#define CbmHoughChiSquareCut 10
#define CbmSttHoughParallelCut 0.00001

#define	limit 0.0000000001
#define doCircular 1

#define rootoutput kFALSE
#define saverootoutput kFALSE
#define DO_TIMING kTRUE
#define PLOT_ALL_STRAWS kTRUE

static Int_t trackThis;

// -----   Default constructor   -------------------------------------------
CbmSttTrackFinderHough::CbmSttTrackFinderHough() 
{ 
  fVerbose      = 1;
  fEventCounter = 0;
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmSttTrackFinderHough::CbmSttTrackFinderHough(Int_t verbose) 
{ 
  fVerbose      = verbose;
  fEventCounter = 0;
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmSttTrackFinderHough::~CbmSttTrackFinderHough() 
{ 
}
// -------------------------------------------------------------------------



// -----   Public method Init   --------------------------------------------
void CbmSttTrackFinderHough::Init() 
{
  fEventCounter = 0;

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

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

  // create the accumulator
  if (doCircular)
  {
      fAccumulator = new CbmSttHoughAccumulatorNew(CbmSttHoughThreshold, kTRUE,
						   CbmSttHoughCDistBins, -75., 75., 
						   CbmSttHoughCPhiBins, -2 * dPi, 2 * dPi,
						   CbmSttHoughCRadBins, -400, 400
	  );
  }
  else
  {
      fAccumulator = new CbmSttHoughAccumulatorNew(CbmSttHoughThreshold, kTRUE,
						   CbmSttHoughThetaBins, -50., 50., 
						   CbmSttHoughPhiBins, -1 * dPi, dPi
	  );
  }
}


Bool_t CbmSttTrackFinderHough::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.)))
    {
	TArc  
	    myArc;
	
	myArc.SetFillStyle(0);
	myArc.SetLineColor(17);
	myArc.DrawArc(xpos, ypos, radius);
    }
    else
    {
	return false;
    }
    
    return true;
}


void CbmSttTrackFinderHough::plotAllStraws()
{
    Double_t
	tubeOuterDiam = 1.032,
	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 CbmSttTrackFinderHough::GetTrackletCircular(Double_t firstX, Double_t firstY, Double_t firstR, 
						 Double_t secondX, Double_t secondY, Double_t secondR, 
						 Double_t thirdX, Double_t thirdY, Double_t thirdR, 
						 Double_t *circleRadii, Double_t *circleCentersX, 
						 Double_t *circleCentersY) const
{
    Int_t
	trackletCounter = 0;

    for (Int_t sign1 = -1; sign1 < 3; sign1 += 2)
    {
	for (Int_t sign2 = -1; sign2 < 3; sign2 += 2)
	{
	    for (Int_t sign3 = -1; sign3 < 3; sign3 += 2)
	    {
   
		Double_t
		    a = -2. * (firstX - secondX),
		    b = -2. * (firstY - secondY),
		    c = 2. * (sign1 * firstR - sign2 * secondR),
		    d = ((firstX * firstX + firstY * firstY - firstR * firstR) - 
			 (secondX * secondX + secondY * secondY - secondR * secondR)),
		    e = -2. * (firstX - thirdX),
		    f = -2. * (firstY - thirdY),
		    g = 2. * (sign1 * firstR - sign3 * thirdR),
		    h = ((firstX * firstX + firstY * firstY - firstR * firstR) - 
			 (thirdX * thirdX + thirdY * thirdY - thirdR * thirdR)),
		    
		    A = -1. * (f * d - b * h) / (a * f - b * e),
		    B = -1. * (f * c - b * g) / (a * f - b * e),
		    C = (e * d - a * h) / (a * f - e * b),
		    D = (e * c - a * g) / (a * f - e * b), 
		    
		    I = B * B + D * D - 1.,
		    II = 2. * A * B - 2. * B * firstX + 2. * C * D - 2. * D * firstY + sign1 * 2. * firstR,
		    III = A * A - 2. * A * firstX + firstX * firstX + C * C - 2. * C * firstY + firstY * firstY - firstR * firstR,
		    
		    r = (-1. * II - sqrt(II * II - 4. * I * III)) / (2. * I),
		    
		    x = A + B * r,
		    y = C + D * r;  

		circleRadii[trackletCounter] = sqrt(r * r);
		circleCentersX[trackletCounter] = x;
		circleCentersY[trackletCounter] = y;
		trackletCounter++;
	    }
	}
    }  
}




void CbmSttTrackFinderHough::GetTracklet(Double_t thisX, Double_t thisY, Double_t thisR, 
					 Double_t otherX, Double_t otherY, Double_t otherR, 
					 Double_t *contactPointsXWithOther, Double_t *contactPointsYWithOther, 
					 Double_t *angleWithOther, 
					 Double_t *contactPointsXAtOther, Double_t *contactPointsYAtOther) const
{
    // TODO: split this up in four sub functions:
    //    calculate tangent angles 
    //    calculate angles
    //    calculate contact points
    //    sort results
    
    // citation: Eric W. Weisstein. "Circle Tangent Line." 
    // From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/CircleTangentLine.html
    
    // calculate the external similitude center
    Double_t
	exX = (thisR * otherX - otherR * thisX) / (thisR - otherR),
	exY = (thisR * otherY - otherR * thisY) / (thisR - otherR);

    // displace the geometry such that external center will be at origin  
    Double_t
	otherXex = otherX - exX,
	otherYex = otherY - exY,
	thisXex = thisX - exX,
	thisYex = thisY - exY;
    
    // calculate the four possible tangent angles from the external similitude points with this circle
    Double_t
	tEx0 = +1 * acos((-thisR * thisXex + thisYex * sqrt(thisXex * thisXex + thisYex * thisYex - thisR * thisR)) / 
			 (thisXex * thisXex + thisYex * thisYex)),
	tEx1 = +1 * acos((-thisR * thisXex - thisYex * sqrt(thisXex * thisXex + thisYex * thisYex - thisR * thisR)) / 
			 (thisXex * thisXex + thisYex * thisYex)),
	tEx2 = -1 * acos((-thisR * thisXex + thisYex * sqrt(thisXex * thisXex + thisYex * thisYex - thisR * thisR)) / 
			 (thisXex * thisXex + thisYex * thisYex)),
	tEx3 = -1 * acos((-thisR * thisXex - thisYex * sqrt(thisXex * thisXex + thisYex * thisYex - thisR * thisR)) / 
			 (thisXex * thisXex + thisYex * thisYex));
    
    // calculate four possible contact points with this circle
    Double_t
	xEx0 = thisX + thisR * cos(tEx0),
	xEx1 = thisX + thisR * cos(tEx1),
	xEx2 = thisX + thisR * cos(tEx2),
	xEx3 = thisX + thisR * cos(tEx3),
	yEx0 = thisY + thisR * sin(tEx0),
	yEx1 = thisY + thisR * sin(tEx1),
	yEx2 = thisY + thisR * sin(tEx2),
	yEx3 = thisY + thisR * sin(tEx3);

    // calculate the four possible tangent angles from the external similitude points with other circle
    Double_t
	tEx0other = +1 * acos((-otherR * otherXex + otherYex * sqrt(otherXex * otherXex + otherYex * otherYex - otherR * otherR)) / 
			      (otherXex * otherXex + otherYex * otherYex)),
	tEx1other = +1 * acos((-otherR * otherXex - otherYex * sqrt(otherXex * otherXex + otherYex * otherYex - otherR * otherR)) / 
			      (otherXex * otherXex + otherYex * otherYex)),
	tEx2other = -1 * acos((-otherR * otherXex + otherYex * sqrt(otherXex * otherXex + otherYex * otherYex - otherR * otherR)) / 
			      (otherXex * otherXex + otherYex * otherYex)),
	tEx3other = -1 * acos((-otherR * otherXex - otherYex * sqrt(otherXex * otherXex + otherYex * otherYex - otherR * otherR)) / 
			      (otherXex * otherXex + otherYex * otherYex));
    
    // calculate four possible contact points with other circle
    Double_t
	xEx0other = otherX + otherR * cos(tEx0other),
	xEx1other = otherX + otherR * cos(tEx1other),
	xEx2other = otherX + otherR * cos(tEx2other),
	xEx3other = otherX + otherR * cos(tEx3other),
	yEx0other = otherY + otherR * sin(tEx0other),
	yEx1other = otherY + otherR * sin(tEx1other),
	yEx2other = otherY + otherR * sin(tEx2other),
	yEx3other = otherY + otherR * sin(tEx3other);

    // calculate line coefficients and offsets
    Double_t
	coefEx0 = (yEx0 - exY) / (xEx0 - exX),
	coefEx1 = (yEx1 - exY) / (xEx1 - exX),
	coefEx2 = (yEx2 - exY) / (xEx2 - exX),
	coefEx3 = (yEx3 - exY) / (xEx3 - exX),
	offEx0 = exY - coefEx0 * exX,
	offEx1 = exY - coefEx1 * exX,
	offEx2 = exY - coefEx2 * exX,
	offEx3 = exY - coefEx3 * exX;

    // and another set
    Double_t
	coefEx0other = (yEx0other - exY) / (xEx0other - exX),
	coefEx1other = (yEx1other - exY) / (xEx1other - exX),
	coefEx2other = (yEx2other - exY) / (xEx2other - exX),
	coefEx3other = (yEx3other - exY) / (xEx3other - exX),
	offEx0other = exY - coefEx0other * exX,
	offEx1other = exY - coefEx1other * exX,
	offEx2other = exY - coefEx2other * exX,
	offEx3other = exY - coefEx3other * exX;

    Int_t
	thisteller = 0;
    
    if (NumberOfLineCircleIntersections(exX, exY, xEx0, yEx0, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xEx0;
	    contactPointsYWithOther[thisteller] = yEx0;
	    angleWithOther[thisteller] = tan((yEx0 - exY) / (xEx0 - exX));
	    
	    if ((NumberOfLineCircleIntersections(exX, exY, xEx0other, yEx0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefEx0 - coefEx0other) < limit) && (fabs(offEx0 - offEx0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx0other;
		contactPointsYAtOther[thisteller] = yEx0other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx1other, yEx1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx0 - coefEx1other) < limit) && (fabs(offEx0 - offEx1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx1other;
		contactPointsYAtOther[thisteller] = yEx1other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx2other, yEx2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx0 - coefEx2other) < limit) && (fabs(offEx0 - offEx2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx2other;
		contactPointsYAtOther[thisteller] = yEx2other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx3other, yEx3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx0 - coefEx3other) < limit) && (fabs(offEx0 - offEx3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx3other;
		contactPointsYAtOther[thisteller] = yEx3other;
	    }
	    
	    thisteller++;
	}
    }
    

    if (NumberOfLineCircleIntersections(exX, exY, xEx1, yEx1, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xEx1;
	    contactPointsYWithOther[thisteller] = yEx1;
	    angleWithOther[thisteller] = tan((yEx1 - exY) / (xEx1 - exX));
	    if ((NumberOfLineCircleIntersections(exX, exY, xEx0other, yEx0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefEx1 - coefEx0other) < limit) && (fabs(offEx1 - offEx0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx0other;
		contactPointsYAtOther[thisteller] = yEx0other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx1other, yEx1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx1 - coefEx1other) < limit) && (fabs(offEx1 - offEx1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx1other;
		contactPointsYAtOther[thisteller] = yEx1other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx2other, yEx2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx1 - coefEx2other) < limit) && (fabs(offEx1 - offEx2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx2other;
		contactPointsYAtOther[thisteller] = yEx2other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx3other, yEx3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx1 - coefEx3other) < limit) && (fabs(offEx1 - offEx3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx3other;
		contactPointsYAtOther[thisteller] = yEx3other;
	    }
	    
	    thisteller++;
	}
    }
    
    
    if (NumberOfLineCircleIntersections(exX, exY, xEx2, yEx2, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xEx2;
	    contactPointsYWithOther[thisteller] = yEx2;
	    angleWithOther[thisteller] = tan((yEx2 - exY) / (xEx2 - exX));
	    
	    if ((NumberOfLineCircleIntersections(exX, exY, xEx0other, yEx0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefEx2 - coefEx0other) < limit) && (fabs(offEx2 - offEx0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx0other;
		contactPointsYAtOther[thisteller] = yEx0other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx1other, yEx1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx2 - coefEx1other) < limit) && (fabs(offEx2 - offEx1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx1other;
		contactPointsYAtOther[thisteller] = yEx1other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx2other, yEx2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx2 - coefEx2other) < limit) && (fabs(offEx2 - offEx2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx2other;
	      contactPointsYAtOther[thisteller] = yEx2other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx3other, yEx3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx2 - coefEx3other) < limit) && (fabs(offEx2 - offEx3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx3other;
		contactPointsYAtOther[thisteller] = yEx3other;
	    }
	    
	    thisteller++;
	}
    }
    
    if (NumberOfLineCircleIntersections(exX, exY, xEx3, yEx3, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xEx3;
	    contactPointsYWithOther[thisteller] = yEx3;
	    angleWithOther[thisteller] = tan((yEx3 - exY) / (xEx3 - exX));
	    
	    if ((NumberOfLineCircleIntersections(exX, exY, xEx0other, yEx0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefEx3 - coefEx0other) < limit) && (fabs(offEx3 - offEx0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx0other;
		contactPointsYAtOther[thisteller] = yEx0other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx1other, yEx1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx3 - coefEx1other) < limit) && (fabs(offEx3 - offEx1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx1other;
		contactPointsYAtOther[thisteller] = yEx1other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx2other, yEx2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx3 - coefEx2other) < limit) && (fabs(offEx3 - offEx2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx2other;
		contactPointsYAtOther[thisteller] = yEx2other;
	    }
	    else if ((NumberOfLineCircleIntersections(exX, exY, xEx3other, yEx3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefEx3 - coefEx3other) < limit) && (fabs(offEx3 - offEx3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xEx3other;
		contactPointsYAtOther[thisteller] = yEx3other;
	    }
	    
	    thisteller++;
	}
    }
    
    
    
    // calculate the internal similitude center
    Double_t
	inX = (thisR * otherX + otherR * thisX) / (thisR + otherR),
	inY = (thisR * otherY + otherR * thisY) / (thisR + otherR);
    
    // displace the geometry such that internal center will be at origin  
    Double_t
	otherXin = otherX - inX,
	otherYin = otherY - inY,
	thisXin = thisX - inX,
	thisYin = thisY - inY;
    
    // calculate the four possible tangent angles from the external similitude points
    Double_t
	tIn0 = +1 * acos((-thisR * thisXin + thisYin * sqrt(thisXin * thisXin + thisYin * thisYin - thisR * thisR)) / 
			 (thisXin * thisXin + thisYin * thisYin)),
	tIn1 = +1 * acos((-thisR * thisXin - thisYin * sqrt(thisXin * thisXin + thisYin * thisYin - thisR * thisR)) / 
			 (thisXin * thisXin + thisYin * thisYin)),
	tIn2 = -1 * acos((-thisR * thisXin + thisYin * sqrt(thisXin * thisXin + thisYin * thisYin - thisR * thisR)) / 
			 (thisXin * thisXin + thisYin * thisYin)),
	tIn3 = -1 * acos((-thisR * thisXin - thisYin * sqrt(thisXin * thisXin + thisYin * thisYin - thisR * thisR)) / 
			 (thisXin * thisXin + thisYin * thisYin));
    
    // calculate four possible contact points with this circle
    Double_t
	xIn0 = thisX + thisR * cos(tIn0),
	xIn1 = thisX + thisR * cos(tIn1),
	xIn2 = thisX + thisR * cos(tIn2),
	xIn3 = thisX + thisR * cos(tIn3),
	yIn0 = thisY + thisR * sin(tIn0),
	yIn1 = thisY + thisR * sin(tIn1),
	yIn2 = thisY + thisR * sin(tIn2),
	yIn3 = thisY + thisR * sin(tIn3);
    
    
    // calculate the four possible tangent angles from the internal similitude points with other circle
    Double_t
	tIn0other = +1 * acos((-otherR * otherXin + otherYin * sqrt(otherXin * otherXin + otherYin * otherYin - otherR * otherR)) / 
			      (otherXin * otherXin + otherYin * otherYin)),
	tIn1other = +1 * acos((-otherR * otherXin - otherYin * sqrt(otherXin * otherXin + otherYin * otherYin - otherR * otherR)) / 
			      (otherXin * otherXin + otherYin * otherYin)),
	tIn2other = -1 * acos((-otherR * otherXin + otherYin * sqrt(otherXin * otherXin + otherYin * otherYin - otherR * otherR)) / 
			      (otherXin * otherXin + otherYin * otherYin)),
	tIn3other = -1 * acos((-otherR * otherXin - otherYin * sqrt(otherXin * otherXin + otherYin * otherYin - otherR * otherR)) / 
			      (otherXin * otherXin + otherYin * otherYin));
    
    // calculate four possible contact points with other circle
    Double_t
	xIn0other = otherX + otherR * cos(tIn0other),
	xIn1other = otherX + otherR * cos(tIn1other),
	xIn2other = otherX + otherR * cos(tIn2other),
	xIn3other = otherX + otherR * cos(tIn3other),
	
	yIn0other = otherY + otherR * sin(tIn0other),
	yIn1other = otherY + otherR * sin(tIn1other),
	yIn2other = otherY + otherR * sin(tIn2other),
	yIn3other = otherY + otherR * sin(tIn3other);
    
    // calculate line coefficients and offsets
    Double_t
	coefIn0 = (yIn0 - inY) / (xIn0 - inX),
	coefIn1 = (yIn1 - inY) / (xIn1 - inX),
	coefIn2 = (yIn2 - inY) / (xIn2 - inX),
	coefIn3 = (yIn3 - inY) / (xIn3 - inX),
	offIn0 = inY - coefIn0 * inX,
	offIn1 = inY - coefIn1 * inX,
	offIn2 = inY - coefIn2 * inX,
	offIn3 = inY - coefIn3 * inX;
    
    // and another set
    Double_t
	coefIn0other = (yIn0other - inY) / (xIn0other - inX),
	coefIn1other = (yIn1other - inY) / (xIn1other - inX),
	coefIn2other = (yIn2other - inY) / (xIn2other - inX),
	coefIn3other = (yIn3other - inY) / (xIn3other - inX),
	offIn0other = inY - coefIn0other * inX,
	offIn1other = inY - coefIn1other * inX,
	offIn2other = inY - coefIn2other * inX,
	offIn3other = inY - coefIn3other * inX;
    
    if (NumberOfLineCircleIntersections(inX, inY, xIn0, yIn0, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xIn0;
	    contactPointsYWithOther[thisteller] = yIn0;
	    angleWithOther[thisteller] = tan((yIn0 - inY) / (xIn0 - inX));
	    
	    if ((NumberOfLineCircleIntersections(inX, inY, xIn0other, yIn0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefIn0 - coefIn0other) < limit) && (fabs(offIn0 - offIn0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn0other;
		contactPointsYAtOther[thisteller] = yIn0other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn1other, yIn1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn0 - coefIn1other) < limit) && (fabs(offIn0 - offIn1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn1other;
		contactPointsYAtOther[thisteller] = yIn1other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn2other, yIn2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn0 - coefIn2other) < limit) && (fabs(offIn0 - offIn2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn2other;
		contactPointsYAtOther[thisteller] = yIn2other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn3other, yIn3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn0 - coefIn3other) < limit) && (fabs(offIn0 - offIn3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn3other;
		contactPointsYAtOther[thisteller] = yIn3other;
	    }
	    
	    thisteller++;
	}
    }
    
    
    if (NumberOfLineCircleIntersections(inX, inY, xIn1, yIn1, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xIn1;
	    contactPointsYWithOther[thisteller] = yIn1;
	    angleWithOther[thisteller] = tan((yIn1 - inY) / (xIn1 - inX));
	    
	    if ((NumberOfLineCircleIntersections(inX, inY, xIn0other, yIn0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefIn1 - coefIn0other) < limit) && (fabs(offIn1 - offIn0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn0other;
		contactPointsYAtOther[thisteller] = yIn0other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn1other, yIn1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn1 - coefIn1other) < limit) && (fabs(offIn1 - offIn1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn1other;
		contactPointsYAtOther[thisteller] = yIn1other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn2other, yIn2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn1 - coefIn2other) < limit) && (fabs(offIn1 - offIn2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn2other;
		contactPointsYAtOther[thisteller] = yIn2other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn3other, yIn3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn1 - coefIn3other) < limit) && (fabs(offIn1 - offIn3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn3other;
		contactPointsYAtOther[thisteller] = yIn3other;
	    }
	    
	    thisteller++;
	}
    }
    
    
    if (NumberOfLineCircleIntersections(inX, inY, xIn2, yIn2, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xIn2;
	    contactPointsYWithOther[thisteller] = yIn2;
	    angleWithOther[thisteller] = tan((yIn2 - inY) / (xIn2 - inX));
	    if ((NumberOfLineCircleIntersections(inX, inY, xIn0other, yIn0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefIn2 - coefIn0other) < limit) && (fabs(offIn2 - offIn0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn0other;
		contactPointsYAtOther[thisteller] = yIn0other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn1other, yIn1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn2 - coefIn1other) < limit) && (fabs(offIn2 - offIn1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn1other;
		contactPointsYAtOther[thisteller] = yIn1other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn2other, yIn2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn2 - coefIn2other) < limit) && (fabs(offIn2 - offIn2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn2other;
		contactPointsYAtOther[thisteller] = yIn2other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn3other, yIn3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn2 - coefIn3other) < limit) && (fabs(offIn2 - offIn3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn3other;
		contactPointsYAtOther[thisteller] = yIn3other;
	    }

	    thisteller++;
	}
    }
    
    
    if (NumberOfLineCircleIntersections(inX, inY, xIn3, yIn3, thisX, thisY, thisR) == 1)
    {
	if (thisteller < 4)
	{
	    contactPointsXWithOther[thisteller] = xIn3;
	    contactPointsYWithOther[thisteller] = yIn3;
	    if ((NumberOfLineCircleIntersections(inX, inY, xIn0other, yIn0other, otherX, otherY, otherR) == 1) &&
		(fabs(coefIn3 - coefIn0other) < limit) && (fabs(offIn3 - offIn0other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn0other;
		contactPointsYAtOther[thisteller] = yIn0other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn1other, yIn1other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn3 - coefIn1other) < limit) && (fabs(offIn3 - offIn1other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn1other;
		contactPointsYAtOther[thisteller] = yIn1other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn2other, yIn2other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn3 - coefIn2other) < limit) && (fabs(offIn3 - offIn2other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn2other;
		contactPointsYAtOther[thisteller] = yIn2other;
	    }
	    else if ((NumberOfLineCircleIntersections(inX, inY, xIn3other, yIn3other, otherX, otherY, otherR) == 1) &&
		     (fabs(coefIn3 - coefIn3other) < limit) && (fabs(offIn3 - offIn3other) < limit))
	    {
		contactPointsXAtOther[thisteller] = xIn3other;
		contactPointsYAtOther[thisteller] = yIn3other;
	    }
	    
	    thisteller++;
	}
    }
}


// -----   Public method DoFind   ------------------------------------------
Int_t CbmSttTrackFinderHough::DoFind(TClonesArray* trackArray) 
{
    Int_t
	retval = 0;

    if (doCircular)
    {
	retval = DoFindCircular(trackArray);
    }
/*
    else
    {
	retval = DoFindLinear(trackArray);
    }
*/    
    return retval;
}

Double_t CbmSttTrackFinderHough::CalculateRadialHitCoordinatesCircular(Double_t xCenter, Double_t yCenter, Double_t radius, CbmSttHit *pMhit)
{
    Double_t 
	angle,
	x0,
	y0,
	wireX = 0.,
	wireY = 0.,
	distance = 0.,
	deltaX = 0., 
	deltaY = 0.,
	isochroneRadius = 0.;
    
    // get XY of wire
    wireX = pMhit->GetX();
    wireY = pMhit->GetY();
    isochroneRadius = pMhit->GetIsochrone();

    deltaX = wireX - xCenter;
    deltaY = wireY - yCenter;

    distance = sqrt(deltaX * deltaX + deltaY * deltaY);

    if (distance > radius)
    {
	x0 = ((deltaX * (distance - isochroneRadius)) / distance) + xCenter;
	y0 = ((deltaY * (distance - isochroneRadius)) / distance) + yCenter;
    }
    else
    {
	x0 = ((deltaX * (distance + isochroneRadius)) / distance) + xCenter;
	y0 = ((deltaY * (distance + isochroneRadius)) / distance) + yCenter;
    }
      
    angle = atan((y0 - yCenter) / (x0 - xCenter));
    return radius * angle;
}

Int_t CbmSttTrackFinderHough::DoFindCircular(TClonesArray* trackArray) 

{
//    cout << "in circular finder" << endl;
    list<CbmSttTrack> 
	tmpTrackArray;

    Double_t
	totaltime = 0.,
	findtime = 0.,
	filltime = 0.,
	createtime = 0.;

    fEventCounter++;
    
    if (rootoutput)
    {
      finderCanvas = new TCanvas("findercanvas", "findercanvas", 800, 600); 
      TH2F myRange("range", "range", 100, -50., 50., 100, -50., 50);
      myRange.DrawCopy();
      if (PLOT_ALL_STRAWS)
	  plotAllStraws();
    }
 
    TH3F
	*hough;
    
    if (rootoutput)
    {
	hough = new TH3F("hough", "hough", CbmSttHoughCDistBins, -75., 75., 
			 CbmSttHoughCPhiBins, -1 * dPi, dPi,
			 CbmSttHoughCRadBins, -200, 200
	    );
    }
    
    TStopwatch 
	timer;
    
    if (DO_TIMING)
    {
	timer.Start(kTRUE);
    }

    fAccumulator->ClearAccumulator();
    fAccumulator->SetThreshold(CbmSttHoughThreshold);

    if (DO_TIMING)
    {
	timer.Stop();
	totaltime += timer.CpuTime();
	createtime = timer.CpuTime();
    }

    if (fHitCollectionList.GetEntries() == 0)
    {
	cout << "-E- CbmSttTrackFinderHough::DoFind: "
	     << "No hit arrays present, call AddHitCollection() first (at least once)! " << endl;
	return -1;
    }
    
    if ( !trackArray ) 
    {
	cout << "-E- CbmSttTrackFinderHough::DoFind: "
	     << "Track array missing! " << endl;
	return -1;
    }
    
    // Create pointers to hit and SttPoint
    CbmSttHit*       pHitFirst  = NULL;
    CbmSttHit*       pHitSecond = NULL;
    CbmSttHit*       pHitThird  = NULL;
    CbmSttTrack*  pTrack     = NULL;

    Int_t 
	nHits = 0;

    for (Int_t hitListCounter = 0; hitListCounter < fHitCollectionList.GetEntries(); hitListCounter++)
    {
	nHits += ((TClonesArray *)fHitCollectionList.At(hitListCounter))->GetEntriesFast();
    }

    Int_t
	ptIndexFirst = 0,
	ptIndexSecond = 0,
	ptIndexThird = 0;

    // guard against ridiculous situations
    if (nHits > 150)
      return 0;

    TVector3
	parallelVector(0., 0., 1.);

    if (DO_TIMING)
    {
	timer.Start(kTRUE);
    }

    // fill accumulator
    for (Int_t iHitFirst = 0; iHitFirst < nHits - 2; iHitFirst++)
    {
	pHitFirst = GetHitFromCollections(iHitFirst);
	
	if ( ! pHitFirst ) 
	    continue;

	if (pHitFirst->GetWireDirection().Angle(parallelVector) > CbmSttHoughParallelCut)
	{
	    continue;
	}

	ptIndexFirst = pHitFirst->GetRefIndex();
	
	if (ptIndexFirst < 0)
	    continue;           // fake or background hit
	
	Double_t
	    firstX = pHitFirst->GetX(),
	    firstY = pHitFirst->GetY(),
	    firstR = pHitFirst->GetIsochrone();

	for (Int_t iHitSecond = iHitFirst + 1; iHitSecond < nHits - 1; iHitSecond++)
	{
	    pHitSecond = GetHitFromCollections(iHitSecond);
	    	    
	    if ( ! pHitSecond ) 
		continue;

	    if (pHitSecond->GetWireDirection().Angle(parallelVector) > CbmSttHoughParallelCut) 
	    {
		continue;
	    }
	
	    ptIndexSecond = pHitSecond->GetRefIndex();
	    
	    if (ptIndexSecond < 0)
		continue;           // fake or background hit
	    
	    Double_t
		secondX = pHitSecond->GetX(),
		secondR = pHitSecond->GetIsochrone(),
		secondY = pHitSecond->GetY();

	    for (Int_t iHitThird = iHitSecond + 1; iHitThird < nHits; iHitThird++)
	    {
		pHitThird = GetHitFromCollections(iHitThird);
		
		if ( ! pHitThird ) 
		    continue;
		
		
		if (pHitThird->GetWireDirection().Angle(parallelVector) > CbmSttHoughParallelCut)
		{
		    continue;
		}

		ptIndexThird = pHitThird->GetRefIndex();
		
		if (ptIndexThird < 0)
		    continue;           // fake or background hit
		
		Double_t
		    thirdX = pHitThird->GetX(),
		    thirdR = pHitThird->GetIsochrone(),
		    thirdY = pHitThird->GetY();

		// all hits regardless of track
		if (rootoutput)
		{
		    TArc  
			myArc;
		    myArc.SetFillStyle(0);
		    myArc.SetLineColor(1);

		    if (firstR == 0)
			myArc.DrawArc(firstX, firstY, 1.0);
		    else
			myArc.DrawArc(firstX, firstY, firstR);

		    if (secondR == 0)
			myArc.DrawArc(secondX, secondY, 1.0);
		    else
			myArc.DrawArc(secondX, secondY, secondR);

		    if (thirdR == 0)
			myArc.DrawArc(thirdX, thirdY, 1.0);
		    else
			myArc.DrawArc(thirdX, thirdY, thirdR);
		}
	
		Double_t
		    deltaR1 = sqrt((thirdY - secondY) * (thirdY - secondY) + (thirdX - secondX) * (thirdX - secondX)),
		    deltaR2 = sqrt((secondY - firstY) * (secondY - firstY) + (secondX - firstX) * (secondX - firstX));

		// take only hits which are close together
		if ((fabs(deltaR1) < 10.)  && (fabs(deltaR2) < 10.))
		{
		    Double_t
			circleRadii[8],
			circleCentersX[8],
			circleCentersY[8];
		    
		    GetTrackletCircular(firstX, firstY, firstR, 
					secondX, secondY, secondR, 
					thirdX, thirdY, thirdR, 
					circleRadii, circleCentersX, circleCentersY);

		    // loop over all tracklets between this and next
		    for (int trackletteller = 0; trackletteller < 8; trackletteller++)
		    {
			Double_t
			    phi = atan(circleCentersY[trackletteller] / circleCentersX[trackletteller]),
			    dist = sqrt(circleCentersX[trackletteller] * circleCentersX[trackletteller] +
					circleCentersY[trackletteller] * circleCentersY[trackletteller])
			    - circleRadii[trackletteller];

			if (circleCentersX[trackletteller] < 0.)
			{
			    if (circleCentersY[trackletteller] < 0.)
				phi -= 3.1415;
			    else
				phi += 3.1415;
			}

			Double_t
			    radial1 = CalculateRadialHitCoordinatesCircular(circleCentersX[trackletteller], 
									    circleCentersY[trackletteller], 
									    circleRadii[trackletteller], 
									    pHitFirst),
			    radial2 = CalculateRadialHitCoordinatesCircular(circleCentersX[trackletteller], 
									    circleCentersY[trackletteller], 
									    circleRadii[trackletteller], 
									    pHitSecond),
			    radial3 = CalculateRadialHitCoordinatesCircular(circleCentersX[trackletteller], 
									    circleCentersY[trackletteller], 
									    circleRadii[trackletteller], 
									    pHitThird);

			Double_t
			    deltaZ1 = (radial1 - radial2) / (pHitFirst->GetZ() - pHitSecond->GetZ()),
			    deltaZ2 = (radial2 - radial3) / (pHitSecond->GetZ() - pHitThird->GetZ());

			if (fabs(deltaZ1 - deltaZ2) < 0.5)
			{
			    fAccumulator->AddHit(dist, iHitFirst, 
						 phi, iHitSecond, 
						 circleRadii[trackletteller], iHitThird);
			}
		    }
		}
	    }
	}
    }

    if (DO_TIMING)
    {
	timer.Stop();
	totaltime += timer.CpuTime();
	filltime = timer.CpuTime();
    }

    Int_t
	nTmpTracks = 0;

    Int_t
	maximumContents;

    if (DO_TIMING)
    {
	timer.Start(kTRUE);
    }
    
    // cluster accumulator
    fAccumulator->Cluster();

    // get all peaks, sort by height
    maximumContents = fAccumulator->GetHighestPeakContents();

    Double_t
	distSeed, phiSeed, rSeed;
    
    distSeed = fAccumulator->GetHighestPeakX();
    phiSeed = fAccumulator->GetHighestPeakY();
    rSeed = fAccumulator->GetHighestPeakZ();
    
    while (maximumContents > 5)
    {
	vector<Int_t>
	    peakContents;
	
	fAccumulator->GetHighestPeakEntryList(peakContents);
	cout << "difsize2: " << fAccumulator->GetHighestPeak()->DifSize() << endl;

	Double_t
	    totalChiSquare = 0.;

	Int_t 
	    totalHits = 0;

	for (Int_t iHitCounter = 0; iHitCounter < peakContents.size(); iHitCounter++)
	{
	    Int_t
		iHit = peakContents[iHitCounter];
	    
	    if (iHit != -1)
	    {
		Double_t
		    myX = (fAccumulator->GetHighestPeakX() + fAccumulator->GetHighestPeakZ()) * cos(fAccumulator->GetHighestPeakY()),
		    myY = (fAccumulator->GetHighestPeakX() + fAccumulator->GetHighestPeakZ()) * sin(fAccumulator->GetHighestPeakY()),
		    myR = fAccumulator->GetHighestPeakZ();

		CbmSttHit
		    *pHit = GetHitFromCollections(iHit);

		totalChiSquare += GetChiSquare(myX, myY, myR, pHit);
		totalHits++;
	    }
	}

	Double_t
	    dSeed = fAccumulator->GetHighestPeakX(),
	    pSeed = fAccumulator->GetHighestPeakY(),
	    rSeed = fAccumulator->GetHighestPeakZ();

	Int_t
	    dSeedBin = fAccumulator->GetHighestPeakXBin(),
	    pSeedBin = fAccumulator->GetHighestPeakYBin(),
	    rSeedBin = fAccumulator->GetHighestPeakZBin();
	
	if ((totalChiSquare / totalHits) < 100.)
	{
	    cout << "chisquare okay" << endl;
	    vector<Int_t>
		possibleHits;

	    for (Int_t iHit = 0; iHit < peakContents.size(); iHit++)
	    {
		Int_t
		    iHitFirst = peakContents[iHit];

		if (iHitFirst != -1)
		{
		    possibleHits.push_back(peakContents[iHit]);
		}
	    }

	    Bool_t
		someOriginalHitsLeft = kFALSE,
		hitsChanged = kTRUE;

	    Int_t
		iterationCounter;

	    while (hitsChanged)
	    {
		ZoomTrack(dSeed, pSeed, rSeed, possibleHits);
	    
		Int_t
		    firstNumber = possibleHits.size();

		AddMissingHits(dSeed, pSeed, rSeed, possibleHits);

		Int_t
		    lastNumber = possibleHits.size();

		someOriginalHitsLeft = kFALSE;

		for (Int_t iOriginalHit = 0; iOriginalHit < peakContents.size(); iOriginalHit++)
		{
		    for (Int_t iHit = 0; iHit < possibleHits.size(); iHit++)
		    {
			if (possibleHits[iHit] == peakContents[iOriginalHit])
			{
			    cout << "found original hit " << possibleHits[iHit] << endl;
			    someOriginalHitsLeft = kTRUE;
			    break;
			}
		    }

		    if (someOriginalHitsLeft)
			break;
		}

		cout << "hit change: " << firstNumber << "  " << lastNumber << endl;
		if (!someOriginalHitsLeft)
		{
		    cout << "no originals left" << endl;
		}

		if ((firstNumber == lastNumber) || (!someOriginalHitsLeft) || (iterationCounter == 10))
		{
		    hitsChanged = kFALSE;
		}
		iterationCounter++;
	    }

	    cout << "fixed hits" << endl;
	    Double_t
		xSeed = (dSeed + rSeed) * cos(pSeed),
		ySeed = (dSeed + rSeed) * sin(pSeed);
	    
	    if ((possibleHits.size() > 3) && (someOriginalHitsLeft))
	    {
		cout << "book em" << endl;
		// create track
		pTrack = new CbmSttTrack();
		nTmpTracks++;

		pTrack->GetParamLast()->SetX(dSeed);                
		pTrack->GetParamLast()->SetY(pSeed);
		pTrack->GetParamLast()->SetZ(0.);                
		pTrack->GetParamLast()->SetTx(rSeed);
		pTrack->GetParamLast()->SetTy(0.);                
		pTrack->GetParamLast()->SetQp(0.);

		// move hits from accumulator to track
		for (int hitTeller = 0; hitTeller < possibleHits.size(); hitTeller++)
		{
		    pHitFirst = GetHitFromCollections(possibleHits[hitTeller]);
		    if (pHitFirst->GetRadial() == 0.)
		    {
			pHitFirst->SetRadial(sqrt(pHitFirst->GetX() * pHitFirst->GetX() + pHitFirst->GetY() * pHitFirst->GetY()));
		    }

		    pTrack->AddHit(possibleHits[hitTeller], pHitFirst);
		    cout << "removing: " << possibleHits[hitTeller] << endl;
		    fAccumulator->RemoveHit(possibleHits[hitTeller]);
		}

		pTrack->SortHits();

		// draw track seed
		if (rootoutput && 0)
		{
		    TArc  
			myArc;
		    
		    myArc.SetFillStyle(0);
		    myArc.SetLineColor(nTmpTracks + 1);
		    myArc.DrawArc(xSeed, ySeed, rSeed);
		}

		// draw track hits
		if (rootoutput && 0)
		{
		    for (Int_t hitCounter = 0; hitCounter < pTrack->GetNofHits(); hitCounter++)
		    {
			CbmSttHit
			    *pHit = GetHitFromCollections(pTrack->GetHitIndex(hitCounter));

			//cout << "drawing: " << pTrack->GetHitIndex(hitCounter) << endl;
			TArc  
			    myArc;
			
			Double_t
			    radius = pHit->GetIsochrone();

			if (radius == 0)
			    radius = 1.;

			myArc.SetFillStyle(0);
			myArc.SetLineColor(nTmpTracks + 1);
			myArc.DrawArc(pHit->GetX(), pHit->GetY(), radius);
		    }

		    finderCanvas->Update();
		    char
			waitchar;
		    
		    cout << "continue:";
		    cin >> waitchar;

		}

		cout << "track done" << endl;
		tmpTrackArray.push_back(*pTrack);
	    }
	    else
	    {
		for (Int_t iHit = 0; iHit < peakContents.size(); iHit++)
		{
		    Int_t
			iHitFirst = peakContents[iHit];
		    
		    if (iHitFirst != -1)
		    {
			fAccumulator->RemoveHit(iHitFirst, kTRUE, dSeedBin, pSeedBin, rSeedBin);
		    }
		}
	    }
	}
	else
	{
	    for (Int_t iHit = 0; iHit < peakContents.size(); iHit++)
	    {
		Int_t
		    iHitFirst = peakContents[iHit];

		if (iHitFirst != -1)
		{
		    fAccumulator->RemoveHit(iHitFirst, kTRUE, dSeedBin, pSeedBin, rSeedBin);
		}
	    }
	}
	maximumContents = fAccumulator->GetHighestPeakContents();
    }

    Bool_t
	tracksLeft = kTRUE;

    Int_t
	nTracks = 0;

    while (tracksLeft)
    {
	tracksLeft = CreateTracks(&tmpTrackArray, trackArray, nTracks);
    }
    
    if (DO_TIMING)
    {
	timer.Stop();
	totaltime += timer.CpuTime();
	findtime = timer.CpuTime();
    }

    if (rootoutput)
    {
	finderCanvas->Update();
	finderCanvas->Show();
	
	if (saverootoutput)
	{
	    stringstream
		conv;
	    
	    string
		eventCounterStr;
	    
	    conv << fEventCounter;
	    conv >> eventCounterStr;
	    
	    string
		fileName = "finderView" + eventCounterStr + ".C";
	    
	    finderCanvas->SaveSource(fileName.c_str()); 
	}
    }


    cout << "total time: " << totaltime << endl;
    cout << "create time: " << createtime << endl;
    cout << "fill time: " << filltime << endl;
    cout << "find time: " << findtime << endl;
    
    if (rootoutput)
    {
	char
	    goOnChar;

	cout << "press any key to continue: " << endl;
	cin >> goOnChar;

	finderCanvas->Clear();
	
	delete hough;
	delete finderCanvas;
    }

    if (fVerbose) 
    {
	cout << endl;
	cout << "-------------------------------------------------------" 
	     << endl;
	cout << "-I-           Hough STT track finding               -I-"  
	     << endl;
	cout << "Hits: " << nHits << endl;
	cout << "-------------------------------------------------------" 
	     << endl;
    }

    return nTracks;
}

Bool_t CbmSttTrackFinderHough::CreateTracks(list<CbmSttTrack> *tmpTrackArray, TClonesArray *realTrackArray, Int_t &nTracks)
{
    cout << "in create hits" << endl;
    Int_t
	maxHits = 0;

    CbmSttTrack
	*realTrack;

    list<CbmSttTrack>::iterator
	trackIter = tmpTrackArray->begin(),
	tmpTrack;

    cout << "tracks: " << tmpTrackArray->size() << endl;

    while (trackIter != tmpTrackArray->end())
    {
	cout << " hits: " << trackIter->GetNofHits() << endl;
	trackIter->SortHits();

	if (trackIter->GetNofHits() > maxHits)
	{
	    tmpTrack = trackIter;
	    maxHits = tmpTrack->GetNofHits();
	}

	trackIter++;
    }

    if (maxHits == 0)
    {
	return kFALSE;
    }

    Int_t
	numberOfHits = 0;

    for (Int_t hitCounter = 0; hitCounter < tmpTrack->GetNofHits(); hitCounter++)
    {
	Int_t
	    hitNumber = tmpTrack->GetHitIndex(hitCounter);
	
	CbmSttHit
	    *pHit = GetHitFromCollections(hitNumber);
	
	if (!pHit->IsAssigned())
	{
	    numberOfHits++;
	}
    }

    if (numberOfHits > 5)
    {
	// create real track
	realTrack = new((*realTrackArray)[nTracks]) CbmSttTrack();
	nTracks++;
	
	cout << "**************** tmptrack radius: " << tmpTrack->GetParamLast()->GetTx() << endl;
	
	realTrack->GetParamLast()->SetX(tmpTrack->GetParamLast()->GetX());                
	realTrack->GetParamLast()->SetY(tmpTrack->GetParamLast()->GetY());
	realTrack->GetParamLast()->SetZ(tmpTrack->GetParamLast()->GetZ());                
	realTrack->GetParamLast()->SetTx(tmpTrack->GetParamLast()->GetTx());
	realTrack->GetParamLast()->SetTy(tmpTrack->GetParamLast()->GetTy());                
	realTrack->GetParamLast()->SetQp(tmpTrack->GetParamLast()->GetQp());
	
	cout << "starting copying" << endl;
	// copy hits to real track
	for (Int_t hitCounter = 0; hitCounter < tmpTrack->GetNofHits(); hitCounter++)
	{
	    cout << "adding hit" << endl;
	    
	    Int_t
		hitNumber = tmpTrack->GetHitIndex(hitCounter);
	    
	    CbmSttHit
		*pHit = GetHitFromCollections(hitNumber);
	    
	    if (!pHit->IsAssigned())
	    {
		// draw hits
		if (rootoutput)
		{
		    TArc  
			myArc;
		    
		    Double_t
			radius = pHit->GetIsochrone();
		    
		    if (radius == 0)
			radius = 1.;
		    
		    myArc.SetFillStyle(0);
		    myArc.SetLineColor(nTracks + 1);
		    myArc.DrawArc(pHit->GetX(), pHit->GetY(), radius);
		}
		
		pHit->SetAssigned();	
		realTrack->AddHit(hitNumber, pHit);		
	    }
	}
    }
    tmpTrackArray->erase(tmpTrack);
    return kTRUE;
}


void CbmSttTrackFinderHough::ZoomTrack(Double_t &dSeed, Double_t &phiSeed, 
				       Double_t &rSeed, vector<int> &hits)
{
    Int_t
	plotCounter = 0;

    if (hits.size() < 3)
    {
	cout << "-W- CbmSttTrackFinderHough::GetSeed2Circular: less than three hits found, not possible to fit!"
	     << endl;
	return;
    }
    else
    {
	Double_t
	    phi = phiSeed,
	    dist = dSeed;
	
	CbmSttHoughAccumulatorNew
	    *accumulator = new CbmSttHoughAccumulatorNew(2, kTRUE,
							 150, dist - 5., dist + 5., 
							 150, phi - 0.5, phi + 0.5, 
							 100, rSeed - 200., rSeed + 200.); 
	
       	for (Int_t i = 0; i < hits.size() - 2; i++)
	{
	    for (Int_t ii = i + 1; ii < hits.size() - 1; ii++)
	    {
		for (Int_t iii = ii + 1; iii < hits.size(); iii++)
		{
		    CbmSttHit
			*iFirstHit = GetHitFromCollections(hits[i]),
			*iSecondHit = GetHitFromCollections(hits[ii]),
			*iThirdHit = GetHitFromCollections(hits[iii]);
		    
		    Double_t
			firstX = iFirstHit->GetX(),
			firstY = iFirstHit->GetY(),
			firstR = iFirstHit->GetIsochrone(),
			secondX = iSecondHit->GetX(),
			secondY = iSecondHit->GetY(),
			secondR = iSecondHit->GetIsochrone(),
			thirdX = iThirdHit->GetX(),
			thirdY = iThirdHit->GetY(),
			thirdR = iThirdHit->GetIsochrone();

		    
		    Double_t
			circleRadii[8],
			circleCentersX[8],
			circleCentersY[8];
		    
		    GetTrackletCircular(firstX, firstY, firstR, 
					secondX, secondY, secondR, 
					thirdX, thirdY, thirdR, 
					circleRadii, circleCentersX, circleCentersY);

		    // loop over all tracklets between this and next
		    for (Int_t trackletteller = 0; trackletteller < 8; trackletteller++)
		    {
			Double_t
			    phi = atan(circleCentersY[trackletteller] / circleCentersX[trackletteller]),
			    dist = sqrt(circleCentersX[trackletteller] * circleCentersX[trackletteller] +
					circleCentersY[trackletteller] * circleCentersY[trackletteller])
			    - circleRadii[trackletteller];
			
			if (circleCentersX[trackletteller] < 0.)
			{
			    if (circleCentersY[trackletteller] < 0.)
				phi -= 3.1415;
			    else
				phi += 3.1415;
			}
			
			
			accumulator->AddHit(dist, i, 
					    phi, ii, 
					    circleRadii[trackletteller], iii);
		    }
		}
	    }
	}

	// cluster accumulator
	accumulator->Cluster();
	
	dSeed = accumulator->GetHighestPeakX();
	phiSeed = accumulator->GetHighestPeakY();
	rSeed = accumulator->GetHighestPeakZ();

	Double_t
	    xSeed = (rSeed + dSeed) * cos(phiSeed),
	    ySeed = (rSeed + dSeed) * sin(phiSeed);

	// the fine track estimation
	if (0)
	{
	    TArc  
		myArc;
	    
	    myArc.SetFillStyle(0);
	    myArc.SetLineColor(5);
	    myArc.DrawArc(xSeed, ySeed, rSeed);

	    finderCanvas->Update();

	    char
		waitchar;

	    cout << "zoomed" << endl;
	    cin >> waitchar;
	}

	vector<int>
	    goodHits;

	for (int hitTeller = 0; hitTeller < hits.size(); hitTeller++)
	{
	    if (GetChiSquare(xSeed, ySeed, rSeed, GetHitFromCollections(hits[hitTeller])) < 20.0)
	    {
		goodHits.push_back(hits[hitTeller]);
	    }
	}

	hits.clear();

	for (int hitTeller = 0; hitTeller < goodHits.size(); hitTeller++)
	{
	    hits.push_back(goodHits[hitTeller]);
	}

	delete accumulator;
    }
}

CbmSttHit* CbmSttTrackFinderHough::GetHitFromCollections(Int_t hitCounter)
{
    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 CbmSttTrackFinderHough::AddMissingHits(Double_t dSeed, Double_t pSeed, Double_t rSeed, vector<int> &possibleHits)
{
    Double_t
	xSeed = (dSeed + rSeed) * cos(pSeed),
	ySeed = (dSeed + rSeed) * sin(pSeed);

    Int_t 
	nHits = 0;

    for (Int_t hitListCounter = 0; hitListCounter < fHitCollectionList.GetEntries(); hitListCounter++)
    {
	nHits += ((TClonesArray *)fHitCollectionList.At(hitListCounter))->GetEntriesFast();
    }

    Double_t
	radialCoordinates[possibleHits.size()],
	zCoordinates[possibleHits.size()];

    sort(possibleHits.begin(), possibleHits.end());

   // close in r-z plane
    for (Int_t hitTeller = 0; hitTeller < possibleHits.size(); hitTeller++)
    {
	CbmSttHit
	    *tmpHit = GetHitFromCollections(possibleHits[hitTeller]);
	
	radialCoordinates[hitTeller] = sqrt(tmpHit->GetX() * tmpHit->GetX() + tmpHit->GetY() * tmpHit->GetY());
	zCoordinates[hitTeller] = tmpHit->GetZ();
    }
    /*
    TCanvas
	tmpCanvas("tmp", "tmp", 800, 600);
    */
    TGraph
	fitGraph(possibleHits.size(), radialCoordinates, zCoordinates);
    /*
    TH2F myRange("range", "range", 100, -100., 100., 100, -100., 100);
    myRange.DrawCopy();
    
    fitGraph.DrawGraph(possibleHits.size(), radialCoordinates, zCoordinates);
    */
    TF1 
	lineFunc("linefunc", "pol1");
   
    if (possibleHits.size() > 3)
    {
	fitGraph.Fit("linefunc", "");
    }

    Double_t
	slope = lineFunc.GetParameter(1),
	offset = lineFunc.GetParameter(0);

/*    
    TLine
	myLine;
    
    myLine.DrawLine(-100, -100 * slope + offset, 100, 100 * slope + offset);

    tmpCanvas.Show();
    tmpCanvas.Update();
	    
    char
	waitchar;
    
    cin >> waitchar;
*/
	
    for (Int_t iHit = 0; iHit < nHits; iHit++)
    {
	CbmSttHit
	    *pHit = GetHitFromCollections(iHit);

	Double_t
	    radial = sqrt(pHit->GetX() * pHit->GetX() + pHit->GetY() * pHit->GetY()),
	    zDifference = radial * slope + offset - pHit->GetZ();

	// close in x-y plane
	if ((GetChiSquare(xSeed, ySeed, rSeed, pHit) < 20.) && (sqrt(zDifference * zDifference) < 20.))
	{
	    Bool_t
		addFlag = kTRUE;

	    for (int posHitTeller = 0; posHitTeller < possibleHits.size(); posHitTeller++)
	    {
		if (possibleHits[posHitTeller] == iHit)
		{
		    addFlag = kFALSE;
		    break;
		}
	    }

	    // check for closeness to track
	    if (addFlag)
	    {
		possibleHits.push_back(iHit);
	    }
	}
    }
}


CbmMCPoint* CbmSttTrackFinderHough::GetPointFromCollections(Int_t hitCounter)
{
    return NULL;
}

Double_t CbmSttTrackFinderHough::GetChiSquare(Double_t xCircleCenter, Double_t yCircleCenter, Double_t rCircle, CbmSttHit *pMhit)
{
    Double_t 
	chiSquare = 0.,
	measuredRadius = 0.,
	radialResolution = 0.,
	distanceOfClosestApproach = 0.,
	x_0 = 0.,
	y_0 = 0.; 
    
    // get isochrone
    measuredRadius = pMhit->GetIsochrone();
    radialResolution = pMhit->GetIsochroneError();

    if (radialResolution == 0)
    {
	radialResolution = sqrt(pMhit->GetDx() * pMhit->GetDx() + pMhit->GetDy() * pMhit->GetDy() + pMhit->GetDz() * pMhit->GetDz());
    }

    // get XY of wire
    x_0 = pMhit->GetX();
    y_0 = pMhit->GetY();
 
    distanceOfClosestApproach = sqrt((x_0 - xCircleCenter) * (x_0 - xCircleCenter) + 
				     (y_0 - yCircleCenter) * (y_0 - yCircleCenter)) - rCircle;
    
    distanceOfClosestApproach = sqrt(distanceOfClosestApproach * distanceOfClosestApproach);
    
    chiSquare = (distanceOfClosestApproach - measuredRadius) / radialResolution;
  
    return chiSquare;
}

/*
Int_t CbmSttTrackFinderHough::DoFindLinear(TClonesArray* mHitArray,
					   TClonesArray* trackArray) 

{
    fEventCounter++;

    if (rootoutput)
    {
      finderCanvas = new TCanvas("findercanvas", "findercanvas", 800, 600); 
      TH2F myRange("range", "range", 100, -50., 50., 100, -50., 50);
      myRange.DrawCopy();
      plotAllStraws();
      //finderCanvas->Update();
    }
 
    TH2F
	*hough;
    
    if (rootoutput)
    {
	cout << "creating hough" << endl;
	hough = new TH2F("hough", "hough", CbmSttHoughThetaBins, -100., 100., CbmSttHoughPhiBins, -1 * dPi, dPi);
    }

    
    fAccumulator->ClearAccumulator();
    fAccumulator->SetThreshold(CbmSttHoughThreshold);
    
    if (!mHitArray)
    {
	cout << "-E- CbmSttTrackFinderHough::DoFind: "
	     << "Hit array missing! " 
	     << mHitArray << endl;
	return -1;
    }
    
    if ( !trackArray ) 
    {
	cout << "-E- CbmSttTrackFinderHough::DoFind: "
	     << "Track array missing! " << endl;
	return -1;
    }
    
    // Create pointers to hit and SttPoint
    CbmSttHit*       pHitFirst = NULL;
    CbmSttHit*       pHitSecond = NULL;
    CbmSttTrack*     pTrack = NULL;
    
    // Number of STT hits
    Int_t nHits = mHitArray->GetEntriesFast();
    
    Int_t
	ptIndexFirst = 0,
	ptIndexSecond = 0;
    
    // fill accumulator
    for (Int_t iHitFirst = 0; iHitFirst < nHits; iHitFirst++)
    {
	pHitFirst = (CbmSttHit*) mHitArray->At(iHitFirst);
	
	if ( ! pHitFirst ) 
	    continue;
	
	ptIndexFirst = pHitFirst->GetRefIndex();
	
	if (ptIndexFirst < 0)
	    continue;           // fake or background hit
	
	Double_t
	    thisX = pHitFirst->GetX(),
	    thisY = pHitFirst->GetY(),
	    thisR = pHitFirst->GetIsochrone();
	
	if (rootoutput)
	{
	    TArc  
		myArc;
	    
	    myArc.SetLineColor(1);
	    myArc.DrawArc(thisX, thisY, thisR);
	}
	
	for (Int_t iHitSecond = iHitFirst + 1; iHitSecond < nHits; iHitSecond++)
	{
	    pHitSecond = (CbmSttHit*) mHitArray->At(iHitSecond);
	    
	    if ( ! pHitSecond ) 
		continue;
	    
	    ptIndexSecond = pHitSecond->GetRefIndex();
	    
	    if (ptIndexSecond < 0)
		continue;           // fake or background hit
	    
	    Double_t
		nextX = pHitSecond->GetX(),
		nextR = pHitSecond->GetIsochrone(),
		nextY = pHitSecond->GetY();
	    
	    Double_t
		contactPointsXWithNext[4],
		contactPointsYWithNext[4],
		contactPointsXAtNext[4],
		contactPointsYAtNext[4],
		anglesWithNext[4];
	    
	    GetTracklet(thisX, thisY, thisR, nextX, nextY, nextR, contactPointsXWithNext, 
			contactPointsYWithNext, anglesWithNext, contactPointsXAtNext, contactPointsYAtNext); 
	    
	    // loop over all tracklets between this and next
	    for (int trackletteller = 0; trackletteller < 4; trackletteller++)
	    {
		double
		    x_1 = contactPointsXWithNext[trackletteller],
		    x_2 = contactPointsXAtNext[trackletteller],
		    y_1 = contactPointsYWithNext[trackletteller],
		    y_2 = contactPointsYAtNext[trackletteller],
		    thetaval = atan((x_2 - x_1) / (y_2 - y_1)),
		    distval = ((( x_2 - x_1 ) * y_1 - x_1 * ( y_2 - y_1)) / 
			       (sqrt((x_2 - x_1) * (x_2 - x_1) + (y_2 - y_1) * (y_2 - y_1))));
		
		distval = sqrt(distval * distval);
		
		if (thetaval < 0)
		    thetaval += dPi;
		
		fAccumulator->AddHit(distval, iHitFirst, thetaval, iHitSecond);
		if (rootoutput)
		{
		    //hough->Fill(distval, thetaval);
		}
	    }
	}
    }
    
    Int_t
	nTracks = 0;

    Int_t
	maximumContents;

    fAccumulator->Cluster();
    // get all peaks, sort by height
    maximumContents = fAccumulator->GetHighestPeakContents();

    while (maximumContents > 0)
    {
	// create track
	pTrack = new((*trackArray)[nTracks]) CbmSttTrack();
	nTracks++;

	cout << "ntrack: " << nTracks << "  " << maximumContents << endl;
	Double_t
	    refYsign,
	    refXsign;

	Bool_t
	    foundFirstTrackHit = kFALSE;

	// loop over houghCell
	for (Int_t iHit= 0; iHit < maximumContents; iHit++)
	{
	    Int_t
		iHitFirst,
		iHitSecond;

	    fAccumulator->GetHighestPeakEntry(iHit, iHitFirst, iHitSecond);
	    if ((iHitFirst != -1) && (iHitSecond!= -1))
	    {
		pHitFirst = (CbmSttHit*) mHitArray->At(iHitFirst);
		pHitSecond = (CbmSttHit*) mHitArray->At(iHitSecond);
		
		// if first trackhit, store phi
		if (!foundFirstTrackHit)
		{
		    refYsign = pHitFirst->GetY();
		    refXsign = pHitFirst->GetX();
		    
		    foundFirstTrackHit = kTRUE;
		}
		
		if ((pHitFirst->GetX() * refXsign >= 0.) && (pHitFirst->GetY() * refYsign >= 0.)) 
		{
		    pTrack->AddHit(iHitFirst, pHitFirst);
		    fAccumulator->RemoveHit(iHitFirst);
		    
		    if (rootoutput)
		    {
			TArc  
			    myArc;
			
			myArc.SetLineColor(nTracks + 1);
			myArc.DrawArc(pHitFirst->GetX(), pHitFirst->GetY(), pHitFirst->GetIsochrone());
		    }
		}
		
		if ((pHitSecond->GetX() * refXsign >= 0.) && (pHitSecond->GetY() * refYsign >= 0.)) 
		{
		    pTrack->AddHit(iHitSecond, pHitSecond);
		    fAccumulator->RemoveHit(iHitSecond);
		    
		    if (rootoutput)
		    {
			TArc  
			    myArc;
			
			myArc.SetLineColor(nTracks + 1);
			myArc.DrawArc(pHitSecond->GetX(), pHitSecond->GetY(), pHitSecond->GetIsochrone());
		    }
		}
	    }
	}

	maximumContents = fAccumulator->GetHighestPeakContents();
    }
	
    cout << "done this" << endl;
    if (rootoutput)
    {
	cout << "getting here1" << endl;
      char
	  goOnChar;

	cout << "getting here2" << endl;
      hough->SetFillColor(2);
	cout << "getting here3" << endl;
	//hough->DrawCopy();
	cout << "getting here4" << endl;
      finderCanvas->Update();
	cout << "getting here5" << endl;
      finderCanvas->Show();
	cout << "getting here6" << endl;

      stringstream
	  conv;
 	cout << "getting here7" << endl;
     
      string
	  eventCounterStr;
	
	cout << "getting here8" << endl;
      conv << fEventCounter;
      conv >> eventCounterStr;
	
      string
	  fileName = "finderView" + eventCounterStr + ".C";

      finderCanvas->SaveSource(fileName.c_str()); 
      cout << "press any key to continue: " << endl;
      cin >> goOnChar;

      finderCanvas->Clear();
      
      delete hough;
      delete finderCanvas;
    }
    cout << "done root" << endl;
    //if (fVerbose) 
    {
	cout << endl;
	cout << "-------------------------------------------------------" 
	     << endl;
	cout << "-I-           Hough STT track finding               -I-"  
	     << endl;
	cout << "Hits: " << nHits << endl;
	cout << "-------------------------------------------------------" 
	     << endl;
    }
    
    return nTracks;
}
*/
Int_t CbmSttTrackFinderHough::NumberOfLineCircleIntersections(Double_t xline1, Double_t yline1, 
							      Double_t xline2, Double_t yline2, 
							      Double_t xcircle, Double_t ycircle, 
							      Double_t rcircle) const
{
    Int_t
	retval = 0;
    
    xcircle -= xline1;
    ycircle -= yline1;
    xline2 -= xline1;
    yline2 -= yline1;
  
    Double_t
	cline = yline2 / xline2;
    
    Double_t
	determinant;

    determinant = 8 * xcircle * cline * ycircle 
	-4 * ycircle * ycircle 
	+4 * rcircle * rcircle 
	-4 * cline * cline * xcircle * xcircle 
	+4 * cline * cline * rcircle * rcircle;

    if (determinant < -1 * limit)
	retval = 0;
    else if (fabs(determinant) < limit)
	retval = 1;
    else
	retval = 2;

    return retval;
}

// -------------------------------------------------------------------------


ClassImp(CbmSttTrackFinderHough)