HLT/src/AliL3Track.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo$fi.uib.no>, Uli Frankenfeld <mailto:franken@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6
   7 #include "AliL3RootTypes.h"
   8
   9 #include "AliL3Logging.h"
  10 #include "AliL3Track.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3Vertex.h"
  13 #include <math.h>
  14
  15 //_____________________________________________________________
  16 // AliL3Track
  17 //
  18 // Track base class
  19 //Begin_Html
  20 /*
  21 <img src="track_coordinates.gif">
  22 */
  23 //End_Html
  24
  25 ClassImp(AliL3Track)
  26
  27 Float_t AliL3Track::BFACT = 0.0029980;
  28 Double_t AliL3Track::pi=3.14159265358979323846;
  29
  30 AliL3Track::AliL3Track()
  31 {
  32   //Constructor
  33
  34   fNHits = 0;
  35   fMCid = -1;
  36   fKappa=0;
  37   fRadius=0;
  38   fCenterX=0;
  39   fCenterY=0;
  40   ComesFromMainVertex(false);
  41   fQ = 0;
  42   fPhi0=0;
  43   fPsi=0;
  44   fR0=0;
  45   fTanl=0;
  46   fZ0=0;
  47   fPt=0;
  48   fLength=0;
  49   fIsLocal=true;
  50   fRowRange[0]=0;
  51   fRowRange[1]=0;
  52   memset(fHitNumbers,0,176*sizeof(UInt_t));
  53 }
  54
  55 void AliL3Track::Set(AliL3Track *tpt){
  56
  57   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  58   SetPhi0(tpt->GetPhi0());
  59   SetKappa(tpt->GetKappa());
  60   SetNHits(tpt->GetNHits());
  61   SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
  62   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
  63   SetPt(tpt->GetPt());
  64   SetPsi(tpt->GetPsi());
  65   SetTgl(tpt->GetTgl());
  66   SetCharge(tpt->GetCharge());
  67   SetHits(tpt->GetNHits(),(UInt_t *)tpt->GetHitNumbers());
  68
  69 }
  70
  71 Int_t AliL3Track::Compare(const AliL3Track *track) const
  72 {
  73   if(track->GetNHits() < GetNHits()) return 1;
  74   if(track->GetNHits() > GetNHits()) return -1;
  75   return 0;
  76 }
  77
  78 AliL3Track::~AliL3Track()
  79 {
  80   //Nothing to do
  81 }
  82
  83 Double_t AliL3Track::GetP() const
  84 {
  85   // Returns total momentum.
  86   return fabs(GetPt())*sqrt(1. + GetTgl()*GetTgl());
  87
  88 }
  89
  90 Double_t AliL3Track::GetPseudoRapidity() const
  91 {
  92   return 0.5 * log((GetP() + GetPz()) / (GetP() - GetPz()));
  93 }
  94 /*
  95 Double_t AliL3Track::GetEta() const
  96 {
  97   return GetPseudoRapidity();
  98 }
  99 */
 100 Double_t AliL3Track::GetRapidity() const
 101 {
 102   Double_t m_pi = 0.13957;
 103   return 0.5 * log((m_pi + GetPz()) / (m_pi - GetPz()));
 104 }
 105
 106 void AliL3Track::Rotate(Int_t slice,Bool_t tolocal)
 107 {
 108   //Rotate track to global parameters
 109   //If flag tolocal is set, the track is rotated
 110   //to local coordinates.
 111
 112
 113   Float_t psi[1] = {GetPsi()};
 114   if(!tolocal)
 115     AliL3Transform::Local2GlobalAngle(psi,slice);
 116   else
 117     AliL3Transform::Global2LocalAngle(psi,slice);
 118   SetPsi(psi[0]);
 119   Float_t first[3];
 120   first[0] = GetFirstPointX();
 121   first[1] = GetFirstPointY();
 122   first[2] = GetFirstPointZ();
 123   if(!tolocal)
 124     AliL3Transform::Local2Global(first,slice);
 125   else
 126     AliL3Transform::Global2Local(first,slice,kTRUE);
 127
 128   SetFirstPoint(first[0],first[1],first[2]);
 129   Float_t last[3];
 130   last[0] = GetLastPointX();
 131   last[1] = GetLastPointY();
 132   last[2] = GetLastPointZ();
 133   if(!tolocal)
 134     AliL3Transform::Local2Global(last,slice);
 135   else
 136     AliL3Transform::Global2Local(last,slice,kTRUE);
 137   SetLastPoint(last[0],last[1],last[2]);
 138
 139   Float_t center[3] = {GetCenterX(),GetCenterY(),0};
 140   if(!tolocal)
 141     AliL3Transform::Local2Global(center,slice);
 142   else
 143     AliL3Transform::Global2Local(center,slice,kTRUE);
 144   SetCenterX(center[0]);
 145   SetCenterY(center[1]);
 146
 147   if(!tolocal)
 148     fIsLocal=kFALSE;
 149   else
 150     fIsLocal=kTRUE;
 151 }
 152
 153 void AliL3Track::CalculateHelix(){
 154   //Calculate Radius, CenterX and Centery from Psi, X0, Y0
 155   //
 156
 157   fRadius = fPt / (BFACT*AliL3Transform::GetBField());
 158   if(fRadius) fKappa = -fQ*1./fRadius;
 159   else fRadius = 999999;  //just zero
 160   Double_t trackPhi0 = fPsi + fQ *0.5 * pi;
 161
 162   fCenterX = fFirstPoint[0] - fRadius *  cos(trackPhi0);
 163   fCenterY = fFirstPoint[1] - fRadius *  sin(trackPhi0);
 164 }
 165
 166 Double_t AliL3Track::GetCrossingAngle(Int_t padrow)
 167 {
 168   //Calculate the crossing angle between track and given padrow.
 169
 170   if(!IsLocal())
 171     {
 172       printf("Track is not given in local coordinates\n");
 173       return 0;
 174     }
 175
 176   Float_t xyz[3];
 177   if(!GetCrossingPoint(padrow,xyz))
 178     printf("AliL3HoughTrack::GetCrossingPoint : Track does not cross line!!\n");
 179
 180   //Take the dot product of the tangent vector of the track, and
 181   //vector perpendicular to the padrow.
 182
 183   Double_t tangent[2];
 184   tangent[1] = (xyz[0] - GetCenterX())/GetRadius();
 185   tangent[0] = -1.*(xyz[1] - GetCenterY())/GetRadius();
 186
 187   Double_t perp_padrow[2] = {1,0}; //locally in slice
 188
 189   Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
 190   return acos(cos_beta);
 191
 192 }
 193
 194 Bool_t AliL3Track::GetCrossingPoint(Int_t padrow,Float_t *xyz)
 195 {
 196   //Assumes the track is given in local coordinates
 197
 198   if(!IsLocal())
 199     {
 200       printf("GetCrossingPoint: Track is given on global coordinates\n");
 201       return false;
 202     }
 203
 204   Double_t xHit = AliL3Transform::Row2X(padrow);
 205
 206   xyz[0] = xHit;
 207   Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
 208   Double_t r2 = GetRadius()*GetRadius();
 209   if(aa > r2)
 210     return false;
 211
 212   Double_t aa2 = sqrt(r2 - aa);
 213   Double_t y1 = GetCenterY() + aa2;
 214   Double_t y2 = GetCenterY() - aa2;
 215   xyz[1] = y1;
 216   if(fabs(y2) < fabs(y1)) xyz[1] = y2;
 217
 218   Double_t yHit = xyz[1];
 219   Double_t angle1 = atan2((yHit - GetCenterY()),(xHit - GetCenterX()));
 220   if(angle1 < 0) angle1 += 2.*AliL3Transform::Pi();
 221   Double_t angle2 = atan2((GetFirstPointY() - GetCenterY()),(GetFirstPointX() - GetCenterX()));
 222   if(angle2 < 0) angle2 += 2.*AliL3Transform::Pi();
 223   Double_t diff_angle = angle1 - angle2;
 224   diff_angle = fmod(diff_angle,2*AliL3Transform::Pi());
 225   if((GetCharge()*diff_angle) > 0) diff_angle = diff_angle - GetCharge()*2.*AliL3Transform::Pi();
 226   Double_t s_tot = fabs(diff_angle)*GetRadius();
 227   Double_t zHit = GetFirstPointZ() + s_tot*GetTgl();
 228   xyz[2] = zHit;
 229
 230   return true;
 231 }
 232
 233
 234 Bool_t AliL3Track::CalculateReferencePoint(Double_t angle,Double_t radius){
 235   // Global coordinate: crossing point with y = ax+ b; a=tan(angle-AliL3Transform::Pi()/2);
 236   //
 237   const Double_t rr=radius;//132; //position of referece plane
 238   const Double_t xr = cos(angle) *rr;
 239   const Double_t yr = sin(angle) *rr;
 240
 241   Double_t a = tan(angle-pi/2);
 242   Double_t b = yr - a * xr;
 243
 244   Double_t pp=(fCenterX+a*fCenterY-a*b)/(1+pow(a,2));
 245   Double_t qq=(pow(fCenterX,2)+pow(fCenterY,2)-2*fCenterY*b+pow(b,2)-pow(fRadius,2))/(1+pow(a,2));
 246
 247   Double_t racine = pp*pp-qq;
 248   if(racine<0) return IsPoint(kFALSE);      //no Point
 249
 250   Double_t rootRacine = sqrt(racine);
 251   Double_t x0 = pp+rootRacine;
 252   Double_t x1 = pp-rootRacine;
 253   Double_t y0 = a*x0 + b;
 254   Double_t y1 = a*x1 + b;
 255
 256   Double_t diff0 = sqrt(pow(x0-xr,2)+pow(y0-yr,2));
 257   Double_t diff1 = sqrt(pow(x1-xr,2)+pow(y1-yr,2));
 258
 259   if(diff0<diff1){
 260     fPoint[0]=x0;
 261     fPoint[1]=y0;
 262   }
 263   else{
 264     fPoint[0]=x1;
 265     fPoint[1]=y1;
 266   }
 267
 268   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 269   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 270   if(fabs(trackPhi0-pointPhi0)>pi){
 271     if(trackPhi0<pointPhi0) trackPhi0 += 2*pi;
 272     else                    pointPhi0 += 2*pi;
 273   }
 274   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 275   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 276
 277   fPointPsi = pointPhi0 - fQ * 0.5 * pi;
 278   if(fPointPsi<0.)  fPointPsi+= 2*pi;
 279   fPointPsi = fmod(fPointPsi, 2*pi);
 280
 281   return IsPoint(kTRUE);
 282 }
 283
 284 Bool_t AliL3Track::CalculateEdgePoint(Double_t angle){
 285   // Global coordinate: crossing point with y = ax; a=tan(angle);
 286   //
 287   Double_t rmin=80;  //min Radius of TPC
 288   Double_t rmax=260; //max Radius of TPC
 289
 290   Double_t a = tan(angle);
 291   Double_t pp=(fCenterX+a*fCenterY)/(1+pow(a,2));
 292   Double_t qq=(pow(fCenterX,2)+pow(fCenterY,2)-pow(fRadius,2))/(1+pow(a,2));
 293   Double_t racine = pp*pp-qq;
 294   if(racine<0) return IsPoint(kFALSE);     //no Point
 295   Double_t rootRacine = sqrt(racine);
 296   Double_t x0 = pp+rootRacine;
 297   Double_t x1 = pp-rootRacine;
 298   Double_t y0 = a*x0;
 299   Double_t y1 = a*x1;
 300
 301   Double_t r0 = sqrt(pow(x0,2)+pow(y0,2));
 302   Double_t r1 = sqrt(pow(x1,2)+pow(y1,2));
 303   //find the right crossing point:
 304   //inside the TPC modules
 305   Bool_t ok0 = kFALSE;
 306   Bool_t ok1 = kFALSE;
 307
 308   if(r0>rmin&&r0<rmax){
 309     Double_t da=atan2(y0,x0);
 310     if(da<0) da+=2*pi;
 311     if(fabs(da-angle)<0.5)
 312       ok0 = kTRUE;
 313   }
 314   if(r1>rmin&&r1<rmax){
 315     Double_t da=atan2(y1,x1);
 316     if(da<0) da+=2*pi;
 317     if(fabs(da-angle)<0.5)
 318       ok1 = kTRUE;
 319   }
 320   if(!(ok0||ok1)) return IsPoint(kFALSE);   //no Point
 321
 322   if(ok0&&ok1){
 323     Double_t diff0 = sqrt(pow(fFirstPoint[0]-x0,2)+pow(fFirstPoint[1]-y0,2));
 324     Double_t diff1 = sqrt(pow(fFirstPoint[0]-x1,2)+pow(fFirstPoint[1]-y1,2));
 325     if(diff0<diff1) ok1 = kFALSE; //use ok0
 326     else ok0 = kFALSE;            //use ok1
 327   }
 328   if(ok0){fPoint[0]=x0; fPoint[1]=y0;}
 329   else   {fPoint[0]=x1; fPoint[1]=y1;}
 330
 331   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 332   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 333   if(fabs(trackPhi0-pointPhi0)>pi){
 334     if(trackPhi0<pointPhi0) trackPhi0 += 2*pi;
 335     else                    pointPhi0 += 2*pi;
 336   }
 337   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 338   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 339
 340   fPointPsi = pointPhi0 - fQ * 0.5 * pi;
 341   if(fPointPsi<0.)  fPointPsi+= 2*pi;
 342   fPointPsi = fmod(fPointPsi, 2*pi);
 343
 344   return IsPoint(kTRUE);
 345 }
 346
 347 Bool_t AliL3Track::CalculatePoint(Double_t xplane){
 348   // Local coordinate: crossing point with x plane
 349   //
 350   Double_t racine = pow(fRadius,2)-pow(xplane-fCenterX,2);
 351   if(racine<0) return IsPoint(kFALSE);
 352   Double_t rootRacine = sqrt(racine);
 353
 354   Double_t y0 = fCenterY + rootRacine;
 355   Double_t y1 = fCenterY - rootRacine;
 356   //Double_t diff0 = sqrt(pow(fFirstPoint[0]-xplane)+pow(fFirstPoint[1]-y0));
 357   //Double_t diff1 = sqrt(pow(fFirstPoint[0]-xplane)+pow(fFirstPoint[1]-y1));
 358   Double_t diff0 = fabs(y0-fFirstPoint[1]);
 359   Double_t diff1 = fabs(y1-fFirstPoint[1]);
 360
 361   fPoint[0]=xplane;
 362   if(diff0<diff1) fPoint[1]=y0;
 363   else            fPoint[1]=y1;
 364
 365   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 366   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 367   if(fabs(trackPhi0-pointPhi0)>pi){
 368     if(trackPhi0<pointPhi0) trackPhi0 += 2*pi;
 369     else                    pointPhi0 += 2*pi;
 370   }
 371   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 372   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 373
 374   fPointPsi = pointPhi0 - fQ * 0.5 * pi;
 375   if(fPointPsi<0.)  fPointPsi+= 2*pi;
 376   fPointPsi = fmod(fPointPsi, 2*pi);
 377
 378   return IsPoint(kTRUE);
 379 }
 380
 381 void AliL3Track::GetClosestPoint(AliL3Vertex *vertex,Double_t &closest_x,Double_t &closest_y,Double_t &closest_z)
 382 {
 383   //Calculate the point of closest approach to the vertex
 384
 385
 386   Double_t xc = GetCenterX() - vertex->GetX();//Shift the center of curvature with respect to the vertex
 387   Double_t yc = GetCenterY() - vertex->GetY();
 388
 389   Double_t dist_x1 = xc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 390   Double_t dist_y1 = yc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 391   Double_t distance1 = sqrt(dist_x1*dist_x1 + dist_y1*dist_y1);
 392
 393   Double_t dist_x2 = xc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 394   Double_t dist_y2 = yc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 395   Double_t distance2 = sqrt(dist_x2*dist_x2 + dist_y2*dist_y2);
 396
 397   //Choose the closest:
 398   if(distance1 < distance2)
 399     {
 400       closest_x = dist_x1 + vertex->GetX();
 401       closest_y = dist_y1 + vertex->GetY();
 402     }
 403   else
 404     {
 405       closest_x = dist_x2 + vertex->GetX();
 406       closest_y = dist_y2 + vertex->GetY();
 407     }
 408
 409   //Get the z coordinate:
 410   Double_t angle1 = atan2((closest_y-GetCenterY()),(closest_x-GetCenterX()));
 411   if(angle1 < 0) angle1 = angle1 + 2*AliL3Transform::Pi();
 412
 413   Double_t angle2 = atan2((GetFirstPointY()-GetCenterY()),(GetFirstPointX()-GetCenterX()));
 414   if(angle2 < 0) angle2 = angle2 + 2*AliL3Transform::Pi();
 415
 416   Double_t diff_angle = angle1 - angle2;
 417   diff_angle = fmod(diff_angle,2*AliL3Transform::Pi());
 418
 419   if((GetCharge()*diff_angle) < 0) diff_angle = diff_angle + GetCharge()*2*AliL3Transform::Pi();
 420   Double_t s_tot = fabs(diff_angle)*GetRadius();
 421
 422   closest_z = GetFirstPointZ() - s_tot*GetTgl();
 423 }