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 "AliL3Defs.h"
  10 #include "AliL3HoughTrack.h"
  11 #include "AliL3Transform.h"
  12
  13 //_____________________________________________________________
  14 // AliL3HoughTrack
  15 //
  16 // Track class for Hough tracklets
  17
  18 ClassImp(AliL3HoughTrack)
  19
  20
  21 AliL3HoughTrack::AliL3HoughTrack()
  22 {
  23   //Constructor
  24
  25   fWeight = 0;
  26   fMinDist=0;
  27   fDLine = 0;
  28   fPsiLine = 0;
  29   fIsHelix = true;
  30   fEtaIndex = -1;
  31 }
  32
  33
  34 AliL3HoughTrack::~AliL3HoughTrack()
  35 {
  36
  37 }
  38
  39 void AliL3HoughTrack::Set(AliL3Track *track)
  40 {
  41
  42   AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
  43   SetTrackParameters(tpt->GetKappa(),tpt->GetPhi0(),tpt->GetWeight());
  44   SetEtaIndex(tpt->GetEtaIndex());
  45   SetEta(tpt->GetEta());
  46   SetTgl(tpt->GetTgl());
  47   SetPsi(tpt->GetPsi());
  48   SetCenterX(tpt->GetCenterX());
  49   SetCenterY(tpt->GetCenterY());
  50   SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
  51   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
  52   SetCharge(tpt->GetCharge());
  53   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  54   SetSlice(tpt->GetSlice());
  55   SetNHits(1);
  56   return;
  57
  58   fWeight = tpt->GetWeight();
  59   fDLine = tpt->GetDLine();
  60   fPsiLine = tpt->GetPsiLine();
  61   SetNHits(tpt->GetWeight());
  62   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  63   fIsHelix = false;
  64
  65
  66 }
  67
  68 Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
  69 {
  70   AliL3HoughTrack *track = (AliL3HoughTrack*)tpt;
  71   if(track->GetWeight() < GetWeight()) return 1;
  72   if(track->GetWeight() > GetWeight()) return -1;
  73   return 0;
  74 }
  75
  76 void AliL3HoughTrack::SetEta(Double_t f)
  77 {
  78   //Set eta, and calculate fTanl, which is the tan of dipangle
  79
  80   fEta = f;
  81   Double_t theta = 2*atan(exp(-1.*fEta));
  82   Double_t dipangle = Pi/2 - theta;
  83   Double_t tgl = tan(dipangle);
  84   SetTgl(tgl);
  85 }
  86
  87 void AliL3HoughTrack::CalculateHelix()
  88 {
  89   return;
  90 }
  91
  92 void AliL3HoughTrack::UpdateToFirstRow()
  93 {
  94   //Update the track parameters to the point where track cross
  95   //its first padrow.
  96
  97   //Get the crossing point with the first padrow:
  98   Float_t xyz[3];
  99   if(!GetCrossingPoint(GetFirstRow(),xyz))
 100     LOG(AliL3Log::kWarning,"AliL3HoughTrack::UpdateToFirstRow()","Track parameters")
 101       <<AliL3Log::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
 102       <<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
 103
 104   //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());
 105   //printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
 106
 107   Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 108
 109   //Get the track parameters
 110
 111   /*
 112     Double_t x0    = GetR0() * cos(GetPhi0()) ;
 113     Double_t y0    = GetR0() * sin(GetPhi0()) ;
 114   */
 115   Double_t rc    = GetRadius();//fabs(GetPt()) / ( BFACT * AliL3Transform::GetBField() )  ;
 116   Double_t tPhi0 = GetPsi() + GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
 117   Double_t xc    = GetCenterX();//x0 - rc * cos(tPhi0) ;
 118   Double_t yc    = GetCenterY();//y0 - rc * sin(tPhi0) ;
 119
 120   //Check helix and cylinder intersect
 121   Double_t fac1 = xc*xc + yc*yc ;
 122   Double_t sfac = sqrt( fac1 ) ;
 123
 124   if ( fabs(sfac-rc) > radius || fabs(sfac+rc) < radius ) {
 125     LOG(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
 126       "Track does not intersect"<<ENDLOG;
 127     return;
 128   }
 129
 130   //Find intersection
 131   Double_t fac2   = ( radius*radius + fac1 - rc*rc) / (2.00 * radius * sfac ) ;
 132   Double_t phi    = atan2(yc,xc) + GetCharge()*acos(fac2) ;
 133   Double_t td     = atan2(radius*sin(phi) - yc,radius*cos(phi) - xc) ;
 134
 135   //Intersection in z
 136   if ( td < 0 ) td = td + 2. * pi ;
 137   Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),2*pi) ;
 138   if ( deltat < 0.      ) deltat += 2. * pi ;
 139   if ( deltat > 2.*pi ) deltat -= 2. * pi ;
 140   Double_t z = GetZ0() + rc * GetTgl() * deltat ;
 141
 142
 143   Double_t xExtra = radius * cos(phi) ;
 144   Double_t yExtra = radius * sin(phi) ;
 145
 146   Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
 147
 148   //if ( tPhi < 0 ) tPhi += 2. * M_PI ;
 149
 150   Double_t tPsi = tPhi - GetCharge() * 0.5 * pi / fabs(GetCharge()) ;
 151   if ( tPsi > 2. * pi ) tPsi -= 2. * pi ;
 152   if ( tPsi < 0.        ) tPsi += 2. * pi ;
 153
 154   //And finally, update the track parameters
 155   SetR0(radius);
 156   SetPhi0(phi);
 157   SetZ0(z);
 158   SetPsi(tPsi);
 159   SetFirstPoint(xyz[0],xyz[1],z);
 160   //printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
 161
 162   //printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
 163
 164   //Also, set the coordinates of the point where track crosses last padrow:
 165   GetCrossingPoint(GetLastRow(),xyz);
 166   SetLastPoint(xyz[0],xyz[1],xyz[2]);
 167   //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
 168 }
 169
 170 void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t phi,Int_t weight)
 171 {
 172
 173   fWeight = weight;
 174   fMinDist = 100000;
 175   SetKappa(kappa);
 176   SetPhi0(phi);
 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(phi); //Psi = Phi when first point is vertex
 183   SetR0(0);
 184   Double_t charge = -1.*kappa;
 185   SetCharge((Int_t)copysign(1.,charge));
 186
 187   Double_t trackPhi0 = GetPsi() + charge*0.5*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 }
 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
 205   /*
 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;
 214   */
 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
 231   /*
 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;
 245   */
 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 = (slice*20)*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(x1*x1+y1*y1);
 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(x2*x2+y2*y2);
 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 */