- extend debug output for GTU simulation
[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     trk->SetAlpha(alpha);
268
269     Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer)); //??? sign?
270     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
271     trk->SetYProj(yproj);
272
273     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
274
275     AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i",
276                       trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
277                       trk->GetYProj(), trk->GetAlpha(), trk->GetPID() ));
278   }
279   return kTRUE;
280 }
281
282 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
283 {
284   // run the z-channel unit
285
286   TIter next(fTrackletsPostInput[layer]);
287
288   while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
289     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
290       if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
291         trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
292
293         TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
294         AliTRDtrackletGTU *t = 0x0;
295         while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
296           if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
297               ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
298             break;
299           }
300         }
301         if (t)
302           fZChannelTracklets[layer][zch].AddAfter(t, trk);
303         else
304           fZChannelTracklets[layer][zch].AddFirst(trk);
305       }
306     }
307     AliDebug(10, Form("stack %d, layer %2d: 0x%08x Z(2..0)=%i/%i/%i",
308                       fStack, layer, trk->GetTrackletWord(),
309                       fGtuParam->IsInZChannel(fStack, layer, 2, trk->GetZbin()) ? trk->GetSubChannel(2) : -1,
310                       fGtuParam->IsInZChannel(fStack, layer, 1, trk->GetZbin()) ? trk->GetSubChannel(1) : -1,
311                       fGtuParam->IsInZChannel(fStack, layer, 0, trk->GetZbin()) ? trk->GetSubChannel(0) : -1
312                       ));
313   }
314   return kTRUE;
315 }
316
317 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
318 {
319   // run the track finding
320
321    Int_t        *notr = new Int_t[fGtuParam->GetNLayers()];
322    Int_t        *ptrA = new Int_t[fGtuParam->GetNLayers()];
323    Int_t        *ptrB = new Int_t[fGtuParam->GetNLayers()];
324
325    Bool_t       *bHitA     = new Bool_t[fGtuParam->GetNLayers()];
326    Bool_t       *bHitB     = new Bool_t[fGtuParam->GetNLayers()];
327    Bool_t       *bAligned  = new Bool_t[fGtuParam->GetNLayers()];
328    Bool_t       *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
329    Bool_t       *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
330    Bool_t       *bDone     = new Bool_t[fGtuParam->GetNLayers()];
331    Bool_t        ready;
332
333    Int_t        *inc      = new Int_t[fGtuParam->GetNLayers()];
334    Int_t        *incprime = new Int_t[fGtuParam->GetNLayers()];
335
336 // ----- signals within current layer -----
337    Int_t        yPlus;
338    Int_t        yMinus;
339    Int_t        ybPlus;
340    Int_t        ybMinus;
341    Int_t        alphaPlus;
342    Int_t        alphaMinus;
343    Int_t        nHits;
344    Int_t        nUnc;
345    Int_t        nWayBeyond;
346
347    AliTRDtrackletGTU    *trkRA  = 0x0;  // reference tracklet A
348    AliTRDtrackletGTU    *trkRB  = 0x0;  // reference tracklet B
349    AliTRDtrackletGTU    *trkA   = 0x0;  // tracklet A in current layer
350    AliTRDtrackletGTU    *trkB   = 0x0;  // tracklet B in current layer
351 /*
352    AliTRDtrackletGTU    *trkEnd = new AliTRDtrackletGTU();
353    for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
354        trkEnd->SetSubChannel(i, 7);
355    trkEnd->SetYProj(0);
356    trkEnd->SetAlpha(0);
357 */
358
359    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
360      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
361      AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
362
363      ready = kFALSE; // ready if all channels done
364
365 //      for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
366 //        for (Int_t iTrkl = 0; iTrkl < fZChannelTracklets[iLayer][zch].GetEntries(); iTrkl++) {
367 //       printf("%4i/%4i(%i,%i) ",
368 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetYProj(),
369 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetAlpha(),
370 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetIndex(),
371 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetSubChannel(zch)
372 //              );
373 //        }
374 //        printf("\n");
375 //      }
376
377      for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
378        notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
379        ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
380        ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
381
382 // not necessary here
383 //       bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
384 //       bDone[layer] = (notr[layer] <= 0);
385 //       ready = ready && bDone[layer];
386
387        if (reflayer == 1)
388          AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
389      }
390
391      if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
392        continue;
393
394      while (!ready) {
395        // ----- reference tracklets -----
396        trkRA = 0x0;
397        trkRB = 0x0;
398        if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
399          trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
400        else  {
401          AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
402          break;
403        }
404
405        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
406          trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
407
408        yPlus      = trkRA->GetYProj() + fGtuParam->GetDeltaY();
409        yMinus     = trkRA->GetYProj() - fGtuParam->GetDeltaY();
410        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
411        alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
412        if (trkRB) {
413          ybPlus           = trkRB->GetYProj() + fGtuParam->GetDeltaY();
414          ybMinus          = trkRB->GetYProj() - fGtuParam->GetDeltaY();
415        }
416        else { // irrelevant (should be, is it?)
417          ybPlus           = trkRA->GetYProj() + fGtuParam->GetDeltaY();
418          ybMinus          = trkRA->GetYProj() - fGtuParam->GetDeltaY();
419        }
420
421        nHits      = 0;
422        nUnc       = 0;
423        nWayBeyond = 0;
424
425        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
426          bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
427          inc[layer] = incprime[layer] = 0;
428
429          if (layer == reflayer) {
430            bHitA[layer] = kTRUE;
431            bAligned[layer] = kTRUE;
432            nHits++;
433            continue;
434          }
435
436          trkA = 0x0;
437          trkB = 0x0;
438          if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
439            trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
440          if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
441            trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
442
443          bAlignedA[layer] = kFALSE;
444          bAlignedB[layer] = kFALSE;
445
446          if (trkA) {
447            bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
448                             !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus)  ) &&
449                             !(trkA->GetAlpha() < alphaMinus) &&
450                             !(trkA->GetAlpha() > alphaPlus) );
451            bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
452          }
453          else {
454            bHitA[layer] = 0;
455            bAlignedA[layer] = kTRUE;
456          }
457
458          if (trkB) {
459            bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
460                             !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
461                             !(trkB->GetAlpha() < alphaMinus) &&
462                             !(trkB->GetAlpha() > alphaPlus) );
463            bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
464          }
465          else {
466            bHitB[layer] = 0;
467            bAlignedB[layer] = kTRUE;
468          }
469
470          bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
471 //       bAligned[layer] = bAlignedA[layer]; //???
472
473          if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
474            nHits++;
475          else if (!bAligned[layer] )
476            nUnc++;
477          if (trkRB) {
478            if (trkA) {
479              if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
480                nWayBeyond++;
481            }
482            else
483              nWayBeyond++;
484          }
485
486          //  pre-calculation for the layer shifting (alignment w. r. t. trkRB)
487          if (trkA) {
488              if(trkRB) {
489                  if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
490                      incprime[layer] = 1;
491                  else
492                      incprime[layer] = 0;
493              }
494              else
495                  incprime[layer] = 1;
496
497              if (trkB) {
498                  if (trkRB) {
499                      if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
500                          incprime[layer] = 2;
501                  }
502                  else
503                      incprime[layer] = 2;
504              }
505          }
506        } // end of loop over layers
507
508        AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
509                        nHits, (nHits >= 4) ? "-> track found" : ""));
510
511        if (nHits >= 4) {
512          // ----- track registration -----
513          AliTRDtrackGTU *track = new AliTRDtrackGTU();
514          track->SetSector(fSector);
515          track->SetStack(fStack);
516          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
517            if (bHitA[layer] || layer == reflayer)
518              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
519            else if (bHitB[layer])
520              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
521          }
522
523          Bool_t registerTrack = kTRUE;
524          for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
525            if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
526              if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
527                AliDebug(1,"Not registered");
528                  registerTrack = kFALSE;
529              }
530            }
531          }
532          if (registerTrack) {
533              track->SetZChannel(zch);
534              track->SetRefLayerIdx(refLayerIdx);
535              fTracks[zch][refLayerIdx].Add(track);
536          }
537        }
538
539        if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
540           inc[reflayer] = 0;
541        else if (nWayBeyond > 2) // no track possible for both reference tracklets
542          inc[reflayer] = 2;
543        else
544          inc[reflayer] = 1;
545
546        if (inc[reflayer] != 0) // reflayer is shifted
547          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
548            if (layer == reflayer)
549              continue;
550            inc[layer] = incprime[layer];
551          }
552        else { // reflayer is not shifted
553          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
554            if (layer == reflayer || notr[layer] == 0)
555              continue;
556            inc[layer] = 0;
557            trkA = 0x0;
558            trkB = 0x0;
559            if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
560              trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
561
562            if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
563              trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
564
565            if (trkA) {
566              if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
567                   !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) )  // trkA could hit trkRA
568              // if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) )
569                    inc[layer] = 0;
570                else if (trkB) {
571                    if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
572                        inc[layer] = 2;
573                    else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
574                        inc[layer] = 1;
575                    else
576                        inc[layer] = incprime[layer];
577                }
578                else
579                    inc[layer] = incprime[layer];
580            }
581          }
582        }
583
584        ready = kTRUE;
585        for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
586
587          bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
588          ready = ready && bDone[layer];
589
590          if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
591            AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
592
593 //       AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
594          AliDebug(10,Form("    Layer: %i   %2i(%2i, %2i, %4i, %3i)%s%s   %2i(%2i, %2i, %4i, %3i)%s%s   +%i   %s  (no: %i)",
595                           layer,
596                           ptrA[layer],
597                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
598                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
599                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
600                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
601                           bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
602                           ptrB[layer],
603                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
604                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
605                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
606                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
607                           bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
608                           inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
609          ptrA[layer] += inc[layer];
610          ptrB[layer] += inc[layer];
611        }
612
613      } // end of while
614    } // end of loop over reflayer
615
616    delete [] bHitA;
617    delete [] bHitB;
618    delete [] bAligned;
619    delete [] bDone;
620    delete [] inc;
621    delete [] incprime;
622    delete [] bAlignedA;
623    delete [] bAlignedB;
624    delete [] notr;
625    delete [] ptrA;
626    delete [] ptrB;
627
628    return kTRUE;
629 }
630
631 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
632 {
633     TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
634     TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
635
636     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
637         tracksRefMerged[zch] = new TList;
638         tracksRefUnique[zch] = new TList;
639     }
640
641     TList *tracksZMergedStage0 = new TList;
642     TList *tracksZUniqueStage0 = new TList;
643
644     TList **tracksZSplitted = new TList*[2];
645     for (Int_t i = 0; i < 2; i++)
646         tracksZSplitted[i] = new TList;
647
648     TList *tracksZMergedStage1 = new TList;
649
650     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
651
652     Bool_t done = kFALSE;
653     Int_t minIdx = 0;
654     AliTRDtrackGTU *trkStage0 = 0x0;
655
656     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
657         // ----- Merging and Unification in Reflayers (seed_merger) -----
658         do {
659             done = kTRUE;
660             trkStage0 = 0x0;
661             for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
662                 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
663                 if (trkInRefLayer[refLayerIdx] == 0) {
664                     continue;
665                 }
666                 else if (trkStage0 == 0x0 ) {
667                     trkStage0 = trkInRefLayer[refLayerIdx];
668                     minIdx = refLayerIdx;
669                     done = kFALSE;
670                 }
671                 else if (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel() ||
672                          (trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel() && trkInRefLayer[refLayerIdx]->GetYapprox() < trkStage0->GetYapprox()) ) {
673                     minIdx = refLayerIdx;
674                     trkStage0 = trkInRefLayer[refLayerIdx];
675                     done = kFALSE;
676                 }
677             }
678             if (!trkStage0)
679               break;
680             tracksRefMerged[zch]->Add(trkStage0);
681             fTracks[zch][minIdx].RemoveFirst();
682         } while (trkStage0 != 0);
683
684         Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
685
686         AliDebug(2, Form("zchannel %i", zch));
687         TIter trackRefMerged(tracksRefMerged[zch]);
688         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
689           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
690                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
691                            trk->GetTrackletIndex(5),
692                            trk->GetTrackletIndex(4),
693                            trk->GetTrackletIndex(3),
694                            trk->GetTrackletIndex(2),
695                            trk->GetTrackletIndex(1),
696                            trk->GetTrackletIndex(0),
697                            trk->GetYapprox() >> 3,
698                            trk->GetZSubChannel()));
699         AliDebug(2, "uniquified:");
700         TIter trackRefUnique(tracksRefUnique[zch]);
701         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
702           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
703                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
704                            trk->GetTrackletIndex(5),
705                            trk->GetTrackletIndex(4),
706                            trk->GetTrackletIndex(3),
707                            trk->GetTrackletIndex(2),
708                            trk->GetTrackletIndex(1),
709                            trk->GetTrackletIndex(0),
710                            trk->GetYapprox() >> 3,
711                            trk->GetZSubChannel()));
712     }
713
714 // ----- Merging in zchannels - 1st stage -----
715
716     do {
717         done = kTRUE;
718         trkStage0 = 0x0;
719         for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
720             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
721             if (trk == 0) {
722                 continue;
723             }
724             else if (trkStage0 == 0x0 ) {
725                 trkStage0 = trk;
726                 minIdx = zch;
727                 done = kFALSE;
728             }
729             else if ( ((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) <  ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) ||
730                       (((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) == ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
731                 minIdx = zch;
732                 trkStage0 = trk;
733                 done = kFALSE;
734             }
735         }
736
737         if (!trkStage0)
738           break;
739         tracksZMergedStage0->Add(trkStage0);
740         tracksRefUnique[minIdx]->RemoveFirst();
741     } while (trkStage0 != 0);
742
743     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
744
745     AliDebug(2, "stage 0:");
746     TIter trackZMergedStage0(tracksZMergedStage0);
747     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
748       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
749                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
750                        trk->GetTrackletIndex(5),
751                        trk->GetTrackletIndex(4),
752                        trk->GetTrackletIndex(3),
753                        trk->GetTrackletIndex(2),
754                        trk->GetTrackletIndex(1),
755                        trk->GetTrackletIndex(0),
756                        trk->GetYapprox() >> 3,
757                        trk->GetZChannel(),
758                        trk->GetZSubChannel()));
759
760     AliDebug(2, "uniquified:");
761     TIter trackZUniqueStage0(tracksZUniqueStage0);
762     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
763       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
764                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
765                        trk->GetTrackletIndex(5),
766                        trk->GetTrackletIndex(4),
767                        trk->GetTrackletIndex(3),
768                        trk->GetTrackletIndex(2),
769                        trk->GetTrackletIndex(1),
770                        trk->GetTrackletIndex(0),
771                        trk->GetYapprox() >> 3,
772                        trk->GetZChannel(),
773                        trk->GetZSubChannel()));
774
775 // ----- Splitting in z -----
776
777     TIter next(tracksZUniqueStage0);
778     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
779         tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
780     }
781
782     for (Int_t i = 0; i < 2; i++) {
783       AliDebug(2, Form("split %i:", i));
784       TIter trackZSplit(tracksZSplitted[i]);
785       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
786         AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
787                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
788                          trk->GetTrackletIndex(5),
789                          trk->GetTrackletIndex(4),
790                          trk->GetTrackletIndex(3),
791                          trk->GetTrackletIndex(2),
792                          trk->GetTrackletIndex(1),
793                          trk->GetTrackletIndex(0),
794                          trk->GetYapprox() >> 3,
795                          trk->GetZChannel(),
796                          trk->GetZSubChannel()));
797     }
798
799 // ----- Merging in zchanels - 2nd stage -----
800
801     do {
802         done = kTRUE;
803         trkStage0 = 0x0;
804         for (Int_t i = 1; i >= 0; i--) {
805             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
806             if (trk == 0) {
807                 continue;
808             }
809             else if (trkStage0 == 0x0 ) {
810                 trkStage0 = trk;
811                 minIdx = i;
812                 done = kFALSE;
813             }
814             else if ( (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) <  ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
815                       ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
816                 minIdx = i;
817                 trkStage0 = trk;
818                 done = kFALSE;
819             }
820         }
821
822         if (!trkStage0)
823           break;
824         tracksZMergedStage1->Add(trkStage0);
825         tracksZSplitted[minIdx]->RemoveFirst();
826     } while (trkStage0 != 0);
827
828     Uniquifier(tracksZMergedStage1, ListOfTracks);
829
830     AliDebug(2, "merged:");
831     TIter trackZMergedStage1(tracksZMergedStage1);
832     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
833       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
834                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
835                        trk->GetTrackletIndex(5),
836                        trk->GetTrackletIndex(4),
837                        trk->GetTrackletIndex(3),
838                        trk->GetTrackletIndex(2),
839                        trk->GetTrackletIndex(1),
840                        trk->GetTrackletIndex(0),
841                        trk->GetYapprox() >> 3,
842                        trk->GetZChannel(),
843                        trk->GetZSubChannel()));
844
845     AliDebug(2, "uniquified:");
846     TIter track(ListOfTracks);
847     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
848       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
849                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
850                        trk->GetTrackletIndex(5),
851                        trk->GetTrackletIndex(4),
852                        trk->GetTrackletIndex(3),
853                        trk->GetTrackletIndex(2),
854                        trk->GetTrackletIndex(1),
855                        trk->GetTrackletIndex(0),
856                        trk->GetYapprox() >> 3,
857                        trk->GetZChannel(),
858                        trk->GetZSubChannel()));
859
860     // cleaning up
861     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
862       delete tracksRefMerged[zch];
863       delete tracksRefUnique[zch];
864     }
865     delete [] tracksRefMerged;
866     delete [] tracksRefUnique;
867
868
869     delete tracksZMergedStage0;
870     delete tracksZUniqueStage0;
871
872     for (Int_t i = 0; i < 2; i++)
873       delete tracksZSplitted[i];
874     delete [] tracksZSplitted;
875
876     delete tracksZMergedStage1;
877
878     delete [] trkInRefLayer;
879
880     return kTRUE;
881 }
882
883 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
884 {
885   // run the track reconstruction for all tracks in the list
886
887   TIter next(ListOfTracks);
888   while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
889     CalculateTrackParams(track);
890     CalculatePID(track);
891     AliDebug(1, Form("track with pid = %i", track->GetPID()));
892   }
893   return kTRUE;
894 }
895
896 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
897 {
898   // calculate PID for the given track
899   if (!track) {
900     AliError("No track to calculate!");
901     return kFALSE;
902   }
903
904   if (AliTRDgtuParam::GetUseGTUconst()) {
905     // averaging as in GTU
906     ULong64_t coeff;
907
908     // selection of coefficient for averaging
909     switch(track->GetTrackletMask()){
910     case 0x3f:
911       // 6 tracklets
912       coeff=0x5558; // approx. 1/6
913       break;
914
915     case 0x3e:
916     case 0x3d:
917     case 0x3b:
918     case 0x37:
919     case 0x2f:
920     case 0x1f:
921       // 5 tracklets
922       coeff=0x6666; // approx. 1/5
923       break;
924
925     default:
926       // 4 tracklets
927       coeff=0x8000; // = 0.25
928     }
929     coeff &= 0x1ffff; // 17-bit constant
930
931     ULong64_t sum = 0;
932     for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
933       if ((track->GetTrackletMask() >> iLayer) & 1) {
934         sum += track->GetTracklet(iLayer)->GetPID();
935       }
936     }
937
938     ULong64_t sumExt = (sum << 1) & 0xfff; // 11 bit -> 12 bit vector
939     ULong64_t prod   = (sumExt * coeff) & 0xfffffffffull; // 18x18 signed -> 36
940     ULong64_t prodFinal = ((prod >> 18) + ((prod >> 17) & 1)) & 0xff; // rounding term is equivalent to adding 5 to sum_ext
941
942     track->SetPID(prodFinal & 0xff);
943
944     return kTRUE;
945   }
946   else {
947
948     // simple averaging
949     Int_t nTracklets = 0;
950     Int_t pidSum = 0;
951     for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
952       if (!track->IsTrackletInLayer(layer)) {
953         continue;
954       }
955       AliTRDtrackletGTU *trk = track->GetTracklet(layer);
956       pidSum += trk->GetPID();
957       nTracklets++;
958     }
959     track->SetPID(pidSum/nTracklets);
960
961     return kTRUE;
962   }
963 }
964
965 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
966 {
967   // calculate the track parameters
968
969   if (!track) {
970     AliError("No track to calculate!");
971     return kFALSE;
972   }
973
974   Int_t a = 0;
975   Float_t b = 0;
976   Float_t c = 0;
977   Float_t x1;
978   Float_t x2;
979
980   AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
981
982   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
983     if (!track->IsTrackletInLayer(layer)) {
984       continue;
985     }
986     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
987     if (!trk) {
988       AliError(Form("Could not get tracklet in layer %i\n", layer));
989       continue;
990     }
991     AliDebug(10,Form("  trk yprime: %6i, aki: %6i", trk->GetYPrime(), (Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))));
992     a += (((Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))) * trk->GetYPrime() + 1) >> 8;
993     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
994     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
995   }
996   if (a < 0)
997       a += 3;
998   a = a >> 2;
999
1000   track->SetFitParams(a, b, c);
1001
1002   fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1003
1004   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()));
1005
1006   return kTRUE;
1007 }
1008
1009 Bool_t AliTRDgtuTMU::WriteTrackletsToTree(TTree *trklTree)
1010 {
1011   if (!trklTree) {
1012     AliError("No tree given");
1013     return kFALSE;
1014   }
1015   AliTRDtrackletGTU *trkl = 0x0;
1016   TBranch *branch = trklTree->GetBranch("gtutracklets");
1017   if (!branch) {
1018       branch = trklTree->Branch("gtutracklets", "AliTRDtrackletGTU", &trkl, 32000, 99);
1019   }
1020
1021   AliDebug(5,Form("---------- Writing tracklets to tree (not yet) ----------"));
1022   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1023     TIter next(fTrackletsPostInput[layer]);
1024     while ((trkl = (AliTRDtrackletGTU*) next())) {
1025         AliDebug(10,Form("InputUnit : GetIndex(): %3i, GetZbin(): %2i, GetY() : %5i, GetdY() : %3i, GetYPrime() : %5i, GetYProj() : %5i, GetAlpha() : %3i, Zidx(2..0): %i  %i  %i", trkl->GetIndex(), trkl->GetZbin(), trkl->GetYbin(), trkl->GetdY(), trkl->GetYPrime(), trkl->GetYProj(), trkl->GetAlpha(), trkl->GetSubChannel(2), trkl->GetSubChannel(1), trkl->GetSubChannel(0) ));
1026         branch->SetAddress(&trkl);
1027         trklTree->Fill();
1028     }
1029   }
1030   return kTRUE;
1031 }
1032
1033 Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
1034 {
1035   // remove multiple occurences of the same track
1036
1037     TIter next(inlist);
1038     AliTRDtrackGTU *trkStage0 = 0x0;
1039     AliTRDtrackGTU *trkStage1 = 0x0;
1040
1041     do {
1042         trkStage0 = (AliTRDtrackGTU*) next();
1043
1044         Bool_t tracksEqual = kFALSE;
1045         if (trkStage0 != 0 && trkStage1 != 0) {
1046             for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1047                 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1048                     tracksEqual = kTRUE;
1049                     break;
1050                 }
1051             }
1052         }
1053
1054         if (tracksEqual) {
1055             if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
1056                 trkStage1 = trkStage0;
1057         }
1058         else {
1059             if (trkStage1 != 0x0)
1060                 outlist->Add(trkStage1);
1061             trkStage1 = trkStage0;
1062         }
1063
1064     } while (trkStage1 != 0x0);
1065     return kTRUE;
1066 }