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 TEveTrack* esd_make_track(TEveTrackPropagator*   rnrStyle,
  11                             Int_t                  index,
  12                             AliESDtrack*           at,
  13                             AliExternalTrackParam* tp=0)
  14 {
  15   // Helper function
  16   Double_t      pbuf[3], vbuf[3];
  17   TEveRecTrack  rt;
  18
  19   if(tp == 0) tp = at;
  20
  21   rt.fLabel  = at->GetLabel();
  22   rt.fIndex  = index;
  23   rt.fStatus = (Int_t) at->GetStatus();
  24   rt.fSign   = tp->GetSign();
  25   tp->GetXYZ(vbuf);
  26   rt.fV.Set(vbuf);
  27   tp->GetPxPyPz(pbuf);
  28   rt.fP.Set(pbuf);
  29   Double_t ep = at->GetP(), mc = at->GetMass();
  30   rt.fBeta = ep/TMath::Sqrt(ep*ep + mc*mc);
  31
  32   TEveTrack* track = new TEveTrack(&rt, rnrStyle);
  33   //PH The line below is replaced waiting for a fix in Root
  34   //PH which permits to use variable siza arguments in CINT
  35   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
  36   //PH  track->SetName(Form("ESDTrack %d", rt.fLabel));
  37   //PH  track->SetTitle(Form("pT=%.3f, pZ=%.3f; V=(%.3f, %.3f, %.3f)",
  38   //PH                 rt.fSign*TMath::Hypot(rt.fP.fX, rt.fP.fY), rt.fP.fZ,
  39   //PH                 rt.fV.fX, rt.fV.fY, rt.fV.fZ));
  40   char form[1000];
  41   sprintf(form,"TEveTrack %d", rt.fIndex);
  42   track->SetName(form);
  43   track->SetStdTitle();
  44   return track;
  45 }
  46
  47 Bool_t gkFixFailedITSExtr = kFALSE;
  48
  49 TEveTrackList* esd_tracks(Double_t min_pt=0.1, Double_t max_pt=100)
  50 {
  51   AliESDEvent* esd = AliEveEventManager::AssertESD();
  52
  53   Double_t minptsq = min_pt*min_pt;
  54   Double_t maxptsq = max_pt*max_pt;
  55   Double_t ptsq;
  56
  57   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
  58   cont->SetMainColor(Color_t(6));
  59   TEveTrackPropagator* rnrStyle = cont->GetPropagator();
  60   rnrStyle->SetMagField( 0.1*esd->GetMagneticField() );
  61
  62   gEve->AddElement(cont);
  63
  64   Int_t    count = 0;
  65   Double_t pbuf[3];
  66   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
  67   {
  68     AliESDtrack* at = esd->GetTrack(n);
  69
  70     // Here would be sweet to have TObjectFormula.
  71     at->GetPxPyPz(pbuf);
  72     ptsq = pbuf[0]*pbuf[0] + pbuf[1]*pbuf[1];
  73     if(ptsq < minptsq || ptsq > maxptsq)
  74       continue;
  75
  76     ++count;
  77
  78     // If gkFixFailedITSExtr is TRUE (FALSE by default) and
  79     // if ITS refit failed, take track parameters at inner TPC radius.
  80     AliExternalTrackParam* tp = at;
  81     if (gkFixFailedITSExtr && !at->IsOn(AliESDtrack::kITSrefit)) {
  82        tp = at->GetInnerParam();
  83     }
  84
  85     TEveTrack* track = esd_make_track(rnrStyle, n, at, tp);
  86     track->SetAttLineAttMarker(cont);
  87     gEve->AddElement(track, cont);
  88   }
  89
  90   //PH The line below is replaced waiting for a fix in Root
  91   //PH which permits to use variable siza arguments in CINT
  92   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
  93   //PH  const Text_t* tooltip = Form("pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
  94   char tooltip[1000];
  95   sprintf(tooltip,"pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
  96   cont->SetTitle(tooltip); // Not broadcasted automatically ...
  97   cont->UpdateItems();
  98
  99   cont->MakeTracks();
 100
 101   gEve->Redraw3D();
 102
 103   return cont;
 104 }
 105
 106 /******************************************************************************/
 107 // esd_tracks_from_array()
 108 /******************************************************************************/
 109
 110 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
 111 {
 112   // Retrieves AliESDTrack's from collection.
 113   // See example usage with AliAnalysisTrackCuts in the next function.
 114
 115   if (esd == 0) esd = AliEveEventManager::AssertESD();
 116
 117   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
 118   cont->SetMainColor(Color_t(6));
 119   TEveTrackPropagator* rnrStyle = cont->GetPropagator();
 120   rnrStyle->SetMagField( 0.1*esd->GetMagneticField() );
 121
 122   gEve->AddElement(cont);
 123
 124   Int_t    count = 0;
 125   TIter    next(col);
 126   TObject *obj;
 127   while ((obj = next()) != 0)
 128   {
 129     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE) {
 130       Warning("Object '%s', '%s' is not an AliESDtrack.",
 131               obj->GetName(), obj->GetTitle());
 132       continue;
 133     }
 134
 135     ++count;
 136     AliESDtrack* at = (AliESDtrack*) obj;
 137
 138     TEveTrack* track = esd_make_track(rnrStyle, count, at);
 139     track->SetAttLineAttMarker(cont);
 140     gEve->AddElement(track, cont);
 141   }
 142
 143   //PH The line below is replaced waiting for a fix in Root
 144   //PH which permits to use variable siza arguments in CINT
 145   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 146   //PH  const Text_t* tooltip = Form("N=%d", count);
 147   const tooltip[1000];
 148   sprintf(tooltip,"N=%d", count);
 149   cont->SetTitle(tooltip); // Not broadcasted automatically ...
 150   cont->UpdateItems();
 151
 152   cont->MakeTracks();
 153
 154   gEve->Redraw3D();
 155
 156   return cont;
 157 }
 158
 159 void esd_tracks_alianalcuts_demo()
 160 {
 161   AliESDEvent* esd = AliEveEventManager::AssertESD();
 162   gSystem->Load("libANALYSIS");
 163
 164   AliAnalysisTrackCuts atc;
 165   atc.SetPtRange(0.1, 5);
 166   atc.SetRapRange(-1, 1);
 167
 168   esd_tracks_from_array(atc.GetAcceptedParticles(esd), esd);
 169 }
 170
 171 /******************************************************************************/
 172 // esd_tracks_vertex_cut
 173 /******************************************************************************/
 174
 175 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
 176 {
 177   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
 178   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 179
 180   Float_t b[2];
 181   Float_t bRes[2];
 182   Float_t bCov[3];
 183   esdTrack->GetImpactParameters(b,bCov);
 184   if (bCov[0]<=0 || bCov[2]<=0) {
 185     printf("Estimated b resolution lower or equal zero!\n");
 186     bCov[0]=0; bCov[2]=0;
 187   }
 188   bRes[0] = TMath::Sqrt(bCov[0]);
 189   bRes[1] = TMath::Sqrt(bCov[2]);
 190
 191   // -----------------------------------
 192   // How to get to a n-sigma cut?
 193   //
 194   // The accumulated statistics from 0 to d is
 195   //
 196   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 197   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 198   //
 199   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 200   // Can this be expressed in a different way?
 201
 202   if (bRes[0] == 0 || bRes[1] ==0)
 203     return -1;
 204
 205   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 206
 207   // stupid rounding problem screws up everything:
 208   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 209   if (TMath::Exp(-d * d / 2) < 1e-10)
 210     return 1000;
 211
 212   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 213   return d;
 214 }
 215
 216 TEveElementList* esd_tracks_vertex_cut()
 217 {
 218   // Import ESD tracks, separate them into five containers according to
 219   // primary-vertex cut and ITS refit status.
 220
 221   AliESDEvent* esd = AliEveEventManager::AssertESD();
 222
 223   TEveElementList* cont = new TEveElementList("ESD Tracks", 0, kTRUE);
 224   gEve->AddElement(cont);
 225   TEveTrackList *tl[5];
 226   Int_t            tc[5];
 227   Int_t            count = 0;
 228
 229   tl[0] = new TEveTrackList("Sigma < 3");
 230   tc[0] = 0;
 231   tl[0]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 232   tl[0]->SetMainColor(Color_t(3));
 233   gEve->AddElement(tl[0], cont);
 234
 235   tl[1] = new TEveTrackList("3 < Sigma < 5");
 236   tc[1] = 0;
 237   tl[1]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 238   tl[1]->SetMainColor(Color_t(7));
 239   gEve->AddElement(tl[1], cont);
 240
 241   tl[2] = new TEveTrackList("5 < Sigma");
 242   tc[2] = 0;
 243   tl[2]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 244   tl[2]->SetMainColor(Color_t(46));
 245   gEve->AddElement(tl[2], cont);
 246
 247   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
 248   tc[3] = 0;
 249   tl[3]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 250   tl[3]->SetMainColor(Color_t(41));
 251   gEve->AddElement(tl[3], cont);
 252
 253   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
 254   tc[4] = 0;
 255   tl[4]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
 256   tl[4]->SetMainColor(Color_t(48));
 257   gEve->AddElement(tl[4], cont);
 258
 259   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
 260   {
 261     AliESDtrack* at = esd->GetTrack(n);
 262
 263     Float_t s  = get_sigma_to_vertex(at);
 264     Int_t   ti;
 265     if      (s <  3) ti = 0;
 266     else if (s <= 5) ti = 1;
 267     else             ti = 2;
 268
 269     AliExternalTrackParam* tp = at;
 270     // If ITS refit failed, optionally take track parameters at inner
 271     // TPC radius and put track in a special container.
 272     // This ignores state of gkFixFailedITSExtr (used in esd_tracks()).
 273     // Use BOTH functions to compare results.
 274     if (!at->IsOn(AliESDtrack::kITSrefit)) {
 275       // tp = at->GetInnerParam();
 276       ti = (ti == 2) ? 4 : 3;
 277     }
 278
 279     TEveTrackList* tlist = tl[ti];
 280     ++tc[ti];
 281     ++count;
 282
 283     TEveTrack* track = esd_make_track(tlist->GetPropagator(), n, at, tp);
 284     track->SetAttLineAttMarker(tlist);
 285
 286     //PH The line below is replaced waiting for a fix in Root
 287     //PH which permits to use variable siza arguments in CINT
 288     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 289     //PH    track->SetName(Form("track %d, sigma=%5.3f", at->GetLabel(), s));
 290     char form[1000];
 291     sprintf(form,"TEveTrack idx=%d, sigma=%5.3f", at->GetID(), s);
 292     track->SetName(form);
 293     gEve->AddElement(track, tlist);
 294   }
 295
 296   for (Int_t ti=0; ti<5; ++ti) {
 297     TEveTrackList* tlist = tl[ti];
 298     //PH The line below is replaced waiting for a fix in Root
 299     //PH which permits to use variable siza arguments in CINT
 300     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 301     //PH    const Text_t* tooltip = Form("N tracks=%d", tc[ti]);
 302     //MT Modified somewhat.
 303     char buff[1000];
 304     sprintf(buff, "%s [%d]", tlist->GetName(), tlist->GetNChildren());
 305     tlist->SetName(buff);
 306     sprintf(buff, "N tracks=%d", tc[ti]);
 307     tlist->SetTitle(buff); // Not broadcasted automatically ...
 308     tlist->UpdateItems();
 309
 310     tlist->MakeTracks();
 311   }
 312   //PH The line below is replaced waiting for a fix in Root
 313   //PH which permits to use variable siza arguments in CINT
 314   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 315   //PH  cont->SetTitle(Form("N all tracks = %d", count));
 316   char form[1000];
 317   sprintf(form,"N all tracks = %d", count);
 318   cont->SetTitle(form);
 319   cont->UpdateItems();
 320   gEve->Redraw3D();
 321
 322   return cont;
 323 }