[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTrack.cxx

//$Id$

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ASV 

#include "AliL3StandardIncludes.h"

#include "AliL3Logging.h"
#include "AliL3HoughTrack.h"
#include "AliL3Transform.h"

#if GCCVERSION == 3
using namespace std;
#endif

//_____________________________________________________________
// AliL3HoughTrack
//
// Track class for Hough tracklets

ClassImp(AliL3HoughTrack)


AliL3HoughTrack::AliL3HoughTrack()
{
  //Constructor
  
  fWeight = 0;
  fMinDist=0;
  fDLine = 0;
  fPsiLine = 0;
  fIsHelix = true;
  fEtaIndex = -1;
  fEta = 0;
  
}


AliL3HoughTrack::~AliL3HoughTrack()
{
  
}

void AliL3HoughTrack::Set(AliL3Track *track)
{
  
  AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
  SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
  SetEtaIndex(tpt->GetEtaIndex());
  SetEta(tpt->GetEta());
  SetTgl(tpt->GetTgl());
  SetPsi(tpt->GetPsi());
  SetCenterX(tpt->GetCenterX());
  SetCenterY(tpt->GetCenterY());
  SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
  SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
  SetCharge(tpt->GetCharge());
  SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  SetSlice(tpt->GetSlice());
  SetNHits(1);
  return;

  fWeight = tpt->GetWeight();
  fDLine = tpt->GetDLine();
  fPsiLine = tpt->GetPsiLine();
  SetNHits(tpt->GetWeight());
  SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  fIsHelix = false;
}

Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
{
  AliL3HoughTrack *track = (AliL3HoughTrack*)tpt;
  if(track->GetWeight() < GetWeight()) return 1;
  if(track->GetWeight() > GetWeight()) return -1;
  return 0;
}

void AliL3HoughTrack::SetEta(Double_t f)
{
  //Set eta, and calculate fTanl, which is the tan of dipangle

  fEta = f;
  Double_t theta = 2*atan(exp(-1.*fEta));
  Double_t dipangle = AliL3Transform::Pi()/2 - theta;
  Double_t tgl = tan(dipangle);
  SetTgl(tgl);
}

void AliL3HoughTrack::CalculateHelix()
{
  return;
}

void AliL3HoughTrack::UpdateToFirstRow()
{
  //Update the track parameters to the point where track cross
  //its first padrow.
  
  //Get the crossing point with the first padrow:
  Float_t xyz[3];
  if(!GetCrossingPoint(GetFirstRow(),xyz))
    LOG(AliL3Log::kWarning,"AliL3HoughTrack::UpdateToFirstRow()","Track parameters")
      <<AliL3Log::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
      <<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
  
  //printf("Track with eta %f tgl %f crosses at x %f y %f z %f on padrow %d\n",GetEta(),GetTgl(),xyz[0],xyz[1],xyz[2],GetFirstRow());
  //printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
  
  Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);

  //Get the track parameters
  
  /*
    Double_t x0    = GetR0() * cos(GetPhi0()) ;
    Double_t y0    = GetR0() * sin(GetPhi0()) ;
  */
  Double_t rc    = GetRadius();//fabs(GetPt()) / ( BFACT * AliL3Transform::GetBField() )  ;
  Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / abs(GetCharge()) ;
  Double_t xc    = GetCenterX();//x0 - rc * cos(tPhi0) ;
  Double_t yc    = GetCenterY();//y0 - rc * sin(tPhi0) ;
  
  //Check helix and cylinder intersect
  Double_t fac1 = xc*xc + yc*yc ;
  Double_t sfac = sqrt( fac1 ) ;
  
  if ( fabs(sfac-rc) > radius || fabs(sfac+rc) < radius ) {
    LOG(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
      "Track does not intersect"<<ENDLOG;
    return;
  }
  
  //Find intersection
  Double_t fac2   = (radius*radius + fac1 - rc*rc) / (2.00 * radius * sfac ) ;
  Double_t phi    = atan2(yc,xc) + GetCharge()*acos(fac2) ;
  Double_t td     = atan2(radius*sin(phi) - yc,radius*cos(phi) - xc) ;
  
  //Intersection in z
  if ( td < 0 ) td = td + 2. * pi ;
  Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),2*pi) ;
  if ( deltat < 0.      ) deltat += 2. * pi ;
  if ( deltat > 2.*pi ) deltat -= 2. * pi ;
  Double_t z = GetZ0() + rc * GetTgl() * deltat ;
  
  
  Double_t xExtra = radius * cos(phi) ;
  Double_t yExtra = radius * sin(phi) ;
  
  Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
  
  //if ( tPhi < 0 ) tPhi += 2. * M_PI ;
  
  Double_t tPsi = tPhi - GetCharge() * 0.5 * pi / abs(GetCharge()) ;
  if ( tPsi > 2. * pi ) tPsi -= 2. * pi ;
  if ( tPsi < 0.        ) tPsi += 2. * pi ;
  
  //And finally, update the track parameters
  SetR0(radius);
  SetPhi0(phi);
  SetZ0(z);
  SetPsi(tPsi);
  SetFirstPoint(xyz[0],xyz[1],z);
  //printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
  
  //printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
  
  //Also, set the coordinates of the point where track crosses last padrow:
  GetCrossingPoint(GetLastRow(),xyz);
  SetLastPoint(xyz[0],xyz[1],xyz[2]);
  //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
}

void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
{

  fWeight = weight;
  fMinDist = 100000;
  SetKappa(kappa);
  Double_t pt = fabs(BFACT*AliL3Transform::GetBField()/kappa);
  SetPt(pt);
  Double_t radius = 1/fabs(kappa);
  SetRadius(radius);
  SetFirstPoint(0,0,0);
  SetPsi(eangle); //Psi = emission angle when first point is vertex
  SetPhi0(0);     //not defined for vertex reference point
  SetR0(0);
  Double_t charge = -1.*kappa;
  SetCharge((Int_t)copysign(1.,charge));
  Double_t trackPhi0 = GetPsi() + charge*0.5*AliL3Transform::Pi()/fabs(charge);
  Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ;
  Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ;
  SetCenterX(xc);
  SetCenterY(yc);
  SetNHits(1); //just for the trackarray IO
  fIsHelix = true;
}

void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row)
{
  //Initialize a track piece, not yet a track
  //Used in case of straight line transformation

  //Transform line parameters to coordinate system of slice:
  
  /*
  D = D + fTransform->Row2X(ref_row)*cos(psi);

  fDLine = D;
  fPsiLine = psi;
  fWeight = weight;
  SetNHits(weight);
  SetRowRange(rowrange[0],rowrange[1]);
  fIsHelix = false;
  */
}

void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t min_dist)
{
  
  if(min_dist < fMinDist)
    {
      fMinDist = min_dist;
      SetMCid(mcid);
    }
  
}

void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Double_t *xy)
{
  
  /*
  if(fIsHelix)
    {
      printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
      return;
    }

  Double_t xhit = fTransform->Row2X(padrow);
  Double_t a = -1/tan(fPsiLine);
  Double_t b = fDLine/sin(fPsiLine);
  
  Double_t yhit = a*xhit + b;
  xy[0] = xhit;
  xy[1] = yhit;
  */
}

/*
Double_t AliL3HoughTrack::GetCrossingAngle(Int_t padrow)
{
  //Calculate the crossing angle between track and given padrow.

  if(!fIsHelix)
    {
      printf("AliL3HoughTrack::GetCrossingAngle : Track is not a helix\n");
      return 0;
    }

  if(!IsLocal())
    {
      printf("Track is not given in local coordinates\n");
      return 0;
    }

  Float_t xyz[3];
  if(!GetCrossingPoint(padrow,xyz))
    printf("AliL3HoughTrack::GetCrossingPoint : Track does not cross line!!\n");
  
  
  //Convert center of curvature to local coordinates:
  //Float_t xyz_coc[3] = {GetCenterX(),GetCenterY(),0};
  //fTransform->Global2Local(xyz_coc,slice);
  
  //Take the dot product of the tangent vector of the track, and
  //vector perpendicular to the padrow.
  
  Double_t tangent[2];
  //tangent[1] = (xyz[0] - xyz_coc[0])/GetRadius();
  //tangent[0] = -1.*(xyz[1] - xyz_coc[1])/GetRadius();
  tangent[1] = (xyz[0] - GetCenterX())/GetRadius();
  tangent[0] = -1.*(xyz[1] - GetCenterY())/GetRadius();

  Double_t perp_padrow[2] = {1,0}; //locally in slice

  Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
  return acos(cos_beta);
  
}

Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t padrow,Float_t *xyz)
{
  //Assumes the track is given in local coordinates

  if(!fIsHelix)
    {
      printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
      return 0;
    }
    

  if(!IsLocal())
    {
      printf("GetCrossingPoint: Track is given on global coordinates\n");
      return false;
    }
  
  Double_t xHit = fTransform->Row2X(padrow);

  //xyz[0] = fTransform->Row2X(padrow);
  xyz[0] = xHit;
  Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
  Double_t r2 = GetRadius()*GetRadius();
  if(aa > r2)
    return false;

  Double_t aa2 = sqrt(r2 - aa);
  Double_t y1 = GetCenterY() + aa2;
  Double_t y2 = GetCenterY() - aa2;
  xyz[1] = y1;
  if(fabs(y2) < fabs(y1)) xyz[1] = y2;
  xyz[2] = 0; //only consider transverse plane
  
  return true;
}


Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t slice,Int_t padrow,Float_t *xyz)
{
  //Calculate the crossing point in transverse plane of this track and given 
  //padrow (y = a*x + b). Point is given in local coordinates in slice.
  //Assumes the track is given in global coordinates.

  if(!fIsHelix)
    {
      printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
      return 0;
    }

  
  if(IsLocal())
    {
      printf("GetCrossingPoint: Track is given in local coordintes!!!\n");
      return false;
    }

  Double_t padrowradii = fTransform->Row2X(padrow);


  Float_t rotation_angle = (slice*20)*ToRad;
  
  Float_t cs,sn;
  cs = cos(rotation_angle);
  sn = sin(rotation_angle);

  Double_t a = -1.*cs/sn;
  Double_t b = padrowradii/sn;

  Double_t ycPrime = GetCenterY() - b ;
  Double_t aa = ( 1. + a * a ) ;
  Double_t bb = -2. * ( GetCenterX() + a * ycPrime ) ;
  Double_t cc = ( GetCenterX() * GetCenterX() + ycPrime * ycPrime - GetRadius() * GetRadius() ) ;

  Double_t racine = bb * bb - 4. * aa * cc ;
  if ( racine < 0 ) return false ;
  Double_t rootRacine = sqrt(racine) ;
  
  Double_t oneOverA = 1./aa;
//
//   First solution
//
   Double_t x1 = 0.5 * oneOverA * ( -1. * bb + rootRacine ) ; 
   Double_t y1 = a * x1 + b ;
   Double_t r1 = sqrt(x1*x1+y1*y1);
//
//   Second solution
//
   Double_t x2 = 0.5 * oneOverA * ( -1. * bb - rootRacine ) ; 
   Double_t y2 = a * x2 + b ;
   Double_t r2 = sqrt(x2*x2+y2*y2);
//
//    Choose close to (0,0) 
//
   Double_t xHit ;
   Double_t yHit ;
   if ( r1 < r2 ) {
      xHit = x1 ;
      yHit = y1 ;
   }
   else {
      xHit = x2 ;
      yHit = y2 ;
   }
  
   xyz[0] = xHit;
   xyz[1] = yHit;
   xyz[2] = 0;
   
   fTransform->Global2Local(xyz,slice);
   
   return true;
}
*/
Commit	Line	Data
b1886074	1	//$Id$
	2
	3	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	4	//*-- Copyright &copy ASV
d712b5b8	5
e06900d5	6	#include "AliL3StandardIncludes.h"
d712b5b8	7
d712b5b8	8	#include "AliL3Logging.h"
d712b5b8	9	#include "AliL3HoughTrack.h"
13a0c3d8	10	#include "AliL3Transform.h"
d712b5b8	11
e06900d5	12	#if GCCVERSION == 3
	13	using namespace std;
	14	#endif
	15
b1886074	16	//_____________________________________________________________
	17	// AliL3HoughTrack
	18	//
	19	// Track class for Hough tracklets
	20
d712b5b8	21	ClassImp(AliL3HoughTrack)
	22
	23
	24	AliL3HoughTrack::AliL3HoughTrack()
	25	{
	26	//Constructor
	27
	28	fWeight = 0;
	29	fMinDist=0;
d712b5b8	30	fDLine = 0;
	31	fPsiLine = 0;
	32	fIsHelix = true;
	33	fEtaIndex = -1;
208b54c5	34	fEta = 0;
208b54c5	35
d712b5b8	36	}
	37
	38
	39	AliL3HoughTrack::~AliL3HoughTrack()
	40	{
237d3f5c	41
d712b5b8	42	}
	43
	44	void AliL3HoughTrack::Set(AliL3Track *track)
	45	{
	46
	47	AliL3HoughTrack tpt = (AliL3HoughTrack)track;
b46b53c1	48	SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
d712b5b8	49	SetEtaIndex(tpt->GetEtaIndex());
d712b5b8	50	SetEta(tpt->GetEta());
ae4d5136	51	SetTgl(tpt->GetTgl());
d712b5b8	52	SetPsi(tpt->GetPsi());
	53	SetCenterX(tpt->GetCenterX());
	54	SetCenterY(tpt->GetCenterY());
	55	SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
	56	SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
	57	SetCharge(tpt->GetCharge());
	58	SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
b1886074	59	SetSlice(tpt->GetSlice());
242c143e	60	SetNHits(1);
d712b5b8	61	return;
	62
	63	fWeight = tpt->GetWeight();
	64	fDLine = tpt->GetDLine();
	65	fPsiLine = tpt->GetPsiLine();
	66	SetNHits(tpt->GetWeight());
	67	SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
	68	fIsHelix = false;
d712b5b8	69	}
	70
	71	Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
	72	{
	73	AliL3HoughTrack track = (AliL3HoughTrack)tpt;
	74	if(track->GetWeight() < GetWeight()) return 1;
	75	if(track->GetWeight() > GetWeight()) return -1;
	76	return 0;
	77	}
	78
	79	void AliL3HoughTrack::SetEta(Double_t f)
	80	{
b1886074	81	//Set eta, and calculate fTanl, which is the tan of dipangle
b1886074	82
d712b5b8	83	fEta = f;
d712b5b8	84	Double_t theta = 2atan(exp(-1.fEta));
26abc209	85	Double_t dipangle = AliL3Transform::Pi()/2 - theta;
8b007a2d	86	Double_t tgl = tan(dipangle);
d712b5b8	87	SetTgl(tgl);
	88	}
	89
	90	void AliL3HoughTrack::CalculateHelix()
	91	{
	92	return;
	93	}
	94
	95	void AliL3HoughTrack::UpdateToFirstRow()
	96	{
	97	//Update the track parameters to the point where track cross
	98	//its first padrow.
	99
	100	//Get the crossing point with the first padrow:
	101	Float_t xyz[3];
	102	if(!GetCrossingPoint(GetFirstRow(),xyz))
	103	LOG(AliL3Log::kWarning,"AliL3HoughTrack::UpdateToFirstRow()","Track parameters")
	104	<<AliL3Log::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
	105	<<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
	106
8b007a2d	107	//printf("Track with eta %f tgl %f crosses at x %f y %f z %f on padrow %d\n",GetEta(),GetTgl(),xyz[0],xyz[1],xyz[2],GetFirstRow());
	108	//printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
	109
d712b5b8	110	Double_t radius = sqrt(xyz[0]xyz[0] + xyz[1]xyz[1]);
	111
	112	//Get the track parameters
b1886074	113
	114	/*
	115	Double_t x0 = GetR0() * cos(GetPhi0()) ;
	116	Double_t y0 = GetR0() * sin(GetPhi0()) ;
	117	*/
13a0c3d8	118	Double_t rc = GetRadius();//fabs(GetPt()) / ( BFACT * AliL3Transform::GetBField() ) ;
e06900d5	119	Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / abs(GetCharge()) ;
b1886074	120	Double_t xc = GetCenterX();//x0 - rc * cos(tPhi0) ;
b1886074	121	Double_t yc = GetCenterY();//y0 - rc * sin(tPhi0) ;
d712b5b8	122
	123	//Check helix and cylinder intersect
	124	Double_t fac1 = xcxc + ycyc ;
	125	Double_t sfac = sqrt( fac1 ) ;
b1886074	126
d712b5b8	127	if ( fabs(sfac-rc) > radius \|\| fabs(sfac+rc) < radius ) {
	128	LOG(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
	129	"Track does not intersect"<<ENDLOG;
	130	return;
	131	}
	132
	133	//Find intersection
b46b53c1	134	Double_t fac2 = (radiusradius + fac1 - rcrc) / (2.00 * radius * sfac ) ;
d712b5b8	135	Double_t phi = atan2(yc,xc) + GetCharge()*acos(fac2) ;
	136	Double_t td = atan2(radiussin(phi) - yc,radiuscos(phi) - xc) ;
	137
	138	//Intersection in z
	139	if ( td < 0 ) td = td + 2. * pi ;
	140	Double_t deltat = fmod((-GetCharge()td + GetCharge()tPhi0),2*pi) ;
	141	if ( deltat < 0. ) deltat += 2. * pi ;
	142	if ( deltat > 2.pi ) deltat -= 2. pi ;
	143	Double_t z = GetZ0() + rc * GetTgl() * deltat ;
8b007a2d	144
d712b5b8	145
	146	Double_t xExtra = radius * cos(phi) ;
	147	Double_t yExtra = radius * sin(phi) ;
	148
	149	Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
	150
	151	//if ( tPhi < 0 ) tPhi += 2. * M_PI ;
	152
e06900d5	153	Double_t tPsi = tPhi - GetCharge() * 0.5 * pi / abs(GetCharge()) ;
d712b5b8	154	if ( tPsi > 2. * pi ) tPsi -= 2. * pi ;
	155	if ( tPsi < 0. ) tPsi += 2. * pi ;
	156
	157	//And finally, update the track parameters
d712b5b8	158	SetR0(radius);
	159	SetPhi0(phi);
	160	SetZ0(z);
	161	SetPsi(tPsi);
	162	SetFirstPoint(xyz[0],xyz[1],z);
8b007a2d	163	//printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
	164
	165	//printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
d712b5b8	166
	167	//Also, set the coordinates of the point where track crosses last padrow:
	168	GetCrossingPoint(GetLastRow(),xyz);
	169	SetLastPoint(xyz[0],xyz[1],xyz[2]);
8b007a2d	170	//printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
d712b5b8	171	}
d712b5b8	172
b46b53c1	173	void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
d712b5b8	174	{
	175
	176	fWeight = weight;
	177	fMinDist = 100000;
	178	SetKappa(kappa);
13a0c3d8	179	Double_t pt = fabs(BFACT*AliL3Transform::GetBField()/kappa);
d712b5b8	180	SetPt(pt);
	181	Double_t radius = 1/fabs(kappa);
	182	SetRadius(radius);
	183	SetFirstPoint(0,0,0);
b46b53c1	184	SetPsi(eangle); //Psi = emission angle when first point is vertex
b46b53c1	185	SetPhi0(0); //not defined for vertex reference point
8b007a2d	186	SetR0(0);
d712b5b8	187	Double_t charge = -1.*kappa;
8b007a2d	188	SetCharge((Int_t)copysign(1.,charge));
26abc209	189	Double_t trackPhi0 = GetPsi() + charge0.5AliL3Transform::Pi()/fabs(charge);
8b007a2d	190	Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ;
8b007a2d	191	Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ;
d712b5b8	192	SetCenterX(xc);
d712b5b8	193	SetCenterY(yc);
242c143e	194	SetNHits(1); //just for the trackarray IO
d712b5b8	195	fIsHelix = true;
	196	}
	197
	198	void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row)
	199	{
	200	//Initialize a track piece, not yet a track
	201	//Used in case of straight line transformation
	202
	203	//Transform line parameters to coordinate system of slice:
	204
237d3f5c	205	/*
d712b5b8	206	D = D + fTransform->Row2X(ref_row)*cos(psi);
	207
	208	fDLine = D;
	209	fPsiLine = psi;
	210	fWeight = weight;
	211	SetNHits(weight);
	212	SetRowRange(rowrange[0],rowrange[1]);
	213	fIsHelix = false;
237d3f5c	214	*/
d712b5b8	215	}
	216
	217	void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t min_dist)
	218	{
	219
	220	if(min_dist < fMinDist)
	221	{
	222	fMinDist = min_dist;
	223	SetMCid(mcid);
	224	}
	225
	226	}
	227
	228	void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Double_t *xy)
	229	{
	230
237d3f5c	231	/*
d712b5b8	232	if(fIsHelix)
	233	{
	234	printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
	235	return;
	236	}
	237
	238	Double_t xhit = fTransform->Row2X(padrow);
	239	Double_t a = -1/tan(fPsiLine);
	240	Double_t b = fDLine/sin(fPsiLine);
	241
	242	Double_t yhit = a*xhit + b;
	243	xy[0] = xhit;
	244	xy[1] = yhit;
237d3f5c	245	*/
d712b5b8	246	}
	247
	248	/*
	249	Double_t AliL3HoughTrack::GetCrossingAngle(Int_t padrow)
	250	{
	251	//Calculate the crossing angle between track and given padrow.
	252
	253	if(!fIsHelix)
	254	{
	255	printf("AliL3HoughTrack::GetCrossingAngle : Track is not a helix\n");
	256	return 0;
	257	}
	258
	259	if(!IsLocal())
	260	{
	261	printf("Track is not given in local coordinates\n");
	262	return 0;
	263	}
	264
	265	Float_t xyz[3];
	266	if(!GetCrossingPoint(padrow,xyz))
	267	printf("AliL3HoughTrack::GetCrossingPoint : Track does not cross line!!\n");
	268
	269
	270	//Convert center of curvature to local coordinates:
	271	//Float_t xyz_coc[3] = {GetCenterX(),GetCenterY(),0};
	272	//fTransform->Global2Local(xyz_coc,slice);
	273
	274	//Take the dot product of the tangent vector of the track, and
	275	//vector perpendicular to the padrow.
	276
	277	Double_t tangent[2];
	278	//tangent[1] = (xyz[0] - xyz_coc[0])/GetRadius();
	279	//tangent[0] = -1.*(xyz[1] - xyz_coc[1])/GetRadius();
	280	tangent[1] = (xyz[0] - GetCenterX())/GetRadius();
	281	tangent[0] = -1.*(xyz[1] - GetCenterY())/GetRadius();
	282
	283	Double_t perp_padrow[2] = {1,0}; //locally in slice
	284
	285	Double_t cos_beta = fabs(tangent[0]perp_padrow[0] + tangent[1]perp_padrow[1]);
	286	return acos(cos_beta);
	287
	288	}
	289
	290	Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t padrow,Float_t *xyz)
	291	{
	292	//Assumes the track is given in local coordinates
	293
	294	if(!fIsHelix)
	295	{
	296	printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
	297	return 0;
	298	}
	299
	300
	301	if(!IsLocal())
	302	{
	303	printf("GetCrossingPoint: Track is given on global coordinates\n");
	304	return false;
	305	}
	306
	307	Double_t xHit = fTransform->Row2X(padrow);
	308
	309	//xyz[0] = fTransform->Row2X(padrow);
310	xyz[0] = xHit;
311	Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
312	Double_t r2 = GetRadius()*GetRadius();
313	if(aa > r2)
314	return false;
315
316	Double_t aa2 = sqrt(r2 - aa);
317	Double_t y1 = GetCenterY() + aa2;
318	Double_t y2 = GetCenterY() - aa2;
319	xyz[1] = y1;
320	if(fabs(y2) < fabs(y1)) xyz[1] = y2;
321	xyz[2] = 0; //only consider transverse plane
322
323	return true;
324	}
325
326
327	Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t slice,Int_t padrow,Float_t *xyz)
328	{
329	//Calculate the crossing point in transverse plane of this track and given
330	//padrow (y = a*x + b). Point is given in local coordinates in slice.
331	//Assumes the track is given in global coordinates.
332
333	if(!fIsHelix)
334	{
335	printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
336	return 0;
337	}
338
339
340	if(IsLocal())
341	{
342	printf("GetCrossingPoint: Track is given in local coordintes!!!\n");
343	return false;
344	}
345
346	Double_t padrowradii = fTransform->Row2X(padrow);
347
348
349	Float_t rotation_angle = (slice20)ToRad;
350
351	Float_t cs,sn;
352	cs = cos(rotation_angle);
353	sn = sin(rotation_angle);
354
355	Double_t a = -1.*cs/sn;
356	Double_t b = padrowradii/sn;
357
358	Double_t ycPrime = GetCenterY() - b ;
359	Double_t aa = ( 1. + a * a ) ;
360	Double_t bb = -2. * ( GetCenterX() + a * ycPrime ) ;
361	Double_t cc = ( GetCenterX() * GetCenterX() + ycPrime * ycPrime - GetRadius() * GetRadius() ) ;
362
363	Double_t racine = bb * bb - 4. * aa * cc ;
364	if ( racine < 0 ) return false ;
365	Double_t rootRacine = sqrt(racine) ;
366
367	Double_t oneOverA = 1./aa;
368	//
369	// First solution
370	//
371	Double_t x1 = 0.5 * oneOverA * ( -1. * bb + rootRacine ) ;
372	Double_t y1 = a * x1 + b ;
373	Double_t r1 = sqrt(x1x1+y1y1);
374	//
375	// Second solution
376	//
377	Double_t x2 = 0.5 * oneOverA * ( -1. * bb - rootRacine ) ;
378	Double_t y2 = a * x2 + b ;
379	Double_t r2 = sqrt(x2x2+y2y2);
380	//
381	// Choose close to (0,0)
382	//
383	Double_t xHit ;
384	Double_t yHit ;
385	if ( r1 < r2 ) {
386	xHit = x1 ;
387	yHit = y1 ;
388	}
389	else {
390	xHit = x2 ;
391	yHit = y2 ;
392	}
393
394	xyz[0] = xHit;
395	xyz[1] = yHit;
396	xyz[2] = 0;
397
398	fTransform->Global2Local(xyz,slice);
399
400	return true;
401	}
402	*/