kine_tracks.C
[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
a6337352 10TEveTrack* esd_make_track(TEveTrackPropagator* trkProp,
a15e6d7d 11 Int_t index,
12 AliESDtrack* at,
13 AliExternalTrackParam* tp=0)
86f12f78 14{
a8600b56 15 // Helper function
84aff7a4 16 Double_t pbuf[3], vbuf[3];
17 TEveRecTrack rt;
86f12f78 18
19 if(tp == 0) tp = at;
20
84aff7a4 21 rt.fLabel = at->GetLabel();
22 rt.fIndex = index;
23 rt.fStatus = (Int_t) at->GetStatus();
24 rt.fSign = tp->GetSign();
86f12f78 25 tp->GetXYZ(vbuf);
84aff7a4 26 rt.fV.Set(vbuf);
86f12f78 27 tp->GetPxPyPz(pbuf);
84aff7a4 28 rt.fP.Set(pbuf);
86f12f78 29 Double_t ep = at->GetP(), mc = at->GetMass();
84aff7a4 30 rt.fBeta = ep/TMath::Sqrt(ep*ep + mc*mc);
51346b82 31
a6337352 32 TEveTrack* track = new TEveTrack(&rt, trkProp);
7be1e8cd 33 //PH The line below is replaced waiting for a fix in Root
34 //PH which permits to use variable siza arguments in CINT
35 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
84aff7a4 36 //PH track->SetName(Form("ESDTrack %d", rt.fLabel));
7be1e8cd 37 //PH track->SetTitle(Form("pT=%.3f, pZ=%.3f; V=(%.3f, %.3f, %.3f)",
84aff7a4 38 //PH rt.fSign*TMath::Hypot(rt.fP.fX, rt.fP.fY), rt.fP.fZ,
39 //PH rt.fV.fX, rt.fV.fY, rt.fV.fZ));
7be1e8cd 40 char form[1000];
84aff7a4 41 sprintf(form,"TEveTrack %d", rt.fIndex);
7be1e8cd 42 track->SetName(form);
fc9514f5 43 track->SetStdTitle();
86f12f78 44 return track;
45}
46
c222eb1d 47void esd_track_add_param(TEveTrack* track, AliExternalTrackParam* tp)
48{
49 if (tp == 0)
50 return;
51
52 Double_t pbuf[3], vbuf[3];
53 tp->GetXYZ(vbuf);
54 tp->GetPxPyPz(pbuf);
55
a15e6d7d 56 TEvePathMark pm(TEvePathMark::kReference);
57 pm.fV.Set(vbuf);
58 pm.fP.Set(pbuf);
c222eb1d 59 track->AddPathMark(pm);
60}
61
a15e6d7d 62// Use inner-tpc track params when its refit failed.
c222eb1d 63Bool_t gkFixFailedITSExtr = kFALSE;
a15e6d7d 64
65// Also show lines as generated by AliESDtrack.
66Bool_t gkMakeTrackParamLines = kFALSE;
86f12f78 67
a6337352 68TEveTrackList* esd_tracks(Double_t min_pt=0, Double_t max_pt=10000)
5a5a1232 69{
d810d0de 70 AliESDEvent* esd = AliEveEventManager::AssertESD();
5a5a1232 71
72 Double_t minptsq = min_pt*min_pt;
73 Double_t maxptsq = max_pt*max_pt;
74 Double_t ptsq;
75
51346b82 76 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3 77 cont->SetMainColor(6);
a6337352 78 TEveTrackPropagator* trkProp = cont->GetPropagator();
79 trkProp->SetMagField( 0.1*esd->GetMagneticField() );
80 trkProp->SetMaxR ( 520 );
5a5a1232 81
84aff7a4 82 gEve->AddElement(cont);
5a5a1232 83
c222eb1d 84 TEveElementList* contLines = 0;
85 if (gkMakeTrackParamLines)
86 {
87 contLines = new TEveElementList("MyTracks");
88 gEve->AddElement(contLines);
89 }
90
86f12f78 91 Int_t count = 0;
92 Double_t pbuf[3];
d31ac39d 93 for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
94 {
86f12f78 95 AliESDtrack* at = esd->GetTrack(n);
5a5a1232 96
97 // Here would be sweet to have TObjectFormula.
98 at->GetPxPyPz(pbuf);
99 ptsq = pbuf[0]*pbuf[0] + pbuf[1]*pbuf[1];
100 if(ptsq < minptsq || ptsq > maxptsq)
101 continue;
102
103 ++count;
104
a8600b56 105 // If gkFixFailedITSExtr is TRUE (FALSE by default) and
106 // if ITS refit failed, take track parameters at inner TPC radius.
c222eb1d 107 Bool_t innerTaken = kFALSE;
86f12f78 108 AliExternalTrackParam* tp = at;
109 if (gkFixFailedITSExtr && !at->IsOn(AliESDtrack::kITSrefit)) {
c222eb1d 110 tp = at->GetInnerParam();
111 innerTaken = kTRUE;
86f12f78 112 }
5a5a1232 113
a6337352 114 TEveTrack* track = esd_make_track(trkProp, n, at, tp);
32e219c2 115 track->SetAttLineAttMarker(cont);
c222eb1d 116
a15e6d7d 117 if (!innerTaken) {
118 esd_track_add_param(track, at->GetInnerParam());
119 }
c222eb1d 120 // esd_track_add_param(track, at->GetOuterParam());
121
84aff7a4 122 gEve->AddElement(track, cont);
c222eb1d 123
124 if (gkMakeTrackParamLines) {
125 TEveLine* l = new TEveLine;
126 l->SetName(Form("Track%d", count));
fbc350a3 127 l->SetLineColor(5);
c222eb1d 128 at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
129 contLines->AddElement(l);
130 }
5a5a1232 131 }
132
7be1e8cd 133 //PH The line below is replaced waiting for a fix in Root
134 //PH which permits to use variable siza arguments in CINT
135 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
136 //PH const Text_t* tooltip = Form("pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
137 char tooltip[1000];
138 sprintf(tooltip,"pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
5a5a1232 139 cont->SetTitle(tooltip); // Not broadcasted automatically ...
5a5a1232 140
141 cont->MakeTracks();
32e219c2 142
84aff7a4 143 gEve->Redraw3D();
5a5a1232 144
145 return cont;
146}
147
57ffa5fb 148/******************************************************************************/
a8600b56 149// esd_tracks_from_array()
57ffa5fb 150/******************************************************************************/
a8600b56 151
84aff7a4 152TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
5a5a1232 153{
154 // Retrieves AliESDTrack's from collection.
155 // See example usage with AliAnalysisTrackCuts in the next function.
156
d810d0de 157 if (esd == 0) esd = AliEveEventManager::AssertESD();
86f12f78 158
51346b82 159 TEveTrackList* cont = new TEveTrackList("ESD Tracks");
fbc350a3 160 cont->SetMainColor(6);
a6337352 161 TEveTrackPropagator* trkProp = cont->GetPropagator();
162 trkProp->SetMagField( 0.1*esd->GetMagneticField() );
163 trkProp->SetMaxR ( 520 );
5a5a1232 164
84aff7a4 165 gEve->AddElement(cont);
5a5a1232 166
5a5a1232 167 Int_t count = 0;
168 TIter next(col);
169 TObject *obj;
84aff7a4 170 while ((obj = next()) != 0)
d31ac39d 171 {
84aff7a4 172 if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE) {
5a5a1232 173 Warning("Object '%s', '%s' is not an AliESDtrack.",
174 obj->GetName(), obj->GetTitle());
175 continue;
176 }
177
178 ++count;
86f12f78 179 AliESDtrack* at = (AliESDtrack*) obj;
5a5a1232 180
a6337352 181 TEveTrack* track = esd_make_track(trkProp, count, at);
32e219c2 182 track->SetAttLineAttMarker(cont);
84aff7a4 183 gEve->AddElement(track, cont);
5a5a1232 184 }
185
7be1e8cd 186 //PH The line below is replaced waiting for a fix in Root
187 //PH which permits to use variable siza arguments in CINT
188 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
189 //PH const Text_t* tooltip = Form("N=%d", count);
190 const tooltip[1000];
191 sprintf(tooltip,"N=%d", count);
5a5a1232 192 cont->SetTitle(tooltip); // Not broadcasted automatically ...
5a5a1232 193
194 cont->MakeTracks();
32e219c2 195
84aff7a4 196 gEve->Redraw3D();
5a5a1232 197
198 return cont;
199}
200
201void esd_tracks_alianalcuts_demo()
202{
d810d0de 203 AliESDEvent* esd = AliEveEventManager::AssertESD();
86f12f78 204 gSystem->Load("libANALYSIS");
5a5a1232 205
206 AliAnalysisTrackCuts atc;
207 atc.SetPtRange(0.1, 5);
208 atc.SetRapRange(-1, 1);
209
86f12f78 210 esd_tracks_from_array(atc.GetAcceptedParticles(esd), esd);
5a5a1232 211}
a8600b56 212
57ffa5fb 213/******************************************************************************/
a8600b56 214// esd_tracks_vertex_cut
57ffa5fb 215/******************************************************************************/
a8600b56 216
217Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
218{
219 // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
220 // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
221
222 Float_t b[2];
223 Float_t bRes[2];
224 Float_t bCov[3];
225 esdTrack->GetImpactParameters(b,bCov);
226 if (bCov[0]<=0 || bCov[2]<=0) {
d31ac39d 227 printf("Estimated b resolution lower or equal zero!\n");
a8600b56 228 bCov[0]=0; bCov[2]=0;
229 }
230 bRes[0] = TMath::Sqrt(bCov[0]);
231 bRes[1] = TMath::Sqrt(bCov[2]);
232
233 // -----------------------------------
234 // How to get to a n-sigma cut?
235 //
236 // The accumulated statistics from 0 to d is
237 //
238 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
239 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
240 //
241 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
242 // Can this be expressed in a different way?
243
244 if (bRes[0] == 0 || bRes[1] ==0)
245 return -1;
246
247 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
248
249 // stupid rounding problem screws up everything:
250 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
251 if (TMath::Exp(-d * d / 2) < 1e-10)
252 return 1000;
253
254 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
255 return d;
256}
257
84aff7a4 258TEveElementList* esd_tracks_vertex_cut()
a8600b56 259{
f206464b 260 // Import ESD tracks, separate them into five containers according to
261 // primary-vertex cut and ITS refit status.
a8600b56 262
d810d0de 263 AliESDEvent* esd = AliEveEventManager::AssertESD();
a8600b56 264
84aff7a4 265 TEveElementList* cont = new TEveElementList("ESD Tracks", 0, kTRUE);
266 gEve->AddElement(cont);
267 TEveTrackList *tl[5];
f206464b 268 Int_t tc[5];
a8600b56 269 Int_t count = 0;
270
84aff7a4 271 tl[0] = new TEveTrackList("Sigma < 3");
a8600b56 272 tc[0] = 0;
daaa6c4d 273 tl[0]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
a6337352 274 tl[0]->GetPropagator()->SetMaxR ( 520 );
fbc350a3 275 tl[0]->SetMainColor(3);
84aff7a4 276 gEve->AddElement(tl[0], cont);
a8600b56 277
84aff7a4 278 tl[1] = new TEveTrackList("3 < Sigma < 5");
a8600b56 279 tc[1] = 0;
daaa6c4d 280 tl[1]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
a6337352 281 tl[1]->GetPropagator()->SetMaxR ( 520 );
fbc350a3 282 tl[1]->SetMainColor(7);
84aff7a4 283 gEve->AddElement(tl[1], cont);
a8600b56 284
84aff7a4 285 tl[2] = new TEveTrackList("5 < Sigma");
a8600b56 286 tc[2] = 0;
daaa6c4d 287 tl[2]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
a6337352 288 tl[2]->GetPropagator()->SetMaxR ( 520 );
fbc350a3 289 tl[2]->SetMainColor(46);
84aff7a4 290 gEve->AddElement(tl[2], cont);
a8600b56 291
84aff7a4 292 tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
a8600b56 293 tc[3] = 0;
daaa6c4d 294 tl[3]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
a6337352 295 tl[3]->GetPropagator()->SetMaxR ( 520 );
fbc350a3 296 tl[3]->SetMainColor(41);
84aff7a4 297 gEve->AddElement(tl[3], cont);
a8600b56 298
84aff7a4 299 tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
f206464b 300 tc[4] = 0;
daaa6c4d 301 tl[4]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
a6337352 302 tl[4]->GetPropagator()->SetMaxR ( 520 );
fbc350a3 303 tl[4]->SetMainColor(48);
84aff7a4 304 gEve->AddElement(tl[4], cont);
f206464b 305
d31ac39d 306 for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
307 {
a8600b56 308 AliESDtrack* at = esd->GetTrack(n);
309
310 Float_t s = get_sigma_to_vertex(at);
311 Int_t ti;
312 if (s < 3) ti = 0;
313 else if (s <= 5) ti = 1;
314 else ti = 2;
315
316 AliExternalTrackParam* tp = at;
f206464b 317 // If ITS refit failed, optionally take track parameters at inner
318 // TPC radius and put track in a special container.
a8600b56 319 // This ignores state of gkFixFailedITSExtr (used in esd_tracks()).
320 // Use BOTH functions to compare results.
321 if (!at->IsOn(AliESDtrack::kITSrefit)) {
f206464b 322 // tp = at->GetInnerParam();
323 ti = (ti == 2) ? 4 : 3;
a8600b56 324 }
325
84aff7a4 326 TEveTrackList* tlist = tl[ti];
a8600b56 327 ++tc[ti];
328 ++count;
329
84aff7a4 330 TEveTrack* track = esd_make_track(tlist->GetPropagator(), n, at, tp);
51346b82 331 track->SetAttLineAttMarker(tlist);
d31ac39d 332
7be1e8cd 333 //PH The line below is replaced waiting for a fix in Root
334 //PH which permits to use variable siza arguments in CINT
335 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
336 //PH track->SetName(Form("track %d, sigma=%5.3f", at->GetLabel(), s));
337 char form[1000];
84aff7a4 338 sprintf(form,"TEveTrack idx=%d, sigma=%5.3f", at->GetID(), s);
7be1e8cd 339 track->SetName(form);
84aff7a4 340 gEve->AddElement(track, tlist);
a8600b56 341 }
342
f206464b 343 for (Int_t ti=0; ti<5; ++ti) {
84aff7a4 344 TEveTrackList* tlist = tl[ti];
7be1e8cd 345 //PH The line below is replaced waiting for a fix in Root
346 //PH which permits to use variable siza arguments in CINT
347 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
348 //PH const Text_t* tooltip = Form("N tracks=%d", tc[ti]);
846b5c55 349 //MT Modified somewhat.
350 char buff[1000];
fbc350a3 351 sprintf(buff, "%s [%d]", tlist->GetName(), tlist->NumChildren());
846b5c55 352 tlist->SetName(buff);
353 sprintf(buff, "N tracks=%d", tc[ti]);
354 tlist->SetTitle(buff); // Not broadcasted automatically ...
a8600b56 355
356 tlist->MakeTracks();
a8600b56 357 }
7be1e8cd 358 //PH The line below is replaced waiting for a fix in Root
359 //PH which permits to use variable siza arguments in CINT
360 //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
361 //PH cont->SetTitle(Form("N all tracks = %d", count));
362 char form[1000];
363 sprintf(form,"N all tracks = %d", count);
364 cont->SetTitle(form);
84aff7a4 365 gEve->Redraw3D();
a8600b56 366
367 return cont;
368}