1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDgtuTMU.cxx 28397 2008-09-02 09:33:00Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Track Matching Unit (TMU) simulation //
22 // Author: J. Klein (Jochen.Klein@cern.ch) //
24 ////////////////////////////////////////////////////////////////////////////
31 #include "AliESDEvent.h"
32 #include "AliESDTrdTrack.h"
35 #include "AliTRDgeometry.h"
36 #include "AliTRDpadPlane.h"
38 #include "AliTRDgtuParam.h"
39 #include "AliTRDgtuTMU.h"
40 #include "AliTRDtrackGTU.h"
42 ClassImp(AliTRDgtuTMU)
44 AliTRDgtuTMU::AliTRDgtuTMU(Int_t stack, Int_t sector) :
47 fTrackletsPostInput(0x0),
48 fZChannelTracklets(0x0),
49 fTrackArray(new TClonesArray("AliTRDtrackGTU", 50)),
55 // constructor which initializes the position information of the TMU
57 fGtuParam = AliTRDgtuParam::Instance();
59 // store tracklets per link
60 fTracklets = new TObjArray*[fGtuParam->GetNLinks()];
61 for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
62 fTracklets[iLink] = new TObjArray();
65 // tracklets after input units per layer
66 fTrackletsPostInput = new TObjArray*[fGtuParam->GetNLayers()];
67 fZChannelTracklets = new TList*[fGtuParam->GetNLayers()];
68 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
69 fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
70 fTrackletsPostInput[layer] = new TObjArray();
73 fTracks = new TList*[fGtuParam->GetNZChannels()];
74 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
75 fTracks[zch] = new TList[fGtuParam->GetNRefLayers()];
84 AliTRDgtuTMU::~AliTRDgtuTMU()
88 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
89 delete [] fTracks[zch];
94 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
95 delete [] fZChannelTracklets[layer];
96 delete fTrackletsPostInput[layer];
98 delete [] fZChannelTracklets;
99 delete [] fTrackletsPostInput;
101 for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
102 delete fTracklets[iLink];
104 delete [] fTracklets;
107 Bool_t AliTRDgtuTMU::SetSector(Int_t sector)
111 if (sector > -1 && sector < fGtuParam->GetGeo()->Nsector() ) {
116 AliError(Form("Invalid sector given: %i", sector));
120 Bool_t AliTRDgtuTMU::SetStack(Int_t stack)
122 // Set the stack (necessary for tracking)
124 if (stack > -1 && stack < fGtuParam->GetGeo()->Nstack() ) {
129 AliError(Form("Invalid stack given: %i", stack));
133 Bool_t AliTRDgtuTMU::Reset()
137 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
138 for (Int_t reflayeridx = 0; reflayeridx < fGtuParam->GetNRefLayers(); reflayeridx++) {
139 fTracks[zch][reflayeridx].Clear();
143 fTrackArray->Delete();
145 // delete all added tracklets
146 for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
147 fTracklets[iLink]->Clear();
149 for (Int_t layer = 0; layer < fGtuParam->GetNLinks()/2; layer++) {
150 fTrackletsPostInput[layer]->Clear();
151 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++)
152 fZChannelTracklets[layer][zch].Clear();
161 Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletGTU *tracklet, Int_t link)
163 // add a tracklet on the given link
165 fTracklets[link]->Add(tracklet);
170 Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd, Int_t outLabel)
172 // performs the analysis as in a TMU module of the GTU, i. e.
174 // calculation of track parameteres (pt, deflection, ???)
176 if (fStack < 0 || fSector < 0) {
177 AliError("No valid stack/sector set for this TMU! No tracking!");
181 // ----- Input units -----
182 AliDebug(1,"--------- Running Input units ----------");
183 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
184 if (!RunInputUnit(layer)) {
185 AliError(Form("Input unit in layer %i failed", layer));
190 // ----- Z-channel units -----
191 AliDebug(1,"--------- Running Z-channel units ----------");
192 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
193 if (!RunZChannelUnit(layer)) {
194 AliError(Form("Z-Channel unit in layer %i failed", layer));
199 // ----- track finding -----
200 AliDebug(1,"--------- Running tracking units ----------");
201 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
202 if (!RunTrackFinder(zch, ListOfTracks)) {
203 AliError(Form("Track Finder in z-channel %i failed", zch));
208 // ----- Track Merging -----
209 if (!RunTrackMerging(ListOfTracks)) {
210 AliError("Track merging failed");
214 // ----- track reconstruction -----
215 if (!RunTrackReconstruction(ListOfTracks)) {
216 AliError("Track reconstruction failed");
220 // ----- label calculation and ESD storage -----
221 TIter next(ListOfTracks);
222 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
226 trk->SetLabel(outLabel);
228 AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
229 esd->AddTrdTrack(trdtrack);
237 Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
239 // precalculations for the tracking and reconstruction
241 Int_t iTrkl0 = 0; // A-side tracklet
242 Int_t iTrkl1 = 0; // B-side tracklet
244 Int_t nTracklets = 0;
245 while ((!fGtuParam->GetLimitNoTracklets() ||
246 (nTracklets < fGtuParam->GetMaxNoTracklets())) &&
247 ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) ||
248 (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast()))) {
249 if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
250 fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
254 if (iTrkl0 >= fTracklets[2*layer + 0]->GetEntriesFast()) {
255 fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
259 if (((AliTRDtrackletGTU*) fTracklets[2*layer + 1]->At(iTrkl1))->GetZbin() <
260 ((AliTRDtrackletGTU*) fTracklets[2*layer + 0]->At(iTrkl0))->GetZbin()) {
261 fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
266 fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
274 TIter next(fTrackletsPostInput[layer]);
276 while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
277 trk->SetIndex(fTrackletsPostInput[layer]->IndexOf(trk));
279 Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
280 alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
281 // wrapping projected alpha as in hardware
282 if ((alpha < -64) || (alpha > 63))
283 AliDebug(1, Form("alpha out of range: %i", alpha));
284 alpha = alpha & 0x7f;
286 trk->SetAlpha(0xffffffc0 | alpha);
288 trk->SetAlpha(alpha);
290 Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer));
291 yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
292 trk->SetYProj(yproj);
294 trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
296 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",
297 trk->GetTrackletWord(), trk->GetIndex(), trk->GetZbin(), trk->GetYbin(), trk->GetdY(), trk->GetYPrime(),
298 trk->GetYProj(), trk->GetAlpha(), trk->GetPID(),
299 fGtuParam->GetCiYProj(layer), fGtuParam->GetYt(fStack, layer, trk->GetZbin()) ));
304 Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
306 // run the z-channel unit
308 TIter next(fTrackletsPostInput[layer]);
310 while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
311 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
312 if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
313 trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
315 TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
316 AliTRDtrackletGTU *t = 0x0;
317 while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
318 if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
319 ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
324 fZChannelTracklets[layer][zch].AddAfter(t, trk);
326 fZChannelTracklets[layer][zch].AddFirst(trk);
329 AliDebug(10, Form("stack %d, layer %2d: 0x%08x Z(2..0)=%i/%i/%i",
330 fStack, layer, trk->GetTrackletWord(),
331 fGtuParam->IsInZChannel(fStack, layer, 2, trk->GetZbin()) ? trk->GetSubChannel(2) : -1,
332 fGtuParam->IsInZChannel(fStack, layer, 1, trk->GetZbin()) ? trk->GetSubChannel(1) : -1,
333 fGtuParam->IsInZChannel(fStack, layer, 0, trk->GetZbin()) ? trk->GetSubChannel(0) : -1
339 Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
341 // run the track finding
343 Int_t *notr = new Int_t[fGtuParam->GetNLayers()];
344 Int_t *ptrA = new Int_t[fGtuParam->GetNLayers()];
345 Int_t *ptrB = new Int_t[fGtuParam->GetNLayers()];
347 Bool_t *bHitA = new Bool_t[fGtuParam->GetNLayers()];
348 Bool_t *bHitB = new Bool_t[fGtuParam->GetNLayers()];
349 Bool_t *bAligned = new Bool_t[fGtuParam->GetNLayers()];
350 Bool_t *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
351 Bool_t *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
352 Bool_t *bDone = new Bool_t[fGtuParam->GetNLayers()];
355 Int_t *inc = new Int_t[fGtuParam->GetNLayers()];
356 Int_t *incprime = new Int_t[fGtuParam->GetNLayers()];
358 // ----- signals within current layer -----
369 AliTRDtrackletGTU *trkRA = 0x0; // reference tracklet A
370 AliTRDtrackletGTU *trkRB = 0x0; // reference tracklet B
371 AliTRDtrackletGTU *trkA = 0x0; // tracklet A in current layer
372 AliTRDtrackletGTU *trkB = 0x0; // tracklet B in current layer
374 AliTRDtrackletGTU *trkEnd = new AliTRDtrackletGTU();
375 for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
376 trkEnd->SetSubChannel(i, 7);
381 for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
382 Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
383 AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
385 ready = kFALSE; // ready if all channels done
387 // for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
388 // for (Int_t iTrkl = 0; iTrkl < fZChannelTracklets[iLayer][zch].GetEntries(); iTrkl++) {
389 // printf("%4i/%4i(%i,%i) ",
390 // ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetYProj(),
391 // ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetAlpha(),
392 // ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetIndex(),
393 // ((AliTRDtrackletGTU*) fZChannelTracklets[iLayer][zch].At(iTrkl))->GetSubChannel(zch)
399 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
400 notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
401 ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
402 ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
404 // not necessary here
405 // bDone[layer] = (ptrB[layer] >= notr[layer] - 1); // trkB is last one
406 // bDone[layer] = (notr[layer] <= 0);
407 // ready = ready && bDone[layer];
410 AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
413 if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
417 // ----- reference tracklets -----
420 if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
421 trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
423 AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
427 if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
428 trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
430 yPlus = trkRA->GetYProj() + fGtuParam->GetDeltaY();
431 yMinus = trkRA->GetYProj() - fGtuParam->GetDeltaY();
432 alphaPlus = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
433 alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
435 ybPlus = trkRB->GetYProj() + fGtuParam->GetDeltaY();
436 ybMinus = trkRB->GetYProj() - fGtuParam->GetDeltaY();
438 else { // irrelevant (should be, is it?)
439 ybPlus = trkRA->GetYProj() + fGtuParam->GetDeltaY();
440 ybMinus = trkRA->GetYProj() - fGtuParam->GetDeltaY();
447 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
448 bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
449 inc[layer] = incprime[layer] = 0;
451 if (layer == reflayer) {
452 bHitA[layer] = kTRUE;
453 bAligned[layer] = kTRUE;
460 if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
461 trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
462 if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
463 trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
465 bAlignedA[layer] = kFALSE;
466 bAlignedB[layer] = kFALSE;
469 bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
470 !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus) ) &&
471 !(trkA->GetAlpha() < alphaMinus) &&
472 !(trkA->GetAlpha() > alphaPlus) );
473 bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
477 bAlignedA[layer] = kTRUE;
481 bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
482 !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
483 !(trkB->GetAlpha() < alphaMinus) &&
484 !(trkB->GetAlpha() > alphaPlus) );
485 bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
489 bAlignedB[layer] = kTRUE;
492 bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
493 // bAligned[layer] = bAlignedA[layer]; //???
495 if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
497 else if (!bAligned[layer] )
501 if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
508 // pre-calculation for the layer shifting (alignment w. r. t. trkRB)
511 if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
521 if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
528 } // end of loop over layers
530 AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
531 nHits, (nHits >= 4) ? "-> track found" : ""));
534 // ----- track registration -----
535 AliTRDtrackGTU *track = new ((*fTrackArray)[fTrackArray->GetEntriesFast()]) AliTRDtrackGTU();
536 track->SetSector(fSector);
537 track->SetStack(fStack);
538 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
539 if (bHitA[layer] || layer == reflayer)
540 track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
541 else if (bHitB[layer])
542 track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
545 Bool_t registerTrack = kTRUE;
546 for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
547 if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
548 if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
549 AliDebug(1,"Not registered");
550 registerTrack = kFALSE;
555 track->SetZChannel(zch);
556 track->SetRefLayerIdx(refLayerIdx);
557 fTracks[zch][refLayerIdx].Add(track);
561 if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
563 else if (nWayBeyond > 2) // no track possible for both reference tracklets
568 if (inc[reflayer] != 0) // reflayer is shifted
569 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
570 if (layer == reflayer)
572 inc[layer] = incprime[layer];
574 else { // reflayer is not shifted
575 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
576 if (layer == reflayer || notr[layer] == 0)
581 if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
582 trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
584 if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
585 trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
588 if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
589 !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus ) ) ) // trkA could hit trkRA
590 // if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) )
593 if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
595 else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
598 inc[layer] = incprime[layer];
601 inc[layer] = incprime[layer];
607 for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
609 bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
610 ready = ready && bDone[layer];
612 if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
613 AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
615 // AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
616 AliDebug(10,Form(" Layer: %i %2i(%2i, %2i, %4i, %3i)%s%s %2i(%2i, %2i, %4i, %3i)%s%s +%i %s (no: %i)",
619 (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetIndex() : -1,
620 (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetSubChannel(zch) : -1,
621 (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetYProj() : -1,
622 (0 <= ptrA[layer] && ptrA[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]))->GetAlpha() : -1,
623 bHitA[layer] ? "*" : " ", bAlignedA[layer] ? "+" : " ",
625 (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetIndex() : -1,
626 (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetSubChannel(zch) : -1,
627 (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetYProj() : -1,
628 (0 <= ptrB[layer] && ptrB[layer] < notr[layer]) ? ((AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]))->GetAlpha() : -1,
629 bHitB[layer] ? "*" : " ", bAlignedB[layer] ? "+" : " ",
630 inc[layer], bDone[layer] ? "done" : " ", notr[layer]));
631 ptrA[layer] += inc[layer];
632 ptrB[layer] += inc[layer];
636 } // end of loop over reflayer
653 Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
655 TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
656 TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
658 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
659 tracksRefMerged[zch] = new TList;
660 tracksRefUnique[zch] = new TList;
663 TList *tracksZMergedStage0 = new TList;
664 TList *tracksZUniqueStage0 = new TList;
666 TList **tracksZSplitted = new TList*[2];
667 for (Int_t i = 0; i < 2; i++)
668 tracksZSplitted[i] = new TList;
670 TList *tracksZMergedStage1 = new TList;
672 AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
674 // Bool_t done = kFALSE;
676 AliTRDtrackGTU *trkStage0 = 0x0;
678 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
679 // ----- Merging and Unification in Reflayers (seed_merger) -----
683 for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
684 trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
685 if (trkInRefLayer[refLayerIdx] == 0) {
688 else if (trkStage0 == 0x0 ) {
689 trkStage0 = trkInRefLayer[refLayerIdx];
690 minIdx = refLayerIdx;
693 else if ( (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel()) ||
694 ((trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel()) &&
695 ((trkInRefLayer[refLayerIdx]->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
696 minIdx = refLayerIdx;
697 trkStage0 = trkInRefLayer[refLayerIdx];
703 tracksRefMerged[zch]->Add(trkStage0);
704 fTracks[zch][minIdx].RemoveFirst();
705 } while (trkStage0 != 0);
707 Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
709 AliDebug(2, Form("zchannel %i", zch));
710 TIter trackRefMerged(tracksRefMerged[zch]);
711 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
712 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
713 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
714 trk->GetTrackletIndex(5),
715 trk->GetTrackletIndex(4),
716 trk->GetTrackletIndex(3),
717 trk->GetTrackletIndex(2),
718 trk->GetTrackletIndex(1),
719 trk->GetTrackletIndex(0),
720 trk->GetYapprox() >> 3,
721 trk->GetZSubChannel()));
722 AliDebug(2, "uniquified:");
723 TIter trackRefUnique(tracksRefUnique[zch]);
724 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
725 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, z_idx=%i",
726 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
727 trk->GetTrackletIndex(5),
728 trk->GetTrackletIndex(4),
729 trk->GetTrackletIndex(3),
730 trk->GetTrackletIndex(2),
731 trk->GetTrackletIndex(1),
732 trk->GetTrackletIndex(0),
733 trk->GetYapprox() >> 3,
734 trk->GetZSubChannel()));
737 // ----- Merging in zchannels - 1st stage -----
739 if (AliTRDgtuParam::GetUseGTUmerge()) {
742 Bool_t lowerThan[3] = { kFALSE, kFALSE, kFALSE };
743 AliTRDtrackGTU *trk[3] = { 0x0, 0x0, 0x0 };
744 for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel)
745 trk[iChannel] = (AliTRDtrackGTU*) tracksRefUnique[iChannel]->First();
747 for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel) {
748 AliTRDtrackGTU *trk1 = trk[iChannel];
749 AliTRDtrackGTU *trk2 = trk[(iChannel + 1) % 3];
751 Int_t sortnum1 = (trk1->GetZChannel() + 3 * trk1->GetZSubChannel()) / 2 - 1;
752 Int_t sortnum2 = (trk2->GetZChannel() + 3 * trk2->GetZSubChannel()) / 2 - 1;
753 AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
754 trk1->GetZChannel(), trk2->GetZChannel(),
755 sortnum1, sortnum2));
756 if ( (sortnum1 < sortnum2) ||
757 ((sortnum1 == sortnum2) &&
758 ((trk1->GetYapprox() >> 3) < (trk2->GetYapprox() >> 3)) ) ) {
759 lowerThan[iChannel] = kTRUE;
764 notEmpty = (trk[2] ? (1 << 2) : 0) |
765 (trk[1] ? (1 << 1) : 0) |
766 (trk[0] ? (1 << 0) : 0);
823 tracksZMergedStage0->Add(trk[pop]);
824 tracksRefUnique[pop]->RemoveFirst();
829 // there is still a difference between this implementation and
830 // the hardware algorithm, only for expert use
835 // compare tracks from all adjacent zchannels
836 // (including comparison of channels 2 and 0)
837 for (Int_t i = fGtuParam->GetNZChannels() - 1; i > -1; i--) {
839 AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
843 else if (trkStage0 == 0x0 ) {
849 Int_t sortnum1 = (trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1;
850 Int_t sortnum2 = (trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 - 1;
851 AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
852 trk->GetZChannel(), trkStage0->GetZChannel(),
853 sortnum1, sortnum2));
854 if ( (sortnum1 < sortnum2) ||
855 ((sortnum1 == sortnum2) &&
856 ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
866 tracksZMergedStage0->Add(trkStage0);
867 tracksRefUnique[minIdx]->RemoveFirst();
868 } while (trkStage0 != 0);
871 Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
873 AliDebug(2, "stage 0:");
874 TIter trackZMergedStage0(tracksZMergedStage0);
875 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
876 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
877 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
878 trk->GetTrackletIndex(5),
879 trk->GetTrackletIndex(4),
880 trk->GetTrackletIndex(3),
881 trk->GetTrackletIndex(2),
882 trk->GetTrackletIndex(1),
883 trk->GetTrackletIndex(0),
884 trk->GetYapprox() >> 3,
886 trk->GetZSubChannel()));
888 AliDebug(2, "uniquified:");
889 TIter trackZUniqueStage0(tracksZUniqueStage0);
890 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
891 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
892 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
893 trk->GetTrackletIndex(5),
894 trk->GetTrackletIndex(4),
895 trk->GetTrackletIndex(3),
896 trk->GetTrackletIndex(2),
897 trk->GetTrackletIndex(1),
898 trk->GetTrackletIndex(0),
899 trk->GetYapprox() >> 3,
901 trk->GetZSubChannel()));
903 // ----- Splitting in z -----
905 TIter next(tracksZUniqueStage0);
906 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
907 if ((trk->GetZChannel() < 0) ||
908 (trk->GetZChannel() > 2) ||
909 (trk->GetZSubChannel() < 0) ||
910 (trk->GetZSubChannel() > 6)) {
911 AliError(Form("track with invalid z-channel information at %p: zch = %i, subchannel = %i",
912 trk, trk->GetZChannel(), trk->GetZSubChannel()));
915 Int_t idx = (trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2;
916 if ((idx < 0) || (idx > 1)) {
917 AliError(Form("invalid index %i null", idx));
921 if (!tracksZSplitted[idx]) {
922 AliError(Form("array pointer %i null", idx));
925 tracksZSplitted[idx]->Add(trk);
928 for (Int_t i = 0; i < 2; i++) {
929 AliDebug(2, Form("split %i:", i));
930 TIter trackZSplit(tracksZSplitted[i]);
931 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
932 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
933 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
934 trk->GetTrackletIndex(5),
935 trk->GetTrackletIndex(4),
936 trk->GetTrackletIndex(3),
937 trk->GetTrackletIndex(2),
938 trk->GetTrackletIndex(1),
939 trk->GetTrackletIndex(0),
940 trk->GetYapprox() >> 3,
942 trk->GetZSubChannel()));
945 // ----- Merging in zchanels - 2nd stage -----
950 for (Int_t i = 1; i >= 0; i--) {
951 AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
955 else if (trkStage0 == 0x0 ) {
960 else if ( (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) < ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
961 ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) &&
962 ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
971 tracksZMergedStage1->Add(trkStage0);
972 tracksZSplitted[minIdx]->RemoveFirst();
973 } while (trkStage0 != 0);
975 Uniquifier(tracksZMergedStage1, ListOfTracks);
977 AliDebug(2, "merged:");
978 TIter trackZMergedStage1(tracksZMergedStage1);
979 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
980 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
981 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
982 trk->GetTrackletIndex(5),
983 trk->GetTrackletIndex(4),
984 trk->GetTrackletIndex(3),
985 trk->GetTrackletIndex(2),
986 trk->GetTrackletIndex(1),
987 trk->GetTrackletIndex(0),
988 trk->GetYapprox() >> 3,
990 trk->GetZSubChannel()));
992 AliDebug(2, "uniquified:");
993 TIter track(ListOfTracks);
994 while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
995 AliDebug(2, Form("track ref layer %i : %i %i %i %i %i %i, y=%i, zch=%i idx=%i",
996 AliTRDgtuParam::GetRefLayer(trk->GetRefLayerIdx()),
997 trk->GetTrackletIndex(5),
998 trk->GetTrackletIndex(4),
999 trk->GetTrackletIndex(3),
1000 trk->GetTrackletIndex(2),
1001 trk->GetTrackletIndex(1),
1002 trk->GetTrackletIndex(0),
1003 trk->GetYapprox() >> 3,
1005 trk->GetZSubChannel()));
1008 for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
1009 delete tracksRefMerged[zch];
1010 delete tracksRefUnique[zch];
1012 delete [] tracksRefMerged;
1013 delete [] tracksRefUnique;
1016 delete tracksZMergedStage0;
1017 delete tracksZUniqueStage0;
1019 for (Int_t i = 0; i < 2; i++)
1020 delete tracksZSplitted[i];
1021 delete [] tracksZSplitted;
1023 delete tracksZMergedStage1;
1025 delete [] trkInRefLayer;
1030 Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
1032 // run the track reconstruction for all tracks in the list
1034 TIter next(ListOfTracks);
1035 while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
1036 CalculateTrackParams(track);
1037 CalculatePID(track);
1038 AliDebug(1, Form("track with pid = %i", track->GetPID()));
1043 Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
1045 // calculate PID for the given track
1047 AliError("No track to calculate!");
1051 if (AliTRDgtuParam::GetUseGTUconst()) {
1052 // averaging as in GTU
1053 AliDebug(1, "using GTU constants for PID calculation");
1056 // selection of coefficient for averaging
1057 switch(track->GetTrackletMask()){
1060 coeff=0x5558; // approx. 1/6
1070 coeff=0x6666; // approx. 1/5
1075 coeff=0x8000; // = 0.25
1077 coeff &= 0x1ffff; // 17-bit constant
1081 for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
1082 if ((track->GetTrackletMask() >> iLayer) & 1) {
1083 sum += track->GetTracklet(iLayer)->GetPID();
1088 Float_t av = 1./i * sum;
1090 ULong64_t prod = (sum * coeff) & 0xfffffffffull;
1091 ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
1093 if (TMath::Abs((prodFinal & 0xff) - av) > 0.5)
1094 AliError(Form("problem with PID averaging (hw <-> ar): %3lli <-> %4.1f", prodFinal & 0xff, av));
1095 track->SetPID(prodFinal & 0xff);
1102 Int_t nTracklets = 0;
1104 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1105 if (!track->IsTrackletInLayer(layer)) {
1108 AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1109 pidSum += trk->GetPID();
1114 track->SetPID(pidSum/nTracklets);
1116 AliError("Track without contributing tracklets, no PID assigned");
1122 Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
1124 // calculate the track parameters
1127 AliError("No track to calculate!");
1137 AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
1139 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1140 if (!track->IsTrackletInLayer(layer)) {
1143 AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1145 AliError(Form("Could not get tracklet in layer %i\n", layer));
1148 AliDebug(10,Form(" layer %i trk yprime: %6i, aki: %6i", layer, trk->GetYPrime(),
1149 fGtuParam->GetAki(track->GetTrackletMask(), layer)));
1150 a += (((fGtuParam->GetAki(track->GetTrackletMask(), layer) * trk->GetYPrime()) >> 7) + 1) >> 1;
1151 b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1152 c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1158 track->SetFitParams(a << 2, b, c);
1160 fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1162 AliDebug(5,Form(" Track parameters: a16 = %i, a18 = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f (trkl mask: %i)",
1163 a, a << 2, b, c, x1, x2, track->GetPt(), track->GetTrackletMask()));
1169 Bool_t AliTRDgtuTMU::Uniquifier(const TList *inlist, TList *outlist)
1171 // remove multiple occurences of the same track
1174 AliTRDtrackGTU *trkStage0 = 0x0;
1175 AliTRDtrackGTU *trkStage1 = 0x0;
1178 trkStage0 = (AliTRDtrackGTU*) next();
1180 Bool_t tracksEqual = kFALSE;
1181 if (trkStage0 != 0 && trkStage1 != 0) {
1182 for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1183 if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1184 tracksEqual = kTRUE;
1191 if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
1192 trkStage1 = trkStage0;
1195 if (trkStage1 != 0x0)
1196 outlist->Add(trkStage1);
1197 trkStage1 = trkStage0;
1200 } while (trkStage1 != 0x0);