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