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