StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AliHLTTPCCAMerger.h
1 //-*- Mode: C++ -*-
2 // ************************************************************************
3 // This file is property of and copyright by the ALICE HLT Project *
4 // ALICE Experiment at CERN, All rights reserved. *
5 // See cxx source for full Copyright notice *
6 // *
7 //*************************************************************************
8 
9 #ifndef ALIHLTTPCCAMERGER_H
10 #define ALIHLTTPCCAMERGER_H
11 
12 #include "AliHLTTPCCADef.h"
13 #include "AliHLTTPCCAParam.h"
14 #include "AliHLTTPCCATrackParam.h"
15 #include "AliHLTTPCCASliceTrack.h"
16 
17 #include <vector>
18 
19 #if !defined(HLTCA_GPUCODE)
20 #include <iostream>
21 #endif
22 
23 class AliHLTTPCCATrackParamVector;
24 class AliHLTTPCCASliceTrack;
25 class AliHLTTPCCASliceOutput;
26 class AliHLTTPCCAMergedTrack;
27 class AliHLTTPCCAMergerOutput;
28 
34 class AliHLTTPCCAClusterInfo;
35 class AliHLTTPCCAMerger
36 {
37  // public:
38  // class AliHLTTPCCAClusterInfo;
39  private:
40 // struct AliHLTTPCCATrackMemory;
41 // struct AliHLTTPCCAHitMemory;
42 
43  class AliHLTTPCCASliceTrackInfo;
44 
45  class AliHLTTPCCABorderTrack
46  {
47  public:
48 
49  int TrackID() const { return fTrackID; }
50  float b() const{ return fb; }
51  float bErr2() const{ return fbErr2; }
52  float p() const{ return fp; }
53  float pErr2() const{ return fpErr2; }
54  unsigned short InnerRow() const {return fInnerRow;}
55  unsigned short OuterRow() const {return fOuterRow;}
56 
57  void SetInnerRow(unsigned short v) {fInnerRow = v;}
58  void SetOuterRow(unsigned short v) {fOuterRow = v;}
59  void SetTrackID ( int v ) { fTrackID = v; }
60  void Setb (float v) { fb = v; }
61  void SetbErr2 (float v) { fbErr2 = v; }
62  void Setp (float v) { fp = v; }
63  void SetpErr2 (float v) { fpErr2 = v; }
64 
65  private:
66  float fb;
67  float fbErr2;
68  float fp;
69  float fpErr2;
70  int fTrackID; // track index
71  unsigned short fInnerRow;
72  unsigned short fOuterRow;
73  };
74 
75  public:
76 
77  AliHLTTPCCAMerger();
78  ~AliHLTTPCCAMerger();
79 
80  void SetSliceParam( const AliHLTTPCCAParam &v ) { fSliceParam = v; }
81 
82  void Clear();
83  void SetSliceData( int index, const AliHLTTPCCASliceOutput *SliceData );
84  void Reconstruct();
85  const AliHLTTPCCAMergerOutput * Output() const { return fOutput; }
86 
87 // bool FitTrack( AliHLTTPCCATrackParam &T, float &Alpha,
88 // AliHLTTPCCATrackParam t0, float Alpha0, int hits[], int &NHits, bool dir = 0 );
89 
90  void SetSlices ( int i, AliHLTTPCCATracker *sl );
91 
92  static void SetDoNotMergeBorders(int i = 0) {fgDoNotMergeBorders = i;}
93 
94  int NTimers() { return fNTimers; }
95  float Timer( int i ) { return fTimers[i]; };
96  private:
97 
98  void ConvertTrackParamToVector( AliHLTTPCCATrackParam t0[ushort_v::Size], AliHLTTPCCATrackParamVector &t, int &nTracksV);
99 
100  sfloat_m FitTrack( AliHLTTPCCATrackParamVector &t, sfloat_v &Alpha0V,
101  int hits[2000][ushort_v::Size], ushort_v &firstHits, ushort_v::Memory &NTrackHits,
102  int &nTracksV, sfloat_m active0 = sfloat_m(true), bool dir = 1 );
103 
104  AliHLTTPCCAMerger( const AliHLTTPCCAMerger& );
105  const AliHLTTPCCAMerger &operator=( const AliHLTTPCCAMerger& ) const;
106 
107  void InvertCholetsky(float a[15]);
108  void MultiplySS(float const C[15], float const V[15], float K[5][5]);
109  void MultiplyMS(float const C[5][5], float const V[15], float K[15]);
110  void MultiplySR(float const C[15], float const r_in[5], float r_out[5]);
111  void FilterTracks(float const r[5], float const C[15], float const m[5], float const V[15], float R[5], float W[15], float &chi2);
112 
113  void InvertCholetsky(sfloat_v a[15]);
114  void MultiplySS(sfloat_v const C[15], sfloat_v const V[15], sfloat_v K[5][5]);
115  void MultiplyMS(sfloat_v const C[5][5], sfloat_v const V[15], sfloat_v K[15]);
116  void MultiplySR(sfloat_v const C[15], sfloat_v const r_in[5], sfloat_v r_out[5]);
117  void FilterTracks(sfloat_v const r[5], sfloat_v const C[15], sfloat_v const m[5], sfloat_v const V[15],
118  sfloat_v R[5], sfloat_v W[15], sfloat_v &chi2, const sfloat_m &mask = sfloat_m( true ));
119 
120  static bool CompareInnerRow (const AliHLTTPCCABorderTrack &b1, const AliHLTTPCCABorderTrack &b2) {
121  return (b1.InnerRow() > b2.InnerRow()) || ( (b1.InnerRow() == b2.InnerRow()) && (b1.b() > b2.b()) ) ;
122  }
123  static bool CompareOuterRow (const AliHLTTPCCABorderTrack &b1, const AliHLTTPCCABorderTrack &b2) {
124  return (b1.OuterRow() < b2.OuterRow())/* || ( (b1.OuterRow() == b2.OuterRow()) && (b1.b() > b2.b()) )*/ ;
125  }
126 
127  void UnpackSlices();
128  void Merging(int number=0);
129 
130  void MakeBorderTracks(AliHLTTPCCABorderTrack B[], unsigned short &nB, unsigned char &iSlice );
131  void MergeBorderTracks( AliHLTTPCCABorderTrack B1[], int N1, int iSlice1, AliHLTTPCCABorderTrack B2[], int N2, int iSlice2, int number,
132  unsigned short FirstTrIR[], unsigned short LastTrIR[]);
133 
134  static const int fgkNSlices = AliHLTTPCCAParameters::NumberOfSlices; //* N slices
135  static int fgDoNotMergeBorders;
136  AliHLTTPCCAParam fSliceParam; //* slice parameters (geometry, calibr, etc.)
137  const AliHLTTPCCASliceOutput *fkSlices[fgkNSlices]; //* array of input slice tracks
138  AliHLTTPCCATracker *slices[fgkNSlices]; //* array of input slice tracks
139  int fMaxClusterInfos; //* booked size of fClusterInfos array
140  AliHLTTPCCAClusterInfo *fClusterInfos; //* information about track clusters
141 
142  int fMaxTrackInfos; //* booked size of fTrackInfos array
143  AliHLTTPCCASliceTrackInfo *fTrackInfos; //* additional information for slice tracks
144  int fSliceTrackInfoStart[fgkNSlices]; //* slice starting index in fTrackInfos array;
145  int fSliceNTrackInfos[fgkNSlices]; //* N of slice track infos in fTrackInfos array;
146 
147  AliHLTTPCCAMergerOutput *fOutput; //* array of output merged tracks
148 
149  static const int fNTimers = 8;
150  float fTimers[fNTimers];
151 
152 // AliHLTTPCCAHitMemory fHitMemory;
153 };
154 
155 #include "AliHLTTPCCATrackParamVector.h"
156 
157 inline void AliHLTTPCCAMerger::ConvertTrackParamToVector( AliHLTTPCCATrackParam t0[ushort_v::Size], AliHLTTPCCATrackParamVector &t, int &nTracksV)
158 {
159  sfloat_v tmpVec;
160  short_v tmpVecShort;
161  sfloat_v::Memory tmpFloat;
162  short_v::Memory tmpShort;
163 
164  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].X();
165  tmpVec.load( tmpFloat );
166  t.SetX(tmpVec);
167  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].SignCosPhi();
168  tmpVec.load( tmpFloat );
169  t.SetSignCosPhi(tmpVec);
170 
171  for(int iP=0; iP<5; iP++)
172  {
173  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Par()[iP];
174  tmpVec.load( tmpFloat );
175  t.SetPar(iP,tmpVec);
176  }
177  for(int iC=0; iC<15; iC++)
178  {
179  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Cov()[iC];
180  tmpVec.load( tmpFloat );
181  t.SetCov(iC,tmpVec);
182  }
183  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Chi2();
184  tmpVec.load( tmpFloat );
185  t.SetChi2(tmpVec);
186  for(int iV=0; iV < nTracksV; iV++) tmpShort[iV] = t0[iV].NDF();
187  tmpVecShort.load( tmpShort );
188  t.SetNDF(tmpVecShort);
189 }
190 
191 inline void AliHLTTPCCAMerger::InvertCholetsky(sfloat_v a[15])
192 {
193  sfloat_v d[5], uud, u[5][5];
194  for(int i=0; i<5; i++)
195  {
196  d[i]=0.f;
197  for(int j=0; j<5; j++)
198  u[i][j]=0.;
199  }
200 
201  for(int i=0; i<5; i++)
202  {
203  uud=0.;
204  for(int j=0; j<i; j++)
205  uud += u[j][i]*u[j][i]*d[j];
206  uud = a[i*(i+3)/2] - uud;
207  sfloat_m small_val = CAMath::Abs(uud)<1.e-12f;
208  uud(small_val) = 1.e-12f;
209  d[i] = uud/CAMath::Abs(uud);
210  u[i][i] = sqrt(CAMath::Abs(uud));
211 
212  for(int j=i+1; j<5; j++)
213  {
214  uud = 0.;
215  for(int k=0; k<i; k++)
216  uud += u[k][i]*u[k][j]*d[k];
217  uud = a[j*(j+1)/2+i] - uud;
218  u[i][j] = d[i]/u[i][i]*uud;
219  }
220  }
221 
222  sfloat_v u1[5];
223 
224  for(int i=0; i<5; i++)
225  {
226  u1[i] = u[i][i];
227  u[i][i] = 1.f/u[i][i];
228  }
229  for(int i=0; i<4; i++)
230  {
231  u[i][i+1] = - u[i][i+1]*u[i][i]*u[i+1][i+1];
232  }
233  for(int i=0; i<3; i++)
234  {
235  u[i][i+2] = u[i][i+1]*u1[i+1]*u[i+1][i+2]-u[i][i+2]*u[i][i]*u[i+2][i+2];
236  }
237  for(int i=0; i<2; i++)
238  {
239  u[i][i+3] = u[i][i+2]*u1[i+2]*u[i+2][i+3] - u[i][i+3]*u[i][i]*u[i+3][i+3];
240  u[i][i+3] -= u[i][i+1]*u1[i+1]*(u[i+1][i+2]*u1[i+2]*u[i+2][i+3] - u[i+1][i+3]);
241  }
242  u[0][4] = u[0][2]*u1[2]*u[2][4] - u[0][4]*u[0][0]*u[4][4];
243  u[0][4] += u[0][1]*u1[1]*(u[1][4] - u[1][3]*u1[3]*u[3][4] - u[1][2]*u1[2]*u[2][4]);
244  u[0][4] += u[3][4]*u1[3]*(u[0][3] - u1[2]*u[2][3]*(u[0][2] - u[0][1]*u1[1]*u[1][2]));
245 
246  for(int i=0; i<5; i++)
247  a[i+10] = u[i][4]*d[4]*u[4][4];
248  for(int i=0; i<4; i++)
249  a[i+6] = u[i][3]*u[3][3]*d[3] + u[i][4]*u[3][4]*d[4];
250  for(int i=0; i<3; i++)
251  a[i+3] = u[i][2]*u[2][2]*d[2] + u[i][3]*u[2][3]*d[3] + u[i][4]*u[2][4]*d[4];
252  for(int i=0; i<2; i++)
253  a[i+1] = u[i][1]*u[1][1]*d[1] + u[i][2]*u[1][2]*d[2] + u[i][3]*u[1][3]*d[3] + u[i][4]*u[1][4]*d[4];
254  a[0] = u[0][0]*u[0][0]*d[0] + u[0][1]*u[0][1]*d[1] + u[0][2]*u[0][2]*d[2] + u[0][3]*u[0][3]*d[3] + u[0][4]*u[0][4]*d[4];
255 }
256 
257 inline void AliHLTTPCCAMerger::MultiplySS(sfloat_v const C[15], sfloat_v const V[15], sfloat_v K[5][5])
258 {
259  K[0][0] = C[0]*V[ 0] + C[1]*V[ 1] + C[3]*V[ 3] + C[6]*V[ 6] + C[10]*V[10];
260  K[0][1] = C[0]*V[ 1] + C[1]*V[ 2] + C[3]*V[ 4] + C[6]*V[ 7] + C[10]*V[11];
261  K[0][2] = C[0]*V[ 3] + C[1]*V[ 4] + C[3]*V[ 5] + C[6]*V[ 8] + C[10]*V[12];
262  K[0][3] = C[0]*V[ 6] + C[1]*V[ 7] + C[3]*V[ 8] + C[6]*V[ 9] + C[10]*V[13];
263  K[0][4] = C[0]*V[10] + C[1]*V[11] + C[3]*V[12] + C[6]*V[13] + C[10]*V[14];
264 
265  K[1][0] = C[1]*V[ 0] + C[2]*V[ 1] + C[4]*V[ 3] + C[7]*V[ 6] + C[11]*V[10];
266  K[1][1] = C[1]*V[ 1] + C[2]*V[ 2] + C[4]*V[ 4] + C[7]*V[ 7] + C[11]*V[11];
267  K[1][2] = C[1]*V[ 3] + C[2]*V[ 4] + C[4]*V[ 5] + C[7]*V[ 8] + C[11]*V[12];
268  K[1][3] = C[1]*V[ 6] + C[2]*V[ 7] + C[4]*V[ 8] + C[7]*V[ 9] + C[11]*V[13];
269  K[1][4] = C[1]*V[10] + C[2]*V[11] + C[4]*V[12] + C[7]*V[13] + C[11]*V[14];
270 
271  K[2][0] = C[3]*V[ 0] + C[4]*V[ 1] + C[5]*V[ 3] + C[8]*V[ 6] + C[12]*V[10];
272  K[2][1] = C[3]*V[ 1] + C[4]*V[ 2] + C[5]*V[ 4] + C[8]*V[ 7] + C[12]*V[11];
273  K[2][2] = C[3]*V[ 3] + C[4]*V[ 4] + C[5]*V[ 5] + C[8]*V[ 8] + C[12]*V[12];
274  K[2][3] = C[3]*V[ 6] + C[4]*V[ 7] + C[5]*V[ 8] + C[8]*V[ 9] + C[12]*V[13];
275  K[2][4] = C[3]*V[10] + C[4]*V[11] + C[5]*V[12] + C[8]*V[13] + C[12]*V[14];
276 
277  K[3][0] = C[6]*V[ 0] + C[7]*V[ 1] + C[8]*V[ 3] + C[9]*V[ 6] + C[13]*V[10];
278  K[3][1] = C[6]*V[ 1] + C[7]*V[ 2] + C[8]*V[ 4] + C[9]*V[ 7] + C[13]*V[11];
279  K[3][2] = C[6]*V[ 3] + C[7]*V[ 4] + C[8]*V[ 5] + C[9]*V[ 8] + C[13]*V[12];
280  K[3][3] = C[6]*V[ 6] + C[7]*V[ 7] + C[8]*V[ 8] + C[9]*V[ 9] + C[13]*V[13];
281  K[3][4] = C[6]*V[10] + C[7]*V[11] + C[8]*V[12] + C[9]*V[13] + C[13]*V[14];
282 
283  K[4][0] = C[10]*V[ 0] + C[11]*V[ 1] + C[12]*V[ 3] + C[13]*V[ 6] + C[14]*V[10];
284  K[4][1] = C[10]*V[ 1] + C[11]*V[ 2] + C[12]*V[ 4] + C[13]*V[ 7] + C[14]*V[11];
285  K[4][2] = C[10]*V[ 3] + C[11]*V[ 4] + C[12]*V[ 5] + C[13]*V[ 8] + C[14]*V[12];
286  K[4][3] = C[10]*V[ 6] + C[11]*V[ 7] + C[12]*V[ 8] + C[13]*V[ 9] + C[14]*V[13];
287  K[4][4] = C[10]*V[10] + C[11]*V[11] + C[12]*V[12] + C[13]*V[13] + C[14]*V[14];
288 }
289 
290 inline void AliHLTTPCCAMerger::MultiplyMS(sfloat_v const C[5][5], sfloat_v const V[15], sfloat_v K[15])
291 {
292  K[0] = C[0][0]*V[0] + C[0][1]*V[1] + C[0][2]*V[3] + C[0][3]*V[6] + C[0][4]*V[10];
293 
294  K[1] = C[1][0]*V[0] + C[1][1]*V[1] + C[1][2]*V[3] + C[1][3]*V[6] + C[1][4]*V[10];
295  K[2] = C[1][0]*V[1] + C[1][1]*V[2] + C[1][2]*V[4] + C[1][3]*V[7] + C[1][4]*V[11];
296 
297  K[3] = C[2][0]*V[0] + C[2][1]*V[1] + C[2][2]*V[3] + C[2][3]*V[6] + C[2][4]*V[10];
298  K[4] = C[2][0]*V[1] + C[2][1]*V[2] + C[2][2]*V[4] + C[2][3]*V[7] + C[2][4]*V[11];
299  K[5] = C[2][0]*V[3] + C[2][1]*V[4] + C[2][2]*V[5] + C[2][3]*V[8] + C[2][4]*V[12];
300 
301  K[6] = C[3][0]*V[0] + C[3][1]*V[1] + C[3][2]*V[3] + C[3][3]*V[6] + C[3][4]*V[10];
302  K[7] = C[3][0]*V[1] + C[3][1]*V[2] + C[3][2]*V[4] + C[3][3]*V[7] + C[3][4]*V[11];
303  K[8] = C[3][0]*V[3] + C[3][1]*V[4] + C[3][2]*V[5] + C[3][3]*V[8] + C[3][4]*V[12];
304  K[9] = C[3][0]*V[6] + C[3][1]*V[7] + C[3][2]*V[8] + C[3][3]*V[9] + C[3][4]*V[13];
305 
306  K[10] = C[4][0]*V[ 0] + C[4][1]*V[ 1] + C[4][2]*V[ 3] + C[4][3]*V[ 6] + C[4][4]*V[10];
307  K[11] = C[4][0]*V[ 1] + C[4][1]*V[ 2] + C[4][2]*V[ 4] + C[4][3]*V[ 7] + C[4][4]*V[11];
308  K[12] = C[4][0]*V[ 3] + C[4][1]*V[ 4] + C[4][2]*V[ 5] + C[4][3]*V[ 8] + C[4][4]*V[12];
309  K[13] = C[4][0]*V[ 6] + C[4][1]*V[ 7] + C[4][2]*V[ 8] + C[4][3]*V[ 9] + C[4][4]*V[13];
310  K[14] = C[4][0]*V[10] + C[4][1]*V[11] + C[4][2]*V[12] + C[4][3]*V[13] + C[4][4]*V[14];
311 }
312 
313 inline void AliHLTTPCCAMerger::MultiplySR(sfloat_v const C[15], sfloat_v const r_in[5], sfloat_v r_out[5])
314 {
315  r_out[0] = r_in[0]*C[ 0] + r_in[1]*C[ 1] + r_in[2]*C[ 3] +r_in[3]*C[ 6] + r_in[4]*C[10];
316  r_out[1] = r_in[0]*C[ 1] + r_in[1]*C[ 2] + r_in[2]*C[ 4] +r_in[3]*C[ 7] + r_in[4]*C[11];
317  r_out[2] = r_in[0]*C[ 3] + r_in[1]*C[ 4] + r_in[2]*C[ 5] +r_in[3]*C[ 8] + r_in[4]*C[12];
318  r_out[3] = r_in[0]*C[ 6] + r_in[1]*C[ 7] + r_in[2]*C[ 8] +r_in[3]*C[ 9] + r_in[4]*C[13];
319  r_out[4] = r_in[0]*C[10] + r_in[1]*C[11] + r_in[2]*C[12] +r_in[3]*C[13] + r_in[4]*C[14];
320 }
321 
322 inline void AliHLTTPCCAMerger::FilterTracks(sfloat_v const r[5], sfloat_v const C[15],
323  sfloat_v const m[5], sfloat_v const V[15],
324  sfloat_v R[5], sfloat_v W[15], sfloat_v &chi2, const sfloat_m &mask)
325 {
326  sfloat_v S[15];
327  for(int i=0; i<15; i++)
328  {
329  W[i] = C[i];
330  S[i] = C[i] + V[i];
331  }
332  for(int i=0; i<5; i++)
333  R[i] = r[i];
334 
335  InvertCholetsky(S);
336 
337  sfloat_v K[5][5];
338  MultiplySS(C,S,K);
339  sfloat_v dzeta[5];
340  for(int i=0; i<5; i++) dzeta[i] = m[i] - r[i];
341  sfloat_v KC[15];
342  MultiplyMS(K,C,KC);
343  for(int i=0; i< 15; i++)
344  W[i](mask) -= KC[i];
345 
346  sfloat_v kd(Vc::Zero);
347  for(int i=0; i<5; i++)
348  {
349  kd = 0.f;
350  for(int j=0; j<5; j++)
351  kd += K[i][j]*dzeta[j];
352  R[i](mask) += kd;
353  }
354  sfloat_v S_dzeta[5];
355  MultiplySR(S, dzeta, S_dzeta);
356  chi2(mask) = dzeta[0]*S_dzeta[0] + dzeta[1]*S_dzeta[1] + dzeta[2]*S_dzeta[2] + dzeta[3]*S_dzeta[3] + dzeta[4]*S_dzeta[4];
357 }
358 
359 class AliHLTTPCCAClusterInfo
360 {
361 
362  public:
363 
364  unsigned char ISlice() const { return fISlice; }
365  unsigned char IRow() const { return fIRow; }
366  unsigned int IClu() const { return fIClu; }
367  float X() const { return fX; }
368  float Y() const { return fY; }
369  float Z() const { return fZ; }
370 
371  void SetISlice ( unsigned char v ) { fISlice = v; }
372  void SetIRow ( unsigned char v ) { fIRow = v; }
373  void SetIClu ( unsigned int v ) { fIClu = v; }
374  void SetX ( float v ) { fX = v; }
375  void SetY ( float v ) { fY = v; }
376  void SetZ ( float v ) { fZ = v; }
377 
378  public:
379 
380  unsigned char fISlice; // slice number
381  unsigned char fIRow; // row number
382  unsigned int fIClu; // cluster number
383  float fX; // x position (slice coord.system)
384  float fY; // y position (slice coord.system)
385  float fZ; // z position (slice coord.system)
386 };
387 
388 class TrackHitsCompare {
389  // typedef AliHLTTPCCAMerger::AliHLTTPCCAClusterInfo AliHLTTPCCAClusterInfo;
390  public:
391  TrackHitsCompare(AliHLTTPCCAClusterInfo *clusterInfos): fClusterInfos(clusterInfos) {};
392 
393  bool operator()(int i, int j) { return fClusterInfos[i].X() < fClusterInfos[j].X(); };
394  private:
395  const AliHLTTPCCAClusterInfo *fClusterInfos;
396 };
397 
398 #endif
Generated by   doxygen 1.8.5