In muon-related macros in EVE:
[u/mrichter/AliRoot.git] / EVE / alice-macros / esd_tracks.C
CommitLineData
5a5a1232 1// $Id$
d810d0de 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 *
51346b82 7 * full copyright notice. *
d810d0de 8 **************************************************************************/
5a5a1232 9
e9ece5b4 10#if !defined(__CINT__) || defined(__MAKECINT__)
11
15c12442 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
b264afb8 32#include <AliPWG0Helper.h>
e9ece5b4 33
34#endif
35
36
29207369 37// Use inner-tpc track params when its refit failed.
38Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
39
6336e5d7 40// Use magnetic-field as retrieved from GRP.
41Bool_t g_esd_tracks_true_field = kTRUE;
42
43// Use Runge-Kutta track stepper.
44Bool_t g_esd_tracks_rk_stepper = kFALSE;
45
46
47//==============================================================================
48
49void 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
29207369 69TString esd_track_title(AliESDtrack* t)
ad5abc55 70{
71 TString s;
29207369 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);
d6962cc4 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);
29207369 85
86 s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
d6962cc4 87 "pT = %.3f + %.3f - %.3f [%.3f]\n"
88 "P = (%.3f, %.3f, %.3f)\n"
89 "V = (%.3f, %.3f, %.3f)\n",
29207369 90 idx.Data(), lbl.Data(), t->Charge(), 0,
d6962cc4 91 pt, ptM - pt, pt - ptm, ptsig*ptsq,
92 p[0], p[1], p[2],
93 v[0], v[1], v[2]);
29207369 94
ad5abc55 95 Int_t o;
d6962cc4 96 s += "Det (in,out,refit,pid):\n";
ad5abc55 97 o = AliESDtrack::kITSin;
6b7fd3a4 98 s += Form("ITS (%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55 99 o = AliESDtrack::kTPCin;
6b7fd3a4 100 s += Form("TPC(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55 101 o = AliESDtrack::kTRDin;
6b7fd3a4 102 s += Form("TRD(%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55 103 o = AliESDtrack::kTOFin;
6b7fd3a4 104 s += Form("TOF(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
ad5abc55 105 o = AliESDtrack::kHMPIDout;
6b7fd3a4 106 s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
ad5abc55 107 s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
86f12f78 108
29207369 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 }
51346b82 115
29207369 116 return s;
86f12f78 117}
118
6336e5d7 119//==============================================================================
120
e9ece5b4 121void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
c222eb1d 122{
29207369 123 // Add additional track parameters as a path-mark to track.
124
c222eb1d 125 if (tp == 0)
126 return;
127
29207369 128 Double_t pbuf[3], vbuf[3];
c222eb1d 129 tp->GetXYZ(vbuf);
130 tp->GetPxPyPz(pbuf);
131
a15e6d7d 132 TEvePathMark pm(TEvePathMark::kReference);
133 pm.fV.Set(vbuf);
134 pm.fP.Set(pbuf);
c222eb1d 135 track->AddPathMark(pm);
136}
137
6336e5d7 138//==============================================================================
139
e9ece5b4 140AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
29207369 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.
e9ece5b4 147
148 const AliExternalTrackParam* tp = at;
29207369 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
50ac85c6 157 AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
29207369 158 track->SetAttLineAttMarker(cont);
0e33c639 159 track->SetName(Form("AliEveTrack %d", at->GetID()));
29207369 160 track->SetElementTitle(esd_track_title(at));
9b6216c0 161 track->SetSourceObject(at);
29207369 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 }
a15e6d7d 172
29207369 173 return track;
174}
86f12f78 175
6336e5d7 176
177//==============================================================================
178// esd_tracks()
179//==============================================================================
180
29207369 181TEveTrackList* esd_tracks()
5a5a1232 182{
d810d0de 183 AliESDEvent* esd = AliEveEventManager::AssertESD();
5a5a1232 184
51346b82 185 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3 186 cont->SetMainColor(6);
6336e5d7 187
188 esd_track_propagator_setup(cont->GetPropagator(),
189 0.1*esd->GetMagneticField(), 520);
5a5a1232 190
84aff7a4 191 gEve->AddElement(cont);
5a5a1232 192
29207369 193 Int_t count = 0;
194 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
d31ac39d 195 {
29207369 196 ++count;
0e33c639 197 AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
5a5a1232 198
29207369 199 cont->AddElement(track);
200 }
201 cont->SetTitle(Form("N=%d", count));
202 cont->MakeTracks();
5a5a1232 203
29207369 204 gEve->Redraw3D();
5a5a1232 205
29207369 206 return cont;
207}
5a5a1232 208
06219f85 209TEveTrackList* esd_tracks_MI()
29207369 210{
211 AliESDEvent* esd = AliEveEventManager::AssertESD();
c222eb1d 212
06219f85 213 TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
214 cont->SetLineColor(5);
29207369 215 gEve->AddElement(cont);
c222eb1d 216
29207369 217 Int_t count = 0;
218 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
219 {
220 ++count;
221 AliESDtrack* at = esd->GetTrack(n);
06219f85 222 AliEveTrack* l = new AliEveTrack(at, cont->GetPropagator());
29207369 223 l->SetName(Form("ESDTrackMI %d", at->GetID()));
224 l->SetElementTitle(esd_track_title(at));
06219f85 225 l->SetAttLineAttMarker(cont);
9b6216c0 226 l->SetSourceObject(at);
c222eb1d 227
06219f85 228 at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
229
230 l->SetLockPoints(kTRUE);
29207369 231 cont->AddElement(l);
5a5a1232 232 }
29207369 233 cont->SetTitle(Form("N=%d", count));
32e219c2 234
84aff7a4 235 gEve->Redraw3D();
5a5a1232 236
237 return cont;
238}
239
6336e5d7 240
241//==============================================================================
a8600b56 242// esd_tracks_from_array()
6336e5d7 243//==============================================================================
a8600b56 244
84aff7a4 245TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
5a5a1232 246{
247 // Retrieves AliESDTrack's from collection.
248 // See example usage with AliAnalysisTrackCuts in the next function.
249
d810d0de 250 if (esd == 0) esd = AliEveEventManager::AssertESD();
86f12f78 251
51346b82 252 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3 253 cont->SetMainColor(6);
6336e5d7 254
255 esd_track_propagator_setup(cont->GetPropagator(),
256 0.1*esd->GetMagneticField(), 520);
5a5a1232 257
84aff7a4 258 gEve->AddElement(cont);
5a5a1232 259
5a5a1232 260 Int_t count = 0;
261 TIter next(col);
262 TObject *obj;
84aff7a4 263 while ((obj = next()) != 0)
d31ac39d 264 {
29207369 265 if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
266 {
267 Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
5a5a1232 268 obj->GetName(), obj->GetTitle());
269 continue;
270 }
271
272 ++count;
e9ece5b4 273 AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
5a5a1232 274
0e33c639 275 AliEveTrack* track = esd_make_track(at, cont);
29207369 276 cont->AddElement(track);
5a5a1232 277 }
6b7fd3a4 278 cont->SetTitle(Form("N=%d", count));
5a5a1232 279 cont->MakeTracks();
32e219c2 280
84aff7a4 281 gEve->Redraw3D();
5a5a1232 282
283 return cont;
284}
285
286void esd_tracks_alianalcuts_demo()
287{
d810d0de 288 AliESDEvent* esd = AliEveEventManager::AssertESD();
5a5a1232 289
fcfc8b01 290 AliESDtrackCuts atc;
5a5a1232 291 atc.SetPtRange(0.1, 5);
292 atc.SetRapRange(-1, 1);
293
fcfc8b01 294 esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
5a5a1232 295}
a8600b56 296
6336e5d7 297
298//==============================================================================
29207369 299// esd_tracks_by_category
6336e5d7 300//==============================================================================
a8600b56 301
302Float_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);
29207369 311 if (bCov[0] <= 0 || bCov[2] <= 0)
312 {
d31ac39d 313 printf("Estimated b resolution lower or equal zero!\n");
29207369 314 bCov[0] = bCov[2] = 0;
a8600b56 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
29207369 330 if (bRes[0] == 0 || bRes[1] == 0)
a8600b56 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
29207369 344TEveElementList* esd_tracks_by_category()
a8600b56 345{
440a8bc6 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;
a8600b56 497}
5320a363 498
499
500//==============================================================================
501// esd_tracks_by_anal_cuts
502//==============================================================================
503
504AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
505
506TEveElementList* 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");
5065a65b 514 Int_t mode = AliPWG0Helper::kTPC;
5320a363 515 if (TMath::Abs(esd->GetMagneticField()) > 0.01)
516 mode |= AliPWG0Helper::kFieldOn;
15c12442 517 gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
5320a363 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}