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