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, Int_t outLabel)
 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     if (!RunZChannelUnit(layer)) {
 189       AliError(Form("Z-Channel unit in layer %i failed", layer));
 190       return kFALSE;
 191     }
 192   }
 193
 194   // ----- track finding -----
 195   AliDebug(1,"--------- Running tracking units ----------");
 196   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 197     if (!RunTrackFinder(zch, ListOfTracks)) {
 198       AliError(Form("Track Finder in z-channel %i failed", zch));
 199       return kFALSE;
 200     }
 201   }
 202
 203   // ----- Track Merging -----
 204   if (!RunTrackMerging(ListOfTracks)) {
 205     AliError("Track merging failed");
 206     return kFALSE;
 207   }
 208
 209   // ----- track reconstruction -----
 210   if (!RunTrackReconstruction(ListOfTracks)) {
 211     AliError("Track reconstruction failed");
 212     return kFALSE;
 213   }
 214
 215   // ----- label calculation and ESD storage -----
 216   TIter next(ListOfTracks);
 217   while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
 218     if (outLabel == -1)
 219       trk->CookLabel();
 220     else
 221       trk->SetLabel(outLabel);
 222     if (esd) {
 223       AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
 224       esd->AddTrdTrack(trdtrack);
 225       delete trdtrack;
 226     }
 227   }
 228
 229   return kTRUE;
 230 }
 231
 232 Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
 233 {
 234   // precalculations for the tracking and reconstruction
 235
 236   Int_t iTrkl0 = 0; // A-side tracklet
 237   Int_t iTrkl1 = 0; // B-side tracklet
 238
 239   Int_t nTracklets = 0;
 240   while ((!fGtuParam->GetLimitNoTracklets() ||
 241           (nTracklets < fGtuParam->GetMaxNoTracklets())) &&
 242          ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) ||
 243           (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast()))) {
 244     if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
 245       fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
 246       iTrkl0++;
 247     }
 248     else {
 249       if (iTrkl0 >= fTracklets[2*layer + 0]->GetEntriesFast()) {
 250         fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
 251         iTrkl1++;
 252       }
 253       else {
 254         if (((AliTRDtrackletGTU*) fTracklets[2*layer + 1]->At(iTrkl1))->GetZbin() <
 255             ((AliTRDtrackletGTU*) fTracklets[2*layer + 0]->At(iTrkl0))->GetZbin()) {
 256           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
 257           iTrkl1++;
 258
 259         }
 260         else {
 261           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
 262           iTrkl0++;
 263         }
 264       }
 265     }
 266     ++nTracklets;
 267   }
 268
 269   TIter next(fTrackletsPostInput[layer]);
 270
 271   while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
 272     trk->SetIndex(fTrackletsPostInput[layer]->IndexOf(trk));
 273
 274     Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
 275     alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
 276     // wrapping projected alpha as in hardware
 277     if ((alpha < -64) || (alpha > 63))
 278       AliDebug(1, Form("alpha out of range: %i", alpha));
 279     alpha = alpha & 0x7f;
 280     if (alpha & 0x40)
 281       trk->SetAlpha(0xffffffc0 | alpha);
 282     else
 283       trk->SetAlpha(alpha);
 284
 285     Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer));
 286     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
 287     trk->SetYProj(yproj);
 288
 289     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
 290
 291     AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i, c: %5i, yt: %5i",
 292                       trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
 293                       trk->GetYProj(), trk->GetAlpha(), trk->GetPID(),
 294                       fGtuParam->GetCiYProj(layer), fGtuParam->GetYt(fStack, layer, trk->GetZbin()) ));
 295   }
 296   return kTRUE;
 297 }
 298
 299 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
 300 {
 301   // run the z-channel unit
 302
 303   TIter next(fTrackletsPostInput[layer]);
 304
 305   while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
 306     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 307       if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
 308         trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
 309
 310         TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
 311         AliTRDtrackletGTU *t = 0x0;
 312         while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
 313           if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
 314               ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
 315             break;
 316           }
 317         }
 318         if (t)
 319           fZChannelTracklets[layer][zch].AddAfter(t, trk);
 320         else
 321           fZChannelTracklets[layer][zch].AddFirst(trk);
 322       }
 323     }
 324     AliDebug(10, Form("stack %d, layer %2d: 0x%08x Z(2..0)=%i/%i/%i",
 325                       fStack, layer, trk->GetTrackletWord(),
 326                       fGtuParam->IsInZChannel(fStack, layer, 2, trk->GetZbin()) ? trk->GetSubChannel(2) : -1,
 327                       fGtuParam->IsInZChannel(fStack, layer, 1, trk->GetZbin()) ? trk->GetSubChannel(1) : -1,
 328                       fGtuParam->IsInZChannel(fStack, layer, 0, trk->GetZbin()) ? trk->GetSubChannel(0) : -1
 329                       ));
 330   }
 331   return kTRUE;
 332 }
 333
 334 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
 335 {
 336   // run the track finding
 337
 338    Int_t        *notr = new Int_t[fGtuParam->GetNLayers()];
 339    Int_t        *ptrA = new Int_t[fGtuParam->GetNLayers()];
 340    Int_t        *ptrB = new Int_t[fGtuParam->GetNLayers()];
 341
 342    Bool_t       *bHitA     = new Bool_t[fGtuParam->GetNLayers()];
 343    Bool_t       *bHitB     = new Bool_t[fGtuParam->GetNLayers()];
 344    Bool_t       *bAligned  = new Bool_t[fGtuParam->GetNLayers()];
 345    Bool_t       *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
 346    Bool_t       *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
 347    Bool_t       *bDone     = new Bool_t[fGtuParam->GetNLayers()];
 348    Bool_t        ready;
 349
 350    Int_t        *inc      = new Int_t[fGtuParam->GetNLayers()];
 351    Int_t        *incprime = new Int_t[fGtuParam->GetNLayers()];
 352
 353 // ----- signals within current layer -----
 354    Int_t        yPlus;
 355    Int_t        yMinus;
 356    Int_t        ybPlus;
 357    Int_t        ybMinus;
 358    Int_t        alphaPlus;
 359    Int_t        alphaMinus;
 360    Int_t        nHits;
 361    Int_t        nUnc;
 362    Int_t        nWayBeyond;
 363
 364    AliTRDtrackletGTU    *trkRA  = 0x0;  // reference tracklet A
 365    AliTRDtrackletGTU    *trkRB  = 0x0;  // reference tracklet B
 366    AliTRDtrackletGTU    *trkA   = 0x0;  // tracklet A in current layer
 367    AliTRDtrackletGTU    *trkB   = 0x0;  // tracklet B in current layer
 368 /*
 369    AliTRDtrackletGTU    *trkEnd = new AliTRDtrackletGTU();
 370    for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
 371        trkEnd->SetSubChannel(i, 7);
 372    trkEnd->SetYProj(0);
 373    trkEnd->SetAlpha(0);
 374 */
 375
 376    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
 377      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
 378      AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
 379
 380      ready = kFALSE; // ready if all channels done
 381
 382 //      for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
 383 //        for (Int_t iTrkl = 0; iTrkl < fZChannelTracklets[iLayer][zch].GetEntries(); iTrkl++) {
 384 //       printf("%4i/%4i(%i,%i) ",
 385 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetYProj(),
 386 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetAlpha(),
 387 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetIndex(),
 388 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetSubChannel(zch)
 389 //              );
 390 //        }
 391 //        printf("\n");
 392 //      }
 393
 394      for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 395        notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
 396        ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
 397        ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
 398
 399 // not necessary here
 400 //       bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
 401 //       bDone[layer] = (notr[layer] <= 0);
 402 //       ready = ready && bDone[layer];
 403
 404        if (reflayer == 1)
 405          AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
 406      }
 407
 408      if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
 409        continue;
 410
 411      while (!ready) {
 412        // ----- reference tracklets -----
 413        trkRA = 0x0;
 414        trkRB = 0x0;
 415        if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
 416          trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
 417        else  {
 418          AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
 419          break;
 420        }
 421
 422        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
 423          trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
 424
 425        yPlus      = trkRA->GetYProj() + fGtuParam->GetDeltaY();
 426        yMinus     = trkRA->GetYProj() - fGtuParam->GetDeltaY();
 427        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
 428        alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
 429        if (trkRB) {
 430          ybPlus           = trkRB->GetYProj() + fGtuParam->GetDeltaY();
 431          ybMinus          = trkRB->GetYProj() - fGtuParam->GetDeltaY();
 432        }
 433        else { // irrelevant (should be, is it?)
 434          ybPlus           = trkRA->GetYProj() + fGtuParam->GetDeltaY();
 435          ybMinus          = trkRA->GetYProj() - fGtuParam->GetDeltaY();
 436        }
 437
 438        nHits      = 0;
 439        nUnc       = 0;
 440        nWayBeyond = 0;
 441
 442        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 443          bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
 444          inc[layer] = incprime[layer] = 0;
 445
 446          if (layer == reflayer) {
 447            bHitA[layer] = kTRUE;
 448            bAligned[layer] = kTRUE;
 449            nHits++;
 450            continue;
 451          }
 452
 453          trkA = 0x0;
 454          trkB = 0x0;
 455          if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
 456            trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
 457          if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
 458            trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
 459
 460          bAlignedA[layer] = kFALSE;
 461          bAlignedB[layer] = kFALSE;
 462
 463          if (trkA) {
 464            bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
 465                             !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus)  ) &&
 466                             !(trkA->GetAlpha() < alphaMinus) &&
 467                             !(trkA->GetAlpha() > alphaPlus) );
 468            bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
 469          }
 470          else {
 471            bHitA[layer] = 0;
 472            bAlignedA[layer] = kTRUE;
 473          }
 474
 475          if (trkB) {
 476            bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
 477                             !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
 478                             !(trkB->GetAlpha() < alphaMinus) &&
 479                             !(trkB->GetAlpha() > alphaPlus) );
 480            bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
 481          }
 482          else {
 483            bHitB[layer] = 0;
 484            bAlignedB[layer] = kTRUE;
 485          }
 486
 487          bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
 488 //       bAligned[layer] = bAlignedA[layer]; //???
 489
 490          if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
 491            nHits++;
 492          else if (!bAligned[layer] )
 493            nUnc++;
 494          if (trkRB) {
 495            if (trkA) {
 496              if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
 497                nWayBeyond++;
 498            }
 499            else
 500              nWayBeyond++;
 501          }
 502
 503          //  pre-calculation for the layer shifting (alignment w. r. t. trkRB)
 504          if (trkA) {
 505              if(trkRB) {
 506                  if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
 507                      incprime[layer] = 1;
 508                  else
 509                      incprime[layer] = 0;
 510              }
 511              else
 512                  incprime[layer] = 1;
 513
 514              if (trkB) {
 515                  if (trkRB) {
 516                      if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
 517                          incprime[layer] = 2;
 518                  }
 519                  else
 520                      incprime[layer] = 2;
 521              }
 522          }
 523        } // end of loop over layers
 524
 525        AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
 526                        nHits, (nHits >= 4) ? "-> track found" : ""));
 527
 528        if (nHits >= 4) {
 529          // ----- track registration -----
 530          AliTRDtrackGTU *track = new AliTRDtrackGTU();
 531          track->SetSector(fSector);
 532          track->SetStack(fStack);
 533          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
 534            if (bHitA[layer] || layer == reflayer)
 535              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
 536            else if (bHitB[layer])
 537              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
 538          }
 539
 540          Bool_t registerTrack = kTRUE;
 541          for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
 542            if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
 543              if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
 544                AliDebug(1,"Not registered");
 545                  registerTrack = kFALSE;
 546              }
 547            }
 548          }
 549          if (registerTrack) {
 550              track->SetZChannel(zch);
 551              track->SetRefLayerIdx(refLayerIdx);
 552              fTracks[zch][refLayerIdx].Add(track);
 553          }
 554          else
 555            delete track;
 556        }
 557
 558        if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
 559           inc[reflayer] = 0;
 560        else if (nWayBeyond > 2) // no track possible for both reference tracklets
 561          inc[reflayer] = 2;
 562        else
 563          inc[reflayer] = 1;
 564
 565        if (inc[reflayer] != 0) // reflayer is shifted
 566          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 567            if (layer == reflayer)
 568              continue;
 569            inc[layer] = incprime[layer];
 570          }
 571        else { // reflayer is not shifted
 572          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 573            if (layer == reflayer || notr[layer] == 0)
 574              continue;
 575            inc[layer] = 0;
 576            trkA = 0x0;
 577            trkB = 0x0;
 578            if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
 579              trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
 580
 581            if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
 582              trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
 583
 584            if (trkA) {
 585              if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
 586                   !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) )  // trkA could hit trkRA
 587              // if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) )
 588                    inc[layer] = 0;
 589                else if (trkB) {
 590                    if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
 591                        inc[layer] = 2;
 592                    else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
 593                        inc[layer] = 1;
 594                    else
 595                        inc[layer] = incprime[layer];
 596                }
 597                else
 598                    inc[layer] = incprime[layer];
 599            }
 600          }
 601        }
 602
 603        ready = kTRUE;
 604        for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
 605
 606          bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
 607          ready = ready && bDone[layer];
 608
 609          if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
 610            AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
 611
 612 //       AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
 613          AliDebug(10,Form("    Layer: %i   %2i(%2i, %2i, %4i, %3i)%s%s   %2i(%2i, %2i, %4i, %3i)%s%s   +%i   %s  (no: %i)",
 614                           layer,
 615                           ptrA[layer],
 616                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
 617                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
 618                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
 619                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
 620                           bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
 621                           ptrB[layer],
 622                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
 623                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
 624                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
 625                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
 626                           bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
 627                           inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
 628          ptrA[layer] += inc[layer];
 629          ptrB[layer] += inc[layer];
 630        }
 631
 632      } // end of while
 633    } // end of loop over reflayer
 634
 635    delete [] bHitA;
 636    delete [] bHitB;
 637    delete [] bAligned;
 638    delete [] bDone;
 639    delete [] inc;
 640    delete [] incprime;
 641    delete [] bAlignedA;
 642    delete [] bAlignedB;
 643    delete [] notr;
 644    delete [] ptrA;
 645    delete [] ptrB;
 646
 647    return kTRUE;
 648 }
 649
 650 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
 651 {
 652     TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
 653     TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
 654
 655     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 656         tracksRefMerged[zch] = new TList;
 657         tracksRefUnique[zch] = new TList;
 658     }
 659
 660     TList *tracksZMergedStage0 = new TList;
 661     TList *tracksZUniqueStage0 = new TList;
 662
 663     TList **tracksZSplitted = new TList*[2];
 664     for (Int_t i = 0; i < 2; i++)
 665         tracksZSplitted[i] = new TList;
 666
 667     TList *tracksZMergedStage1 = new TList;
 668
 669     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
 670
 671     // Bool_t done = kFALSE;
 672     Int_t minIdx = 0;
 673     AliTRDtrackGTU *trkStage0 = 0x0;
 674
 675     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 676         // ----- Merging and Unification in Reflayers (seed_merger) -----
 677         do {
 678           // done = kTRUE;
 679             trkStage0 = 0x0;
 680             for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
 681                 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
 682                 if (trkInRefLayer[refLayerIdx] == 0) {
 683                     continue;
 684                 }
 685                 else if (trkStage0 == 0x0 ) {
 686                     trkStage0 = trkInRefLayer[refLayerIdx];
 687                     minIdx = refLayerIdx;
 688                     // done = kFALSE;
 689                 }
 690                 else if ( (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel()) ||
 691                           ((trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel()) &&
 692                            ((trkInRefLayer[refLayerIdx]->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
 693                     minIdx = refLayerIdx;
 694                     trkStage0 = trkInRefLayer[refLayerIdx];
 695                     // done = kFALSE;
 696                 }
 697             }
 698             if (!trkStage0)
 699               break;
 700             tracksRefMerged[zch]->Add(trkStage0);
 701             fTracks[zch][minIdx].RemoveFirst();
 702         } while (trkStage0 != 0);
 703
 704         Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
 705
 706         AliDebug(2, Form("zchannel %i", zch));
 707         TIter trackRefMerged(tracksRefMerged[zch]);
 708         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
 709           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
 710                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 711                            trk->GetTrackletIndex(5),
 712                            trk->GetTrackletIndex(4),
 713                            trk->GetTrackletIndex(3),
 714                            trk->GetTrackletIndex(2),
 715                            trk->GetTrackletIndex(1),
 716                            trk->GetTrackletIndex(0),
 717                            trk->GetYapprox() >> 3,
 718                            trk->GetZSubChannel()));
 719         AliDebug(2, "uniquified:");
 720         TIter trackRefUnique(tracksRefUnique[zch]);
 721         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
 722           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
 723                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 724                            trk->GetTrackletIndex(5),
 725                            trk->GetTrackletIndex(4),
 726                            trk->GetTrackletIndex(3),
 727                            trk->GetTrackletIndex(2),
 728                            trk->GetTrackletIndex(1),
 729                            trk->GetTrackletIndex(0),
 730                            trk->GetYapprox() >> 3,
 731                            trk->GetZSubChannel()));
 732     }
 733
 734 // ----- Merging in zchannels - 1st stage -----
 735
 736     if (AliTRDgtuParam::GetUseGTUmerge()) {
 737       Int_t notEmpty;
 738       do {
 739         Bool_t lowerThan[3] = { kFALSE, kFALSE, kFALSE };
 740         AliTRDtrackGTU *trk[3] = { 0x0, 0x0, 0x0 };
 741         for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel)
 742           trk[iChannel] = (AliTRDtrackGTU*) tracksRefUnique[iChannel]->First();
 743
 744         for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel) {
 745           AliTRDtrackGTU *trk1 = trk[iChannel];
 746           AliTRDtrackGTU *trk2 = trk[(iChannel + 1) % 3];
 747           if (trk1 && trk2) {
 748             Int_t sortnum1 = (trk1->GetZChannel() + 3 * trk1->GetZSubChannel()) / 2 - 1;
 749             Int_t sortnum2 = (trk2->GetZChannel() + 3 * trk2->GetZSubChannel()) / 2 - 1;
 750             AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
 751                              trk1->GetZChannel(), trk2->GetZChannel(),
 752                              sortnum1, sortnum2));
 753             if ( (sortnum1 < sortnum2) ||
 754                  ((sortnum1 == sortnum2) &&
 755                   ((trk1->GetYapprox() >> 3) < (trk2->GetYapprox() >> 3)) ) ) {
 756               lowerThan[iChannel] = kTRUE;
 757             }
 758           }
 759         }
 760
 761         notEmpty = (trk[2] ? (1 << 2) : 0) |
 762           (trk[1] ? (1 << 1) : 0) |
 763           (trk[0] ? (1 << 0) : 0);
 764         Int_t pop = -1;
 765
 766         switch (notEmpty) {
 767           // one track only
 768         case 0x1:
 769           pop = 0;
 770           break;
 771         case 0x2:
 772           pop = 1;
 773           break;
 774         case 0x4:
 775           pop = 2;
 776           break;
 777
 778           // two tracks
 779         case 0x3:
 780           if (lowerThan[0])
 781             pop = 0;
 782           else
 783             pop = 1;
 784           break;
 785         case 0x5:
 786           if (lowerThan[2])
 787             pop = 2;
 788           else
 789             pop = 0;
 790           break;
 791         case 0x6:
 792           if (lowerThan[1])
 793             pop = 1;
 794           else
 795             pop = 2;
 796           break;
 797
 798           // three tracks
 799         case 0x7:
 800           if (lowerThan[0]) {
 801             if (lowerThan[2])
 802               pop = 2;
 803             else
 804               pop = 0;
 805           } else {
 806             if (lowerThan[1])
 807               pop = 1;
 808             else
 809               pop = 2;
 810           }
 811           break;
 812
 813           // no tracks
 814         default:
 815           // nop
 816           ;
 817         }
 818
 819         if (pop > -1) {
 820           tracksZMergedStage0->Add(trk[pop]);
 821           tracksRefUnique[pop]->RemoveFirst();
 822         }
 823       } while (notEmpty);
 824     }
 825     else {
 826       // there is still a difference between this implementation and
 827       // the hardware algorithm, only for expert use
 828
 829       do {
 830         // done = kTRUE;
 831         trkStage0 = 0x0;
 832         // compare tracks from all adjacent zchannels
 833         // (including comparison of channels 2 and 0)
 834         for (Int_t i = fGtuParam->GetNZChannels() - 1; i > -1; i--) {
 835           Int_t zch = i % 3;
 836           AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
 837           if (trk == 0) {
 838             continue;
 839           }
 840           else if (trkStage0 == 0x0 ) {
 841             trkStage0 = trk;
 842             minIdx = zch;
 843             // done = kFALSE;
 844           }
 845           else {
 846             Int_t sortnum1 = (trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1;
 847             Int_t sortnum2 = (trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 - 1;
 848             AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
 849                              trk->GetZChannel(), trkStage0->GetZChannel(),
 850                              sortnum1, sortnum2));
 851             if ( (sortnum1 < sortnum2) ||
 852                  ((sortnum1 == sortnum2) &&
 853                   ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
 854               minIdx = zch;
 855               trkStage0 = trk;
 856               // done = kFALSE;
 857             }
 858           }
 859         }
 860
 861         if (!trkStage0)
 862           break;
 863         tracksZMergedStage0->Add(trkStage0);
 864         tracksRefUnique[minIdx]->RemoveFirst();
 865       } while (trkStage0 != 0);
 866     }
 867
 868     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
 869
 870     AliDebug(2, "stage 0:");
 871     TIter trackZMergedStage0(tracksZMergedStage0);
 872     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
 873       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 874                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 875                        trk->GetTrackletIndex(5),
 876                        trk->GetTrackletIndex(4),
 877                        trk->GetTrackletIndex(3),
 878                        trk->GetTrackletIndex(2),
 879                        trk->GetTrackletIndex(1),
 880                        trk->GetTrackletIndex(0),
 881                        trk->GetYapprox() >> 3,
 882                        trk->GetZChannel(),
 883                        trk->GetZSubChannel()));
 884
 885     AliDebug(2, "uniquified:");
 886     TIter trackZUniqueStage0(tracksZUniqueStage0);
 887     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
 888       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 889                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 890                        trk->GetTrackletIndex(5),
 891                        trk->GetTrackletIndex(4),
 892                        trk->GetTrackletIndex(3),
 893                        trk->GetTrackletIndex(2),
 894                        trk->GetTrackletIndex(1),
 895                        trk->GetTrackletIndex(0),
 896                        trk->GetYapprox() >> 3,
 897                        trk->GetZChannel(),
 898                        trk->GetZSubChannel()));
 899
 900 // ----- Splitting in z -----
 901
 902     TIter next(tracksZUniqueStage0);
 903     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
 904         tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
 905     }
 906
 907     for (Int_t i = 0; i < 2; i++) {
 908       AliDebug(2, Form("split %i:", i));
 909       TIter trackZSplit(tracksZSplitted[i]);
 910       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
 911         AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 912                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 913                          trk->GetTrackletIndex(5),
 914                          trk->GetTrackletIndex(4),
 915                          trk->GetTrackletIndex(3),
 916                          trk->GetTrackletIndex(2),
 917                          trk->GetTrackletIndex(1),
 918                          trk->GetTrackletIndex(0),
 919                          trk->GetYapprox() >> 3,
 920                          trk->GetZChannel(),
 921                          trk->GetZSubChannel()));
 922     }
 923
 924 // ----- Merging in zchanels - 2nd stage -----
 925
 926     do {
 927       // done = kTRUE;
 928         trkStage0 = 0x0;
 929         for (Int_t i = 1; i >= 0; i--) {
 930             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
 931             if (trk == 0) {
 932                 continue;
 933             }
 934             else if (trkStage0 == 0x0 ) {
 935                 trkStage0 = trk;
 936                 minIdx = i;
 937                 // done = kFALSE;
 938             }
 939             else if (  (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) <  ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
 940                       ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) &&
 941                        ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
 942                 minIdx = i;
 943                 trkStage0 = trk;
 944                 // done = kFALSE;
 945             }
 946         }
 947
 948         if (!trkStage0)
 949           break;
 950         tracksZMergedStage1->Add(trkStage0);
 951         tracksZSplitted[minIdx]->RemoveFirst();
 952     } while (trkStage0 != 0);
 953
 954     Uniquifier(tracksZMergedStage1, ListOfTracks);
 955
 956     AliDebug(2, "merged:");
 957     TIter trackZMergedStage1(tracksZMergedStage1);
 958     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
 959       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 960                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 961                        trk->GetTrackletIndex(5),
 962                        trk->GetTrackletIndex(4),
 963                        trk->GetTrackletIndex(3),
 964                        trk->GetTrackletIndex(2),
 965                        trk->GetTrackletIndex(1),
 966                        trk->GetTrackletIndex(0),
 967                        trk->GetYapprox() >> 3,
 968                        trk->GetZChannel(),
 969                        trk->GetZSubChannel()));
 970
 971     AliDebug(2, "uniquified:");
 972     TIter track(ListOfTracks);
 973     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
 974       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 975                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 976                        trk->GetTrackletIndex(5),
 977                        trk->GetTrackletIndex(4),
 978                        trk->GetTrackletIndex(3),
 979                        trk->GetTrackletIndex(2),
 980                        trk->GetTrackletIndex(1),
 981                        trk->GetTrackletIndex(0),
 982                        trk->GetYapprox() >> 3,
 983                        trk->GetZChannel(),
 984                        trk->GetZSubChannel()));
 985
 986     // cleaning up
 987     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 988       delete tracksRefMerged[zch];
 989       delete tracksRefUnique[zch];
 990     }
 991     delete [] tracksRefMerged;
 992     delete [] tracksRefUnique;
 993
 994
 995     delete tracksZMergedStage0;
 996     delete tracksZUniqueStage0;
 997
 998     for (Int_t i = 0; i < 2; i++)
 999       delete tracksZSplitted[i];
1000     delete [] tracksZSplitted;
1001
1002     delete tracksZMergedStage1;
1003
1004     delete [] trkInRefLayer;
1005
1006     return kTRUE;
1007 }
1008
1009 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
1010 {
1011   // run the track reconstruction for all tracks in the list
1012
1013   TIter next(ListOfTracks);
1014   while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
1015     CalculateTrackParams(track);
1016     CalculatePID(track);
1017     AliDebug(1, Form("track with pid = %i", track->GetPID()));
1018   }
1019   return kTRUE;
1020 }
1021
1022 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
1023 {
1024   // calculate PID for the given track
1025   if (!track) {
1026     AliError("No track to calculate!");
1027     return kFALSE;
1028   }
1029
1030   if (AliTRDgtuParam::GetUseGTUconst()) {
1031     // averaging as in GTU
1032     AliDebug(1, "using GTU constants for PID calculation");
1033     ULong64_t coeff;
1034
1035     // selection of coefficient for averaging
1036     switch(track->GetTrackletMask()){
1037     case 0x3f:
1038       // 6 tracklets
1039       coeff=0x5558; // approx. 1/6
1040       break;
1041
1042     case 0x3e:
1043     case 0x3d:
1044     case 0x3b:
1045     case 0x37:
1046     case 0x2f:
1047     case 0x1f:
1048       // 5 tracklets
1049       coeff=0x6666; // approx. 1/5
1050       break;
1051
1052     default:
1053       // 4 tracklets
1054       coeff=0x8000; // = 0.25
1055     }
1056     coeff &= 0x1ffff; // 17-bit constant
1057
1058     ULong64_t sum = 0;
1059     Int_t i = 0;
1060     for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
1061       if ((track->GetTrackletMask() >> iLayer) & 1) {
1062         sum += track->GetTracklet(iLayer)->GetPID();
1063         ++i;
1064       }
1065     }
1066
1067     Float_t av = 1./i * sum;
1068     sum = sum & 0x7ff;
1069     ULong64_t prod   = (sum * coeff) & 0xfffffffffull;
1070     ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
1071
1072     if (TMath::Abs((prodFinal & 0xff) - av) > 0.5)
1073       AliError(Form("problem with PID averaging (hw <-> ar): %3lli <-> %4.1f", prodFinal & 0xff, av));
1074     track->SetPID(prodFinal & 0xff);
1075
1076     return kTRUE;
1077   }
1078   else {
1079
1080     // simple averaging
1081     Int_t nTracklets = 0;
1082     Int_t pidSum = 0;
1083     for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1084       if (!track->IsTrackletInLayer(layer)) {
1085         continue;
1086       }
1087       AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1088       pidSum += trk->GetPID();
1089       nTracklets++;
1090     }
1091
1092     if (nTracklets>0)
1093       track->SetPID(pidSum/nTracklets);
1094     else
1095       AliError("Track without contributing tracklets, no PID assigned");
1096
1097     return kTRUE;
1098   }
1099 }
1100
1101 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
1102 {
1103   // calculate the track parameters
1104
1105   if (!track) {
1106     AliError("No track to calculate!");
1107     return kFALSE;
1108   }
1109
1110   Int_t a = 0;
1111   Float_t b = 0;
1112   Float_t c = 0;
1113   Float_t x1;
1114   Float_t x2;
1115
1116   AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
1117
1118   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1119     if (!track->IsTrackletInLayer(layer)) {
1120       continue;
1121     }
1122     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1123     if (!trk) {
1124       AliError(Form("Could not get tracklet in layer %i\n", layer));
1125       continue;
1126     }
1127     AliDebug(10,Form("  layer %i trk yprime: %6i, aki: %6i", layer, trk->GetYPrime(),
1128                      fGtuParam->GetAki(track->GetTrackletMask(), layer)));
1129     a += (((fGtuParam->GetAki(track->GetTrackletMask(), layer) * trk->GetYPrime()) >> 7) + 1) >> 1;
1130     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1131     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1132   }
1133   if (a < 0)
1134       a += 3;
1135   a = a >> 2;
1136
1137   track->SetFitParams(a << 2, b, c);
1138
1139   fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1140
1141   AliDebug(5,Form("  Track parameters: a16 = %i, a18 = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f (trkl mask: %i)",
1142                   a, a << 2, b, c, x1, x2, track->GetPt(), track->GetTrackletMask()));
1143
1144   return kTRUE;
1145 }
1146
1147
1148 Bool_t AliTRDgtuTMU::Uniquifier(const TList *inlist, TList *outlist)
1149 {
1150   // remove multiple occurences of the same track
1151
1152     TIter next(inlist);
1153     AliTRDtrackGTU *trkStage0 = 0x0;
1154     AliTRDtrackGTU *trkStage1 = 0x0;
1155
1156     do {
1157       if (trkStage0 != trkStage1)
1158         delete trkStage0;
1159
1160         trkStage0 = (AliTRDtrackGTU*) next();
1161
1162         Bool_t tracksEqual = kFALSE;
1163         if (trkStage0 != 0 && trkStage1 != 0) {
1164             for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1165                 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1166                     tracksEqual = kTRUE;
1167                     break;
1168                 }
1169             }
1170         }
1171
1172         if (tracksEqual) {
1173           if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets()) {
1174             delete trkStage1;
1175             trkStage1 = trkStage0;
1176           }
1177         }
1178         else {
1179             if (trkStage1 != 0x0)
1180                 outlist->Add(trkStage1);
1181             trkStage1 = trkStage0;
1182         }
1183
1184     } while (trkStage1 != 0x0);
1185     return kTRUE;
1186 }