EVE/alice-macros/esd_tracks.C

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