Coverity
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuTMU.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliTRDgtuTMU.cxx 28397 2008-09-02 09:33:00Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  Track Matching Unit (TMU) simulation                                  //
21 //                                                                        //
22 //  Author: J. Klein (Jochen.Klein@cern.ch)                               //
23 //                                                                        //
24 ////////////////////////////////////////////////////////////////////////////
25
26 #include "TTree.h"
27 #include "TList.h"
28 #include "TVectorD.h"
29 #include "TMath.h"
30
31 #include "AliESDEvent.h"
32 #include "AliESDTrdTrack.h"
33
34 #include "AliLog.h"
35 #include "AliTRDgeometry.h"
36 #include "AliTRDpadPlane.h"
37
38 #include "AliTRDgtuParam.h"
39 #include "AliTRDgtuTMU.h"
40 #include "AliTRDtrackGTU.h"
41
42 ClassImp(AliTRDgtuTMU)
43
44 AliTRDgtuTMU::AliTRDgtuTMU(Int_t stack, Int_t sector) :
45   TObject(),
46   fTracklets(0x0),
47   fTrackletsPostInput(0x0),
48   fZChannelTracklets(0x0),
49   fTracks(0x0),
50   fGtuParam(0x0),
51   fStack(-1),
52   fSector(-1)
53 {
54   // constructor which initializes the position information of the TMU
55
56   fGtuParam = AliTRDgtuParam::Instance();
57
58   // store tracklets per link
59   fTracklets = new TObjArray*[fGtuParam->GetNLinks()];
60   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
61     fTracklets[iLink] = new TObjArray();
62   }
63
64   // tracklets after input units per layer
65   fTrackletsPostInput = new TObjArray*[fGtuParam->GetNLayers()];
66   fZChannelTracklets = new TList*[fGtuParam->GetNLayers()];
67   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
68     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
69     fTrackletsPostInput[layer] = new TObjArray();
70   }
71
72   fTracks = new TList*[fGtuParam->GetNZChannels()];
73   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
74       fTracks[zch] = new TList[fGtuParam->GetNRefLayers()];
75   }
76
77   if (stack > -1)
78       SetStack(stack);
79   if (sector > -1)
80       SetSector(sector);
81 }
82
83 AliTRDgtuTMU::~AliTRDgtuTMU()
84 {
85   // destructor
86
87   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
88     delete [] fTracks[zch];
89   }
90   delete [] fTracks;
91   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
92     delete [] fZChannelTracklets[layer];
93     //    delete [] fTrackletsPostInput[layer];
94   }
95   delete [] fZChannelTracklets;
96   //  delete [] fTrackletsPostInput;
97
98   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
99     delete fTracklets[iLink];
100   }
101   delete [] fTracklets;
102 }
103
104 Bool_t AliTRDgtuTMU::SetSector(Int_t sector)
105 {
106   // set the sector
107
108   if (sector > -1 && sector < fGtuParam->GetGeo()->Nsector() ) {
109     fSector = sector;
110     return kTRUE;
111   }
112
113   AliError(Form("Invalid sector given: %i", sector));
114   return kFALSE;
115 }
116
117 Bool_t AliTRDgtuTMU::SetStack(Int_t stack)
118 {
119   // Set the stack (necessary for tracking)
120
121   if (stack > -1 && stack < fGtuParam->GetGeo()->Nstack() ) {
122     fStack = stack;
123     return kTRUE;
124   }
125
126   AliError(Form("Invalid stack given: %i", stack));
127   return kFALSE;
128 }
129
130 Bool_t AliTRDgtuTMU::Reset()
131 {
132   // delete all tracks
133
134   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
135     for (Int_t reflayeridx = 0; reflayeridx < fGtuParam->GetNRefLayers(); reflayeridx++) {
136       fTracks[zch][reflayeridx].Clear();
137     }
138   }
139
140   // delete all added tracklets
141   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
142     fTracklets[iLink]->Clear();
143   }
144   for (Int_t layer = 0; layer < fGtuParam->GetNLinks()/2; layer++) {
145     fTrackletsPostInput[layer]->Clear();
146     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++)
147       fZChannelTracklets[layer][zch].Clear();
148   }
149
150   fSector = -1;
151   fStack = -1;
152
153   return kTRUE;
154 }
155
156 Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletGTU *tracklet, Int_t link)
157 {
158   // add a tracklet on the given link
159
160   fTracklets[link]->Add(tracklet);
161   return kTRUE;
162 }
163
164
165 Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
166 {
167   // performs the analysis as in a TMU module of the GTU, i. e.
168   // track matching
169   // calculation of track parameteres (pt, deflection, ???)
170
171   if (fStack < 0 || fSector < 0) {
172     AliError("No valid stack/sector set for this TMU! No tracking!");
173     return kFALSE;
174   }
175
176   // ----- Input units -----
177   AliDebug(1,"--------- Running Input units ----------");
178   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
179     if (!RunInputUnit(layer)) {
180       AliError(Form("Input unit in layer %i failed", layer));
181       return kFALSE;
182     }
183   }
184
185   // ----- Z-channel units -----
186   AliDebug(1,"--------- Running Z-channel units ----------");
187   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
188     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
189     if (!RunZChannelUnit(layer)) {
190       AliError(Form("Z-Channel unit in layer %i failed", layer));
191       return kFALSE;
192     }
193   }
194
195   // ----- track finding -----
196   AliDebug(1,"--------- Running tracking units ----------");
197   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
198     if (!RunTrackFinder(zch, ListOfTracks)) {
199       AliError(Form("Track Finder in z-channel %i failed", zch));
200       return kFALSE;
201     }
202   }
203
204   // ----- Track Merging -----
205   if (!RunTrackMerging(ListOfTracks)) {
206     AliError("Track merging failed");
207     return kFALSE;
208   }
209
210   // ----- track reconstruction -----
211   if (!RunTrackReconstruction(ListOfTracks)) {
212     AliError("Track reconstruction failed");
213     return kFALSE;
214   }
215
216   if (esd) {
217       TIter next(ListOfTracks);
218       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
219           AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
220           esd->AddTrdTrack(trdtrack);
221           delete trdtrack;
222       }
223   }
224
225   return kTRUE;
226 }
227
228 Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
229 {
230   // precalculations for the tracking and reconstruction
231
232   Int_t iTrkl0 = 0; // A-side tracklet
233   Int_t iTrkl1 = 0; // B-side tracklet
234
235   while ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) || (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast())) {
236     if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
237       fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
238       iTrkl0++;
239     }
240     else {
241       if (iTrkl0 >= fTracklets[2*layer + 0]->GetEntriesFast()) {
242         fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
243         iTrkl1++;
244       }
245       else {
246         if (((AliTRDtrackletGTU*) fTracklets[2*layer + 1]->At(iTrkl1))->GetZbin() <
247             ((AliTRDtrackletGTU*) fTracklets[2*layer + 0]->At(iTrkl0))->GetZbin()) {
248           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
249           iTrkl1++;
250
251         }
252         else {
253           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
254           iTrkl0++;
255         }
256       }
257     }
258   }
259
260   TIter next(fTrackletsPostInput[layer]);
261
262   while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
263     trk->SetIndex(fTrackletsPostInput[layer]->IndexOf(trk));
264
265     Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
266     alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
267     // wrapping projected alpha as in hardware
268     if ((alpha < -64) || (alpha > 63))
269       AliError(Form("alpha out of range: %i", alpha));
270     alpha = alpha & 0x7f;
271     if (alpha & 0x40)
272       trk->SetAlpha(0xffffffc0 | alpha);
273     else
274       trk->SetAlpha(alpha);
275
276     Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer)); //??? sign?
277     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
278     trk->SetYProj(yproj);
279
280     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
281
282     AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i",
283                       trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
284                       trk->GetYProj(), trk->GetAlpha(), trk->GetPID() ));
285   }
286   return kTRUE;
287 }
288
289 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
290 {
291   // run the z-channel unit
292
293   TIter next(fTrackletsPostInput[layer]);
294
295   while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
296     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
297       if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
298         trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
299
300         TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
301         AliTRDtrackletGTU *t = 0x0;
302         while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
303           if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
304               ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
305             break;
306           }
307         }
308         if (t)
309           fZChannelTracklets[layer][zch].AddAfter(t, trk);
310         else
311           fZChannelTracklets[layer][zch].AddFirst(trk);
312       }
313     }
314     AliDebug(10, Form("stack %d, layer %2d: 0x%08x Z(2..0)=%i/%i/%i",
315                       fStack, layer, trk->GetTrackletWord(),
316                       fGtuParam->IsInZChannel(fStack, layer, 2, trk->GetZbin()) ? trk->GetSubChannel(2) : -1,
317                       fGtuParam->IsInZChannel(fStack, layer, 1, trk->GetZbin()) ? trk->GetSubChannel(1) : -1,
318                       fGtuParam->IsInZChannel(fStack, layer, 0, trk->GetZbin()) ? trk->GetSubChannel(0) : -1
319                       ));
320   }
321   return kTRUE;
322 }
323
324 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
325 {
326   // run the track finding
327
328    Int_t        *notr = new Int_t[fGtuParam->GetNLayers()];
329    Int_t        *ptrA = new Int_t[fGtuParam->GetNLayers()];
330    Int_t        *ptrB = new Int_t[fGtuParam->GetNLayers()];
331
332    Bool_t       *bHitA     = new Bool_t[fGtuParam->GetNLayers()];
333    Bool_t       *bHitB     = new Bool_t[fGtuParam->GetNLayers()];
334    Bool_t       *bAligned  = new Bool_t[fGtuParam->GetNLayers()];
335    Bool_t       *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
336    Bool_t       *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
337    Bool_t       *bDone     = new Bool_t[fGtuParam->GetNLayers()];
338    Bool_t        ready;
339
340    Int_t        *inc      = new Int_t[fGtuParam->GetNLayers()];
341    Int_t        *incprime = new Int_t[fGtuParam->GetNLayers()];
342
343 // ----- signals within current layer -----
344    Int_t        yPlus;
345    Int_t        yMinus;
346    Int_t        ybPlus;
347    Int_t        ybMinus;
348    Int_t        alphaPlus;
349    Int_t        alphaMinus;
350    Int_t        nHits;
351    Int_t        nUnc;
352    Int_t        nWayBeyond;
353
354    AliTRDtrackletGTU    *trkRA  = 0x0;  // reference tracklet A
355    AliTRDtrackletGTU    *trkRB  = 0x0;  // reference tracklet B
356    AliTRDtrackletGTU    *trkA   = 0x0;  // tracklet A in current layer
357    AliTRDtrackletGTU    *trkB   = 0x0;  // tracklet B in current layer
358 /*
359    AliTRDtrackletGTU    *trkEnd = new AliTRDtrackletGTU();
360    for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
361        trkEnd->SetSubChannel(i, 7);
362    trkEnd->SetYProj(0);
363    trkEnd->SetAlpha(0);
364 */
365
366    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
367      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
368      AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
369
370      ready = kFALSE; // ready if all channels done
371
372 //      for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
373 //        for (Int_t iTrkl = 0; iTrkl < fZChannelTracklets[iLayer][zch].GetEntries(); iTrkl++) {
374 //       printf("%4i/%4i(%i,%i) ",
375 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetYProj(),
376 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetAlpha(),
377 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetIndex(),
378 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetSubChannel(zch)
379 //              );
380 //        }
381 //        printf("\n");
382 //      }
383
384      for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
385        notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
386        ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
387        ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
388
389 // not necessary here
390 //       bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
391 //       bDone[layer] = (notr[layer] <= 0);
392 //       ready = ready && bDone[layer];
393
394        if (reflayer == 1)
395          AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
396      }
397
398      if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
399        continue;
400
401      while (!ready) {
402        // ----- reference tracklets -----
403        trkRA = 0x0;
404        trkRB = 0x0;
405        if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
406          trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
407        else  {
408          AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
409          break;
410        }
411
412        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
413          trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
414
415        yPlus      = trkRA->GetYProj() + fGtuParam->GetDeltaY();
416        yMinus     = trkRA->GetYProj() - fGtuParam->GetDeltaY();
417        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
418        alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
419        if (trkRB) {
420          ybPlus           = trkRB->GetYProj() + fGtuParam->GetDeltaY();
421          ybMinus          = trkRB->GetYProj() - fGtuParam->GetDeltaY();
422        }
423        else { // irrelevant (should be, is it?)
424          ybPlus           = trkRA->GetYProj() + fGtuParam->GetDeltaY();
425          ybMinus          = trkRA->GetYProj() - fGtuParam->GetDeltaY();
426        }
427
428        nHits      = 0;
429        nUnc       = 0;
430        nWayBeyond = 0;
431
432        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
433          bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
434          inc[layer] = incprime[layer] = 0;
435
436          if (layer == reflayer) {
437            bHitA[layer] = kTRUE;
438            bAligned[layer] = kTRUE;
439            nHits++;
440            continue;
441          }
442
443          trkA = 0x0;
444          trkB = 0x0;
445          if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
446            trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
447          if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
448            trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
449
450          bAlignedA[layer] = kFALSE;
451          bAlignedB[layer] = kFALSE;
452
453          if (trkA) {
454            bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
455                             !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus)  ) &&
456                             !(trkA->GetAlpha() < alphaMinus) &&
457                             !(trkA->GetAlpha() > alphaPlus) );
458            bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
459          }
460          else {
461            bHitA[layer] = 0;
462            bAlignedA[layer] = kTRUE;
463          }
464
465          if (trkB) {
466            bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
467                             !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
468                             !(trkB->GetAlpha() < alphaMinus) &&
469                             !(trkB->GetAlpha() > alphaPlus) );
470            bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
471          }
472          else {
473            bHitB[layer] = 0;
474            bAlignedB[layer] = kTRUE;
475          }
476
477          bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
478 //       bAligned[layer] = bAlignedA[layer]; //???
479
480          if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
481            nHits++;
482          else if (!bAligned[layer] )
483            nUnc++;
484          if (trkRB) {
485            if (trkA) {
486              if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
487                nWayBeyond++;
488            }
489            else
490              nWayBeyond++;
491          }
492
493          //  pre-calculation for the layer shifting (alignment w. r. t. trkRB)
494          if (trkA) {
495              if(trkRB) {
496                  if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
497                      incprime[layer] = 1;
498                  else
499                      incprime[layer] = 0;
500              }
501              else
502                  incprime[layer] = 1;
503
504              if (trkB) {
505                  if (trkRB) {
506                      if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
507                          incprime[layer] = 2;
508                  }
509                  else
510                      incprime[layer] = 2;
511              }
512          }
513        } // end of loop over layers
514
515        AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
516                        nHits, (nHits >= 4) ? "-> track found" : ""));
517
518        if (nHits >= 4) {
519          // ----- track registration -----
520          AliTRDtrackGTU *track = new AliTRDtrackGTU();
521          track->SetSector(fSector);
522          track->SetStack(fStack);
523          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
524            if (bHitA[layer] || layer == reflayer)
525              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
526            else if (bHitB[layer])
527              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
528          }
529
530          Bool_t registerTrack = kTRUE;
531          for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
532            if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
533              if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
534                AliDebug(1,"Not registered");
535                  registerTrack = kFALSE;
536              }
537            }
538          }
539          if (registerTrack) {
540              track->SetZChannel(zch);
541              track->SetRefLayerIdx(refLayerIdx);
542              fTracks[zch][refLayerIdx].Add(track);
543          }
544        }
545
546        if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
547           inc[reflayer] = 0;
548        else if (nWayBeyond > 2) // no track possible for both reference tracklets
549          inc[reflayer] = 2;
550        else
551          inc[reflayer] = 1;
552
553        if (inc[reflayer] != 0) // reflayer is shifted
554          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
555            if (layer == reflayer)
556              continue;
557            inc[layer] = incprime[layer];
558          }
559        else { // reflayer is not shifted
560          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
561            if (layer == reflayer || notr[layer] == 0)
562              continue;
563            inc[layer] = 0;
564            trkA = 0x0;
565            trkB = 0x0;
566            if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
567              trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
568
569            if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
570              trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
571
572            if (trkA) {
573              if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
574                   !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) )  // trkA could hit trkRA
575              // if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) )
576                    inc[layer] = 0;
577                else if (trkB) {
578                    if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
579                        inc[layer] = 2;
580                    else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
581                        inc[layer] = 1;
582                    else
583                        inc[layer] = incprime[layer];
584                }
585                else
586                    inc[layer] = incprime[layer];
587            }
588          }
589        }
590
591        ready = kTRUE;
592        for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
593
594          bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
595          ready = ready && bDone[layer];
596
597          if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
598            AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
599
600 //       AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
601          AliDebug(10,Form("    Layer: %i   %2i(%2i, %2i, %4i, %3i)%s%s   %2i(%2i, %2i, %4i, %3i)%s%s   +%i   %s  (no: %i)",
602                           layer,
603                           ptrA[layer],
604                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
605                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
606                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
607                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
608                           bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
609                           ptrB[layer],
610                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
611                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
612                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
613                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
614                           bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
615                           inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
616          ptrA[layer] += inc[layer];
617          ptrB[layer] += inc[layer];
618        }
619
620      } // end of while
621    } // end of loop over reflayer
622
623    delete [] bHitA;
624    delete [] bHitB;
625    delete [] bAligned;
626    delete [] bDone;
627    delete [] inc;
628    delete [] incprime;
629    delete [] bAlignedA;
630    delete [] bAlignedB;
631    delete [] notr;
632    delete [] ptrA;
633    delete [] ptrB;
634
635    return kTRUE;
636 }
637
638 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
639 {
640     TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
641     TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
642
643     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
644         tracksRefMerged[zch] = new TList;
645         tracksRefUnique[zch] = new TList;
646     }
647
648     TList *tracksZMergedStage0 = new TList;
649     TList *tracksZUniqueStage0 = new TList;
650
651     TList **tracksZSplitted = new TList*[2];
652     for (Int_t i = 0; i < 2; i++)
653         tracksZSplitted[i] = new TList;
654
655     TList *tracksZMergedStage1 = new TList;
656
657     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
658
659     Bool_t done = kFALSE;
660     Int_t minIdx = 0;
661     AliTRDtrackGTU *trkStage0 = 0x0;
662
663     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
664         // ----- Merging and Unification in Reflayers (seed_merger) -----
665         do {
666             done = kTRUE;
667             trkStage0 = 0x0;
668             for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
669                 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
670                 if (trkInRefLayer[refLayerIdx] == 0) {
671                     continue;
672                 }
673                 else if (trkStage0 == 0x0 ) {
674                     trkStage0 = trkInRefLayer[refLayerIdx];
675                     minIdx = refLayerIdx;
676                     done = kFALSE;
677                 }
678                 else if (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel() ||
679                          (trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel() && trkInRefLayer[refLayerIdx]->GetYapprox() < trkStage0->GetYapprox()) ) {
680                     minIdx = refLayerIdx;
681                     trkStage0 = trkInRefLayer[refLayerIdx];
682                     done = kFALSE;
683                 }
684             }
685             if (!trkStage0)
686               break;
687             tracksRefMerged[zch]->Add(trkStage0);
688             fTracks[zch][minIdx].RemoveFirst();
689         } while (trkStage0 != 0);
690
691         Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
692
693         AliDebug(2, Form("zchannel %i", zch));
694         TIter trackRefMerged(tracksRefMerged[zch]);
695         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
696           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
697                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
698                            trk->GetTrackletIndex(5),
699                            trk->GetTrackletIndex(4),
700                            trk->GetTrackletIndex(3),
701                            trk->GetTrackletIndex(2),
702                            trk->GetTrackletIndex(1),
703                            trk->GetTrackletIndex(0),
704                            trk->GetYapprox() >> 3,
705                            trk->GetZSubChannel()));
706         AliDebug(2, "uniquified:");
707         TIter trackRefUnique(tracksRefUnique[zch]);
708         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
709           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
710                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
711                            trk->GetTrackletIndex(5),
712                            trk->GetTrackletIndex(4),
713                            trk->GetTrackletIndex(3),
714                            trk->GetTrackletIndex(2),
715                            trk->GetTrackletIndex(1),
716                            trk->GetTrackletIndex(0),
717                            trk->GetYapprox() >> 3,
718                            trk->GetZSubChannel()));
719     }
720
721 // ----- Merging in zchannels - 1st stage -----
722
723     do {
724         done = kTRUE;
725         trkStage0 = 0x0;
726         for (Int_t zch = fGtuParam->GetNZChannels() - 1; zch > -1; zch--) {
727             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
728             if (trk == 0) {
729                 continue;
730             }
731             else if (trkStage0 == 0x0 ) {
732                 trkStage0 = trk;
733                 minIdx = zch;
734                 done = kFALSE;
735             }
736             else if ( ((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) <  ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) ||
737                       (((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) == ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
738                 minIdx = zch;
739                 trkStage0 = trk;
740                 done = kFALSE;
741             }
742         }
743
744         if (!trkStage0)
745           break;
746         tracksZMergedStage0->Add(trkStage0);
747         tracksRefUnique[minIdx]->RemoveFirst();
748     } while (trkStage0 != 0);
749
750     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
751
752     AliDebug(2, "stage 0:");
753     TIter trackZMergedStage0(tracksZMergedStage0);
754     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
755       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
756                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
757                        trk->GetTrackletIndex(5),
758                        trk->GetTrackletIndex(4),
759                        trk->GetTrackletIndex(3),
760                        trk->GetTrackletIndex(2),
761                        trk->GetTrackletIndex(1),
762                        trk->GetTrackletIndex(0),
763                        trk->GetYapprox() >> 3,
764                        trk->GetZChannel(),
765                        trk->GetZSubChannel()));
766
767     AliDebug(2, "uniquified:");
768     TIter trackZUniqueStage0(tracksZUniqueStage0);
769     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
770       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
771                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
772                        trk->GetTrackletIndex(5),
773                        trk->GetTrackletIndex(4),
774                        trk->GetTrackletIndex(3),
775                        trk->GetTrackletIndex(2),
776                        trk->GetTrackletIndex(1),
777                        trk->GetTrackletIndex(0),
778                        trk->GetYapprox() >> 3,
779                        trk->GetZChannel(),
780                        trk->GetZSubChannel()));
781
782 // ----- Splitting in z -----
783
784     TIter next(tracksZUniqueStage0);
785     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
786         tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
787     }
788
789     for (Int_t i = 0; i < 2; i++) {
790       AliDebug(2, Form("split %i:", i));
791       TIter trackZSplit(tracksZSplitted[i]);
792       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
793         AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
794                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
795                          trk->GetTrackletIndex(5),
796                          trk->GetTrackletIndex(4),
797                          trk->GetTrackletIndex(3),
798                          trk->GetTrackletIndex(2),
799                          trk->GetTrackletIndex(1),
800                          trk->GetTrackletIndex(0),
801                          trk->GetYapprox() >> 3,
802                          trk->GetZChannel(),
803                          trk->GetZSubChannel()));
804     }
805
806 // ----- Merging in zchanels - 2nd stage -----
807
808     do {
809         done = kTRUE;
810         trkStage0 = 0x0;
811         for (Int_t i = 1; i >= 0; i--) {
812             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
813             if (trk == 0) {
814                 continue;
815             }
816             else if (trkStage0 == 0x0 ) {
817                 trkStage0 = trk;
818                 minIdx = i;
819                 done = kFALSE;
820             }
821             else if ( (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) <  ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
822                       ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
823                 minIdx = i;
824                 trkStage0 = trk;
825                 done = kFALSE;
826             }
827         }
828
829         if (!trkStage0)
830           break;
831         tracksZMergedStage1->Add(trkStage0);
832         tracksZSplitted[minIdx]->RemoveFirst();
833     } while (trkStage0 != 0);
834
835     Uniquifier(tracksZMergedStage1, ListOfTracks);
836
837     AliDebug(2, "merged:");
838     TIter trackZMergedStage1(tracksZMergedStage1);
839     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
840       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
841                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
842                        trk->GetTrackletIndex(5),
843                        trk->GetTrackletIndex(4),
844                        trk->GetTrackletIndex(3),
845                        trk->GetTrackletIndex(2),
846                        trk->GetTrackletIndex(1),
847                        trk->GetTrackletIndex(0),
848                        trk->GetYapprox() >> 3,
849                        trk->GetZChannel(),
850                        trk->GetZSubChannel()));
851
852     AliDebug(2, "uniquified:");
853     TIter track(ListOfTracks);
854     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
855       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
856                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
857                        trk->GetTrackletIndex(5),
858                        trk->GetTrackletIndex(4),
859                        trk->GetTrackletIndex(3),
860                        trk->GetTrackletIndex(2),
861                        trk->GetTrackletIndex(1),
862                        trk->GetTrackletIndex(0),
863                        trk->GetYapprox() >> 3,
864                        trk->GetZChannel(),
865                        trk->GetZSubChannel()));
866
867     // cleaning up
868     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
869       delete tracksRefMerged[zch];
870       delete tracksRefUnique[zch];
871     }
872     delete [] tracksRefMerged;
873     delete [] tracksRefUnique;
874
875
876     delete tracksZMergedStage0;
877     delete tracksZUniqueStage0;
878
879     for (Int_t i = 0; i < 2; i++)
880       delete tracksZSplitted[i];
881     delete [] tracksZSplitted;
882
883     delete tracksZMergedStage1;
884
885     delete [] trkInRefLayer;
886
887     return kTRUE;
888 }
889
890 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
891 {
892   // run the track reconstruction for all tracks in the list
893
894   TIter next(ListOfTracks);
895   while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
896     CalculateTrackParams(track);
897     CalculatePID(track);
898     AliDebug(1, Form("track with pid = %i", track->GetPID()));
899   }
900   return kTRUE;
901 }
902
903 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
904 {
905   // calculate PID for the given track
906   if (!track) {
907     AliError("No track to calculate!");
908     return kFALSE;
909   }
910
911   if (AliTRDgtuParam::GetUseGTUconst()) {
912     // averaging as in GTU
913     ULong64_t coeff;
914
915     // selection of coefficient for averaging
916     switch(track->GetTrackletMask()){
917     case 0x3f:
918       // 6 tracklets
919       coeff=0x5558; // approx. 1/6
920       break;
921
922     case 0x3e:
923     case 0x3d:
924     case 0x3b:
925     case 0x37:
926     case 0x2f:
927     case 0x1f:
928       // 5 tracklets
929       coeff=0x6666; // approx. 1/5
930       break;
931
932     default:
933       // 4 tracklets
934       coeff=0x8000; // = 0.25
935     }
936     coeff &= 0x1ffff; // 17-bit constant
937
938     ULong64_t sum = 0;
939     for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
940       if ((track->GetTrackletMask() >> iLayer) & 1) {
941         sum += track->GetTracklet(iLayer)->GetPID();
942       }
943     }
944
945     sum = sum & 0x7ff;
946     ULong64_t prod   = (sum * coeff) & 0xfffffffffull;
947     ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
948
949     track->SetPID(prodFinal & 0xff);
950
951     return kTRUE;
952   }
953   else {
954
955     // simple averaging
956     Int_t nTracklets = 0;
957     Int_t pidSum = 0;
958     for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
959       if (!track->IsTrackletInLayer(layer)) {
960         continue;
961       }
962       AliTRDtrackletGTU *trk = track->GetTracklet(layer);
963       pidSum += trk->GetPID();
964       nTracklets++;
965     }
966     track->SetPID(pidSum/nTracklets);
967
968     return kTRUE;
969   }
970 }
971
972 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
973 {
974   // calculate the track parameters
975
976   if (!track) {
977     AliError("No track to calculate!");
978     return kFALSE;
979   }
980
981   Int_t a = 0;
982   Float_t b = 0;
983   Float_t c = 0;
984   Float_t x1;
985   Float_t x2;
986
987   AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
988
989   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
990     if (!track->IsTrackletInLayer(layer)) {
991       continue;
992     }
993     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
994     if (!trk) {
995       AliError(Form("Could not get tracklet in layer %i\n", layer));
996       continue;
997     }
998     AliDebug(10,Form("  layer %i trk yprime: %6i, aki: %6i", layer, trk->GetYPrime(), (Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))));
999     a += (((Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))) * trk->GetYPrime() + 1) >> 8;
1000     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1001     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1002   }
1003   if (a < 0)
1004       a += 3;
1005   a = a >> 2;
1006
1007   track->SetFitParams(a, b, c);
1008
1009   fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1010
1011   AliDebug(5,Form("  Track parameters: a = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f (trkl mask: %i)", a, b, c, x1, x2, track->GetPt(), track->GetTrackletMask()));
1012
1013   return kTRUE;
1014 }
1015
1016
1017 Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
1018 {
1019   // remove multiple occurences of the same track
1020
1021     TIter next(inlist);
1022     AliTRDtrackGTU *trkStage0 = 0x0;
1023     AliTRDtrackGTU *trkStage1 = 0x0;
1024
1025     do {
1026         trkStage0 = (AliTRDtrackGTU*) next();
1027
1028         Bool_t tracksEqual = kFALSE;
1029         if (trkStage0 != 0 && trkStage1 != 0) {
1030             for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1031                 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1032                     tracksEqual = kTRUE;
1033                     break;
1034                 }
1035             }
1036         }
1037
1038         if (tracksEqual) {
1039             if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
1040                 trkStage1 = trkStage0;
1041         }
1042         else {
1043             if (trkStage1 != 0x0)
1044                 outlist->Add(trkStage1);
1045             trkStage1 = trkStage0;
1046         }
1047
1048     } while (trkStage1 != 0x0);
1049     return kTRUE;
1050 }