HLT/hough/AliL3HoughTrack.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Logging.h"
   9 #include "AliL3Track.h"
  10 #include "AliL3HoughTrack.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3HoughTransformerRow.h"
  13
  14 #if __GNUC__ >= 3
  15 using namespace std;
  16 #endif
  17
  18 /** \class AliL3HoughTrack
  19 <pre>
  20 //_____________________________________________________________
  21 // AliL3HoughTrack
  22 //
  23 // Track class for Hough tracklets
  24 //
  25 </pre>
  26 */
  27
  28 ClassImp(AliL3HoughTrack)
  29
  30
  31     AliL3HoughTrack::AliL3HoughTrack() : AliL3Track()
  32 {
  33   //Constructor
  34
  35   fWeight = 0;
  36   fMinDist=0;
  37   fDLine = 0;
  38   fPsiLine = 0;
  39   fIsHelix = true;
  40   fEtaIndex = -1;
  41   fEta = 0;
  42   ComesFromMainVertex(kTRUE);
  43 }
  44
  45 AliL3HoughTrack::~AliL3HoughTrack()
  46 {
  47   //dtor
  48 }
  49
  50 void AliL3HoughTrack::Set(AliL3Track *track)
  51 {
  52   //Basically copy constructor
  53   AliL3HoughTrack *tpt = (AliL3HoughTrack*)track;
  54   SetTrackParameters(tpt->GetKappa(),tpt->GetPsi(),tpt->GetWeight());
  55   SetEtaIndex(tpt->GetEtaIndex());
  56   SetEta(tpt->GetEta());
  57   SetTgl(tpt->GetTgl());
  58   SetPsi(tpt->GetPsi());
  59   SetPterr(tpt->GetPterr());
  60   SetTglerr(tpt->GetTglerr());
  61   SetPsierr(tpt->GetPsierr());
  62   SetCenterX(tpt->GetCenterX());
  63   SetCenterY(tpt->GetCenterY());
  64   SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
  65   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
  66   SetCharge(tpt->GetCharge());
  67   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  68   SetSlice(tpt->GetSlice());
  69   SetHits(tpt->GetNHits(),(UInt_t *)tpt->GetHitNumbers());
  70   SetMCid(tpt->GetMCid());
  71   SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
  72   SetSector(tpt->GetSector());
  73   return;
  74
  75 //    fWeight = tpt->GetWeight();
  76 //    fDLine = tpt->GetDLine();
  77 //    fPsiLine = tpt->GetPsiLine();
  78 //    SetNHits(tpt->GetWeight());
  79 //    SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  80 //    fIsHelix = false;
  81 }
  82
  83 Int_t AliL3HoughTrack::Compare(const AliL3Track *tpt) const
  84 {
  85   //Compare 2 hough tracks according to their weight
  86   AliL3HoughTrack *track = (AliL3HoughTrack*)tpt;
  87   if(track->GetWeight() < GetWeight()) return 1;
  88   if(track->GetWeight() > GetWeight()) return -1;
  89   return 0;
  90 }
  91
  92 void AliL3HoughTrack::SetEta(Double_t f)
  93 {
  94   //Set eta, and calculate fTanl, which is the tan of dipangle
  95
  96   fEta = f;
  97   Double_t theta = 2*atan(exp(-1.*fEta));
  98   Double_t dipangle = AliL3Transform::PiHalf() - theta;
  99   Double_t tgl = tan(dipangle);
 100   SetTgl(tgl);
 101 }
 102
 103 void AliL3HoughTrack::UpdateToFirstRow()
 104 {
 105   //Update the track parameters to the point where track cross
 106   //its first padrow.`
 107
 108   //Get the crossing point with the first padrow:
 109   Float_t xyz[3];
 110   if(!GetCrossingPoint(GetFirstRow(),xyz))
 111     LOG(AliL3Log::kWarning,"AliL3HoughTrack::UpdateToFirstRow()","Track parameters")
 112       <<AliL3Log::kDec<<"Track does not cross padrow "<<GetFirstRow()<<" centerx "
 113       <<GetCenterX()<<" centery "<<GetCenterY()<<" Radius "<<GetRadius()<<" tgl "<<GetTgl()<<ENDLOG;
 114
 115   //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());
 116   //printf("Before: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
 117
 118   Double_t radius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 119
 120   //Get the track parameters
 121
 122   /*
 123     Double_t x0    = GetR0() * cos(GetPhi0()) ;
 124     Double_t y0    = GetR0() * sin(GetPhi0()) ;
 125   */
 126   Double_t rc    = GetRadius();//fabs(GetPt()) / AliL3Transform::GetBFieldValue();
 127   Double_t tPhi0 = GetPsi() + GetCharge() * AliL3Transform::PiHalf() / abs(GetCharge()) ;
 128   Double_t xc    = GetCenterX();//x0 - rc * cos(tPhi0) ;
 129   Double_t yc    = GetCenterY();//y0 - rc * sin(tPhi0) ;
 130
 131   //Check helix and cylinder intersect
 132   Double_t fac1 = xc*xc + yc*yc ;
 133   Double_t sfac = sqrt( fac1 ) ;
 134
 135   if ( fabs(sfac-rc) > radius || fabs(sfac+rc) < radius ) {
 136     LOG(AliL3Log::kError,"AliL3HoughTrack::UpdateToFirstRow","Tracks")<<AliL3Log::kDec<<
 137       "Track does not intersect"<<ENDLOG;
 138     return;
 139   }
 140
 141   //Find intersection
 142   Double_t fac2   = (radius*radius + fac1 - rc*rc) / (2.00 * radius * sfac ) ;
 143   Double_t phi    = atan2(yc,xc) + GetCharge()*acos(fac2) ;
 144   Double_t td     = atan2(radius*sin(phi) - yc,radius*cos(phi) - xc) ;
 145
 146   //Intersection in z
 147   if ( td < 0 ) td = td + AliL3Transform::TwoPi();
 148   Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),AliL3Transform::TwoPi());
 149   if ( deltat < 0. ) deltat += AliL3Transform::TwoPi();
 150   else if ( deltat > AliL3Transform::TwoPi() ) deltat -= AliL3Transform::TwoPi();
 151   Double_t z = GetZ0() + rc * GetTgl() * deltat ;
 152
 153   Double_t xExtra = radius * cos(phi) ;
 154   Double_t yExtra = radius * sin(phi) ;
 155
 156   Double_t tPhi = atan2(yExtra-yc,xExtra-xc);
 157
 158   //if ( tPhi < 0 ) tPhi += 2. * M_PI ;
 159   Double_t tPsi = tPhi - GetCharge() * AliL3Transform::PiHalf() / abs(GetCharge()) ;
 160   if ( tPsi > AliL3Transform::TwoPi() ) tPsi -= AliL3Transform::TwoPi() ;
 161   else if ( tPsi < 0. ) tPsi += AliL3Transform::TwoPi();
 162
 163   //And finally, update the track parameters
 164   SetR0(radius);
 165   SetPhi0(phi);
 166   SetZ0(z);
 167   SetPsi(tPsi);
 168   SetFirstPoint(xyz[0],xyz[1],z);
 169   //printf("After: first %f %f %f tgl %f center %f %f charge %d\n",GetFirstPointX(),GetFirstPointY(),GetFirstPointZ(),GetTgl(),GetCenterX(),GetCenterY(),GetCharge());
 170
 171   //printf("First point set %f %f %f\n",xyz[0],xyz[1],z);
 172
 173   //Also, set the coordinates of the point where track crosses last padrow:
 174   GetCrossingPoint(GetLastRow(),xyz);
 175   SetLastPoint(xyz[0],xyz[1],xyz[2]);
 176   //printf("last point %f %f %f\n",xyz[0],xyz[1],xyz[2]);
 177 }
 178
 179 void AliL3HoughTrack::SetTrackParameters(Double_t kappa,Double_t eangle,Int_t weight)
 180 {
 181   //Set track parameters - sort of ctor
 182   fWeight = weight;
 183   fMinDist = 100000;
 184   SetKappa(kappa);
 185   Double_t pt = fabs(AliL3Transform::GetBFieldValue()/kappa);
 186   SetPt(pt);
 187   Double_t radius = 1/fabs(kappa);
 188   SetRadius(radius);
 189   SetFirstPoint(0,0,0);
 190   SetPsi(eangle); //Psi = emission angle when first point is vertex
 191   SetPhi0(0);     //not defined for vertex reference point
 192   SetR0(0);
 193   Double_t charge = -1.*kappa;
 194   SetCharge((Int_t)copysign(1.,charge));
 195   Double_t trackPhi0 = GetPsi() + charge*0.5*AliL3Transform::Pi()/fabs(charge);
 196   Double_t xc = GetFirstPointX() - GetRadius() * cos(trackPhi0) ;
 197   Double_t yc = GetFirstPointY() - GetRadius() * sin(trackPhi0) ;
 198   SetCenterX(xc);
 199   SetCenterY(yc);
 200   SetNHits(1); //just for the trackarray IO
 201   fIsHelix = true;
 202 }
 203
 204 void AliL3HoughTrack::SetTrackParametersRow(Double_t alpha1,Double_t alpha2,Double_t eta,Int_t weight)
 205 {
 206   //Set track parameters for HoughTransformerRow
 207   //This includes curvature,emission angle and eta
 208   Double_t psi = atan((alpha1-alpha2)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
 209   Double_t kappa = 2.0*(alpha1*cos(psi)-AliL3HoughTransformerRow::GetBeta1()*sin(psi));
 210   SetTrackParameters(kappa,psi,weight);
 211
 212   Double_t zovr;
 213   Double_t etaparam1 = AliL3HoughTransformerRow::GetEtaCalcParam1();
 214   Double_t etaparam2 = AliL3HoughTransformerRow::GetEtaCalcParam2();
 215   if(eta>0)
 216     zovr = (etaparam1 - sqrt(etaparam1*etaparam1 - 4.*etaparam2*eta))/(2.*etaparam2);
 217   else
 218     zovr = -1.*(etaparam1 - sqrt(etaparam1*etaparam1 + 4.*etaparam2*eta))/(2.*etaparam2);
 219   Double_t r = sqrt(1.+zovr*zovr);
 220   Double_t exacteta = 0.5*log((1+zovr/r)/(1-zovr/r));
 221   SetEta(exacteta);
 222 }
 223
 224 void AliL3HoughTrack::SetLineParameters(Double_t psi,Double_t D,Int_t weight,Int_t *rowrange,Int_t /*ref_row*/)
 225 {
 226   //Initialize a track piece, not yet a track
 227   //Used in case of straight line transformation
 228
 229   //Transform line parameters to coordinate system of slice:
 230
 231   //  D = D + fTransform->Row2X(ref_row)*cos(psi);
 232
 233   fDLine = D;
 234   fPsiLine = psi;
 235   fWeight = weight;
 236   SetNHits(1);
 237   SetRowRange(rowrange[0],rowrange[1]);
 238   fIsHelix = false;
 239 }
 240
 241 void AliL3HoughTrack::SetBestMCid(Int_t mcid,Double_t mindist)
 242 {
 243   //Finds and set the closest mc label
 244   if(mindist < fMinDist)
 245     {
 246       fMinDist = mindist;
 247       SetMCid(mcid);
 248     }
 249 }
 250
 251 void AliL3HoughTrack::GetLineCrossingPoint(Int_t padrow,Float_t *xy)
 252 {
 253   //Returns the crossing point of the track with a given padrow
 254   if(fIsHelix)
 255     {
 256       printf("AliL3HoughTrack::GetLineCrossingPoint : Track is not a line\n");
 257       return;
 258     }
 259
 260   Float_t xhit = AliL3Transform::Row2X(padrow) - AliL3Transform::Row2X(GetFirstRow());
 261   Float_t a = -1/tan(fPsiLine);
 262   Float_t b = fDLine/sin(fPsiLine);
 263   Float_t yhit = a*xhit + b;
 264   xy[0] = xhit;
 265   xy[1] = yhit;
 266 }