* esd_tracks.C
[u/mrichter/AliRoot.git] / 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 #if !defined(__CINT__) || defined(__MAKECINT__)
11
12 #include "TString.h"
13 #include "TMath.h"
14 #include "TGListTree.h"
15 #include "TEveVSDStructs.h"
16 #include "TEveManager.h"
17 #include "TEveTrackPropagator.h"
18
19 #include "AliESDEvent.h"
20 #include "AliESDtrackCuts.h"
21 #include "AliESDtrack.h"
22 #include "AliExternalTrackParam.h"
23
24 #include "EVE/EveBase/AliEveTrack.h"
25 #include "EVE/EveBase/AliEveMagField.h"
26 #include "EVE/EveBase/AliEveEventManager.h"
27
28 #endif
29
30
31 // Use inner-tpc track params when its refit failed.
32 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
33
34 // Use magnetic-field as retrieved from GRP.
35 Bool_t g_esd_tracks_true_field = kTRUE;
36
37 // Use Runge-Kutta track stepper.
38 Bool_t g_esd_tracks_rk_stepper = kFALSE;
39
40
41 //==============================================================================
42
43 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
44                                 Float_t magF, Float_t maxR)
45 {
46   if (g_esd_tracks_true_field)
47   {
48     trkProp->SetMagFieldObj(new AliEveMagField);
49   }
50   else
51   {
52     trkProp->SetMagField(magF);
53   }
54   if (g_esd_tracks_rk_stepper)
55   {
56     trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
57   }
58   trkProp->SetMaxR(maxR);
59 }
60
61 //==============================================================================
62
63 TString esd_track_title(AliESDtrack* t)
64 {
65   TString s;
66
67   Int_t label = t->GetLabel(), index = t->GetID();
68   TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
69   TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));
70
71   Double_t p[3], v[3];
72   t->GetXYZ(v);
73   t->GetPxPyPz(p);
74   Double_t pt    = t->Pt();
75   Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
76   Double_t ptsq  = pt*pt;
77   Double_t ptm   = pt / (1.0 + pt*ptsig);
78   Double_t ptM   = pt / (1.0 - pt*ptsig);
79
80   s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
81            "pT = %.3f + %.3f - %.3f [%.3f]\n"
82            "P  = (%.3f, %.3f, %.3f)\n"
83            "V  = (%.3f, %.3f, %.3f)\n",
84            idx.Data(), lbl.Data(), t->Charge(), 0,
85            pt, ptM - pt, pt - ptm, ptsig*ptsq,
86            p[0], p[1], p[2],
87            v[0], v[1], v[2]);
88
89   Int_t   o;
90   s += "Det (in,out,refit,pid):\n";
91   o  = AliESDtrack::kITSin;
92   s += Form("ITS (%d,%d,%d,%d)  ",  t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
93   o  = AliESDtrack::kTPCin;
94   s += Form("TPC(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
95   o  = AliESDtrack::kTRDin;
96   s += Form("TRD(%d,%d,%d,%d) ",    t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
97   o  = AliESDtrack::kTOFin;
98   s += Form("TOF(%d,%d,%d,%d)\n",   t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
99   o  = AliESDtrack::kHMPIDout;
100   s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
101   s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
102
103   if (t->IsOn(AliESDtrack::kESDpid))
104   {
105     Double_t pid[5];
106     t->GetESDpid(pid);
107     s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
108   }
109
110   return s;
111 }
112
113 //==============================================================================
114
115 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
116 {
117   // Add additional track parameters as a path-mark to track.
118
119   if (tp == 0)
120     return;
121
122   Double_t pbuf[3], vbuf[3];
123   tp->GetXYZ(vbuf);
124   tp->GetPxPyPz(pbuf);
125
126   TEvePathMark pm(TEvePathMark::kReference);
127   pm.fV.Set(vbuf);
128   pm.fP.Set(pbuf);
129   track->AddPathMark(pm);
130 }
131
132 //==============================================================================
133
134 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
135 {
136   // Make a standard track representation and put it into given container.
137
138   // Choose which parameters to use a track's starting point.
139   // If gkFixFailedITSExtr is TRUE (FALSE by default) and
140   // if ITS refit failed, take track parameters at inner TPC radius.
141
142   const AliExternalTrackParam* tp = at;
143
144   Bool_t innerTaken = kFALSE;
145   if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
146   {
147     tp = at->GetInnerParam();
148     innerTaken = kTRUE;
149   }
150
151   AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
152   track->SetAttLineAttMarker(cont);
153   track->SetName(Form("AliEveTrack %d", at->GetID()));
154   track->SetElementTitle(esd_track_title(at));
155   track->SetSourceObject(at);
156
157   // Add inner/outer track parameters as path-marks.
158   if (at->IsOn(AliESDtrack::kTPCrefit))
159   {
160     if ( ! innerTaken)
161     {
162       esd_track_add_param(track, at->GetInnerParam());
163     }
164     esd_track_add_param(track, at->GetOuterParam());
165   }
166
167   return track;
168 }
169
170
171 //==============================================================================
172 // esd_tracks()
173 //==============================================================================
174
175 TEveTrackList* esd_tracks()
176 {
177   AliESDEvent* esd = AliEveEventManager::AssertESD();
178
179   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
180   cont->SetMainColor(6);
181
182   esd_track_propagator_setup(cont->GetPropagator(),
183                              0.1*esd->GetMagneticField(), 520);
184
185   gEve->AddElement(cont);
186
187   Int_t count = 0;
188   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
189   {
190     ++count;
191     AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
192
193     cont->AddElement(track);
194   }
195   cont->SetTitle(Form("N=%d", count));
196   cont->MakeTracks();
197
198   gEve->Redraw3D();
199
200   return cont;
201 }
202
203 TEveElementList* esd_tracks_MI()
204 {
205   AliESDEvent* esd = AliEveEventManager::AssertESD();
206
207   TEveElementList* cont = new TEveElementList("ESD Tracks MI");
208   gEve->AddElement(cont);
209
210   Int_t count = 0;
211   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
212   {
213     ++count;
214     AliESDtrack* at = esd->GetTrack(n);
215     TEveLine* l = new TEveLine; 
216     l->SetLineColor(5);
217     at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
218     l->SetName(Form("ESDTrackMI %d", at->GetID()));
219     l->SetElementTitle(esd_track_title(at));
220     l->SetSourceObject(at);
221
222     cont->AddElement(l);
223   }
224   cont->SetTitle(Form("N=%d", count));
225
226   gEve->Redraw3D();
227
228   return cont;
229 }
230
231
232 //==============================================================================
233 // esd_tracks_from_array()
234 //==============================================================================
235
236 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
237 {
238   // Retrieves AliESDTrack's from collection.
239   // See example usage with AliAnalysisTrackCuts in the next function.
240
241   if (esd == 0) esd = AliEveEventManager::AssertESD();
242
243   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
244   cont->SetMainColor(6);
245
246   esd_track_propagator_setup(cont->GetPropagator(),
247                              0.1*esd->GetMagneticField(), 520);
248
249   gEve->AddElement(cont);
250
251   Int_t    count = 0;
252   TIter    next(col);
253   TObject *obj;
254   while ((obj = next()) != 0)
255   {
256     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
257     {
258       Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
259               obj->GetName(), obj->GetTitle());
260       continue;
261     }
262
263     ++count;
264     AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
265
266     AliEveTrack* track = esd_make_track(at, cont);
267     cont->AddElement(track);
268   }
269   cont->SetTitle(Form("N=%d", count));
270   cont->MakeTracks();
271
272   gEve->Redraw3D();
273
274   return cont;
275 }
276
277 void esd_tracks_alianalcuts_demo()
278 {
279   AliESDEvent* esd = AliEveEventManager::AssertESD();
280
281   AliESDtrackCuts atc;
282   atc.SetPtRange(0.1, 5);
283   atc.SetRapRange(-1, 1);
284
285   esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
286 }
287
288
289 //==============================================================================
290 // esd_tracks_by_category
291 //==============================================================================
292
293 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
294 {
295   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
296   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
297
298   Float_t b[2];
299   Float_t bRes[2];
300   Float_t bCov[3];
301   esdTrack->GetImpactParameters(b,bCov);
302   if (bCov[0] <= 0 || bCov[2] <= 0)
303   {
304     printf("Estimated b resolution lower or equal zero!\n");
305     bCov[0] = bCov[2] = 0;
306   }
307   bRes[0] = TMath::Sqrt(bCov[0]);
308   bRes[1] = TMath::Sqrt(bCov[2]);
309
310   // -----------------------------------
311   // How to get to a n-sigma cut?
312   //
313   // The accumulated statistics from 0 to d is
314   //
315   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
316   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
317   //
318   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
319   // Can this be expressed in a different way?
320
321   if (bRes[0] == 0 || bRes[1] == 0)
322     return -1;
323
324   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
325
326   // stupid rounding problem screws up everything:
327   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
328   if (TMath::Exp(-d * d / 2) < 1e-10)
329     return 1000;
330
331   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
332   return d;
333 }
334
335 TEveElementList* esd_tracks_by_category()
336 {
337   // Import ESD tracks, separate them into several containers
338   // according to primary-vertex cut and ITS&TPC refit status.
339
340   AliESDEvent* esd = AliEveEventManager::AssertESD();
341
342   TEveElementList* cont = new TEveElementList("ESD Tracks by category");
343   gEve->AddElement(cont);
344
345   const Int_t   nCont = 7;
346   const Float_t maxR  = 520;
347   const Float_t magF  = 0.1*esd->GetMagneticField();
348
349   TEveTrackList *tl[nCont];
350   Int_t          tc[nCont];
351   Int_t          count = 0;
352
353   tl[0] = new TEveTrackList("Sigma < 3");
354   tc[0] = 0;
355   esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
356   tl[0]->SetMainColor(3);
357   cont->AddElement(tl[0]);
358
359   tl[1] = new TEveTrackList("3 < Sigma < 5");
360   tc[1] = 0;
361   esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
362   tl[1]->SetMainColor(7);
363   cont->AddElement(tl[1]);
364
365   tl[2] = new TEveTrackList("5 < Sigma");
366   tc[2] = 0;
367   esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
368   tl[2]->SetMainColor(46);
369   cont->AddElement(tl[2]);
370
371   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
372   tc[3] = 0;
373   esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
374   tl[3]->SetMainColor(41);
375   cont->AddElement(tl[3]);
376
377   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
378   tc[4] = 0;
379   esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
380   tl[4]->SetMainColor(48);
381   cont->AddElement(tl[4]);
382
383   tl[5] = new TEveTrackList("no TPC refit");
384   tc[5] = 0;
385   esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
386   tl[5]->SetMainColor(kRed);
387   cont->AddElement(tl[5]);
388
389   tl[6] = new TEveTrackList("ITS stand-alone");
390   tc[6] = 0;
391   esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
392   tl[6]->SetMainColor(kMagenta - 9);
393   cont->AddElement(tl[6]);
394
395   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
396   {
397     AliESDtrack* at = esd->GetTrack(n);
398
399     Float_t s  = get_sigma_to_vertex(at);
400     Int_t   ti;
401     if      (s <  3) ti = 0;
402     else if (s <= 5) ti = 1;
403     else             ti = 2;
404
405     if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
406     {
407       ti = 6;
408     }
409     else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
410     {
411       ti = 5;
412     }
413     else if ( ! at->IsOn(AliESDtrack::kITSrefit))
414     {
415       ti = (ti == 2) ? 4 : 3;
416     }
417
418     TEveTrackList* tlist = tl[ti];
419     ++tc[ti];
420     ++count;
421
422     AliEveTrack* track = esd_make_track(at, tlist);
423     track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
424     tlist->AddElement(track);
425   }
426
427   for (Int_t ti = 0; ti < nCont; ++ti)
428   {
429     TEveTrackList* tlist = tl[ti];
430     tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
431     tlist->SetTitle(Form("N tracks=%d", tc[ti]));
432
433     tlist->MakeTracks();
434
435     Bool_t good_cont = ti <= 1;
436     if (AliEveTrackCounter::IsActive())
437     {
438       AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
439     }
440     else
441     {
442       if ( ! good_cont)
443         tlist->SetLineStyle(6);
444     }
445   }
446   cont->SetTitle(Form("N all tracks = %d", count));
447   // ??? The following does not always work:
448   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
449
450   gEve->Redraw3D();
451
452   return cont;
453 }
454
455
456 //==============================================================================
457 // esd_tracks_by_anal_cuts
458 //==============================================================================
459
460 AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
461
462 TEveElementList* esd_tracks_by_anal_cuts()
463 {
464   AliESDEvent* esd = AliEveEventManager::AssertESD();
465
466   if (g_esd_tracks_anal_cuts == 0)
467   {
468     gSystem->Load("libPWG0base");
469     gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/CreateStandardCuts.C");
470     AliPWG0Helper::AnalysisMode mode = AliPWG0Helper::kTPCITS;
471     if (TMath::Abs(esd->GetMagneticField()) > 0.01)
472       mode |= AliPWG0Helper::kFieldOn;
473     g_esd_tracks_anal_cuts = CreateTrackCuts(mode, kFALSE);
474   }
475
476   TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
477   gEve->AddElement(cont);
478
479   const Int_t   nCont = 2;
480   const Float_t maxR  = 520;
481   const Float_t magF  = 0.1*esd->GetMagneticField();
482
483   TEveTrackList *tl[nCont];
484   Int_t          tc[nCont];
485   Int_t          count = 0;
486
487   tl[0] = new TEveTrackList("Passed");
488   tc[0] = 0;
489   esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
490   tl[0]->SetMainColor(3);
491   cont->AddElement(tl[0]);
492
493   tl[1] = new TEveTrackList("Rejected");
494   tc[1] = 0;
495   esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
496   tl[1]->SetMainColor(kRed);
497   cont->AddElement(tl[1]);
498
499   for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
500   {
501     AliESDtrack* at = esd->GetTrack(n);
502
503     Float_t s  = get_sigma_to_vertex(at);
504     Int_t   ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
505
506     TEveTrackList* tlist = tl[ti];
507     ++tc[ti];
508     ++count;
509
510     AliEveTrack* track = esd_make_track(at, tlist);
511     track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
512     tlist->AddElement(track);
513   }
514
515   for (Int_t ti = 0; ti < nCont; ++ti)
516   {
517     TEveTrackList* tlist = tl[ti];
518     tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
519     tlist->SetTitle(Form("N tracks=%d", tc[ti]));
520
521     tlist->MakeTracks();
522
523     Bool_t good_cont = ti < 1;
524     if (AliEveTrackCounter::IsActive())
525     {
526       AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
527     }
528     else
529     {
530       if ( ! good_cont)
531         tlist->SetLineStyle(6);
532     }
533   }
534   cont->SetTitle(Form("N all tracks = %d", count));
535   // ??? The following does not always work:
536   cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
537
538   gEve->Redraw3D();
539
540   return cont;
541 }