CbmRoot
L1CAMergeClones.cxx
Go to the documentation of this file.
1 /*
2  *=====================================================
3  *
4  * CBM Level 1 Reconstruction
5  *
6  * Authors: M.Zyzak
7  *
8  * e-mail :
9  *
10  *=====================================================
11  *
12  * Merge Clones
13  *
14  */
15 
16 #include "L1Algo.h"
17 #include "L1Extrapolation.h"
18 
19 #include <iostream>
20 
21 // using namespace std;
22 using std::cout;
23 using std::endl;
24 using std::vector;
25 
26 void L1Algo::InvertCholetsky(fvec a[15]) {
27  fvec d[5], uud, u[5][5];
28  for (int i = 0; i < 5; i++) {
29  d[i] = 0.f;
30  for (int j = 0; j < 5; j++)
31  u[i][j] = 0.f;
32  }
33 
34  for (int i = 0; i < 5; i++) {
35  uud = 0.f;
36  for (int j = 0; j < i; j++)
37  uud += u[j][i] * u[j][i] * d[j];
38  uud = a[i * (i + 3) / 2] - uud;
39 
40  fvec smallval = 1.e-12;
41  fvec initialised = fvec(uud < smallval);
42  uud = ((!initialised) & uud) + (smallval & initialised);
43 
44  d[i] = uud / fabs(uud);
45  u[i][i] = sqrt(fabs(uud));
46 
47  for (int j = i + 1; j < 5; j++) {
48  uud = 0.f;
49  for (int k = 0; k < i; k++)
50  uud += u[k][i] * u[k][j] * d[k];
51  uud = a[j * (j + 1) / 2 + i] /*a[i][j]*/ - uud;
52  u[i][j] = d[i] / u[i][i] * uud;
53  }
54  }
55 
56  fvec u1[5];
57 
58  for (int i = 0; i < 5; i++) {
59  u1[i] = u[i][i];
60  u[i][i] = 1.f / u[i][i];
61  }
62  for (int i = 0; i < 4; i++) {
63  u[i][i + 1] = -u[i][i + 1] * u[i][i] * u[i + 1][i + 1];
64  }
65  for (int i = 0; i < 3; i++) {
66  u[i][i + 2] = u[i][i + 1] * u1[i + 1] * u[i + 1][i + 2]
67  - u[i][i + 2] * u[i][i] * u[i + 2][i + 2];
68  }
69  for (int i = 0; i < 2; i++) {
70  u[i][i + 3] = u[i][i + 2] * u1[i + 2] * u[i + 2][i + 3]
71  - u[i][i + 3] * u[i][i] * u[i + 3][i + 3];
72  u[i][i + 3] -=
73  u[i][i + 1] * u1[i + 1]
74  * (u[i + 1][i + 2] * u1[i + 2] * u[i + 2][i + 3] - u[i + 1][i + 3]);
75  }
76  u[0][4] = u[0][2] * u1[2] * u[2][4] - u[0][4] * u[0][0] * u[4][4];
77  u[0][4] +=
78  u[0][1] * u1[1]
79  * (u[1][4] - u[1][3] * u1[3] * u[3][4] - u[1][2] * u1[2] * u[2][4]);
80  u[0][4] +=
81  u[3][4] * u1[3]
82  * (u[0][3] - u1[2] * u[2][3] * (u[0][2] - u[0][1] * u1[1] * u[1][2]));
83 
84  for (int i = 0; i < 5; i++)
85  a[i + 10] = u[i][4] * d[4] * u[4][4];
86  for (int i = 0; i < 4; i++)
87  a[i + 6] = u[i][3] * u[3][3] * d[3] + u[i][4] * u[3][4] * d[4];
88  for (int i = 0; i < 3; i++)
89  a[i + 3] = u[i][2] * u[2][2] * d[2] + u[i][3] * u[2][3] * d[3]
90  + u[i][4] * u[2][4] * d[4];
91  for (int i = 0; i < 2; i++)
92  a[i + 1] = u[i][1] * u[1][1] * d[1] + u[i][2] * u[1][2] * d[2]
93  + u[i][3] * u[1][3] * d[3] + u[i][4] * u[1][4] * d[4];
94  a[0] = u[0][0] * u[0][0] * d[0] + u[0][1] * u[0][1] * d[1]
95  + u[0][2] * u[0][2] * d[2] + u[0][3] * u[0][3] * d[3]
96  + u[0][4] * u[0][4] * d[4];
97 }
98 
99 void L1Algo::MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]) {
100  K[0][0] =
101  C[0] * V[0] + C[1] * V[1] + C[3] * V[3] + C[6] * V[6] + C[10] * V[10];
102  K[0][1] =
103  C[0] * V[1] + C[1] * V[2] + C[3] * V[4] + C[6] * V[7] + C[10] * V[11];
104  K[0][2] =
105  C[0] * V[3] + C[1] * V[4] + C[3] * V[5] + C[6] * V[8] + C[10] * V[12];
106  K[0][3] =
107  C[0] * V[6] + C[1] * V[7] + C[3] * V[8] + C[6] * V[9] + C[10] * V[13];
108  K[0][4] =
109  C[0] * V[10] + C[1] * V[11] + C[3] * V[12] + C[6] * V[13] + C[10] * V[14];
110 
111  K[1][0] =
112  C[1] * V[0] + C[2] * V[1] + C[4] * V[3] + C[7] * V[6] + C[11] * V[10];
113  K[1][1] =
114  C[1] * V[1] + C[2] * V[2] + C[4] * V[4] + C[7] * V[7] + C[11] * V[11];
115  K[1][2] =
116  C[1] * V[3] + C[2] * V[4] + C[4] * V[5] + C[7] * V[8] + C[11] * V[12];
117  K[1][3] =
118  C[1] * V[6] + C[2] * V[7] + C[4] * V[8] + C[7] * V[9] + C[11] * V[13];
119  K[1][4] =
120  C[1] * V[10] + C[2] * V[11] + C[4] * V[12] + C[7] * V[13] + C[11] * V[14];
121 
122  K[2][0] =
123  C[3] * V[0] + C[4] * V[1] + C[5] * V[3] + C[8] * V[6] + C[12] * V[10];
124  K[2][1] =
125  C[3] * V[1] + C[4] * V[2] + C[5] * V[4] + C[8] * V[7] + C[12] * V[11];
126  K[2][2] =
127  C[3] * V[3] + C[4] * V[4] + C[5] * V[5] + C[8] * V[8] + C[12] * V[12];
128  K[2][3] =
129  C[3] * V[6] + C[4] * V[7] + C[5] * V[8] + C[8] * V[9] + C[12] * V[13];
130  K[2][4] =
131  C[3] * V[10] + C[4] * V[11] + C[5] * V[12] + C[8] * V[13] + C[12] * V[14];
132 
133  K[3][0] =
134  C[6] * V[0] + C[7] * V[1] + C[8] * V[3] + C[9] * V[6] + C[13] * V[10];
135  K[3][1] =
136  C[6] * V[1] + C[7] * V[2] + C[8] * V[4] + C[9] * V[7] + C[13] * V[11];
137  K[3][2] =
138  C[6] * V[3] + C[7] * V[4] + C[8] * V[5] + C[9] * V[8] + C[13] * V[12];
139  K[3][3] =
140  C[6] * V[6] + C[7] * V[7] + C[8] * V[8] + C[9] * V[9] + C[13] * V[13];
141  K[3][4] =
142  C[6] * V[10] + C[7] * V[11] + C[8] * V[12] + C[9] * V[13] + C[13] * V[14];
143 
144  K[4][0] =
145  C[10] * V[0] + C[11] * V[1] + C[12] * V[3] + C[13] * V[6] + C[14] * V[10];
146  K[4][1] =
147  C[10] * V[1] + C[11] * V[2] + C[12] * V[4] + C[13] * V[7] + C[14] * V[11];
148  K[4][2] =
149  C[10] * V[3] + C[11] * V[4] + C[12] * V[5] + C[13] * V[8] + C[14] * V[12];
150  K[4][3] =
151  C[10] * V[6] + C[11] * V[7] + C[12] * V[8] + C[13] * V[9] + C[14] * V[13];
152  K[4][4] = C[10] * V[10] + C[11] * V[11] + C[12] * V[12] + C[13] * V[13]
153  + C[14] * V[14];
154 }
155 
156 void L1Algo::MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]) {
157  K[0] = C[0][0] * V[0] + C[0][1] * V[1] + C[0][2] * V[3] + C[0][3] * V[6]
158  + C[0][4] * V[10];
159 
160  K[1] = C[1][0] * V[0] + C[1][1] * V[1] + C[1][2] * V[3] + C[1][3] * V[6]
161  + C[1][4] * V[10];
162  K[2] = C[1][0] * V[1] + C[1][1] * V[2] + C[1][2] * V[4] + C[1][3] * V[7]
163  + C[1][4] * V[11];
164 
165  K[3] = C[2][0] * V[0] + C[2][1] * V[1] + C[2][2] * V[3] + C[2][3] * V[6]
166  + C[2][4] * V[10];
167  K[4] = C[2][0] * V[1] + C[2][1] * V[2] + C[2][2] * V[4] + C[2][3] * V[7]
168  + C[2][4] * V[11];
169  K[5] = C[2][0] * V[3] + C[2][1] * V[4] + C[2][2] * V[5] + C[2][3] * V[8]
170  + C[2][4] * V[12];
171 
172  K[6] = C[3][0] * V[0] + C[3][1] * V[1] + C[3][2] * V[3] + C[3][3] * V[6]
173  + C[3][4] * V[10];
174  K[7] = C[3][0] * V[1] + C[3][1] * V[2] + C[3][2] * V[4] + C[3][3] * V[7]
175  + C[3][4] * V[11];
176  K[8] = C[3][0] * V[3] + C[3][1] * V[4] + C[3][2] * V[5] + C[3][3] * V[8]
177  + C[3][4] * V[12];
178  K[9] = C[3][0] * V[6] + C[3][1] * V[7] + C[3][2] * V[8] + C[3][3] * V[9]
179  + C[3][4] * V[13];
180 
181  K[10] = C[4][0] * V[0] + C[4][1] * V[1] + C[4][2] * V[3] + C[4][3] * V[6]
182  + C[4][4] * V[10];
183  K[11] = C[4][0] * V[1] + C[4][1] * V[2] + C[4][2] * V[4] + C[4][3] * V[7]
184  + C[4][4] * V[11];
185  K[12] = C[4][0] * V[3] + C[4][1] * V[4] + C[4][2] * V[5] + C[4][3] * V[8]
186  + C[4][4] * V[12];
187  K[13] = C[4][0] * V[6] + C[4][1] * V[7] + C[4][2] * V[8] + C[4][3] * V[9]
188  + C[4][4] * V[13];
189  K[14] = C[4][0] * V[10] + C[4][1] * V[11] + C[4][2] * V[12] + C[4][3] * V[13]
190  + C[4][4] * V[14];
191 }
192 
193 void L1Algo::MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]) {
194  r_out[0] = r_in[0] * C[0] + r_in[1] * C[1] + r_in[2] * C[3] + r_in[3] * C[6]
195  + r_in[4] * C[10];
196  r_out[1] = r_in[0] * C[1] + r_in[1] * C[2] + r_in[2] * C[4] + r_in[3] * C[7]
197  + r_in[4] * C[11];
198  r_out[2] = r_in[0] * C[3] + r_in[1] * C[4] + r_in[2] * C[5] + r_in[3] * C[8]
199  + r_in[4] * C[12];
200  r_out[3] = r_in[0] * C[6] + r_in[1] * C[7] + r_in[2] * C[8] + r_in[3] * C[9]
201  + r_in[4] * C[13];
202  r_out[4] = r_in[0] * C[10] + r_in[1] * C[11] + r_in[2] * C[12]
203  + r_in[3] * C[13] + r_in[4] * C[14];
204 }
205 
206 void L1Algo::FilterTracks(fvec const r[5],
207  fvec const C[15],
208  fvec const m[5],
209  fvec const V[15],
210  fvec R[5],
211  fvec W[15],
212  fvec* chi2) {
213  fvec S[15];
214  for (int i = 0; i < 15; i++) {
215  if (W) W[i] = C[i];
216  S[i] = C[i] + V[i];
217  }
218 
219  InvertCholetsky(S);
220 
221  fvec dzeta[5];
222 
223  for (int i = 0; i < 5; i++)
224  dzeta[i] = m[i] - r[i];
225 
226  if (W && R) {
227  for (int i = 0; i < 5; i++)
228  R[i] = r[i];
229 
230  fvec K[5][5];
231  MultiplySS(C, S, K);
232 
233  fvec KC[15];
234  MultiplyMS(K, C, KC);
235  for (int i = 0; i < 15; i++)
236  W[i] -= KC[i];
237 
238  fvec kd;
239  for (int i = 0; i < 5; i++) {
240  kd = 0.f;
241  for (int j = 0; j < 5; j++)
242  kd += K[i][j] * dzeta[j];
243  R[i] += kd;
244  }
245  }
246  if (chi2) {
247  fvec S_dzeta[5];
248  MultiplySR(S, dzeta, S_dzeta);
249  *chi2 = dzeta[0] * S_dzeta[0] + dzeta[1] * S_dzeta[1]
250  + dzeta[2] * S_dzeta[2] + dzeta[3] * S_dzeta[3]
251  + dzeta[4] * S_dzeta[4];
252  }
253 }
254 
255 void L1Algo::CAMergeClones() {
256  // vector<unsigned short> FirstHit;
257  // vector<unsigned short> LastHit;
258  // vector<THitI> FirstHitIndex;
259  // vector<THitI> LastHitIndex;
260  // vector<unsigned short> Neighbour;
261  // vector<float> TrackChi2;
262  vector<bool> IsNext;
263  vector<bool> IsUsed;
264 
265  // vector< L1Track > vTracksNew;
266  vTracksNew.clear();
267  vTracksNew.reserve(NTracksIsecAll);
268  // vector< THitI > vRecoHitsNew;
269  vRecoHitsNew.clear();
270  vRecoHitsNew.reserve(vRecoHits.size());
271 
272  FirstHit.resize(NTracksIsecAll);
273  LastHit.resize(NTracksIsecAll);
274  FirstHitIndex.resize(NTracksIsecAll);
275  LastHitIndex.resize(NTracksIsecAll);
276  IsUsed.resize(NTracksIsecAll);
277  TrackChi2.resize(NTracksIsecAll);
278  Neighbour.resize(NTracksIsecAll);
279  IsNext.resize(NTracksIsecAll);
280 
281  THitI start_hit = 0;
282  unsigned short ista = 0;
283 #ifdef OMP
284 #pragma omp parallel for
285 #endif
286 
287  for (unsigned short iTr = 0; iTr < NTracksIsecAll; iTr++) {
288  FirstHitIndex[iTr] = start_hit;
289  ista = (*vSFlag)[(*vStsHits)[vRecoHits[start_hit]].f] / 4;
290  FirstHit[iTr] = ista;
291  start_hit += vTracks[iTr].NHits - 1;
292  LastHitIndex[iTr] = start_hit;
293  ista = (*vSFlag)[(*vStsHits)[vRecoHits[start_hit]].f] / 4;
294  LastHit[iTr] = ista;
295  start_hit++;
296 
297  IsUsed[iTr] = 0;
298  Neighbour[iTr] = 50000;
299  TrackChi2[iTr] = 100000;
300  IsNext[iTr] = 0;
301  }
302 
303  L1KFTrackFitter();
304  //KFTrackFitter_simple();
305 
306  L1TrackPar Tb;
307  L1TrackPar Tf;
308  L1FieldValue fBm, fBb, fBf _fvecalignment;
309  L1FieldRegion fld _fvecalignment;
310 #ifdef OMP
311 #pragma omp parallel for
312 #endif
313  for (int iTr = 0; iTr < static_cast<unsigned short>(NTracksIsecAll); iTr++) {
314  if (static_cast<int>(vTracks[iTr].NHits) > (NStations - 3)) continue;
315  for (int jTr = 0; jTr < static_cast<unsigned short>(NTracksIsecAll);
316  jTr++) {
317  if (static_cast<int>(vTracks[jTr].NHits) > (NStations - 3)) continue;
318 
319  if (iTr == jTr) continue;
320  // if(vTracks[iTr].n != vTracks[jTr].n) continue;
321 
322 
323  // if (fabs(vTracks[iTr].fTrackTime - vTracks[jTr].fTrackTime) > 3) continue;
324 
325  unsigned short dist = 0;
326  unsigned short stab = 0, staf = 0;
327  bool IsNextTemp = 0;
328  //if((vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])*(vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])
329  // > 9*(vTracks[iTr].CFirst[14]+vTracks[jTr].CFirst[14]) ) continue;
330  if (FirstHit[iTr] > LastHit[jTr]) {
331  dist = FirstHit[iTr] - LastHit[jTr];
332 
333  stab = FirstHit[iTr];
334  staf = LastHit[jTr];
335  IsNextTemp = 1;
336 
337  Tb.x = vTracks[iTr].TFirst[0];
338  Tb.y = vTracks[iTr].TFirst[1];
339  Tb.tx = vTracks[iTr].TFirst[2];
340  Tb.ty = vTracks[iTr].TFirst[3];
341  Tb.qp = vTracks[iTr].TFirst[4];
342  Tb.z = vTracks[iTr].TFirst[5];
343  Tb.C00 = vTracks[iTr].CFirst[0];
344  Tb.C10 = vTracks[iTr].CFirst[1];
345  Tb.C11 = vTracks[iTr].CFirst[2];
346  Tb.C20 = vTracks[iTr].CFirst[3];
347  Tb.C21 = vTracks[iTr].CFirst[4];
348  Tb.C22 = vTracks[iTr].CFirst[5];
349  Tb.C30 = vTracks[iTr].CFirst[6];
350  Tb.C31 = vTracks[iTr].CFirst[7];
351  Tb.C32 = vTracks[iTr].CFirst[8];
352  Tb.C33 = vTracks[iTr].CFirst[9];
353  Tb.C40 = vTracks[iTr].CFirst[10];
354  Tb.C41 = vTracks[iTr].CFirst[11];
355  Tb.C42 = vTracks[iTr].CFirst[12];
356  Tb.C43 = vTracks[iTr].CFirst[13];
357  Tb.C44 = vTracks[iTr].CFirst[14];
358 
359  Tf.x = vTracks[jTr].TLast[0];
360  Tf.y = vTracks[jTr].TLast[1];
361  Tf.tx = vTracks[jTr].TLast[2];
362  Tf.ty = vTracks[jTr].TLast[3];
363  Tf.qp = vTracks[jTr].TLast[4];
364  Tf.z = vTracks[jTr].TLast[5];
365  Tf.C00 = vTracks[jTr].CLast[0];
366  Tf.C10 = vTracks[jTr].CLast[1];
367  Tf.C11 = vTracks[jTr].CLast[2];
368  Tf.C20 = vTracks[jTr].CLast[3];
369  Tf.C21 = vTracks[jTr].CLast[4];
370  Tf.C22 = vTracks[jTr].CLast[5];
371  Tf.C30 = vTracks[jTr].CLast[6];
372  Tf.C31 = vTracks[jTr].CLast[7];
373  Tf.C32 = vTracks[jTr].CLast[8];
374  Tf.C33 = vTracks[jTr].CLast[9];
375  Tf.C40 = vTracks[jTr].CLast[10];
376  Tf.C41 = vTracks[jTr].CLast[11];
377  Tf.C42 = vTracks[jTr].CLast[12];
378  Tf.C43 = vTracks[jTr].CLast[13];
379  Tf.C44 = vTracks[jTr].CLast[14];
380  //std::cout << "!!!!!!! Chi2 !!!!!! "<<vTracks[iTr].TFirst[0]<<" "<<vTracks[jTr].TLast[0]<<std::endl;
381  }
382  // if(FirstHit[jTr] > LastHit[iTr])
383  // {
384  // dist = FirstHit[jTr] - LastHit[iTr];
385  //
386  // stab = FirstHit[jTr];
387  // staf = LastHit[iTr];
388  //
389  // Tb.x = vTracks[jTr].TFirst[0];
390  // Tb.y = vTracks[jTr].TFirst[1];
391  // Tb.tx = vTracks[jTr].TFirst[2];
392  // Tb.ty = vTracks[jTr].TFirst[3];
393  // Tb.qp = vTracks[jTr].TFirst[4];
394  // Tb.z = vTracks[jTr].TFirst[5];
395  // Tb.C00 = vTracks[jTr].CFirst[0];
396  // Tb.C10 = vTracks[jTr].CFirst[1];
397  // Tb.C11 = vTracks[jTr].CFirst[2];
398  // Tb.C20 = vTracks[jTr].CFirst[3];
399  // Tb.C21 = vTracks[jTr].CFirst[4];
400  // Tb.C22 = vTracks[jTr].CFirst[5];
401  // Tb.C30 = vTracks[jTr].CFirst[6];
402  // Tb.C31 = vTracks[jTr].CFirst[7];
403  // Tb.C32 = vTracks[jTr].CFirst[8];
404  // Tb.C33 = vTracks[jTr].CFirst[9];
405  // Tb.C40 = vTracks[jTr].CFirst[10];
406  // Tb.C41 = vTracks[jTr].CFirst[11];
407  // Tb.C42 = vTracks[jTr].CFirst[12];
408  // Tb.C43 = vTracks[jTr].CFirst[13];
409  // Tb.C44 = vTracks[jTr].CFirst[14];
410  //
411  // Tf.x = vTracks[iTr].TLast[0];
412  // Tf.y = vTracks[iTr].TLast[1];
413  // Tf.tx = vTracks[iTr].TLast[2];
414  // Tf.ty = vTracks[iTr].TLast[3];
415  // Tf.qp = vTracks[iTr].TLast[4];
416  // Tf.z = vTracks[iTr].TLast[5];
417  // Tf.C00 = vTracks[iTr].CLast[0];
418  // Tf.C10 = vTracks[iTr].CLast[1];
419  // Tf.C11 = vTracks[iTr].CLast[2];
420  // Tf.C20 = vTracks[iTr].CLast[3];
421  // Tf.C21 = vTracks[iTr].CLast[4];
422  // Tf.C22 = vTracks[iTr].CLast[5];
423  // Tf.C30 = vTracks[iTr].CLast[6];
424  // Tf.C31 = vTracks[iTr].CLast[7];
425  // Tf.C32 = vTracks[iTr].CLast[8];
426  // Tf.C33 = vTracks[iTr].CLast[9];
427  // Tf.C40 = vTracks[iTr].CLast[10];
428  // Tf.C41 = vTracks[iTr].CLast[11];
429  // Tf.C42 = vTracks[iTr].CLast[12];
430  // Tf.C43 = vTracks[iTr].CLast[13];
431  // Tf.C44 = vTracks[iTr].CLast[14];
432  // }
433 
434  if (dist == 0) continue;
435  //if(((Tf.qp - Tb.qp)*(Tf.qp - Tb.qp)/(Tb.C44+Tf.C44))[0] > 25*10*7) continue;
436  if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
437  // if (fabs (Tf.time[0] - Tb.time[0]) > 500000) continue;
438  unsigned short stam;
439 
440  vStations[staf].fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
441  vStations[stab].fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
442  if (dist > 1)
443  stam = staf + 1;
444  else
445  stam = staf - 1;
446  fvec zm = vStations[stam].z;
447  fvec xm = 0.5 * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
448  fvec ym = 0.5 * (Tb.y + Tb.ty * (zm - Tb.z) + Tb.y + Tb.ty * (zm - Tb.z));
449  vStations[stam].fieldSlice.GetFieldValue(xm, ym, fBm);
450  fld.Set(fBb, Tb.z, fBm, zm, fBf, Tf.z);
451 
452  fvec zMiddle = 0.5 * (Tb.z + Tf.z);
453 
454  L1Extrapolate(Tf, zMiddle, Tf.qp, fld);
455  L1Extrapolate(Tb, zMiddle, Tb.qp, fld);
456 
457  fvec Chi2Tracks = 0.f;
458  FilterTracks(&(Tf.x), &(Tf.C00), &(Tb.x), &(Tb.C00), 0, 0, &Chi2Tracks);
459  if (Chi2Tracks[0] > 50) continue;
460  if (Chi2Tracks[0] < TrackChi2[iTr] || Chi2Tracks[0] < TrackChi2[jTr]) {
461  if (Neighbour[iTr] < static_cast<unsigned short>(50000)) {
462  Neighbour[Neighbour[iTr]] = 50000;
463  TrackChi2[Neighbour[iTr]] = 100000;
464  IsNext[Neighbour[iTr]] = 0;
465  }
466  if (Neighbour[jTr] < static_cast<unsigned short>(50000)) {
467  Neighbour[Neighbour[jTr]] = 50000;
468  TrackChi2[Neighbour[jTr]] = 100000;
469  IsNext[Neighbour[jTr]] = 0;
470  }
471  Neighbour[iTr] = jTr;
472  Neighbour[jTr] = iTr;
473  TrackChi2[iTr] = Chi2Tracks[0];
474  TrackChi2[jTr] = Chi2Tracks[0];
475  IsNext[iTr] = IsNextTemp;
476  IsNext[jTr] = (!IsNextTemp);
477  }
478  }
479  }
480  for (int iTr = 0; iTr < static_cast<unsigned short>(NTracksIsecAll); iTr++) {
481  if (IsUsed[iTr]) continue;
482 
483  vTracksNew.push_back(vTracks[iTr]);
484  if (!IsNext[iTr])
485  for (THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++)
486  vRecoHitsNew.push_back(vRecoHits[HI]);
487 
488  if (Neighbour[iTr] < 50000) {
489  IsUsed[Neighbour[iTr]] = 1;
490  vTracksNew.back().NHits += vTracks[Neighbour[iTr]].NHits;
491  for (THitI HI = FirstHitIndex[Neighbour[iTr]];
492  HI <= LastHitIndex[Neighbour[iTr]];
493  HI++)
494  vRecoHitsNew.push_back(vRecoHits[HI]);
495  }
496 
497  if (IsNext[iTr])
498  for (THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++)
499  vRecoHitsNew.push_back(vRecoHits[HI]);
500  }
501  //vTracks.resize(vTracksNew.size());
502  for (unsigned short iTr = 0; iTr < vTracksNew.size(); iTr++)
503  vTracks[iTr] = vTracksNew[iTr];
504  // vRecoHits.resize(vRecoHitsNew.size());
505  for (THitI iHi = 0; iHi < vRecoHitsNew.size(); iHi++)
506  vRecoHits[iHi] = vRecoHitsNew[iHi];
507 
508 
509  NHitsIsecAll = vRecoHitsNew.size();
510  NTracksIsecAll = vTracksNew.size();
511 
512 
513  //std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!! new "<<vTracksNew.size()<<" old "<< vTracks.size()<<std::endl;
514 }
L1TrackPar::C10
fvec C10
Definition: L1TrackPar.h:9
L1Algo::LastHit
L1Vector< unsigned short > LastHit
Definition: L1Algo.h:275
L1TrackPar::qp
fvec qp
Definition: L1TrackPar.h:9
L1TrackPar::t
fvec t
Definition: L1TrackPar.h:9
L1Algo::CAMergeClones
void CAMergeClones()
Definition: L1CAMergeClones.cxx:255
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
L1Algo.h
L1TrackPar::C20
fvec C20
Definition: L1TrackPar.h:9
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
L1TrackPar::C41
fvec C41
Definition: L1TrackPar.h:10
L1Algo::FirstHitIndex
L1Vector< THitI > FirstHitIndex
Definition: L1Algo.h:276
L1Algo::MultiplySS
void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5])
Definition: L1CAMergeClones.cxx:99
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
L1TrackPar::C30
fvec C30
Definition: L1TrackPar.h:9
L1TrackPar::C31
fvec C31
Definition: L1TrackPar.h:9
fvec
F32vec4 fvec
Definition: L1/vectors/P4_F32vec4.h:249
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, const L1FieldRegion &F, fvec *w=0)
Definition: L1Extrapolation.h:314
L1Algo::Neighbour
L1Vector< unsigned short > Neighbour
Definition: L1Algo.h:278
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
L1Algo::MultiplySR
void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5])
Definition: L1CAMergeClones.cxx:193
L1TrackPar::C42
fvec C42
Definition: L1TrackPar.h:10
L1TrackPar::z
fvec z
Definition: L1TrackPar.h:9
L1Algo::TrackChi2
L1Vector< float > TrackChi2
Definition: L1Algo.h:279
L1Algo::vRecoHitsNew
L1Vector< THitI > vRecoHitsNew
Definition: L1Algo.h:282
L1TrackPar::ty
fvec ty
Definition: L1TrackPar.h:9
L1TrackPar::C32
fvec C32
Definition: L1TrackPar.h:9
L1Algo::FirstHit
L1Vector< unsigned short > FirstHit
Definition: L1Algo.h:274
L1TrackPar::y
fvec y
Definition: L1TrackPar.h:9
L1Algo::NTracksIsecAll
unsigned int NTracksIsecAll
Definition: L1Algo.h:364
THitI
unsigned int THitI
Definition: L1StsHit.h:8
L1Algo::vTracksNew
L1Vector< L1Track > vTracksNew
Definition: L1Algo.h:283
L1Algo::NStations
int NStations
Definition: L1Algo.h:333
L1TrackPar::C44
fvec C44
Definition: L1TrackPar.h:10
L1FieldRegion
Definition: L1Field.h:85
L1Algo::NHitsIsecAll
unsigned int NHitsIsecAll
Definition: L1Algo.h:363
L1Algo::InvertCholetsky
void InvertCholetsky(fvec a[15])
--— Subroutines used by L1Algo::CAMergeClones() ---—
Definition: L1CAMergeClones.cxx:26
d
double d
Definition: P4_F64vec2.h:24
L1TrackPar::tx
fvec tx
Definition: L1TrackPar.h:9
L1TrackPar::C21
fvec C21
Definition: L1TrackPar.h:9
L1Algo::vTracks
L1Vector< L1Track > vTracks
Definition: L1Algo.h:355
L1TrackPar::x
fvec x
Definition: L1TrackPar.h:9
L1TrackPar::C00
fvec C00
Definition: L1TrackPar.h:9
L1Extrapolation.h
L1TrackPar::C33
fvec C33
Definition: L1TrackPar.h:9
L1Algo::vRecoHits
L1Vector< THitI > vRecoHits
Definition: L1Algo.h:356
L1Algo::_fvecalignment
L1Station vStations[MaxNStations] _fvecalignment
Definition: L1Algo.h:336
L1Algo::L1KFTrackFitter
void L1KFTrackFitter()
Definition: L1TrackFitter.cxx:335
L1Algo::FilterTracks
void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5], fvec W[15], fvec *chi2)
Definition: L1CAMergeClones.cxx:206
L1TrackPar
Definition: L1TrackPar.h:6
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
L1TrackPar::C55
fvec C55
Definition: L1TrackPar.h:10
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
L1TrackPar::C40
fvec C40
Definition: L1TrackPar.h:10
L1TrackPar::C11
fvec C11
Definition: L1TrackPar.h:9
L1Algo::MultiplyMS
void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15])
Definition: L1CAMergeClones.cxx:156
L1TrackPar::C43
fvec C43
Definition: L1TrackPar.h:10
L1FieldValue
Definition: L1Field.h:11
L1Algo::LastHitIndex
L1Vector< THitI > LastHitIndex
Definition: L1Algo.h:277
L1TrackPar::C22
fvec C22
Definition: L1TrackPar.h:9