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