fix error in setting the number of dEdx slices to be saved in ESD
[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   fZChannelTracklets(0x0),
48   fTracks(0x0),
49   fGtuParam(0x0),
50   fStack(-1),
51   fSector(-1)
52 {
53   // constructor which initializes the position information of the TMU
54
55   fGtuParam = AliTRDgtuParam::Instance();
56   fTracklets = new TObjArray*[fGtuParam->GetNLayers()];
57   fZChannelTracklets = new TList*[fGtuParam->GetNLayers()];
58   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
59     fTracklets[layer] = new TObjArray();
60     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
61   }
62   fTracks = new TList*[fGtuParam->GetNZChannels()];
63   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
64       fTracks[zch] = new TList[fGtuParam->GetNRefLayers()];
65   }
66   if (stack > -1) 
67       SetStack(stack);
68   if (sector > -1)
69       SetSector(sector);
70 }
71
72 AliTRDgtuTMU::~AliTRDgtuTMU() 
73 {
74   // destructor
75
76   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
77     delete [] fTracks[zch];
78   }
79   delete [] fTracks;
80   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
81     fTracklets[layer]->Delete(); 
82     delete [] fZChannelTracklets[layer];
83     delete fTracklets[layer];
84   }
85   delete [] fZChannelTracklets;
86   delete [] fTracklets;
87 }
88
89 Bool_t AliTRDgtuTMU::SetSector(Int_t sector)
90 {
91   // set the sector
92
93   if (sector > -1 && sector < fGtuParam->GetGeo()->Nsector() ) {
94     fSector = sector; 
95     return kTRUE;
96   }
97
98   AliError(Form("Invalid sector given: %i", sector));
99   return kFALSE;
100 }
101
102 Bool_t AliTRDgtuTMU::SetStack(Int_t stack) 
103 {
104   // Set the stack (necessary for tracking)
105
106   if (stack > -1 && stack < fGtuParam->GetGeo()->Nstack() ) {
107     fStack = stack;
108     return kTRUE;
109   }
110
111   AliError(Form("Invalid stack given: %i", stack));
112   return kFALSE;
113 }
114
115 Bool_t AliTRDgtuTMU::Reset() 
116 {
117   // delete all tracks
118
119   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
120     for (Int_t reflayeridx = 0; reflayeridx < fGtuParam->GetNRefLayers(); reflayeridx++) {
121       fTracks[zch][reflayeridx].Clear();
122     }
123   }
124
125   // delete all added tracklets
126   for (Int_t layer = 0; layer < fGtuParam->GetNLinks()/2; layer++) {
127     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++)
128       fZChannelTracklets[layer][zch].Clear();
129     fTracklets[layer]->Delete();
130   }
131
132   fSector = -1;
133   fStack = -1;
134
135   return kTRUE;
136 }
137
138 Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletBase *tracklet, Int_t link) 
139 {
140   // add a tracklet on the given link
141
142   AliTRDtrackletGTU *trkl = new AliTRDtrackletGTU(tracklet); 
143   fTracklets[(Int_t) link/2]->Add(trkl);
144   return kTRUE;
145 }
146
147
148 Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd) 
149 {
150   // performs the analysis as in a TMU module of the GTU, i. e.
151   // track matching
152   // calculation of track parameteres (pt, deflection, ???)
153
154   if (fStack < 0 || fSector < 0) {
155     AliError("No valid stack/sector set for this TMU! No tracking!");
156     return kFALSE;
157   }
158
159   // ----- Input units -----
160   AliDebug(1,"--------- Running Input units ----------");
161   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
162     if (!RunInputUnit(layer)) {
163       AliError(Form("Input unit in layer %i failed", layer));
164       return kFALSE;
165     }
166   }
167
168   // ----- Z-channel units -----
169   AliDebug(1,"--------- Running Z-channel units ----------");
170   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
171     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
172     if (!RunZChannelUnit(layer)) {
173       AliError(Form("Z-Channel unit in layer %i failed", layer));
174       return kFALSE;
175     }
176   }
177
178   // ----- track finding -----
179   AliDebug(1,"--------- Running tracking units ----------");
180   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
181     AliDebug(2,Form("Track finder for Zchannel: %i", zch));
182     if (!RunTrackFinder(zch, ListOfTracks)) {
183       AliError(Form("Track Finder in z-channel %i failed", zch));
184       return kFALSE;
185     }
186   } 
187
188   // ----- Track Merging -----
189   if (!RunTrackMerging(ListOfTracks)) {
190     AliError("Track merging failed");
191     return kFALSE;
192   }
193
194   // ----- track reconstruction -----
195   if (!RunTrackReconstruction(ListOfTracks)) {
196     AliError("Track reconstruction failed");
197     return kFALSE;
198   }
199
200   if (esd) {
201       TIter next(ListOfTracks);
202       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
203           AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
204           esd->AddTrdTrack(trdtrack);
205           delete trdtrack;
206       }
207   }
208
209   return kTRUE;
210 }
211
212 Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer) 
213 {
214   // precalculations for the tracking and reconstruction
215
216   fTracklets[layer]->Sort();
217   TIter next(fTracklets[layer]);
218
219   while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
220     trk->SetIndex(fTracklets[layer]->IndexOf(trk));
221
222     Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer); 
223     alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
224     trk->SetAlpha(alpha);
225
226     Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer)); //??? sign?
227     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
228     trk->SetYProj(yproj);
229
230     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
231
232     AliDebug(10, Form("InputUnit : GetIndex(): %3i, GetZbin(): %2i, GetY() : %5i, GetdY() : %3i, GetYPrime() : %5i, GetYProj() : %5i, GetAlpha() : %3i", 
233                       trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(), trk->GetYProj(), trk->GetAlpha() ));
234   }
235   return kTRUE;
236 }
237
238 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer) 
239 {
240   // run the z-channel unit
241
242   TIter next(fTracklets[layer]);
243
244   while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
245     AliDebug(10,Form("*TMU* Tracklet in stack %d, layer %2d: 0x%08x ", fStack, layer, trk->GetTrackletWord()));
246     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
247       if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
248         trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
249 //      printf("Z%i(%i) ", zch, trk->GetSubChannel(zch));
250
251         TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
252         AliTRDtrackletGTU *t = 0x0;
253         while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
254           if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) || 
255               (t->GetSubChannel(zch) == trk->GetSubChannel(zch) && t->GetYProj() < trk->GetYProj()) )
256             break;
257         }
258         fZChannelTracklets[layer][zch].AddAfter(t, trk);
259       }
260 //      else 
261 //        printf("      ");
262     }
263 //    printf("\n");
264   }
265   return kTRUE;
266 }
267
268 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */) 
269 {
270   // run the track finding
271
272    Int_t        *notr = new Int_t[fGtuParam->GetNLayers()];
273    Int_t        *ptrA = new Int_t[fGtuParam->GetNLayers()];
274    Int_t        *ptrB = new Int_t[fGtuParam->GetNLayers()];
275
276    Bool_t       *bHitA     = new Bool_t[fGtuParam->GetNLayers()]; 
277    Bool_t       *bHitB     = new Bool_t[fGtuParam->GetNLayers()];
278    Bool_t       *bAligned  = new Bool_t[fGtuParam->GetNLayers()];
279    Bool_t       *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
280    Bool_t       *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
281    Bool_t       *bDone     = new Bool_t[fGtuParam->GetNLayers()];
282    Bool_t        ready;
283
284    Int_t        *inc      = new Int_t[fGtuParam->GetNLayers()];
285    Int_t        *incprime = new Int_t[fGtuParam->GetNLayers()];
286
287 // ----- signals within current layer -----
288    Int_t        yPlus;
289    Int_t        yMinus;            
290    Int_t        ybPlus;    
291    Int_t        ybMinus;
292    Int_t        alphaPlus;
293    Int_t        alphaMinus; 
294    Int_t        nHits;
295    Int_t        nUnc;   
296    Int_t        nWayBeyond;
297
298    AliTRDtrackletGTU    *trkRA  = 0x0;  // reference tracklet A
299    AliTRDtrackletGTU    *trkRB  = 0x0;  // reference tracklet B
300    AliTRDtrackletGTU    *trkA   = 0x0;  // tracklet A in current layer
301    AliTRDtrackletGTU    *trkB   = 0x0;  // tracklet B in current layer
302 /*
303    AliTRDtrackletGTU    *trkEnd = new AliTRDtrackletGTU();
304    for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++) 
305        trkEnd->SetSubChannel(i, 7);
306    trkEnd->SetYProj(0);
307    trkEnd->SetAlpha(0);
308 */
309
310    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
311      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
312      AliDebug(5,Form("~~~~~ Reflayer: %i", reflayer));
313
314      ready = kFALSE; // ready if all channels done
315
316      for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
317        notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
318        ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
319        ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
320
321 // not necessary here
322 //       bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
323 //       bDone[layer] = (notr[layer] <= 0);
324 //       ready = ready && bDone[layer];
325
326        if (reflayer == 1)
327          AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
328      }
329      
330      if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0) 
331        continue;
332
333      while (!ready) {
334        // ----- reference tracklets -----
335        trkRA = 0x0;
336        trkRB = 0x0;
337        if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
338          trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
339        else  {
340          AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
341          break; 
342        }
343
344        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
345          trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
346
347        AliDebug(10,Form("ptrRA: %i, ptrRB: %i", ptrA[reflayer], ptrB[reflayer]));
348        yPlus      = trkRA->GetYProj() + fGtuParam->GetDeltaY();
349        yMinus     = trkRA->GetYProj() - fGtuParam->GetDeltaY();
350        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
351        alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
352        if (trkRB) {
353          ybPlus           = trkRB->GetYProj() + fGtuParam->GetDeltaY();
354          ybMinus          = trkRB->GetYProj() - fGtuParam->GetDeltaY();
355        }
356        else { // irrelevant (should be, is it?)
357          ybPlus           = trkRA->GetYProj() + fGtuParam->GetDeltaY();
358          ybMinus          = trkRA->GetYProj() - fGtuParam->GetDeltaY();
359        }
360
361        nHits      = 0;
362        nUnc       = 0;
363        nWayBeyond = 0;
364        
365        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
366          bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
367          inc[layer] = incprime[layer] = 0;
368
369          if (layer == reflayer) {
370            bHitA[layer] = kTRUE;
371            bAligned[layer] = kTRUE;
372            nHits++;
373            continue; 
374          }
375
376          trkA = 0x0;
377          trkB = 0x0;
378          if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
379            trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
380          if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
381            trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
382
383          bAlignedA[layer] = kFALSE;
384          bAlignedB[layer] = kFALSE;
385
386          if (trkA) { 
387            bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
388                             !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus)  ) &&
389                             !(trkA->GetAlpha() < alphaMinus) &&
390                             !(trkA->GetAlpha() > alphaPlus) );
391            bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ); 
392          }
393          else {
394            bHitA[layer] = 0;
395            bAlignedA[layer] = kTRUE;
396          }
397
398          if (trkB) {
399            bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
400                             !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
401                             !(alphaMinus > trkB->GetAlpha()) &&
402                             !(alphaPlus  > trkB->GetAlpha()) );
403            bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
404          } 
405          else {
406            bHitB[layer] = 0;
407            bAlignedB[layer] = kTRUE;
408          }
409           
410          bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
411 //       bAligned[layer] = bAlignedA[layer]; //???
412           
413          if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
414            nHits++;
415          else if (!bAligned[layer] )
416            nUnc++;
417          if (trkRB) {
418            if (trkA) {
419              if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
420                nWayBeyond++;
421            }
422            else 
423              nWayBeyond++;
424          }
425
426          //  pre-calculation for the layer shifting (alignment w. r. t. trkRB)
427          if (trkA) {
428              if(trkRB) {
429                  if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
430                      incprime[layer] = 1;
431                  else 
432                      incprime[layer] = 0;
433              }
434              else 
435                  incprime[layer] = 1;
436
437              if (trkB) { 
438                  if (trkRB) {
439                      if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
440                          incprime[layer] = 2;
441                  }
442                  else 
443                      incprime[layer] = 2;
444              }
445          }
446        } // end of loop over layers
447
448        AliDebug(5,Form("logic calculation finished, Nhits: %i", nHits));
449
450        if (nHits >= 4) {
451          // ----- track registration -----
452          AliDebug(1,"***** TMU: Track found *****");
453          AliTRDtrackGTU *track = new AliTRDtrackGTU();
454          track->SetSector(fSector);
455          track->SetStack(fStack);
456          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
457            if (bHitA[layer] || layer == reflayer) 
458              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
459            else if (bHitB[layer]) 
460              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
461          }
462
463          Bool_t registerTrack = kTRUE;
464          for (Int_t layerIdx = refLayerIdx; layerIdx > 0; layerIdx--) {
465            if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
466              if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
467                AliDebug(1,"Not registered");
468                  registerTrack = kFALSE;
469              }
470            }
471          }
472          if (registerTrack) {
473              track->SetZChannel(zch);
474              track->SetRefLayerIdx(refLayerIdx);
475              fTracks[zch][refLayerIdx].Add(track);
476          }
477        }
478         
479        if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
480           inc[reflayer] = 0;
481        else if (nWayBeyond > 2) // no track possible for both reference tracklets
482          inc[reflayer] = 2;
483        else 
484          inc[reflayer] = 1;
485         
486        if (inc[reflayer] != 0) // reflayer is shifted
487          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
488            if (layer == reflayer)
489              continue;
490            inc[layer] = incprime[layer]; 
491          }
492        else { // reflayer is not shifted
493          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
494            if (layer == reflayer || notr[layer] == 0)
495              continue;
496            inc[layer] = 0;
497            trkA = 0x0;
498            trkB = 0x0;
499            if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
500              trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
501
502            if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
503              trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
504
505            if (trkA) {
506                if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) && 
507                     !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) )  // trkA could hit trkRA
508                    inc[layer] = 0;
509                else if (trkB) {
510                    if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
511                        inc[layer] = 2;
512                    else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
513                        inc[layer] = 1;
514                    else 
515                        inc[layer] = incprime[layer];
516                }
517                else 
518                    inc[layer] = incprime[layer];
519            }
520          }
521        }
522        
523        ready = kTRUE;
524        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
525
526          bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
527          ready = ready && bDone[layer];
528
529          if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
530            AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
531
532 //       AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
533          AliDebug(10,Form(" -- Layer: %i   %i   %i   +%i   %s  (no: %i)", layer, ptrA[layer], ptrB[layer], inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
534          ptrA[layer] += inc[layer];
535          ptrB[layer] += inc[layer];
536        }
537
538      } // end of while
539    } // end of loop over reflayer
540
541    delete [] bHitA;
542    delete [] bHitB;
543    delete [] bAligned;
544    delete [] bDone;
545    delete [] inc;
546    delete [] incprime;
547    delete [] bAlignedA;
548    delete [] bAlignedB;
549    delete [] notr;
550    delete [] ptrA;
551    delete [] ptrB;
552
553    return kTRUE;
554 }
555
556 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks) 
557 {
558     TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
559     TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
560
561     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
562         tracksRefMerged[zch] = new TList;
563         tracksRefUnique[zch] = new TList;
564     }
565
566     TList *tracksZMergedStage0 = new TList;
567     TList *tracksZUniqueStage0 = new TList;
568
569     TList **tracksZSplitted = new TList*[2];
570     for (Int_t i = 0; i < 2; i++)
571         tracksZSplitted[i] = new TList;
572
573     TList *tracksZMergedStage1 = new TList;
574
575     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
576
577     Bool_t done = kFALSE;
578     Int_t minIdx = 0;
579     AliTRDtrackGTU *trkStage0 = 0x0;
580
581     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
582         // ----- Merging and Unification in Reflayers -----
583         do {
584             done = kTRUE;
585             trkStage0 = 0x0;
586             for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
587                 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
588                 if (trkInRefLayer[refLayerIdx] == 0) {
589                     continue;
590                 }
591                 else if (trkStage0 == 0x0 ) {
592                     trkStage0 = trkInRefLayer[refLayerIdx];
593                     minIdx = refLayerIdx;
594                     done = kFALSE;
595                 }
596                 else if (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel() || 
597                          (trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel() && trkInRefLayer[refLayerIdx]->GetYapprox() < trkStage0->GetYapprox()) ) {
598                     minIdx = refLayerIdx;
599                     trkStage0 = trkInRefLayer[refLayerIdx];
600                     done = kFALSE;
601                 }
602             }
603             if (!trkStage0)
604               break;
605             tracksRefMerged[zch]->Add(trkStage0);
606             fTracks[zch][minIdx].RemoveFirst();
607         } while (trkStage0 != 0);
608
609         Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
610     }
611
612 // ----- Merging in zchannels - 1st stage -----
613     
614     do {
615         done = kTRUE;
616         trkStage0 = 0x0;
617         for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
618             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
619             if (trk == 0) {
620                 continue;
621             }
622             else if (trkStage0 == 0x0 ) {
623                 trkStage0 = trk;
624                 minIdx = zch;
625                 done = kFALSE;
626             }
627             else if ( ((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) < ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) ) {
628                 minIdx = zch;
629                 trkStage0 = trk;
630                 done = kFALSE;
631             }
632         }
633         
634         if (!trkStage0)
635           break;
636         tracksZMergedStage0->Add(trkStage0);
637         tracksRefUnique[minIdx]->RemoveFirst();
638     } while (trkStage0 != 0);
639     
640     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
641     
642 // ----- Splitting in z -----
643     
644     TIter next(tracksZUniqueStage0);
645     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
646         tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
647     }
648     
649 // ----- Merging in zchanels - 2nd stage -----
650     
651     do {
652         done = kTRUE;
653         trkStage0 = 0x0;
654         for (Int_t i = 0; i < 2; i++) {
655             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
656             if (trk == 0) {
657                 continue;
658             }
659             else if (trkStage0 == 0x0 ) {
660                 trkStage0 = trk;
661                 minIdx = i;
662                 done = kFALSE;
663             }
664             else if ( ((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) < ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2) ) {
665                 minIdx = i;
666                 trkStage0 = trk;
667                 done = kFALSE;
668             }
669         }
670         
671         if (!trkStage0)
672           break;
673         tracksZMergedStage1->Add(trkStage0);
674         tracksZSplitted[minIdx]->RemoveFirst();
675     } while (trkStage0 != 0);
676     
677     Uniquifier(tracksZMergedStage1, ListOfTracks);
678     
679     // cleaning up                                                                                                                         
680     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
681       delete tracksRefMerged[zch];
682       delete tracksRefUnique[zch];
683     }
684     delete [] tracksRefMerged;
685     delete [] tracksRefUnique;
686
687
688     delete tracksZMergedStage0;
689     delete tracksZUniqueStage0;
690
691     for (Int_t i = 0; i < 2; i++)
692       delete tracksZSplitted[i];
693     delete [] tracksZSplitted;
694
695     delete tracksZMergedStage1;
696
697     delete [] trkInRefLayer;
698
699     return kTRUE;
700 }
701
702 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks) 
703 {
704   // run the track reconstruction for all tracks in the list
705
706   TIter next(ListOfTracks);
707   while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
708     CalculateTrackParams(track);
709     CalculatePID(track);
710   }
711   return kTRUE;
712 }
713
714 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
715 {
716   // calculate PID for the given track
717   if (!track) {
718     AliError("No track to calculate!");
719     return kFALSE;
720   }
721
722   Int_t nTracklets = 0;
723   Int_t pidSum = 0;
724   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
725     if (!track->IsTrackletInLayer(layer)) {
726       continue;
727     }
728     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
729     pidSum += trk->GetPID();
730     nTracklets++;
731   }
732   track->SetPID(pidSum/nTracklets);
733   return kTRUE;
734 }
735
736 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track) 
737 {
738   // calculate the track parameters
739
740   if (!track) {
741     AliError("No track to calculate!");
742     return kFALSE;
743   }
744
745   Int_t a = 0;
746   Float_t b = 0;
747   Float_t c = 0;
748   Float_t x1;
749   Float_t x2;
750
751   AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
752
753   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
754     if (!track->IsTrackletInLayer(layer)) {
755       continue;
756     }
757     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
758     if (!trk) {
759       AliError(Form("Could not get tracklet in layer %i\n", layer));
760       continue;
761     }
762     AliDebug(10,Form("trk yprime: %i", trk->GetYPrime()));
763     a += (((Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))) * trk->GetYPrime() + 1) >> 8; 
764     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
765     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
766   }
767   if (a < 0)
768       a += 3;
769   a = a >> 2;
770
771   fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
772   AliDebug(10,Form("Intersection points: %f, %f", x1, x2));
773   AliDebug(10,Form("Sum: a = %5i, b = %9.2f, c = %9.2f\n", a, b, c));
774   track->SetFitParams(a, b, c);
775
776   Float_t r = fGtuParam->GetPt(a, b, x1, x2);
777   Int_t pt = (Int_t) (2 * r);
778   if (pt >= 0) 
779       pt += 32;
780   else
781       pt -= 29;
782   pt /= 2;
783   track->SetPtInt(pt);
784   AliDebug(5,Form("Track parameters: a = %i, b = %f, c = %f, x1 = %f, x2 = %f, r = %f, pt = %f (trkl mask: %i)", a, b, c, x1, x2, r, track->GetPt(), track->GetTrackletMask()));
785   return kTRUE;
786 }
787
788 Bool_t AliTRDgtuTMU::WriteTrackletsToTree(TTree *trklTree)
789 {
790   if (!trklTree) {
791     AliError("No tree given");
792     return kFALSE;
793   }
794   AliTRDtrackletGTU *trkl = 0x0;
795   TBranch *branch = trklTree->GetBranch("gtutracklets");
796   if (!branch) {
797       branch = trklTree->Branch("gtutracklets", "AliTRDtrackletGTU", &trkl, 32000, 99);
798   }
799
800   AliDebug(5,Form("---------- Writing tracklets to tree (not yet) ----------"));
801   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
802     TIter next(fTracklets[layer]);
803     while ((trkl = (AliTRDtrackletGTU*) next())) {
804         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) ));
805         branch->SetAddress(&trkl);
806         trklTree->Fill();
807     }
808   }
809   return kTRUE;
810 }
811
812 Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
813 {
814   // remove multiple occurences of the same track
815
816     TIter next(inlist);
817     AliTRDtrackGTU *trkStage0 = 0x0;
818     AliTRDtrackGTU *trkStage1 = 0x0;
819
820     do {
821         trkStage0 = (AliTRDtrackGTU*) next();
822
823         Bool_t tracksEqual = kFALSE;
824         if (trkStage0 != 0 && trkStage1 != 0) {
825             for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
826                 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
827                     tracksEqual = kTRUE;
828                     break;
829                 }
830             }
831         }
832
833         if (tracksEqual) {
834             if (trkStage0->GetNTracklets() > trkStage1->GetNTracklets())
835                 trkStage1 = trkStage0;
836         } 
837         else {
838             if (trkStage1 != 0x0) 
839                 outlist->Add(trkStage1);
840             trkStage1 = trkStage0;
841         } 
842         
843     } while (trkStage1 != 0x0);
844     return kTRUE;
845 }