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