HLT/src/AliL3HoughTrack.cxx

   1
   2 #include <math.h>
   3
   4 #include "AliL3Transform.h"
   5 #include "AliL3Defs.h"
   6 #include "AliL3HoughTrack.h"
   7
   8 ClassImp(AliL3HoughTrack)
   9
  10
  11 AliL3HoughTrack::AliL3HoughTrack()
  12 {
  13   //Constructor
  14
  15   fWeight = 0;
  16   fMinDist=0;
  17   fTransform = new AliL3Transform();
  18   fDLine = 0;
  19   fPsiLine = 0;
  20   fIsHelix = true;
  21 }
  22
  23
  24 AliL3HoughTrack::~AliL3HoughTrack()
  25 {
  26   //Destructor
  27   if(fTransform)
  28     delete fTransform;
  29 }
  30
  31 void AliL3HoughTrack::Set(AliL3Track *track)
  32 {
  33
  34   AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
  35
  36   fWeight = tpt->GetWeight();
  37   fDLine = tpt->GetDLine();
  38   fPsiLine = tpt->GetPsiLine();
  39   SetNHits(tpt->GetWeight());
  40   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  41   fIsHelix = false;
  42
  43
  44 }
  45
  46 void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t phi,Int_t weight)
  47 {
  48
  49   fWeight = weight;
  50   fMinDist = 100000;
  51   SetKappa(kappa);
  52   SetPhi0(phi);
  53   Double_t pt = fabs(BFACT*bField/kappa);
  54   SetPt(pt);
  55   Double_t radius = 1/fabs(kappa);
  56   SetRadius(radius);
  57
  58   //set nhits for sorting.
  59   SetNHits(weight);
  60
  61   SetCharge(copysign(1,kappa));
  62   Double_t charge = -1.*kappa;
  63   Double_t trackPhi0 = GetPhi0() + charge*0.5*Pi/fabs(charge);
  64
  65   //The first point on track is origo:
  66   Double_t x0=0;
  67   Double_t y0=0;
  68
  69   Double_t xc = x0 - GetRadius() * cos(trackPhi0) ;
  70   Double_t yc = y0 - GetRadius() * sin(trackPhi0) ;
  71   SetCenterX(xc);
  72   SetCenterY(yc);
  73   fIsHelix = true;
  74 }
  75
  76 void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t ref_row)
  77 {
  78   //Initialize a track piece, not yet a track
  79   //Used in case of straight line transformation
  80
  81   //Transform line parameters to coordinate system of slice:
  82
  83   D = D + fTransform->Row2X(ref_row)*cos(psi);
  84
  85   fDLine = D;
  86   fPsiLine = psi;
  87   fWeight = weight;
  88   SetNHits(weight);
  89   SetRowRange(rowrange[0],rowrange[1]);
  90   fIsHelix = false;
  91
  92 }
  93
  94 void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t min_dist)
  95 {
  96
  97   if(min_dist < fMinDist)
  98     {
  99       fMinDist = min_dist;
 100       SetMCid(mcid);
 101     }
 102
 103 }
 104
 105 void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Double_t *xy)
 106 {
 107
 108
 109   if(fIsHelix)
 110     {
 111       printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
 112       return;
 113     }
 114
 115   Double_t xhit = fTransform->Row2X(padrow);
 116   Double_t a = -1/tan(fPsiLine);
 117   Double_t b = fDLine/sin(fPsiLine);
 118
 119   Double_t yhit = a*xhit + b;
 120   xy[0] = xhit;
 121   xy[1] = yhit;
 122
 123 }
 124
 125 /*
 126 Double_t AliL3HoughTrack::GetCrossingAngle(Int_t padrow)
 127 {
 128   //Calculate the crossing angle between track and given padrow.
 129
 130   if(!fIsHelix)
 131     {
 132       printf("AliL3HoughTrack::GetCrossingAngle : Track is not a helix\n");
 133       return 0;
 134     }
 135
 136   if(!IsLocal())
 137     {
 138       printf("Track is not given in local coordinates\n");
 139       return 0;
 140     }
 141
 142   Float_t xyz[3];
 143   if(!GetCrossingPoint(padrow,xyz))
 144     printf("AliL3HoughTrack::GetCrossingPoint : Track does not cross line!!\n");
 145
 146
 147   //Convert center of curvature to local coordinates:
 148   //Float_t xyz_coc[3] = {GetCenterX(),GetCenterY(),0};
 149   //fTransform->Global2Local(xyz_coc,slice);
 150
 151   //Take the dot product of the tangent vector of the track, and
 152   //vector perpendicular to the padrow.
 153
 154   Double_t tangent[2];
 155   //tangent[1] = (xyz[0] - xyz_coc[0])/GetRadius();
 156   //tangent[0] = -1.*(xyz[1] - xyz_coc[1])/GetRadius();
 157   tangent[1] = (xyz[0] - GetCenterX())/GetRadius();
 158   tangent[0] = -1.*(xyz[1] - GetCenterY())/GetRadius();
 159
 160   Double_t perp_padrow[2] = {1,0}; //locally in slice
 161
 162   Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
 163   return acos(cos_beta);
 164
 165 }
 166
 167 Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t padrow,Float_t *xyz)
 168 {
 169   //Assumes the track is given in local coordinates
 170
 171   if(!fIsHelix)
 172     {
 173       printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
 174       return 0;
 175     }
 176
 177
 178   if(!IsLocal())
 179     {
 180       printf("GetCrossingPoint: Track is given on global coordinates\n");
 181       return false;
 182     }
 183
 184   Double_t xHit = fTransform->Row2X(padrow);
 185
 186   //xyz[0] = fTransform->Row2X(padrow);
 187   xyz[0] = xHit;
 188   Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
 189   Double_t r2 = GetRadius()*GetRadius();
 190   if(aa > r2)
 191     return false;
 192
 193   Double_t aa2 = sqrt(r2 - aa);
 194   Double_t y1 = GetCenterY() + aa2;
 195   Double_t y2 = GetCenterY() - aa2;
 196   xyz[1] = y1;
 197   if(fabs(y2) < fabs(y1)) xyz[1] = y2;
 198   xyz[2] = 0; //only consider transverse plane
 199
 200   return true;
 201 }
 202
 203
 204 Bool_t AliL3HoughTrack::GetCrossingPoint(Int_t slice,Int_t padrow,Float_t *xyz)
 205 {
 206   //Calculate the crossing point in transverse plane of this track and given
 207   //padrow (y = a*x + b). Point is given in local coordinates in slice.
 208   //Assumes the track is given in global coordinates.
 209
 210   if(!fIsHelix)
 211     {
 212       printf("AliL3HoughTrack::GetCrossingPoint : Track is not a helix\n");
 213       return 0;
 214     }
 215
 216
 217   if(IsLocal())
 218     {
 219       printf("GetCrossingPoint: Track is given in local coordintes!!!\n");
 220       return false;
 221     }
 222
 223   Double_t padrowradii = fTransform->Row2X(padrow);
 224
 225
 226   Float_t rotation_angle = (slice*20)*ToRad;
 227
 228   Float_t cs,sn;
 229   cs = cos(rotation_angle);
 230   sn = sin(rotation_angle);
 231
 232   Double_t a = -1.*cs/sn;
 233   Double_t b = padrowradii/sn;
 234
 235   Double_t ycPrime = GetCenterY() - b ;
 236   Double_t aa = ( 1. + a * a ) ;
 237   Double_t bb = -2. * ( GetCenterX() + a * ycPrime ) ;
 238   Double_t cc = ( GetCenterX() * GetCenterX() + ycPrime * ycPrime - GetRadius() * GetRadius() ) ;
 239
 240   Double_t racine = bb * bb - 4. * aa * cc ;
 241   if ( racine < 0 ) return false ;
 242   Double_t rootRacine = sqrt(racine) ;
 243
 244   Double_t oneOverA = 1./aa;
 245 //
 246 //   First solution
 247 //
 248    Double_t x1 = 0.5 * oneOverA * ( -1. * bb + rootRacine ) ;
 249    Double_t y1 = a * x1 + b ;
 250    Double_t r1 = sqrt(x1*x1+y1*y1);
 251 //
 252 //   Second solution
 253 //
 254    Double_t x2 = 0.5 * oneOverA * ( -1. * bb - rootRacine ) ;
 255    Double_t y2 = a * x2 + b ;
 256    Double_t r2 = sqrt(x2*x2+y2*y2);
 257 //
 258 //    Choose close to (0,0)
 259 //
 260    Double_t xHit ;
 261    Double_t yHit ;
 262    if ( r1 < r2 ) {
 263       xHit = x1 ;
 264       yHit = y1 ;
 265    }
 266    else {
 267       xHit = x2 ;
 268       yHit = y2 ;
 269    }
 270
 271    xyz[0] = xHit;
 272    xyz[1] = yHit;
 273    xyz[2] = 0;
 274
 275    fTransform->Global2Local(xyz,slice);
 276
 277    return true;
 278 }
 279 */