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__)
15 #include <TGListTree.h>
16 #include <TEveVSDStructs.h>
17 #include <TEveManager.h>
18 #include <TEveTrackPropagator.h>
20 #include <AliESDtrackCuts.h>
21 #include <AliPWG0Helper.h>
22 #include <AliESDEvent.h>
23 #include <AliESDfriend.h>
24 #include <AliESDtrack.h>
25 #include <AliESDfriendTrack.h>
26 #include <AliExternalTrackParam.h>
27 #include <AliEveTrack.h>
28 #include <AliEveTrackCounter.h>
29 #include <AliEveMagField.h>
30 #include <AliEveEventManager.h>
34 // Use inner-tpc track params when its refit failed.
35 Bool_t g_esd_tracks_use_ip_on_failed_its_refit = kFALSE;
37 // Use magnetic-field as retrieved from GRP.
38 Bool_t g_esd_tracks_true_field = kTRUE;
40 // Use Runge-Kutta track stepper.
41 Bool_t g_esd_tracks_rk_stepper = kFALSE;
44 //==============================================================================
46 void esd_track_propagator_setup(TEveTrackPropagator* trkProp,
47 Float_t magF, Float_t maxR)
49 if (g_esd_tracks_true_field)
51 trkProp->SetMagFieldObj(new AliEveMagField);
55 trkProp->SetMagField(magF);
57 if (g_esd_tracks_rk_stepper)
59 trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
61 trkProp->SetMaxR(maxR);
64 //==============================================================================
66 TString esd_track_title(AliESDtrack* t)
70 Int_t label = t->GetLabel(), index = t->GetID();
71 TString idx(index == kMinInt ? "<undef>" : Form("%d", index));
72 TString lbl(label == kMinInt ? "<undef>" : Form("%d", label));
77 Double_t pt = t->Pt();
78 Double_t ptsig = TMath::Sqrt(t->GetSigma1Pt2());
79 Double_t ptsq = pt*pt;
80 Double_t ptm = pt / (1.0 + pt*ptsig);
81 Double_t ptM = pt / (1.0 - pt*ptsig);
83 s = Form("Index=%s, Label=%s\nChg=%d, Pdg=%d\n"
84 "pT = %.3f + %.3f - %.3f [%.3f]\n"
85 "P = (%.3f, %.3f, %.3f)\n"
86 "V = (%.3f, %.3f, %.3f)\n",
87 idx.Data(), lbl.Data(), t->Charge(), 0,
88 pt, ptM - pt, pt - ptm, ptsig*ptsq,
93 s += "Det (in,out,refit,pid):\n";
94 o = AliESDtrack::kITSin;
95 s += Form("ITS (%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
96 o = AliESDtrack::kTPCin;
97 s += Form("TPC(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
98 o = AliESDtrack::kTRDin;
99 s += Form("TRD(%d,%d,%d,%d) ", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
100 o = AliESDtrack::kTOFin;
101 s += Form("TOF(%d,%d,%d,%d)\n", t->IsOn(o), t->IsOn(o<<1), t->IsOn(o<<2), t->IsOn(o<<3));
102 o = AliESDtrack::kHMPIDout;
103 s += Form("HMPID(out=%d,pid=%d)\n", t->IsOn(o), t->IsOn(o<<1));
104 s += Form("ESD pid=%d", t->IsOn(AliESDtrack::kESDpid));
106 if (t->IsOn(AliESDtrack::kESDpid))
110 s += Form("\n[%.2f %.2f %.2f %.2f %.2f]", pid[0], pid[1], pid[2], pid[3], pid[4]);
116 //==============================================================================
118 void esd_track_add_param(AliEveTrack* track, const AliExternalTrackParam* tp)
120 // Add additional track parameters as a path-mark to track.
125 Double_t pbuf[3], vbuf[3];
129 TEvePathMark pm(TEvePathMark::kReference);
132 track->AddPathMark(pm);
135 //==============================================================================
137 AliEveTrack* esd_make_track_TPC(AliESDtrack *at, AliESDfriendTrack* aft, TEveTrackList* cont)
139 // Make a TPC track representation and put it into given container.
141 AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
142 track->SetAttLineAttMarker(cont);
143 track->SetName(Form("AliEveTrack %d", at->GetID()));
144 track->SetElementTitle(esd_track_title(at));
145 track->SetSourceObject(at);
147 // Add inner/outer track parameters as start point and pathmark.
148 if (at->GetInnerParam()) track->SetStartParams(at->GetInnerParam());
150 if (aft->GetTPCOut()) esd_track_add_param(track, aft->GetTPCOut());
156 AliEveTrack* esd_make_track_ITS_standalone(AliESDtrack *at, AliESDfriendTrack* aft, TEveTrackList* cont)
158 // Make a ITS standalone track representation and put it into given container.
160 if ( !(!at->IsOn(AliESDtrack::kTPCin) &&
161 at->IsOn(AliESDtrack::kITSout)) ) return NULL; //only ITS standalone
162 AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
163 track->SetAttLineAttMarker(cont);
164 track->SetName(Form("AliEveTrack %d", at->GetID()));
165 track->SetElementTitle(esd_track_title(at));
166 track->SetSourceObject(at);
168 // Add inner/outer track parameters as path-marks.
169 if (aft->GetITSOut())
171 esd_track_add_param(track, aft->GetITSOut());
178 AliEveTrack* esd_make_track_ITS(AliESDtrack *at, AliESDfriendTrack* aft, TEveTrackList* cont)
180 // Make a ITS track representation and put it into given container.
182 if ( (!at->IsOn(AliESDtrack::kTPCin) &&
183 at->IsOn(AliESDtrack::kITSout)) ) return NULL; //ignore ITS standalone
184 AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
185 track->SetAttLineAttMarker(cont);
186 track->SetName(Form("AliEveTrack %d", at->GetID()));
187 track->SetElementTitle(esd_track_title(at));
188 track->SetSourceObject(at);
190 // Add inner/outer track parameters as path-marks.
191 if (aft->GetITSOut())
193 esd_track_add_param(track, aft->GetITSOut());
200 AliEveTrack* esd_make_track(AliESDtrack *at, TEveTrackList* cont)
202 // Make a standard track representation and put it into given container.
204 // Choose which parameters to use a track's starting point.
205 // If gkFixFailedITSExtr is TRUE (FALSE by default) and
206 // if ITS refit failed, take track parameters at inner TPC radius.
208 const AliExternalTrackParam* tp = at;
210 Bool_t innerTaken = kFALSE;
211 if ( ! at->IsOn(AliESDtrack::kITSrefit) && g_esd_tracks_use_ip_on_failed_its_refit)
213 tp = at->GetInnerParam();
217 AliEveTrack* track = new AliEveTrack(at, cont->GetPropagator());
218 track->SetAttLineAttMarker(cont);
219 track->SetName(Form("AliEveTrack %d", at->GetID()));
220 track->SetElementTitle(esd_track_title(at));
221 track->SetSourceObject(at);
223 // Add inner/outer track parameters as path-marks.
224 if (at->IsOn(AliESDtrack::kTPCrefit))
228 esd_track_add_param(track, at->GetInnerParam());
230 esd_track_add_param(track, at->GetOuterParam());
237 //==============================================================================
239 //==============================================================================
241 TEveTrackList* esd_tracks_TPC()
243 AliESDEvent* esd = AliEveEventManager::AssertESD();
244 AliESDfriend* esd_friend = AliEveEventManager::AssertESDfriend();
246 TEveTrackList* cont = new TEveTrackList("TPC Tracks");
247 cont->SetMainColor(kMagenta);
249 esd_track_propagator_setup(cont->GetPropagator(),
250 0.1*esd->GetMagneticField(), 520);
252 gEve->AddElement(cont);
255 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
258 if (!esd->GetTrack(n)) continue;
259 if (!esd_friend->GetTrack(n)) continue;
260 AliEveTrack* track = esd_make_track_TPC(esd->GetTrack(n), esd_friend->GetTrack(n), cont);
261 if (!track) continue;
263 cont->AddElement(track);
265 cont->SetTitle(Form("N=%d", count));
273 TEveTrackList* esd_tracks_ITS()
275 AliESDEvent* esd = AliEveEventManager::AssertESD();
276 AliESDfriend* esd_friend = AliEveEventManager::AssertESDfriend();
278 TEveTrackList* cont = new TEveTrackList("ITS Tracks");
279 cont->SetMainColor(kMagenta+3);
281 esd_track_propagator_setup(cont->GetPropagator(),
282 0.1*esd->GetMagneticField(), 520);
283 cont->GetPropagator()->SetMaxR(85.0);
284 cont->SetLineWidth(1);
286 gEve->AddElement(cont);
289 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
292 if (!esd->GetTrack(n)) continue;
293 if (!esd_friend->GetTrack(n)) continue;
294 AliEveTrack* track = esd_make_track_ITS(esd->GetTrack(n), esd_friend->GetTrack(n), cont);
295 if (!track) continue;
297 cont->AddElement(track);
299 cont->SetTitle(Form("N=%d", count));
307 TEveTrackList* esd_tracks_ITS_standalone()
309 AliESDEvent* esd = AliEveEventManager::AssertESD();
310 AliESDfriend* esd_friend = AliEveEventManager::AssertESDfriend();
312 TEveTrackList* cont = new TEveTrackList("ITS Standalone Tracks");
313 cont->SetMainColor(kBlue);
315 esd_track_propagator_setup(cont->GetPropagator(),
316 0.1*esd->GetMagneticField(), 520);
317 cont->GetPropagator()->SetMaxR(85.0);
318 cont->SetLineWidth(1);
320 gEve->AddElement(cont);
323 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
326 if (!esd->GetTrack(n)) continue;
327 if (!esd_friend->GetTrack(n)) continue;
328 AliEveTrack* track = esd_make_track_ITS_standalone(esd->GetTrack(n), esd_friend->GetTrack(n), cont);
329 if (!track) continue;
331 cont->AddElement(track);
333 cont->SetTitle(Form("N=%d", count));
341 TEveTrackList* esd_tracks()
343 AliESDEvent* esd = AliEveEventManager::AssertESD();
345 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
346 cont->SetMainColor(6);
348 esd_track_propagator_setup(cont->GetPropagator(),
349 0.1*esd->GetMagneticField(), 520);
351 gEve->AddElement(cont);
354 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
357 AliEveTrack* track = esd_make_track(esd->GetTrack(n), cont);
359 cont->AddElement(track);
361 cont->SetTitle(Form("N=%d", count));
369 TEveTrackList* esd_tracks_MI()
371 AliESDEvent* esd = AliEveEventManager::AssertESD();
373 TEveTrackList* cont = new TEveTrackList("ESD Tracks MI");
374 cont->SetLineColor(5);
375 gEve->AddElement(cont);
378 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
381 AliESDtrack* at = esd->GetTrack(n);
382 AliEveTrack* l = new AliEveTrack(at, cont->GetPropagator());
383 l->SetName(Form("ESDTrackMI %d", at->GetID()));
384 l->SetElementTitle(esd_track_title(at));
385 l->SetAttLineAttMarker(cont);
386 l->SetSourceObject(at);
388 at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
390 l->SetLockPoints(kTRUE);
393 cont->SetTitle(Form("N=%d", count));
401 //==============================================================================
402 // esd_tracks_from_array()
403 //==============================================================================
405 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
407 // Retrieves AliESDTrack's from collection.
408 // See example usage with AliAnalysisTrackCuts in the next function.
410 if (esd == 0) esd = AliEveEventManager::AssertESD();
412 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
413 cont->SetMainColor(6);
415 esd_track_propagator_setup(cont->GetPropagator(),
416 0.1*esd->GetMagneticField(), 520);
418 gEve->AddElement(cont);
423 while ((obj = next()) != 0)
425 if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE)
427 Warning("esd_tracks_from_array", "Object '%s', '%s' is not an AliESDtrack.",
428 obj->GetName(), obj->GetTitle());
433 AliESDtrack* at = reinterpret_cast<AliESDtrack*>(obj);
435 AliEveTrack* track = esd_make_track(at, cont);
436 cont->AddElement(track);
438 cont->SetTitle(Form("N=%d", count));
446 void esd_tracks_alianalcuts_demo()
448 AliESDEvent* esd = AliEveEventManager::AssertESD();
451 atc.SetPtRange(0.1, 5);
452 atc.SetRapRange(-1, 1);
454 esd_tracks_from_array(atc.GetAcceptedTracks(esd), esd);
458 //==============================================================================
459 // esd_tracks_by_category
460 //==============================================================================
462 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
464 // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
465 // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
470 esdTrack->GetImpactParameters(b,bCov);
471 if (bCov[0] <= 0 || bCov[2] <= 0)
473 printf("Estimated b resolution lower or equal zero!\n");
474 bCov[0] = bCov[2] = 0;
476 bRes[0] = TMath::Sqrt(bCov[0]);
477 bRes[1] = TMath::Sqrt(bCov[2]);
479 // -----------------------------------
480 // How to get to a n-sigma cut?
482 // The accumulated statistics from 0 to d is
484 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
485 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
487 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
488 // Can this be expressed in a different way?
490 if (bRes[0] == 0 || bRes[1] == 0)
493 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
495 // stupid rounding problem screws up everything:
496 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
497 if (TMath::Exp(-d * d / 2) < 1e-10)
500 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
504 TEveElementList* esd_tracks_by_category()
506 // Import ESD tracks, separate them into several containers
507 // according to primary-vertex cut and ITS&TPC refit status.
509 AliESDEvent* esd = AliEveEventManager::AssertESD();
511 TEveElementList* cont = new TEveElementList("ESD Tracks by category");
512 gEve->AddElement(cont);
514 const Int_t nCont = 9;
515 const Float_t maxR = 520;
516 const Float_t magF = 0.1*esd->GetMagneticField();
518 TEveTrackList *tl[nCont];
522 tl[0] = new TEveTrackList("Sigma < 3");
524 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
525 tl[0]->SetMainColor(3);
526 cont->AddElement(tl[0]);
528 tl[1] = new TEveTrackList("3 < Sigma < 5");
530 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
531 tl[1]->SetMainColor(7);
532 cont->AddElement(tl[1]);
534 tl[2] = new TEveTrackList("5 < Sigma");
536 esd_track_propagator_setup(tl[2]->GetPropagator(), magF, maxR);
537 tl[2]->SetMainColor(46);
538 cont->AddElement(tl[2]);
540 tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
542 esd_track_propagator_setup(tl[3]->GetPropagator(), magF, maxR);
543 tl[3]->SetMainColor(41);
544 cont->AddElement(tl[3]);
546 tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
548 esd_track_propagator_setup(tl[4]->GetPropagator(), magF, maxR);
549 tl[4]->SetMainColor(48);
550 cont->AddElement(tl[4]);
552 tl[5] = new TEveTrackList("no TPC refit");
554 esd_track_propagator_setup(tl[5]->GetPropagator(), magF, maxR);
555 tl[5]->SetMainColor(kRed);
556 cont->AddElement(tl[5]);
558 tl[6] = new TEveTrackList("ITS ncl>=3 & SPD Inner");
560 esd_track_propagator_setup(tl[6]->GetPropagator(), magF, maxR);
561 tl[6]->SetMainColor(kGreen);
562 cont->AddElement(tl[6]);
564 tl[7] = new TEveTrackList("ITS ncl>=3 & b<3 cm");
566 esd_track_propagator_setup(tl[7]->GetPropagator(), magF, maxR);
567 tl[7]->SetMainColor(kMagenta - 9);
568 cont->AddElement(tl[7]);
570 tl[8] = new TEveTrackList("ITS others");
572 esd_track_propagator_setup(tl[8]->GetPropagator(), magF, maxR);
573 tl[8]->SetMainColor(kRed + 2);
574 cont->AddElement(tl[8]);
576 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
578 AliESDtrack* at = esd->GetTrack(n);
580 Float_t s = get_sigma_to_vertex(at);
583 else if (s <= 5) ti = 1;
589 if (at->IsOn(AliESDtrack::kITSin) && ! at->IsOn(AliESDtrack::kTPCin))
591 UChar_t itsclmap = at->GetITSClusterMap();
592 Bool_t spdinner = (itsclmap & 3) != 0;
595 for (Int_t iter = 0; iter < 6; ++iter)
596 if (itsclmap & (1 << iter)) nclits ++;
600 dtobeam = TMath::Hypot(xyz[0], xyz[1]);
602 if ((nclits >= 3) && (spdinner))
604 else if ((nclits >= 3) && (dtobeam < 3.0))
609 else if (at->IsOn(AliESDtrack::kTPCin) && ! at->IsOn(AliESDtrack::kTPCrefit))
613 else if ( ! at->IsOn(AliESDtrack::kITSrefit))
615 ti = (ti == 2) ? 4 : 3;
618 TEveTrackList* tlist = tl[ti];
622 AliEveTrack* track = esd_make_track(at, tlist);
624 track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
626 track->SetName(Form("ESD Track idx=%d, dxy=%5.3f cl=%i", at->GetID(), dtobeam, nclits));
627 tlist->AddElement(track);
630 for (Int_t ti = 0; ti < nCont; ++ti)
632 TEveTrackList* tlist = tl[ti];
633 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
634 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
638 // Bool_t good_cont = ti <= 1;
639 Bool_t good_cont = ((ti == 6) || (ti == 7));
640 if (AliEveTrackCounter::IsActive())
642 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
647 tlist->SetLineStyle(6);
650 cont->SetTitle(Form("N all tracks = %d", count));
651 // ??? The following does not always work:
652 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);
660 //==============================================================================
661 // esd_tracks_by_anal_cuts
662 //==============================================================================
664 AliESDtrackCuts* g_esd_tracks_anal_cuts = 0;
666 TEveElementList* esd_tracks_by_anal_cuts()
668 AliESDEvent* esd = AliEveEventManager::AssertESD();
670 if (g_esd_tracks_anal_cuts == 0)
672 gSystem->Load("libPWGUDbase");
673 gROOT->ProcessLine(".L $ALICE_ROOT/PWGUD/CreateStandardCuts.C");
674 Int_t mode = AliPWG0Helper::kTPC;
675 if (TMath::Abs(esd->GetMagneticField()) > 0.01)
676 mode |= AliPWG0Helper::kFieldOn;
677 gROOT->ProcessLine(Form("g_esd_tracks_anal_cuts = CreateTrackCuts(%d, kFALSE)", mode));
680 TEveElementList* cont = new TEveElementList("ESD Tracks by Analysis Cuts");
681 gEve->AddElement(cont);
683 const Int_t nCont = 2;
684 const Float_t maxR = 520;
685 const Float_t magF = 0.1*esd->GetMagneticField();
687 TEveTrackList *tl[nCont];
691 tl[0] = new TEveTrackList("Passed");
693 esd_track_propagator_setup(tl[0]->GetPropagator(), magF, maxR);
694 tl[0]->SetMainColor(3);
695 cont->AddElement(tl[0]);
697 tl[1] = new TEveTrackList("Rejected");
699 esd_track_propagator_setup(tl[1]->GetPropagator(), magF, maxR);
700 tl[1]->SetMainColor(kRed);
701 cont->AddElement(tl[1]);
703 for (Int_t n = 0; n < esd->GetNumberOfTracks(); ++n)
705 AliESDtrack* at = esd->GetTrack(n);
707 Float_t s = get_sigma_to_vertex(at);
708 Int_t ti = (g_esd_tracks_anal_cuts->AcceptTrack(at)) ? 0 : 1;
710 TEveTrackList* tlist = tl[ti];
714 AliEveTrack* track = esd_make_track(at, tlist);
715 track->SetName(Form("ESD Track idx=%d, sigma=%5.3f", at->GetID(), s));
716 tlist->AddElement(track);
719 for (Int_t ti = 0; ti < nCont; ++ti)
721 TEveTrackList* tlist = tl[ti];
722 tlist->SetName(Form("%s [%d]", tlist->GetName(), tlist->NumChildren()));
723 tlist->SetTitle(Form("N tracks=%d", tc[ti]));
727 Bool_t good_cont = ti < 1;
728 if (AliEveTrackCounter::IsActive())
730 AliEveTrackCounter::fgInstance->RegisterTracks(tlist, good_cont);
735 tlist->SetLineStyle(6);
738 cont->SetTitle(Form("N all tracks = %d", count));
739 // ??? The following does not always work:
740 cont->FindListTreeItem(gEve->GetListTree())->SetOpen(kTRUE);