L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / hough / AliHLTHoughTrack.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliHLTStandardIncludes.h"
7
8 #include "AliHLTLogging.h"
9 #include "AliHLTTrack.h"
10 #include "AliHLTHoughTrack.h"
11 #include "AliHLTTransform.h"
12 #include "AliHLTHoughTransformerRow.h"
13
14 #if __GNUC__ >= 3
15 using namespace std;
16 #endif
17
18 /** \class AliHLTHoughTrack
19 <pre>
20 //_____________________________________________________________
21 // AliHLTHoughTrack
22 //
23 // Track class for Hough tracklets
24 //
25 </pre>
26 */
27
28 ClassImp(AliHLTHoughTrack)
29
30
31     AliHLTHoughTrack::AliHLTHoughTrack() : AliHLTTrack()
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 AliHLTHoughTrack::~AliHLTHoughTrack()
46 {
47   //dtor
48 }
49
50 void AliHLTHoughTrack::Set(AliHLTTrack *track)
51 {
52   //Basically copy constructor
53   AliHLTHoughTrack *tpt = (AliHLTHoughTrack*)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 AliHLTHoughTrack::Compare(const AliHLTTrack *tpt) const
84 {
85   //Compare 2 hough tracks according to their weight
86   AliHLTHoughTrack *track = (AliHLTHoughTrack*)tpt;
87   if(track->GetWeight() < GetWeight()) return 1;
88   if(track->GetWeight() > GetWeight()) return -1;
89   return 0;
90 }
91
92 void AliHLTHoughTrack::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 = AliHLTTransform::PiHalf() - theta;
99   Double_t tgl = tan(dipangle);
100   SetTgl(tgl);
101 }
102
103 void AliHLTHoughTrack::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(AliHLTLog::kWarning,"AliHLTHoughTrack::UpdateToFirstRow()","Track parameters")
112       <<AliHLTLog::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()) / AliHLTTransform::GetBFieldValue();
127   Double_t tPhi0 = GetPsi() + GetCharge() * AliHLTTransform::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(AliHLTLog::kError,"AliHLTHoughTrack::UpdateToFirstRow","Tracks")<<AliHLTLog::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 + AliHLTTransform::TwoPi();
148   Double_t deltat = fmod((-GetCharge()*td + GetCharge()*tPhi0),AliHLTTransform::TwoPi());
149   if ( deltat < 0. ) deltat += AliHLTTransform::TwoPi();
150   else if ( deltat > AliHLTTransform::TwoPi() ) deltat -= AliHLTTransform::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() * AliHLTTransform::PiHalf() / abs(GetCharge()) ;
160   if ( tPsi > AliHLTTransform::TwoPi() ) tPsi -= AliHLTTransform::TwoPi() ;
161   else if ( tPsi < 0. ) tPsi += AliHLTTransform::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 AliHLTHoughTrack::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(AliHLTTransform::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*AliHLTTransform::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 AliHLTHoughTrack::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)/(AliHLTHoughTransformerRow::GetBeta1()-AliHLTHoughTransformerRow::GetBeta2()));
209   Double_t kappa = 2.0*(alpha1*cos(psi)-AliHLTHoughTransformerRow::GetBeta1()*sin(psi));
210   SetTrackParameters(kappa,psi,weight);
211
212   Double_t zovr;
213   Double_t etaparam1 = AliHLTHoughTransformerRow::GetEtaCalcParam1();
214   Double_t etaparam2 = AliHLTHoughTransformerRow::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 AliHLTHoughTrack::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 AliHLTHoughTrack::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 AliHLTHoughTrack::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("AliHLTHoughTrack::GetLineCrossingPoint : Track is not a line\n");
257       return;
258     }
259
260   Float_t xhit = AliHLTTransform::Row2X(padrow) - AliHLTTransform::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 }