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