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