kine_tracks.C
[u/mrichter/AliRoot.git] / EVE / alice-macros / esd_tracks.C
1 // $Id$
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          *
7  * full copyright notice.                                                 *
8  **************************************************************************/
9
10 TEveTrack* esd_make_track(TEveTrackPropagator*   trkProp,
11                           Int_t                  index,
12                           AliESDtrack*           at,
13                           AliExternalTrackParam* tp=0)
14 {
15   // Helper function
16   Double_t      pbuf[3], vbuf[3];
17   TEveRecTrack  rt;
18
19   if(tp == 0) tp = at;
20
21   rt.fLabel  = at->GetLabel();
22   rt.fIndex  = index;
23   rt.fStatus = (Int_t) at->GetStatus();
24   rt.fSign   = tp->GetSign();
25   tp->GetXYZ(vbuf);
26   rt.fV.Set(vbuf);
27   tp->GetPxPyPz(pbuf);
28   rt.fP.Set(pbuf);
29   Double_t ep = at->GetP(), mc = at->GetMass();
30   rt.fBeta = ep/TMath::Sqrt(ep*ep + mc*mc);
31
32   TEveTrack* track = new TEveTrack(&rt, trkProp);
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.)
36   //PH  track->SetName(Form("ESDTrack %d", rt.fLabel));
37   //PH  track->SetTitle(Form("pT=%.3f, pZ=%.3f; V=(%.3f, %.3f, %.3f)",
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));
40   char form[1000];
41   sprintf(form,"TEveTrack %d", rt.fIndex);
42   track->SetName(form);
43   track->SetStdTitle();
44   return track;
45 }
46
47 void 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
56   TEvePathMark pm(TEvePathMark::kReference);
57   pm.fV.Set(vbuf);
58   pm.fP.Set(pbuf);
59   track->AddPathMark(pm);
60 }
61
62 // Use inner-tpc track params when its refit failed.
63 Bool_t gkFixFailedITSExtr    = kFALSE;
64
65 // Also show lines as generated by AliESDtrack.
66 Bool_t gkMakeTrackParamLines = kFALSE;
67
68 TEveTrackList* esd_tracks(Double_t min_pt=0, Double_t max_pt=10000)
69 {
70   AliESDEvent* esd = AliEveEventManager::AssertESD();
71
72   Double_t minptsq = min_pt*min_pt;
73   Double_t maxptsq = max_pt*max_pt;
74   Double_t ptsq;
75
76   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
77   cont->SetMainColor(6);
78   TEveTrackPropagator* trkProp = cont->GetPropagator();
79   trkProp->SetMagField( 0.1*esd->GetMagneticField() );
80   trkProp->SetMaxR    ( 520 );
81
82   gEve->AddElement(cont);
83
84   TEveElementList* contLines = 0;
85   if (gkMakeTrackParamLines)
86   {
87     contLines = new TEveElementList("MyTracks");  
88     gEve->AddElement(contLines);
89   }
90
91   Int_t    count = 0;
92   Double_t pbuf[3];
93   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
94   {
95     AliESDtrack* at = esd->GetTrack(n);
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
105     // If gkFixFailedITSExtr is TRUE (FALSE by default) and
106     // if ITS refit failed, take track parameters at inner TPC radius.
107     Bool_t innerTaken = kFALSE;
108     AliExternalTrackParam* tp = at;
109     if (gkFixFailedITSExtr && !at->IsOn(AliESDtrack::kITSrefit)) {
110       tp = at->GetInnerParam();
111       innerTaken = kTRUE;
112     }
113
114     TEveTrack* track = esd_make_track(trkProp, n, at, tp);
115     track->SetAttLineAttMarker(cont);
116
117     if (!innerTaken) {
118       esd_track_add_param(track, at->GetInnerParam());
119     }
120     // esd_track_add_param(track, at->GetOuterParam());
121
122     gEve->AddElement(track, cont);
123
124     if (gkMakeTrackParamLines) {
125       TEveLine* l = new TEveLine; 
126       l->SetName(Form("Track%d", count));
127       l->SetLineColor(5);
128       at->FillPolymarker(l, esd->GetMagneticField(), 0, 250, 5);
129       contLines->AddElement(l);
130     }
131   }
132
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);
139   cont->SetTitle(tooltip); // Not broadcasted automatically ...
140
141   cont->MakeTracks();
142
143   gEve->Redraw3D();
144
145   return cont;
146 }
147
148 /******************************************************************************/
149 // esd_tracks_from_array()
150 /******************************************************************************/
151
152 TEveTrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
153 {
154   // Retrieves AliESDTrack's from collection.
155   // See example usage with AliAnalysisTrackCuts in the next function.
156
157   if (esd == 0) esd = AliEveEventManager::AssertESD();
158
159   TEveTrackList* cont = new TEveTrackList("ESD Tracks");
160   cont->SetMainColor(6);
161   TEveTrackPropagator* trkProp = cont->GetPropagator();
162   trkProp->SetMagField( 0.1*esd->GetMagneticField() );
163   trkProp->SetMaxR    ( 520 );
164
165   gEve->AddElement(cont);
166
167   Int_t    count = 0;
168   TIter    next(col);
169   TObject *obj;
170   while ((obj = next()) != 0)
171   {
172     if (obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE) {
173       Warning("Object '%s', '%s' is not an AliESDtrack.",
174               obj->GetName(), obj->GetTitle());
175       continue;
176     }
177
178     ++count;
179     AliESDtrack* at = (AliESDtrack*) obj;
180
181     TEveTrack* track = esd_make_track(trkProp, count, at);
182     track->SetAttLineAttMarker(cont);
183     gEve->AddElement(track, cont);
184   }
185
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);
192   cont->SetTitle(tooltip); // Not broadcasted automatically ...
193
194   cont->MakeTracks();
195
196   gEve->Redraw3D();
197
198   return cont;
199 }
200
201 void esd_tracks_alianalcuts_demo()
202 {
203   AliESDEvent* esd = AliEveEventManager::AssertESD();
204   gSystem->Load("libANALYSIS");
205
206   AliAnalysisTrackCuts atc;
207   atc.SetPtRange(0.1, 5);
208   atc.SetRapRange(-1, 1);
209
210   esd_tracks_from_array(atc.GetAcceptedParticles(esd), esd);
211 }
212
213 /******************************************************************************/
214 // esd_tracks_vertex_cut
215 /******************************************************************************/
216
217 Float_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) {
227     printf("Estimated b resolution lower or equal zero!\n");
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
258 TEveElementList* esd_tracks_vertex_cut()
259 {
260   // Import ESD tracks, separate them into five containers according to
261   // primary-vertex cut and ITS refit status.
262
263   AliESDEvent* esd = AliEveEventManager::AssertESD();
264
265   TEveElementList* cont = new TEveElementList("ESD Tracks", 0, kTRUE);
266   gEve->AddElement(cont);
267   TEveTrackList *tl[5];
268   Int_t            tc[5];
269   Int_t            count = 0;
270
271   tl[0] = new TEveTrackList("Sigma < 3");
272   tc[0] = 0;
273   tl[0]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
274   tl[0]->GetPropagator()->SetMaxR    ( 520 );
275   tl[0]->SetMainColor(3);
276   gEve->AddElement(tl[0], cont);
277
278   tl[1] = new TEveTrackList("3 < Sigma < 5");
279   tc[1] = 0;
280   tl[1]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
281   tl[1]->GetPropagator()->SetMaxR    ( 520 );
282   tl[1]->SetMainColor(7);
283   gEve->AddElement(tl[1], cont);
284
285   tl[2] = new TEveTrackList("5 < Sigma");
286   tc[2] = 0;
287   tl[2]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
288   tl[2]->GetPropagator()->SetMaxR    ( 520 );
289   tl[2]->SetMainColor(46);
290   gEve->AddElement(tl[2], cont);
291
292   tl[3] = new TEveTrackList("no ITS refit; Sigma < 5");
293   tc[3] = 0;
294   tl[3]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
295   tl[3]->GetPropagator()->SetMaxR    ( 520 );
296   tl[3]->SetMainColor(41);
297   gEve->AddElement(tl[3], cont);
298
299   tl[4] = new TEveTrackList("no ITS refit; Sigma > 5");
300   tc[4] = 0;
301   tl[4]->GetPropagator()->SetMagField( 0.1*esd->GetMagneticField() );
302   tl[4]->GetPropagator()->SetMaxR    ( 520 );
303   tl[4]->SetMainColor(48);
304   gEve->AddElement(tl[4], cont);
305
306   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
307   {
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;
317     // If ITS refit failed, optionally take track parameters at inner
318     // TPC radius and put track in a special container.
319     // This ignores state of gkFixFailedITSExtr (used in esd_tracks()).
320     // Use BOTH functions to compare results.
321     if (!at->IsOn(AliESDtrack::kITSrefit)) {
322       // tp = at->GetInnerParam();
323       ti = (ti == 2) ? 4 : 3;
324     }
325
326     TEveTrackList* tlist = tl[ti];
327     ++tc[ti];
328     ++count;
329
330     TEveTrack* track = esd_make_track(tlist->GetPropagator(), n, at, tp);
331     track->SetAttLineAttMarker(tlist);
332
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];
338     sprintf(form,"TEveTrack idx=%d, sigma=%5.3f", at->GetID(), s);
339     track->SetName(form);
340     gEve->AddElement(track, tlist);
341   }
342
343   for (Int_t ti=0; ti<5; ++ti) {
344     TEveTrackList* tlist = tl[ti];
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]);
349     //MT Modified somewhat.
350     char buff[1000];
351     sprintf(buff, "%s [%d]", tlist->GetName(), tlist->NumChildren());
352     tlist->SetName(buff);
353     sprintf(buff, "N tracks=%d", tc[ti]);
354     tlist->SetTitle(buff); // Not broadcasted automatically ...
355
356     tlist->MakeTracks();
357   }
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);
365   gEve->Redraw3D();
366
367   return cont;
368 }