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 "TMath.h"
  14 #include "TGListTree.h"
  15 #include "TEveVSDStructs.h"
  16 #include "TEveManager.h"
  17 #include "TEveTrackPropagator.h"
  18
  19 #include "AliESDEvent.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliESDtrack.h"
  22 #include "AliExternalTrackParam.h"
  23
  24 #include "EVE/EveBase/AliEveTrack.h"
  25 #include "EVE/EveBase/AliEveMagField.h"
  26 #include "EVE/EveBase/AliEveEventManager.h"
  27
  28 #endif
  29
  30
  31 // Use inner-tpc track params when its refit failed.
  32 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
  33
  34 // Use magnetic-field as retrieved from GRP.
  35 Bool_t g_esd_tracks_true_field = kTRUE;
  36
  37 // Use Runge-Kutta track stepper.
  38 Bool_t g_esd_tracks_rk_stepper = kFALSE;
  39
  40
  41 //==============================================================================
  42
  43 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
  44                                 Float_t magF, Float_t maxR)
  45 {
  46   if (g_esd_tracks_true_field)
  47   {
  48     trkProp->SetMagFieldObj(new AliEveMagField);
  49   }
  50   else
  51   {
  52     trkProp->SetMagField(magF);
  53   }
  54   if (g_esd_tracks_rk_stepper)
  55   {
  56     trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
  57   }
  58   trkProp->SetMaxR(maxR);
  59 }
  60
  61 //==============================================================================
  62
  63 TString esd_track_title(AliESDtrack* t)
  64 {
  65   TString s;
  66
  67   Int_t label = t->GetLabel(), index = t->GetID();
  68   TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
  69   TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));
  70
  71   Double_t p[3], v[3];
  72   t->GetXYZ(v);
  73   t->GetPxPyPz(p);
  74   Double_t pt    = t->Pt();
  75   Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
  76   Double_t ptsq  = pt*pt;
  77   Double_t ptm   = pt / (1.0 + pt*ptsig);
  78   Double_t ptM   = pt / (1.0 - pt*ptsig);
  79
  80   s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
  81            "pT = %.3f + %.3f - %.3f [%.3f]\n"
  82            "P  = (%.3f, %.3f, %.3f)\n"
  83            "V  = (%.3f, %.3f, %.3f)\n",
  84            idx.Data(), lbl.Data(), t->Charge(), 0,
  85            pt, ptM - pt, pt - ptm, ptsig*ptsq,
  86            p[0], p[1], p[2],
  87            v[0], v[1], v[2]);
  88
  89   Int_t   o;
  90   s += "Det (in,out,refit,pid):\n";
  91   o  = AliESDtrack::kITSin;
  92   s += Form("ITS (%d,%d,%d,%d)  ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  93   o  = AliESDtrack::kTPCin;
  94   s += Form("TPC(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  95   o  = AliESDtrack::kTRDin;
  96   s += Form("TRD(%d,%d,%d,%d) ",    t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  97   o  = AliESDtrack::kTOFin;
  98   s += Form("TOF(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  99   o  = AliESDtrack::kHMPIDout;
 100   s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
 101   s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
 102
 103   if (t->IsOn(AliESDtrack::kESDpid))
 104   {
 105     Double_t pid[5];
 106     t->GetESDpid(pid);
 107     s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
 108   }
 109
 110   return s;
 111 }
 112
 113 //==============================================================================
 114
 115 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
 116 {
 117   // Add additional track parameters as a path-mark to track.
 118
 119   if (tp == 0)
 120     return;
 121
 122   Double_t pbuf[3], vbuf[3];
 123   tp->GetXYZ(vbuf);
 124   tp->GetPxPyPz(pbuf);
 125
 126   TEvePathMark pm(TEvePathMark::kReference);
 127   pm.fV.Set(vbuf);
 128   pm.fP.Set(pbuf);
 129   track->AddPathMark(pm);
 130 }
 131
 132 //==============================================================================
 133
 134 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
 135 {
 136   // Make a standard track representation and put it into given container.
 137
 138   // Choose which parameters to use a track's starting point.
 139   // If gkFixFailedITSExtr is TRUE (FALSE by default) and
 140   // if ITS refit failed, take track parameters at inner TPC radius.
 141
 142   const AliExternalTrackParam* tp = at;
 143
 144   Bool_t innerTaken = kFALSE;
 145   if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
 146   {
 147     tp = at->GetInnerParam();
 148     innerTaken = kTRUE;
 149   }
 150
 151   AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
 152   track->SetAttLineAttMarker(cont);
 153   track->SetName(Form("AliEveTrack %d", at->GetID()));
 154   track->SetElementTitle(esd_track_title(at));
 155   track->SetSourceObject(at);
 156
 157   // Add inner/outer track parameters as path-marks.
 158   if (at->IsOn(AliESDtrack::kTPCrefit))
 159   {
 160     if ( ! innerTaken)
 161     {
 162       esd_track_add_param(track, at->GetInnerParam());
 163     }
 164     esd_track_add_param(track, at->GetOuterParam());
 165   }
 166
 167   return track;
 168 }
 169
 170
 171 //==============================================================================
 172 // esd_tracks()
 173 //==============================================================================
 174
 175 TEveTrackList* esd_tracks()
 176 {
 177   AliESDEvent* esd = AliEveEventManager::AssertESD();
 178
 179   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 180   cont->SetMainColor(6);
 181
 182   esd_track_propagator_setup(cont->GetPropagator(),
 183                              0.1*esd->GetMagneticField(), 520);
 184
 185   gEve->AddElement(cont);
 186
 187   Int_t count = 0;
 188   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 189   {
 190     ++count;
 191     AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
 192
 193     cont->AddElement(track);
 194   }
 195   cont->SetTitle(Form("N=%d", count));
 196   cont->MakeTracks();
 197
 198   gEve->Redraw3D();
 199
 200   return cont;
 201 }
 202
 203 TEveElementList* esd_tracks_MI()
 204 {
 205   AliESDEvent* esd = AliEveEventManager::AssertESD();
 206
 207   TEveElementList* cont = new TEveElementList("ESD Tracks MI");
 208   gEve->AddElement(cont);
 209
 210   Int_t count = 0;
 211   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 212   {
 213     ++count;
 214     AliESDtrack* at = esd->GetTrack(n);
 215     TEveLine* l = new TEveLine;
 216     l->SetLineColor(5);
 217     at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
 218     l->SetName(Form("ESDTrackMI %d", at->GetID()));
 219     l->SetElementTitle(esd_track_title(at));
 220     l->SetSourceObject(at);
 221
 222     cont->AddElement(l);
 223   }
 224   cont->SetTitle(Form("N=%d", count));
 225
 226   gEve->Redraw3D();
 227
 228   return cont;
 229 }
 230
 231
 232 //==============================================================================
 233 // esd_tracks_from_array()
 234 //==============================================================================
 235
 236 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
 237 {
 238   // Retrieves AliESDTrack's from collection.
 239   // See example usage with AliAnalysisTrackCuts in the next function.
 240
 241   if (esd == 0) esd = AliEveEventManager::AssertESD();
 242
 243   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 244   cont->SetMainColor(6);
 245
 246   esd_track_propagator_setup(cont->GetPropagator(),
 247                              0.1*esd->GetMagneticField(), 520);
 248
 249   gEve->AddElement(cont);
 250
 251   Int_t    count = 0;
 252   TIter    next(col);
 253   TObject *obj;
 254   while ((obj = next()) != 0)
 255   {
 256     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
 257     {
 258       Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
 259               obj->GetName(), obj->GetTitle());
 260       continue;
 261     }
 262
 263     ++count;
 264     AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
 265
 266     AliEveTrack* track = esd_make_track(at, cont);
 267     cont->AddElement(track);
 268   }
 269   cont->SetTitle(Form("N=%d", count));
 270   cont->MakeTracks();
 271
 272   gEve->Redraw3D();
 273
 274   return cont;
 275 }
 276
 277 void esd_tracks_alianalcuts_demo()
 278 {
 279   AliESDEvent* esd = AliEveEventManager::AssertESD();
 280
 281   AliESDtrackCuts atc;
 282   atc.SetPtRange(0.1, 5);
 283   atc.SetRapRange(-1, 1);
 284
 285   esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
 286 }
 287
 288
 289 //==============================================================================
 290 // esd_tracks_by_category
 291 //==============================================================================
 292
 293 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
 294 {
 295   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
 296   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 297
 298   Float_t b[2];
 299   Float_t bRes[2];
 300   Float_t bCov[3];
 301   esdTrack->GetImpactParameters(b,bCov);
 302   if (bCov[0] <= 0 || bCov[2] <= 0)
 303   {
 304     printf("Estimated b resolution lower or equal zero!\n");
 305     bCov[0] = bCov[2] = 0;
 306   }
 307   bRes[0] = TMath::Sqrt(bCov[0]);
 308   bRes[1] = TMath::Sqrt(bCov[2]);
 309
 310   // -----------------------------------
 311   // How to get to a n-sigma cut?
 312   //
 313   // The accumulated statistics from 0 to d is
 314   //
 315   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 316   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 317   //
 318   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 319   // Can this be expressed in a different way?
 320
 321   if (bRes[0] == 0 || bRes[1] == 0)
 322     return -1;
 323
 324   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 325
 326   // stupid rounding problem screws up everything:
 327   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 328   if (TMath::Exp(-d * d / 2) < 1e-10)
 329     return 1000;
 330
 331   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 332   return d;
 333 }
 334
 335 TEveElementList* g_esd_tracks_by_category_container = 0;
 336
 337 TEveElementList* esd_tracks_by_category()
 338 {
 339   // Import ESD tracks, separate them into several containers
 340   // according to primary-vertex cut and ITS&TPC refit status.
 341
 342   AliESDEvent* esd = AliEveEventManager::AssertESD();
 343
 344   TEveElementList* cont = new TEveElementList("ESD Tracks by category");
 345   gEve->AddElement(cont);
 346
 347   g_esd_tracks_by_category_container = cont;
 348
 349   const Int_t   nCont = 7;
 350   const Float_t maxR  = 520;
 351   const Float_t magF  = 0.1*esd->GetMagneticField();
 352
 353   TEveTrackList *tl[nCont];
 354   Int_t          tc[nCont];
 355   Int_t          count = 0;
 356
 357   tl[0] = new TEveTrackList("Sigma < 3");
 358   tc[0] = 0;
 359   esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
 360   tl[0]->SetMainColor(3);
 361   cont->AddElement(tl[0]);
 362
 363   tl[1] = new TEveTrackList("3 < Sigma < 5");
 364   tc[1] = 0;
 365   esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
 366   tl[1]->SetMainColor(7);
 367   cont->AddElement(tl[1]);
 368
 369   tl[2] = new TEveTrackList("5 < Sigma");
 370   tc[2] = 0;
 371   esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
 372   tl[2]->SetMainColor(46);
 373   cont->AddElement(tl[2]);
 374
 375   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
 376   tc[3] = 0;
 377   esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
 378   tl[3]->SetMainColor(41);
 379   cont->AddElement(tl[3]);
 380
 381   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
 382   tc[4] = 0;
 383   esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
 384   tl[4]->SetMainColor(48);
 385   cont->AddElement(tl[4]);
 386
 387   tl[5] = new TEveTrackList("no TPC refit");
 388   tc[5] = 0;
 389   esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
 390   tl[5]->SetMainColor(kRed);
 391   cont->AddElement(tl[5]);
 392
 393   tl[6] = new TEveTrackList("ITS stand-alone");
 394   tc[6] = 0;
 395   esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
 396   tl[6]->SetMainColor(kMagenta - 9);
 397   cont->AddElement(tl[6]);
 398
 399   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 400   {
 401     AliESDtrack* at = esd->GetTrack(n);
 402
 403     Float_t s  = get_sigma_to_vertex(at);
 404     Int_t   ti;
 405     if      (s <  3) ti = 0;
 406     else if (s <= 5) ti = 1;
 407     else             ti = 2;
 408
 409     if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
 410     {
 411       ti = 6;
 412     }
 413     else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
 414     {
 415       ti = 5;
 416     }
 417     else if ( ! at->IsOn(AliESDtrack::kITSrefit))
 418     {
 419       ti = (ti == 2) ? 4 : 3;
 420     }
 421
 422     TEveTrackList* tlist = tl[ti];
 423     ++tc[ti];
 424     ++count;
 425
 426     AliEveTrack* track = esd_make_track(at, tlist);
 427     track->SetName(Form("AliEveTrack idx=%d, sigma=%5.3f", at->GetID(), s));
 428     tlist->AddElement(track);
 429   }
 430
 431   for (Int_t ti = 0; ti < nCont; ++ti)
 432   {
 433     TEveTrackList* tlist = tl[ti];
 434     tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
 435     tlist->SetTitle(Form("N tracks=%d", tc[ti]));
 436
 437     tlist->MakeTracks();
 438   }
 439   cont->SetTitle(Form("N all tracks = %d", count));
 440   // ??? The following does not always work:
 441   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
 442
 443   gEve->Redraw3D();
 444
 445   return cont;
 446 }