EVE/alice-macros/esd_tracks.C

   1 // $Id$
   2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
   6  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
   7  * full copyright notice.                                                 *
   8  **************************************************************************/
   9
  10 #if !defined(__CINT__) || defined(__MAKECINT__)
  11
  12 #include <TString.h>
  13 #include <TSystem.h>
  14 #include <TROOT.h>
  15 #include <TMath.h>
  16 #include <TGListTree.h>
  17
  18 #include <TEveVSDStructs.h>
  19 #include <TEveManager.h>
  20 #include <TEveTrackPropagator.h>
  21
  22 #include <EveBase/AliEveTrack.h>
  23 #include <EveBase/AliEveTrackCounter.h>
  24 #include <EveBase/AliEveMagField.h>
  25 #include <EveBase/AliEveEventManager.h>
  26
  27 #include <AliESDEvent.h>
  28 #include <AliESDtrackCuts.h>
  29 #include <AliESDtrack.h>
  30 #include <AliExternalTrackParam.h>
  31
  32 #include <AliPWG0Helper.h>
  33
  34 #endif
  35
  36
  37 // Use inner-tpc track params when its refit failed.
  38 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
  39
  40 // Use magnetic-field as retrieved from GRP.
  41 Bool_t g_esd_tracks_true_field = kTRUE;
  42
  43 // Use Runge-Kutta track stepper.
  44 Bool_t g_esd_tracks_rk_stepper = kFALSE;
  45
  46
  47 //==============================================================================
  48
  49 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
  50                                 Float_t magF, Float_t maxR)
  51 {
  52   if (g_esd_tracks_true_field)
  53   {
  54     trkProp->SetMagFieldObj(new AliEveMagField);
  55   }
  56   else
  57   {
  58     trkProp->SetMagField(magF);
  59   }
  60   if (g_esd_tracks_rk_stepper)
  61   {
  62     trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
  63   }
  64   trkProp->SetMaxR(maxR);
  65 }
  66
  67 //==============================================================================
  68
  69 TString esd_track_title(AliESDtrack* t)
  70 {
  71   TString s;
  72
  73   Int_t label = t->GetLabel(), index = t->GetID();
  74   TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
  75   TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));
  76
  77   Double_t p[3], v[3];
  78   t->GetXYZ(v);
  79   t->GetPxPyPz(p);
  80   Double_t pt    = t->Pt();
  81   Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
  82   Double_t ptsq  = pt*pt;
  83   Double_t ptm   = pt / (1.0 + pt*ptsig);
  84   Double_t ptM   = pt / (1.0 - pt*ptsig);
  85
  86   s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
  87            "pT = %.3f + %.3f - %.3f [%.3f]\n"
  88            "P  = (%.3f, %.3f, %.3f)\n"
  89            "V  = (%.3f, %.3f, %.3f)\n",
  90            idx.Data(), lbl.Data(), t->Charge(), 0,
  91            pt, ptM - pt, pt - ptm, ptsig*ptsq,
  92            p[0], p[1], p[2],
  93            v[0], v[1], v[2]);
  94
  95   Int_t   o;
  96   s += "Det (in,out,refit,pid):\n";
  97   o  = AliESDtrack::kITSin;
  98   s += Form("ITS (%d,%d,%d,%d)  ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  99   o  = AliESDtrack::kTPCin;
 100   s += Form("TPC(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
 101   o  = AliESDtrack::kTRDin;
 102   s += Form("TRD(%d,%d,%d,%d) ",    t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
 103   o  = AliESDtrack::kTOFin;
 104   s += Form("TOF(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
 105   o  = AliESDtrack::kHMPIDout;
 106   s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
 107   s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
 108
 109   if (t->IsOn(AliESDtrack::kESDpid))
 110   {
 111     Double_t pid[5];
 112     t->GetESDpid(pid);
 113     s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
 114   }
 115
 116   return s;
 117 }
 118
 119 //==============================================================================
 120
 121 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
 122 {
 123   // Add additional track parameters as a path-mark to track.
 124
 125   if (tp == 0)
 126     return;
 127
 128   Double_t pbuf[3], vbuf[3];
 129   tp->GetXYZ(vbuf);
 130   tp->GetPxPyPz(pbuf);
 131
 132   TEvePathMark pm(TEvePathMark::kReference);
 133   pm.fV.Set(vbuf);
 134   pm.fP.Set(pbuf);
 135   track->AddPathMark(pm);
 136 }
 137
 138 //==============================================================================
 139
 140 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
 141 {
 142   // Make a standard track representation and put it into given container.
 143
 144   // Choose which parameters to use a track's starting point.
 145   // If gkFixFailedITSExtr is TRUE (FALSE by default) and
 146   // if ITS refit failed, take track parameters at inner TPC radius.
 147
 148   const AliExternalTrackParam* tp = at;
 149
 150   Bool_t innerTaken = kFALSE;
 151   if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
 152   {
 153     tp = at->GetInnerParam();
 154     innerTaken = kTRUE;
 155   }
 156
 157   AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
 158   track->SetAttLineAttMarker(cont);
 159   track->SetName(Form("AliEveTrack %d", at->GetID()));
 160   track->SetElementTitle(esd_track_title(at));
 161   track->SetSourceObject(at);
 162
 163   // Add inner/outer track parameters as path-marks.
 164   if (at->IsOn(AliESDtrack::kTPCrefit))
 165   {
 166     if ( ! innerTaken)
 167     {
 168       esd_track_add_param(track, at->GetInnerParam());
 169     }
 170     esd_track_add_param(track, at->GetOuterParam());
 171   }
 172
 173   return track;
 174 }
 175
 176
 177 //==============================================================================
 178 // esd_tracks()
 179 //==============================================================================
 180
 181 TEveTrackList* esd_tracks()
 182 {
 183   AliESDEvent* esd = AliEveEventManager::AssertESD();
 184
 185   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 186   cont->SetMainColor(6);
 187
 188   esd_track_propagator_setup(cont->GetPropagator(),
 189                              0.1*esd->GetMagneticField(), 520);
 190
 191   gEve->AddElement(cont);
 192
 193   Int_t count = 0;
 194   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 195   {
 196     ++count;
 197     AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
 198
 199     cont->AddElement(track);
 200   }
 201   cont->SetTitle(Form("N=%d", count));
 202   cont->MakeTracks();
 203
 204   gEve->Redraw3D();
 205
 206   return cont;
 207 }
 208
 209 TEveTrackList* esd_tracks_MI()
 210 {
 211   AliESDEvent* esd = AliEveEventManager::AssertESD();
 212
 213   TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
 214   cont->SetLineColor(5);
 215   gEve->AddElement(cont);
 216
 217   Int_t count = 0;
 218   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 219   {
 220     ++count;
 221     AliESDtrack* at = esd->GetTrack(n);
 222     AliEveTrack* l = new AliEveTrack(at, cont->GetPropagator());
 223     l->SetName(Form("ESDTrackMI %d", at->GetID()));
 224     l->SetElementTitle(esd_track_title(at));
 225     l->SetAttLineAttMarker(cont);
 226     l->SetSourceObject(at);
 227
 228     at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
 229
 230     l->SetLockPoints(kTRUE);
 231     cont->AddElement(l);
 232   }
 233   cont->SetTitle(Form("N=%d", count));
 234
 235   gEve->Redraw3D();
 236
 237   return cont;
 238 }
 239
 240
 241 //==============================================================================
 242 // esd_tracks_from_array()
 243 //==============================================================================
 244
 245 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
 246 {
 247   // Retrieves AliESDTrack's from collection.
 248   // See example usage with AliAnalysisTrackCuts in the next function.
 249
 250   if (esd == 0) esd = AliEveEventManager::AssertESD();
 251
 252   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 253   cont->SetMainColor(6);
 254
 255   esd_track_propagator_setup(cont->GetPropagator(),
 256                              0.1*esd->GetMagneticField(), 520);
 257
 258   gEve->AddElement(cont);
 259
 260   Int_t    count = 0;
 261   TIter    next(col);
 262   TObject *obj;
 263   while ((obj = next()) != 0)
 264   {
 265     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
 266     {
 267       Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
 268               obj->GetName(), obj->GetTitle());
 269       continue;
 270     }
 271
 272     ++count;
 273     AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
 274
 275     AliEveTrack* track = esd_make_track(at, cont);
 276     cont->AddElement(track);
 277   }
 278   cont->SetTitle(Form("N=%d", count));
 279   cont->MakeTracks();
 280
 281   gEve->Redraw3D();
 282
 283   return cont;
 284 }
 285
 286 void esd_tracks_alianalcuts_demo()
 287 {
 288   AliESDEvent* esd = AliEveEventManager::AssertESD();
 289
 290   AliESDtrackCuts atc;
 291   atc.SetPtRange(0.1, 5);
 292   atc.SetRapRange(-1, 1);
 293
 294   esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
 295 }
 296
 297
 298 //==============================================================================
 299 // esd_tracks_by_category
 300 //==============================================================================
 301
 302 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
 303 {
 304   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
 305   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 306
 307   Float_t b[2];
 308   Float_t bRes[2];
 309   Float_t bCov[3];
 310   esdTrack->GetImpactParameters(b,bCov);
 311   if (bCov[0] <= 0 || bCov[2] <= 0)
 312   {
 313     printf("Estimated b resolution lower or equal zero!\n");
 314     bCov[0] = bCov[2] = 0;
 315   }
 316   bRes[0] = TMath::Sqrt(bCov[0]);
 317   bRes[1] = TMath::Sqrt(bCov[2]);
 318
 319   // -----------------------------------
 320   // How to get to a n-sigma cut?
 321   //
 322   // The accumulated statistics from 0 to d is
 323   //
 324   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 325   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 326   //
 327   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 328   // Can this be expressed in a different way?
 329
 330   if (bRes[0] == 0 || bRes[1] == 0)
 331     return -1;
 332
 333   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 334
 335   // stupid rounding problem screws up everything:
 336   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 337   if (TMath::Exp(-d * d / 2) < 1e-10)
 338     return 1000;
 339
 340   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 341   return d;
 342 }
 343
 344 TEveElementList* esd_tracks_by_category()
 345 {
 346    // Import ESD tracks, separate them into several containers
 347    // according to primary-vertex cut and ITS&TPC refit status.
 348
 349    AliESDEvent* esd = AliEveEventManager::AssertESD();
 350
 351    TEveElementList* cont = new TEveElementList("ESD Tracks by category");
 352    gEve->AddElement(cont);
 353
 354    const Int_t   nCont = 9;
 355    const Float_t maxR  = 520;
 356    const Float_t magF  = 0.1*esd->GetMagneticField();
 357
 358    TEveTrackList *tl[nCont];
 359    Int_t          tc[nCont];
 360    Int_t          count = 0;
 361
 362    tl[0] = new TEveTrackList("Sigma < 3");
 363    tc[0] = 0;
 364    esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
 365    tl[0]->SetMainColor(3);
 366    cont->AddElement(tl[0]);
 367
 368    tl[1] = new TEveTrackList("3 < Sigma < 5");
 369    tc[1] = 0;
 370    esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
 371    tl[1]->SetMainColor(7);
 372    cont->AddElement(tl[1]);
 373
 374    tl[2] = new TEveTrackList("5 < Sigma");
 375    tc[2] = 0;
 376    esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
 377    tl[2]->SetMainColor(46);
 378    cont->AddElement(tl[2]);
 379
 380    tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
 381    tc[3] = 0;
 382    esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
 383    tl[3]->SetMainColor(41);
 384    cont->AddElement(tl[3]);
 385
 386    tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
 387    tc[4] = 0;
 388    esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
 389    tl[4]->SetMainColor(48);
 390    cont->AddElement(tl[4]);
 391
 392    tl[5] = new TEveTrackList("no TPC refit");
 393    tc[5] = 0;
 394    esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
 395    tl[5]->SetMainColor(kRed);
 396    cont->AddElement(tl[5]);
 397
 398    tl[6] = new TEveTrackList("ITS ncl>=3 & SPD Inner");
 399    tc[6] = 0;
 400    esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
 401    tl[6]->SetMainColor(kGreen);
 402    cont->AddElement(tl[6]);
 403
 404    tl[7] = new TEveTrackList("ITS ncl>=3 & b<3 cm");
 405    tc[7] = 0;
 406    esd_track_propagator_setup(tl[7]->GetPropagator(), magF, maxR);
 407    tl[7]->SetMainColor(kMagenta - 9);
 408    cont->AddElement(tl[7]);
 409
 410    tl[8] = new TEveTrackList("ITS others");
 411    tc[8] = 0;
 412    esd_track_propagator_setup(tl[8]->GetPropagator(), magF, maxR);
 413    tl[8]->SetMainColor(kRed + 2);
 414    cont->AddElement(tl[8]);
 415
 416    for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 417    {
 418       AliESDtrack* at = esd->GetTrack(n);
 419
 420       Float_t s  = get_sigma_to_vertex(at);
 421       Int_t   ti;
 422       if      (s <  3) ti = 0;
 423       else if (s <= 5) ti = 1;
 424       else             ti = 2;
 425
 426       Int_t    nclits;
 427       Double_t dtobeam;
 428
 429       if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
 430       {
 431          UChar_t itsclmap = at->GetITSClusterMap();
 432          Bool_t  spdinner = (itsclmap & 3) != 0;
 433
 434          nclits = 0;
 435          for (Int_t iter = 0; iter < 6; ++iter)
 436             if (itsclmap & (1 << iter)) nclits ++;
 437
 438          Double_t xyz[3];
 439          at->GetXYZ(xyz);
 440          dtobeam = TMath::Hypot(xyz[0], xyz[1]);
 441
 442          if ((nclits >= 3) && (spdinner))
 443             ti = 6;
 444          else if ((nclits >= 3) && (dtobeam < 3.0))
 445             ti = 7;
 446          else
 447             ti = 8;
 448       }
 449       else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
 450       {
 451          ti = 5;
 452       }
 453       else if ( ! at->IsOn(AliESDtrack::kITSrefit))
 454       {
 455          ti = (ti == 2) ? 4 : 3;
 456       }
 457
 458       TEveTrackList* tlist = tl[ti];
 459       ++tc[ti];
 460       ++count;
 461
 462       AliEveTrack* track = esd_make_track(at, tlist);
 463       if (ti<6)
 464          track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
 465       else
 466          track->SetName(Form("ESD Track idx=%d, dxy=%5.3f cl=%i", at->GetID(), dtobeam, nclits));
 467       tlist->AddElement(track);
 468    }
 469
 470    for (Int_t ti = 0; ti < nCont; ++ti)
 471    {
 472       TEveTrackList* tlist = tl[ti];
 473       tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
 474       tlist->SetTitle(Form("N tracks=%d", tc[ti]));
 475
 476       tlist->MakeTracks();
 477
 478       //    Bool_t good_cont = ti <= 1;
 479       Bool_t good_cont = ((ti == 6) || (ti == 7));
 480       if (AliEveTrackCounter::IsActive())
 481       {
 482          AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
 483       }
 484       else
 485       {
 486          if ( ! good_cont)
 487             tlist->SetLineStyle(6);
 488       }
 489    }
 490    cont->SetTitle(Form("N all tracks = %d", count));
 491    // ??? The following does not always work:
 492    cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
 493
 494    gEve->Redraw3D();
 495
 496    return cont;
 497 }
 498
 499
 500 //==============================================================================
 501 // esd_tracks_by_anal_cuts
 502 //==============================================================================
 503
 504 AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
 505
 506 TEveElementList* esd_tracks_by_anal_cuts()
 507 {
 508   AliESDEvent* esd = AliEveEventManager::AssertESD();
 509
 510   if (g_esd_tracks_anal_cuts == 0)
 511   {
 512     gSystem->Load("libPWG0base");
 513     gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/CreateStandardCuts.C");
 514     Int_t mode = AliPWG0Helper::kTPC;
 515     if (TMath::Abs(esd->GetMagneticField()) > 0.01)
 516       mode |= AliPWG0Helper::kFieldOn;
 517     gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
 518   }
 519
 520   TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
 521   gEve->AddElement(cont);
 522
 523   const Int_t   nCont = 2;
 524   const Float_t maxR  = 520;
 525   const Float_t magF  = 0.1*esd->GetMagneticField();
 526
 527   TEveTrackList *tl[nCont];
 528   Int_t          tc[nCont];
 529   Int_t          count = 0;
 530
 531   tl[0] = new TEveTrackList("Passed");
 532   tc[0] = 0;
 533   esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
 534   tl[0]->SetMainColor(3);
 535   cont->AddElement(tl[0]);
 536
 537   tl[1] = new TEveTrackList("Rejected");
 538   tc[1] = 0;
 539   esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
 540   tl[1]->SetMainColor(kRed);
 541   cont->AddElement(tl[1]);
 542
 543   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 544   {
 545     AliESDtrack* at = esd->GetTrack(n);
 546
 547     Float_t s  = get_sigma_to_vertex(at);
 548     Int_t   ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
 549
 550     TEveTrackList* tlist = tl[ti];
 551     ++tc[ti];
 552     ++count;
 553
 554     AliEveTrack* track = esd_make_track(at, tlist);
 555     track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
 556     tlist->AddElement(track);
 557   }
 558
 559   for (Int_t ti = 0; ti < nCont; ++ti)
 560   {
 561     TEveTrackList* tlist = tl[ti];
 562     tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
 563     tlist->SetTitle(Form("N tracks=%d", tc[ti]));
 564
 565     tlist->MakeTracks();
 566
 567     Bool_t good_cont = ti < 1;
 568     if (AliEveTrackCounter::IsActive())
 569     {
 570       AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
 571     }
 572     else
 573     {
 574       if ( ! good_cont)
 575         tlist->SetLineStyle(6);
 576     }
 577   }
 578   cont->SetTitle(Form("N all tracks = %d", count));
 579   // ??? The following does not always work:
 580   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
 581
 582   gEve->Redraw3D();
 583
 584   return cont;
 585 }