]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDcheckESD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fixed bugs with the rule checker.
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckESD.cxx
CommitLineData
695019eb 1/**************************************************************************\r
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3* *\r
4* Author: The ALICE Off-line Project. *\r
5* Contributors are mentioned in the code where appropriate. *\r
6* *\r
7* Permission to use, copy, modify and distribute this software and its *\r
8* documentation strictly for non-commercial purposes is hereby granted *\r
9* without fee, provided that the above copyright notice appears in all *\r
10* copies and that both the copyright notice and this permission notice *\r
11* appear in the supporting documentation. The authors make no claims *\r
12* about the suitability of this software for any purpose. It is *\r
13* provided "as is" without express or implied warranty. *\r
14**************************************************************************/\r
15\r
16/////////////////////////////////////////////////////\r
17//\r
18// Check basic detector results at ESD level\r
19// - Geometrical efficiency \r
20// - Tracking efficiency \r
21// - PID efficiency \r
22// - Refit efficiency \r
23//\r
24// Author\r
25// Alex Bercuci <A.Bercuci@gsi.de>\r
26//\r
27//////////////////////////////////////////////////////\r
28\r
29#include <TClonesArray.h>\r
30#include <TCanvas.h>\r
31#include <TObjArray.h>\r
32#include <TPad.h>\r
33#include <TLegend.h>\r
34#include <TLatex.h>\r
35#include <TLine.h>\r
36#include <TF1.h>\r
37#include <TH2I.h>\r
38#include <TH2F.h>\r
39#include <TH3S.h>\r
40#include <TH3F.h>\r
41#include <TProfile2D.h>\r
42#include <TProfile.h>\r
43#include <TGraphErrors.h>\r
44#include <TGraphAsymmErrors.h>\r
45#include <TFile.h>\r
46#include <TTree.h>\r
47#include <TROOT.h>\r
48#include <TChain.h>\r
49#include <TParticle.h>\r
50#include <TTimeStamp.h>\r
51\r
52#include "AliLog.h"\r
53#include "AliAnalysisManager.h"\r
54#include "AliESDEvent.h"\r
55#include "AliESDkink.h"\r
56#include "AliMCEvent.h"\r
57#include "AliESDInputHandler.h"\r
58#include "AliMCEventHandler.h"\r
59\r
60#include "AliESDtrack.h"\r
61#include "AliMCParticle.h"\r
62#include "AliPID.h"\r
63#include "AliStack.h"\r
64#include "AliTrackReference.h"\r
65\r
66#include "AliTRDcheckESD.h"\r
67\r
68ClassImp(AliTRDcheckESD)\r
69\r
70const Float_t AliTRDcheckESD::fgkxTPC = 290.;\r
71const Float_t AliTRDcheckESD::fgkxTOF = 365.;\r
72const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={\r
738, 4, 2, 20};\r
74FILE* AliTRDcheckESD::fgFile = NULL;\r
75\r
76const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;\r
77const Int_t AliTRDcheckESD::fgkEvVertexN = 1;\r
78const Float_t AliTRDcheckESD::fgkTrkDCAxy = 40.;\r
79const Float_t AliTRDcheckESD::fgkTrkDCAz = 15.;\r
80const Int_t AliTRDcheckESD::fgkNclTPC = 100;\r
81const Float_t AliTRDcheckESD::fgkPt = 0.2;\r
82const Float_t AliTRDcheckESD::fgkEta = 0.9;\r
83const Float_t AliTRDcheckESD::fgkQs = 0.002;\r
84\r
85//____________________________________________________________________\r
86AliTRDcheckESD::AliTRDcheckESD():\r
87 AliAnalysisTaskSE()\r
88 ,fStatus(0)\r
89 ,fNRefFigures(0)\r
90 ,fESD(NULL)\r
91 ,fMC(NULL)\r
92 ,fHistos(NULL)\r
93 ,fResults(NULL)\r
94{\r
95 //\r
96 // Default constructor\r
97 //\r
76106bcc 98 SetNameTitle("TRDcheckESD", "Check TRD @ ESD level");\r
695019eb 99 SetMC(kTRUE);\r
100}\r
101\r
102//____________________________________________________________________\r
103AliTRDcheckESD::AliTRDcheckESD(char* name):\r
104 AliAnalysisTaskSE(name)\r
105 ,fStatus(0)\r
106 ,fNRefFigures(0)\r
107 ,fESD(NULL)\r
108 ,fMC(NULL)\r
109 ,fHistos(NULL)\r
110 ,fResults(NULL)\r
111{\r
112 //\r
113 // Default constructor\r
114 //\r
115 SetMC(kTRUE);\r
116 SetTitle("Check TRD @ ESD level");\r
117 DefineOutput(1, TObjArray::Class());\r
118}\r
119\r
120//____________________________________________________________________\r
121AliTRDcheckESD::~AliTRDcheckESD()\r
122{\r
123// Destructor\r
124 if(fHistos){\r
125 //fHistos->Delete();\r
126 delete fHistos;\r
127 }\r
128 if(fResults){\r
129 fResults->Delete();\r
130 delete fResults;\r
131 }\r
132}\r
133\r
134//____________________________________________________________________\r
135void AliTRDcheckESD::UserCreateOutputObjects()\r
136{ \r
137 //\r
138 // Create Output Containers (TObjectArray containing 1D histograms)\r
139 //\r
140 Histos();\r
068e2c00 141 PostData(1, fHistos);\r
695019eb 142}\r
143\r
144//____________________________________________________________________\r
145void AliTRDcheckESD::MakeSummary(){\r
146 TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);\r
147 cOut->cd();\r
148 GetRefFigure(4);\r
149 cOut->SaveAs(Form("TRDsummary%s1.gif", GetName()));\r
150\r
151 cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 768);\r
152 cOut->cd();\r
153 GetRefFigure(5);\r
154 cOut->SaveAs(Form("TRDsummary%s2.gif", GetName()));\r
155}\r
156\r
157//____________________________________________________________________\r
158Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)\r
159{\r
160 if(ifig>=fNRefFigures){\r
161 AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));\r
162 return kFALSE;\r
163 }\r
164 if(!gPad){\r
165 AliWarning("Please provide a canvas to draw results.");\r
166 return kFALSE;\r
167 } else {\r
168 gPad->SetLogx(0);gPad->SetLogy(0);\r
169 gPad->SetMargin(0.125, 0.015, 0.1, 0.015);\r
170 }\r
171\r
172 const Char_t *title[20];\r
173 Float_t nada(0.0);\r
174 TH1 *hF(NULL);\r
175 TH1 *hFeffP(NULL); TH1 *hFeffN(NULL);\r
176 TH2 *h2F(NULL); TH2 *h2Feff(NULL);\r
177 TH2 *h2FtpcP(NULL); TH2 *h2FtpcN(NULL);\r
178 TH2 *h2FtrdP(NULL); TH2 *h2FtrdN(NULL);\r
179 TH3 *h3F(NULL);\r
180 if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;\r
181 TLegend *leg(NULL);\r
182 TList *l(NULL); TVirtualPad *pad(NULL);\r
183 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);\r
184 TObjArray *arr(NULL);\r
185 TProfile2D *hProf2D(NULL);\r
186 TProfile *hProf(NULL);\r
187 TLatex *lat=new TLatex();\r
188 lat->SetTextSize(0.07);\r
189 lat->SetTextColor(2);\r
190 TLine line;\r
191 TTimeStamp now;\r
192 TF1* fitFunc(NULL);\r
193 switch(ifig){\r
194 case kNCl: // number of clusters/track\r
195 if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;\r
196\r
197 leg = new TLegend(.83, .7, .99, .96);\r
198 leg->SetHeader("Species");\r
199 leg->SetBorderSize(0); leg->SetFillStyle(0);\r
200 for(Int_t ig(0); ig<fgkNgraph[kNCl]; ig++){\r
201 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;\r
202 if(!g->GetN()) continue;\r
203 g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");\r
204 if(ig) continue;\r
205 hF=g->GetHistogram();\r
206 hF->SetXTitle("no of clusters");\r
207 hF->SetYTitle("entries"); \r
208 hF->GetYaxis()->CenterTitle(1);\r
209 hF->GetYaxis()->SetTitleOffset(1.2);\r
210 hF->SetMinimum(5);\r
211 }\r
212 leg->Draw(); gPad->SetLogy();\r
213 break;\r
214 case kTRDstat: // Efficiency\r
215 if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;\r
216 leg = new TLegend(.62, .77, .98, .98);\r
217 leg->SetHeader("TRD Efficiency");\r
218 leg->SetBorderSize(0); leg->SetFillStyle(0);\r
219 title[0] = "Geometrical (TRDin/TPCout)";\r
220 title[1] = "Tracking (TRDout/TRDin)";\r
221 title[2] = "PID (TRDpid/TRDin)";\r
222 title[3] = "Refit (TRDrefit/TRDin)";\r
223 hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);\r
224 hF->SetMaximum(1.4);\r
225 hF->GetXaxis()->SetMoreLogLabels();\r
226 hF->GetYaxis()->CenterTitle(1);\r
227 hF->Draw("p");\r
228 for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){\r
229 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;\r
230 g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");\r
231 //PutTrendValue(name[id], g->GetMean(2));\r
232 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));\r
233 }\r
234 leg->Draw(); gPad->SetLogx();\r
235 break;\r
236 case kTRDmom: // Energy loss\r
237 if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;\r
238 leg = new TLegend(.65, .7, .95, .99);\r
239 leg->SetHeader("Energy Loss");\r
240 leg->SetBorderSize(1); leg->SetFillColor(0);\r
241 title[0] = "Max & 90% quantile";\r
242 title[1] = "Mean & 60% quantile";\r
243 hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);\r
244 hF->SetMaximum(1.3);hF->SetMinimum(-.3);\r
245 hF->Draw("p");\r
246 for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){\r
247 if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;\r
248 ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");\r
249 //PutTrendValue(name[id], g->GetMean(2));\r
250 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));\r
251 }\r
252 leg->Draw();gPad->SetLogx(kFALSE);\r
253 break;\r
254 case kPtRes: // Pt resolution @ vertex\r
255 if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;\r
256 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); \r
257 pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();\r
258 pad->SetMargin(0.1, 0.022, 0.1, 0.023);\r
259 hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);\r
260 hF->SetMaximum(10.);hF->SetMinimum(-3.);\r
261 hF->GetXaxis()->SetMoreLogLabels();\r
262 hF->GetXaxis()->SetTitleOffset(1.2);\r
263 hF->GetYaxis()->CenterTitle();\r
264 hF->Draw("p");\r
265 //for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){\r
266 for(Int_t ig(2); ig<6; ig++){\r
267 if(!(g = (TGraphErrors*)arr->At(ig))) continue;\r
268 if(!g->GetN()) continue;\r
269 g->Draw("pl");\r
270 //PutTrendValue(name[id], g->GetMean(2));\r
271 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));\r
272 }\r
273 pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();\r
274 pad->SetMargin(0.1, 0.22, 0.1, 0.023);\r
275 hF = (TH1*)hF->Clone("hFcheckESD1");\r
276 hF->SetTitle("ITS+TPC");\r
277 hF->SetMaximum(10.);hF->SetMinimum(-3.);\r
278 hF->Draw("p");\r
279 leg = new TLegend(.78, .1, .99, .98);\r
280 leg->SetHeader("P_{t} @ DCA");\r
281 leg->SetBorderSize(1); leg->SetFillColor(0);\r
282 leg->SetTextAlign(22);\r
283 leg->SetTextFont(12);\r
284 leg->SetTextSize(0.03813559);\r
285 {\r
286 Int_t nPlots(0);\r
287 //for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){\r
288 for(Int_t ig(12); ig<16; ig++){\r
289 if(!(g = (TGraphErrors*)arr->At(ig))) continue;\r
290 if(!g->GetN()) continue;\r
291 nPlots++;\r
292 g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");\r
293 //PutTrendValue(name[id], g->GetMean(2));\r
294 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));\r
295 }\r
296 if(nPlots) leg->Draw();\r
297 }\r
298 break;\r
299 case 4: // plot a 3x3 canvas with tracking related histograms\r
300 gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);\r
301 gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);\r
302 gPad->Divide(3,3,0.,0.);\r
303 l=gPad->GetListOfPrimitives();\r
304 // eta-phi distr. for positive TPC tracks\r
305 pad = ((TVirtualPad*)l->At(0)); pad->cd();\r
306 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
307 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
308 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE); \r
309 h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos));\r
310 h2FtpcP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();\r
311 h2FtpcP->GetXaxis()->SetTitle("#eta");\r
312 h2FtpcP->GetXaxis()->CenterTitle();\r
313 h2FtpcP->GetXaxis()->SetTitleSize(0.07);\r
314 h2FtpcP->GetXaxis()->SetTitleOffset(0.8);\r
315 h2FtpcP->GetXaxis()->SetLabelSize(0.05);\r
316 h2FtpcP->GetYaxis()->SetTitle("detector #varphi");\r
317 h2FtpcP->GetYaxis()->CenterTitle();\r
318 h2FtpcP->GetYaxis()->SetTitleSize(0.07);\r
319 h2FtpcP->GetYaxis()->SetTitleOffset(0.8);\r
320 h2FtpcP->GetYaxis()->SetLabelSize(0.05);\r
321 h2FtpcP->SetTitle("");\r
322 h2FtpcP->Draw("colz");\r
323 lat->DrawLatex(-0.9, 3.6, "TPC positive ref. tracks");\r
324 //-----------------\r
325 // eta-phi distr. for negative TPC tracks\r
326 pad = ((TVirtualPad*)l->At(1)); pad->cd();\r
327 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
328 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
329 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
330 h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg));\r
331 h2FtpcN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();\r
332 h2FtpcN->GetXaxis()->SetTitle("#eta");\r
333 h2FtpcN->GetXaxis()->CenterTitle();\r
334 h2FtpcN->GetXaxis()->SetTitleSize(0.07);\r
335 h2FtpcN->GetXaxis()->SetTitleOffset(0.8);\r
336 h2FtpcN->GetXaxis()->SetLabelSize(0.05);\r
337 h2FtpcN->GetYaxis()->SetTitle("detector #varphi");\r
338 h2FtpcN->GetYaxis()->CenterTitle();\r
339 h2FtpcN->GetYaxis()->SetTitleSize(0.07);\r
340 h2FtpcN->GetYaxis()->SetTitleOffset(0.8);\r
341 h2FtpcN->GetYaxis()->SetLabelSize(0.05);\r
342 h2FtpcN->SetTitle("");\r
343 h2FtpcN->Draw("colz");\r
344 lat->DrawLatex(-0.9, 3.6, "TPC negative ref. tracks");\r
345 // eta-phi distr. for positive TRD tracks\r
346 pad = ((TVirtualPad*)l->At(3)); pad->cd();\r
347 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
348 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
349 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
350 h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos));\r
351 h2FtrdP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();\r
352 h2FtrdP->GetXaxis()->SetTitle("#eta");\r
353 h2FtrdP->GetXaxis()->CenterTitle();\r
354 h2FtrdP->GetXaxis()->SetTitleSize(0.07);\r
355 h2FtrdP->GetXaxis()->SetTitleOffset(0.8);\r
356 h2FtrdP->GetXaxis()->SetLabelSize(0.05);\r
357 h2FtrdP->GetYaxis()->SetTitle("detector #varphi");\r
358 h2FtrdP->GetYaxis()->CenterTitle();\r
359 h2FtrdP->GetYaxis()->SetTitleSize(0.07);\r
360 h2FtrdP->GetYaxis()->SetTitleOffset(0.8);\r
361 h2FtrdP->GetYaxis()->SetLabelSize(0.05);\r
362 h2FtrdP->SetMaximum(h2FtpcP->GetMaximum());\r
363 h2FtrdP->SetTitle("");\r
364 h2FtrdP->Draw("colz");\r
365 lat->DrawLatex(-0.9, 3.6, "TRD positive ref. tracks");\r
366 //-----------------\r
367 // eta-phi distr. for negative TRD tracks\r
368 pad = ((TVirtualPad*)l->At(4)); pad->cd();\r
369 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
370 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
371 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
372 h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg));\r
373 h2FtrdN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();\r
374 h2FtrdN->GetXaxis()->SetTitle("#eta");\r
375 h2FtrdN->GetXaxis()->CenterTitle();\r
376 h2FtrdN->GetXaxis()->SetTitleSize(0.07);\r
377 h2FtrdN->GetXaxis()->SetTitleOffset(0.8);\r
378 h2FtrdN->GetXaxis()->SetLabelSize(0.05);\r
379 h2FtrdN->GetYaxis()->SetTitle("detector #varphi");\r
380 h2FtrdN->GetYaxis()->CenterTitle();\r
381 h2FtrdN->GetYaxis()->SetTitleSize(0.07);\r
382 h2FtrdN->GetYaxis()->SetTitleOffset(0.8);\r
383 h2FtrdN->GetYaxis()->SetLabelSize(0.05);\r
384 h2FtrdN->SetMaximum(h2FtpcN->GetMaximum());\r
385 h2FtrdN->SetTitle("");\r
386 h2FtrdN->Draw("colz");\r
387 lat->DrawLatex(-0.9, 3.6, "TRD negative ref. tracks");\r
388 // eta-phi efficiency for positive TRD tracks\r
389 pad = ((TVirtualPad*)l->At(6)); pad->cd();\r
390 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
391 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
392 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
393 h2Feff = (TH2F*)h2FtrdP->Clone();\r
394 h2Feff->Reset();\r
395 h2Feff->Divide(h2FtrdP, h2FtpcP);\r
396 h2Feff->GetXaxis()->SetTitle("#eta");\r
397 h2Feff->GetXaxis()->CenterTitle();\r
398 h2Feff->GetXaxis()->SetTitleSize(0.07);\r
399 h2Feff->GetXaxis()->SetTitleOffset(0.8);\r
400 h2Feff->GetXaxis()->SetLabelSize(0.05);\r
401 h2Feff->GetYaxis()->SetTitle("detector #varphi");\r
402 h2Feff->GetYaxis()->CenterTitle();\r
403 h2Feff->GetYaxis()->SetTitleSize(0.07);\r
404 h2Feff->GetYaxis()->SetTitleOffset(0.8);\r
405 h2Feff->GetYaxis()->SetLabelSize(0.05);\r
406 h2Feff->SetMaximum(1.0);\r
407 h2Feff->SetTitle("");\r
408 h2Feff->Draw("colz");\r
409 lat->DrawLatex(-0.9, 3.6, "Efficiency positive tracks");\r
410 // eta-phi efficiency for negative TRD tracks\r
411 pad = ((TVirtualPad*)l->At(7)); pad->cd();\r
412 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
413 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
414 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
415 h2Feff = (TH2F*)h2FtrdN->Clone();\r
416 h2Feff->Reset();\r
417 h2Feff->Divide(h2FtrdN, h2FtpcN);\r
418 h2Feff->GetXaxis()->SetTitle("#eta");\r
419 h2Feff->GetXaxis()->CenterTitle();\r
420 h2Feff->GetXaxis()->SetTitleSize(0.07);\r
421 h2Feff->GetXaxis()->SetTitleOffset(0.8);\r
422 h2Feff->GetXaxis()->SetLabelSize(0.05);\r
423 h2Feff->GetYaxis()->SetTitle("detector #varphi");\r
424 h2Feff->GetYaxis()->CenterTitle();\r
425 h2Feff->GetYaxis()->SetTitleSize(0.07);\r
426 h2Feff->GetYaxis()->SetTitleOffset(0.8);\r
427 h2Feff->GetYaxis()->SetLabelSize(0.05);\r
428 h2Feff->SetMaximum(1.0);\r
429 h2Feff->SetTitle("");\r
430 h2Feff->Draw("colz");\r
431 lat->DrawLatex(-0.9, 3.6, "Efficiency negative tracks");\r
432 \r
433 // <ntracklets> vs (phi,eta)\r
434 pad = ((TVirtualPad*)l->At(2)); pad->cd();\r
435 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
436 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
437 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
438 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvNtrkl));\r
439 hProf2D->SetStats(kFALSE);\r
440 hProf2D->SetTitle("");\r
441 hProf2D->GetXaxis()->SetTitle("#eta");\r
442 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
443 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
444 hProf2D->GetXaxis()->CenterTitle();\r
445 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
446 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
447 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
448 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
449 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
450 hProf2D->GetYaxis()->CenterTitle();\r
451 hProf2D->SetMinimum(0.);\r
452 hProf2D->SetMaximum(6.);\r
453 hProf2D->Draw("colz");\r
454 lat->DrawLatex(-0.9, 3.6, "TRD <N_{tracklets}>");\r
455 // TPC-TRD matching efficiency vs pt\r
456 pad = ((TVirtualPad*)l->At(5)); pad->cd();\r
457 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);\r
458 pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);\r
459 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
460 hFeffP = TRDEfficiency(+1);\r
461 hFeffN = TRDEfficiency(-1);\r
462 h2F=new TH2F("rangeEffPt", "",10,0.,10.,10,0.,1.1);\r
463 h2F->SetStats(kFALSE);\r
464 h2F->GetXaxis()->SetTitle("p_{T} [GeV/c]");\r
465 h2F->GetXaxis()->SetTitleOffset(0.8); \r
466 h2F->GetXaxis()->SetTitleSize(0.07);\r
467 h2F->GetXaxis()->CenterTitle();\r
468 h2F->GetXaxis()->SetLabelSize(0.05);\r
469 h2F->GetYaxis()->SetTitle("TRD-TPC matching efficiency");\r
470 h2F->GetYaxis()->SetTitleOffset(0.8); \r
471 h2F->GetYaxis()->SetTitleSize(0.07);\r
472 h2F->GetYaxis()->SetLabelSize(0.05);\r
473 h2F->GetYaxis()->CenterTitle();\r
474 h2F->Draw();\r
475 line.SetLineStyle(2);\r
476 line.SetLineWidth(2);\r
477 line.DrawLine(h2F->GetXaxis()->GetXmin(), 0.7, h2F->GetXaxis()->GetXmax(), 0.7);\r
478 line.DrawLine(h2F->GetXaxis()->GetXmin(), 0.9, h2F->GetXaxis()->GetXmax(), 0.9);\r
479 hFeffP->SetMarkerStyle(20);\r
480 hFeffP->SetMarkerColor(2);\r
481 hFeffN->SetMarkerStyle(22);\r
482 hFeffN->SetMarkerColor(4);\r
483 hFeffP->Draw("same");\r
484 hFeffN->Draw("same");\r
485 leg=new TLegend(0.65, 0.2, 0.95, 0.4);\r
486 leg->SetFillColor(0);\r
487 leg->AddEntry(hFeffP, "positives", "p");\r
488 leg->AddEntry(hFeffN, "negatives", "p");\r
489 leg->Draw();\r
490 // create trending values for the TPC-TRD matching efficiency\r
491 // fit the efficiency histos with a constant in the range [1.0,1.5] GeV/c\r
492 fitFunc = new TF1("constantFunc","[0]",1.0,1.5);\r
493 hFeffP->Fit(fitFunc,"Q0","",1.0,1.5);\r
494 PutTrendValue("TrackingEffPos1GeV", fitFunc->GetParameter(0));\r
495 PutTrendValue("TrackingEffPos1GeVErr", fitFunc->GetParError(0));\r
496 hFeffN->Fit(fitFunc,"Q0","",1.0,1.5);\r
497 PutTrendValue("TrackingEffNeg1GeV", fitFunc->GetParameter(0));\r
498 PutTrendValue("TrackingEffNeg1GeVErr", fitFunc->GetParError(0));\r
499 \r
500 // Nclusters per TRD track\r
501 pad = ((TVirtualPad*)l->At(8)); pad->cd();\r
502 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.12);\r
503 pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);\r
504 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
505 pad->SetLogz();\r
506 h2F = dynamic_cast<TH2F*>(fHistos->At(kNClsTrackTRD));\r
507 h2F->SetStats(kFALSE);\r
508 h2F->SetTitle("");\r
509 h2F->GetXaxis()->SetTitle("p [GeV/c]");\r
510 h2F->GetXaxis()->SetTitleOffset(0.8); \r
511 h2F->GetXaxis()->SetTitleSize(0.07);\r
512 h2F->GetXaxis()->CenterTitle();\r
513 h2F->GetXaxis()->SetLabelSize(0.05);\r
514 h2F->GetYaxis()->SetTitle("#clusters per TRD track");\r
515 h2F->GetYaxis()->SetTitleOffset(0.8); \r
516 h2F->GetYaxis()->SetTitleSize(0.07);\r
517 h2F->GetYaxis()->CenterTitle();\r
518 h2F->GetYaxis()->SetLabelSize(0.05);\r
519 h2F->Draw("colz");\r
520 break;\r
521 case 5: // plot a 3x3 canvas with PID related histograms\r
522 gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);\r
523 gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);\r
524 gPad->Divide(3,3,0.,0.);\r
525 l=gPad->GetListOfPrimitives();\r
526 // eta-phi distr. for <Qtot> in layer 0\r
527 pad = ((TVirtualPad*)l->At(0)); pad->cd();\r
528 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
529 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
530 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
531 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+0));\r
532 hProf2D->SetStats(kFALSE);\r
533 hProf2D->SetTitle("");\r
534 hProf2D->GetXaxis()->SetTitle("#eta");\r
535 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
536 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
537 hProf2D->GetXaxis()->CenterTitle();\r
538 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
539 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
540 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
541 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
542 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
543 hProf2D->GetYaxis()->CenterTitle();\r
544 hProf2D->SetMinimum(0.);\r
545 hProf2D->SetMaximum(25.);\r
546 hProf2D->Draw("colz");\r
547 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 0");\r
548 // eta-phi distr. for <Qtot> in layer 1\r
549 pad = ((TVirtualPad*)l->At(3)); pad->cd();\r
550 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
551 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
552 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
553 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+1));\r
554 hProf2D->SetStats(kFALSE);\r
555 hProf2D->SetTitle("");\r
556 hProf2D->GetXaxis()->SetTitle("#eta");\r
557 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
558 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
559 hProf2D->GetXaxis()->CenterTitle();\r
560 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
561 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
562 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
563 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
564 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
565 hProf2D->GetYaxis()->CenterTitle();\r
566 hProf2D->SetMinimum(0.);\r
567 hProf2D->SetMaximum(25.);\r
568 hProf2D->Draw("colz");\r
569 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 1");\r
570 // eta-phi distr. for <Qtot> in layer 2\r
571 pad = ((TVirtualPad*)l->At(6)); pad->cd();\r
572 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
573 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
574 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
575 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+2));\r
576 hProf2D->SetStats(kFALSE);\r
577 hProf2D->SetTitle("");\r
578 hProf2D->GetXaxis()->SetTitle("#eta");\r
579 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
580 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
581 hProf2D->GetXaxis()->CenterTitle();\r
582 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
583 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
584 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
585 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
586 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
587 hProf2D->GetYaxis()->CenterTitle();\r
588 hProf2D->SetMinimum(0.);\r
589 hProf2D->SetMaximum(25.);\r
590 hProf2D->Draw("colz");\r
591 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 2");\r
592 // eta-phi distr. for <Qtot> in layer 3\r
593 pad = ((TVirtualPad*)l->At(1)); pad->cd();\r
594 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
595 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
596 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
597 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+3));\r
598 hProf2D->SetStats(kFALSE);\r
599 hProf2D->SetTitle("");\r
600 hProf2D->GetXaxis()->SetTitle("#eta");\r
601 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
602 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
603 hProf2D->GetXaxis()->CenterTitle();\r
604 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
605 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
606 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
607 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
608 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
609 hProf2D->GetYaxis()->CenterTitle();\r
610 hProf2D->SetMinimum(0.);\r
611 hProf2D->SetMaximum(25.);\r
612 hProf2D->Draw("colz");\r
613 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 3");\r
614 // eta-phi distr. for <Qtot> in layer 4\r
615 pad = ((TVirtualPad*)l->At(4)); pad->cd();\r
616 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
617 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
618 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
619 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+4));\r
620 hProf2D->SetStats(kFALSE);\r
621 hProf2D->SetTitle("");\r
622 hProf2D->GetXaxis()->SetTitle("#eta");\r
623 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
624 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
625 hProf2D->GetXaxis()->CenterTitle();\r
626 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
627 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
628 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
629 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
630 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
631 hProf2D->GetYaxis()->CenterTitle();\r
632 hProf2D->SetMinimum(0.);\r
633 hProf2D->SetMaximum(25.);\r
634 hProf2D->Draw("colz");\r
635 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 4");\r
636 // eta-phi distr. for <Qtot> in layer 5\r
637 pad = ((TVirtualPad*)l->At(7)); pad->cd();\r
638 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
639 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
640 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
641 hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+5));\r
642 hProf2D->SetStats(kFALSE);\r
643 hProf2D->SetTitle("");\r
644 hProf2D->GetXaxis()->SetTitle("#eta");\r
645 hProf2D->GetXaxis()->SetTitleOffset(0.8); \r
646 hProf2D->GetXaxis()->SetTitleSize(0.07);\r
647 hProf2D->GetXaxis()->CenterTitle();\r
648 hProf2D->GetXaxis()->SetLabelSize(0.05);\r
649 hProf2D->GetYaxis()->SetTitle("detector #varphi");\r
650 hProf2D->GetYaxis()->SetTitleOffset(0.8); \r
651 hProf2D->GetYaxis()->SetTitleSize(0.07);\r
652 hProf2D->GetYaxis()->SetLabelSize(0.05);\r
653 hProf2D->GetYaxis()->CenterTitle();\r
654 hProf2D->SetMinimum(0.);\r
655 hProf2D->SetMaximum(25.);\r
656 hProf2D->Draw("colz");\r
657 lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 5");\r
658 // PH versus slice number\r
659 pad = ((TVirtualPad*)l->At(2)); pad->cd();\r
660 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
661 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
662 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
663 h2F = dynamic_cast<TH2F*>(fHistos->At(kPHSlice));\r
664 h2F->SetStats(kFALSE);\r
665 h2F->SetTitle("");\r
666 h2F->GetXaxis()->SetTitle("slice");\r
667 h2F->GetXaxis()->SetTitleOffset(0.8); \r
668 h2F->GetXaxis()->SetTitleSize(0.07);\r
669 h2F->GetXaxis()->CenterTitle();\r
670 h2F->GetXaxis()->SetLabelSize(0.05);\r
671 h2F->GetYaxis()->SetTitle("PH");\r
672 h2F->GetYaxis()->SetTitleOffset(0.8); \r
673 h2F->GetYaxis()->SetTitleSize(0.07);\r
674 h2F->GetYaxis()->SetLabelSize(0.05);\r
675 h2F->GetYaxis()->CenterTitle();\r
676 h2F->Draw("colz");\r
677 //hProf = h2F->ProfileX("profileX");\r
678 //hProf->SetLineWidth(2);\r
679 //hProf->Draw("same");\r
680 // Qtot vs P\r
681 pad = ((TVirtualPad*)l->At(5)); pad->cd();\r
682 pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);\r
683 pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);\r
684 pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);\r
685 pad->SetLogz();\r
686 h2F = dynamic_cast<TH2F*>(fHistos->At(kQtotP));\r
687 h2F->SetStats(kFALSE);\r
688 h2F->SetTitle("");\r
689 h2F->GetXaxis()->SetTitle("P [GeV/c]");\r
690 h2F->GetXaxis()->SetTitleOffset(0.8); \r
691 h2F->GetXaxis()->SetTitleSize(0.07);\r
692 h2F->GetXaxis()->CenterTitle();\r
693 h2F->GetXaxis()->SetLabelSize(0.05);\r
694 h2F->GetYaxis()->SetRangeUser(0.0,100.0);\r
695 h2F->GetYaxis()->SetTitle("Q_{tot}");\r
696 h2F->GetYaxis()->SetTitleOffset(0.8); \r
697 h2F->GetYaxis()->SetTitleSize(0.07);\r
698 h2F->GetYaxis()->SetLabelSize(0.05);\r
699 h2F->GetYaxis()->CenterTitle();\r
700 h2F->Draw("colz");\r
701 // create trending value for the average Qtot at 1 GeV/c\r
702 hProf = h2F->ProfileX("profileQtot",1,h2F->GetYaxis()->FindBin(40.));\r
703 PutTrendValue("AvQtot1GeV", hProf->GetBinContent(hProf->GetXaxis()->FindBin(1.)));\r
704 PutTrendValue("AvQtot1GeVErr", hProf->GetBinError(hProf->GetXaxis()->FindBin(1.)));\r
705 break;\r
706 }\r
707 return kTRUE;\r
708}\r
709\r
710//____________________________________________________________________\r
711void AliTRDcheckESD::UserExec(Option_t *){\r
712 //\r
713 // Run the Analysis\r
714 //\r
715 fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
716 fMC = MCEvent();\r
717\r
718 if(!fESD){\r
719 AliError("ESD event missing.");\r
720 return;\r
721 }\r
722 \r
723 // Get MC information if available\r
724 AliStack * fStack = NULL;\r
725 if(HasMC()){\r
726 if(!fMC){ \r
727 AliWarning("MC event missing");\r
728 SetMC(kFALSE);\r
729 } else {\r
730 if(!(fStack = fMC->Stack())){\r
731 AliWarning("MC stack missing");\r
732 SetMC(kFALSE);\r
733 }\r
734 }\r
735 }\r
736 TH1 *h(NULL);\r
737 \r
738 // fill event vertex histos\r
739 h = (TH1F*)fHistos->At(kTPCVertex);\r
740 if(fESD->GetPrimaryVertexTPC()) h->Fill(fESD->GetPrimaryVertexTPC()->GetZv());\r
741 h = (TH1F*)fHistos->At(kEventVertex);\r
742 if(fESD->GetPrimaryVertex()) h->Fill(fESD->GetPrimaryVertex()->GetZv());\r
743 // fill the uncutted number of tracks\r
744 h = (TH1I*)fHistos->At(kNTracksAll);\r
745 h->Fill(fESD->GetNumberOfTracks());\r
746 \r
747 // counters for number of tracks in acceptance&DCA and for those with a minimum of TPC clusters\r
748 Int_t nTracksAcc=0;\r
749 Int_t nTracksTPC=0;\r
750 \r
751 AliESDtrack *esdTrack(NULL);\r
752 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){\r
753 esdTrack = fESD->GetTrack(itrk);\r
754\r
755 // track status\r
756 ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);\r
757\r
758 // track selection\r
759 Bool_t selected(kTRUE);\r
760 if(esdTrack->Pt() < fgkPt){ \r
761 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESD->GetEventNumberInFile(), itrk, esdTrack->Pt()));\r
762 selected = kFALSE;\r
763 }\r
764 if(TMath::Abs(esdTrack->Eta()) > fgkEta){\r
765 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));\r
766 selected = kFALSE;\r
767 }\r
768 if(!Bool_t(status & AliESDtrack::kTPCout)){\r
769 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] !TPCout", fESD->GetEventNumberInFile(), itrk));\r
770 selected = kFALSE;\r
771 }\r
772 if(esdTrack->GetKinkIndex(0) > 0){\r
773 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Kink", fESD->GetEventNumberInFile(), itrk));\r
774 selected = kFALSE;\r
775 }\r
776 \r
777 Float_t par[2], cov[3];\r
778 esdTrack->GetImpactParameters(par, cov);\r
779 if(selected && esdTrack->GetTPCNcls()>=10) {\r
780 // fill DCA histograms\r
781 h = (TH1F*)fHistos->At(kDCAxy); h->Fill(par[0]);\r
782 h = (TH1F*)fHistos->At(kDCAz); h->Fill(par[1]);\r
783 // fill pt distribution at this stage\r
784 h = (TH1F*)fHistos->At(kPt1); h->Fill(esdTrack->Pt());\r
785 }\r
786 if(IsCollision()){ // cuts on DCA\r
787 if(TMath::Abs(par[0]) > fgkTrkDCAxy){ \r
788 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));\r
789 selected = kFALSE;\r
790 }\r
791 if(TMath::Abs(par[1]) > fgkTrkDCAz){ \r
792 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));\r
793 selected = kFALSE;\r
794 }\r
795 }\r
796 Float_t theta=esdTrack->Theta();\r
797 Float_t phi=esdTrack->Phi();\r
798 Int_t nClustersTPC = esdTrack->GetTPCNcls();\r
799 Float_t eta=-TMath::Log(TMath::Tan(theta/2.));\r
800 if(selected) {\r
801 nTracksAcc++; // number of tracks in acceptance and DCA cut\r
802 // fill pt distribution at this stage\r
803 h = (TH1F*)fHistos->At(kPt2); h->Fill(esdTrack->Pt());\r
804 // TPC nclusters distribution\r
805 h = (TH1I*)fHistos->At(kNTPCCl); h->Fill(nClustersTPC);\r
806 if(esdTrack->Pt()>1.0) {\r
807 h = (TH1I*)fHistos->At(kNTPCCl2); h->Fill(nClustersTPC);\r
808 }\r
809 // (eta,nclustersTPC) distrib of TPC ref. tracks\r
810 h = (TH2F*)fHistos->At(kEtaNclsTPC); h->Fill(eta, nClustersTPC);\r
811 // (phi,nclustersTPC) distrib of TPC ref. tracks\r
812 h = (TH2F*)fHistos->At(kPhiNclsTPC); h->Fill(phi, nClustersTPC);\r
813 \r
814 }\r
815 \r
816 if(nClustersTPC < fgkNclTPC){ \r
817 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESD->GetEventNumberInFile(), itrk, nClustersTPC));\r
818 selected = kFALSE;\r
819 }\r
820 if(!selected) continue;\r
821 \r
822 // number of TPC reference tracks\r
823 nTracksTPC++;\r
824 \r
825 Int_t nTRD(esdTrack->GetNcls(2));\r
826 Double_t pt(esdTrack->Pt());\r
827 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);\r
828 // pid quality\r
829 Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);\r
830\r
831 TH3F *hhh;\r
832 // find position and momentum of the track at entrance in TRD\r
833 Double_t localCoord[3] = {0., 0., 0.};\r
834 Bool_t localCoordGood = esdTrack->GetXYZAt(298., fESD->GetMagneticField(), localCoord);\r
835 if(localCoordGood) {\r
836 hhh = (TH3F*)fHistos->At(kPropagXYvsP); hhh->Fill(localCoord[0], localCoord[1], esdTrack->GetP());\r
837 hhh = (TH3F*)fHistos->At(kPropagRZvsP); hhh->Fill(localCoord[2], TMath::Sqrt(localCoord[0]*localCoord[0]+localCoord[1]*localCoord[1]), esdTrack->GetP());\r
838 }\r
839 Double_t localMom[3] = {0., 0., 0.};\r
840 Bool_t localMomGood = esdTrack->GetPxPyPzAt(298., fESD->GetMagneticField(), localMom);\r
841 Double_t localPhi = (localMomGood ? TMath::ATan2(localMom[1], localMom[0]) : 0.0);\r
842 Double_t localSagitaPhi = (localCoordGood ? TMath::ATan2(localCoord[1], localCoord[0]) : 0.0);\r
843\r
844 // fill pt distribution at this stage\r
845 if(esdTrack->Charge()>0) {\r
846 h = (TH1F*)fHistos->At(kPt3pos); h->Fill(pt);\r
847 // fill eta-phi map of TPC positive ref. tracks\r
848 if(localCoordGood && localMomGood) {\r
849 hhh = (TH3F*)fHistos->At(kTPCRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);\r
850 }\r
851 }\r
852 if(esdTrack->Charge()<0) {\r
853 h = (TH1F*)fHistos->At(kPt3neg); h->Fill(pt);\r
854 // fill eta-phi map of TPC negative ref. tracks\r
855 if(localCoordGood && localMomGood) {\r
856 hhh = (TH3F*)fHistos->At(kTPCRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);\r
857 }\r
858 }\r
859 // TPC dE/dx vs P\r
860 h = (TH2F*)fHistos->At(kTPCDedx); h->Fill(esdTrack->GetP(), esdTrack->GetTPCsignal());\r
861 // (eta,phi) distrib of TPC ref. tracks\r
862 h = (TH2F*)fHistos->At(kEtaPhi); h->Fill(eta, phi);\r
863 \r
864 Int_t nTRDtrkl = esdTrack->GetTRDntracklets();\r
865 // TRD reference tracks\r
866 if(nTRDtrkl>=1) {\r
867 // fill pt distribution at this stage\r
868 if(esdTrack->Charge()>0) {\r
869 h = (TH1F*)fHistos->At(kPt4pos); h->Fill(pt);\r
870 // fill eta-phi map of TRD positive ref. tracks\r
871 if(localCoordGood && localMomGood) {\r
872 hhh = (TH3F*)fHistos->At(kTRDRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);\r
873 }\r
874 }\r
875 if(esdTrack->Charge()<0) {\r
876 h = (TH1F*)fHistos->At(kPt4neg); h->Fill(pt);\r
877 // fill eta-phi map of TRD negative ref. tracks\r
878 if(localCoordGood && localMomGood) {\r
879 hhh = (TH3F*)fHistos->At(kTRDRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);\r
880 }\r
881 }\r
882 TProfile2D *h2d;\r
883 // fill eta-phi map of TRD negative ref. tracks\r
884 if(localCoordGood && localMomGood) {\r
885 h2d = (TProfile2D*)fHistos->At(kTRDEtaPhiAvNtrkl); h2d->Fill(eta, localSagitaPhi, (Float_t)nTRDtrkl);\r
886 h2d = (TProfile2D*)fHistos->At(kTRDEtaDeltaPhiAvNtrkl); h2d->Fill(eta, localPhi-localSagitaPhi, (Float_t)nTRDtrkl);\r
887 }\r
888 // ntracklets/track vs P\r
889 h = (TH2F*)fHistos->At(kNTrackletsTRD); h->Fill(esdTrack->GetP(), nTRDtrkl);\r
890 // ntracklets/track vs P\r
891 h = (TH2F*)fHistos->At(kNClsTrackTRD); h->Fill(esdTrack->GetP(), esdTrack->GetTRDncls());\r
892 // (slicePH,sliceNo) distribution and Qtot from slices\r
893 for(Int_t iPlane=0; iPlane<6; iPlane++) {\r
894 Float_t Qtot=0;\r
895 for(Int_t iSlice=0; iSlice<8; iSlice++) {\r
896 if(esdTrack->GetTRDslice(iPlane, iSlice)>20.) {\r
897 h = (TH2F*)fHistos->At(kPHSlice); h->Fill(iSlice, esdTrack->GetTRDslice(iPlane, iSlice));\r
898 Qtot += esdTrack->GetTRDslice(iPlane, iSlice);\r
899 }\r
900 }\r
901 // Qtot>100 to avoid noise\r
902 if(Qtot>100.) {\r
903 h = (TH2F*)fHistos->At(kQtotP); h->Fill(esdTrack->GetTRDmomentum(iPlane), fgkQs*Qtot);\r
904 }\r
905 // Qtot>100 to avoid noise\r
906 // fgkQs*Qtot<40. so that the average will give a value close to the peak\r
907 if(localCoordGood && localMomGood && Qtot>100. && fgkQs*Qtot<40.) {\r
908 h2d = (TProfile2D*)fHistos->At(kTRDEtaPhiAvQtot+iPlane);\r
909 h2d->Fill(eta, localSagitaPhi, fgkQs*Qtot);\r
910 }\r
911 }\r
912 // theta distribution\r
913 h = (TH1F*)fHistos->At(kTheta); h->Fill(theta);\r
914 h = (TH1F*)fHistos->At(kPhi); h->Fill(phi);\r
915 } // end if nTRDtrkl>=1\r
916 \r
917 // look at external track param\r
918 const AliExternalTrackParam *op = esdTrack->GetOuterParam();\r
919 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();\r
920\r
921 Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.); \r
922 // read MC info if available\r
923 Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);\r
924 AliMCParticle *mcParticle(NULL);\r
925 if(HasMC()){\r
926 AliTrackReference *ref(NULL); \r
927 Int_t fLabel(esdTrack->GetLabel());\r
928 Int_t fIdx(TMath::Abs(fLabel));\r
929 if(fIdx > fStack->GetNtrack()) continue; \r
930 \r
931 // read MC particle \r
932 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {\r
933 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));\r
934 continue;\r
935 }\r
936 pt0 = mcParticle->Pt();\r
937 eta0 = mcParticle->Eta();\r
938 phi0 = mcParticle->Phi();\r
939 kPhysPrim = fMC->IsPhysicalPrimary(fIdx);\r
940\r
941 // read track references\r
942 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();\r
943 if(!nRefs){\r
944 AliWarning(Form("No TR found for track @ Label[%d].", fLabel));\r
945 continue;\r
946 }\r
947 Int_t iref = 0;\r
948 while(iref<nRefs){\r
949 ref = mcParticle->GetTrackReference(iref);\r
950 if(ref->LocalX() > fgkxTPC) break;\r
951 ref=NULL; iref++;\r
952 }\r
953 if(ref){ \r
954 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume\r
955 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));\r
956 }\r
957 } else { // track stopped in TPC \r
958 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));\r
959 }\r
960 ptTRD = ref->Pt();kFOUND=kTRUE;\r
961 } else { // use reconstructed values\r
962 if(op){\r
963 Double_t x(op->GetX());\r
964 if(x<fgkxTOF && x>fgkxTPC){\r
965 ptTRD=op->Pt();\r
966 kFOUND=kTRUE;\r
967 }\r
968 }\r
969\r
970 if(!kFOUND && ip){\r
971 ptTRD=ip->Pt();\r
972 kFOUND=kTRUE;\r
973 }\r
974 } // end if(HasMC())\r
975\r
976 if(kFOUND){\r
977 h = (TH2I*)fHistos->At(kTRDstat);\r
978 if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);\r
979 if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);\r
980 if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);\r
981 if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);\r
982 if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);\r
983 }\r
984 Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)\r
985 ,sgn(esdTrack->Charge()<0?0:1);\r
986 if(kBarrel && kPhysPrim) {\r
987 TH3 *h3 = (TH3S*)fHistos->At(kPtRes);\r
988 Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10; \r
989 h3->Fill(pt0, 1.e2*(pt/pt0-1.), \r
990 offset + 2*idx + sgn);\r
991 }\r
992 ((TH1*)fHistos->At(kNCl))->Fill(nTRD, 2*idx + sgn);\r
993 if(ip){\r
994 h = (TH2I*)fHistos->At(kTRDmom);\r
995 Float_t pTRD(0.);\r
996 for(Int_t ily=6; ily--;){\r
997 if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;\r
998 h->Fill(ip->GetP()-pTRD, ily);\r
999 }\r
1000 }\r
1001 } // end loop over tracks\r
1002 \r
1003 // fill the number of tracks histograms\r
1004 h = (TH1I*)fHistos->At(kNTracksAcc);\r
1005 h->Fill(nTracksAcc);\r
1006 h = (TH1I*)fHistos->At(kNTracksTPC);\r
1007 h->Fill(nTracksTPC);\r
695019eb 1008}\r
1009\r
1010//____________________________________________________________________\r
1011TObjArray* AliTRDcheckESD::Histos()\r
1012{\r
1013// Retrieve histograms array if already build or build it\r
1014\r
1015 if(fHistos) return fHistos;\r
1016\r
1017 fHistos = new TObjArray(kNhistos);\r
1018 //fHistos->SetOwner(kTRUE);\r
1019\r
1020 TH1 *h = NULL;\r
1021\r
1022 // clusters per track\r
1023 const Int_t kNpt(30);\r
1024 Float_t Pt(0.2);\r
1025 Float_t binsPt[kNpt+1];\r
1026 for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;\r
1027 if(!(h = (TH2S*)gROOT->FindObject("hNCl"))){\r
1028 h = new TH2S("hNCl", "Clusters per TRD track;N_{cl}^{TRD};SPECIES;entries", 60, 0., 180., 10, -0.5, 9.5);\r
1029 TAxis *ay(h->GetYaxis());\r
1030 ay->SetLabelOffset(0.015);\r
1031 for(Int_t i(0); i<AliPID::kSPECIES; i++){\r
1032 ay->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));\r
1033 ay->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));\r
1034 }\r
1035 } else h->Reset();\r
1036 fHistos->AddAt(h, kNCl); fNRefFigures++;\r
1037\r
1038 // status bits histogram\r
1039 const Int_t kNbits(5);\r
1040 Float_t Bits(.5);\r
1041 Float_t binsBits[kNbits+1];\r
1042 for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;\r
1043 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){\r
1044 h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);\r
1045 TAxis *ay(h->GetYaxis());\r
1046 ay->SetBinLabel(1, "kTPCout");\r
1047 ay->SetBinLabel(2, "kTRDin");\r
1048 ay->SetBinLabel(3, "kTRDout");\r
1049 ay->SetBinLabel(4, "kTRDpid");\r
1050 ay->SetBinLabel(5, "kTRDrefit");\r
1051 } else h->Reset();\r
1052 fHistos->AddAt(h, kTRDstat);\r
1053\r
1054 // energy loss\r
1055 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){\r
1056 h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);\r
1057 } else h->Reset();\r
1058 fHistos->AddAt(h, kTRDmom);\r
1059 //if(!HasMC()) return fHistos;\r
1060\r
1061 // pt resolution\r
1062 const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);\r
1063 Float_t DPt(-3.), Spec(-0.5);\r
1064 Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];\r
1065 for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;\r
1066 for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;\r
1067 if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){\r
1068 h = new TH3S("hPtRes", "P_{t} resolution @ DCA;p_{t}^{MC} [GeV/c];#Delta p_{t}/p_{t}^{MC} [%];SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspec, binsSpec);\r
1069 TAxis *az(h->GetZaxis());\r
1070 az->SetLabelOffset(0.015);\r
1071 for(Int_t i(0); i<AliPID::kSPECIES; i++){\r
1072 az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));\r
1073 az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));\r
1074 az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));\r
1075 az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));\r
1076 }\r
1077 } else h->Reset();\r
1078 fHistos->AddAt(h, kPtRes);\r
1079\r
1080 // TPC event vertex distribution\r
1081 if(!(h = (TH1F*)gROOT->FindObject("hTPCVertex"))){\r
1082 h = new TH1F("hTPCVertex", "Event vertex Z coord. from TPC tracks", 100, -25., 25.);\r
1083 } else h->Reset();\r
1084 fHistos->AddAt(h, kTPCVertex);\r
1085 \r
1086 // Event vertex\r
1087 if(!(h = (TH1F*)gROOT->FindObject("hEventVertex"))){\r
1088 h = new TH1F("hEventVertex", "Event vertex Z coord.", 100, -25., 25.);\r
1089 } else h->Reset();\r
1090 fHistos->AddAt(h, kEventVertex);\r
1091 \r
1092 // Number of all tracks\r
1093 if(!(h = (TH1I*)gROOT->FindObject("hNTracksAll"))){\r
1094 h = new TH1I("hNTracksAll", "Number of tracks per event, event vertex cuts", 5000, 0, 5000);\r
1095 } else h->Reset();\r
1096 fHistos->AddAt(h, kNTracksAll);\r
1097 \r
1098 // Number of tracks in acceptance and DCA cut\r
1099 if(!(h = (TH1I*)gROOT->FindObject("hNTracksAcc"))){\r
1100 h = new TH1I("hNTracksAcc", Form("Number of tracks per event, |#eta|<%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",\r
1101 fgkEta, fgkTrkDCAxy, fgkTrkDCAz), 5000, 0, 5000);\r
1102 } else h->Reset();\r
1103 fHistos->AddAt(h, kNTracksAcc);\r
1104 \r
1105 // Number of tracks in TPC (Ncls>10)\r
1106 if(!(h = (TH1I*)gROOT->FindObject("hNTracksTPC"))){\r
1107 h = new TH1I("hNTracksTPC", Form("Number of tracks per event, |#eta|<%.1f, pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1108 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 5000, 0, 5000);\r
1109 } else h->Reset();\r
1110 fHistos->AddAt(h, kNTracksTPC);\r
1111 \r
1112 // Distribution of DCA-xy\r
1113 if(!(h = (TH1F*)gROOT->FindObject("hDCAxy"))){\r
1114 h = new TH1F("hDCAxy", "Distribution of transverse DCA", 100, -100., 100.);\r
1115 } else h->Reset();\r
1116 fHistos->AddAt(h, kDCAxy);\r
1117 \r
1118 // Distribution of DCA-z\r
1119 if(!(h = (TH1F*)gROOT->FindObject("hDCAz"))){\r
1120 h = new TH1F("hDCAz", "Distribution of longitudinal DCA", 100, -100., 100.);\r
1121 } else h->Reset();\r
1122 fHistos->AddAt(h, kDCAz);\r
1123 \r
1124 Float_t binPtLimits[33] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,\r
1125 1.0, 1.1, 1.2, 1.3, 1.4, \r
1126 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,\r
1127 3.4, 3.8, 4.2, 4.6, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};\r
1128 // Pt distributions\r
1129 if(!(h = (TH1F*)gROOT->FindObject("hPt1"))){\r
1130 h = new TH1F("hPt1", Form("dN/dpt, |#eta|<%.1f and pt>%.1f", fgkEta, fgkPt), 32, binPtLimits);\r
1131 } else h->Reset();\r
1132 fHistos->AddAt(h, kPt1);\r
1133 \r
1134 if(!(h = (TH1F*)gROOT->FindObject("hPt2"))){\r
1135 h = new TH1F("hPt2", Form("dN/dpt, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",\r
1136 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 32, binPtLimits);\r
1137 } else h->Reset();\r
1138 fHistos->AddAt(h, kPt2);\r
1139 \r
1140 if(!(h = (TH1F*)gROOT->FindObject("hPt3pos"))){\r
1141 h = new TH1F("hPt3pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1142 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);\r
1143 } else h->Reset();\r
1144 fHistos->AddAt(h, kPt3pos);\r
1145 \r
1146 if(!(h = (TH1F*)gROOT->FindObject("hPt3neg"))){\r
1147 h = new TH1F("hPt3neg", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1148 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);\r
1149 } else h->Reset();\r
1150 fHistos->AddAt(h, kPt3neg);\r
1151 \r
1152 if(!(h = (TH1F*)gROOT->FindObject("hPt4pos"))){\r
1153 h = new TH1F("hPt4pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",\r
1154 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);\r
1155 } else h->Reset();\r
1156 fHistos->AddAt(h, kPt4pos);\r
1157 \r
1158 if(!(h = (TH1F*)gROOT->FindObject("hPt4neg"))){\r
1159 h = new TH1F("hPt4pos", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",\r
1160 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);\r
1161 } else h->Reset();\r
1162 fHistos->AddAt(h, kPt4neg);\r
1163 \r
1164 // theta distribution of TRD tracks\r
1165 if(!(h = (TH1F*)gROOT->FindObject("hTheta"))){\r
1166 h = new TH1F("hTheta", Form("dN/d#theta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",\r
1167 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 220,.5,2.7);\r
1168 } else h->Reset();\r
1169 fHistos->AddAt(h, kTheta);\r
1170 \r
1171 // phi distribution of TRD tracks\r
1172 if(!(h = (TH1F*)gROOT->FindObject("hPhi"))){\r
1173 h = new TH1F("hPhi", Form("dN/d#varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",\r
1174 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 157,0,6.28);\r
1175 } else h->Reset();\r
1176 fHistos->AddAt(h, kPhi);\r
1177 \r
1178 // TPC cluster distribution\r
1179 if(!(h = (TH1F*)gROOT->FindObject("hNTPCCl"))){\r
1180 h = new TH1I("hNTPCCl", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",\r
1181 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);\r
1182 } else h->Reset();\r
1183 fHistos->AddAt(h, kNTPCCl);\r
1184 \r
1185 if(!(h = (TH1I*)gROOT->FindObject("hNTPCCl2"))){\r
1186 h = new TH1F("hNTPCCl2", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, pt>1.0 GeV/c",\r
1187 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);\r
1188 } else h->Reset();\r
1189 fHistos->AddAt(h, kNTPCCl2);\r
1190 \r
1191 // dE/dx vs P for TPC reference tracks\r
1192 if(!(h = (TH2F*)gROOT->FindObject("hTPCDedx"))){\r
1193 h = new TH2F("hTPCDedx", Form("TPC dE/dx vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1194 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 120, 0,600.);\r
1195 } else h->Reset();\r
1196 fHistos->AddAt(h, kTPCDedx);\r
1197 \r
1198 // eta,phi distribution of TPC reference tracks\r
1199 if(!(h = (TH2F*)gROOT->FindObject("hEtaPhi"))){\r
1200 h = new TH2F("hEtaPhi", Form("TPC (#eta,#varphi), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1201 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 50, -1, 1, 157, 0, 6.28);\r
1202 } else h->Reset();\r
1203 fHistos->AddAt(h, kEtaPhi);\r
1204 \r
1205 // Nclusters vs eta distribution for TPC tracks\r
1206 if(!(h = (TH2F*)gROOT->FindObject("hEtaNclsTPC"))){\r
1207 h = new TH2F("hEtaNclsTPC", Form("TPC Nclusters vs. #eta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",\r
1208 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 50, -1, 1, 160, 0, 160.);\r
1209 } else h->Reset();\r
1210 fHistos->AddAt(h, kEtaNclsTPC);\r
1211 \r
1212 // Nclusters vs phi distribution for TPC reference tracks\r
1213 if(!(h = (TH2F*)gROOT->FindObject("hPhiNclsTPC"))){\r
1214 h = new TH2F("hPhiNclsTPC", Form("TPC Nclusters vs. #varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",\r
1215 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 157, 0, 6.28, 160, 0, 160.);\r
1216 } else h->Reset();\r
1217 fHistos->AddAt(h, kPhiNclsTPC);\r
1218 \r
1219 // Ntracklets/track vs P for TRD reference tracks\r
1220 Double_t binsP[19] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.7, 2.0,\r
1221 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 9.0, 12.0};\r
1222 if(!(h = (TH2F*)gROOT->FindObject("hNTrackletsTRD"))){\r
1223 h = new TH2F("hNTrackletsTRD", Form("TRD Ntracklets/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1224 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 7, -0.5, 6.5);\r
1225 } else h->Reset();\r
1226 fHistos->AddAt(h, kNTrackletsTRD);\r
1227 \r
1228 // Nclusters/track vs P for TRD reference tracks\r
1229 if(!(h = (TH2F*)gROOT->FindObject("hNClsTrackTRD"))){\r
1230 h = new TH2F("hNClsTrackTRD", Form("TRD Nclusters/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1231 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 180, 0., 180.);\r
1232 } else h->Reset();\r
1233 fHistos->AddAt(h, kNClsTrackTRD);\r
1234 \r
1235 // <PH> vs slice number for TRD reference tracklets\r
1236 if(!(h = (TH2F*)gROOT->FindObject("hPHSlice"))){\r
1237 h = new TH2F("hPHSlice", Form("<PH> vs sliceNo, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1238 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 2000, 0., 2000.);\r
1239 } else h->Reset();\r
1240 fHistos->AddAt(h, kPHSlice);\r
1241 \r
1242 // Qtot vs P for TRD reference tracklets\r
1243 if(!(h = (TH2F*)gROOT->FindObject("hQtotP"))){\r
1244 h = new TH2F("hQtotP", Form("Qtot(from slices) vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1245 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 400, 0., 200);\r
1246 } else h->Reset();\r
1247 fHistos->AddAt(h, kQtotP);\r
1248 \r
1249 // (X,Y,momentum) distribution after AliESDtrack::PropagateTo(r=300.)\r
1250 if(!(h = (TH3F*)gROOT->FindObject("hPropagXYvsP"))){\r
1251 h = new TH3F("hPropagXYvsP", Form("(x,y) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1252 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-500,500, 100,-500,500, 10, 0.,10.);\r
1253 } else h->Reset();\r
1254 fHistos->AddAt(h, kPropagXYvsP);\r
1255 \r
1256 // (R,Z,momentum) distribution after AliESDtrack::PropagateTo(r=300.)\r
1257 if(!(h = (TH3F*)gROOT->FindObject("hPropagRZvsP"))){\r
1258 h = new TH3F("hPropagRZvsP", Form("(r,z) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1259 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-350., 350., 100,0.,500., 10, 0.,10.);\r
1260 } else h->Reset();\r
1261 fHistos->AddAt(h, kPropagRZvsP);\r
1262 \r
1263 Float_t etaBinLimits[101]; \r
1264 for(Int_t i=0; i<101; i++) etaBinLimits[i] = -1.0 + i*2.0/100.;\r
1265 Float_t phiBinLimits[151];\r
1266 for(Int_t i=0; i<151; i++) phiBinLimits[i] = -1.1*TMath::Pi() + i*2.2*TMath::Pi()/150.;\r
1267 // (eta,detector phi,P) distribution of reference TPC positive tracks\r
1268 if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksPos"))){\r
1269 h = new TH3F("hTPCRefTracksPos", Form("(#eta,detector #varphi,p) for TPC positive reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1270 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);\r
1271 } else h->Reset();\r
1272 fHistos->AddAt(h, kTPCRefTracksPos);\r
1273 \r
1274 // (eta,detector phi,P) distribution of reference TPC negative tracks\r
1275 if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksNeg"))){\r
1276 h = new TH3F("hTPCRefTracksNeg", Form("(#eta,detector #varphi,p) for TPC negative reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",\r
1277 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);\r
1278 } else h->Reset();\r
1279 fHistos->AddAt(h, kTPCRefTracksNeg);\r
1280 \r
1281 // (eta,detector phi,P) distribution of reference TRD positive tracks\r
1282 if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksPos"))){\r
1283 h = new TH3F("hTRDRefTracksPos", Form("(#eta,detector #varphi,p) for TRD positive reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",\r
1284 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);\r
1285 } else h->Reset();\r
1286 fHistos->AddAt(h, kTRDRefTracksPos);\r
1287 \r
1288 // (eta,detector phi,P) distribution of reference TRD negative tracks\r
1289 if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksNeg"))){\r
1290 h = new TH3F("hTRDRefTracksNeg", Form("(#eta,detector #varphi,p) for TRD negative reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",\r
1291 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);\r
1292 } else h->Reset();\r
1293 fHistos->AddAt(h, kTRDRefTracksNeg);\r
1294 \r
1295 // (eta,detector phi) profile of average number of TRD tracklets/track\r
1296 if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaPhiAvNtrkl"))){\r
1297 h = new TProfile2D("hTRDEtaPhiAvNtrkl", Form("<Ntracklets/track> vs (#eta,detector #varphi) for TRD reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",\r
1298 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 150, -1.1*TMath::Pi(), 1.1*TMath::Pi());\r
1299 } else h->Reset();\r
1300 fHistos->AddAt(h, kTRDEtaPhiAvNtrkl);\r
1301\r
1302 // (eta,delta phi) profile of average number of TRD tracklets/track\r
1303 if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaDeltaPhiAvNtrkl"))){\r
1304 h = new TProfile2D("hTRDEtaDeltaPhiAvNtrkl", Form("<Ntracklets/track> vs (#eta, #Delta#varphi) for TRD reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",\r
1305 fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -0.4*TMath::Pi(), 0.4*TMath::Pi());\r
1306 } else h->Reset();\r
1307 fHistos->AddAt(h, kTRDEtaDeltaPhiAvNtrkl);\r
1308 \r
1309 // (eta, detector phi) profile of average tracklet Qtot from slices\r
1310 for(Int_t iLayer=0;iLayer<6;iLayer++) {\r
1311 if(!(h = (TProfile2D*)gROOT->FindObject(Form("hTRDEtaPhiAvQtot_Layer%d",iLayer)))) {\r
1312 h = new TProfile2D(Form("hTRDEtaPhiAvQtot_Layer%d",iLayer),\r
1313 Form("<Q_{tot}> vs (#eta, detector #varphi) for TRD reference tracks (layer %d), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",\r
1314 iLayer, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -1.1*TMath::Pi(), 1.1*TMath::Pi());\r
1315 } else h->Reset();\r
1316 fHistos->AddAt(h, kTRDEtaPhiAvQtot+iLayer);\r
1317 }\r
1318 \r
1319 return fHistos;\r
1320}\r
1321\r
1322//____________________________________________________________________\r
1323Bool_t AliTRDcheckESD::Load(const Char_t *file, const Char_t *dir, const Char_t *name)\r
1324{\r
1325// Load data from performance file\r
1326\r
1327 if(!TFile::Open(file)){\r
1328 AliWarning(Form("Couldn't open file %s.", file));\r
1329 return kFALSE;\r
1330 }\r
1331 if(dir){\r
1332 if(!gFile->cd(dir)){\r
1333 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
1334 return kFALSE;\r
1335 }\r
1336 }\r
1337 TObjArray *o(NULL);\r
1338 const Char_t *tn=(name ? name : GetName());\r
1339 if(!(o = (TObjArray*)gDirectory->Get(tn))){\r
1340 AliWarning(Form("Missing histogram container %s.", tn));\r
1341 return kFALSE;\r
1342 }\r
1343 fHistos = (TObjArray*)o->Clone(GetName());\r
1344 gFile->Close();\r
1345 return kTRUE;\r
1346}\r
1347\r
1348//_______________________________________________________\r
1349Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)\r
1350{\r
1351// Dump trending value to default file\r
1352\r
1353 if(!fgFile){\r
1354 fgFile = fopen("TRD.Performance.txt", "at");\r
1355 }\r
1356 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);\r
1357 return kTRUE;\r
1358}\r
1359\r
1360//____________________________________________________________________\r
1361void AliTRDcheckESD::Terminate(Option_t *)\r
1362{\r
1363// Steer post-processing \r
1364 if(!fHistos){\r
1365 fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));\r
1366 if(!fHistos){\r
1367 AliError("Histogram container not found in output");\r
1368 return;\r
1369 }\r
1370 }\r
1371\r
1372 const Char_t *name[kNrefs] = {\r
1373 "Ncl", "Eff", "Eloss", "PtResDCA"\r
1374 };\r
1375 TObjArray *arr(NULL); TGraph *g(NULL);\r
1376 if(!fResults){\r
1377 fResults = new TObjArray(kNrefs);\r
1378 fResults->SetOwner();\r
1379 fResults->SetName("results");\r
1380 for(Int_t iref(0); iref<kNrefs; iref++){\r
1381 fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);\r
1382 arr->SetName(name[iref]); arr->SetOwner();\r
1383 switch(iref){\r
1384 case kNCl:\r
1385 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){\r
1386 arr->AddAt(g = new TGraphErrors(), ig);\r
1387 g->SetLineColor(ig+1); \r
1388 g->SetMarkerColor(ig+1); \r
1389 g->SetMarkerStyle(ig+20); \r
1390 g->SetName(Form("s%d", ig));\r
1391 switch(ig){\r
1392 case 0: g->SetTitle("ALL"); break;\r
1393 case 1: g->SetTitle("NEG"); break;\r
1394 case 2: g->SetTitle("POS"); break;\r
1395 default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;\r
1396 };\r
1397 }\r
1398 break;\r
1399 case kTRDmom:\r
1400 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){\r
1401 arr->AddAt(g = new TGraphAsymmErrors(), ig);\r
1402 g->SetLineColor(ig+1); \r
1403 g->SetMarkerColor(ig+1); \r
1404 g->SetMarkerStyle(ig+20); \r
1405 }\r
1406 break;\r
1407 case kPtRes:\r
1408 for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){\r
1409 Int_t ig(2*idx);\r
1410 arr->AddAt(g = new TGraphErrors(), ig);\r
1411 g->SetLineColor(kRed-idx); \r
1412 g->SetMarkerColor(kRed-idx); \r
1413 g->SetMarkerStyle(20+idx); \r
1414 g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));\r
1415 arr->AddAt(g = new TGraphErrors(), ig+1);\r
1416 g->SetLineColor(kBlue-idx); \r
1417 g->SetMarkerColor(kBlue-idx); \r
1418 g->SetMarkerStyle(20+idx); \r
1419 g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));\r
1420\r
1421 ig+=10;\r
1422 arr->AddAt(g = new TGraphErrors(), ig);\r
1423 g->SetLineColor(kRed-idx); \r
1424 g->SetMarkerColor(kRed-idx); \r
1425 g->SetMarkerStyle(20+idx); \r
1426 g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));\r
1427 arr->AddAt(g = new TGraphErrors(), ig+1);\r
1428 g->SetLineColor(kBlue-idx); \r
1429 g->SetMarkerColor(kBlue-idx); \r
1430 g->SetMarkerStyle(20+idx); \r
1431 g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));\r
1432 }\r
1433 break;\r
1434 default:\r
1435 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){\r
1436 arr->AddAt(g = new TGraphErrors(), ig);\r
1437 g->SetLineColor(ig+1); \r
1438 g->SetMarkerColor(ig+1); \r
1439 g->SetMarkerStyle(ig+20); \r
1440 }\r
1441 break;\r
1442 }\r
1443 }\r
1444 }\r
1445 TH1 *h1[2] = {NULL, NULL};\r
1446 TH2I *h2(NULL);\r
1447 TAxis *ax(NULL);\r
1448\r
1449 // No of clusters\r
1450 if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;\r
1451 ax = h2->GetXaxis();\r
1452 arr = (TObjArray*)fResults->At(kNCl);\r
1453 // All tracks\r
1454 h1[0] = h2->ProjectionX("Ncl_px");\r
1455 TGraphErrors *ge=(TGraphErrors*)arr->At(0);\r
1456 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){\r
1457 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));\r
1458 }\r
1459 // All charged tracks\r
1460 TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};\r
1461 hNclCh[0]->Reset();hNclCh[1]->Reset();\r
1462 for(Int_t is(1); is<=AliPID::kSPECIES; is++){\r
1463 hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg\r
1464 hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is)); // pos\r
1465 }\r
1466 if(Int_t(hNclCh[0]->GetEntries())){\r
1467 ge=(TGraphErrors*)arr->At(1);\r
1468 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){\r
1469 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));\r
1470 }\r
1471 }\r
1472 if(Int_t(hNclCh[1]->GetEntries())){\r
1473 ge=(TGraphErrors*)arr->At(2);\r
1474 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){\r
1475 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));\r
1476 }\r
1477 }\r
1478 // Species wise\r
1479 for(Int_t is(1); is<=AliPID::kSPECIES; is++){\r
1480 h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);\r
1481 if(!Int_t(h1[0]->GetEntries())) continue;\r
1482 ge=(TGraphErrors*)arr->At(2+is);\r
1483 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){\r
1484 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));\r
1485 }\r
1486 }\r
1487 fNRefFigures = 1;\r
1488\r
1489 // EFFICIENCY\r
1490 // geometrical efficiency\r
1491 if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;\r
1492 arr = (TObjArray*)fResults->At(kTRDstat);\r
1493 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);\r
1494 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);\r
1495 Process(h1, (TGraphErrors*)arr->At(0));\r
1496 delete h1[0];delete h1[1];\r
1497 // tracking efficiency\r
1498 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);\r
1499 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);\r
1500 Process(h1, (TGraphErrors*)arr->At(1));\r
1501 delete h1[1];\r
1502 // PID efficiency\r
1503 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);\r
1504 Process(h1, (TGraphErrors*)arr->At(2));\r
1505 delete h1[1];\r
1506 // Refit efficiency\r
1507 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);\r
1508 Process(h1, (TGraphErrors*)arr->At(3));\r
1509 delete h1[1];\r
1510 fNRefFigures++;\r
1511\r
1512 // ENERGY LOSS\r
1513 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;\r
1514 arr = (TObjArray*)fResults->At(kTRDmom);\r
1515 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);\r
1516 ax=h2->GetXaxis();\r
1517 const Int_t nq(4);\r
1518 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};\r
1519 Double_t yq[nq];\r
1520 for(Int_t ily=6; ily--;){\r
1521 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);\r
1522 h1[0]->GetQuantiles(nq,yq,xq);\r
1523 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));\r
1524 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);\r
1525 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());\r
1526 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);\r
1527\r
1528 //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);\r
1529 delete h1[0];\r
1530 }\r
1531 fNRefFigures++;\r
1532// if(!HasMC()) return;\r
1533\r
1534 // Pt RESOLUTION @ DCA\r
1535 TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};\r
1536 if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;\r
1537 arr = (TObjArray*)fResults->At(kPtRes);\r
1538 TAxis *az(h3->GetZaxis());\r
1539 for(Int_t i(0); i<AliPID::kSPECIES; i++){\r
1540 Int_t idx(2*i);\r
1541 az->SetRange(idx+1, idx+2); \r
1542 gg[1] = (TGraphErrors*)arr->At(idx);\r
1543 gg[0] = (TGraphErrors*)arr->At(idx+1);\r
1544 Process2D((TH2*)h3->Project3D("yx"), gg);\r
1545\r
1546 idx+=10;\r
1547 az->SetRange(idx+1, idx+2); \r
1548 gg[1] = (TGraphErrors*)arr->At(idx);\r
1549 gg[0] = (TGraphErrors*)arr->At(idx+1);\r
1550 Process2D((TH2*)h3->Project3D("yx"), gg);\r
1551 }\r
1552 fNRefFigures++;\r
1553 \r
1554 // 3x3 tracking summary canvas\r
1555 fNRefFigures++;\r
1556 // 3x3 PID summary canvas\r
1557 fNRefFigures++;\r
1558}\r
1559\r
1560//____________________________________________________________________\r
1561Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)\r
1562{\r
1563 switch(pdg){\r
1564 case kElectron: \r
1565 case kPositron: return AliPID::kElectron; \r
1566 case kMuonPlus:\r
1567 case kMuonMinus: return AliPID::kMuon; \r
1568 case kPiPlus: \r
1569 case kPiMinus: return AliPID::kPion; \r
1570 case kKPlus: \r
1571 case kKMinus: return AliPID::kKaon;\r
1572 case kProton: \r
1573 case kProtonBar: return AliPID::kProton;\r
1574 } \r
1575 return -1;\r
1576}\r
1577\r
1578//____________________________________________________________________\r
1579void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)\r
1580{\r
1581// Generic function to process one reference plot\r
1582\r
1583 Int_t n1 = 0, n2 = 0, ip=0;\r
1584 Double_t eff = 0.;\r
1585\r
1586 TAxis *ax = h1[0]->GetXaxis();\r
1587 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){\r
1588 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;\r
1589 n2 = (Int_t)h1[1]->GetBinContent(ib);\r
1590 eff = n2/Float_t(n1);\r
1591\r
1592 ip=g->GetN();\r
1593 g->SetPoint(ip, ax->GetBinCenter(ib), eff);\r
1594 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);\r
1595 }\r
1596} \r
1597//________________________________________________________\r
1598void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)\r
1599{\r
1600 //\r
1601 // Do the processing\r
1602 //\r
1603\r
1604 Int_t n = 0;\r
1605 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);\r
1606 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);\r
1607 TF1 f("fg", "gaus", -3.,3.);\r
1608 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){\r
1609 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);\r
1610 TH1D *h = h2->ProjectionY("py", ibin, ibin);\r
1611 if(h->GetEntries()<100) continue;\r
1612 //AdjustF1(h, f);\r
1613\r
1614 h->Fit(&f, "QN");\r
1615 Int_t ip = g[0]->GetN();\r
1616 g[0]->SetPoint(ip, x, f.GetParameter(1));\r
1617 g[0]->SetPointError(ip, 0., f.GetParError(1));\r
1618 g[1]->SetPoint(ip, x, f.GetParameter(2));\r
1619 g[1]->SetPointError(ip, 0., f.GetParError(2));\r
1620 }\r
1621 return;\r
1622}\r
1623//____________________________________________________________________\r
1624void AliTRDcheckESD::PrintStatus(ULong_t status)\r
1625{\r
1626// Dump track status to stdout\r
1627\r
1628 printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"\r
1629 ,Bool_t(status & AliESDtrack::kITSin)\r
1630 ,Bool_t(status & AliESDtrack::kITSout)\r
1631 ,Bool_t(status & AliESDtrack::kITSrefit)\r
1632 ,Bool_t(status & AliESDtrack::kTPCin)\r
1633 ,Bool_t(status & AliESDtrack::kTPCout)\r
1634 ,Bool_t(status & AliESDtrack::kTPCrefit)\r
1635 ,Bool_t(status & AliESDtrack::kTPCpid)\r
1636 ,Bool_t(status & AliESDtrack::kTRDin)\r
1637 ,Bool_t(status & AliESDtrack::kTRDout)\r
1638 ,Bool_t(status & AliESDtrack::kTRDrefit)\r
1639 ,Bool_t(status & AliESDtrack::kTRDpid)\r
1640 ,Bool_t(status & AliESDtrack::kTRDStop)\r
1641 ,Bool_t(status & AliESDtrack::kHMPIDout)\r
1642 ,Bool_t(status & AliESDtrack::kHMPIDpid)\r
1643 );\r
1644}\r
1645//____________________________________________________________________\r
1646TH2F* AliTRDcheckESD::Proj3D(TH3F* hist, TH2F* accMap, Int_t zbinLow, Int_t zbinHigh, Float_t &entries) {\r
1647 //\r
1648 // Project a 3D histogram to a 2D histogram in the Z axis interval [zbinLow,zbinHigh] \r
1649 // Return the 2D histogram and also the number of entries into this projection (entries)\r
1650\r
1651 Int_t nBinsX = hist->GetXaxis()->GetNbins(); // X and Y axis bins are assumed to be all equal\r
1652 Float_t minX = hist->GetXaxis()->GetXmin();\r
1653 Float_t maxX = hist->GetXaxis()->GetXmax();\r
1654 Int_t nBinsY = hist->GetYaxis()->GetNbins();\r
1655 Float_t minY = hist->GetYaxis()->GetXmin();\r
1656 Float_t maxY = hist->GetYaxis()->GetXmax();\r
1657 Int_t nBinsZ = hist->GetZaxis()->GetNbins(); // Z axis bins (pt) might have different widths\r
1658 //Float_t minZ = hist->GetZaxis()->GetXmin();\r
1659 //Float_t maxZ = hist->GetZaxis()->GetXmax();\r
1660\r
1661 TH2F* projHisto = (TH2F*)gROOT->FindObject("projHisto");\r
1662 if(projHisto) \r
1663 projHisto->Reset();\r
1664 else\r
1665 projHisto = new TH2F("projHisto", "projection", nBinsX, minX, maxX, nBinsY, minY, maxY);\r
1666\r
1667 const Double_t kMinAccFraction = 0.1;\r
1668 entries = 0.0;\r
1669 Double_t maxAcc = (accMap ? accMap->GetMaximum() : -1);\r
1670 \r
1671 for(Int_t iZ=1; iZ<=nBinsZ; iZ++) {\r
1672 if(iZ<zbinLow) continue;\r
1673 if(iZ>zbinHigh) continue;\r
1674 for(Int_t iX=1; iX<=nBinsX; iX++) {\r
1675 for(Int_t iY=1; iY<=nBinsY; iY++) {\r
1676 if(accMap && maxAcc>0) {\r
1677 if(accMap->GetBinContent(iX,iY)/maxAcc>kMinAccFraction)\r
1678 projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));\r
1679 }\r
1680 else // no acc. cut \r
1681 projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));\r
1682 // count only the entries which are inside the acceptance map\r
1683 if(accMap && maxAcc>0) { \r
1684 if(accMap->GetBinContent(iX,iY)/maxAcc>kMinAccFraction)\r
1685 entries+=hist->GetBinContent(iX,iY,iZ);\r
1686 }\r
1687 else // no acc. cut\r
1688 entries+=hist->GetBinContent(iX,iY,iZ);\r
1689 }\r
1690 }\r
1691 }\r
1692 return projHisto;\r
1693}\r
1694//____________________________________________________________________\r
1695TH1F* AliTRDcheckESD::TRDEfficiency(Short_t positives) {\r
1696 //\r
1697 // Calculate the TRD-TPC matching efficiency as function of pt\r
1698 //\r
1699 TH3F* tpc3D(NULL); TH3F* trd3D(NULL);\r
1700 if(positives>0) { // get the histos for positive tracks \r
1701 if(!(tpc3D = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos)))) return 0x0;\r
1702 if(!(trd3D = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos)))) return 0x0;\r
1703 }\r
1704 else { // get the histos for positive tracks\r
1705 if(!(tpc3D = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg)))) return 0x0;\r
1706 if(!(trd3D = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg)))) return 0x0;\r
1707 }\r
1708 \r
1709 Int_t nBinsZ = trd3D->GetZaxis()->GetNbins();\r
1710 // project everything on the eta-phi map to obtain an acceptance map (make sure there is enough statistics)\r
1711 Float_t nada = 0.;\r
1712 TH2F *trdAcc = (TH2F*)Proj3D(trd3D, 0x0, 1, nBinsZ, nada)->Clone();\r
1713 // get the bin limits from the Z axis of 3D histos\r
1714 Float_t *ptBinLimits = new Float_t[nBinsZ+1];\r
1715 for(Int_t i=1; i<=nBinsZ; i++) {\r
1716 ptBinLimits[i-1] = trd3D->GetZaxis()->GetBinLowEdge(i);\r
1717 }\r
1718 ptBinLimits[nBinsZ] = trd3D->GetZaxis()->GetBinUpEdge(nBinsZ);\r
1719 TH1F *efficiency = new TH1F("eff", "TRD-TPC matching efficiency", nBinsZ, ptBinLimits);\r
1720 // loop over Z bins\r
1721 for(Int_t i=1; i<=nBinsZ; i++) {\r
1722 Float_t tpcEntries = 0.0; Float_t trdEntries = 0.0;\r
1723 Proj3D(tpc3D, trdAcc, i, i, tpcEntries);\r
1724 Proj3D(trd3D, trdAcc, i, i, trdEntries);\r
1725 Float_t ratio = 0;\r
1726 if(tpcEntries>0) ratio = trdEntries/tpcEntries;\r
1727 Float_t error = 0;\r
1728 if(tpcEntries>0 && trdEntries>0) \r
1729 error = ratio*TMath::Sqrt((1.0/tpcEntries)+(1.0/trdEntries));\r
1730 if(ratio>0.0) {\r
1731 efficiency->SetBinContent(i,ratio);\r
1732 efficiency->SetBinError(i,error);\r
1733 }\r
1734 } // end loop over Z bins\r
1735 \r
1736 efficiency->SetLineColor((positives>0 ? 2 : 4));\r
1737 return efficiency;\r
1738}\r