Follow rename of method in AliESDtrack.
[u/mrichter/AliRoot.git] / EVE / alice-macros / esd_tracks.C
1 // $Id$
2
3 Reve::Track* esd_make_track(Reve::TrackRnrStyle*   rnrStyle,
4                             Int_t                  index,
5                             AliESDtrack*           at,
6                             AliExternalTrackParam* tp=0)
7 {
8   // Helper function
9   Double_t        pbuf[3], vbuf[3];
10   Reve::RecTrack  rt;
11
12   if(tp == 0) tp = at;
13
14   rt.label  = at->GetLabel();
15   rt.index  = index;
16   rt.status = (Int_t) at->GetStatus();
17   rt.sign   = tp->GetSign();
18   tp->GetXYZ(vbuf);
19   rt.V.Set(vbuf);
20   tp->GetPxPyPz(pbuf);
21   rt.P.Set(pbuf);
22   Double_t ep = at->GetP(), mc = at->GetMass();
23   rt.beta = ep/TMath::Sqrt(ep*ep + mc*mc);
24  
25   Reve::Track* track = new Reve::Track(&rt, rnrStyle);
26   //PH The line below is replaced waiting for a fix in Root
27   //PH which permits to use variable siza arguments in CINT
28   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
29   //PH  track->SetName(Form("ESDTrack %d", rt.label));
30   //PH  track->SetTitle(Form("pT=%.3f, pZ=%.3f; V=(%.3f, %.3f, %.3f)",
31   //PH                 rt.sign*TMath::Hypot(rt.P.x, rt.P.y), rt.P.z,
32   //PH                 rt.V.x, rt.V.y, rt.V.z));
33   char form[1000];
34   sprintf(form,"Track %d", rt.index);
35   track->SetName(form);
36   track->SetStdTitle();
37   return track;
38 }
39
40 Bool_t gkFixFailedITSExtr = kFALSE;
41
42 Reve::TrackList* esd_tracks(Double_t min_pt=0.1, Double_t max_pt=100)
43 {
44   AliESDEvent* esd = Alieve::Event::AssertESD();
45
46   Double_t minptsq = min_pt*min_pt;
47   Double_t maxptsq = max_pt*max_pt;
48   Double_t ptsq;
49
50   Reve::TrackList* cont = new Reve::TrackList("ESD Tracks"); 
51   cont->SetMainColor(Color_t(6));
52   Reve::TrackRnrStyle* rnrStyle = cont->GetRnrStyle();
53   rnrStyle->SetMagField( esd->GetMagneticField() );
54
55   gReve->AddRenderElement(cont);
56
57   Int_t    count = 0;
58   Double_t pbuf[3];
59   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
60   {
61     AliESDtrack* at = esd->GetTrack(n);
62
63     // Here would be sweet to have TObjectFormula.
64     at->GetPxPyPz(pbuf);
65     ptsq = pbuf[0]*pbuf[0] + pbuf[1]*pbuf[1];
66     if(ptsq < minptsq || ptsq > maxptsq)
67       continue;
68
69     ++count;
70
71     // If gkFixFailedITSExtr is TRUE (FALSE by default) and
72     // if ITS refit failed, take track parameters at inner TPC radius.
73     AliExternalTrackParam* tp = at;
74     if (gkFixFailedITSExtr && !at->IsOn(AliESDtrack::kITSrefit)) {
75        tp = at->GetInnerParam();
76     }
77
78     Reve::Track* track = esd_make_track(rnrStyle, n, at, tp);
79     track->SetAttLineAttMarker(cont);
80     gReve->AddRenderElement(track, cont);
81   }
82
83   //PH The line below is replaced waiting for a fix in Root
84   //PH which permits to use variable siza arguments in CINT
85   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
86   //PH  const Text_t* tooltip = Form("pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
87   char tooltip[1000];
88   sprintf(tooltip,"pT ~ (%.2lf, %.2lf), N=%d", min_pt, max_pt, count);
89   cont->SetTitle(tooltip); // Not broadcasted automatically ...
90   cont->UpdateItems();
91
92   cont->MakeTracks();
93
94   gReve->Redraw3D();
95
96   return cont;
97 }
98
99 /**************************************************************************/
100 // esd_tracks_from_array()
101 /**************************************************************************/
102
103 Reve::TrackList* esd_tracks_from_array(TCollection* col, AliESDEvent* esd=0)
104 {
105   // Retrieves AliESDTrack's from collection.
106   // See example usage with AliAnalysisTrackCuts in the next function.
107
108   if(esd == 0) esd = Alieve::Event::AssertESD();
109
110   Reve::TrackList* cont = new Reve::TrackList("ESD Tracks"); 
111   cont->SetMainColor(Color_t(6));
112   Reve::TrackRnrStyle* rnrStyle = cont->GetRnrStyle();
113   rnrStyle->SetMagField( esd->GetMagneticField() );
114
115   gReve->AddRenderElement(cont);
116
117   Int_t    count = 0;
118   TIter    next(col);
119   TObject *obj;
120   while((obj = next()) != 0)
121   {
122     if(obj->IsA()->InheritsFrom("AliESDtrack") == kFALSE) {
123       Warning("Object '%s', '%s' is not an AliESDtrack.",
124               obj->GetName(), obj->GetTitle());
125       continue;
126     }
127
128     ++count;
129     AliESDtrack* at = (AliESDtrack*) obj;
130
131     Reve::Track* track = esd_make_track(rnrStyle, count, at);
132     track->SetAttLineAttMarker(cont);
133     gReve->AddRenderElement(track, cont);
134   }
135
136   //PH The line below is replaced waiting for a fix in Root
137   //PH which permits to use variable siza arguments in CINT
138   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
139   //PH  const Text_t* tooltip = Form("N=%d", count);
140   const tooltip[1000];
141   sprintf(tooltip,"N=%d", count);
142   cont->SetTitle(tooltip); // Not broadcasted automatically ...
143   cont->UpdateItems();
144
145   cont->MakeTracks();
146
147   gReve->Redraw3D();
148
149   return cont;
150 }
151
152 void esd_tracks_alianalcuts_demo()
153 {
154   AliESDEvent* esd = Alieve::Event::AssertESD();
155   gSystem->Load("libANALYSIS");
156
157   AliAnalysisTrackCuts atc;
158   atc.SetPtRange(0.1, 5);
159   atc.SetRapRange(-1, 1);
160
161   esd_tracks_from_array(atc.GetAcceptedParticles(esd), esd);
162 }
163
164 /**************************************************************************/
165 // esd_tracks_vertex_cut
166 /**************************************************************************/
167
168 Float_t get_sigma_to_vertex(AliESDtrack* esdTrack)
169 {
170   // Taken from: PWG0/esdTrackCuts/AliESDtrackCuts.cxx
171   // Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
172
173   Float_t b[2];
174   Float_t bRes[2];
175   Float_t bCov[3];
176   esdTrack->GetImpactParameters(b,bCov);
177   if (bCov[0]<=0 || bCov[2]<=0) {
178     printf("Estimated b resolution lower or equal zero!\n");
179     bCov[0]=0; bCov[2]=0;
180   }
181   bRes[0] = TMath::Sqrt(bCov[0]);
182   bRes[1] = TMath::Sqrt(bCov[2]);
183
184   // -----------------------------------
185   // How to get to a n-sigma cut?
186   //
187   // The accumulated statistics from 0 to d is
188   //
189   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
190   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
191   //
192   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
193   // Can this be expressed in a different way?
194
195   if (bRes[0] == 0 || bRes[1] ==0)
196     return -1;
197
198   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
199
200   // stupid rounding problem screws up everything:
201   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
202   if (TMath::Exp(-d * d / 2) < 1e-10)
203     return 1000;
204
205   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
206   return d;
207 }
208
209 Reve::RenderElementList* esd_tracks_vertex_cut()
210 {
211   // Import ESD tracks, separate them into five containers according to
212   // primary-vertex cut and ITS refit status.
213
214   AliESDEvent* esd = Alieve::Event::AssertESD();
215
216   Reve::RenderElementList* cont = new Reve::RenderElementList("ESD Tracks", 0, kTRUE);
217   gReve->AddRenderElement(cont);
218   Reve::TrackList *tl[5];
219   Int_t            tc[5];
220   Int_t            count = 0;
221
222   tl[0] = new Reve::TrackList("Sigma < 3");
223   tc[0] = 0;
224   tl[0]->GetRnrStyle()->SetMagField( esd->GetMagneticField() );
225   tl[0]->SetMainColor(Color_t(3));
226   gReve->AddRenderElement(tl[0], cont);
227
228   tl[1] = new Reve::TrackList("3 < Sigma < 5");
229   tc[1] = 0;
230   tl[1]->GetRnrStyle()->SetMagField( esd->GetMagneticField() );
231   tl[1]->SetMainColor(Color_t(7));
232   gReve->AddRenderElement(tl[1], cont);
233
234   tl[2] = new Reve::TrackList("5 < Sigma");
235   tc[2] = 0;
236   tl[2]->GetRnrStyle()->SetMagField( esd->GetMagneticField() );
237   tl[2]->SetMainColor(Color_t(46));
238   gReve->AddRenderElement(tl[2], cont);
239
240   tl[3] = new Reve::TrackList("no ITS refit; Sigma < 5");
241   tc[3] = 0;
242   tl[3]->GetRnrStyle()->SetMagField( esd->GetMagneticField() );
243   tl[3]->SetMainColor(Color_t(41));
244   gReve->AddRenderElement(tl[3], cont);
245
246   tl[4] = new Reve::TrackList("no ITS refit; Sigma > 5");
247   tc[4] = 0;
248   tl[4]->GetRnrStyle()->SetMagField( esd->GetMagneticField() );
249   tl[4]->SetMainColor(Color_t(48));
250   gReve->AddRenderElement(tl[4], cont);
251
252   for (Int_t n=0; n<esd->GetNumberOfTracks(); n++)
253   {
254     AliESDtrack* at = esd->GetTrack(n);
255
256     Float_t s  = get_sigma_to_vertex(at);
257     Int_t   ti;
258     if      (s <  3) ti = 0;
259     else if (s <= 5) ti = 1;
260     else             ti = 2;
261
262     AliExternalTrackParam* tp = at;
263     // If ITS refit failed, optionally take track parameters at inner
264     // TPC radius and put track in a special container.
265     // This ignores state of gkFixFailedITSExtr (used in esd_tracks()).
266     // Use BOTH functions to compare results.
267     if (!at->IsOn(AliESDtrack::kITSrefit)) {
268       // tp = at->GetInnerParam();
269       ti = (ti == 2) ? 4 : 3;
270     }
271
272     Reve::TrackList* tlist = tl[ti];
273     ++tc[ti];
274     ++count;
275
276     Reve::Track* track = esd_make_track(tlist->GetRnrStyle(), n, at, tp);
277     track->SetAttLineAttMarker(tlist);    
278
279     //PH The line below is replaced waiting for a fix in Root
280     //PH which permits to use variable siza arguments in CINT
281     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
282     //PH    track->SetName(Form("track %d, sigma=%5.3f", at->GetLabel(), s));
283     char form[1000];
284     sprintf(form,"Track idx=%d, sigma=%5.3f", at->GetID(), s);
285     track->SetName(form);
286     gReve->AddRenderElement(track, tlist);
287   }
288
289   for (Int_t ti=0; ti<5; ++ti) {
290     Reve::TrackList* tlist = tl[ti];
291     //PH The line below is replaced waiting for a fix in Root
292     //PH which permits to use variable siza arguments in CINT
293     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
294     //PH    const Text_t* tooltip = Form("N tracks=%d", tc[ti]);
295     //MT Modified somewhat.
296     char buff[1000];
297     sprintf(buff, "%s [%d]", tlist->GetName(), tlist->GetNChildren());
298     tlist->SetName(buff);
299     sprintf(buff, "N tracks=%d", tc[ti]);
300     tlist->SetTitle(buff); // Not broadcasted automatically ...
301     tlist->UpdateItems();
302
303     tlist->MakeTracks();
304   }
305   //PH The line below is replaced waiting for a fix in Root
306   //PH which permits to use variable siza arguments in CINT
307   //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
308   //PH  cont->SetTitle(Form("N all tracks = %d", count));
309   char form[1000];
310   sprintf(form,"N all tracks = %d", count);
311   cont->SetTitle(form);
312   cont->UpdateItems();
313   gReve->Redraw3D();
314
315   return cont;
316 }