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