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