TRD/AliTRDgtuTrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //                                                                           //
  21 //  TRD module global track (GTU)                                            //
  22 //                                                                           //
  23 //                                                                           //
  24 ///////////////////////////////////////////////////////////////////////////////
  25
  26 #include <TMath.h>
  27 #include <TMatrixD.h>
  28 #include <TObjArray.h>
  29
  30 #include "AliLog.h"
  31
  32 #include "AliTRDgeometry.h"
  33 #include "AliTRDcalibDB.h"
  34 #include "AliTRDltuTracklet.h"
  35 #include "AliTRDgtuTrack.h"
  36 #include "Cal/AliTRDCalPID.h"
  37
  38 ClassImp(AliTRDgtuTrack)
  39
  40 //_____________________________________________________________________________
  41 AliTRDgtuTrack::AliTRDgtuTrack()
  42   :TObject()
  43   ,fTracklets(new TObjArray(400))
  44   ,fYproj(0)
  45   ,fZproj(0)
  46   ,fSlope(0)
  47   ,fDetector(-1)
  48   ,fNtracklets(0)
  49   ,fNplanes(0)
  50   ,fNclusters(0)
  51   ,fPt(0)
  52   ,fPhi(0)
  53   ,fEta(0)
  54   ,fLabel(-1)
  55   ,fPID(0)
  56   ,fIsElectron(kFALSE)
  57 {
  58   //
  59   // Default constructor
  60   //
  61
  62 }
  63
  64 //_____________________________________________________________________________
  65 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
  66   :TObject(t)
  67   ,fTracklets(NULL)
  68   ,fYproj(t.fYproj)
  69   ,fZproj(t.fZproj)
  70   ,fSlope(t.fSlope)
  71   ,fDetector(t.fDetector)
  72   ,fNtracklets(t.fNtracklets)
  73   ,fNplanes(t.fNplanes)
  74   ,fNclusters(t.fNclusters)
  75   ,fPt(t.fPt)
  76   ,fPhi(t.fPhi)
  77   ,fEta(t.fEta)
  78   ,fLabel(t.fLabel)
  79   ,fPID(t.fPID)
  80   ,fIsElectron(t.fIsElectron)
  81 {
  82   //
  83   // Copy contructor
  84   //
  85
  86 }
  87
  88 //_____________________________________________________________________________
  89 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
  90 {
  91   //
  92   // Assignment operator
  93   //
  94
  95   if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this);
  96
  97   return *this;
  98
  99 }
 100
 101 //_____________________________________________________________________________
 102 void AliTRDgtuTrack::Copy(TObject &t) const
 103 {
 104   //
 105   // Copy function
 106   //
 107
 108   ((AliTRDgtuTrack &) t).fTracklets  = NULL;
 109   ((AliTRDgtuTrack &) t).fYproj      = fYproj;
 110   ((AliTRDgtuTrack &) t).fZproj      = fZproj;
 111   ((AliTRDgtuTrack &) t).fSlope      = fSlope;
 112   ((AliTRDgtuTrack &) t).fDetector   = fDetector;
 113   ((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
 114   ((AliTRDgtuTrack &) t).fNplanes    = fNplanes;
 115   ((AliTRDgtuTrack &) t).fNclusters  = fNclusters;
 116   ((AliTRDgtuTrack &) t).fPt         = fPt;
 117   ((AliTRDgtuTrack &) t).fPhi        = fPhi;
 118   ((AliTRDgtuTrack &) t).fEta        = fEta;
 119   ((AliTRDgtuTrack &) t).fLabel      = fLabel;
 120   ((AliTRDgtuTrack &) t).fPID        = fPID;
 121   ((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;
 122
 123 }
 124
 125 //_____________________________________________________________________________
 126 AliTRDgtuTrack::~AliTRDgtuTrack()
 127 {
 128   //
 129   // Destructor
 130   //
 131
 132   if (fTracklets) {
 133     //fTracklets->Delete();
 134     fTracklets->Clear();
 135     delete fTracklets;
 136     fTracklets = 0;
 137   }
 138
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
 143 {
 144   //
 145   // Add a LTU tracklet to this track
 146   //
 147
 148   Tracklets()->Add(trk);
 149
 150 }
 151
 152 //_____________________________________________________________________________
 153 AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
 154 {
 155   //
 156   // Return LTU tracklet at position "pos"
 157   //
 158
 159   if (fTracklets == 0) {
 160     return 0;
 161   }
 162   void *trk = fTracklets->UncheckedAt(pos);
 163   if (trk == 0) {
 164     return 0;
 165   }
 166
 167   return (AliTRDltuTracklet *) trk;
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 Int_t AliTRDgtuTrack::Compare(const TObject *o) const
 173 {
 174   //
 175   // Compare function for sorting the tracks
 176   //
 177
 178   AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
 179
 180   if (fYproj <  gtutrack->GetYproj()) {
 181     return -1;
 182   }
 183   if (fYproj == gtutrack->GetYproj()) {
 184     return  0;
 185   }
 186
 187   return +1;
 188
 189 }
 190
 191 //_____________________________________________________________________________
 192 void AliTRDgtuTrack::Reset()
 193 {
 194   //
 195   // Reset the track information
 196   //
 197
 198   fYproj      = 0.0;
 199   fZproj      = 0.0;
 200   fSlope      = 0.0;
 201   fDetector   = -1;
 202   fNtracklets = 0;
 203   fNplanes    = 0;
 204   fNclusters  = 0;
 205   fPt         = 0.0;
 206   fPhi        = 0.0;
 207   fEta        = 0.0;
 208   fLabel      = -1;
 209   fPID        = 0.0;
 210   fIsElectron = kFALSE;
 211
 212 }
 213
 214 //_____________________________________________________________________________
 215 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
 216 {
 217   //
 218   // Calculate the kinematics of the found track
 219   //
 220
 221   Int_t  i    = 0;
 222   Int_t  iDet = 0;
 223   Int_t  nDet = 0;
 224   Bool_t newDetector;
 225
 226   AliTRDltuTracklet *trk;
 227   Int_t nTracklets = GetNtracklets();
 228   Float_t fC[kNmaxTrk][3];            // X, Y, Z  coordinates of segments
 229
 230   fYproj      = 0.0;
 231   fZproj      = 0.0;
 232   fSlope      = 0.0;
 233   fNclusters  = 0;
 234   fNplanes    = 0;
 235   fNtracklets = GetNtracklets();
 236   Int_t inDetector[kNplan];
 237   for (i = 0; i < kNplan; i++) {
 238     inDetector[i] = -1;
 239   }
 240
 241   for (i = 0; i < nTracklets; i++) {
 242
 243     trk         = GetTracklet(i);
 244     fYproj     += trk->GetYproj(xpl);
 245     fZproj     += trk->GetZproj(xpl);
 246     fSlope     += trk->GetSlope();
 247     fNclusters += trk->GetNclusters();
 248     iDet        = trk->GetDetector();
 249
 250     newDetector = kTRUE;
 251     for (Int_t id = 0; id < nDet; id++) {
 252       if (iDet == inDetector[id]) {
 253         newDetector = kFALSE;
 254         break;
 255       }
 256     }
 257     if (newDetector) {
 258       inDetector[nDet++] = iDet;
 259       fNplanes++;
 260     }
 261
 262     fC[i][0] = trk->GetTime0();
 263     fC[i][1] = trk->GetOffset();
 264     fC[i][2] = trk->GetRowz();
 265
 266   }
 267   fYproj /= (Float_t) nTracklets;
 268   fZproj /= (Float_t) nTracklets;
 269   fSlope /= (Float_t) nTracklets;
 270
 271   Float_t x[kNmaxTrk+1];
 272   Float_t y[kNmaxTrk+1];
 273   Float_t z[kNmaxTrk+1];
 274   Bool_t  count[kNmaxTrk];
 275   for (i = 0; i < kNmaxTrk; i++) {
 276     count[i] = kFALSE;
 277   }
 278
 279   Int_t iXmin = -1;
 280   Int_t j     =  0;
 281   x[0] = y[0] = z[0] = 0.0;
 282   while (j < nTracklets) {
 283     iXmin = -1;
 284     for (i = 0; i < nTracklets; i++) {
 285       if (count[i]) continue;
 286       if (iXmin == -1) {
 287         iXmin = i;
 288         continue;
 289       }
 290       if (fC[i][0] < fC[iXmin][0]) {
 291         iXmin = i;
 292       }
 293     }
 294     x[j+1] = fC[iXmin][0];
 295     y[j+1] = fC[iXmin][1];
 296     z[j+1] = fC[iXmin][2];
 297     j++;
 298     count[iXmin] = kTRUE;
 299   }
 300
 301   TMatrixD smatrix(2,2);
 302   TMatrixD sums(2,1);
 303   TMatrixD res(2,1);
 304   Double_t xv, yv;
 305   Float_t a, b;
 306
 307   smatrix.Zero();
 308   sums.Zero();
 309   for (i = 0; i < nTracklets; i++) {
 310     xv = (Double_t) x[i+1];
 311     yv = (Double_t) y[i+1];
 312     smatrix(0,0) += 1.0;
 313     smatrix(1,1) += xv*xv;
 314     smatrix(0,1) += xv;
 315     smatrix(1,0) += xv;
 316     sums(0,0)    += yv;
 317     sums(1,0)    += xv*yv;
 318   }
 319   res = smatrix.Invert() * sums;
 320   a   = res(0,0);
 321   b   = res(1,0);
 322
 323   Float_t dist  = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
 324   Float_t fx1   = x[1]          + dist * (Float_t) (nTracklets-1) / 6.0;
 325   Float_t fy1   = a + b * fx1;
 326   Float_t fx2   = x[nTracklets] - dist * (Float_t) (nTracklets-1) / 6.0;
 327   Float_t fy2   = a + b * fx2;
 328   Float_t d12   = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
 329   Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
 330   Float_t r     = (d12 / 2.0) / TMath::Sin(alpha);
 331
 332   fPt = 0.3 * field * 0.01 * r;
 333
 334   Float_t d1 = fx1*fx1 + fy1*fy1;
 335   Float_t d2 = fx2*fx2 + fy2*fy2;
 336   Float_t d  = fx1*fy2 - fx2*fy1;
 337
 338   Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
 339   Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
 340
 341   if (yc != 0.0) {
 342     fPhi = TMath::ATan(xc/yc);
 343   }
 344   else {
 345     fPhi = TMath::PiOver2();
 346   }
 347
 348   fPhi *= 180.0/TMath::Pi();
 349
 350   smatrix.Zero();
 351   sums.Zero();
 352   for (i = 0; i < nTracklets+1; i++) {
 353     xv = (Double_t) z[i];
 354     yv = (Double_t) x[i];
 355     smatrix(0,0) += 1.0;
 356     smatrix(1,1) += xv*xv;
 357     smatrix(0,1) += xv;
 358     smatrix(1,0) += xv;
 359     sums(0,0)    += yv;
 360     sums(1,0)    += xv*yv;
 361   }
 362   res = smatrix.Invert() * sums;
 363   a   = res(0,0);
 364   b   = res(1,0);
 365   Float_t theta = TMath::ATan(b);
 366
 367   if (theta < 0.0) {
 368     theta = TMath::Pi() + theta;
 369   }
 370   if (theta == 0.0) {
 371     fEta = 0.0;
 372   }
 373   else {
 374     fEta = -TMath::Log(TMath::Tan(theta/2.0));
 375   }
 376
 377 }
 378
 379 //_____________________________________________________________________________
 380 void AliTRDgtuTrack::MakePID()
 381 {
 382   //
 383   // Electron likelihood signal
 384   //
 385
 386   Int_t i = 0;
 387
 388   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 389   if (!calibration) {
 390     AliError("No instance of AliTRDcalibDB.");
 391     return;
 392   }
 393   const AliTRDCalPID *pd = calibration->GetPIDObject(1);
 394
 395   AliTRDltuTracklet *trk;
 396   Int_t   nTracklets = GetNtracklets();
 397   Int_t   det;
 398   Int_t   pla;
 399   Float_t sl;
 400   Float_t th;
 401   Float_t q;
 402   Float_t probPio = 1.0;
 403   Float_t probEle = 1.0;
 404
 405   for (i = 0; i < nTracklets; i++) {
 406
 407     trk = GetTracklet(i);
 408
 409     sl  = TMath::Abs(trk->GetSlope());     // Tracklet inclination in X-y plane
 410     th  = trk->GetRowz()/trk->GetTime0();  // Tracklet inclination in X-z plane
 411     th  = TMath::ATan(TMath::Abs(th));
 412
 413     q   = trk->GetQ()
 414         * TMath::Cos(sl/180.0*TMath::Pi())
 415         * TMath::Cos(th/180.0*TMath::Pi());
 416
 417     det = trk->GetDetector();
 418     pla = trk->GetPlane(det);
 419
 420     // Unclear yet factor to match the PS distributions = 5.8
 421     // not explained only by the tail filter ...
 422     // AliRoot v4-03-07 , v4-03-Release
 423     //q = q * 5.8;
 424
 425     // HEAD28Mar06
 426     // Temporary (B. Vulpescu):
 427     // The charge distributions do not match the new changes in simulation (A. Bercuci),
 428     // which are nevertheless now in agreement with the beam tests.
 429     // Some tricks will be used to still have reasonable results
 430     // To match the existing charge distributions, the charge per layer has to be modified
 431     // as folows:
 432     /*
 433       if (k == 0) {
 434       // electrons
 435       q = 4.3 * q + 95.0;
 436       } else {
 437       // others
 438       q = 4.2 * q + 70.0;
 439       }
 440     */
 441     // Since at tracking time we have no information on the particle type, we will modify
 442     // instead the charge distributions accordingly. This will slow down the sampling.
 443     // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
 444     // been regenerated with AliTRDCreateDummyCDB.C
 445     // The new PIDLQ data base has the version :
 446     // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
 447     // version v0_s1
 448
 449     // HEAD28Mar06 new distributions: factor = 6.0
 450
 451                 // A.Bercuci on 2nd May 2007
 452                 // Temporary modification to deal with more energy slices. The values
 453                 // attached to charge slices and track length are dummy
 454                 Float_t dedx[3];
 455                 dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
 456                 Float_t length = 3.7;
 457
 458     probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
 459     probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
 460
 461   }
 462
 463   if ((probEle + probPio) > 0.0) {
 464     fPID = probEle/(probEle + probPio);
 465   }
 466   else {
 467     fPID = 0.0;
 468   }
 469
 470   // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
 471   // P [GeV/c]  fPIDcut (between 0 and 1)
 472   // 2          0.925
 473   // 3          0.915
 474   // 4          0.875
 475   // 5          0.855
 476   // 6          0.845
 477   // 8          0.785
 478   // 10         0.735
 479   //
 480   // PIDcut = 0.978 - 0.024 * P[GeV/c]
 481   //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
 482
 483   // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
 484   //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
 485
 486   // HEAD28Mar06 with new generated distributions (pol2 fit)
 487   Float_t absPt   = TMath::Abs(fPt);
 488   //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt;  // 800 bins in dEdx
 489   Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt;    // 200 bins in dEdx
 490
 491   fIsElectron = kFALSE;
 492   if (fPID >= fPIDcut) {
 493     fIsElectron = kTRUE;
 494   }
 495
 496 }
 497
 498 //_____________________________________________________________________________
 499 void AliTRDgtuTrack::CookLabel()
 500 {
 501   //
 502   // Cook the track label from tracklets labels
 503   //
 504
 505   AliTRDltuTracklet *trk;
 506
 507   const Int_t kMaxTracks = 10;
 508   Int_t trackLabel[kMaxTracks];
 509   Int_t trackCount[kMaxTracks];
 510   for (Int_t it = 0; it < kMaxTracks; it++) {
 511     trackLabel[it] = -1;
 512     trackCount[it] =  0;
 513   }
 514
 515   Bool_t counted;
 516   Int_t  label;
 517   Int_t  nTracks = 0;
 518   for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
 519
 520     trk = GetTracklet(itrk);
 521
 522     if (trk->GetLabel() == -1) {
 523       continue;
 524     }
 525
 526     label = trk->GetLabel();
 527
 528     counted = kFALSE;
 529     for (Int_t it = 0; it < nTracks; it++) {
 530       if (label == trackLabel[it]) {
 531         trackCount[it]++;
 532         counted = kTRUE;
 533         break;
 534       }
 535     }
 536     if (!counted) {
 537       trackLabel[nTracks] = label;
 538       trackCount[nTracks]++;
 539       nTracks++;
 540       if (nTracks == kMaxTracks) {
 541         Warning("CookLabel","Too many tracks for this tracklet.");
 542         nTracks--;
 543         break;
 544       }
 545     }
 546
 547   }
 548
 549   Float_t frac = 4.0 / 5.0;
 550   for (Int_t it = 0; it < kMaxTracks; it++) {
 551     if (trackCount[it] >= (Int_t) (frac*fNtracklets)) {
 552       fLabel = trackLabel[it];
 553       break;
 554     }
 555   }
 556
 557 }
 558
 559 //_____________________________________________________________________________
 560 void AliTRDgtuTrack::ResetTracklets()
 561 {
 562   //
 563   // Resets the list of tracklets
 564   //
 565
 566   if (fTracklets) {
 567     //fTracklets->Delete();
 568     fTracklets->Clear();
 569   }
 570
 571 }
 572
 573 //_____________________________________________________________________________
 574 TObjArray *AliTRDgtuTrack::Tracklets()
 575 {
 576   //
 577   // Returns the list of tracklets
 578   //
 579
 580   if (!fTracklets) {
 581     fTracklets = new TObjArray(400);
 582   }
 583
 584   return fTracklets;
 585
 586 }
 587
 588 //_____________________________________________________________________________
 589 Int_t AliTRDgtuTrack::GetNtracklets() const
 590 {
 591   //
 592   // Returns the number of tracklets
 593   //
 594
 595   if (fTracklets) {
 596     return fTracklets->GetEntriesFast();
 597   }
 598
 599   return 0;
 600
 601 }