Line data Source code
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 48 : ClassImp(AliTRDgtuTMU)
43 :
44 : AliTRDgtuTMU::AliTRDgtuTMU(Int_t stack, Int_t sector) :
45 4 : TObject(),
46 4 : fTracklets(0x0),
47 4 : fTrackletsPostInput(0x0),
48 4 : fZChannelTracklets(0x0),
49 12 : fTrackArray(new TClonesArray("AliTRDtrackGTU", 50)),
50 4 : fTracks(0x0),
51 4 : fGtuParam(0x0),
52 4 : fStack(-1),
53 4 : fSector(-1)
54 20 : {
55 : // constructor which initializes the position information of the TMU
56 :
57 8 : fGtuParam = AliTRDgtuParam::Instance();
58 :
59 : // store tracklets per link
60 8 : fTracklets = new TObjArray*[fGtuParam->GetNLinks()];
61 104 : for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
62 144 : fTracklets[iLink] = new TObjArray();
63 : }
64 :
65 : // tracklets after input units per layer
66 8 : fTrackletsPostInput = new TObjArray*[fGtuParam->GetNLayers()];
67 8 : fZChannelTracklets = new TList*[fGtuParam->GetNLayers()];
68 56 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
69 240 : fZChannelTracklets[layer] = new TList[fGtuParam->GetNZChannels()];
70 72 : fTrackletsPostInput[layer] = new TObjArray();
71 : }
72 :
73 8 : fTracks = new TList*[fGtuParam->GetNZChannels()];
74 32 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
75 120 : fTracks[zch] = new TList[fGtuParam->GetNRefLayers()];
76 : }
77 :
78 4 : if (stack > -1)
79 0 : SetStack(stack);
80 4 : if (sector > -1)
81 0 : SetSector(sector);
82 8 : }
83 :
84 : AliTRDgtuTMU::~AliTRDgtuTMU()
85 24 : {
86 : // destructor
87 :
88 32 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
89 72 : delete [] fTracks[zch];
90 : }
91 8 : delete [] fTracks;
92 8 : delete fTrackArray;
93 :
94 56 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
95 144 : delete [] fZChannelTracklets[layer];
96 48 : delete fTrackletsPostInput[layer];
97 : }
98 8 : delete [] fZChannelTracklets;
99 8 : delete [] fTrackletsPostInput;
100 :
101 104 : for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
102 96 : delete fTracklets[iLink];
103 : }
104 8 : delete [] fTracklets;
105 12 : }
106 :
107 : Bool_t AliTRDgtuTMU::SetSector(Int_t sector)
108 : {
109 : // set the sector
110 :
111 117 : if (sector > -1 && sector < fGtuParam->GetGeo()->Nsector() ) {
112 39 : fSector = sector;
113 39 : return kTRUE;
114 : }
115 :
116 0 : AliError(Form("Invalid sector given: %i", sector));
117 0 : return kFALSE;
118 39 : }
119 :
120 : Bool_t AliTRDgtuTMU::SetStack(Int_t stack)
121 : {
122 : // Set the stack (necessary for tracking)
123 :
124 117 : if (stack > -1 && stack < fGtuParam->GetGeo()->Nstack() ) {
125 39 : fStack = stack;
126 39 : return kTRUE;
127 : }
128 :
129 0 : AliError(Form("Invalid stack given: %i", stack));
130 0 : return kFALSE;
131 39 : }
132 :
133 : Bool_t AliTRDgtuTMU::Reset()
134 : {
135 : // delete all tracks
136 :
137 387 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
138 1032 : for (Int_t reflayeridx = 0; reflayeridx < fGtuParam->GetNRefLayers(); reflayeridx++) {
139 387 : fTracks[zch][reflayeridx].Clear();
140 : }
141 : }
142 :
143 43 : fTrackArray->Delete();
144 :
145 : // delete all added tracklets
146 1118 : for (Int_t iLink = 0; iLink < fGtuParam->GetNLinks(); iLink++) {
147 516 : fTracklets[iLink]->Clear();
148 : }
149 602 : for (Int_t layer = 0; layer < fGtuParam->GetNLinks()/2; layer++) {
150 258 : fTrackletsPostInput[layer]->Clear();
151 2064 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++)
152 774 : fZChannelTracklets[layer][zch].Clear();
153 : }
154 :
155 43 : fSector = -1;
156 43 : fStack = -1;
157 :
158 43 : return kTRUE;
159 : }
160 :
161 : Bool_t AliTRDgtuTMU::AddTracklet(AliTRDtrackletGTU *tracklet, Int_t link)
162 : {
163 : // add a tracklet on the given link
164 :
165 720 : fTracklets[link]->Add(tracklet);
166 360 : return kTRUE;
167 : }
168 :
169 :
170 : Bool_t AliTRDgtuTMU::RunTMU(TList *ListOfTracks, AliESDEvent *esd, Int_t outLabel)
171 : {
172 : // performs the analysis as in a TMU module of the GTU, i. e.
173 : // track matching
174 : // calculation of track parameteres (pt, deflection, ???)
175 :
176 117 : if (fStack < 0 || fSector < 0) {
177 0 : AliError("No valid stack/sector set for this TMU! No tracking!");
178 0 : return kFALSE;
179 : }
180 :
181 : // ----- Input units -----
182 117 : AliDebug(1,"--------- Running Input units ----------");
183 585 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
184 234 : if (!RunInputUnit(layer)) {
185 0 : AliError(Form("Input unit in layer %i failed", layer));
186 0 : return kFALSE;
187 : }
188 : }
189 :
190 : // ----- Z-channel units -----
191 117 : AliDebug(1,"--------- Running Z-channel units ----------");
192 585 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
193 234 : if (!RunZChannelUnit(layer)) {
194 0 : AliError(Form("Z-Channel unit in layer %i failed", layer));
195 0 : return kFALSE;
196 : }
197 : }
198 :
199 : // ----- track finding -----
200 117 : AliDebug(1,"--------- Running tracking units ----------");
201 351 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
202 117 : if (!RunTrackFinder(zch, ListOfTracks)) {
203 0 : AliError(Form("Track Finder in z-channel %i failed", zch));
204 0 : return kFALSE;
205 : }
206 : }
207 :
208 : // ----- Track Merging -----
209 39 : if (!RunTrackMerging(ListOfTracks)) {
210 0 : AliError("Track merging failed");
211 0 : return kFALSE;
212 : }
213 :
214 : // ----- track reconstruction -----
215 39 : if (!RunTrackReconstruction(ListOfTracks)) {
216 0 : AliError("Track reconstruction failed");
217 0 : return kFALSE;
218 : }
219 :
220 : // ----- label calculation and ESD storage -----
221 39 : TIter next(ListOfTracks);
222 186 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
223 15 : if (outLabel == -1)
224 15 : trk->CookLabel();
225 : else
226 0 : trk->SetLabel(outLabel);
227 15 : if (esd) {
228 0 : AliESDTrdTrack *trdtrack = trk->CreateTrdTrack();
229 0 : esd->AddTrdTrack(trdtrack);
230 0 : delete trdtrack;
231 0 : }
232 15 : }
233 :
234 : return kTRUE;
235 78 : }
236 :
237 : Bool_t AliTRDgtuTMU::RunInputUnit(Int_t layer)
238 : {
239 : // precalculations for the tracking and reconstruction
240 :
241 : Int_t iTrkl0 = 0; // A-side tracklet
242 : Int_t iTrkl1 = 0; // B-side tracklet
243 :
244 : Int_t nTracklets = 0;
245 1808 : while ((!fGtuParam->GetLimitNoTracklets() ||
246 594 : (nTracklets < fGtuParam->GetMaxNoTracklets())) &&
247 594 : ((iTrkl0 < fTracklets[2*layer + 0]->GetEntriesFast()) ||
248 386 : (iTrkl1 < fTracklets[2*layer + 1]->GetEntriesFast()))) {
249 360 : if (iTrkl1 >= fTracklets[2*layer + 1]->GetEntriesFast()) {
250 134 : fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
251 134 : iTrkl0++;
252 134 : }
253 : else {
254 226 : if (iTrkl0 >= fTracklets[2*layer + 0]->GetEntriesFast()) {
255 152 : fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
256 152 : iTrkl1++;
257 152 : }
258 : else {
259 222 : if (((AliTRDtrackletGTU*) fTracklets[2*layer + 1]->At(iTrkl1))->GetZbin() <
260 74 : ((AliTRDtrackletGTU*) fTracklets[2*layer + 0]->At(iTrkl0))->GetZbin()) {
261 125 : fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 1]->At(iTrkl1));
262 51 : iTrkl1++;
263 :
264 51 : }
265 : else {
266 23 : fTrackletsPostInput[layer]->AddLast(fTracklets[2*layer + 0]->At(iTrkl0));
267 23 : iTrkl0++;
268 : }
269 : }
270 : }
271 360 : ++nTracklets;
272 : }
273 :
274 234 : TIter next(fTrackletsPostInput[layer]);
275 :
276 1656 : while ( AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next() ) {
277 720 : trk->SetIndex(fTrackletsPostInput[layer]->IndexOf(trk));
278 :
279 1080 : Int_t alpha = (trk->GetYbin() >> fGtuParam->GetBitExcessY()) * fGtuParam->GetCiAlpha(layer);
280 720 : alpha = ( 2 * trk->GetdY() - (alpha >> fGtuParam->GetBitExcessAlpha()) + 1 ) >> 1;
281 : // wrapping projected alpha as in hardware
282 360 : if ((alpha < -64) || (alpha > 63))
283 185 : AliDebug(1, Form("alpha out of range: %i", alpha));
284 360 : alpha = alpha & 0x7f;
285 360 : if (alpha & 0x40)
286 193 : trk->SetAlpha(0xffffffc0 | alpha);
287 : else
288 167 : trk->SetAlpha(alpha);
289 :
290 1080 : Int_t yproj = trk->GetdY() * (fGtuParam->GetCiYProj(layer));
291 720 : yproj = ( ( ( (yproj >> fGtuParam->GetBitExcessYProj()) + trk->GetYbin() ) >> 2) + 1) >> 1;
292 360 : trk->SetYProj(yproj);
293 :
294 1440 : trk->SetYPrime(trk->GetYbin() + fGtuParam->GetYt(fStack, layer, trk->GetZbin()));
295 :
296 1800 : 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()) ));
300 360 : }
301 : return kTRUE;
302 234 : }
303 :
304 : Bool_t AliTRDgtuTMU::RunZChannelUnit(Int_t layer)
305 : {
306 : // run the z-channel unit
307 :
308 468 : TIter next(fTrackletsPostInput[layer]);
309 :
310 1656 : while (AliTRDtrackletGTU *trk = (AliTRDtrackletGTU*) next()) {
311 2880 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
312 3240 : if (fGtuParam->IsInZChannel(fStack, layer, zch, trk->GetZbin()) ) {
313 2313 : trk->SetSubChannel(zch, fGtuParam->GetZSubchannel(fStack, layer, zch, trk->GetZbin()) );
314 :
315 771 : TIter nexttrkl(&fZChannelTracklets[layer][zch], kIterBackward);
316 : AliTRDtrackletGTU *t = 0x0;
317 3259 : while ((t = (AliTRDtrackletGTU*) nexttrkl.Next())) {
318 3372 : if (t->GetSubChannel(zch) < trk->GetSubChannel(zch) ||
319 3168 : ((t->GetSubChannel(zch) == trk->GetSubChannel(zch)) && (t->GetYProj() < trk->GetYProj())) ) {
320 : break;
321 : }
322 : }
323 771 : if (t)
324 1158 : fZChannelTracklets[layer][zch].AddAfter(t, trk);
325 : else
326 384 : fZChannelTracklets[layer][zch].AddFirst(trk);
327 771 : }
328 : }
329 1800 : 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
334 : ));
335 360 : }
336 : return kTRUE;
337 234 : }
338 :
339 : Bool_t AliTRDgtuTMU::RunTrackFinder(Int_t zch, TList* /* ListOfTracks */)
340 : {
341 : // run the track finding
342 :
343 234 : Int_t *notr = new Int_t[fGtuParam->GetNLayers()];
344 117 : Int_t *ptrA = new Int_t[fGtuParam->GetNLayers()];
345 117 : Int_t *ptrB = new Int_t[fGtuParam->GetNLayers()];
346 :
347 117 : Bool_t *bHitA = new Bool_t[fGtuParam->GetNLayers()];
348 117 : Bool_t *bHitB = new Bool_t[fGtuParam->GetNLayers()];
349 117 : Bool_t *bAligned = new Bool_t[fGtuParam->GetNLayers()];
350 117 : Bool_t *bAlignedA = new Bool_t[fGtuParam->GetNLayers()];
351 117 : Bool_t *bAlignedB = new Bool_t[fGtuParam->GetNLayers()];
352 117 : Bool_t *bDone = new Bool_t[fGtuParam->GetNLayers()];
353 : Bool_t ready;
354 :
355 117 : Int_t *inc = new Int_t[fGtuParam->GetNLayers()];
356 117 : Int_t *incprime = new Int_t[fGtuParam->GetNLayers()];
357 :
358 : // ----- signals within current layer -----
359 : Int_t yPlus;
360 : Int_t yMinus;
361 : Int_t ybPlus;
362 : Int_t ybMinus;
363 : Int_t alphaPlus;
364 : Int_t alphaMinus;
365 : Int_t nHits;
366 : Int_t nUnc;
367 : Int_t nWayBeyond;
368 :
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
373 : /*
374 : AliTRDtrackletGTU *trkEnd = new AliTRDtrackletGTU();
375 : for (Int_t i = 0; i < fGtuParam->GetNZChannels(); i++)
376 : trkEnd->SetSubChannel(i, 7);
377 : trkEnd->SetYProj(0);
378 : trkEnd->SetAlpha(0);
379 : */
380 :
381 936 : for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
382 351 : Int_t reflayer = fGtuParam->GetRefLayer(refLayerIdx);
383 1053 : AliDebug(5,Form("Tracking for z-channel: %i, reflayer: %i", zch, reflayer));
384 :
385 : ready = kFALSE; // ready if all channels done
386 :
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)
394 : // );
395 : // }
396 : // printf("\n");
397 : // }
398 :
399 4914 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
400 2106 : notr[layer] = fZChannelTracklets[layer][zch].GetEntries();
401 2106 : ptrA[layer] = 0; // notr[layer] > 0 ? 0 : -1;
402 2106 : ptrB[layer] = 1; // notr[layer] > 1 ? 1 : -1;
403 :
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];
408 :
409 2106 : if (reflayer == 1)
410 2106 : AliDebug(5,Form("in layer: %i (zchannel = %i) there are: %i tracklets", layer, zch, notr[layer]));
411 : }
412 :
413 351 : if (ptrA[reflayer] < 0 && ptrB[reflayer] < 0)
414 0 : continue;
415 :
416 1005 : while (!ready) {
417 : // ----- reference tracklets -----
418 : trkRA = 0x0;
419 : trkRB = 0x0;
420 1224 : if (0 <= ptrA[reflayer] && ptrA[reflayer] < notr[reflayer])
421 327 : trkRA = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrA[reflayer]);
422 : else {
423 855 : AliDebug(10,Form("No valid tracklet in the reference at ptr: %i! Nothing done!", ptrA[reflayer]));
424 : break;
425 : }
426 :
427 654 : if (0 <= ptrB[reflayer] && ptrB[reflayer] < notr[reflayer])
428 196 : trkRB = (AliTRDtrackletGTU*) fZChannelTracklets[reflayer][zch].At(ptrB[reflayer]);
429 :
430 327 : yPlus = trkRA->GetYProj() + fGtuParam->GetDeltaY();
431 327 : yMinus = trkRA->GetYProj() - fGtuParam->GetDeltaY();
432 327 : alphaPlus = trkRA->GetAlpha() + fGtuParam->GetDeltaAlpha();
433 327 : alphaMinus = trkRA->GetAlpha() - fGtuParam->GetDeltaAlpha();
434 327 : if (trkRB) {
435 196 : ybPlus = trkRB->GetYProj() + fGtuParam->GetDeltaY();
436 196 : ybMinus = trkRB->GetYProj() - fGtuParam->GetDeltaY();
437 196 : }
438 : else { // irrelevant (should be, is it?)
439 131 : ybPlus = trkRA->GetYProj() + fGtuParam->GetDeltaY();
440 131 : ybMinus = trkRA->GetYProj() - fGtuParam->GetDeltaY();
441 : }
442 :
443 : nHits = 0;
444 : nUnc = 0;
445 : nWayBeyond = 0;
446 :
447 4578 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
448 1962 : bHitA[layer] = bHitB[layer] = bAligned[layer] = kFALSE;
449 1962 : inc[layer] = incprime[layer] = 0;
450 :
451 1962 : if (layer == reflayer) {
452 327 : bHitA[layer] = kTRUE;
453 327 : bAligned[layer] = kTRUE;
454 327 : nHits++;
455 327 : continue;
456 : }
457 :
458 : trkA = 0x0;
459 : trkB = 0x0;
460 3270 : if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
461 1029 : trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
462 3270 : if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
463 600 : trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
464 :
465 1635 : bAlignedA[layer] = kFALSE;
466 1635 : bAlignedB[layer] = kFALSE;
467 :
468 1635 : if (trkA) {
469 3810 : bHitA[layer] = ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
470 2112 : !(trkA->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() > yPlus) ) &&
471 465 : !(trkA->GetAlpha() < alphaMinus) &&
472 411 : !(trkA->GetAlpha() > alphaPlus) );
473 3810 : bAlignedA[layer] = !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) );
474 1029 : }
475 : else {
476 606 : bHitA[layer] = 0;
477 606 : bAlignedA[layer] = kTRUE;
478 : }
479 :
480 1635 : if (trkB) {
481 2149 : bHitB[layer] = ( !(trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) ) &&
482 1138 : !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) &&
483 120 : !(trkB->GetAlpha() < alphaMinus) &&
484 102 : !(trkB->GetAlpha() > alphaPlus) );
485 2147 : bAlignedB[layer] = (trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) );
486 600 : }
487 : else {
488 1035 : bHitB[layer] = 0;
489 1035 : bAlignedB[layer] = kTRUE;
490 : }
491 :
492 3547 : bAligned[layer] = bAlignedA[layer] || bAlignedB[layer]; //???
493 : // bAligned[layer] = bAlignedA[layer]; //???
494 :
495 4190 : if (bAligned[layer] && (bHitA[layer] || bHitB[layer]) )
496 354 : nHits++;
497 1281 : else if (!bAligned[layer] )
498 186 : nUnc++;
499 1635 : if (trkRB) {
500 980 : if (trkA) {
501 1556 : if ((trkA->GetSubChannel(zch) > trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() > ybPlus) )
502 161 : nWayBeyond++;
503 : }
504 : else
505 388 : nWayBeyond++;
506 : }
507 :
508 : // pre-calculation for the layer shifting (alignment w. r. t. trkRB)
509 1635 : if (trkA) {
510 1029 : if(trkRB) {
511 1483 : if ((trkA->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkA->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkA->GetYProj() < ybMinus )) // could trkA be aligned for trkRB
512 310 : incprime[layer] = 1;
513 : else
514 282 : incprime[layer] = 0;
515 : }
516 : else
517 437 : incprime[layer] = 1;
518 :
519 1029 : if (trkB) {
520 600 : if (trkRB) {
521 1170 : if ((trkB->GetSubChannel(zch) < trkRB->GetSubChannel(zch)) || (trkB->GetSubChannel(zch) == trkRB->GetSubChannel(zch) && trkB->GetYProj() < ybMinus )) // could trkB be aligned for trkRB
522 200 : incprime[layer] = 2;
523 : }
524 : else
525 134 : incprime[layer] = 2;
526 : }
527 : }
528 : } // end of loop over layers
529 :
530 981 : AliDebug(5,Form("logic calculation finished, Nhits: %i %s",
531 : nHits, (nHits >= 4) ? "-> track found" : ""));
532 :
533 327 : if (nHits >= 4) {
534 : // ----- track registration -----
535 71 : AliTRDtrackGTU *track = new ((*fTrackArray)[fTrackArray->GetEntriesFast()]) AliTRDtrackGTU();
536 71 : track->SetSector(fSector);
537 71 : track->SetStack(fStack);
538 994 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++ ) {
539 525 : if (bHitA[layer] || layer == reflayer)
540 327 : track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrA[layer]), layer );
541 99 : else if (bHitB[layer])
542 7 : track->AddTracklet((AliTRDtrackletGTU* ) fZChannelTracklets[layer][zch].At(ptrB[layer]), layer );
543 : }
544 :
545 : Bool_t registerTrack = kTRUE;
546 288 : for (Int_t layerIdx = refLayerIdx-1; layerIdx >= 0; layerIdx--) {
547 73 : if (track->IsTrackletInLayer(fGtuParam->GetRefLayer(layerIdx))) {
548 50 : if ((track->GetTracklet(fGtuParam->GetRefLayer(layerIdx)))->GetSubChannel(zch) > 0) {
549 150 : AliDebug(1,"Not registered");
550 : registerTrack = kFALSE;
551 50 : }
552 : }
553 : }
554 71 : if (registerTrack) {
555 32 : track->SetZChannel(zch);
556 32 : track->SetRefLayerIdx(refLayerIdx);
557 32 : fTracks[zch][refLayerIdx].Add(track);
558 32 : }
559 71 : }
560 :
561 431 : if ( (nUnc != 0) && (nUnc + nHits >= 4) ) // could this position of the reference layer give some track //??? special check in case of hit?
562 40 : inc[reflayer] = 0;
563 574 : else if (nWayBeyond > 2) // no track possible for both reference tracklets
564 394 : inc[reflayer] = 2;
565 : else
566 180 : inc[reflayer] = 1;
567 :
568 327 : if (inc[reflayer] != 0) // reflayer is shifted
569 4018 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
570 1722 : if (layer == reflayer)
571 : continue;
572 1435 : inc[layer] = incprime[layer];
573 1722 : }
574 : else { // reflayer is not shifted
575 560 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
576 440 : if (layer == reflayer || notr[layer] == 0)
577 : continue;
578 161 : inc[layer] = 0;
579 : trkA = 0x0;
580 : trkB = 0x0;
581 322 : if (0 <= ptrA[layer] && ptrA[layer] < notr[layer])
582 154 : trkA = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrA[layer]);
583 :
584 322 : if (0 <= ptrB[layer] && ptrB[layer] < notr[layer])
585 145 : trkB = (AliTRDtrackletGTU*) fZChannelTracklets[layer][zch].At(ptrB[layer]);
586 :
587 161 : if (trkA) {
588 359 : if ( !(trkA->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkA->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkA->GetYProj() < yMinus) ) &&
589 139 : !(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) ) )
591 36 : inc[layer] = 0;
592 118 : else if (trkB) {
593 194 : if ( trkB->GetSubChannel(zch) < trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() < yMinus) )
594 94 : inc[layer] = 2;
595 46 : else if ( !(trkB->GetSubChannel(zch) > trkRA->GetSubChannel(zch) || (trkB->GetSubChannel(zch) == trkRA->GetSubChannel(zch) && trkB->GetYProj() > yPlus) ) )
596 8 : inc[layer] = 1;
597 : else
598 10 : inc[layer] = incprime[layer];
599 : }
600 : else
601 6 : inc[layer] = incprime[layer];
602 : }
603 : }
604 : }
605 :
606 : ready = kTRUE;
607 4578 : for (Int_t layer = fGtuParam->GetNLayers()-1; layer >= 0; layer--) {
608 :
609 5886 : bDone[layer] = ptrB[layer] < 0 || ptrB[layer] >= notr[layer];
610 4911 : ready = ready && bDone[layer];
611 :
612 2970 : if (inc[layer] != 0 && ptrA[layer] >= notr[layer])
613 0 : AliError(Form("Invalid increment: %i at ptrA: %i, notr: %i", inc[layer], ptrA[layer], notr[layer]));
614 :
615 : // AliInfo(Form("Shifting layer: %i, notr: %i, ptrA: %i, ptrB: %i, inc: %i", layer, notr[layer], ptrA[layer], ptrB[layer], inc[layer]));
616 5886 : AliDebug(10,Form(" Layer: %i %2i(%2i, %2i, %4i, %3i)%s%s %2i(%2i, %2i, %4i, %3i)%s%s +%i %s (no: %i)",
617 : layer,
618 : ptrA[layer],
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] ? "+" : " ",
624 : ptrB[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 1962 : ptrA[layer] += inc[layer];
632 1962 : ptrB[layer] += inc[layer];
633 : }
634 :
635 : } // end of while
636 351 : } // end of loop over reflayer
637 :
638 234 : delete [] bHitA;
639 234 : delete [] bHitB;
640 234 : delete [] bAligned;
641 234 : delete [] bDone;
642 234 : delete [] inc;
643 234 : delete [] incprime;
644 234 : delete [] bAlignedA;
645 234 : delete [] bAlignedB;
646 234 : delete [] notr;
647 234 : delete [] ptrA;
648 234 : delete [] ptrB;
649 :
650 117 : return kTRUE;
651 0 : }
652 :
653 : Bool_t AliTRDgtuTMU::RunTrackMerging(TList* ListOfTracks)
654 : {
655 78 : TList **tracksRefMerged = new TList*[fGtuParam->GetNZChannels()];
656 39 : TList **tracksRefUnique = new TList*[fGtuParam->GetNZChannels()];
657 :
658 312 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
659 234 : tracksRefMerged[zch] = new TList;
660 234 : tracksRefUnique[zch] = new TList;
661 : }
662 :
663 39 : TList *tracksZMergedStage0 = new TList;
664 39 : TList *tracksZUniqueStage0 = new TList;
665 :
666 39 : TList **tracksZSplitted = new TList*[2];
667 234 : for (Int_t i = 0; i < 2; i++)
668 156 : tracksZSplitted[i] = new TList;
669 :
670 39 : TList *tracksZMergedStage1 = new TList;
671 :
672 39 : AliTRDtrackGTU **trkInRefLayer = new AliTRDtrackGTU*[fGtuParam->GetNRefLayers()];
673 :
674 : // Bool_t done = kFALSE;
675 : Int_t minIdx = 0;
676 : AliTRDtrackGTU *trkStage0 = 0x0;
677 :
678 312 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
679 : // ----- Merging and Unification in Reflayers (seed_merger) -----
680 : do {
681 : // done = kTRUE;
682 : trkStage0 = 0x0;
683 1192 : for (Int_t refLayerIdx = 0; refLayerIdx < fGtuParam->GetNRefLayers(); refLayerIdx++) {
684 447 : trkInRefLayer[refLayerIdx] = (AliTRDtrackGTU*) fTracks[zch][refLayerIdx].First();
685 447 : if (trkInRefLayer[refLayerIdx] == 0) {
686 : continue;
687 : }
688 33 : else if (trkStage0 == 0x0 ) {
689 : trkStage0 = trkInRefLayer[refLayerIdx];
690 : minIdx = refLayerIdx;
691 : // done = kFALSE;
692 32 : }
693 2 : else if ( (trkInRefLayer[refLayerIdx]->GetZSubChannel() < trkStage0->GetZSubChannel()) ||
694 1 : ((trkInRefLayer[refLayerIdx]->GetZSubChannel() == trkStage0->GetZSubChannel()) &&
695 1 : ((trkInRefLayer[refLayerIdx]->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
696 : minIdx = refLayerIdx;
697 0 : trkStage0 = trkInRefLayer[refLayerIdx];
698 : // done = kFALSE;
699 0 : }
700 : }
701 149 : if (!trkStage0)
702 : break;
703 32 : tracksRefMerged[zch]->Add(trkStage0);
704 32 : fTracks[zch][minIdx].RemoveFirst();
705 32 : } while (trkStage0 != 0);
706 :
707 117 : Uniquifier(tracksRefMerged[zch], tracksRefUnique[zch]);
708 :
709 351 : AliDebug(2, Form("zchannel %i", zch));
710 117 : TIter trackRefMerged(tracksRefMerged[zch]);
711 532 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefMerged())
712 192 : 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 585 : AliDebug(2, "uniquified:");
723 117 : TIter trackRefUnique(tracksRefUnique[zch]);
724 399 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackRefUnique())
725 144 : 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()));
735 117 : }
736 :
737 : // ----- Merging in zchannels - 1st stage -----
738 :
739 39 : if (AliTRDgtuParam::GetUseGTUmerge()) {
740 : Int_t notEmpty;
741 39 : do {
742 63 : Bool_t lowerThan[3] = { kFALSE, kFALSE, kFALSE };
743 63 : AliTRDtrackGTU *trk[3] = { 0x0, 0x0, 0x0 };
744 504 : for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel)
745 189 : trk[iChannel] = (AliTRDtrackGTU*) tracksRefUnique[iChannel]->First();
746 :
747 504 : for (Int_t iChannel = 0; iChannel < fGtuParam->GetNZChannels(); ++iChannel) {
748 189 : AliTRDtrackGTU *trk1 = trk[iChannel];
749 189 : AliTRDtrackGTU *trk2 = trk[(iChannel + 1) % 3];
750 189 : if (trk1 && trk2) {
751 11 : Int_t sortnum1 = (trk1->GetZChannel() + 3 * trk1->GetZSubChannel()) / 2 - 1;
752 11 : Int_t sortnum2 = (trk2->GetZChannel() + 3 * trk2->GetZSubChannel()) / 2 - 1;
753 33 : AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
754 : trk1->GetZChannel(), trk2->GetZChannel(),
755 : sortnum1, sortnum2));
756 15 : if ( (sortnum1 < sortnum2) ||
757 5 : ((sortnum1 == sortnum2) &&
758 4 : ((trk1->GetYapprox() >> 3) < (trk2->GetYapprox() >> 3)) ) ) {
759 7 : lowerThan[iChannel] = kTRUE;
760 7 : }
761 11 : }
762 : }
763 :
764 174 : notEmpty = (trk[2] ? (1 << 2) : 0) |
765 174 : (trk[1] ? (1 << 1) : 0) |
766 87 : (trk[0] ? (1 << 0) : 0);
767 : Int_t pop = -1;
768 :
769 87 : switch (notEmpty) {
770 : // one track only
771 : case 0x1:
772 : pop = 0;
773 7 : break;
774 : case 0x2:
775 : pop = 1;
776 4 : break;
777 : case 0x4:
778 : pop = 2;
779 4 : break;
780 :
781 : // two tracks
782 : case 0x3:
783 2 : if (lowerThan[0])
784 1 : pop = 0;
785 : else
786 : pop = 1;
787 : break;
788 : case 0x5:
789 3 : if (lowerThan[2])
790 3 : pop = 2;
791 : else
792 : pop = 0;
793 : break;
794 : case 0x6:
795 3 : if (lowerThan[1])
796 2 : pop = 1;
797 : else
798 : pop = 2;
799 : break;
800 :
801 : // three tracks
802 : case 0x7:
803 1 : if (lowerThan[0]) {
804 1 : if (lowerThan[2])
805 0 : pop = 2;
806 : else
807 : pop = 0;
808 : } else {
809 0 : if (lowerThan[1])
810 0 : pop = 1;
811 : else
812 : pop = 2;
813 : }
814 : break;
815 :
816 : // no tracks
817 : default:
818 : // nop
819 : ;
820 : }
821 :
822 63 : if (pop > -1) {
823 24 : tracksZMergedStage0->Add(trk[pop]);
824 24 : tracksRefUnique[pop]->RemoveFirst();
825 24 : }
826 63 : } while (notEmpty);
827 39 : }
828 : else {
829 : // there is still a difference between this implementation and
830 : // the hardware algorithm, only for expert use
831 :
832 : do {
833 : // done = kTRUE;
834 : trkStage0 = 0x0;
835 : // compare tracks from all adjacent zchannels
836 : // (including comparison of channels 2 and 0)
837 0 : for (Int_t i = fGtuParam->GetNZChannels() - 1; i > -1; i--) {
838 0 : Int_t zch = i % 3;
839 0 : AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksRefUnique[zch]->First();
840 0 : if (trk == 0) {
841 0 : continue;
842 : }
843 0 : else if (trkStage0 == 0x0 ) {
844 : trkStage0 = trk;
845 : minIdx = zch;
846 : // done = kFALSE;
847 0 : }
848 : else {
849 0 : Int_t sortnum1 = (trk->GetZChannel() + 3 * trk->GetZSubChannel()) / 2 - 1;
850 0 : Int_t sortnum2 = (trkStage0->GetZChannel() + 3 * trkStage0->GetZSubChannel()) / 2 - 1;
851 0 : AliDebug(5, Form("comparing tracks %i - %i: %i - %i",
852 : trk->GetZChannel(), trkStage0->GetZChannel(),
853 : sortnum1, sortnum2));
854 0 : if ( (sortnum1 < sortnum2) ||
855 0 : ((sortnum1 == sortnum2) &&
856 0 : ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
857 : minIdx = zch;
858 : trkStage0 = trk;
859 : // done = kFALSE;
860 0 : }
861 : }
862 0 : }
863 :
864 0 : if (!trkStage0)
865 : break;
866 0 : tracksZMergedStage0->Add(trkStage0);
867 0 : tracksRefUnique[minIdx]->RemoveFirst();
868 0 : } while (trkStage0 != 0);
869 : }
870 :
871 39 : Uniquifier(tracksZMergedStage0, tracksZUniqueStage0);
872 :
873 117 : AliDebug(2, "stage 0:");
874 39 : TIter trackZMergedStage0(tracksZMergedStage0);
875 204 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage0())
876 144 : 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,
885 : trk->GetZChannel(),
886 : trk->GetZSubChannel()));
887 :
888 195 : AliDebug(2, "uniquified:");
889 39 : TIter trackZUniqueStage0(tracksZUniqueStage0);
890 147 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZUniqueStage0())
891 90 : 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,
900 : trk->GetZChannel(),
901 : trk->GetZSubChannel()));
902 :
903 : // ----- Splitting in z -----
904 :
905 39 : TIter next(tracksZUniqueStage0);
906 147 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) next()) {
907 30 : if ((trk->GetZChannel() < 0) ||
908 15 : (trk->GetZChannel() > 2) ||
909 30 : (trk->GetZSubChannel() < 0) ||
910 30 : (trk->GetZSubChannel() > 6)) {
911 0 : AliError(Form("track with invalid z-channel information at %p: zch = %i, subchannel = %i",
912 : trk, trk->GetZChannel(), trk->GetZSubChannel()));
913 0 : trk->Dump();
914 : }
915 30 : Int_t idx = (trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) % 2;
916 15 : if ((idx < 0) || (idx > 1)) {
917 0 : AliError(Form("invalid index %i null", idx));
918 0 : trk->Dump();
919 0 : continue;
920 : }
921 15 : if (!tracksZSplitted[idx]) {
922 0 : AliError(Form("array pointer %i null", idx));
923 0 : continue;
924 : }
925 15 : tracksZSplitted[idx]->Add(trk);
926 30 : }
927 :
928 234 : for (Int_t i = 0; i < 2; i++) {
929 390 : AliDebug(2, Form("split %i:", i));
930 78 : TIter trackZSplit(tracksZSplitted[i]);
931 264 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZSplit())
932 90 : 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,
941 : trk->GetZChannel(),
942 : trk->GetZSubChannel()));
943 78 : }
944 :
945 : // ----- Merging in zchanels - 2nd stage -----
946 :
947 39 : do {
948 : // done = kTRUE;
949 : trkStage0 = 0x0;
950 324 : for (Int_t i = 1; i >= 0; i--) {
951 216 : AliTRDtrackGTU *trk = (AliTRDtrackGTU*) tracksZSplitted[i]->First();
952 108 : if (trk == 0) {
953 93 : continue;
954 : }
955 15 : else if (trkStage0 == 0x0 ) {
956 : trkStage0 = trk;
957 : minIdx = i;
958 : // done = kFALSE;
959 15 : }
960 0 : else if ( (((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) < ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) ||
961 0 : ((((trk->GetZChannel() + 3 * (trk->GetZSubChannel() - 1)) / 2) == ((trkStage0->GetZChannel() + 3 * (trkStage0->GetZSubChannel() - 1)) / 2)) &&
962 0 : ((trk->GetYapprox() >> 3) < (trkStage0->GetYapprox() >> 3)) ) ) {
963 : minIdx = i;
964 : trkStage0 = trk;
965 : // done = kFALSE;
966 0 : }
967 15 : }
968 :
969 54 : if (!trkStage0)
970 : break;
971 15 : tracksZMergedStage1->Add(trkStage0);
972 15 : tracksZSplitted[minIdx]->RemoveFirst();
973 15 : } while (trkStage0 != 0);
974 :
975 39 : Uniquifier(tracksZMergedStage1, ListOfTracks);
976 :
977 195 : AliDebug(2, "merged:");
978 39 : TIter trackZMergedStage1(tracksZMergedStage1);
979 147 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) trackZMergedStage1())
980 90 : 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,
989 : trk->GetZChannel(),
990 : trk->GetZSubChannel()));
991 :
992 195 : AliDebug(2, "uniquified:");
993 39 : TIter track(ListOfTracks);
994 147 : while (AliTRDtrackGTU *trk = (AliTRDtrackGTU*) track())
995 90 : 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,
1004 : trk->GetZChannel(),
1005 : trk->GetZSubChannel()));
1006 :
1007 : // cleaning up
1008 312 : for (Int_t zch = 0; zch < fGtuParam->GetNZChannels(); zch++) {
1009 234 : delete tracksRefMerged[zch];
1010 234 : delete tracksRefUnique[zch];
1011 : }
1012 78 : delete [] tracksRefMerged;
1013 78 : delete [] tracksRefUnique;
1014 :
1015 :
1016 78 : delete tracksZMergedStage0;
1017 78 : delete tracksZUniqueStage0;
1018 :
1019 234 : for (Int_t i = 0; i < 2; i++)
1020 156 : delete tracksZSplitted[i];
1021 78 : delete [] tracksZSplitted;
1022 :
1023 78 : delete tracksZMergedStage1;
1024 :
1025 78 : delete [] trkInRefLayer;
1026 :
1027 : return kTRUE;
1028 39 : }
1029 :
1030 : Bool_t AliTRDgtuTMU::RunTrackReconstruction(TList* ListOfTracks)
1031 : {
1032 : // run the track reconstruction for all tracks in the list
1033 :
1034 78 : TIter next(ListOfTracks);
1035 186 : while (AliTRDtrackGTU *track = (AliTRDtrackGTU*) next()) {
1036 15 : CalculateTrackParams(track);
1037 15 : CalculatePID(track);
1038 75 : AliDebug(1, Form("track with pid = %i", track->GetPID()));
1039 15 : }
1040 : return kTRUE;
1041 39 : }
1042 :
1043 : Bool_t AliTRDgtuTMU::CalculatePID(AliTRDtrackGTU *track)
1044 : {
1045 : // calculate PID for the given track
1046 30 : if (!track) {
1047 0 : AliError("No track to calculate!");
1048 0 : return kFALSE;
1049 : }
1050 :
1051 15 : if (AliTRDgtuParam::GetUseGTUconst()) {
1052 : // averaging as in GTU
1053 45 : AliDebug(1, "using GTU constants for PID calculation");
1054 : ULong64_t coeff;
1055 :
1056 : // selection of coefficient for averaging
1057 15 : switch(track->GetTrackletMask()){
1058 : case 0x3f:
1059 : // 6 tracklets
1060 : coeff=0x5558; // approx. 1/6
1061 3 : break;
1062 :
1063 : case 0x3e:
1064 : case 0x3d:
1065 : case 0x3b:
1066 : case 0x37:
1067 : case 0x2f:
1068 : case 0x1f:
1069 : // 5 tracklets
1070 : coeff=0x6666; // approx. 1/5
1071 7 : break;
1072 :
1073 : default:
1074 : // 4 tracklets
1075 : coeff=0x8000; // = 0.25
1076 5 : }
1077 15 : coeff &= 0x1ffff; // 17-bit constant
1078 :
1079 : ULong64_t sum = 0;
1080 : Int_t i = 0;
1081 210 : for (Int_t iLayer = 0; iLayer < fGtuParam->GetNLayers(); iLayer++) {
1082 90 : if ((track->GetTrackletMask() >> iLayer) & 1) {
1083 73 : sum += track->GetTracklet(iLayer)->GetPID();
1084 73 : ++i;
1085 73 : }
1086 : }
1087 :
1088 15 : Float_t av = 1./i * sum;
1089 15 : sum = sum & 0x7ff;
1090 15 : ULong64_t prod = (sum * coeff) & 0xfffffffffull;
1091 15 : ULong64_t prodFinal = ((prod >> 17) + ((prod >> 16) & 1)) & 0xff;
1092 :
1093 15 : if (TMath::Abs((prodFinal & 0xff) - av) > 0.5)
1094 0 : AliError(Form("problem with PID averaging (hw <-> ar): %3lli <-> %4.1f", prodFinal & 0xff, av));
1095 15 : track->SetPID(prodFinal & 0xff);
1096 :
1097 : return kTRUE;
1098 : }
1099 : else {
1100 :
1101 : // simple averaging
1102 : Int_t nTracklets = 0;
1103 : Int_t pidSum = 0;
1104 0 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1105 0 : if (!track->IsTrackletInLayer(layer)) {
1106 : continue;
1107 : }
1108 0 : AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1109 0 : pidSum += trk->GetPID();
1110 0 : nTracklets++;
1111 0 : }
1112 :
1113 0 : if (nTracklets>0)
1114 0 : track->SetPID(pidSum/nTracklets);
1115 : else
1116 0 : AliError("Track without contributing tracklets, no PID assigned");
1117 :
1118 : return kTRUE;
1119 : }
1120 15 : }
1121 :
1122 : Bool_t AliTRDgtuTMU::CalculateTrackParams(AliTRDtrackGTU *track)
1123 : {
1124 30 : Int_t corrMode = fGtuParam->GetCorrectionMode(); //correct y position for sagitta method only
1125 :
1126 : // calculate the track parameters
1127 :
1128 15 : if (!track) {
1129 0 : AliError("No track to calculate!");
1130 0 : return kFALSE;
1131 : }
1132 :
1133 : Int_t a = 0;
1134 : Float_t b = 0;
1135 : Float_t c = 0;
1136 15 : Float_t x1;
1137 15 : Float_t x2;
1138 : Int_t s = 0;
1139 : Int_t sTmp = 0; // only temporary until optimization methods are chosen
1140 : Int_t d = 0;
1141 : Int_t invPtDev = 0;
1142 :
1143 : //additional params for sagitta optimization
1144 15 : Int_t tmpTrackletMask = track->GetTrackletMask();
1145 15 : Int_t trklYpos[6] = { 0 };
1146 : Int_t invXpos = 0;
1147 : Int_t stack = 0;
1148 : Int_t zBin = 0;
1149 : Int_t zPos = 0;
1150 : Int_t fitParam = 0;
1151 : Int_t zDiff = 0;
1152 : Int_t padTiltingCorrection = 0;
1153 : Int_t refLayer = 0; //reference layer for pad tilting correction
1154 26 : if (track->GetTracklet(3)!=0x0) refLayer = 3;
1155 8 : else if (track->GetTracklet(4)!=0x0) refLayer = 4;
1156 0 : else if (track->GetTracklet(2)!=0x0) refLayer = 2;
1157 15 : invXpos = fGtuParam->GetLayerInvXpos(refLayer);
1158 15 : stack = fGtuParam->GetGeo()->GetStack(track->GetTracklet(refLayer)->GetDetector());
1159 15 : zBin = track->GetTracklet(refLayer)->GetZbin();
1160 15 : zPos = fGtuParam->GetZpos(stack, refLayer, zBin);
1161 15 : fitParam = zPos*invXpos;
1162 15 : fitParam *= 1e-4;
1163 :
1164 45 : AliDebug(5,Form("There are %i tracklets in this track.", track->GetNTracklets()));
1165 :
1166 15 : if ( (corrMode == 1) || (corrMode == 3) )
1167 : {
1168 : // tracklet removal
1169 0 : if (track->GetTrackletMask() == 31) tmpTrackletMask = 30;
1170 0 : else if (track->GetTrackletMask() == 47) tmpTrackletMask = 45;
1171 0 : else if (track->GetTrackletMask() == 55) tmpTrackletMask = 39;
1172 0 : else if (track->GetTrackletMask() == 59) tmpTrackletMask = 57;
1173 0 : else if (track->GetTrackletMask() == 61) tmpTrackletMask = 45;
1174 0 : else if (track->GetTrackletMask() == 62) tmpTrackletMask = 30;
1175 : }
1176 :
1177 210 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1178 90 : if (!track->IsTrackletInLayer(layer)) {
1179 : continue;
1180 : }
1181 73 : AliTRDtrackletGTU *trk = track->GetTracklet(layer);
1182 73 : if (!trk) {
1183 0 : AliError(Form("Could not get tracklet in layer %i\n", layer));
1184 0 : continue;
1185 : }
1186 73 : Int_t stackTmp = fGtuParam->GetGeo()->GetStack(trk->GetDetector());
1187 73 : Int_t zBinTmp = trk->GetZbin();
1188 131 : if (layer != refLayer) zDiff = fGtuParam->GetZpos(stackTmp, layer, zBinTmp)*1e3 - fitParam*fGtuParam->GetLayerXpos(layer);
1189 146 : trklYpos[layer] = trk->GetYbin() - ((Int_t) (1./fGtuParam->GetBinWidthY()*fGtuParam->CorrectYforAlignmentOCDB(trk->GetDetector(),
1190 73 : fGtuParam->GetZpos(fGtuParam->GetGeo()->GetStack(trk->GetDetector()), layer, trk->GetZbin())))); //ocdb alignment correction
1191 73 : padTiltingCorrection = zDiff * fGtuParam->GetTanOfTiltingAngle(layer);
1192 73 : padTiltingCorrection *= 6.25e-7; // *=1e-8 to account for previous shifts and /=160e-4 for the bin width of yPos
1193 73 : if ( (corrMode == 2) || (corrMode == 3) ) trklYpos[layer] = trklYpos[layer] + padTiltingCorrection;
1194 :
1195 :
1196 219 : AliDebug(10,Form(" layer %i trk yprime: %6i, aki: %6i, pki: %i, ybin: %6i", layer, trk->GetYPrime(),
1197 : fGtuParam->GetAki(track->GetTrackletMask(), layer), fGtuParam->GetPki(track->GetTrackletMask(), layer), trk->GetYbin()));
1198 73 : a += (((fGtuParam->GetAki(track->GetTrackletMask(), layer) * trk->GetYPrime()) >> 7) + 1) >> 1;
1199 73 : b += fGtuParam->GetBki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1200 73 : c += fGtuParam->GetCki(track->GetTrackletMask(), layer) * trk->GetYPrime() * fGtuParam->GetBinWidthY();
1201 73 : s += fGtuParam->GetPki(track->GetTrackletMask(), layer) * trk->GetYbin() ;
1202 73 : sTmp += fGtuParam->GetPki(tmpTrackletMask, layer) * trklYpos[layer];
1203 73 : }
1204 :
1205 15 : s *= fGtuParam->GetLengthNorm(track->GetTrackletMask());
1206 15 : sTmp *= fGtuParam->GetLengthNorm(tmpTrackletMask);
1207 :
1208 15 : if (a < 0)
1209 3 : a += 3;
1210 15 : a = a >> 2;
1211 :
1212 15 : invPtDev = a * fGtuParam->Getc1Inv(track->GetTrackletMask()) - s;
1213 :
1214 15 : if (!fGtuParam->GetWriteSagittaOutputToTrackWordBC())
1215 15 : track->SetFitParams(a << 2, TMath::Nint(128. * b), TMath::Nint(256. * c));
1216 : else {
1217 0 : c = ((invPtDev & 0xfff) ^ 0x800) - 0x800;
1218 0 : b = (invPtDev >> 12) & 0xfff;
1219 0 : if (TMath::Abs(invPtDev) < fGtuParam->GetInvPtDevCut()) b += 1 << 12;
1220 0 : track->SetFitParams(a << 2, b, c);
1221 : }
1222 30 : if (corrMode == 0) track->SetInvPtDev(invPtDev);
1223 0 : else if (corrMode != 0) track->SetInvPtDev( a * fGtuParam->Getc1Inv(track->GetTrackletMask()) - ((Int_t) sTmp) );
1224 : //following lines are for debugging purposes only
1225 : //if (corrMode == 0) track->SetInvPtDev( s );
1226 : //else if (corrMode == 1) track->SetInvPtDev( (Int_t) sTmp ) ;
1227 :
1228 15 : fGtuParam->GetIntersectionPoints(track->GetTrackletMask(), x1, x2);
1229 :
1230 45 : AliDebug(5,Form(" Track parameters: a16 = %i, a18 = %i, b = %f, c = %f, x1 = %f, x2 = %f, pt = %f, s = %i, invPtDev = %i (trkl mask: %i)",
1231 : a, a << 2, b, c, x1, x2, track->GetPt(), s, invPtDev, track->GetTrackletMask()));
1232 :
1233 : return kTRUE;
1234 30 : }
1235 :
1236 : Bool_t AliTRDgtuTMU::Uniquifier(const TList *inlist, TList *outlist)
1237 : {
1238 : // remove multiple occurences of the same track
1239 :
1240 390 : TIter next(inlist);
1241 : AliTRDtrackGTU *trkStage0 = 0x0;
1242 : AliTRDtrackGTU *trkStage1 = 0x0;
1243 :
1244 195 : do {
1245 532 : trkStage0 = (AliTRDtrackGTU*) next();
1246 :
1247 : Bool_t tracksEqual = kFALSE;
1248 266 : if (trkStage0 != 0 && trkStage1 != 0) {
1249 42 : for (Int_t layer = 0; layer < fGtuParam->GetNLayers(); layer++) {
1250 93 : if (trkStage0->GetTrackletIndex(layer) != -1 && trkStage0->GetTrackletIndex(layer) == trkStage1->GetTrackletIndex(layer)) {
1251 : tracksEqual = kTRUE;
1252 17 : break;
1253 : }
1254 : }
1255 17 : }
1256 :
1257 266 : if (tracksEqual) {
1258 51 : if (trkStage0->GetNTracklets() >= trkStage1->GetNTracklets())
1259 16 : trkStage1 = trkStage0;
1260 : }
1261 : else {
1262 249 : if (trkStage1 != 0x0)
1263 54 : outlist->Add(trkStage1);
1264 : trkStage1 = trkStage0;
1265 : }
1266 :
1267 266 : } while (trkStage1 != 0x0);
1268 : return kTRUE;
1269 195 : }
|
