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