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