TRD/AliTRDgtuTMU.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliTRDgtuTMU.cxx 28397 2008-09-02 09:33:00Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Track Matching Unit (TMU) simulation                                  //
  21 //                                                                        //
  22 //  Author: J. Klein (Jochen.Klein@cern.ch)                               //
  23 //                                                                        //
  24 ////////////////////////////////////////////////////////////////////////////
  25
  26 #include "TTree.h"
  27 #include "TList.h"
  28 #include "TVectorD.h"
  29 #include "TMath.h"
  30
  31 #include "AliESDEvent.h"
  32 #include "AliESDTrdTrack.h"
  33
  34 #include "AliLog.h"
  35 #include "AliTRDgeometry.h"
  36 #include "AliTRDpadPlane.h"
  37
  38 #include "AliTRDgtuParam.h"
  39 #include "AliTRDgtuTMU.h"
  40 #include "AliTRDtrackGTU.h"
  41
  42 ClassImp(AliTRDgtuTMU)
  43
  44 AliTRDgtuTMU::AliTRDgtuTMU(Int_t stack, Int_t sector) :
  45   TObject(),
  46   fTracklets(0x0),
  47   fTrackletsPostInput(0x0),
  48   fZChannelTracklets(0x0),
  49   fTracks(0x0),
  50   fGtuParam(0x0),
  51   fStack(-1),
  52   fSector(-1)
  53 {
  54   // constructor which initializes the position information of the TMU
  55
  56   fGtuParam = AliTRDgtuParam::Instance();
  57
  58   // store tracklets per link
  59   fTracklets = new TObjArray*[fGtuParam->GetNLinks()];
  60   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
  61     fTracklets[iLink] = new TObjArray();
  62   }
  63
  64   // tracklets after input units per layer
  65   fTrackletsPostInput = new TObjArray*[fGtuParam->GetNLayers()];
  66   fZChannelTracklets = new TList*[fGtuParam->GetNLayers()];
  67   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
  68     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
  69     fTrackletsPostInput[layer] = new TObjArray();
  70   }
  71
  72   fTracks = new TList*[fGtuParam->GetNZChannels()];
  73   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
  74       fTracks[zch] = new TList[fGtuParam->GetNRefLayers()];
  75   }
  76
  77   if (stack > -1)
  78       SetStack(stack);
  79   if (sector > -1)
  80       SetSector(sector);
  81 }
  82
  83 AliTRDgtuTMU::~AliTRDgtuTMU()
  84 {
  85   // destructor
  86
  87   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
  88     delete [] fTracks[zch];
  89   }
  90   delete [] fTracks;
  91   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
  92     delete [] fZChannelTracklets[layer];
  93     //    delete [] fTrackletsPostInput[layer];
  94   }
  95   delete [] fZChannelTracklets;
  96   //  delete [] fTrackletsPostInput;
  97
  98   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
  99     delete fTracklets[iLink];
 100   }
 101   delete [] fTracklets;
 102 }
 103
 104 Bool_t AliTRDgtuTMU::SetSector(Int_t sector)
 105 {
 106   // set the sector
 107
 108   if (sector > -1 && sector < fGtuParam->GetGeo()->Nsector() ) {
 109     fSector = sector;
 110     return kTRUE;
 111   }
 112
 113   AliError(Form("Invalid sector given: %i", sector));
 114   return kFALSE;
 115 }
 116
 117 Bool_t AliTRDgtuTMU::SetStack(Int_t stack)
 118 {
 119   // Set the stack (necessary for tracking)
 120
 121   if (stack > -1 && stack < fGtuParam->GetGeo()->Nstack() ) {
 122     fStack = stack;
 123     return kTRUE;
 124   }
 125
 126   AliError(Form("Invalid stack given: %i", stack));
 127   return kFALSE;
 128 }
 129
 130 Bool_t AliTRDgtuTMU::Reset()
 131 {
 132   // delete all tracks
 133
 134   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 135     for (Int_t reflayeridx = 0; reflayeridx < fGtuParam->GetNRefLayers(); reflayeridx++) {
 136       fTracks[zch][reflayeridx].Clear();
 137     }
 138   }
 139
 140   // delete all added tracklets
 141   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
 142     fTracklets[iLink]->Clear();
 143   }
 144   for (Int_t layer = 0; layer < fGtuParam->GetNLinks()/2; layer++) {
 145     fTrackletsPostInput[layer]->Clear();
 146     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++)
 147       fZChannelTracklets[layer][zch].Clear();
 148   }
 149
 150   fSector = -1;
 151   fStack = -1;
 152
 153   return kTRUE;
 154 }
 155
 156 Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletGTU *tracklet, Int_t link)
 157 {
 158   // add a tracklet on the given link
 159
 160   fTracklets[link]->Add(tracklet);
 161   return kTRUE;
 162 }
 163
 164
 165 Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
 166 {
 167   // performs the analysis as in a TMU module of the GTU, i. e.
 168   // track matching
 169   // calculation of track parameteres (pt, deflection, ???)
 170
 171   if (fStack < 0 || fSector < 0) {
 172     AliError("No valid stack/sector set for this TMU! No tracking!");
 173     return kFALSE;
 174   }
 175
 176   // ----- Input units -----
 177   AliDebug(1,"--------- Running Input units ----------");
 178   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 179     if (!RunInputUnit(layer)) {
 180       AliError(Form("Input unit in layer %i failed", layer));
 181       return kFALSE;
 182     }
 183   }
 184
 185   // ----- Z-channel units -----
 186   AliDebug(1,"--------- Running Z-channel units ----------");
 187   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
 188     fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
 189     if (!RunZChannelUnit(layer)) {
 190       AliError(Form("Z-Channel unit in layer %i failed", layer));
 191       return kFALSE;
 192     }
 193   }
 194
 195   // ----- track finding -----
 196   AliDebug(1,"--------- Running tracking units ----------");
 197   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 198     if (!RunTrackFinder(zch, ListOfTracks)) {
 199       AliError(Form("Track Finder in z-channel %i failed", zch));
 200       return kFALSE;
 201     }
 202   }
 203
 204   // ----- Track Merging -----
 205   if (!RunTrackMerging(ListOfTracks)) {
 206     AliError("Track merging failed");
 207     return kFALSE;
 208   }
 209
 210   // ----- track reconstruction -----
 211   if (!RunTrackReconstruction(ListOfTracks)) {
 212     AliError("Track reconstruction failed");
 213     return kFALSE;
 214   }
 215
 216   // ----- label calculation and ESD storage -----
 217   TIter next(ListOfTracks);
 218   while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
 219     trk->CookLabel();
 220     if (esd) {
 221       AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
 222       esd->AddTrdTrack(trdtrack);
 223       delete trdtrack;
 224     }
 225   }
 226
 227   return kTRUE;
 228 }
 229
 230 Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
 231 {
 232   // precalculations for the tracking and reconstruction
 233
 234   Int_t iTrkl0 = 0; // A-side tracklet
 235   Int_t iTrkl1 = 0; // B-side tracklet
 236
 237   while ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) || (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast())) {
 238     if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
 239       fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
 240       iTrkl0++;
 241     }
 242     else {
 243       if (iTrkl0 >= fTracklets[2*layer + 0]->GetEntriesFast()) {
 244         fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
 245         iTrkl1++;
 246       }
 247       else {
 248         if (((AliTRDtrackletGTU*) fTracklets[2*layer + 1]->At(iTrkl1))->GetZbin() <
 249             ((AliTRDtrackletGTU*) fTracklets[2*layer + 0]->At(iTrkl0))->GetZbin()) {
 250           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
 251           iTrkl1++;
 252
 253         }
 254         else {
 255           fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
 256           iTrkl0++;
 257         }
 258       }
 259     }
 260   }
 261
 262   TIter next(fTrackletsPostInput[layer]);
 263
 264   while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
 265     trk->SetIndex(fTrackletsPostInput[layer]->IndexOf(trk));
 266
 267     Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
 268     alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
 269     // wrapping projected alpha as in hardware
 270     if ((alpha < -64) || (alpha > 63))
 271       AliError(Form("alpha out of range: %i", alpha));
 272     alpha = alpha & 0x7f;
 273     if (alpha & 0x40)
 274       trk->SetAlpha(0xffffffc0 | alpha);
 275     else
 276       trk->SetAlpha(alpha);
 277
 278     Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer)); //??? sign?
 279     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
 280     trk->SetYProj(yproj);
 281
 282     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
 283
 284     AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i",
 285                       trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
 286                       trk->GetYProj(), trk->GetAlpha(), trk->GetPID() ));
 287   }
 288   return kTRUE;
 289 }
 290
 291 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
 292 {
 293   // run the z-channel unit
 294
 295   TIter next(fTrackletsPostInput[layer]);
 296
 297   while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
 298     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 299       if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
 300         trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
 301
 302         TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
 303         AliTRDtrackletGTU *t = 0x0;
 304         while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
 305           if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
 306               ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
 307             break;
 308           }
 309         }
 310         if (t)
 311           fZChannelTracklets[layer][zch].AddAfter(t, trk);
 312         else
 313           fZChannelTracklets[layer][zch].AddFirst(trk);
 314       }
 315     }
 316     AliDebug(10, Form("stack %d, layer %2d: 0x%08x Z(2..0)=%i/%i/%i",
 317                       fStack, layer, trk->GetTrackletWord(),
 318                       fGtuParam->IsInZChannel(fStack, layer, 2, trk->GetZbin()) ? trk->GetSubChannel(2) : -1,
 319                       fGtuParam->IsInZChannel(fStack, layer, 1, trk->GetZbin()) ? trk->GetSubChannel(1) : -1,
 320                       fGtuParam->IsInZChannel(fStack, layer, 0, trk->GetZbin()) ? trk->GetSubChannel(0) : -1
 321                       ));
 322   }
 323   return kTRUE;
 324 }
 325
 326 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
 327 {
 328   // run the track finding
 329
 330    Int_t        *notr = new Int_t[fGtuParam->GetNLayers()];
 331    Int_t        *ptrA = new Int_t[fGtuParam->GetNLayers()];
 332    Int_t        *ptrB = new Int_t[fGtuParam->GetNLayers()];
 333
 334    Bool_t       *bHitA     = new Bool_t[fGtuParam->GetNLayers()];
 335    Bool_t       *bHitB     = new Bool_t[fGtuParam->GetNLayers()];
 336    Bool_t       *bAligned  = new Bool_t[fGtuParam->GetNLayers()];
 337    Bool_t       *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
 338    Bool_t       *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
 339    Bool_t       *bDone     = new Bool_t[fGtuParam->GetNLayers()];
 340    Bool_t        ready;
 341
 342    Int_t        *inc      = new Int_t[fGtuParam->GetNLayers()];
 343    Int_t        *incprime = new Int_t[fGtuParam->GetNLayers()];
 344
 345 // ----- signals within current layer -----
 346    Int_t        yPlus;
 347    Int_t        yMinus;
 348    Int_t        ybPlus;
 349    Int_t        ybMinus;
 350    Int_t        alphaPlus;
 351    Int_t        alphaMinus;
 352    Int_t        nHits;
 353    Int_t        nUnc;
 354    Int_t        nWayBeyond;
 355
 356    AliTRDtrackletGTU    *trkRA  = 0x0;  // reference tracklet A
 357    AliTRDtrackletGTU    *trkRB  = 0x0;  // reference tracklet B
 358    AliTRDtrackletGTU    *trkA   = 0x0;  // tracklet A in current layer
 359    AliTRDtrackletGTU    *trkB   = 0x0;  // tracklet B in current layer
 360 /*
 361    AliTRDtrackletGTU    *trkEnd = new AliTRDtrackletGTU();
 362    for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
 363        trkEnd->SetSubChannel(i, 7);
 364    trkEnd->SetYProj(0);
 365    trkEnd->SetAlpha(0);
 366 */
 367
 368    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
 369      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
 370      AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
 371
 372      ready = kFALSE; // ready if all channels done
 373
 374 //      for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
 375 //        for (Int_t iTrkl = 0; iTrkl < fZChannelTracklets[iLayer][zch].GetEntries(); iTrkl++) {
 376 //       printf("%4i/%4i(%i,%i) ",
 377 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetYProj(),
 378 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetAlpha(),
 379 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetIndex(),
 380 //              ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetSubChannel(zch)
 381 //              );
 382 //        }
 383 //        printf("\n");
 384 //      }
 385
 386      for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 387        notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
 388        ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
 389        ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
 390
 391 // not necessary here
 392 //       bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
 393 //       bDone[layer] = (notr[layer] <= 0);
 394 //       ready = ready && bDone[layer];
 395
 396        if (reflayer == 1)
 397          AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
 398      }
 399
 400      if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
 401        continue;
 402
 403      while (!ready) {
 404        // ----- reference tracklets -----
 405        trkRA = 0x0;
 406        trkRB = 0x0;
 407        if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
 408          trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
 409        else  {
 410          AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
 411          break;
 412        }
 413
 414        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
 415          trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
 416
 417        yPlus      = trkRA->GetYProj() + fGtuParam->GetDeltaY();
 418        yMinus     = trkRA->GetYProj() - fGtuParam->GetDeltaY();
 419        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
 420        alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
 421        if (trkRB) {
 422          ybPlus           = trkRB->GetYProj() + fGtuParam->GetDeltaY();
 423          ybMinus          = trkRB->GetYProj() - fGtuParam->GetDeltaY();
 424        }
 425        else { // irrelevant (should be, is it?)
 426          ybPlus           = trkRA->GetYProj() + fGtuParam->GetDeltaY();
 427          ybMinus          = trkRA->GetYProj() - fGtuParam->GetDeltaY();
 428        }
 429
 430        nHits      = 0;
 431        nUnc       = 0;
 432        nWayBeyond = 0;
 433
 434        for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 435          bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
 436          inc[layer] = incprime[layer] = 0;
 437
 438          if (layer == reflayer) {
 439            bHitA[layer] = kTRUE;
 440            bAligned[layer] = kTRUE;
 441            nHits++;
 442            continue;
 443          }
 444
 445          trkA = 0x0;
 446          trkB = 0x0;
 447          if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
 448            trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
 449          if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
 450            trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
 451
 452          bAlignedA[layer] = kFALSE;
 453          bAlignedB[layer] = kFALSE;
 454
 455          if (trkA) {
 456            bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
 457                             !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus)  ) &&
 458                             !(trkA->GetAlpha() < alphaMinus) &&
 459                             !(trkA->GetAlpha() > alphaPlus) );
 460            bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
 461          }
 462          else {
 463            bHitA[layer] = 0;
 464            bAlignedA[layer] = kTRUE;
 465          }
 466
 467          if (trkB) {
 468            bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
 469                             !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
 470                             !(trkB->GetAlpha() < alphaMinus) &&
 471                             !(trkB->GetAlpha() > alphaPlus) );
 472            bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
 473          }
 474          else {
 475            bHitB[layer] = 0;
 476            bAlignedB[layer] = kTRUE;
 477          }
 478
 479          bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
 480 //       bAligned[layer] = bAlignedA[layer]; //???
 481
 482          if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
 483            nHits++;
 484          else if (!bAligned[layer] )
 485            nUnc++;
 486          if (trkRB) {
 487            if (trkA) {
 488              if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
 489                nWayBeyond++;
 490            }
 491            else
 492              nWayBeyond++;
 493          }
 494
 495          //  pre-calculation for the layer shifting (alignment w. r. t. trkRB)
 496          if (trkA) {
 497              if(trkRB) {
 498                  if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
 499                      incprime[layer] = 1;
 500                  else
 501                      incprime[layer] = 0;
 502              }
 503              else
 504                  incprime[layer] = 1;
 505
 506              if (trkB) {
 507                  if (trkRB) {
 508                      if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
 509                          incprime[layer] = 2;
 510                  }
 511                  else
 512                      incprime[layer] = 2;
 513              }
 514          }
 515        } // end of loop over layers
 516
 517        AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
 518                        nHits, (nHits >= 4) ? "-> track found" : ""));
 519
 520        if (nHits >= 4) {
 521          // ----- track registration -----
 522          AliTRDtrackGTU *track = new AliTRDtrackGTU();
 523          track->SetSector(fSector);
 524          track->SetStack(fStack);
 525          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
 526            if (bHitA[layer] || layer == reflayer)
 527              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
 528            else if (bHitB[layer])
 529              track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
 530          }
 531
 532          Bool_t registerTrack = kTRUE;
 533          for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
 534            if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
 535              if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
 536                AliDebug(1,"Not registered");
 537                  registerTrack = kFALSE;
 538              }
 539            }
 540          }
 541          if (registerTrack) {
 542              track->SetZChannel(zch);
 543              track->SetRefLayerIdx(refLayerIdx);
 544              fTracks[zch][refLayerIdx].Add(track);
 545          }
 546        }
 547
 548        if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
 549           inc[reflayer] = 0;
 550        else if (nWayBeyond > 2) // no track possible for both reference tracklets
 551          inc[reflayer] = 2;
 552        else
 553          inc[reflayer] = 1;
 554
 555        if (inc[reflayer] != 0) // reflayer is shifted
 556          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 557            if (layer == reflayer)
 558              continue;
 559            inc[layer] = incprime[layer];
 560          }
 561        else { // reflayer is not shifted
 562          for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 563            if (layer == reflayer || notr[layer] == 0)
 564              continue;
 565            inc[layer] = 0;
 566            trkA = 0x0;
 567            trkB = 0x0;
 568            if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
 569              trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
 570
 571            if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
 572              trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
 573
 574            if (trkA) {
 575              if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
 576                   !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) )  // trkA could hit trkRA
 577              // if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) )
 578                    inc[layer] = 0;
 579                else if (trkB) {
 580                    if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
 581                        inc[layer] = 2;
 582                    else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
 583                        inc[layer] = 1;
 584                    else
 585                        inc[layer] = incprime[layer];
 586                }
 587                else
 588                    inc[layer] = incprime[layer];
 589            }
 590          }
 591        }
 592
 593        ready = kTRUE;
 594        for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
 595
 596          bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
 597          ready = ready && bDone[layer];
 598
 599          if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
 600            AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
 601
 602 //       AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
 603          AliDebug(10,Form("    Layer: %i   %2i(%2i, %2i, %4i, %3i)%s%s   %2i(%2i, %2i, %4i, %3i)%s%s   +%i   %s  (no: %i)",
 604                           layer,
 605                           ptrA[layer],
 606                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
 607                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
 608                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
 609                           (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
 610                           bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
 611                           ptrB[layer],
 612                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
 613                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
 614                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
 615                           (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
 616                           bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
 617                           inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
 618          ptrA[layer] += inc[layer];
 619          ptrB[layer] += inc[layer];
 620        }
 621
 622      } // end of while
 623    } // end of loop over reflayer
 624
 625    delete [] bHitA;
 626    delete [] bHitB;
 627    delete [] bAligned;
 628    delete [] bDone;
 629    delete [] inc;
 630    delete [] incprime;
 631    delete [] bAlignedA;
 632    delete [] bAlignedB;
 633    delete [] notr;
 634    delete [] ptrA;
 635    delete [] ptrB;
 636
 637    return kTRUE;
 638 }
 639
 640 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
 641 {
 642     TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
 643     TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
 644
 645     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 646         tracksRefMerged[zch] = new TList;
 647         tracksRefUnique[zch] = new TList;
 648     }
 649
 650     TList *tracksZMergedStage0 = new TList;
 651     TList *tracksZUniqueStage0 = new TList;
 652
 653     TList **tracksZSplitted = new TList*[2];
 654     for (Int_t i = 0; i < 2; i++)
 655         tracksZSplitted[i] = new TList;
 656
 657     TList *tracksZMergedStage1 = new TList;
 658
 659     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
 660
 661     Bool_t done = kFALSE;
 662     Int_t minIdx = 0;
 663     AliTRDtrackGTU *trkStage0 = 0x0;
 664
 665     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 666         // ----- Merging and Unification in Reflayers (seed_merger) -----
 667         do {
 668             done = kTRUE;
 669             trkStage0 = 0x0;
 670             for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
 671                 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
 672                 if (trkInRefLayer[refLayerIdx] == 0) {
 673                     continue;
 674                 }
 675                 else if (trkStage0 == 0x0 ) {
 676                     trkStage0 = trkInRefLayer[refLayerIdx];
 677                     minIdx = refLayerIdx;
 678                     done = kFALSE;
 679                 }
 680                 else if (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel() ||
 681                          (trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel() && trkInRefLayer[refLayerIdx]->GetYapprox() < trkStage0->GetYapprox()) ) {
 682                     minIdx = refLayerIdx;
 683                     trkStage0 = trkInRefLayer[refLayerIdx];
 684                     done = kFALSE;
 685                 }
 686             }
 687             if (!trkStage0)
 688               break;
 689             tracksRefMerged[zch]->Add(trkStage0);
 690             fTracks[zch][minIdx].RemoveFirst();
 691         } while (trkStage0 != 0);
 692
 693         Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
 694
 695         AliDebug(2, Form("zchannel %i", zch));
 696         TIter trackRefMerged(tracksRefMerged[zch]);
 697         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
 698           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
 699                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 700                            trk->GetTrackletIndex(5),
 701                            trk->GetTrackletIndex(4),
 702                            trk->GetTrackletIndex(3),
 703                            trk->GetTrackletIndex(2),
 704                            trk->GetTrackletIndex(1),
 705                            trk->GetTrackletIndex(0),
 706                            trk->GetYapprox() >> 3,
 707                            trk->GetZSubChannel()));
 708         AliDebug(2, "uniquified:");
 709         TIter trackRefUnique(tracksRefUnique[zch]);
 710         while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
 711           AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
 712                            AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 713                            trk->GetTrackletIndex(5),
 714                            trk->GetTrackletIndex(4),
 715                            trk->GetTrackletIndex(3),
 716                            trk->GetTrackletIndex(2),
 717                            trk->GetTrackletIndex(1),
 718                            trk->GetTrackletIndex(0),
 719                            trk->GetYapprox() >> 3,
 720                            trk->GetZSubChannel()));
 721     }
 722
 723 // ----- Merging in zchannels - 1st stage -----
 724
 725     do {
 726         done = kTRUE;
 727         trkStage0 = 0x0;
 728         for (Int_t zch = fGtuParam->GetNZChannels() - 1; zch > -1; zch--) {
 729             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
 730             if (trk == 0) {
 731                 continue;
 732             }
 733             else if (trkStage0 == 0x0 ) {
 734                 trkStage0 = trk;
 735                 minIdx = zch;
 736                 done = kFALSE;
 737             }
 738             else if ( ((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) <  ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) ||
 739                       (((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) == ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
 740                 minIdx = zch;
 741                 trkStage0 = trk;
 742                 done = kFALSE;
 743             }
 744         }
 745
 746         if (!trkStage0)
 747           break;
 748         tracksZMergedStage0->Add(trkStage0);
 749         tracksRefUnique[minIdx]->RemoveFirst();
 750     } while (trkStage0 != 0);
 751
 752     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
 753
 754     AliDebug(2, "stage 0:");
 755     TIter trackZMergedStage0(tracksZMergedStage0);
 756     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
 757       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 758                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 759                        trk->GetTrackletIndex(5),
 760                        trk->GetTrackletIndex(4),
 761                        trk->GetTrackletIndex(3),
 762                        trk->GetTrackletIndex(2),
 763                        trk->GetTrackletIndex(1),
 764                        trk->GetTrackletIndex(0),
 765                        trk->GetYapprox() >> 3,
 766                        trk->GetZChannel(),
 767                        trk->GetZSubChannel()));
 768
 769     AliDebug(2, "uniquified:");
 770     TIter trackZUniqueStage0(tracksZUniqueStage0);
 771     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
 772       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 773                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 774                        trk->GetTrackletIndex(5),
 775                        trk->GetTrackletIndex(4),
 776                        trk->GetTrackletIndex(3),
 777                        trk->GetTrackletIndex(2),
 778                        trk->GetTrackletIndex(1),
 779                        trk->GetTrackletIndex(0),
 780                        trk->GetYapprox() >> 3,
 781                        trk->GetZChannel(),
 782                        trk->GetZSubChannel()));
 783
 784 // ----- Splitting in z -----
 785
 786     TIter next(tracksZUniqueStage0);
 787     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
 788         tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
 789     }
 790
 791     for (Int_t i = 0; i < 2; i++) {
 792       AliDebug(2, Form("split %i:", i));
 793       TIter trackZSplit(tracksZSplitted[i]);
 794       while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
 795         AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 796                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 797                          trk->GetTrackletIndex(5),
 798                          trk->GetTrackletIndex(4),
 799                          trk->GetTrackletIndex(3),
 800                          trk->GetTrackletIndex(2),
 801                          trk->GetTrackletIndex(1),
 802                          trk->GetTrackletIndex(0),
 803                          trk->GetYapprox() >> 3,
 804                          trk->GetZChannel(),
 805                          trk->GetZSubChannel()));
 806     }
 807
 808 // ----- Merging in zchanels - 2nd stage -----
 809
 810     do {
 811         done = kTRUE;
 812         trkStage0 = 0x0;
 813         for (Int_t i = 1; i >= 0; i--) {
 814             AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
 815             if (trk == 0) {
 816                 continue;
 817             }
 818             else if (trkStage0 == 0x0 ) {
 819                 trkStage0 = trk;
 820                 minIdx = i;
 821                 done = kFALSE;
 822             }
 823             else if ( (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) <  ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
 824                       ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) && (trk->GetYapprox() < trkStage0->GetYapprox()) ) ) {
 825                 minIdx = i;
 826                 trkStage0 = trk;
 827                 done = kFALSE;
 828             }
 829         }
 830
 831         if (!trkStage0)
 832           break;
 833         tracksZMergedStage1->Add(trkStage0);
 834         tracksZSplitted[minIdx]->RemoveFirst();
 835     } while (trkStage0 != 0);
 836
 837     Uniquifier(tracksZMergedStage1, ListOfTracks);
 838
 839     AliDebug(2, "merged:");
 840     TIter trackZMergedStage1(tracksZMergedStage1);
 841     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
 842       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 843                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 844                        trk->GetTrackletIndex(5),
 845                        trk->GetTrackletIndex(4),
 846                        trk->GetTrackletIndex(3),
 847                        trk->GetTrackletIndex(2),
 848                        trk->GetTrackletIndex(1),
 849                        trk->GetTrackletIndex(0),
 850                        trk->GetYapprox() >> 3,
 851                        trk->GetZChannel(),
 852                        trk->GetZSubChannel()));
 853
 854     AliDebug(2, "uniquified:");
 855     TIter track(ListOfTracks);
 856     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
 857       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
 858                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
 859                        trk->GetTrackletIndex(5),
 860                        trk->GetTrackletIndex(4),
 861                        trk->GetTrackletIndex(3),
 862                        trk->GetTrackletIndex(2),
 863                        trk->GetTrackletIndex(1),
 864                        trk->GetTrackletIndex(0),
 865                        trk->GetYapprox() >> 3,
 866                        trk->GetZChannel(),
 867                        trk->GetZSubChannel()));
 868
 869     // cleaning up
 870     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
 871       delete tracksRefMerged[zch];
 872       delete tracksRefUnique[zch];
 873     }
 874     delete [] tracksRefMerged;
 875     delete [] tracksRefUnique;
 876
 877
 878     delete tracksZMergedStage0;
 879     delete tracksZUniqueStage0;
 880
 881     for (Int_t i = 0; i < 2; i++)
 882       delete tracksZSplitted[i];
 883     delete [] tracksZSplitted;
 884
 885     delete tracksZMergedStage1;
 886
 887     delete [] trkInRefLayer;
 888
 889     return kTRUE;
 890 }
 891
 892 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
 893 {
 894   // run the track reconstruction for all tracks in the list
 895
 896   TIter next(ListOfTracks);
 897   while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
 898     CalculateTrackParams(track);
 899     CalculatePID(track);
 900     AliDebug(1, Form("track with pid = %i", track->GetPID()));
 901   }
 902   return kTRUE;
 903 }
 904
 905 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
 906 {
 907   // calculate PID for the given track
 908   if (!track) {
 909     AliError("No track to calculate!");
 910     return kFALSE;
 911   }
 912
 913   if (AliTRDgtuParam::GetUseGTUconst()) {
 914     // averaging as in GTU
 915     ULong64_t coeff;
 916
 917     // selection of coefficient for averaging
 918     switch(track->GetTrackletMask()){
 919     case 0x3f:
 920       // 6 tracklets
 921       coeff=0x5558; // approx. 1/6
 922       break;
 923
 924     case 0x3e:
 925     case 0x3d:
 926     case 0x3b:
 927     case 0x37:
 928     case 0x2f:
 929     case 0x1f:
 930       // 5 tracklets
 931       coeff=0x6666; // approx. 1/5
 932       break;
 933
 934     default:
 935       // 4 tracklets
 936       coeff=0x8000; // = 0.25
 937     }
 938     coeff &= 0x1ffff; // 17-bit constant
 939
 940     ULong64_t sum = 0;
 941     for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
 942       if ((track->GetTrackletMask() >> iLayer) & 1) {
 943         sum += track->GetTracklet(iLayer)->GetPID();
 944       }
 945     }
 946
 947     sum = sum & 0x7ff;
 948     ULong64_t prod   = (sum * coeff) & 0xfffffffffull;
 949     ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
 950
 951     track->SetPID(prodFinal & 0xff);
 952
 953     return kTRUE;
 954   }
 955   else {
 956
 957     // simple averaging
 958     Int_t nTracklets = 0;
 959     Int_t pidSum = 0;
 960     for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 961       if (!track->IsTrackletInLayer(layer)) {
 962         continue;
 963       }
 964       AliTRDtrackletGTU *trk = track->GetTracklet(layer);
 965       pidSum += trk->GetPID();
 966       nTracklets++;
 967     }
 968
 969     if (nTracklets>0)
 970       track->SetPID(pidSum/nTracklets);
 971     else
 972       AliError("Track without contributing tracklets, no PID assigned");
 973
 974     return kTRUE;
 975   }
 976 }
 977
 978 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
 979 {
 980   // calculate the track parameters
 981
 982   if (!track) {
 983     AliError("No track to calculate!");
 984     return kFALSE;
 985   }
 986
 987   Int_t a = 0;
 988   Float_t b = 0;
 989   Float_t c = 0;
 990   Float_t x1;
 991   Float_t x2;
 992
 993   AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
 994
 995   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
 996     if (!track->IsTrackletInLayer(layer)) {
 997       continue;
 998     }
 999     AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1000     if (!trk) {
1001       AliError(Form("Could not get tracklet in layer %i\n", layer));
1002       continue;
1003     }
1004     AliDebug(10,Form("  layer %i trk yprime: %6i, aki: %6i", layer, trk->GetYPrime(), (Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))));
1005     a += (((Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))) * trk->GetYPrime() + 1) >> 8;
1006     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1007     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1008   }
1009   if (a < 0)
1010       a += 3;
1011   a = a >> 2;
1012
1013   track->SetFitParams(a, b, c);
1014
1015   fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1016
1017   AliDebug(5,Form("  Track parameters: a = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f (trkl mask: %i)", a, b, c, x1, x2, track->GetPt(), track->GetTrackletMask()));
1018
1019   return kTRUE;
1020 }
1021
1022
1023 Bool_t AliTRDgtuTMU::Uniquifier(const TList *inlist, TList *outlist)
1024 {
1025   // remove multiple occurences of the same track
1026
1027     TIter next(inlist);
1028     AliTRDtrackGTU *trkStage0 = 0x0;
1029     AliTRDtrackGTU *trkStage1 = 0x0;
1030
1031     do {
1032         trkStage0 = (AliTRDtrackGTU*) next();
1033
1034         Bool_t tracksEqual = kFALSE;
1035         if (trkStage0 != 0 && trkStage1 != 0) {
1036             for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1037                 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1038                     tracksEqual = kTRUE;
1039                     break;
1040                 }
1041             }
1042         }
1043
1044         if (tracksEqual) {
1045             if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
1046                 trkStage1 = trkStage0;
1047         }
1048         else {
1049             if (trkStage1 != 0x0)
1050                 outlist->Add(trkStage1);
1051             trkStage1 = trkStage0;
1052         }
1053
1054     } while (trkStage1 != 0x0);
1055     return kTRUE;
1056 }