]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDgtuTMU.cxx
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuTMU.cxx
index 661382db81a5b892685b2e54b51a2051472c9ffe..abb69af89dc314ed6f5c04012fe9382eae14631d 100644 (file)
@@ -46,6 +46,7 @@ AliTRDgtuTMU::AliTRDgtuTMU(Int_t stack, Int_t sector) :
   fTracklets(0x0),
   fTrackletsPostInput(0x0),
   fZChannelTracklets(0x0),
+  fTrackArray(new TClonesArray("AliTRDtrackGTU", 50)),
   fTracks(0x0),
   fGtuParam(0x0),
   fStack(-1),
@@ -88,12 +89,14 @@ AliTRDgtuTMU::~AliTRDgtuTMU()
     delete [] fTracks[zch];
   }
   delete [] fTracks;
+  delete fTrackArray;
+
   for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
     delete [] fZChannelTracklets[layer];
-    //    delete [] fTrackletsPostInput[layer];
+    delete fTrackletsPostInput[layer];
   }
   delete [] fZChannelTracklets;
-  //  delete [] fTrackletsPostInput;
+  delete [] fTrackletsPostInput;
 
   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
     delete fTracklets[iLink];
@@ -137,6 +140,8 @@ Bool_t AliTRDgtuTMU::Reset()
     }
   }
 
+  fTrackArray->Delete();
+
   // delete all added tracklets
   for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
     fTracklets[iLink]->Clear();
@@ -162,7 +167,7 @@ Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletGTU *tracklet, Int_t link)
 }
 
 
-Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
+Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd, Int_t outLabel)
 {
   // performs the analysis as in a TMU module of the GTU, i. e.
   // track matching
@@ -185,7 +190,6 @@ Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
   // ----- Z-channel units -----
   AliDebug(1,"--------- Running Z-channel units ----------");
   for (Int_t layer = 0;  layer <  fGtuParam->GetNLayers(); layer++) {
-    fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
     if (!RunZChannelUnit(layer)) {
       AliError(Form("Z-Channel unit in layer %i failed", layer));
       return kFALSE;
@@ -195,7 +199,6 @@ Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
   // ----- track finding -----
   AliDebug(1,"--------- Running tracking units ----------");
   for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
-    AliDebug(2,Form("Track finder for Zchannel: %i", zch));
     if (!RunTrackFinder(zch, ListOfTracks)) {
       AliError(Form("Track Finder in z-channel %i failed", zch));
       return kFALSE;
@@ -214,13 +217,18 @@ Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd)
     return kFALSE;
   }
 
-  if (esd) {
-      TIter next(ListOfTracks);
-      while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
-         AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
-         esd->AddTrdTrack(trdtrack);
-         delete trdtrack;
-      }
+  // ----- label calculation and ESD storage -----
+  TIter next(ListOfTracks);
+  while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
+    if (outLabel == -1)
+      trk->CookLabel();
+    else
+      trk->SetLabel(outLabel);
+    if (esd) {
+      AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
+      esd->AddTrdTrack(trdtrack);
+      delete trdtrack;
+    }
   }
 
   return kTRUE;
@@ -233,7 +241,11 @@ Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
   Int_t iTrkl0 = 0; // A-side tracklet
   Int_t iTrkl1 = 0; // B-side tracklet
 
-  while ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) || (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast())) {
+  Int_t nTracklets = 0;
+  while ((!fGtuParam->GetLimitNoTracklets() ||
+         (nTracklets < fGtuParam->GetMaxNoTracklets())) &&
+        ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) ||
+         (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast()))) {
     if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
       fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
       iTrkl0++;
@@ -256,6 +268,7 @@ Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
        }
       }
     }
+    ++nTracklets;
   }
 
   TIter next(fTrackletsPostInput[layer]);
@@ -265,17 +278,25 @@ Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
 
     Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
     alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
-    trk->SetAlpha(alpha);
-
-    Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer)); //??? sign?
+    // wrapping projected alpha as in hardware
+    if ((alpha < -64) || (alpha > 63))
+      AliDebug(1, Form("alpha out of range: %i", alpha));
+    alpha = alpha & 0x7f;
+    if (alpha & 0x40)
+      trk->SetAlpha(0xffffffc0 | alpha);
+    else
+      trk->SetAlpha(alpha);
+
+    Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer));
     yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
     trk->SetYProj(yproj);
 
     trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
 
-    AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i",
+    AliDebug(10, Form("0x%08x: idx: %3i, z: %2i, y: %5i, dy: %3i, y': %5i, y_proj: %5i, alpha: %3i, pid: %3i, c: %5i, yt: %5i",
                      trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
-                     trk->GetYProj(), trk->GetAlpha(), trk->GetPID() ));
+                     trk->GetYProj(), trk->GetAlpha(), trk->GetPID(),
+                     fGtuParam->GetCiYProj(layer), fGtuParam->GetYt(fStack, layer, trk->GetZbin()) ));
   }
   return kTRUE;
 }
@@ -359,7 +380,7 @@ Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
 
    for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
      Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
-     AliDebug(5,Form("~~~~~ Reflayer: %i", reflayer));
+     AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
 
      ready = kFALSE; // ready if all channels done
 
@@ -406,7 +427,6 @@ Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
        if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
         trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
 
-       AliDebug(10,Form("ptrRA: %i, ptrRB: %i", ptrA[reflayer], ptrB[reflayer]));
        yPlus     = trkRA->GetYProj() + fGtuParam->GetDeltaY();
        yMinus    = trkRA->GetYProj() - fGtuParam->GetDeltaY();
        alphaPlus  = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
@@ -507,12 +527,12 @@ Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
         }
        } // end of loop over layers
 
-       AliDebug(5,Form("logic calculation finished, Nhits: %i", nHits));
+       AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
+                      nHits, (nHits >= 4) ? "-> track found" : ""));
 
        if (nHits >= 4) {
         // ----- track registration -----
-        AliDebug(1,"***** TMU: Track found *****");
-        AliTRDtrackGTU *track = new AliTRDtrackGTU();
+        AliTRDtrackGTU *track = new ((*fTrackArray)[fTrackArray->GetEntriesFast()]) AliTRDtrackGTU();
         track->SetSector(fSector);
         track->SetStack(fStack);
         for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
@@ -584,7 +604,7 @@ Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
        }
 
        ready = kTRUE;
-       for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
+       for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
 
         bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
         ready = ready && bDone[layer];
@@ -593,13 +613,19 @@ Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
           AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
 
 //      AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
-        AliDebug(10,Form(" -- Layer: %i   %2i(%2i)%s%s   %2i(%2i)%s%s   +%i   %s  (no: %i)",
+        AliDebug(10,Form("    Layer: %i   %2i(%2i, %2i, %4i, %3i)%s%s   %2i(%2i, %2i, %4i, %3i)%s%s   +%i   %s  (no: %i)",
                          layer,
                          ptrA[layer],
                          (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
+                         (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
+                         (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
+                         (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
                          bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
                          ptrB[layer],
                          (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
+                         (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
+                         (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
+                         (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
                          bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
                          inc[layer], bDone[layer] ? "done" : "    ", notr[layer]));
         ptrA[layer] += inc[layer];
@@ -645,14 +671,14 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
 
     AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
 
-    Bool_t done = kFALSE;
+    // Bool_t done = kFALSE;
     Int_t minIdx = 0;
     AliTRDtrackGTU *trkStage0 = 0x0;
 
     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
         // ----- Merging and Unification in Reflayers (seed_merger) -----
        do {
-           done = kTRUE;
+         // done = kTRUE;
            trkStage0 = 0x0;
            for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
                trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
@@ -662,13 +688,14 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
                else if (trkStage0 == 0x0 ) {
                    trkStage0 = trkInRefLayer[refLayerIdx];
                    minIdx = refLayerIdx;
-                   done = kFALSE;
+                   // done = kFALSE;
                }
-               else if (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel() ||
-                        (trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel() && trkInRefLayer[refLayerIdx]->GetYapprox() < trkStage0->GetYapprox()) ) {
+               else if ( (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel()) ||
+                         ((trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel()) &&
+                          ((trkInRefLayer[refLayerIdx]->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
                    minIdx = refLayerIdx;
                    trkStage0 = trkInRefLayer[refLayerIdx];
-                   done = kFALSE;
+                   // done = kFALSE;
                }
            }
             if (!trkStage0)
@@ -683,7 +710,7 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
        TIter trackRefMerged(tracksRefMerged[zch]);
        while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
          AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
-                          trk->GetRefLayerIdx(),
+                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
                           trk->GetTrackletIndex(5),
                           trk->GetTrackletIndex(4),
                           trk->GetTrackletIndex(3),
@@ -696,7 +723,7 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
        TIter trackRefUnique(tracksRefUnique[zch]);
        while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
          AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
-                          trk->GetRefLayerIdx(),
+                          AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
                           trk->GetTrackletIndex(5),
                           trk->GetTrackletIndex(4),
                           trk->GetTrackletIndex(3),
@@ -709,45 +736,216 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
 
 // ----- Merging in zchannels - 1st stage -----
 
-    do {
-       done = kTRUE;
-       trkStage0 = 0x0;
-       for (Int_t zch = fGtuParam->GetNZChannels() - 1; zch > -1; zch--) {
-           AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
-           if (trk == 0) {
-               continue;
-           }
-           else if (trkStage0 == 0x0 ) {
-               trkStage0 = trk;
-               minIdx = zch;
-               done = kFALSE;
+    if (AliTRDgtuParam::GetUseGTUmerge()) {
+      Int_t notEmpty;
+      do {
+       Bool_t lowerThan[3] = { kFALSE, kFALSE, kFALSE };
+       AliTRDtrackGTU *trk[3] = { 0x0, 0x0, 0x0 };
+       for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel)
+         trk[iChannel] = (AliTRDtrackGTU*) tracksRefUnique[iChannel]->First();
+
+       for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel) {
+         AliTRDtrackGTU *trk1 = trk[iChannel];
+         AliTRDtrackGTU *trk2 = trk[(iChannel + 1) % 3];
+         if (trk1 && trk2) {
+           Int_t sortnum1 = (trk1->GetZChannel() + 3 * trk1->GetZSubChannel()) / 2 - 1;
+           Int_t sortnum2 = (trk2->GetZChannel() + 3 * trk2->GetZSubChannel()) / 2 - 1;
+           AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
+                            trk1->GetZChannel(), trk2->GetZChannel(),
+                            sortnum1, sortnum2));
+           if ( (sortnum1 < sortnum2) ||
+                ((sortnum1 == sortnum2) &&
+                 ((trk1->GetYapprox() >> 3) < (trk2->GetYapprox() >> 3)) ) ) {
+             lowerThan[iChannel] = kTRUE;
            }
-           else if ( ((trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1) < ((trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 -1 ) ) {
-               minIdx = zch;
-               trkStage0 = trk;
-               done = kFALSE;
+         }
+       }
+
+       notEmpty = (trk[2] ? (1 << 2) : 0) |
+         (trk[1] ? (1 << 1) : 0) |
+         (trk[0] ? (1 << 0) : 0);
+       Int_t pop = -1;
+
+       switch (notEmpty) {
+         // one track only
+       case 0x1:
+         pop = 0;
+         break;
+       case 0x2:
+         pop = 1;
+         break;
+       case 0x4:
+         pop = 2;
+         break;
+
+         // two tracks
+       case 0x3:
+         if (lowerThan[0])
+           pop = 0;
+         else
+           pop = 1;
+         break;
+       case 0x5:
+         if (lowerThan[2])
+           pop = 2;
+         else
+           pop = 0;
+         break;
+       case 0x6:
+         if (lowerThan[1])
+           pop = 1;
+         else
+           pop = 2;
+         break;
+
+         // three tracks
+       case 0x7:
+         if (lowerThan[0]) {
+           if (lowerThan[2])
+             pop = 2;
+           else
+             pop = 0;
+         } else {
+           if (lowerThan[1])
+             pop = 1;
+           else
+             pop = 2;
+         }
+         break;
+
+         // no tracks
+       default:
+         // nop
+         ;
+       }
+
+       if (pop > -1) {
+         tracksZMergedStage0->Add(trk[pop]);
+         tracksRefUnique[pop]->RemoveFirst();
+       }
+      } while (notEmpty);
+    }
+    else {
+      // there is still a difference between this implementation and
+      // the hardware algorithm, only for expert use
+
+      do {
+       // done = kTRUE;
+       trkStage0 = 0x0;
+       // compare tracks from all adjacent zchannels
+       // (including comparison of channels 2 and 0)
+        for (Int_t i = fGtuParam->GetNZChannels() - 1; i > -1; i--) {
+         Int_t zch = i % 3;
+         AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
+         if (trk == 0) {
+           continue;
+         }
+         else if (trkStage0 == 0x0 ) {
+           trkStage0 = trk;
+           minIdx = zch;
+           // done = kFALSE;
+         }
+         else {
+           Int_t sortnum1 = (trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1;
+           Int_t sortnum2 = (trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 - 1;
+           AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
+                            trk->GetZChannel(), trkStage0->GetZChannel(),
+                            sortnum1, sortnum2));
+           if ( (sortnum1 < sortnum2) ||
+                ((sortnum1 == sortnum2) &&
+                 ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
+             minIdx = zch;
+             trkStage0 = trk;
+             // done = kFALSE;
            }
+         }
        }
 
        if (!trkStage0)
           break;
        tracksZMergedStage0->Add(trkStage0);
        tracksRefUnique[minIdx]->RemoveFirst();
-    } while (trkStage0 != 0);
+      } while (trkStage0 != 0);
+    }
 
     Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
 
+    AliDebug(2, "stage 0:");
+    TIter trackZMergedStage0(tracksZMergedStage0);
+    while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
+      AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
+                      AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
+                      trk->GetTrackletIndex(5),
+                      trk->GetTrackletIndex(4),
+                      trk->GetTrackletIndex(3),
+                      trk->GetTrackletIndex(2),
+                      trk->GetTrackletIndex(1),
+                      trk->GetTrackletIndex(0),
+                      trk->GetYapprox() >> 3,
+                      trk->GetZChannel(),
+                      trk->GetZSubChannel()));
+
+    AliDebug(2, "uniquified:");
+    TIter trackZUniqueStage0(tracksZUniqueStage0);
+    while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
+      AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
+                      AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
+                      trk->GetTrackletIndex(5),
+                      trk->GetTrackletIndex(4),
+                      trk->GetTrackletIndex(3),
+                      trk->GetTrackletIndex(2),
+                      trk->GetTrackletIndex(1),
+                      trk->GetTrackletIndex(0),
+                      trk->GetYapprox() >> 3,
+                      trk->GetZChannel(),
+                      trk->GetZSubChannel()));
+
 // ----- Splitting in z -----
 
     TIter next(tracksZUniqueStage0);
     while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
-       tracksZSplitted[(trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2]->Add(trk);
+      if ((trk->GetZChannel() < 0) ||
+         (trk->GetZChannel() > 2) ||
+         (trk->GetZSubChannel() < 0) ||
+         (trk->GetZSubChannel() > 6)) {
+       AliError(Form("track with invalid z-channel information at %p: zch = %i, subchannel = %i",
+                     trk, trk->GetZChannel(), trk->GetZSubChannel()));
+       trk->Dump();
+      }
+      Int_t idx = (trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2;
+      if ((idx < 0) || (idx > 1)) {
+       AliError(Form("invalid index %i null", idx));
+       trk->Dump();
+       continue;
+      }
+      if (!tracksZSplitted[idx]) {
+       AliError(Form("array pointer %i null", idx));
+       continue;
+      }
+      tracksZSplitted[idx]->Add(trk);
+    }
+
+    for (Int_t i = 0; i < 2; i++) {
+      AliDebug(2, Form("split %i:", i));
+      TIter trackZSplit(tracksZSplitted[i]);
+      while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
+       AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
+                        AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
+                        trk->GetTrackletIndex(5),
+                        trk->GetTrackletIndex(4),
+                        trk->GetTrackletIndex(3),
+                        trk->GetTrackletIndex(2),
+                        trk->GetTrackletIndex(1),
+                        trk->GetTrackletIndex(0),
+                        trk->GetYapprox() >> 3,
+                        trk->GetZChannel(),
+                        trk->GetZSubChannel()));
     }
 
 // ----- Merging in zchanels - 2nd stage -----
 
     do {
-       done = kTRUE;
+      // done = kTRUE;
        trkStage0 = 0x0;
        for (Int_t i = 1; i >= 0; i--) {
            AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
@@ -757,12 +955,14 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
            else if (trkStage0 == 0x0 ) {
                trkStage0 = trk;
                minIdx = i;
-               done = kFALSE;
+               // done = kFALSE;
            }
-           else if ( ((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) < ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2) ) {
+           else if (  (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) <  ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
+                     ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) &&
+                      ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
                minIdx = i;
                trkStage0 = trk;
-               done = kFALSE;
+               // done = kFALSE;
            }
        }
 
@@ -774,6 +974,36 @@ Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
 
     Uniquifier(tracksZMergedStage1, ListOfTracks);
 
+    AliDebug(2, "merged:");
+    TIter trackZMergedStage1(tracksZMergedStage1);
+    while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
+      AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
+                      AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
+                      trk->GetTrackletIndex(5),
+                      trk->GetTrackletIndex(4),
+                      trk->GetTrackletIndex(3),
+                      trk->GetTrackletIndex(2),
+                      trk->GetTrackletIndex(1),
+                      trk->GetTrackletIndex(0),
+                      trk->GetYapprox() >> 3,
+                      trk->GetZChannel(),
+                      trk->GetZSubChannel()));
+
+    AliDebug(2, "uniquified:");
+    TIter track(ListOfTracks);
+    while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
+      AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
+                      AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
+                      trk->GetTrackletIndex(5),
+                      trk->GetTrackletIndex(4),
+                      trk->GetTrackletIndex(3),
+                      trk->GetTrackletIndex(2),
+                      trk->GetTrackletIndex(1),
+                      trk->GetTrackletIndex(0),
+                      trk->GetYapprox() >> 3,
+                      trk->GetZChannel(),
+                      trk->GetZSubChannel()));
+
     // cleaning up
     for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
       delete tracksRefMerged[zch];
@@ -820,6 +1050,7 @@ Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
 
   if (AliTRDgtuParam::GetUseGTUconst()) {
     // averaging as in GTU
+    AliDebug(1, "using GTU constants for PID calculation");
     ULong64_t coeff;
 
     // selection of coefficient for averaging
@@ -846,16 +1077,21 @@ Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
     coeff &= 0x1ffff; // 17-bit constant
 
     ULong64_t sum = 0;
+    Int_t i = 0;
     for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
       if ((track->GetTrackletMask() >> iLayer) & 1) {
        sum += track->GetTracklet(iLayer)->GetPID();
+       ++i;
       }
     }
 
-    ULong64_t sumExt = (sum << 1) & 0xfff; // 11 bit -> 12 bit vector
-    ULong64_t prod   = (sumExt * coeff) & 0xfffffffffull; // 18x18 signed -> 36
-    ULong64_t prodFinal = ((prod >> 18) + ((prod >> 17) & 1)) & 0xff; // rounding term is equivalent to adding 5 to sum_ext
+    Float_t av = 1./i * sum;
+    sum = sum & 0x7ff;
+    ULong64_t prod   = (sum * coeff) & 0xfffffffffull;
+    ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
 
+    if (TMath::Abs((prodFinal & 0xff) - av) > 0.5)
+      AliError(Form("problem with PID averaging (hw <-> ar): %3lli <-> %4.1f", prodFinal & 0xff, av));
     track->SetPID(prodFinal & 0xff);
 
     return kTRUE;
@@ -873,7 +1109,11 @@ Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
       pidSum += trk->GetPID();
       nTracklets++;
     }
-    track->SetPID(pidSum/nTracklets);
+
+    if (nTracklets>0)
+      track->SetPID(pidSum/nTracklets);
+    else
+      AliError("Track without contributing tracklets, no PID assigned");
 
     return kTRUE;
   }
@@ -905,8 +1145,9 @@ Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
       AliError(Form("Could not get tracklet in layer %i\n", layer));
       continue;
     }
-    AliDebug(10,Form("trk yprime: %i", trk->GetYPrime()));
-    a += (((Int_t) (2048 * fGtuParam->GetAki(track->GetTrackletMask(), layer))) * trk->GetYPrime() + 1) >> 8;
+    AliDebug(10,Form("  layer %i trk yprime: %6i, aki: %6i", layer, trk->GetYPrime(),
+                    fGtuParam->GetAki(track->GetTrackletMask(), layer)));
+    a += (((fGtuParam->GetAki(track->GetTrackletMask(), layer) * trk->GetYPrime()) >> 7) + 1) >> 1;
     b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
     c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
   }
@@ -914,42 +1155,18 @@ Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
       a += 3;
   a = a >> 2;
 
-  fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
-  AliDebug(10,Form("Intersection points: %f, %f", x1, x2));
-  AliDebug(10,Form("Sum: a = %5i, b = %9.2f, c = %9.2f\n", a, b, c));
-  track->SetFitParams(a, b, c);
+  track->SetFitParams(a << 2, b, c);
 
-  Int_t pt = fGtuParam->GetPt(track->GetTrackletMask(), a, b, x1, x2);
-  track->SetPtInt(pt);
-  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()));
-  return kTRUE;
-}
+  fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
 
-Bool_t AliTRDgtuTMU::WriteTrackletsToTree(TTree *trklTree)
-{
-  if (!trklTree) {
-    AliError("No tree given");
-    return kFALSE;
-  }
-  AliTRDtrackletGTU *trkl = 0x0;
-  TBranch *branch = trklTree->GetBranch("gtutracklets");
-  if (!branch) {
-      branch = trklTree->Branch("gtutracklets", "AliTRDtrackletGTU", &trkl, 32000, 99);
-  }
+  AliDebug(5,Form("  Track parameters: a16 = %i, a18 = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f (trkl mask: %i)",
+                 a, a << 2, b, c, x1, x2, track->GetPt(), track->GetTrackletMask()));
 
-  AliDebug(5,Form("---------- Writing tracklets to tree (not yet) ----------"));
-  for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
-    TIter next(fTrackletsPostInput[layer]);
-    while ((trkl = (AliTRDtrackletGTU*) next())) {
-       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) ));
-       branch->SetAddress(&trkl);
-       trklTree->Fill();
-    }
-  }
   return kTRUE;
 }
 
-Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
+
+Bool_t AliTRDgtuTMU::Uniquifier(const TList *inlist, TList *outlist)
 {
   // remove multiple occurences of the same track
 
@@ -958,7 +1175,7 @@ Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
     AliTRDtrackGTU *trkStage1 = 0x0;
 
     do {
-       trkStage0 = (AliTRDtrackGTU*) next();
+        trkStage0 = (AliTRDtrackGTU*) next();
 
        Bool_t tracksEqual = kFALSE;
        if (trkStage0 != 0 && trkStage1 != 0) {
@@ -971,8 +1188,8 @@ Bool_t AliTRDgtuTMU::Uniquifier(TList *inlist, TList *outlist)
        }
 
        if (tracksEqual) {
-           if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
-               trkStage1 = trkStage0;
+         if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
+           trkStage1 = trkStage0;
        }
        else {
            if (trkStage1 != 0x0)