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