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