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