TRD/AliTRDtrackV1.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 #include "AliESDtrack.h"
  19 #include "AliTracker.h"
  20
  21 #include "AliTRDtrackV1.h"
  22 #include "AliTRDcluster.h"
  23 #include "AliTRDcalibDB.h"
  24 #include "AliTRDrecoParam.h"
  25
  26 ClassImp(AliTRDtrackV1)
  27
  28 ///////////////////////////////////////////////////////////////////////////////
  29 //                                                                           //
  30 //  Represents a reconstructed TRD track                                     //
  31 //  Local TRD Kalman track                                                   //
  32 //                                                                           //
  33 //  Authors:                                                                 //
  34 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  35 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  36 //                                                                           //
  37 ///////////////////////////////////////////////////////////////////////////////
  38
  39 //_______________________________________________________________
  40 AliTRDtrackV1::AliTRDtrackV1()
  41   :AliTRDtrack()
  42 {
  43   //
  44   // Default constructor
  45   //
  46
  47         for(int ip=0; ip<6; ip++){
  48                 fTrackletIndex[ip] = -1;
  49                 fTracklet[ip].Reset();
  50         }
  51 }
  52
  53 //_______________________________________________________________
  54 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t)
  55   :AliTRDtrack(t)
  56 {
  57   //
  58   // Constructor from AliESDtrack
  59   //
  60
  61         //AliInfo(Form("alpha %f", GetAlpha()));
  62         t.GetTRDtracklets(&fTrackletIndex[0]);
  63         for(int ip=0; ip<6; ip++) fTracklet[ip].Reset();
  64 }
  65
  66 //_______________________________________________________________
  67 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref)
  68   :AliTRDtrack(ref)
  69 {
  70   //
  71   // Copy constructor
  72   //
  73
  74         for(int ip=0; ip<6; ip++){
  75                 fTrackletIndex[ip] = ref.fTrackletIndex[ip];
  76                 fTracklet[ip]      = ref.fTracklet[ip];
  77         }
  78 }
  79
  80 //_______________________________________________________________
  81 // AliTRDtrackV1::~AliTRDtrackV1()
  82 // {
  83 //
  84 // }
  85
  86 //_______________________________________________________________
  87 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
  88              , Double_t x, Double_t alpha)
  89   :AliTRDtrack()
  90 {
  91   //
  92   // The stand alone tracking constructor
  93   // TEMPORARY !!!!!!!!!!!
  94   // to check :
  95   // 1. covariance matrix
  96   // 2. dQdl calculation
  97   //
  98
  99   Double_t cnv   = 1.0 / (GetBz() * kB2C);
 100
 101   Double_t pp[5] = { p[0]
 102                    , p[1]
 103                    , x*p[4] - p[2]
 104                    , p[3]
 105                    , p[4]*cnv      };
 106
 107   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
 108   Double_t c32 =   x*cov[13] -     cov[ 8];
 109   Double_t c20 =   x*cov[10] -     cov[ 3];
 110   Double_t c21 =   x*cov[11] -     cov[ 4];
 111   Double_t c42 =   x*cov[14] -     cov[12];
 112
 113   Double_t cc[15] = { cov[ 0]
 114                     , cov[ 1],     cov[ 2]
 115                     , c20,         c21,         c22
 116                     , cov[ 6],     cov[ 7],     c32,     cov[ 9]
 117                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
 118
 119         Set(x,alpha,pp,cc);
 120   Double_t s = GetSnp();
 121   Double_t t = GetTgl();
 122
 123         Int_t ncl = 0, nclPlane; AliTRDcluster *c = 0x0;
 124         for(int iplane=0; iplane<6; iplane++){
 125                 fTrackletIndex[iplane] = -1;
 126                 fTracklet[iplane] = trklts[iplane];
 127                 nclPlane = 0;
 128                 for(int ic = 0; ic<AliTRDseed::knTimebins; ic++){
 129                         if(!fTracklet[iplane].IsUsable(ic)) continue;
 130                         if(!(c = fTracklet[iplane].GetClusters(ic))) continue;
 131
 132                         fIndex[ncl]      = fTracklet[iplane].GetIndexes(ic);
 133                         Double_t q = TMath::Abs(c->GetQ());
 134                         fClusters[ncl]   = c;
 135                         // temporary !!!
 136                         // This is not correctly. Has to be updated in the AliTRDtrackerV1::FollowBackProlonagation()
 137                         fdQdl[ncl]       = q * (s*s < 1.) ? TMath::Sqrt((1-s*s)/(1+t*t)) : 1.;
 138                         ncl++;
 139                         nclPlane++;
 140                 }
 141                 //printf("%d N clusters plane %d [%d %d].\n", iplane, nclPlane, fTracklet[iplane].GetN2(), trklts[iplane].GetN());
 142         }
 143         //printf("N clusters in AliTRDtrackV1 %d.\n", ncl);
 144         SetNumberOfClusters(/*ncl*/);
 145 }
 146
 147 //_______________________________________________________________
 148 Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
 149 {
 150         // set MC label for this tracklet
 151
 152   Int_t s[kMAXCLUSTERSPERTRACK][2];
 153   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 154     s[i][0] = -1;
 155     s[i][1] =  0;
 156   }
 157
 158   Bool_t labelAdded;
 159         Int_t label;
 160         AliTRDcluster *c    = 0x0;
 161   for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
 162                 if(fTrackletIndex[ip] < 0) continue;
 163                 for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
 164                         if(!(c = fTracklet[ip].GetClusters(ic))) continue;
 165                         for (Int_t k = 0; k < 3; k++) {
 166                                 label      = c->GetLabel(k);
 167                                 labelAdded = kFALSE;
 168                                 Int_t j = 0;
 169                                 if (label >= 0) {
 170                                         while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
 171                                                 if ((s[j][0] == label) ||
 172                                                                 (s[j][1] ==     0)) {
 173                                                         s[j][0] = label;
 174                                                         s[j][1]++;
 175                                                         labelAdded = kTRUE;
 176                                                 }
 177                                                 j++;
 178                                         }
 179                                 }
 180                         }
 181                 }
 182         }
 183
 184   Int_t max = 0;
 185   label = -123456789;
 186   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 187     if (s[i][1] <= max) continue;
 188                 max   = s[i][1];
 189                 label = s[i][0];
 190   }
 191
 192   if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
 193
 194   SetLabel(label);
 195
 196         return kTRUE;
 197 }
 198
 199 //_______________________________________________________________
 200 Bool_t AliTRDtrackV1::CookPID()
 201 {
 202   //
 203   // Cook the PID information
 204   //
 205
 206 // CookdEdx();  // truncated mean ... do we still need it ?
 207
 208 // CookdEdxTimBin(seed->GetID());
 209
 210   // Sets the a priori probabilities
 211   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 212     fPID[ispec] = 1.0 / AliPID::kSPECIES;
 213   }
 214
 215         // steer PID calculation @ tracklet level
 216         Double_t *prob = 0x0;
 217         fPIDquality = 0;
 218         for(int itrklt=0; itrklt<AliESDtrack::kNPlane; itrklt++){
 219     //for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
 220
 221                 if(fTrackletIndex[itrklt]<0) continue;
 222                 if(!fTracklet[itrklt].IsOK()) continue;
 223                 if(!(prob = fTracklet[itrklt].GetProbability())) return kFALSE;
 224
 225                 Int_t nspec = 0; // quality check of tracklet dEdx
 226                 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 227                         if(prob[ispec] < 0.) continue;
 228                         fPID[ispec] *= prob[ispec];
 229                         nspec++;
 230                 }
 231                 if(!nspec) continue;
 232
 233                 fPIDquality++;
 234         }
 235
 236   // no tracklet found for PID calculations
 237   if(!fPIDquality) return kTRUE;
 238
 239         // slot for PID calculation @ track level
 240
 241
 242   // normalize probabilities
 243   Double_t probTotal = 0.0;
 244   for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
 245
 246
 247   if (probTotal <= 0.0) {
 248     AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
 249     return kFALSE;
 250   }
 251
 252   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
 253
 254   return kTRUE;
 255 }
 256
 257 //_______________________________________________________________
 258 Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
 259 {
 260   //
 261   // Get the momentum at a given plane
 262   //
 263
 264         return plane >=0 && plane < 6 && fTrackletIndex[plane] >= 0 ? fTracklet[plane].GetMomentum() : -1.;
 265 }
 266
 267 //_______________________________________________________________
 268 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 269 {
 270   //
 271   // Get the predicted chi2
 272   //
 273
 274   Double_t x      = trklt->GetX0();
 275   Double_t p[2]   = { trklt->GetYat(x)
 276                     , trklt->GetZat(x) };
 277   Double_t cov[3];
 278         trklt->GetCovAt(x, cov);
 279
 280   return AliExternalTrackParam::GetPredictedChi2(p, cov);
 281
 282 }
 283
 284 //_______________________________________________________________
 285 Bool_t AliTRDtrackV1::IsOwner() const
 286 {
 287   //
 288   // Check whether track owns the tracklets
 289   //
 290
 291         for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
 292                 if(fTrackletIndex[ip] < 0) continue;
 293                 if(!fTracklet[ip].IsOwner()) return kFALSE;
 294         }
 295         return kTRUE;
 296 }
 297
 298 //___________________________________________________________
 299 void AliTRDtrackV1::SetNumberOfClusters()
 300 {
 301 // Calculate the number of clusters attached to this track
 302
 303         Int_t ncls = 0;
 304         for(int ip=0; ip<6; ip++){
 305                 if(fTrackletIndex[ip] >= 0) ncls += fTracklet[ip].GetN();
 306         }
 307         AliKalmanTrack::SetNumberOfClusters(ncls);
 308 }
 309
 310
 311 //_______________________________________________________________
 312 void AliTRDtrackV1::SetOwner(Bool_t own)
 313 {
 314   //
 315   // Toggle ownership of tracklets
 316   //
 317
 318         for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
 319                 if(fTrackletIndex[ip] < 0) continue;
 320                 //AliInfo(Form("p[%d] index[%d]", ip, fTrackletIndex[ip]));
 321                 fTracklet[ip].SetOwner(own);
 322         }
 323 }
 324
 325 //_______________________________________________________________
 326 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
 327 {
 328   //
 329   // Set the tracklets
 330   //
 331
 332         if(plane < 0 || plane >= AliESDtrack::kNPlane) return;
 333         fTracklet[plane]      = (*trklt);
 334         fTrackletIndex[plane] = index;
 335 }
 336
 337 //_______________________________________________________________
 338 Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 339 {
 340   //
 341   // Update track parameters
 342   //
 343
 344   Double_t x      = trklt->GetX0();
 345   Double_t p[2]   = { trklt->GetYat(x)
 346                     , trklt->GetZat(x) };
 347   Double_t cov[3];
 348         trklt->GetCovAt(x, cov);
 349
 350 //      Print();
 351 //      AliInfo(Form("cov[%f %f %f]", cov[0], cov[1], cov[2]));
 352
 353   if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
 354         //Print();
 355
 356         AliTRDcluster *c = 0x0;
 357         Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
 358         AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
 359
 360
 361   // Register info to track
 362 //   Int_t n      = GetNumberOfClusters();
 363 //   fIndex[n]    = index;
 364 //   fClusters[n] = c;
 365   SetNumberOfClusters(/*GetNumberOfClusters()+trklt->GetN()*/);
 366   SetChi2(GetChi2() + chisq);
 367
 368         // update tracklet
 369         trklt->SetMomentum(GetP());
 370         trklt->SetSnp(GetSnp());
 371         trklt->SetTgl(GetTgl());
 372         return kTRUE;
 373 }
 374
 375 //_______________________________________________________________
 376 void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
 377 {
 378   //
 379   // Update the ESD track
 380   //
 381
 382         // copy dEdx to ESD
 383         Float_t *dedx = 0x0;
 384         for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
 385                 if(fTrackletIndex[ip] < 0) continue;
 386                 fTracklet[ip].CookdEdx(AliESDtrack::kNSlice);
 387                 dedx = fTracklet[ip].GetdEdx();
 388                 for (Int_t js = 0; js < AliESDtrack::kNSlice; js++) track->SetTRDsignals(dedx[js], ip, js);
 389                 //track->SetTRDTimBin(fTimBinPlane[i], i);
 390         }
 391
 392         // copy PID to ESD
 393         track->SetTRDpid(fPID);
 394         track->SetTRDpidQuality(fPIDquality);
 395 }