HLT/hough/AliL3HoughTrack.cxx

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