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(TEveTrack* 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   TEveTrack* track = new TEveTrack(&rt, cont->GetPropagator());
 102   track->SetAttLineAttMarker(cont);
 103   track->SetName(Form("TEveTrack %d", at->GetID()));
 104   track->SetElementTitle(esd_track_title(at));
 105
 106   // Add inner/outer track parameters as path-marks.
 107   if (at->IsOn(AliESDtrack::kTPCrefit))
 108   {
 109     if ( ! innerTaken)
 110     {
 111       esd_track_add_param(track, at->GetInnerParam());
 112     }
 113     esd_track_add_param(track, at->GetOuterParam());
 114   }
 115
 116   return track;
 117 }
 118
 119 TEveTrackList* esd_tracks()
 120 {
 121   AliESDEvent* esd = AliEveEventManager::AssertESD();
 122
 123   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 124   cont->SetMainColor(6);
 125   TEveTrackPropagator* trkProp = cont->GetPropagator();
 126   trkProp->SetMagField(0.1*esd->GetMagneticField());
 127   trkProp->SetMaxR    (520);
 128
 129   gEve->AddElement(cont);
 130
 131   Int_t count = 0;
 132   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 133   {
 134     ++count;
 135     TEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
 136
 137     cont->AddElement(track);
 138   }
 139   cont->SetTitle(Form("N=%d", count));
 140   cont->MakeTracks();
 141
 142   gEve->Redraw3D();
 143
 144   return cont;
 145 }
 146
 147 TEveElementList* esd_tracks_MI()
 148 {
 149   AliESDEvent* esd = AliEveEventManager::AssertESD();
 150
 151   TEveElementList* cont = new TEveElementList("ESD Tracks MI");
 152   gEve->AddElement(cont);
 153
 154   Int_t count = 0;
 155   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 156   {
 157     ++count;
 158     AliESDtrack* at = esd->GetTrack(n);
 159     TEveLine* l = new TEveLine;
 160     l->SetLineColor(5);
 161     at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
 162     l->SetName(Form("ESDTrackMI %d", at->GetID()));
 163     l->SetElementTitle(esd_track_title(at));
 164
 165     cont->AddElement(l);
 166   }
 167   cont->SetTitle(Form("N=%d", count));
 168
 169   gEve->Redraw3D();
 170
 171   return cont;
 172 }
 173
 174 /******************************************************************************/
 175 // esd_tracks_from_array()
 176 /******************************************************************************/
 177
 178 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
 179 {
 180   // Retrieves AliESDTrack's from collection.
 181   // See example usage with AliAnalysisTrackCuts in the next function.
 182
 183   if (esd == 0) esd = AliEveEventManager::AssertESD();
 184
 185   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 186   cont->SetMainColor(6);
 187   TEveTrackPropagator* trkProp = cont->GetPropagator();
 188   trkProp->SetMagField(0.1*esd->GetMagneticField());
 189   trkProp->SetMaxR    (520);
 190
 191   gEve->AddElement(cont);
 192
 193   Int_t    count = 0;
 194   TIter    next(col);
 195   TObject *obj;
 196   while ((obj = next()) != 0)
 197   {
 198     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
 199     {
 200       Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
 201               obj->GetName(), obj->GetTitle());
 202       continue;
 203     }
 204
 205     ++count;
 206     AliESDtrack* at = (AliESDtrack*) obj;
 207
 208     TEveTrack* track = esd_make_track(at, cont);
 209     cont->AddElement(track);
 210   }
 211   cont->SetTitle(Form("N=%d", count));
 212   cont->MakeTracks();
 213
 214   gEve->Redraw3D();
 215
 216   return cont;
 217 }
 218
 219 void esd_tracks_alianalcuts_demo()
 220 {
 221   AliESDEvent* esd = AliEveEventManager::AssertESD();
 222   gSystem->Load("libANALYSIS");
 223
 224   AliAnalysisTrackCuts atc;
 225   atc.SetPtRange(0.1, 5);
 226   atc.SetRapRange(-1, 1);
 227
 228   esd_tracks_from_array(atc.GetAcceptedParticles(esd), esd);
 229 }
 230
 231 /******************************************************************************/
 232 // esd_tracks_by_category
 233 /******************************************************************************/
 234
 235 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
 236 {
 237   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
 238   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 239
 240   Float_t b[2];
 241   Float_t bRes[2];
 242   Float_t bCov[3];
 243   esdTrack->GetImpactParameters(b,bCov);
 244   if (bCov[0] <= 0 || bCov[2] <= 0)
 245   {
 246     printf("Estimated b resolution lower or equal zero!\n");
 247     bCov[0] = bCov[2] = 0;
 248   }
 249   bRes[0] = TMath::Sqrt(bCov[0]);
 250   bRes[1] = TMath::Sqrt(bCov[2]);
 251
 252   // -----------------------------------
 253   // How to get to a n-sigma cut?
 254   //
 255   // The accumulated statistics from 0 to d is
 256   //
 257   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 258   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 259   //
 260   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 261   // Can this be expressed in a different way?
 262
 263   if (bRes[0] == 0 || bRes[1] == 0)
 264     return -1;
 265
 266   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 267
 268   // stupid rounding problem screws up everything:
 269   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 270   if (TMath::Exp(-d * d / 2) < 1e-10)
 271     return 1000;
 272
 273   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 274   return d;
 275 }
 276
 277 TEveElementList* g_esd_tracks_by_category_container = 0;
 278
 279 TEveElementList* esd_tracks_by_category()
 280 {
 281   // Import ESD tracks, separate them into several containers
 282   // according to primary-vertex cut and ITS&TPC refit status.
 283
 284   AliESDEvent* esd = AliEveEventManager::AssertESD();
 285
 286   TEveElementList* cont = new TEveElementList("ESD Tracks by category");
 287   gEve->AddElement(cont);
 288
 289   g_esd_tracks_by_category_container = cont;
 290
 291   const Int_t   nCont = 7;
 292   const Float_t maxR  = 520;
 293   const Float_t magF  = 0.1*esd->GetMagneticField();
 294
 295   TEveTrackList *tl[nCont];
 296   Int_t          tc[nCont];
 297   Int_t          count = 0;
 298
 299   tl[0] = new TEveTrackList("Sigma < 3");
 300   tc[0] = 0;
 301   tl[0]->GetPropagator()->SetMagField(magF);
 302   tl[0]->GetPropagator()->SetMaxR    (maxR);
 303   tl[0]->SetMainColor(3);
 304   cont->AddElement(tl[0]);
 305
 306   tl[1] = new TEveTrackList("3 < Sigma < 5");
 307   tc[1] = 0;
 308   tl[1]->GetPropagator()->SetMagField(magF);
 309   tl[1]->GetPropagator()->SetMaxR    (maxR);
 310   tl[1]->SetMainColor(7);
 311   cont->AddElement(tl[1]);
 312
 313   tl[2] = new TEveTrackList("5 < Sigma");
 314   tc[2] = 0;
 315   tl[2]->GetPropagator()->SetMagField(magF);
 316   tl[2]->GetPropagator()->SetMaxR    (maxR);
 317   tl[2]->SetMainColor(46);
 318   cont->AddElement(tl[2]);
 319
 320   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
 321   tc[3] = 0;
 322   tl[3]->GetPropagator()->SetMagField(magF);
 323   tl[3]->GetPropagator()->SetMaxR    (maxR);
 324   tl[3]->SetMainColor(41);
 325   cont->AddElement(tl[3]);
 326
 327   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
 328   tc[4] = 0;
 329   tl[4]->GetPropagator()->SetMagField(magF);
 330   tl[4]->GetPropagator()->SetMaxR    (maxR);
 331   tl[4]->SetMainColor(48);
 332   cont->AddElement(tl[4]);
 333
 334   tl[5] = new TEveTrackList("no TPC refit");
 335   tc[5] = 0;
 336   tl[5]->GetPropagator()->SetMagField(magF);
 337   tl[5]->GetPropagator()->SetMaxR    (maxR);
 338   tl[5]->SetMainColor(kRed);
 339   cont->AddElement(tl[5]);
 340
 341   tl[6] = new TEveTrackList("ITS stand-alone");
 342   tc[6] = 0;
 343   tl[6]->GetPropagator()->SetMagField(magF);
 344   tl[6]->GetPropagator()->SetMaxR    (maxR);
 345   tl[6]->SetMainColor(kMagenta - 9);
 346   cont->AddElement(tl[6]);
 347
 348   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
 349   {
 350     AliESDtrack* at = esd->GetTrack(n);
 351
 352     Float_t s  = get_sigma_to_vertex(at);
 353     Int_t   ti;
 354     if      (s <  3) ti = 0;
 355     else if (s <= 5) ti = 1;
 356     else             ti = 2;
 357
 358     if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
 359     {
 360       ti = 6;
 361     }
 362     else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
 363     {
 364       ti = 5;
 365     }
 366     else if ( ! at->IsOn(AliESDtrack::kITSrefit))
 367     {
 368       ti = (ti == 2) ? 4 : 3;
 369     }
 370
 371     TEveTrackList* tlist = tl[ti];
 372     ++tc[ti];
 373     ++count;
 374
 375     TEveTrack* track = esd_make_track(at, tlist);
 376     track->SetName(Form("TEveTrack idx=%d, sigma=%5.3f", at->GetID(), s));
 377     tlist->AddElement(track);
 378   }
 379
 380   for (Int_t ti = 0; ti < nCont; ++ti)
 381   {
 382     TEveTrackList* tlist = tl[ti];
 383     tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
 384     tlist->SetTitle(Form("N tracks=%d", tc[ti]));
 385
 386     tlist->MakeTracks();
 387   }
 388   cont->SetTitle(Form("N all tracks = %d", count));
 389   // ??? The following does not always work:
 390   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
 391
 392   gEve->Redraw3D();
 393
 394   return cont;
 395 }