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__)
14 #include "TGListTree.h"
15 #include "TEveVSDStructs.h"
16 #include "TEveManager.h"
17 #include "TEveTrackPropagator.h"
19 #include "AliESDEvent.h"
20 #include "AliESDtrackCuts.h"
21 #include "AliESDtrack.h"
22 #include "AliExternalTrackParam.h"
24 #include "EVE/EveBase/AliEveTrack.h"
25 #include "EVE/EveBase/AliEveMagField.h"
26 #include "EVE/EveBase/AliEveEventManager.h"
31 // Use inner-tpc track params when its refit failed.
32 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
34 // Use magnetic-field as retrieved from GRP.
35 Bool_t g_esd_tracks_true_field = kTRUE;
37 // Use Runge-Kutta track stepper.
38 Bool_t g_esd_tracks_rk_stepper = kFALSE;
41 //==============================================================================
43 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
44 Float_t magF, Float_t maxR)
46 if (g_esd_tracks_true_field)
48 trkProp->SetMagFieldObj(new AliEveMagField);
52 trkProp->SetMagField(magF);
54 if (g_esd_tracks_rk_stepper)
56 trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
58 trkProp->SetMaxR(maxR);
61 //==============================================================================
63 TString esd_track_title(AliESDtrack* t)
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));
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);
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,
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));
103 if (t->IsOn(AliESDtrack::kESDpid))
107 s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
113 //==============================================================================
115 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
117 // Add additional track parameters as a path-mark to track.
122 Double_t pbuf[3], vbuf[3];
126 TEvePathMark pm(TEvePathMark::kReference);
129 track->AddPathMark(pm);
132 //==============================================================================
134 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
136 // Make a standard track representation and put it into given container.
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.
142 const AliExternalTrackParam* tp = at;
144 Bool_t innerTaken = kFALSE;
145 if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
147 tp = at->GetInnerParam();
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);
157 // Add inner/outer track parameters as path-marks.
158 if (at->IsOn(AliESDtrack::kTPCrefit))
162 esd_track_add_param(track, at->GetInnerParam());
164 esd_track_add_param(track, at->GetOuterParam());
171 //==============================================================================
173 //==============================================================================
175 TEveTrackList* esd_tracks()
177 AliESDEvent* esd = AliEveEventManager::AssertESD();
179 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
180 cont->SetMainColor(6);
182 esd_track_propagator_setup(cont->GetPropagator(),
183 0.1*esd->GetMagneticField(), 520);
185 gEve->AddElement(cont);
188 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
191 AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
193 cont->AddElement(track);
195 cont->SetTitle(Form("N=%d", count));
203 TEveElementList* esd_tracks_MI()
205 AliESDEvent* esd = AliEveEventManager::AssertESD();
207 TEveElementList* cont = new TEveElementList("ESD Tracks MI");
208 gEve->AddElement(cont);
211 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
214 AliESDtrack* at = esd->GetTrack(n);
215 TEveLine* l = new TEveLine;
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);
224 cont->SetTitle(Form("N=%d", count));
232 //==============================================================================
233 // esd_tracks_from_array()
234 //==============================================================================
236 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
238 // Retrieves AliESDTrack's from collection.
239 // See example usage with AliAnalysisTrackCuts in the next function.
241 if (esd == 0) esd = AliEveEventManager::AssertESD();
243 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
244 cont->SetMainColor(6);
246 esd_track_propagator_setup(cont->GetPropagator(),
247 0.1*esd->GetMagneticField(), 520);
249 gEve->AddElement(cont);
254 while ((obj = next()) != 0)
256 if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
258 Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
259 obj->GetName(), obj->GetTitle());
264 AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
266 AliEveTrack* track = esd_make_track(at, cont);
267 cont->AddElement(track);
269 cont->SetTitle(Form("N=%d", count));
277 void esd_tracks_alianalcuts_demo()
279 AliESDEvent* esd = AliEveEventManager::AssertESD();
282 atc.SetPtRange(0.1, 5);
283 atc.SetRapRange(-1, 1);
285 esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
289 //==============================================================================
290 // esd_tracks_by_category
291 //==============================================================================
293 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
295 // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
296 // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
301 esdTrack->GetImpactParameters(b,bCov);
302 if (bCov[0] <= 0 || bCov[2] <= 0)
304 printf("Estimated b resolution lower or equal zero!\n");
305 bCov[0] = bCov[2] = 0;
307 bRes[0] = TMath::Sqrt(bCov[0]);
308 bRes[1] = TMath::Sqrt(bCov[2]);
310 // -----------------------------------
311 // How to get to a n-sigma cut?
313 // The accumulated statistics from 0 to d is
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)
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?
321 if (bRes[0] == 0 || bRes[1] == 0)
324 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
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)
331 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
335 TEveElementList* esd_tracks_by_category()
337 // Import ESD tracks, separate them into several containers
338 // according to primary-vertex cut and ITS&TPC refit status.
340 AliESDEvent* esd = AliEveEventManager::AssertESD();
342 TEveElementList* cont = new TEveElementList("ESD Tracks by category");
343 gEve->AddElement(cont);
345 const Int_t nCont = 7;
346 const Float_t maxR = 520;
347 const Float_t magF = 0.1*esd->GetMagneticField();
349 TEveTrackList *tl[nCont];
353 tl[0] = new TEveTrackList("Sigma < 3");
355 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
356 tl[0]->SetMainColor(3);
357 cont->AddElement(tl[0]);
359 tl[1] = new TEveTrackList("3 < Sigma < 5");
361 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
362 tl[1]->SetMainColor(7);
363 cont->AddElement(tl[1]);
365 tl[2] = new TEveTrackList("5 < Sigma");
367 esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
368 tl[2]->SetMainColor(46);
369 cont->AddElement(tl[2]);
371 tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
373 esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
374 tl[3]->SetMainColor(41);
375 cont->AddElement(tl[3]);
377 tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
379 esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
380 tl[4]->SetMainColor(48);
381 cont->AddElement(tl[4]);
383 tl[5] = new TEveTrackList("no TPC refit");
385 esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
386 tl[5]->SetMainColor(kRed);
387 cont->AddElement(tl[5]);
389 tl[6] = new TEveTrackList("ITS stand-alone");
391 esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
392 tl[6]->SetMainColor(kMagenta - 9);
393 cont->AddElement(tl[6]);
395 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
397 AliESDtrack* at = esd->GetTrack(n);
399 Float_t s = get_sigma_to_vertex(at);
402 else if (s <= 5) ti = 1;
405 if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
409 else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
413 else if ( ! at->IsOn(AliESDtrack::kITSrefit))
415 ti = (ti == 2) ? 4 : 3;
418 TEveTrackList* tlist = tl[ti];
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);
427 for (Int_t ti = 0; ti < nCont; ++ti)
429 TEveTrackList* tlist = tl[ti];
430 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
431 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
435 Bool_t good_cont = ti <= 1;
436 if (AliEveTrackCounter::IsActive())
438 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
443 tlist->SetLineStyle(6);
446 cont->SetTitle(Form("N all tracks = %d", count));
447 // ??? The following does not always work:
448 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
456 //==============================================================================
457 // esd_tracks_by_anal_cuts
458 //==============================================================================
460 AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
462 TEveElementList* esd_tracks_by_anal_cuts()
464 AliESDEvent* esd = AliEveEventManager::AssertESD();
466 if (g_esd_tracks_anal_cuts == 0)
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);
476 TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
477 gEve->AddElement(cont);
479 const Int_t nCont = 2;
480 const Float_t maxR = 520;
481 const Float_t magF = 0.1*esd->GetMagneticField();
483 TEveTrackList *tl[nCont];
487 tl[0] = new TEveTrackList("Passed");
489 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
490 tl[0]->SetMainColor(3);
491 cont->AddElement(tl[0]);
493 tl[1] = new TEveTrackList("Rejected");
495 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
496 tl[1]->SetMainColor(kRed);
497 cont->AddElement(tl[1]);
499 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
501 AliESDtrack* at = esd->GetTrack(n);
503 Float_t s = get_sigma_to_vertex(at);
504 Int_t ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
506 TEveTrackList* tlist = tl[ti];
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);
515 for (Int_t ti = 0; ti < nCont; ++ti)
517 TEveTrackList* tlist = tl[ti];
518 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
519 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
523 Bool_t good_cont = ti < 1;
524 if (AliEveTrackCounter::IsActive())
526 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
531 tlist->SetLineStyle(6);
534 cont->SetTitle(Form("N all tracks = %d", count));
535 // ??? The following does not always work:
536 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);