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