HLT/hough/AliL3HoughTrack.cxx

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