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