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