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