- add protection against invalid sector/stack information in on-line track matching
[u/mrichter/AliRoot.git] / TRD / AliTRDonlineTrackMatching.cxx
CommitLineData
30bd9a5a 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
34const Float_t AliTRDonlineTrackMatching::fgkSaveInnerRadius = 290.5;
35const Float_t AliTRDonlineTrackMatching::fgkSaveOuterRadius = 364.5;
36
37Float_t AliTRDonlineTrackMatching::fEsdTrackCutMinTPCrows = 0.;
38Float_t AliTRDonlineTrackMatching::fEsdTrackCutMinRatioRowsFindableClusters = 0.;
39Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxChi2TPCclusters = 0.;
40Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxChi2ITSclusters = 0.;
41Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxDCAtoVertexXY = 0.;
42Float_t AliTRDonlineTrackMatching::fEsdTrackCutMaxDCAtoVertexZ = 0.;
43UShort_t AliTRDonlineTrackMatching::fEsdTrackCutsITSlayerMask = 0; // similar to 2011 default cut: 0x3
44Float_t AliTRDonlineTrackMatching::fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 0.;
45Float_t AliTRDonlineTrackMatching::fEsdTrackCutPtDCAOfs = 0.;
46Float_t AliTRDonlineTrackMatching::fEsdTrackCutPtDCACoeff = 0.;
47Bool_t AliTRDonlineTrackMatching::fEsdTrackCutMinimal = kFALSE;
48Bool_t AliTRDonlineTrackMatching::fEsdTrackCutRequireTPCrefit = kTRUE;
49Bool_t AliTRDonlineTrackMatching::fEsdTrackCutRequireITSrefit = kFALSE;
50Bool_t AliTRDonlineTrackMatching::fEsdTrackCutPrim = kFALSE;
51
52AliTRDonlineTrackMatching::AliTRDonlineTrackMatching() :
633f3ec4 53 TObject(),
30bd9a5a 54 fTRDgeo(NULL),
55 fMinMatchRating(0.25),
56 fHistMatchRating(NULL)
57{
58 // default ctor
59 SetEsdTrackDefaultCuts("minimal");
60}
61
62AliTRDonlineTrackMatching::AliTRDonlineTrackMatching(const AliTRDonlineTrackMatching &c) :
633f3ec4 63 TObject(c),
30bd9a5a 64 fTRDgeo(c.fTRDgeo),
65 fMinMatchRating(c.fMinMatchRating),
66 fHistMatchRating(c.fHistMatchRating)
67{
68 // copy ctor
69}
70
71AliTRDonlineTrackMatching::~AliTRDonlineTrackMatching() {
72
73 // dtor
74
75 delete fTRDgeo;
76 fTRDgeo = NULL;
77}
78
79Short_t AliTRDonlineTrackMatching::EstimateSector(const Double_t globalCoords[3]) {
80
81 // estimates sector by phi angle in x-y plane
82
83 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 return -1;
86 } else {
87 Double_t ang = TMath::ATan2(globalCoords[1], globalCoords[0]);
88 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 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 return 17 - TMath::FloorNint(TMath::Abs(ang)/(20./180.*TMath::Pi()));
102 }
103
104 }
105}
106
6403ed35 107Short_t AliTRDonlineTrackMatching::EstimateLayer(Double_t radius) {
30bd9a5a 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 for (UShort_t i = 0; i < 7; ++i){
114 if (radius < rBoundaries[i])
115 return rLayers[i];
116 }
117 return -2; // radius larger than outmost layer
118}
119
120Short_t AliTRDonlineTrackMatching::EstimateLocalStack(const Double_t globalCoords[3]) {
121
122 // determines stack within sector by z position
123
124 Double_t absZ = TMath::Abs(globalCoords[2]);
125 Short_t signZ = (globalCoords[2] > 0.) ? 1 : -1;
126 Double_t r = TMath::Sqrt(globalCoords[0]*globalCoords[0] + globalCoords[1]*globalCoords[1]);
127 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 if (layer < 0)
135 return -1;
136
137 Double_t innerStackHalfLength = AliTRDgeometry::GetChamberLength(0, 2) / 2.; // same for all layers
138 if (absZ < innerStackHalfLength)
139 return 2;
140
141 Double_t outerStackLength = AliTRDgeometry::GetChamberLength(layer, 1);
142
143 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 if (absZ > 2.05*outerStackLength)
151 return (signZ > 0) ? -2 : -1; // outside supermodule in z direction
152
153 if (absZ < outerStackLength)
154 return (signZ > 0) ? 1 : 3;
155 else
156 return (signZ > 0) ? 0 : 4;
157
158}
159
160Short_t AliTRDonlineTrackMatching::EstimateStack(const Double_t globalCoords[3]) {
161
162 // returns the closest TRD stack to a 3D position in global coordinates
163
164 Short_t sec = EstimateSector(globalCoords);
165 Short_t st = EstimateLocalStack(globalCoords);
166#ifdef TRD_TM_DEBUG
167 printf("EstimateStack sec %d st %d\n", sec, st);
168#endif
169 if ((sec < 0) || (st < 0))
170 return -1;
171 else
172 return 5*sec + st;
173}
174
6403ed35 175Bool_t AliTRDonlineTrackMatching::StackToTrack(const AliExternalTrackParam *track, Short_t &stack, UShort_t &layersWithTracklet, Double_t magFieldinKiloGauss){
30bd9a5a 176
177 // returns stack to track param
178
179 stack = -1;
180 layersWithTracklet = 0;
181
182 UInt_t stackHits[fgkTrdStacks];
5d889c14 183 Double_t x[3] = { 0. };
30bd9a5a 184 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 while (r < fgkSaveOuterRadius){
5d889c14 192 if (track->GetXYZAt(r, magFieldinKiloGauss, x)) {
193 stack = EstimateStack(x);
194 if (stack >= 0){
195 stackHits[stack]++;
196 if (stackHits[stack] > 16) // experimental
197 break;
30bd9a5a 198#ifdef TRD_TM_DEBUG
5d889c14 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);
30bd9a5a 201#endif
5d889c14 202 }
30bd9a5a 203 }
204 r += 1.;
205 }
206
207 // find stack with most hits
208 UInt_t bestHits = 0;
209 for (UShort_t iStack = 0; iStack < fgkTrdStacks; ++iStack){
210 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 if (stackHits[iStack] > bestHits){
216 bestHits = stackHits[iStack];
217 stack = iStack;
218 }
219 }
220
221 if (stack >= 0){
222#ifdef TRD_TM_DEBUG
223 printf("best stack: S%02d-%d\n", TrdLsiSec(stack), TrdLsiSi(stack));
224#endif
225 return kTRUE;
226 }
227
228 return kFALSE;
229}
230
6403ed35 231Bool_t AliTRDonlineTrackMatching::StackToTrack(const AliESDtrack* track, Short_t &stack, UShort_t &layersWithTracklet, Double_t magFieldinKiloGauss){
30bd9a5a 232
233 // returns stack to ESD track
234
235 if (track->GetOuterParam())
236 return StackToTrack(track->GetOuterParam(), stack, layersWithTracklet, magFieldinKiloGauss);
237 else if (track->GetInnerParam())
238 return StackToTrack(track->GetInnerParam(), stack, layersWithTracklet, magFieldinKiloGauss);
239 else
240 return StackToTrack(track, stack, layersWithTracklet, magFieldinKiloGauss);
241}
242
243Bool_t AliTRDonlineTrackMatching::AcceptTrack(const AliESDtrack* esdTrack, const AliESDEvent* esdEvent){
244
245 // returns result ESD track cuts
246
ca892663 247 if (!esdTrack)
248 return kFALSE;
249
250 UInt_t status = esdTrack->GetStatus();
30bd9a5a 251
252 if (fEsdTrackCutMinimal){
253 return ((status & AliESDtrack::kTPCout) > 0);
254 }
255
256 // require TPC fit
257 if ((fEsdTrackCutRequireTPCrefit) && (!(status & AliESDtrack::kTPCrefit)))
258 return kFALSE;
259
260 // require ITS re-fit
261 if ((fEsdTrackCutRequireITSrefit) && (!(status & AliESDtrack::kITSrefit)))
262 return kFALSE;
263
264 // TPC requirements
265 Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
266 Float_t ratioCrossedRowsOverFindableClustersTPC =
267 (esdTrack->GetTPCNclsF() > 0) ? (nCrossedRowsTPC / esdTrack->GetTPCNclsF()) : 1.0;
268 Float_t chi2PerClusterTPC =
269 (esdTrack->GetTPCclusters(0) > 0) ? (esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCclusters(0))) : 100.;
270
271 if (
272 (nCrossedRowsTPC < fEsdTrackCutMinTPCrows) ||
273 (ratioCrossedRowsOverFindableClustersTPC < fEsdTrackCutMinRatioRowsFindableClusters) ||
274 (chi2PerClusterTPC > fEsdTrackCutMaxChi2TPCclusters)
275 )
276 return kFALSE;
277
278 // ITS requirements
279 Float_t chi2PerClusterITS = (esdTrack->GetITSclusters(0) > 0) ? esdTrack->GetITSchi2()/Float_t(esdTrack->GetITSclusters(0)) : 1000.;
280 UShort_t clustersInAnyITSlayer = kFALSE;
281 for (UShort_t layer = 0; layer < 6; ++layer)
282 clustersInAnyITSlayer += (esdTrack->HasPointOnITSLayer(layer) & ((fEsdTrackCutsITSlayerMask >> layer) & 1));
283
284 if ((fEsdTrackCutsITSlayerMask != 0) &&
285 ((clustersInAnyITSlayer == 0) || (chi2PerClusterITS >= fEsdTrackCutMaxChi2ITSclusters))
286 )
287 return kFALSE;
288
289 // geometric requirements
290 Float_t impactPos[2], impactCov[3];
291 esdTrack->GetImpactParameters(impactPos, impactCov);
292
293 if (TMath::Abs(impactPos[0]) > fEsdTrackCutMaxDCAtoVertexXY)
294 return kFALSE;
295
296 if (TMath::Abs(impactPos[1]) > fEsdTrackCutMaxDCAtoVertexZ)
297 return kFALSE;
298
299 if (fEsdTrackCutPrim){
300 // additional requirements for primary tracks
301
302 const AliESDVertex* vertex = esdEvent->GetPrimaryVertexTracks();
303 if ((!vertex) || (!vertex->GetStatus()))
304 vertex = esdEvent->GetPrimaryVertexSPD();
305
306 Float_t chi2TPCConstrainedVsGlobal =
307 (vertex->GetStatus()) ? esdTrack->GetChi2TPCConstrainedVsGlobal(vertex) : (fEsdTrackVCutsChi2TPCconstrainedVsGlobal + 10.);
308
309 if (chi2TPCConstrainedVsGlobal > fEsdTrackVCutsChi2TPCconstrainedVsGlobal)
310 return kFALSE;
311
312 Float_t cutDCAToVertexXYPtDep =
313 fEsdTrackCutPtDCAOfs + fEsdTrackCutPtDCACoeff/((TMath::Abs(esdTrack->Pt()) > 0.0001) ? esdTrack->Pt() : 0.0001);
314
315 if (TMath::Abs(impactPos[0]) >= cutDCAToVertexXYPtDep)
316 return kFALSE;
317
318 }
319
320 return kTRUE;
321}
322
34d48b8f 323Bool_t AliTRDonlineTrackMatching::ProcessEvent(AliESDEvent *esdEvent, Bool_t updateRef, Int_t label) {
30bd9a5a 324
325 // performs track matching for all TRD online tracks of the ESD event
326
327 UInt_t numTrdTracks = esdEvent->GetNumberOfTrdTracks();
328 if (numTrdTracks <= 0)
329 return kTRUE;
330
331 if (!AliGeomManager::GetGeometry()){
633f3ec4 332 AliError("Geometry not available! Skipping TRD track matching.");
30bd9a5a 333 return kFALSE;
334 }
335
336 if (!fTRDgeo){
337 fTRDgeo = new AliTRDgeometry();
338 }
339
340 //
341 // ESD track selection and sorting by TRD stack
342 //
343
344 UInt_t esdTracksByStack[fgkTrdStacks][fgkMaxEsdTracksPerStack];
345 UInt_t esdTrackNumByStack[fgkTrdStacks];
346 memset(esdTrackNumByStack, 0, fgkTrdStacks*sizeof(UInt_t));
347
348 UInt_t numEsdTracks = esdEvent->GetNumberOfTracks();
349#ifdef TRD_TM_DEBUG
350 UInt_t numEsdTracksAccepted = 0;
351#endif
352 Short_t stack;
353 UShort_t layers;
354 AliESDtrack* esdTrack;
355
356 for (UInt_t iEsdTrack = 0; iEsdTrack < numEsdTracks; ++iEsdTrack){
357 esdTrack = esdEvent->GetTrack(iEsdTrack);
358
359 if (!esdTrack){
633f3ec4 360 AliError("invalid ESD track!");
30bd9a5a 361 continue;
362 }
363
364 // track filter here
365 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 if (StackToTrack(esdTrack, stack, layers, esdEvent->GetMagneticField())){
374
375 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 Int_t stacksForReg[9] = {-1, -1, -1, -1, -1, -1, -1, -1, -1};
384 stacksForReg[0] = stack; // stack hit
385 stacksForReg[1] = (stack + 5) % 90; // same stack in next supermodule
386 stacksForReg[2] = (stack - 5); // same stack in previous supermodule
387 if (stacksForReg[2] < 0)
388 stacksForReg[2] += 90;
389
390 switch(TrdLsiSi(stack)){
391 case 0:
392 // stack 0
393 stacksForReg[3] = stack + 1; // next stack in same supermodule
394 stacksForReg[4] = stacksForReg[1] + 1; // next stack in next supermodule
395 stacksForReg[5] = stacksForReg[2] + 1; // next stack in previous supermodule
396 break;
397 case 1:
398 case 2:
399 case 3:
400 stacksForReg[3] = stack + 1; // next stack in same supermodule
401 stacksForReg[4] = stacksForReg[1] + 1; // next stack in next supermodule
402 stacksForReg[5] = stacksForReg[2] + 1; // next stack in previous supermodule
403 stacksForReg[6] = stack - 1; // previous stack in same supermodule
404 stacksForReg[7] = stacksForReg[1] - 1; // previous stack in next supermodule
405 stacksForReg[8] = stacksForReg[2] - 1; // previous stack in previous supermodule
406 break;
407 case 4:
408 stacksForReg[3] = stack - 1; // previous stack in same supermodule
409 stacksForReg[4] = stacksForReg[1] - 1; // previous stack in next supermodule
410 stacksForReg[5] = stacksForReg[2] - 1; // previous stack in previous supermodule
411 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 for (UShort_t iReg = 0; iReg < 9; ++iReg){
422 if (stacksForReg[iReg] < 0)
423 break;
424
425 if (stacksForReg[iReg] >= 90){
633f3ec4 426 AliError(Form("invalid stack for registration: %i", stacksForReg[iReg]));
30bd9a5a 427 continue;
428 }
429
430 if (esdTrackNumByStack[stacksForReg[iReg]] < fgkMaxEsdTracksPerStack - 1)
431 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 }
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 Double_t distY, distZ;
476
477 for (UInt_t iTrdTrack = 0; iTrdTrack < numTrdTracks; ++iTrdTrack){
478
479 trdTrack = esdEvent->GetTrdTrack(iTrdTrack);
34d48b8f 480 if ((label != -1) &&
481 (trdTrack->GetLabel() != label))
482 continue;
483
7bee2115 484 if ((trdTrack->GetSector() < 0) || (trdTrack->GetSector() > 17) ||
485 (trdTrack->GetStack() < 0) || (trdTrack->GetStack() > 4))
486 continue;
487
30bd9a5a 488 stack = TrdSecSiLsi(trdTrack->GetSector(), trdTrack->GetStack());
489 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 for (UInt_t iEsdTrack = 0; iEsdTrack < esdTrackNumByStack[stack]; ++iEsdTrack){
502 matchCandidate = esdEvent->GetTrack(esdTracksByStack[stack][iEsdTrack]);
503
504 if (EstimateTrackDistance(matchCandidate, trdTrack, esdEvent->GetMagneticField(), &distY, &distZ) == 0){
505 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 if (rating > 0.){
513 // possibly matching pair found
514 matchCandidateCount++;
515 if ((matchTrack == NULL) || (rating > matchRating)){
516 // new best match
517 matchTrack = matchCandidate;
518 matchEsdTrackIndexInStack = iEsdTrack;
519 matchRating = rating;
520 }
521 }
522
523 } 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 if (fHistMatchRating){
533 fHistMatchRating->Fill(matchRating);
534 }
535
536 if ((matchTrack) && (matchRating >= fMinMatchRating)){
1d94ad26 537 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()));
30bd9a5a 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
34d48b8f 545 if (updateRef)
546 trdTrack->SetTrackMatchReference(matchTrack);
547 } else {
548 if (updateRef)
549 trdTrack->SetTrackMatchReference(NULL);
550 }
30bd9a5a 551
552 } // loop over TRD online tracks
553
554 return kTRUE;
555}
556
557Bool_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 Double_t x[3] = {0., 0., 0.};
567
568 dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) *norm[1] + (x[2] - pnt[2]) * norm[2];
569
570 while(TMath::Abs(dist) > 0.1) {
571
572 trk->GetXYZAt(r, mag, x);
573
574 if ((x[0] * x[0] + x[1] * x[1]) < 100.) // extrapolation to radius failed
575 return kFALSE;
576
577 //distance between current track position and plane
578 dist_prev = TMath::Abs(dist);
579 dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) * norm[1];
580 if ((flag) && (TMath::Abs(dist) > dist_prev)){
581 step /= -2.;
582 }
583 flag=1;
584 r += step;
585 its++;
586 if ((r > 380.) || (r < 100.) || (its > 100) || (TMath::Abs(step) < 0.00001)){
587 break;
588 }
589 }
590 for (Int_t i=0; i<3; i++)
591 pnt[i] = x[i];
592
593 return kTRUE;
594}
595
596Int_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 if ((!esd_track) || (!gtu_track))
601 return -3;
602
603 // AssertTRDGeometry();
604 if (!fTRDgeo)
605 fTRDgeo = new AliTRDgeometry();
606
607 Float_t diff_y = 0;
608 Float_t diff_z = 0;
609 Int_t nLayers = 0;
610 Double_t xtrkl[3];
611 Double_t ptrkl[3];
612 Double_t ptrkl2[3];
613 UInt_t trklDet;
614 UShort_t trklLayer;
615 UInt_t stack_gtu;
616 UShort_t stackInSector;
617
618 for (UShort_t iLayer = 0; iLayer < 6; iLayer++){
619 AliESDTrdTracklet* trkl = gtu_track->GetTracklet(iLayer);
620 if (trkl){
621 trklDet = trkl->GetDetector();
622 trklLayer = TrdDetLyr(trklDet);
623 stack_gtu = TrdDetLsi(trklDet);
624 stackInSector = TrdDetSi(trklDet);
625
626 // local coordinates of the outer end point of the tracklet
627 xtrkl[0] = AliTRDgeometry::AnodePos();
628 xtrkl[1] = trkl->GetLocalY();
629
630 if(stackInSector == 2){ // corrected version by Felix Muecke
631 xtrkl[2] = fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(trkl->GetBinZ()) -
632 (fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowSize(trkl->GetBinZ()))/2. -
633 fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(6);
634 } else {
635 xtrkl[2] = fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowPos(trkl->GetBinZ()) -
636 (fTRDgeo->GetPadPlane(trklLayer, stackInSector)->GetRowSize(trkl->GetBinZ()))/2. -
637 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 TGeoHMatrix *matrix = fTRDgeo->GetClusterMatrix(trklDet);
647 if (!matrix){
633f3ec4 648 if ((stack_gtu != 13*5+2) && (stack_gtu != 14*5+2) && (stack_gtu != 15*5+2))
aa67162b 649 AliDebug(1, Form("invalid TRD cluster matrix in EstimateTrackDistance for detector %i", trklDet));
30bd9a5a 650 return -5;
651 }
652 matrix->LocalToMaster(xtrkl, ptrkl);
653 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 Double_t layer_zero_local[3] = {0., 0., 0.};
657 Double_t layer_zero_global[3], layer_zero_global2[3];
658
659 matrix->LocalToMaster(layer_zero_local, layer_zero_global);
660 fTRDgeo->RotateBack(trklDet, layer_zero_global, layer_zero_global2); // layer_zero_global2 points to chamber origin in global coords
661
662 Double_t layer_ref_local[3] = {AliTRDgeometry::AnodePos(), 0., 0.};
663 Double_t layer_ref_global[3], layer_ref_global2[3];
664
665 matrix->LocalToMaster(layer_ref_local, layer_ref_global);
666 fTRDgeo->RotateBack(trklDet, layer_ref_global, layer_ref_global2); // layer_ref_global2 points to center anode pos within plane in global coords
667
668 Double_t n0[3] = {layer_ref_global2[0]-layer_zero_global2[0],
669 layer_ref_global2[1]-layer_zero_global2[1],
670 layer_ref_global2[2]-layer_zero_global2[2]};
671
672 Double_t n_len = TMath::Sqrt(n0[0]*n0[0] + n0[1]*n0[1] + n0[2]*n0[2]);
673 if (n_len == 0.){ // This should never happen
633f3ec4 674 AliError("divison by zero in estimate_track_distance!");
30bd9a5a 675 n_len = 1.;
676 }
677 Double_t n[3] = {n0[0]/n_len, n0[1]/n_len, n0[2]/n_len}; // normal vector of plane
678
c433ad04 679 const AliExternalTrackParam *trackParam = esd_track->GetOuterParam();
680 if (!trackParam) {
681 trackParam = esd_track->GetInnerParam();
682 if (!trackParam)
683 trackParam = esd_track;
684 }
685
686 AliExternalTrackParam *outerTPC = new AliExternalTrackParam(*trackParam);
30bd9a5a 687 Bool_t isects = TrackPlaneIntersect(outerTPC, layer_ref_global2, n, mag); // find intersection point between track and TRD layer
688 delete outerTPC;
689 outerTPC = NULL;
690
691 if (isects == kFALSE){ // extrapolation fails, because track never reaches the TRD radius
692 return -1;
693 }
694
695 Double_t m[2] = {ptrkl2[0] - layer_ref_global2[0], ptrkl2[1] - layer_ref_global2[1]};
696 Double_t len_m = TMath::Sqrt(m[0]*m[0] + m[1]*m[1]);
697 diff_y += len_m;
698 diff_z += TMath::Abs(ptrkl2[2] - layer_ref_global2[2]);
699 nLayers++;
700 }
701 }
702
703 if (nLayers > 0){
704 *ydist = diff_y / nLayers;
705 *zdist = diff_z / nLayers;
706 return 0;
707 }
708 else
709 return -4;
710}
711
6403ed35 712Double_t AliTRDonlineTrackMatching::PtDiffRel(Double_t refPt, Double_t gtuPt){
30bd9a5a 713
714 // return relative pt difference
715
716 if (TMath::Abs(refPt) > 0.000001){
717 return (gtuPt - refPt) / refPt;
718 } else
719 return 0.;
720}
721
722
6403ed35 723Double_t AliTRDonlineTrackMatching::RateTrackMatch(Double_t distY, Double_t distZ, Double_t rpt, Double_t gpt){
30bd9a5a 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 if ((distY > 5.) || (distZ > 20.))
729 return 0.;
730
731 // same pt sign required
732 if ((rpt * gpt) < 0.)
733 return 0.;
734
735 Double_t rating_distY = -0.1 * distY + 1.;
736 Double_t rating_distZ = -0.025 * distZ + 1.;
737 Double_t rating_ptDiff = 1. - TMath::Abs(PtDiffRel(rpt, gpt));
738
739 if (rating_ptDiff < 0.)
740 rating_ptDiff = 0.2;
741
742 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}
752
753
754void AliTRDonlineTrackMatching::SetEsdTrackDefaultCuts(const char* cutIdent) {
755
756 if (strcmp(cutIdent, "strict") == 0){
757
758#ifdef TRD_TM_DEBUG
759 printf("AliTRDonlineTrackMatching -- default track cuts selected");
760#endif
761
762 fEsdTrackCutMinimal = kFALSE;
763 fEsdTrackCutPrim = kFALSE;
764
765 fEsdTrackCutMinTPCrows = 70;
766 fEsdTrackCutRequireTPCrefit = kTRUE;
767 fEsdTrackCutMinRatioRowsFindableClusters = 0.8;
768 fEsdTrackCutMaxChi2TPCclusters = 4.;
769 fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 36.;
770
771 fEsdTrackCutRequireITSrefit = kFALSE;
772 fEsdTrackCutMaxChi2ITSclusters = 36.;
773
774 fEsdTrackCutMaxDCAtoVertexXY = 1000.;
775 fEsdTrackCutMaxDCAtoVertexZ = 2.;
776 fEsdTrackCutsITSlayerMask = 0x0;
777
778 fEsdTrackCutPtDCAOfs = 0.0105;
779 fEsdTrackCutPtDCACoeff = 0.0350;
780 } else if (strcmp(cutIdent, "minimal") == 0){
781
782#ifdef TRD_TM_DEBUG
783 printf("AliTRDonlineTrackMatching -- minimal track cuts selected\n");
784#endif
785
786 fEsdTrackCutMinimal = kFALSE;
787 fEsdTrackCutPrim = kFALSE;
788
789 fEsdTrackCutMinTPCrows = 70;
790 fEsdTrackCutRequireTPCrefit = kTRUE;
791 fEsdTrackCutMinRatioRowsFindableClusters = 0.;
792 fEsdTrackCutMaxChi2TPCclusters = 100.;
793 fEsdTrackVCutsChi2TPCconstrainedVsGlobal = 1000.;
794
795 fEsdTrackCutRequireITSrefit = kFALSE;
796 fEsdTrackCutMaxChi2ITSclusters = 0.;
797
798 fEsdTrackCutMaxDCAtoVertexXY = 1000.;
799 fEsdTrackCutMaxDCAtoVertexZ = 1000.;
800 fEsdTrackCutsITSlayerMask = 0x0;
801 } else
633f3ec4 802 AliErrorClass("invalid cut set");
803
30bd9a5a 804}