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