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 }