TRD/AliTRDgtuTMU.cxx

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