2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
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 **************************************************************************/
10 #if !defined(__CINT__) || defined(__MAKECINT__)
16 #include <TGListTree.h>
18 #include <TEveVSDStructs.h>
19 #include <TEveManager.h>
20 #include <TEveTrackPropagator.h>
22 #include <EveBase/AliEveTrack.h>
23 #include <EveBase/AliEveTrackCounter.h>
24 #include <EveBase/AliEveMagField.h>
25 #include <EveBase/AliEveEventManager.h>
27 #include <AliESDEvent.h>
28 #include <AliESDtrackCuts.h>
29 #include <AliESDtrack.h>
30 #include <AliExternalTrackParam.h>
32 #include <AliPWG0Helper.h>
37 // Use inner-tpc track params when its refit failed.
38 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
40 // Use magnetic-field as retrieved from GRP.
41 Bool_t g_esd_tracks_true_field = kTRUE;
43 // Use Runge-Kutta track stepper.
44 Bool_t g_esd_tracks_rk_stepper = kFALSE;
47 //==============================================================================
49 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
50 Float_t magF, Float_t maxR)
52 if (g_esd_tracks_true_field)
54 trkProp->SetMagFieldObj(new AliEveMagField);
58 trkProp->SetMagField(magF);
60 if (g_esd_tracks_rk_stepper)
62 trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
64 trkProp->SetMaxR(maxR);
67 //==============================================================================
69 TString esd_track_title(AliESDtrack* t)
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));
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);
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,
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));
109 if (t->IsOn(AliESDtrack::kESDpid))
113 s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
119 //==============================================================================
121 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
123 // Add additional track parameters as a path-mark to track.
128 Double_t pbuf[3], vbuf[3];
132 TEvePathMark pm(TEvePathMark::kReference);
135 track->AddPathMark(pm);
138 //==============================================================================
140 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
142 // Make a standard track representation and put it into given container.
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.
148 const AliExternalTrackParam* tp = at;
150 Bool_t innerTaken = kFALSE;
151 if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
153 tp = at->GetInnerParam();
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);
163 // Add inner/outer track parameters as path-marks.
164 if (at->IsOn(AliESDtrack::kTPCrefit))
168 esd_track_add_param(track, at->GetInnerParam());
170 esd_track_add_param(track, at->GetOuterParam());
177 //==============================================================================
179 //==============================================================================
181 TEveTrackList* esd_tracks()
183 AliESDEvent* esd = AliEveEventManager::AssertESD();
185 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
186 cont->SetMainColor(6);
188 esd_track_propagator_setup(cont->GetPropagator(),
189 0.1*esd->GetMagneticField(), 520);
191 gEve->AddElement(cont);
194 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
197 AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
199 cont->AddElement(track);
201 cont->SetTitle(Form("N=%d", count));
209 TEveTrackList* esd_tracks_MI()
211 AliESDEvent* esd = AliEveEventManager::AssertESD();
213 TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
214 cont->SetLineColor(5);
215 gEve->AddElement(cont);
218 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
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);
228 at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
230 l->SetLockPoints(kTRUE);
233 cont->SetTitle(Form("N=%d", count));
241 //==============================================================================
242 // esd_tracks_from_array()
243 //==============================================================================
245 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
247 // Retrieves AliESDTrack's from collection.
248 // See example usage with AliAnalysisTrackCuts in the next function.
250 if (esd == 0) esd = AliEveEventManager::AssertESD();
252 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
253 cont->SetMainColor(6);
255 esd_track_propagator_setup(cont->GetPropagator(),
256 0.1*esd->GetMagneticField(), 520);
258 gEve->AddElement(cont);
263 while ((obj = next()) != 0)
265 if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
267 Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
268 obj->GetName(), obj->GetTitle());
273 AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
275 AliEveTrack* track = esd_make_track(at, cont);
276 cont->AddElement(track);
278 cont->SetTitle(Form("N=%d", count));
286 void esd_tracks_alianalcuts_demo()
288 AliESDEvent* esd = AliEveEventManager::AssertESD();
291 atc.SetPtRange(0.1, 5);
292 atc.SetRapRange(-1, 1);
294 esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
298 //==============================================================================
299 // esd_tracks_by_category
300 //==============================================================================
302 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
304 // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
305 // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
310 esdTrack->GetImpactParameters(b,bCov);
311 if (bCov[0] <= 0 || bCov[2] <= 0)
313 printf("Estimated b resolution lower or equal zero!\n");
314 bCov[0] = bCov[2] = 0;
316 bRes[0] = TMath::Sqrt(bCov[0]);
317 bRes[1] = TMath::Sqrt(bCov[2]);
319 // -----------------------------------
320 // How to get to a n-sigma cut?
322 // The accumulated statistics from 0 to d is
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)
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?
330 if (bRes[0] == 0 || bRes[1] == 0)
333 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
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)
340 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
344 TEveElementList* esd_tracks_by_category()
346 // Import ESD tracks, separate them into several containers
347 // according to primary-vertex cut and ITS&TPC refit status.
349 AliESDEvent* esd = AliEveEventManager::AssertESD();
351 TEveElementList* cont = new TEveElementList("ESD Tracks by category");
352 gEve->AddElement(cont);
354 const Int_t nCont = 7;
355 const Float_t maxR = 520;
356 const Float_t magF = 0.1*esd->GetMagneticField();
358 TEveTrackList *tl[nCont];
362 tl[0] = new TEveTrackList("Sigma < 3");
364 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
365 tl[0]->SetMainColor(3);
366 cont->AddElement(tl[0]);
368 tl[1] = new TEveTrackList("3 < Sigma < 5");
370 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
371 tl[1]->SetMainColor(7);
372 cont->AddElement(tl[1]);
374 tl[2] = new TEveTrackList("5 < Sigma");
376 esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
377 tl[2]->SetMainColor(46);
378 cont->AddElement(tl[2]);
380 tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
382 esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
383 tl[3]->SetMainColor(41);
384 cont->AddElement(tl[3]);
386 tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
388 esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
389 tl[4]->SetMainColor(48);
390 cont->AddElement(tl[4]);
392 tl[5] = new TEveTrackList("no TPC refit");
394 esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
395 tl[5]->SetMainColor(kRed);
396 cont->AddElement(tl[5]);
398 tl[6] = new TEveTrackList("ITS stand-alone");
400 esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
401 tl[6]->SetMainColor(kMagenta - 9);
402 cont->AddElement(tl[6]);
404 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
406 AliESDtrack* at = esd->GetTrack(n);
408 Float_t s = get_sigma_to_vertex(at);
411 else if (s <= 5) ti = 1;
414 if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
418 else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
422 else if ( ! at->IsOn(AliESDtrack::kITSrefit))
424 ti = (ti == 2) ? 4 : 3;
427 TEveTrackList* tlist = tl[ti];
431 AliEveTrack* track = esd_make_track(at, tlist);
432 track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
433 tlist->AddElement(track);
436 for (Int_t ti = 0; ti < nCont; ++ti)
438 TEveTrackList* tlist = tl[ti];
439 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
440 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
444 Bool_t good_cont = ti <= 1;
445 if (AliEveTrackCounter::IsActive())
447 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
452 tlist->SetLineStyle(6);
455 cont->SetTitle(Form("N all tracks = %d", count));
456 // ??? The following does not always work:
457 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
465 //==============================================================================
466 // esd_tracks_by_anal_cuts
467 //==============================================================================
469 AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
471 TEveElementList* esd_tracks_by_anal_cuts()
473 AliESDEvent* esd = AliEveEventManager::AssertESD();
475 if (g_esd_tracks_anal_cuts == 0)
477 gSystem->Load("libPWG0base");
478 gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/CreateStandardCuts.C");
479 Int_t mode = AliPWG0Helper::kTPC;
480 if (TMath::Abs(esd->GetMagneticField()) > 0.01)
481 mode |= AliPWG0Helper::kFieldOn;
482 gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
485 TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
486 gEve->AddElement(cont);
488 const Int_t nCont = 2;
489 const Float_t maxR = 520;
490 const Float_t magF = 0.1*esd->GetMagneticField();
492 TEveTrackList *tl[nCont];
496 tl[0] = new TEveTrackList("Passed");
498 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
499 tl[0]->SetMainColor(3);
500 cont->AddElement(tl[0]);
502 tl[1] = new TEveTrackList("Rejected");
504 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
505 tl[1]->SetMainColor(kRed);
506 cont->AddElement(tl[1]);
508 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
510 AliESDtrack* at = esd->GetTrack(n);
512 Float_t s = get_sigma_to_vertex(at);
513 Int_t ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
515 TEveTrackList* tlist = tl[ti];
519 AliEveTrack* track = esd_make_track(at, tlist);
520 track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
521 tlist->AddElement(track);
524 for (Int_t ti = 0; ti < nCont; ++ti)
526 TEveTrackList* tlist = tl[ti];
527 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
528 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
532 Bool_t good_cont = ti < 1;
533 if (AliEveTrackCounter::IsActive())
535 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
540 tlist->SetLineStyle(6);
543 cont->SetTitle(Form("N all tracks = %d", count));
544 // ??? The following does not always work:
545 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);