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