PWG1/TRD/info/AliTRDtrackInfo.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: AliTRDtrackInfo.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Reconstruction QA                                                     //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  25 //                                                                        //
  26 ////////////////////////////////////////////////////////////////////////////
  27
  28 #include "TDatabasePDG.h"
  29 #include "TPDGCode.h"
  30 #include "TVectorT.h"
  31
  32 #include "AliTrackReference.h"
  33 #include "AliExternalTrackParam.h"
  34
  35 #include "AliTRDseedV1.h"
  36 #include "AliTRDtrackV1.h"
  37 #include "AliTRDgeometry.h"
  38 #include "AliTRDtrackerV1.h"
  39
  40 #include "AliTRDtrackInfo.h"
  41
  42 ClassImp(AliTRDtrackInfo)
  43 ClassImp(AliTRDtrackInfo::AliMCinfo)
  44 ClassImp(AliTRDtrackInfo::AliESDinfo)
  45 Double_t AliTRDtrackInfo::AliMCinfo::fgKalmanStep = 2.;
  46
  47 //___________________________________________________
  48 AliTRDtrackInfo::AliTRDtrackInfo():
  49   TObject()
  50   ,fNClusters(0)
  51   ,fTRDtrack(NULL)
  52   ,fMC(NULL)
  53   ,fESD()
  54 {
  55   //
  56   // Default constructor
  57   //
  58 }
  59
  60
  61 //___________________________________________________
  62 AliTRDtrackInfo::AliTRDtrackInfo(const AliTRDtrackInfo &trdInfo):
  63   TObject((const TObject&)trdInfo)
  64   ,fNClusters(trdInfo.fNClusters)
  65   ,fTRDtrack(NULL)
  66   ,fMC(NULL)
  67   ,fESD(trdInfo.fESD)
  68 {
  69   //
  70   // copy Entries
  71   //
  72
  73   if(trdInfo.fMC) fMC = new AliMCinfo(*trdInfo.fMC);
  74
  75   if(trdInfo.fTRDtrack){
  76     fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
  77     if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
  78   }
  79 }
  80
  81 //___________________________________________________
  82 AliTRDtrackInfo::AliMCinfo::AliMCinfo()
  83   :fLabel(0)
  84   ,fPDG(0)
  85   ,fNTrackRefs(0)
  86 {
  87   // Set 0-Pointers
  88   memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
  89 }
  90
  91 //___________________________________________________
  92 AliTRDtrackInfo::AliMCinfo::AliMCinfo(const AliMCinfo &mc)
  93   :fLabel(mc.fLabel)
  94   ,fPDG(mc.fPDG)
  95   ,fNTrackRefs(mc.fNTrackRefs)
  96 {
  97   //
  98   // Constructor
  99   //
 100
 101   memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
 102   for(Int_t ien = 0; ien < 12; ien++){
 103     if(mc.fTrackRefs[ien])
 104       fTrackRefs[ien] = new AliTrackReference(*(mc.fTrackRefs[ien]));
 105   }
 106 }
 107
 108 //________________________________________________________
 109 Int_t AliTRDtrackInfo::AliMCinfo::GetPID() const
 110 {
 111   switch(fPDG){
 112   case kElectron:
 113   case kPositron: return AliPID::kElectron;
 114   case kMuonPlus:
 115   case kMuonMinus: return AliPID::kMuon;
 116   case kPiPlus:
 117   case kPiMinus: return AliPID::kPion;
 118   case kKPlus:
 119   case kKMinus: return AliPID::kKaon;
 120   case kProton:
 121   case kProtonBar: return AliPID::kProton;
 122   }
 123   return -1;
 124 }
 125
 126 //___________________________________________________
 127 AliTRDtrackInfo::AliESDinfo::AliESDinfo()
 128   :fHasV0(0)
 129   ,fId(-1)
 130   ,fStatus(0)
 131   ,fKinkIndex(0)
 132   ,fTPCncls(0)
 133   ,fTRDpidQuality(0)
 134   ,fTRDnSlices(0)
 135   ,fTRDslices(NULL)
 136   ,fOP(NULL)
 137 {
 138   //
 139   // Constructor
 140   //
 141
 142   memset(fTRDr, 0, AliPID::kSPECIES*sizeof(Double32_t));
 143 }
 144
 145 //___________________________________________________
 146 AliTRDtrackInfo::AliESDinfo::AliESDinfo(const AliESDinfo &esd)
 147   :fHasV0(esd.fHasV0)
 148   ,fId(esd.fId)
 149   ,fStatus(esd.fStatus)
 150   ,fKinkIndex(esd.fKinkIndex)
 151   ,fTPCncls(esd.fTPCncls)
 152   ,fTRDpidQuality(esd.fTRDpidQuality)
 153   ,fTRDnSlices(esd.fTRDnSlices)
 154   ,fTRDslices(NULL)
 155   ,fOP(NULL)
 156 {
 157   //
 158   // Constructor
 159   //
 160
 161   memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
 162
 163   if(fTRDnSlices){
 164     fTRDslices = new Double32_t[fTRDnSlices];
 165     memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
 166   }
 167   if(esd.fOP) fOP = new AliExternalTrackParam(*esd.fOP);
 168 }
 169
 170
 171 //___________________________________________________
 172 AliTRDtrackInfo::~AliTRDtrackInfo()
 173 {
 174   //
 175   // Destructor
 176   //
 177
 178   if(fMC) delete fMC;
 179   if(fTRDtrack) delete fTRDtrack;
 180 }
 181
 182 //___________________________________________________
 183 AliTRDtrackInfo::AliMCinfo::~AliMCinfo()
 184 {
 185   //
 186   // Destructor
 187   //
 188
 189   fNTrackRefs = 0;
 190   for(Int_t ien = 0; ien < 12; ien++){
 191     if(fTrackRefs[ien]) delete fTrackRefs[ien];
 192     fTrackRefs[ien] = NULL;
 193   }
 194 }
 195
 196 //___________________________________________________
 197 AliTRDtrackInfo::AliESDinfo::~AliESDinfo()
 198 {
 199   //
 200   // Destructor
 201   //
 202
 203   if(fTRDnSlices){
 204     delete [] fTRDslices;
 205     fTRDslices = NULL;
 206     fTRDnSlices = 0;
 207   }
 208   if(fOP) delete fOP; fOP = NULL;
 209 }
 210
 211 //___________________________________________________
 212 void AliTRDtrackInfo::AliESDinfo::Delete(const Option_t *){
 213   //
 214   // Delete Pointer members
 215   //
 216   if(fTRDnSlices){
 217     delete [] fTRDslices;
 218     fTRDslices = NULL;
 219     fTRDnSlices = 0;
 220   }
 221   if(fOP) delete fOP; fOP = NULL;
 222 }
 223
 224
 225 //___________________________________________________
 226 AliTRDtrackInfo& AliTRDtrackInfo::operator=(const AliTRDtrackInfo &trdInfo)
 227 {
 228   //
 229   // = Operator
 230   //
 231
 232   fNClusters  = trdInfo.fNClusters;
 233   fESD = trdInfo.fESD;
 234
 235   if(trdInfo.fMC){
 236     if(!fMC)
 237       fMC = new AliMCinfo(*trdInfo.fMC);
 238     else
 239       new(fMC) AliMCinfo(*trdInfo.fMC);
 240   }
 241
 242   if(trdInfo.fTRDtrack){
 243     if(!fTRDtrack)
 244       fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
 245     else
 246       new(fTRDtrack) AliTRDtrackV1(*trdInfo.fTRDtrack);
 247     if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
 248   }
 249
 250   return *this;
 251 }
 252
 253 //___________________________________________________
 254 AliTRDtrackInfo::AliMCinfo& AliTRDtrackInfo::AliMCinfo::operator=(const AliMCinfo &mc)
 255 {
 256   //
 257   // Assignment operator
 258   //
 259
 260   fLabel      = mc.fLabel;
 261   fPDG        = mc.fPDG;
 262   fNTrackRefs = mc.fNTrackRefs;
 263
 264   AliTrackReference **itr = &fTrackRefs[0];
 265   AliTrackReference* const *jtr = &mc.fTrackRefs[0];
 266   for(Int_t ien = 0; ien < 12; ien++, itr++, jtr++){
 267     if((*jtr)){
 268       if(!(*itr)) (*itr) = new AliTrackReference(*(*jtr));
 269       else new(&(*itr)) AliTrackReference(*(*jtr));
 270     } else (*itr) = NULL;
 271   }
 272   return *this;
 273 }
 274
 275 //___________________________________________________
 276 AliTRDtrackInfo::AliESDinfo& AliTRDtrackInfo::AliESDinfo::operator=(const AliESDinfo &esd)
 277 {
 278   //
 279   // Assignment operator
 280   //
 281
 282   fId          = esd.fId;
 283   fStatus      = esd.fStatus;
 284   fTRDpidQuality = esd.fTRDpidQuality;
 285   fTRDnSlices  = esd.fTRDnSlices;
 286   fTRDslices   = NULL;
 287
 288   memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
 289
 290   if(fTRDnSlices){
 291     fTRDslices = new Double32_t[fTRDnSlices];
 292     memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
 293   }
 294   if(esd.fOP){
 295     if(fOP) new(fOP) AliExternalTrackParam(esd.fOP);
 296     else fOP = new AliExternalTrackParam(esd.fOP);
 297   } else fOP = NULL;
 298
 299   return *this;
 300 }
 301
 302 //___________________________________________________
 303 void AliTRDtrackInfo::Delete(const Option_t *)
 304 {
 305   //
 306   // Delete
 307   //
 308
 309   fNClusters  = 0;
 310   if(fMC) delete fMC; fMC = NULL;
 311   fESD.Delete(NULL);
 312   if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
 313 }
 314
 315 //___________________________________________________
 316 void AliTRDtrackInfo::SetTrack(const AliTRDtrackV1 *track)
 317 {
 318   //
 319   // Set the TRD track
 320   //
 321
 322   if(!fTRDtrack)
 323     fTRDtrack = new AliTRDtrackV1(*track);
 324   else
 325     new(fTRDtrack)AliTRDtrackV1(*track);
 326   fTRDtrack->SetOwner();
 327   // Make a copy for the object in order to avoid ownership problems
 328 }
 329
 330 //___________________________________________________
 331 void AliTRDtrackInfo::AddTrackRef(const AliTrackReference *tref)
 332 {
 333   //
 334   // Add track reference
 335   //
 336
 337   if(fMC->fNTrackRefs >= 12){
 338     SetCurved();
 339     return;
 340   }
 341   // Make a copy for the object in order to avoid ownership problems
 342   fMC->fTrackRefs[fMC->fNTrackRefs++] = new AliTrackReference(*tref);
 343 }
 344
 345 //___________________________________________________
 346 AliTrackReference* AliTRDtrackInfo::GetTrackRef(Int_t idx) const
 347 {
 348 //
 349 // Returns a track reference
 350 //
 351   if(!fMC) return NULL;
 352   return (idx>=0 && idx < 12) ? fMC->fTrackRefs[idx] : NULL;
 353 }
 354
 355 //___________________________________________________
 356 AliTrackReference* AliTRDtrackInfo::GetTrackRef(AliTRDseedV1* const tracklet) const
 357 {
 358 //
 359 // Returns a track reference
 360 //
 361   if(!fMC) return NULL;
 362   Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
 363   AliTrackReference * const* jtr = &(fMC->fTrackRefs[0]);
 364   for(Int_t itr = 0; itr < fMC->fNTrackRefs; itr++, ++jtr){
 365     if(!(*jtr)) break;
 366     if(TMath::Abs(tracklet->GetX0() - (*jtr)->LocalX()) < cw) return (*jtr);
 367   }
 368   return NULL;
 369 }
 370
 371 //___________________________________________________
 372 Int_t AliTRDtrackInfo::GetNumberOfClusters() const
 373 {
 374   //
 375   // Returns the number of clusters
 376   //
 377
 378   Int_t n = 0;
 379   if(!fTRDtrack) return 0;
 380   if(fTRDtrack->GetNumberOfTracklets() == 0) return n;
 381   AliTRDseedV1 *tracklet = NULL;
 382   for(Int_t ip=0; ip<6; ip++){
 383     if(!(tracklet = const_cast<AliTRDseedV1 *>(fTRDtrack->GetTracklet(ip)))) continue;
 384     n+=tracklet->GetN();
 385   }
 386   return n;
 387 }
 388
 389
 390 //___________________________________________________
 391 void  AliTRDtrackInfo::SetOuterParam(const AliExternalTrackParam *op)
 392 {
 393   //
 394   // Set outer track parameters
 395   //
 396
 397   if(!op) return;
 398   if(fESD.fOP) new(fESD.fOP) AliExternalTrackParam(*op);
 399   else fESD.fOP = new AliExternalTrackParam(*op);
 400 }
 401
 402 //___________________________________________________
 403 Int_t AliTRDtrackInfo::GetNTracklets() const
 404 {
 405   //
 406   // Return the number of tracklets
 407   //
 408
 409   if(!fTRDtrack) return 0;
 410   return fTRDtrack->GetNumberOfTracklets();
 411 }
 412
 413 //___________________________________________________
 414 void AliTRDtrackInfo::SetSlices(Int_t n, Double32_t *s)
 415 {
 416   //
 417   // Set the slices
 418   //
 419   if(fESD.fTRDnSlices != n){
 420     fESD.fTRDnSlices = 0;
 421     delete [] fESD.fTRDslices;
 422     fESD.fTRDslices = NULL;
 423   }
 424
 425   if(!fESD.fTRDnSlices){
 426     fESD.fTRDnSlices = n;
 427     fESD.fTRDslices = new Double32_t[fESD.fTRDnSlices];
 428   }
 429
 430   memcpy(fESD.fTRDslices, s, n*sizeof(Double32_t));
 431 }
 432
 433 //___________________________________________________
 434 Bool_t AliTRDtrackInfo::AliMCinfo::GetDirections(Float_t &x0, Float_t &y0, Float_t &z0, Float_t &dydx, Float_t &dzdx, Float_t &pt, UChar_t &status) const
 435 {
 436 // Check for 2 track ref for the tracklet defined bythe radial position x0
 437 // The "status" is a bit map and gives a more informative output in case of failure:
 438 //   - 0 : everything is OK
 439 //   - BIT(0) : 0 track refs found
 440 //   - BIT(1) : 1 track reference found
 441 //   - BIT(2) : dx <= 0 between track references
 442 //   - BIT(3) : dx > 0 && dx < 3.7 - tangent tracks
 443
 444   status = 0;
 445   Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
 446   Int_t nFound = 0;
 447   AliTrackReference *tr[2] = {NULL, NULL};
 448   AliTrackReference * const* jtr = &fTrackRefs[0];
 449   for(Int_t itr = 0; itr < fNTrackRefs; itr++, ++jtr){
 450     if(!(*jtr)) break;
 451 /*
 452     if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, (*jtr)->LocalX());*/
 453     if(TMath::Abs(x0 - (*jtr)->LocalX()) > cw) continue;
 454     tr[nFound++] = (*jtr);
 455     if(nFound == 2) break;
 456   }
 457   if(nFound < 2){
 458     AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Missing track ref x0[%6.3f] nref[%d]", x0, nFound));
 459     if(!nFound) SETBIT(status, 0);
 460     else SETBIT(status, 1);
 461     return kFALSE;
 462   }
 463   if((pt=tr[1]->Pt()) < 1.e-3) return kFALSE;
 464
 465   Double_t dx = tr[1]->LocalX() - tr[0]->LocalX();
 466   if(dx <= 0.){
 467     AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Track ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]", tr[0]->LocalX(), tr[1]->LocalX()));
 468     SETBIT(status, 2);
 469     return kFALSE;
 470   }
 471   if(TMath::Abs(dx-AliTRDgeometry::CamHght()-AliTRDgeometry::CdrHght())>1.E-3) SETBIT(status, 3);
 472
 473   dydx = (tr[1]->LocalY() - tr[0]->LocalY()) / dx;
 474   dzdx = (tr[1]->Z() - tr[0]->Z()) / dx;
 475   //Float_t dx0 = tr[1]->LocalX() - x0;
 476   y0   =  tr[1]->LocalY()/* - dydx*dx0*/;
 477   z0   =  tr[1]->Z()/* - dzdx*dx0*/;
 478   x0   =  tr[1]->LocalX();
 479   return kTRUE;
 480 }
 481
 482 //___________________________________________________
 483 void AliTRDtrackInfo::AliMCinfo::PropagateKalman(TVectorD *dx, TVectorD *dy, TVectorD *dz, TVectorD *pt, TVectorD *dpt, TVectorD *c) const
 484 {
 485 // Propagate Kalman from the first TRD track reference to
 486 // last one and save residuals in the y, z and pt.
 487 //
 488 // This is to calibrate the dEdx and MS corrections
 489
 490   if(!fNTrackRefs) return;
 491   for(Int_t itr=kNTrackRefs; itr--;){
 492     (*dx)[itr] = -1.; (*dy)[itr] = 100.; (*dz)[itr] = 100.; (*dpt)[itr] = 100.;
 493   }
 494
 495   // Initialize TRD track to the first track reference
 496   AliTrackReference *tr(NULL);
 497   Int_t itr(0); while(!(tr = fTrackRefs[itr])) itr++;
 498   if(tr->Pt()<1.e-3) return;
 499
 500   AliTRDtrackV1 tt;
 501   Double_t xyz[3]={tr->X(),tr->Y(),tr->Z()};
 502   Double_t pxyz[3]={tr->Px(),tr->Py(),tr->Pz()};
 503   Double_t var[6] = {1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4};
 504   Double_t cov[21]={
 505     var[0],  0.,  0.,  0.,  0.,  0.,
 506          var[1],  0.,  0.,  0.,  0.,
 507               var[2],  0.,  0.,  0.,
 508                    var[3],  0.,  0.,
 509                         var[4],  0.,
 510                              var[5]
 511   };
 512   TDatabasePDG db;
 513   const TParticlePDG *pdg=db.GetParticle(fPDG);
 514   if(!pdg){
 515     AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::PropagateKalman()", Form("PDG entry missing for code %d. References for track %d", fPDG, fNTrackRefs));
 516     return;
 517   }
 518   tt.Set(xyz, pxyz, cov, Short_t(pdg->Charge()));
 519   tt.SetMass(pdg->Mass());
 520
 521   Double_t x0(tr->LocalX());
 522   const Double_t *cc(NULL);
 523   for(Int_t ip=0; itr<fNTrackRefs; itr++){
 524     if(!(tr = fTrackRefs[itr])) continue;
 525     if(!AliTRDtrackerV1::PropagateToX(tt, tr->LocalX(), fgKalmanStep)) continue;
 526
 527     //if(update) ...
 528     (*dx)[ip]  = tt.GetX() - x0;
 529     (*dy)[ip]  = tt.GetY() - tr->LocalY();
 530     (*dz)[ip]  = tt.GetZ() - tr->Z();
 531     (*pt)[ip]  = tr->Pt();
 532     (*dpt)[ip] = tt.Pt()- tr->Pt();
 533     cc = tt.GetCovariance();
 534     c->Use(ip*15, (ip+1)*15, cc);
 535     ip++;
 536   }
 537 }