Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-2012, 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 : ///////////////////////////////////////////////////////////////////////////////
17 : //
18 : // Track matching between TRD online tracks and ESD tracks.
19 : //
20 : // Author: Felix Rettig <rettig@compeng.uni-frankfurt.de>
21 : //
22 : ///////////////////////////////////////////////////////////////////////////////
23 :
24 : #include <TH1.h>
25 : #include <AliESDEvent.h>
26 : #include <AliExternalTrackParam.h>
27 : #include "AliESDtrack.h"
28 : #include "AliESDTrdTrack.h"
29 : #include <AliGeomManager.h>
30 : #include "AliTRDgeometry.h"
31 : #include "AliTRDpadPlane.h"
32 : #include "AliTRDonlineTrackMatching.h"
33 :
34 : const Float_t AliTRDonlineTrackMatching::fgkSaveInnerRadius = 290.5;
35 : const Float_t AliTRDonlineTrackMatching::fgkSaveOuterRadius = 364.5;
36 :
37 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMinTPCrows = 0.;
38 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMinRatioRowsFindableClusters = 0.;
39 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxChi2TPCclusters = 0.;
40 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxChi2ITSclusters = 0.;
41 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxDCAtoVertexXY = 0.;
42 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxDCAtoVertexZ = 0.;
43 : UShort_t AliTRDonlineTrackMatching::fEsdTrackCutsITSlayerMask = 0; // similar to 2011 default cut: 0x3
44 : Float_t AliTRDonlineTrackMatching::fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 0.;
45 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutPtDCAOfs = 0.;
46 : Float_t AliTRDonlineTrackMatching::fEsdTrackCutPtDCACoeff = 0.;
47 : Bool_t AliTRDonlineTrackMatching::fEsdTrackCutMinimal = kFALSE;
48 : Bool_t AliTRDonlineTrackMatching::fEsdTrackCutRequireTPCrefit = kTRUE;
49 : Bool_t AliTRDonlineTrackMatching::fEsdTrackCutRequireITSrefit = kFALSE;
50 : Bool_t AliTRDonlineTrackMatching::fEsdTrackCutPrim = kFALSE;
51 :
52 : AliTRDonlineTrackMatching::AliTRDonlineTrackMatching() :
53 24 : TObject(),
54 24 : fTRDgeo(NULL),
55 24 : fMinMatchRating(0.25),
56 24 : fHistMatchRating(NULL)
57 120 : {
58 : // default ctor
59 24 : SetEsdTrackDefaultCuts("minimal");
60 48 : }
61 :
62 : AliTRDonlineTrackMatching::AliTRDonlineTrackMatching(const AliTRDonlineTrackMatching &c) :
63 0 : TObject(c),
64 0 : fTRDgeo(c.fTRDgeo),
65 0 : fMinMatchRating(c.fMinMatchRating),
66 0 : fHistMatchRating(c.fHistMatchRating)
67 0 : {
68 : // copy ctor
69 0 : }
70 :
71 96 : AliTRDonlineTrackMatching::~AliTRDonlineTrackMatching() {
72 :
73 : // dtor
74 :
75 26 : delete fTRDgeo;
76 24 : fTRDgeo = NULL;
77 48 : }
78 :
79 : Short_t AliTRDonlineTrackMatching::EstimateSector(const Double_t globalCoords[3]) {
80 :
81 : // estimates sector by phi angle in x-y plane
82 :
83 7288 : if ((TMath::Abs(globalCoords[0]) > 600) || (TMath::Abs(globalCoords[0]) > 600) || (TMath::Sqrt(globalCoords[0]*globalCoords[0] + globalCoords[1]*globalCoords[1]) < 0.01)){
84 : //printf("GGG %.3f/%.3f\n", globalCoords[0], globalCoords[1]);
85 0 : return -1;
86 : } else {
87 1822 : Double_t ang = TMath::ATan2(globalCoords[1], globalCoords[0]);
88 1822 : if (ang > 0){
89 : #ifdef TRD_TM_DEBUG
90 : printf(" es: %.2f/%.2f -> phi: %.2fdeg -> Sec %02d (A)\n",
91 : globalCoords[0], globalCoords[1], TMath::ATan2(globalCoords[1], globalCoords[0])*180./TMath::Pi(),
92 : TMath::FloorNint(ang/(20./180.*TMath::Pi())));
93 : #endif
94 1426 : return TMath::FloorNint(ang/(20./180.*TMath::Pi()));
95 : } else {
96 : #ifdef TRD_TM_DEBUG
97 : printf(" es: %.2f/%.2f -> phi: %.2fdeg -> Sec %02d (B)\n",
98 : globalCoords[0], globalCoords[1], TMath::ATan2(globalCoords[1], globalCoords[0])*180./TMath::Pi(),
99 : 17 - TMath::FloorNint(TMath::Abs(ang)/(20./180.*TMath::Pi())));
100 : #endif
101 396 : return 17 - TMath::FloorNint(TMath::Abs(ang)/(20./180.*TMath::Pi()));
102 : }
103 :
104 : }
105 1822 : }
106 :
107 : Short_t AliTRDonlineTrackMatching::EstimateLayer(Double_t radius) {
108 :
109 : // estimates layer by radial distance (for virtual stack at phi = 0)
110 :
111 : const Float_t rBoundaries[7] = {290.80, 302.20, 315.06, 327.55, 340.3, 352.80, 364.15}; // radial border lines centered between anode plane and successing radiator
112 : const Short_t rLayers[7] = {-1, 0, 1, 2, 3, 4, 5};
113 10422 : for (UShort_t i = 0; i < 7; ++i){
114 4300 : if (radius < rBoundaries[i])
115 1822 : return rLayers[i];
116 : }
117 0 : return -2; // radius larger than outmost layer
118 1822 : }
119 :
120 : Short_t AliTRDonlineTrackMatching::EstimateLocalStack(const Double_t globalCoords[3]) {
121 :
122 : // determines stack within sector by z position
123 :
124 3644 : Double_t absZ = TMath::Abs(globalCoords[2]);
125 1822 : Short_t signZ = (globalCoords[2] > 0.) ? 1 : -1;
126 1822 : Double_t r = TMath::Sqrt(globalCoords[0]*globalCoords[0] + globalCoords[1]*globalCoords[1]);
127 1822 : Short_t layer = EstimateLayer(r);
128 :
129 : #ifdef TRD_TM_DEBUG
130 : printf("EstimateLocalStack A r: %.2f x: %.2f/%.2f/%.2f -> layer: %i absZ = %.2f\n",
131 : r, globalCoords[0], globalCoords[1], globalCoords[2], layer, absZ);
132 : #endif
133 :
134 1822 : if (layer < 0)
135 34 : return -1;
136 :
137 1788 : Double_t innerStackHalfLength = AliTRDgeometry::GetChamberLength(0, 2) / 2.; // same for all layers
138 1788 : if (absZ < innerStackHalfLength)
139 1142 : return 2;
140 :
141 646 : Double_t outerStackLength = AliTRDgeometry::GetChamberLength(layer, 1);
142 :
143 646 : absZ -= innerStackHalfLength;
144 :
145 : #ifdef TRD_TM_DEBUG
146 : printf("EstimateLocalStack B r: %.2f x: %.2f/%.2f/%.2f -> layer: %i absZ = %.2f il: %.2f ol: %.2f\n",
147 : r, globalCoords[0], globalCoords[1], globalCoords[2], layer, absZ, 2.*innerStackHalfLength, outerStackLength);
148 : #endif
149 :
150 646 : if (absZ > 2.05*outerStackLength)
151 0 : return (signZ > 0) ? -2 : -1; // outside supermodule in z direction
152 :
153 1292 : if (absZ < outerStackLength)
154 1190 : return (signZ > 0) ? 1 : 3;
155 : else
156 102 : return (signZ > 0) ? 0 : 4;
157 :
158 1822 : }
159 :
160 : Short_t AliTRDonlineTrackMatching::EstimateStack(const Double_t globalCoords[3]) {
161 :
162 : // returns the closest TRD stack to a 3D position in global coordinates
163 :
164 3644 : Short_t sec = EstimateSector(globalCoords);
165 1822 : Short_t st = EstimateLocalStack(globalCoords);
166 : #ifdef TRD_TM_DEBUG
167 : printf("EstimateStack sec %d st %d\n", sec, st);
168 : #endif
169 3644 : if ((sec < 0) || (st < 0))
170 34 : return -1;
171 : else
172 1788 : return 5*sec + st;
173 1822 : }
174 :
175 : Bool_t AliTRDonlineTrackMatching::StackToTrack(const AliExternalTrackParam *track, Short_t &stack, UShort_t &layersWithTracklet, Double_t magFieldinKiloGauss){
176 :
177 : // returns stack to track param
178 :
179 260 : stack = -1;
180 130 : layersWithTracklet = 0;
181 :
182 130 : UInt_t stackHits[fgkTrdStacks];
183 130 : Double_t x[3] = { 0. };
184 130 : memset(stackHits, 0, fgkTrdStacks*sizeof(UInt_t));
185 :
186 : #ifdef TRD_TM_DEBUG
187 : printf("STACK-TO-TRACK\n");
188 : #endif
189 :
190 : Double_t r = fgkSaveInnerRadius;
191 7504 : while (r < fgkSaveOuterRadius){
192 3726 : if (track->GetXYZAt(r, magFieldinKiloGauss, x)) {
193 1822 : stack = EstimateStack(x);
194 1822 : if (stack >= 0){
195 1788 : stackHits[stack]++;
196 1788 : if (stackHits[stack] > 16) // experimental
197 : break;
198 : #ifdef TRD_TM_DEBUG
199 : printf(" r=%.3fcm %.2f/%.2f - %d hits for stack %d S%02d-%d (mag=%.1f)\n",
200 : r, x[0], x[1], stackHits[stack], stack, stack/5, stack%5, magFieldinKiloGauss);
201 : #endif
202 : }
203 : }
204 3622 : r += 1.;
205 : }
206 :
207 : // find stack with most hits
208 : UInt_t bestHits = 0;
209 23660 : for (UShort_t iStack = 0; iStack < fgkTrdStacks; ++iStack){
210 11700 : if (stackHits[iStack] == 0)
211 : continue;
212 : #ifdef TRD_TM_DEBUG
213 : printf(" finally %d hits in stack S%02d-%d\n", stackHits[iStack], iStack/5, iStack%5);
214 : #endif
215 108 : if (stackHits[iStack] > bestHits){
216 : bestHits = stackHits[iStack];
217 108 : stack = iStack;
218 108 : }
219 : }
220 :
221 130 : if (stack >= 0){
222 : #ifdef TRD_TM_DEBUG
223 : printf("best stack: S%02d-%d\n", TrdLsiSec(stack), TrdLsiSi(stack));
224 : #endif
225 106 : return kTRUE;
226 : }
227 :
228 24 : return kFALSE;
229 130 : }
230 :
231 : Bool_t AliTRDonlineTrackMatching::StackToTrack(const AliESDtrack* track, Short_t &stack, UShort_t &layersWithTracklet, Double_t magFieldinKiloGauss){
232 :
233 : // returns stack to ESD track
234 :
235 390 : if (track->GetOuterParam())
236 260 : return StackToTrack(track->GetOuterParam(), stack, layersWithTracklet, magFieldinKiloGauss);
237 0 : else if (track->GetInnerParam())
238 0 : return StackToTrack(track->GetInnerParam(), stack, layersWithTracklet, magFieldinKiloGauss);
239 : else
240 0 : return StackToTrack(track, stack, layersWithTracklet, magFieldinKiloGauss);
241 130 : }
242 :
243 : Bool_t AliTRDonlineTrackMatching::AcceptTrack(const AliESDtrack* esdTrack, const AliESDEvent* esdEvent){
244 :
245 : // returns result ESD track cuts
246 :
247 304 : if (!esdTrack)
248 0 : return kFALSE;
249 :
250 152 : UInt_t status = esdTrack->GetStatus();
251 :
252 152 : if (fEsdTrackCutMinimal){
253 0 : return ((status & AliESDtrack::kTPCout) > 0);
254 : }
255 :
256 : // require TPC fit
257 304 : if ((fEsdTrackCutRequireTPCrefit) && (!(status & AliESDtrack::kTPCrefit)))
258 16 : return kFALSE;
259 :
260 : // require ITS re-fit
261 136 : if ((fEsdTrackCutRequireITSrefit) && (!(status & AliESDtrack::kITSrefit)))
262 0 : return kFALSE;
263 :
264 : // TPC requirements
265 136 : Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
266 : Float_t ratioCrossedRowsOverFindableClustersTPC =
267 408 : (esdTrack->GetTPCNclsF() > 0) ? (nCrossedRowsTPC / esdTrack->GetTPCNclsF()) : 1.0;
268 : Float_t chi2PerClusterTPC =
269 408 : (esdTrack->GetTPCclusters(0) > 0) ? (esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCclusters(0))) : 100.;
270 :
271 : if (
272 266 : (nCrossedRowsTPC < fEsdTrackCutMinTPCrows) ||
273 130 : (ratioCrossedRowsOverFindableClustersTPC < fEsdTrackCutMinRatioRowsFindableClusters) ||
274 130 : (chi2PerClusterTPC > fEsdTrackCutMaxChi2TPCclusters)
275 : )
276 6 : return kFALSE;
277 :
278 : // ITS requirements
279 342 : Float_t chi2PerClusterITS = (esdTrack->GetITSclusters(0) > 0) ? esdTrack->GetITSchi2()/Float_t(esdTrack->GetITSclusters(0)) : 1000.;
280 : UShort_t clustersInAnyITSlayer = kFALSE;
281 1820 : for (UShort_t layer = 0; layer < 6; ++layer)
282 780 : clustersInAnyITSlayer += (esdTrack->HasPointOnITSLayer(layer) & ((fEsdTrackCutsITSlayerMask >> layer) & 1));
283 :
284 130 : if ((fEsdTrackCutsITSlayerMask != 0) &&
285 0 : ((clustersInAnyITSlayer == 0) || (chi2PerClusterITS >= fEsdTrackCutMaxChi2ITSclusters))
286 : )
287 0 : return kFALSE;
288 :
289 : // geometric requirements
290 130 : Float_t impactPos[2], impactCov[3];
291 130 : esdTrack->GetImpactParameters(impactPos, impactCov);
292 :
293 130 : if (TMath::Abs(impactPos[0]) > fEsdTrackCutMaxDCAtoVertexXY)
294 0 : return kFALSE;
295 :
296 130 : if (TMath::Abs(impactPos[1]) > fEsdTrackCutMaxDCAtoVertexZ)
297 0 : return kFALSE;
298 :
299 130 : if (fEsdTrackCutPrim){
300 : // additional requirements for primary tracks
301 :
302 0 : const AliESDVertex* vertex = esdEvent->GetPrimaryVertexTracks();
303 0 : if ((!vertex) || (!vertex->GetStatus()))
304 0 : vertex = esdEvent->GetPrimaryVertexSPD();
305 :
306 : Float_t chi2TPCConstrainedVsGlobal =
307 0 : (vertex->GetStatus()) ? esdTrack->GetChi2TPCConstrainedVsGlobal(vertex) : (fEsdTrackVCutsChi2TPCconstrainedVsGlobal + 10.);
308 :
309 0 : if (chi2TPCConstrainedVsGlobal > fEsdTrackVCutsChi2TPCconstrainedVsGlobal)
310 0 : return kFALSE;
311 :
312 : Float_t cutDCAToVertexXYPtDep =
313 0 : fEsdTrackCutPtDCAOfs + fEsdTrackCutPtDCACoeff/((TMath::Abs(esdTrack->Pt()) > 0.0001) ? esdTrack->Pt() : 0.0001);
314 :
315 0 : if (TMath::Abs(impactPos[0]) >= cutDCAToVertexXYPtDep)
316 0 : return kFALSE;
317 :
318 0 : }
319 :
320 130 : return kTRUE;
321 282 : }
322 :
323 : Bool_t AliTRDonlineTrackMatching::ProcessEvent(AliESDEvent *esdEvent, Bool_t updateRef, Int_t label) {
324 :
325 : // performs track matching for all TRD online tracks of the ESD event
326 :
327 8 : UInt_t numTrdTracks = esdEvent->GetNumberOfTrdTracks();
328 8 : if (numTrdTracks <= 0)
329 0 : return kTRUE;
330 :
331 8 : if (!AliGeomManager::GetGeometry()){
332 0 : AliError("Geometry not available! Skipping TRD track matching.");
333 0 : return kFALSE;
334 : }
335 :
336 8 : if (!fTRDgeo){
337 4 : fTRDgeo = new AliTRDgeometry();
338 2 : }
339 :
340 : //
341 : // ESD track selection and sorting by TRD stack
342 : //
343 :
344 8 : UInt_t esdTracksByStack[fgkTrdStacks][fgkMaxEsdTracksPerStack];
345 8 : UInt_t esdTrackNumByStack[fgkTrdStacks];
346 8 : memset(esdTrackNumByStack, 0, fgkTrdStacks*sizeof(UInt_t));
347 :
348 8 : UInt_t numEsdTracks = esdEvent->GetNumberOfTracks();
349 : #ifdef TRD_TM_DEBUG
350 : UInt_t numEsdTracksAccepted = 0;
351 : #endif
352 8 : Short_t stack;
353 8 : UShort_t layers;
354 : AliESDtrack* esdTrack;
355 :
356 320 : for (UInt_t iEsdTrack = 0; iEsdTrack < numEsdTracks; ++iEsdTrack){
357 152 : esdTrack = esdEvent->GetTrack(iEsdTrack);
358 :
359 152 : if (!esdTrack){
360 0 : AliError("invalid ESD track!");
361 0 : continue;
362 : }
363 :
364 : // track filter here
365 152 : if (!AcceptTrack(esdTrack, esdEvent))
366 : continue;
367 : #ifdef TRD_TM_DEBUG
368 : else
369 : numEsdTracksAccepted++;
370 : #endif
371 :
372 : // assign ESD track to TRD stack
373 130 : if (StackToTrack(esdTrack, stack, layers, esdEvent->GetMagneticField())){
374 :
375 106 : if (stack < 0){
376 : #ifdef TRD_TM_DEBUG
377 : printf("#TRACKMATCHING - invalid stack for ESD track\n");
378 : #endif
379 : continue;
380 : }
381 :
382 : // register track in relevant stacks
383 106 : Int_t stacksForReg[9] = {-1, -1, -1, -1, -1, -1, -1, -1, -1};
384 106 : stacksForReg[0] = stack; // stack hit
385 106 : stacksForReg[1] = (stack + 5) % 90; // same stack in next supermodule
386 106 : stacksForReg[2] = (stack - 5); // same stack in previous supermodule
387 106 : if (stacksForReg[2] < 0)
388 16 : stacksForReg[2] += 90;
389 :
390 212 : switch(TrdLsiSi(stack)){
391 : case 0:
392 : // stack 0
393 4 : stacksForReg[3] = stack + 1; // next stack in same supermodule
394 4 : stacksForReg[4] = stacksForReg[1] + 1; // next stack in next supermodule
395 4 : stacksForReg[5] = stacksForReg[2] + 1; // next stack in previous supermodule
396 4 : break;
397 : case 1:
398 : case 2:
399 : case 3:
400 100 : stacksForReg[3] = stack + 1; // next stack in same supermodule
401 100 : stacksForReg[4] = stacksForReg[1] + 1; // next stack in next supermodule
402 100 : stacksForReg[5] = stacksForReg[2] + 1; // next stack in previous supermodule
403 100 : stacksForReg[6] = stack - 1; // previous stack in same supermodule
404 100 : stacksForReg[7] = stacksForReg[1] - 1; // previous stack in next supermodule
405 100 : stacksForReg[8] = stacksForReg[2] - 1; // previous stack in previous supermodule
406 100 : break;
407 : case 4:
408 2 : stacksForReg[3] = stack - 1; // previous stack in same supermodule
409 2 : stacksForReg[4] = stacksForReg[1] - 1; // previous stack in next supermodule
410 2 : stacksForReg[5] = stacksForReg[2] - 1; // previous stack in previous supermodule
411 2 : break;
412 : default:
413 : break;
414 : }
415 :
416 : #ifdef TRD_TM_DEBUG
417 : printf("#TRACKMATCHING - assigned ESD track %d to following TRD stacks:", iEsdTrack);
418 : #endif
419 :
420 : // register for stacks
421 2184 : for (UShort_t iReg = 0; iReg < 9; ++iReg){
422 942 : if (stacksForReg[iReg] < 0)
423 6 : break;
424 :
425 936 : if (stacksForReg[iReg] >= 90){
426 0 : AliError(Form("invalid stack for registration: %i", stacksForReg[iReg]));
427 0 : continue;
428 : }
429 :
430 936 : if (esdTrackNumByStack[stacksForReg[iReg]] < fgkMaxEsdTracksPerStack - 1)
431 936 : esdTracksByStack[stacksForReg[iReg]][esdTrackNumByStack[stacksForReg[iReg]]++] = iEsdTrack;
432 : #ifdef TRD_TM_DEBUG
433 : else
434 : printf("#TRACKMATCHING - maximum number (%d) of ESD tracks per stack reached for S%02d-%d (%d tracks total). Skipping track!\n",
435 : fgkMaxEsdTracksPerStack, TrdLsiSec(stacksForReg[iReg]), TrdLsiSi(stacksForReg[iReg]), numEsdTracks);
436 : printf(" S%02d-%d", TrdLsiSec(stacksForReg[iReg]), TrdLsiSi(stacksForReg[iReg]));
437 : #endif
438 : }
439 : #ifdef TRD_TM_DEBUG
440 : printf(" (ESD-ASSIGN)\n");
441 : #endif
442 :
443 : // if (esdTrackNumByStack[stack] >= fgkMaxEsdTracksPerStack){
444 : //#ifdef TRD_TM_DEBUG
445 : // printf("#TRACKMATCHING - maximum number (%d) of ESD tracks per stack reached for S%02d-%d (%d tracks total). Skipping track!\n",
446 : // fgkMaxEsdTracksPerStack, TrdLsiSec(stack), TrdLsiSi(stack), numEsdTracks);
447 : //#endif
448 : // continue;
449 : // }
450 : //
451 : // esdTracksByStack[stack][esdTrackNumByStack[stack]++] = iEsdTrack;
452 : //#ifdef TRD_TM_DEBUG
453 : // printf("#TRACKMATCHING - assigned ESD track %d to TRD stack S%02d-%d\n",
454 : // iEsdTrack, TrdLsiSec(stack), TrdLsiSi(stack));
455 : //#endif
456 106 : }
457 :
458 : } // loop over esd tracks
459 :
460 : #ifdef TRD_TM_DEBUG
461 : printf("#TRACKMATCHING - %d ESD tracks accepted, %d rejected\n",
462 : numEsdTracksAccepted, numEsdTracks - numEsdTracksAccepted);
463 : #endif
464 :
465 : //
466 : // search matching ESD track for each TRD online track
467 : //
468 : AliESDTrdTrack* trdTrack;
469 : Double_t trdPt;
470 : AliESDtrack* matchCandidate;
471 : AliESDtrack* matchTrack;
472 : Int_t matchEsdTrackIndexInStack;
473 : Double_t matchRating;
474 : Int_t matchCandidateCount;
475 8 : Double_t distY, distZ;
476 :
477 4654 : for (UInt_t iTrdTrack = 0; iTrdTrack < numTrdTracks; ++iTrdTrack){
478 :
479 2319 : trdTrack = esdEvent->GetTrdTrack(iTrdTrack);
480 2319 : if ((label != -1) &&
481 0 : (trdTrack->GetLabel() != label))
482 : continue;
483 :
484 6957 : if ((trdTrack->GetSector() < 0) || (trdTrack->GetSector() > 17) ||
485 4638 : (trdTrack->GetStack() < 0) || (trdTrack->GetStack() > 4))
486 : continue;
487 :
488 2319 : stack = TrdSecSiLsi(trdTrack->GetSector(), trdTrack->GetStack());
489 2319 : trdPt = (esdEvent->GetMagneticField() > 0.) ? (-1.*trdTrack->Pt()) : trdTrack->Pt();
490 : matchTrack = NULL;
491 : matchEsdTrackIndexInStack = -1;
492 : matchRating = 0.;
493 : matchCandidateCount = 0;
494 :
495 : #ifdef TRD_TM_DEBUG
496 : printf("#TRACKMATCHING - trying to match TRD online track %d in S%02d-%d\n",
497 : iTrdTrack, trdTrack->GetSector(), trdTrack->GetStack());
498 : #endif
499 :
500 : // loop over all esd tracks in the same stack and check distance
501 6678 : for (UInt_t iEsdTrack = 0; iEsdTrack < esdTrackNumByStack[stack]; ++iEsdTrack){
502 1020 : matchCandidate = esdEvent->GetTrack(esdTracksByStack[stack][iEsdTrack]);
503 :
504 1020 : if (EstimateTrackDistance(matchCandidate, trdTrack, esdEvent->GetMagneticField(), &distY, &distZ) == 0){
505 60 : Double_t rating = RateTrackMatch(distY, distZ, matchCandidate->GetSignedPt(), trdPt);
506 : #ifdef TRD_TM_DEBUG
507 : printf("#TRACKMATCHING S%02d-%d trd %d - esd %d dy: %.3f dz: %.3f r: %.3f pt e: %.2f t: %.2f match: %d\n",
508 : trdTrack->GetSector(), trdTrack->GetStack(), iTrdTrack, iEsdTrack,
509 : distY, distZ, rating, matchCandidate->GetSignedPt(), trdPt,
510 : (rating >= fMinMatchRating) ? 1 : 0);
511 : #endif
512 60 : if (rating > 0.){
513 : // possibly matching pair found
514 14 : matchCandidateCount++;
515 14 : if ((matchTrack == NULL) || (rating > matchRating)){
516 : // new best match
517 : matchTrack = matchCandidate;
518 : matchEsdTrackIndexInStack = iEsdTrack;
519 : matchRating = rating;
520 14 : }
521 : }
522 :
523 60 : } else {
524 : // estimation of distance failed
525 : #ifdef TRD_TM_DEBUG
526 : printf("TRACKMATCHING S%02d-%d trd %d - esd %d failed\n",
527 : trdTrack->GetSector(), trdTrack->GetStack(), iTrdTrack, iEsdTrack);
528 : #endif
529 : }
530 : } // loop over esd tracks in same stack
531 :
532 2319 : if (fHistMatchRating){
533 0 : fHistMatchRating->Fill(matchRating);
534 0 : }
535 :
536 2333 : if ((matchTrack) && (matchRating >= fMinMatchRating)){
537 39 : AliDebug(1, Form("S%02d-%d trd %d - esd %d match! pt: %.2f %.2f",
538 : trdTrack->GetSector(), trdTrack->GetStack(), iTrdTrack, matchEsdTrackIndexInStack,
539 : trdPt, matchTrack->GetSignedPt()));
540 : #ifdef TRD_TM_DEBUG
541 : printf("#TRACKMATCHING S%02d-%d trd %d - esd %d match! pt: %.2f %.2f\n",
542 : trdTrack->GetSector(), trdTrack->GetStack(), iTrdTrack, matchEsdTrackIndexInStack,
543 : trdPt, matchTrack->GetSignedPt());
544 : #endif
545 13 : if (updateRef)
546 13 : trdTrack->SetTrackMatchReference(matchTrack);
547 : } else {
548 2306 : if (updateRef)
549 2306 : trdTrack->SetTrackMatchReference(NULL);
550 : }
551 :
552 : } // loop over TRD online tracks
553 :
554 : return kTRUE;
555 16 : }
556 :
557 : Bool_t AliTRDonlineTrackMatching::TrackPlaneIntersect(AliExternalTrackParam *trk, Double_t pnt[3], Double_t norm[3], Double_t mag){
558 :
559 : // calculates the intersection point of a track param and a plane defined by point pnt and normal vector norm
560 :
561 : UInt_t its = 0;
562 : Double_t r = 290.;
563 : Double_t step = 10.;
564 : Int_t flag = 0;
565 : Double_t dist = 0, dist_prev = 0;
566 584 : Double_t x[3] = {0., 0., 0.};
567 :
568 292 : dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) *norm[1] + (x[2] - pnt[2]) * norm[2];
569 :
570 5544 : while(TMath::Abs(dist) > 0.1) {
571 :
572 5001 : trk->GetXYZAt(r, mag, x);
573 :
574 5001 : if ((x[0] * x[0] + x[1] * x[1]) < 100.) // extrapolation to radius failed
575 0 : return kFALSE;
576 :
577 : //distance between current track position and plane
578 5001 : dist_prev = TMath::Abs(dist);
579 5001 : dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) * norm[1];
580 9710 : if ((flag) && (TMath::Abs(dist) > dist_prev)){
581 1080 : step /= -2.;
582 1080 : }
583 : flag=1;
584 5001 : r += step;
585 5001 : its++;
586 9961 : if ((r > 380.) || (r < 100.) || (its > 100) || (TMath::Abs(step) < 0.00001)){
587 : break;
588 : }
589 : }
590 2336 : for (Int_t i=0; i<3; i++)
591 876 : pnt[i] = x[i];
592 :
593 292 : return kTRUE;
594 292 : }
595 :
596 : Int_t AliTRDonlineTrackMatching::EstimateTrackDistance(AliESDtrack *esd_track, AliESDTrdTrack* gtu_track, Double_t mag, Double_t *ydist, Double_t *zdist){
597 :
598 : // returns an estimate for the spatial distance between TPC offline track and GTU online track
599 :
600 2040 : if ((!esd_track) || (!gtu_track))
601 0 : return -3;
602 :
603 : // AssertTRDGeometry();
604 1020 : if (!fTRDgeo)
605 0 : fTRDgeo = new AliTRDgeometry();
606 :
607 : Float_t diff_y = 0;
608 : Float_t diff_z = 0;
609 : Int_t nLayers = 0;
610 1020 : Double_t xtrkl[3];
611 1020 : Double_t ptrkl[3];
612 1020 : Double_t ptrkl2[3];
613 : UInt_t trklDet;
614 : UShort_t trklLayer;
615 : UInt_t stack_gtu;
616 : UShort_t stackInSector;
617 :
618 15300 : for (UShort_t iLayer = 0; iLayer < 6; iLayer++){
619 6120 : AliESDTrdTracklet* trkl = gtu_track->GetTracklet(iLayer);
620 6120 : if (trkl){
621 292 : trklDet = trkl->GetDetector();
622 292 : trklLayer = TrdDetLyr(trklDet);
623 292 : stack_gtu = TrdDetLsi(trklDet);
624 292 : stackInSector = TrdDetSi(trklDet);
625 :
626 : // local coordinates of the outer end point of the tracklet
627 292 : xtrkl[0] = AliTRDgeometry::AnodePos();
628 292 : xtrkl[1] = trkl->GetLocalY();
629 :
630 584 : if(stackInSector == 2){ // corrected version by Felix Muecke
631 800 : xtrkl[2] = fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(trkl->GetBinZ()) -
632 508 : (fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowSize(trkl->GetBinZ()))/2. -
633 508 : fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(6);
634 216 : } else {
635 76 : xtrkl[2] = fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(trkl->GetBinZ()) -
636 76 : (fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowSize(trkl->GetBinZ()))/2. -
637 76 : fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(8);
638 : }
639 :
640 : // old draft version
641 : // xtrkl[2] = fTRDgeo->GetPadPlane(trklLayer, (trklDet/6) % 5)->GetRowPos(trkl->GetBinZ()) -
642 : // fTRDgeo->GetPadPlane(trklLayer, (trklDet/6) % 5)->GetRowSize(trkl->GetBinZ()) -
643 : // fTRDgeo->GetPadPlane(trklLayer, (trklDet/6) % 5)->GetRowPos(8);
644 :
645 : // transform to global coordinates
646 292 : TGeoHMatrix *matrix = fTRDgeo->GetClusterMatrix(trklDet);
647 292 : if (!matrix){
648 0 : if ((stack_gtu != 13*5+2) && (stack_gtu != 14*5+2) && (stack_gtu != 15*5+2))
649 0 : AliDebug(1, Form("invalid TRD cluster matrix in EstimateTrackDistance for detector %i", trklDet));
650 0 : return -5;
651 : }
652 292 : matrix->LocalToMaster(xtrkl, ptrkl);
653 292 : fTRDgeo->RotateBack(gtu_track->GetSector() * 30, ptrkl, ptrkl2); // ptrkl2 now contains the global position of the outer end point of the tracklet
654 :
655 : // calculate parameterization of plane representing the tracklets layer
656 292 : Double_t layer_zero_local[3] = {0., 0., 0.};
657 292 : Double_t layer_zero_global[3], layer_zero_global2[3];
658 :
659 292 : matrix->LocalToMaster(layer_zero_local, layer_zero_global);
660 292 : fTRDgeo->RotateBack(trklDet, layer_zero_global, layer_zero_global2); // layer_zero_global2 points to chamber origin in global coords
661 :
662 292 : Double_t layer_ref_local[3] = {AliTRDgeometry::AnodePos(), 0., 0.};
663 292 : Double_t layer_ref_global[3], layer_ref_global2[3];
664 :
665 292 : matrix->LocalToMaster(layer_ref_local, layer_ref_global);
666 292 : fTRDgeo->RotateBack(trklDet, layer_ref_global, layer_ref_global2); // layer_ref_global2 points to center anode pos within plane in global coords
667 :
668 292 : Double_t n0[3] = {layer_ref_global2[0]-layer_zero_global2[0],
669 292 : layer_ref_global2[1]-layer_zero_global2[1],
670 292 : layer_ref_global2[2]-layer_zero_global2[2]};
671 :
672 292 : Double_t n_len = TMath::Sqrt(n0[0]*n0[0] + n0[1]*n0[1] + n0[2]*n0[2]);
673 292 : if (n_len == 0.){ // This should never happen
674 0 : AliError("divison by zero in estimate_track_distance!");
675 : n_len = 1.;
676 0 : }
677 292 : Double_t n[3] = {n0[0]/n_len, n0[1]/n_len, n0[2]/n_len}; // normal vector of plane
678 :
679 292 : const AliExternalTrackParam *trackParam = esd_track->GetOuterParam();
680 292 : if (!trackParam) {
681 0 : trackParam = esd_track->GetInnerParam();
682 0 : if (!trackParam)
683 0 : trackParam = esd_track;
684 : }
685 :
686 292 : AliExternalTrackParam *outerTPC = new AliExternalTrackParam(*trackParam);
687 292 : Bool_t isects = TrackPlaneIntersect(outerTPC, layer_ref_global2, n, mag); // find intersection point between track and TRD layer
688 584 : delete outerTPC;
689 : outerTPC = NULL;
690 :
691 292 : if (isects == kFALSE){ // extrapolation fails, because track never reaches the TRD radius
692 0 : return -1;
693 : }
694 :
695 292 : Double_t m[2] = {ptrkl2[0] - layer_ref_global2[0], ptrkl2[1] - layer_ref_global2[1]};
696 292 : Double_t len_m = TMath::Sqrt(m[0]*m[0] + m[1]*m[1]);
697 292 : diff_y += len_m;
698 292 : diff_z += TMath::Abs(ptrkl2[2] - layer_ref_global2[2]);
699 292 : nLayers++;
700 584 : }
701 6120 : }
702 :
703 1020 : if (nLayers > 0){
704 60 : *ydist = diff_y / nLayers;
705 60 : *zdist = diff_z / nLayers;
706 60 : return 0;
707 : }
708 : else
709 960 : return -4;
710 2040 : }
711 :
712 : Double_t AliTRDonlineTrackMatching::PtDiffRel(Double_t refPt, Double_t gtuPt){
713 :
714 : // return relative pt difference
715 :
716 28 : if (TMath::Abs(refPt) > 0.000001){
717 14 : return (gtuPt - refPt) / refPt;
718 : } else
719 0 : return 0.;
720 14 : }
721 :
722 :
723 : Double_t AliTRDonlineTrackMatching::RateTrackMatch(Double_t distY, Double_t distZ, Double_t rpt, Double_t gpt){
724 :
725 : // returns a match rating derived from Y and Z distance as well as pt difference
726 :
727 : // maximum limits for spatial distance
728 120 : if ((distY > 5.) || (distZ > 20.))
729 46 : return 0.;
730 :
731 : // same pt sign required
732 14 : if ((rpt * gpt) < 0.)
733 0 : return 0.;
734 :
735 14 : Double_t rating_distY = -0.1 * distY + 1.;
736 14 : Double_t rating_distZ = -0.025 * distZ + 1.;
737 14 : Double_t rating_ptDiff = 1. - TMath::Abs(PtDiffRel(rpt, gpt));
738 :
739 14 : if (rating_ptDiff < 0.)
740 0 : rating_ptDiff = 0.2;
741 :
742 14 : Double_t total = rating_distY * rating_distZ * rating_ptDiff;
743 :
744 : #ifdef TRD_TM_DEBUG
745 : if (total > 1.){
746 : printf("<ERROR> track match rating exceeds limit of 1.0: %.3f", total);
747 : }
748 : #endif
749 :
750 : return total;
751 60 : }
752 :
753 :
754 : void AliTRDonlineTrackMatching::SetEsdTrackDefaultCuts(const char* cutIdent) {
755 :
756 48 : if (strcmp(cutIdent, "strict") == 0){
757 :
758 : #ifdef TRD_TM_DEBUG
759 : printf("AliTRDonlineTrackMatching -- default track cuts selected");
760 : #endif
761 :
762 0 : fEsdTrackCutMinimal = kFALSE;
763 0 : fEsdTrackCutPrim = kFALSE;
764 :
765 0 : fEsdTrackCutMinTPCrows = 70;
766 0 : fEsdTrackCutRequireTPCrefit = kTRUE;
767 0 : fEsdTrackCutMinRatioRowsFindableClusters = 0.8;
768 0 : fEsdTrackCutMaxChi2TPCclusters = 4.;
769 0 : fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 36.;
770 :
771 0 : fEsdTrackCutRequireITSrefit = kFALSE;
772 0 : fEsdTrackCutMaxChi2ITSclusters = 36.;
773 :
774 0 : fEsdTrackCutMaxDCAtoVertexXY = 1000.;
775 0 : fEsdTrackCutMaxDCAtoVertexZ = 2.;
776 0 : fEsdTrackCutsITSlayerMask = 0x0;
777 :
778 0 : fEsdTrackCutPtDCAOfs = 0.0105;
779 0 : fEsdTrackCutPtDCACoeff = 0.0350;
780 24 : } else if (strcmp(cutIdent, "minimal") == 0){
781 :
782 : #ifdef TRD_TM_DEBUG
783 : printf("AliTRDonlineTrackMatching -- minimal track cuts selected\n");
784 : #endif
785 :
786 24 : fEsdTrackCutMinimal = kFALSE;
787 24 : fEsdTrackCutPrim = kFALSE;
788 :
789 24 : fEsdTrackCutMinTPCrows = 70;
790 24 : fEsdTrackCutRequireTPCrefit = kTRUE;
791 24 : fEsdTrackCutMinRatioRowsFindableClusters = 0.;
792 24 : fEsdTrackCutMaxChi2TPCclusters = 100.;
793 24 : fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 1000.;
794 :
795 24 : fEsdTrackCutRequireITSrefit = kFALSE;
796 24 : fEsdTrackCutMaxChi2ITSclusters = 0.;
797 :
798 24 : fEsdTrackCutMaxDCAtoVertexXY = 1000.;
799 24 : fEsdTrackCutMaxDCAtoVertexZ = 1000.;
800 24 : fEsdTrackCutsITSlayerMask = 0x0;
801 24 : } else
802 0 : AliErrorClass("invalid cut set");
803 :
804 24 : }
|
