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