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