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(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, 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, 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   if (!esdTrack)
 245     return kFALSE;
 246
 247   UInt_t status = esdTrack->GetStatus();
 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
 320 Bool_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
 540 Bool_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
 579 Int_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
 687 Double_t AliTRDonlineTrackMatching::PtDiffRel(Double_t refPt, Double_t gtuPt){
 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
 698 Double_t AliTRDonlineTrackMatching::RateTrackMatch(Double_t distY, Double_t distZ, Double_t rpt, Double_t gpt){
 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
 729 void 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 }