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