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 TString esd_track_flags(AliESDtrack* t)
  11 {
  12   TString s;
  13   Int_t   o;
  14   s += "det(in,out,refit,pid):\n";
  15   o  = AliESDtrack::kITSin;
  16   s += Form("its(%d,%d,%d,%d) ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  17   o  = AliESDtrack::kTPCin;
  18   s += Form("tpc(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  19   o  = AliESDtrack::kTRDin;
  20   s += Form("trd(%d,%d,%d,%d) ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  21   o  = AliESDtrack::kTOFin;
  22   s += Form("tof(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
  23   o  = AliESDtrack::kHMPIDout;
  24   s += Form("hmpid(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
  25   s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
  26   return s;
  27 }
  28
  29 TEveTrack* esd_make_track(TEveTrackPropagator*   trkProp,
  30                           Int_t                  index,
  31                           AliESDtrack*           at,
  32                           AliExternalTrackParam* tp=0)
  33 {
  34   // Helper function
  35   Double_t      pbuf[3], vbuf[3];
  36   TEveRecTrack  rt;
  37
  38   if(tp == 0) tp = at;
  39
  40   rt.fLabel  = at->GetLabel();
  41   rt.fIndex  = index;
  42   rt.fStatus = (Int_t) at->GetStatus();
  43   rt.fSign   = tp->GetSign();
  44   tp->GetXYZ(vbuf);
  45   rt.fV.Set(vbuf);
  46   tp->GetPxPyPz(pbuf);
  47   rt.fP.Set(pbuf);
  48   Double_t ep = at->GetP(), mc = at->GetMass();
  49   rt.fBeta = ep/TMath::Sqrt(ep*ep + mc*mc);
  50
  51   TEveTrack* track = new TEveTrack(&rt, trkProp);
  52   //PH The line below is replaced waiting for a fix in Root
  53   //PH which permits to use variable siza arguments in CINT
  54   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
  55   //PH  track->SetName(Form("ESDTrack %d", rt.fLabel));
  56   //PH  track->SetTitle(Form("pT=%.3f, pZ=%.3f; V=(%.3f, %.3f, %.3f)",
  57   //PH                 rt.fSign*TMath::Hypot(rt.fP.fX, rt.fP.fY), rt.fP.fZ,
  58   //PH                 rt.fV.fX, rt.fV.fY, rt.fV.fZ));
  59   char form[1000];
  60   sprintf(form,"TEveTrack %d", rt.fIndex);
  61   track->SetName(form);
  62   track->SetStdTitle();
  63   track->SetElementTitle(Form("%s\n%s", track->GetElementTitle(), esd_track_flags(at).Data()));
  64   return track;
  65 }
  66
  67 void esd_track_add_param(TEveTrack* track, AliExternalTrackParam* tp)
  68 {
  69   if (tp == 0)
  70     return;
  71
  72   Double_t      pbuf[3], vbuf[3];
  73   tp->GetXYZ(vbuf);
  74   tp->GetPxPyPz(pbuf);
  75
  76   TEvePathMark pm(TEvePathMark::kReference);
  77   pm.fV.Set(vbuf);
  78   pm.fP.Set(pbuf);
  79   track->AddPathMark(pm);
  80 }
  81
  82 // Use inner-tpc track params when its refit failed.
  83 Bool_t gkFixFailedITSExtr    = kFALSE;
  84
  85 // Also show lines as generated by AliESDtrack.
  86 Bool_t gkMakeTrackParamLines = kFALSE;
  87
  88 TEveTrackList* esd_tracks(Double_t min_pt=0, Double_t max_pt=10000)
  89 {
  90   AliESDEvent* esd = AliEveEventManager::AssertESD();
  91
  92   Double_t minptsq = min_pt*min_pt;
  93   Double_t maxptsq = max_pt*max_pt;
  94   Double_t ptsq;
  95
  96   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
  97   cont->SetMainColor(6);
  98   TEveTrackPropagator* trkProp = cont->GetPropagator();
  99   trkProp->SetMagField( 0.1*esd->GetMagneticField() );
 100   trkProp->SetMaxR    ( 520 );
 101
 102   gEve->AddElement(cont);
 103
 104   TEveElementList* contLines = 0;
 105   if (gkMakeTrackParamLines)
 106   {
 107     contLines = new TEveElementList("MyTracks");
 108     gEve->AddElement(contLines);
 109   }
 110
 111   Int_t    count = 0;
 112   Double_t pbuf[3];
 113   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
 114   {
 115     AliESDtrack* at = esd->GetTrack(n);
 116
 117     // Here would be sweet to have TObjectFormula.
 118     at->GetPxPyPz(pbuf);
 119     ptsq = pbuf[0]*pbuf[0] + pbuf[1]*pbuf[1];
 120     if(ptsq < minptsq || ptsq > maxptsq)
 121       continue;
 122
 123     ++count;
 124
 125     // If gkFixFailedITSExtr is TRUE (FALSE by default) and
 126     // if ITS refit failed, take track parameters at inner TPC radius.
 127     Bool_t innerTaken = kFALSE;
 128     AliExternalTrackParam* tp = at;
 129     if (gkFixFailedITSExtr && !at->IsOn(AliESDtrack::kITSrefit)) {
 130       tp = at->GetInnerParam();
 131       innerTaken = kTRUE;
 132     }
 133
 134     TEveTrack* track = esd_make_track(trkProp, n, at, tp);
 135     track->SetAttLineAttMarker(cont);
 136
 137     if (!innerTaken) {
 138       esd_track_add_param(track, at->GetInnerParam());
 139     }
 140     // esd_track_add_param(track, at->GetOuterParam());
 141
 142     gEve->AddElement(track, cont);
 143
 144     if (gkMakeTrackParamLines) {
 145       TEveLine* l = new TEveLine;
 146       l->SetName(Form("Track%d", count));
 147       l->SetLineColor(5);
 148       at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
 149       contLines->AddElement(l);
 150     }
 151   }
 152
 153   //PH The line below is replaced waiting for a fix in Root
 154   //PH which permits to use variable siza arguments in CINT
 155   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 156   //PH  const Text_t* tooltip = Form("pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
 157   char tooltip[1000];
 158   sprintf(tooltip,"pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
 159   cont->SetTitle(tooltip); // Not broadcasted automatically ...
 160
 161   cont->MakeTracks();
 162
 163   gEve->Redraw3D();
 164
 165   return cont;
 166 }
 167
 168 /******************************************************************************/
 169 // esd_tracks_from_array()
 170 /******************************************************************************/
 171
 172 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
 173 {
 174   // Retrieves AliESDTrack's from collection.
 175   // See example usage with AliAnalysisTrackCuts in the next function.
 176
 177   if (esd == 0) esd = AliEveEventManager::AssertESD();
 178
 179   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 180   cont->SetMainColor(6);
 181   TEveTrackPropagator* trkProp = cont->GetPropagator();
 182   trkProp->SetMagField( 0.1*esd->GetMagneticField() );
 183   trkProp->SetMaxR    ( 520 );
 184
 185   gEve->AddElement(cont);
 186
 187   Int_t    count = 0;
 188   TIter    next(col);
 189   TObject *obj;
 190   while ((obj = next()) != 0)
 191   {
 192     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE) {
 193       Warning("Object '%s', '%s' is not an AliESDtrack.",
 194               obj->GetName(), obj->GetTitle());
 195       continue;
 196     }
 197
 198     ++count;
 199     AliESDtrack* at = (AliESDtrack*) obj;
 200
 201     TEveTrack* track = esd_make_track(trkProp, count, at);
 202     track->SetAttLineAttMarker(cont);
 203     gEve->AddElement(track, cont);
 204   }
 205
 206   //PH The line below is replaced waiting for a fix in Root
 207   //PH which permits to use variable siza arguments in CINT
 208   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 209   //PH  const Text_t* tooltip = Form("N=%d", count);
 210   const tooltip[1000];
 211   sprintf(tooltip,"N=%d", count);
 212   cont->SetTitle(tooltip); // Not broadcasted automatically ...
 213
 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_vertex_cut
 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     printf("Estimated b resolution lower or equal zero!\n");
 248     bCov[0]=0; bCov[2]=0;
 249   }
 250   bRes[0] = TMath::Sqrt(bCov[0]);
 251   bRes[1] = TMath::Sqrt(bCov[2]);
 252
 253   // -----------------------------------
 254   // How to get to a n-sigma cut?
 255   //
 256   // The accumulated statistics from 0 to d is
 257   //
 258   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 259   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 260   //
 261   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 262   // Can this be expressed in a different way?
 263
 264   if (bRes[0] == 0 || bRes[1] ==0)
 265     return -1;
 266
 267   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 268
 269   // stupid rounding problem screws up everything:
 270   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 271   if (TMath::Exp(-d * d / 2) < 1e-10)
 272     return 1000;
 273
 274   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 275   return d;
 276 }
 277
 278 TEveElementList* esd_tracks_vertex_cut()
 279 {
 280   // Import ESD tracks, separate them into five containers according to
 281   // primary-vertex cut and ITS refit status.
 282
 283   AliESDEvent* esd = AliEveEventManager::AssertESD();
 284
 285   TEveElementList* cont = new TEveElementList("ESD Tracks");
 286
 287   gEve->AddElement(cont);
 288   TEveTrackList *tl[5];
 289   Int_t            tc[5];
 290   Int_t            count = 0;
 291
 292   tl[0] = new TEveTrackList("Sigma < 3");
 293   tc[0] = 0;
 294   tl[0]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 295   tl[0]->GetPropagator()->SetMaxR    ( 520 );
 296   tl[0]->SetMainColor(3);
 297   gEve->AddElement(tl[0], cont);
 298
 299   tl[1] = new TEveTrackList("3 < Sigma < 5");
 300   tc[1] = 0;
 301   tl[1]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 302   tl[1]->GetPropagator()->SetMaxR    ( 520 );
 303   tl[1]->SetMainColor(7);
 304   gEve->AddElement(tl[1], cont);
 305
 306   tl[2] = new TEveTrackList("5 < Sigma");
 307   tc[2] = 0;
 308   tl[2]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 309   tl[2]->GetPropagator()->SetMaxR    ( 520 );
 310   tl[2]->SetMainColor(46);
 311   gEve->AddElement(tl[2], cont);
 312
 313   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
 314   tc[3] = 0;
 315   tl[3]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 316   tl[3]->GetPropagator()->SetMaxR    ( 520 );
 317   tl[3]->SetMainColor(41);
 318   gEve->AddElement(tl[3], cont);
 319
 320   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
 321   tc[4] = 0;
 322   tl[4]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 323   tl[4]->GetPropagator()->SetMaxR    ( 520 );
 324   tl[4]->SetMainColor(48);
 325   gEve->AddElement(tl[4], cont);
 326
 327   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
 328   {
 329     AliESDtrack* at = esd->GetTrack(n);
 330
 331     Float_t s  = get_sigma_to_vertex(at);
 332     Int_t   ti;
 333     if      (s <  3) ti = 0;
 334     else if (s <= 5) ti = 1;
 335     else             ti = 2;
 336
 337     AliExternalTrackParam* tp = at;
 338     // If ITS refit failed, optionally take track parameters at inner
 339     // TPC radius and put track in a special container.
 340     // This ignores state of gkFixFailedITSExtr (used in esd_tracks()).
 341     // Use BOTH functions to compare results.
 342     if (!at->IsOn(AliESDtrack::kITSrefit)) {
 343       // tp = at->GetInnerParam();
 344       ti = (ti == 2) ? 4 : 3;
 345     }
 346
 347     TEveTrackList* tlist = tl[ti];
 348     ++tc[ti];
 349     ++count;
 350
 351     TEveTrack* track = esd_make_track(tlist->GetPropagator(), n, at, tp);
 352     track->SetAttLineAttMarker(tlist);
 353
 354     //PH The line below is replaced waiting for a fix in Root
 355     //PH which permits to use variable siza arguments in CINT
 356     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 357     //PH    track->SetName(Form("track %d, sigma=%5.3f", at->GetLabel(), s));
 358     char form[1000];
 359     sprintf(form,"TEveTrack idx=%d, sigma=%5.3f", at->GetID(), s);
 360     track->SetName(form);
 361     gEve->AddElement(track, tlist);
 362   }
 363
 364   for (Int_t ti=0; ti<5; ++ti) {
 365     TEveTrackList* tlist = tl[ti];
 366     //PH The line below is replaced waiting for a fix in Root
 367     //PH which permits to use variable siza arguments in CINT
 368     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 369     //PH    const Text_t* tooltip = Form("N tracks=%d", tc[ti]);
 370     //MT Modified somewhat.
 371     char buff[1000];
 372     sprintf(buff, "%s [%d]", tlist->GetName(), tlist->NumChildren());
 373     tlist->SetName(buff);
 374     sprintf(buff, "N tracks=%d", tc[ti]);
 375     tlist->SetTitle(buff); // Not broadcasted automatically ...
 376
 377     tlist->MakeTracks();
 378   }
 379   //PH The line below is replaced waiting for a fix in Root
 380   //PH which permits to use variable siza arguments in CINT
 381   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 382   //PH  cont->SetTitle(Form("N all tracks = %d", count));
 383   char form[1000];
 384   sprintf(form,"N all tracks = %d", count);
 385   cont->SetTitle(form);
 386   gEve->Redraw3D();
 387
 388   return cont;
 389 }