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