]> 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
update task to read results from AnalysisResults.root file
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckESD.cxx
CommitLineData
1ee39b3a 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 <TObjArray.h>
629ae9b8 31#include <TPad.h>
32#include <TLegend.h>
33#include <TF1.h>
1ee39b3a 34#include <TH2I.h>
fb09668f 35#include <TH3S.h>
1ee39b3a 36#include <TGraphErrors.h>
37#include <TGraphAsymmErrors.h>
38#include <TFile.h>
39#include <TTree.h>
40#include <TROOT.h>
41#include <TChain.h>
42#include <TParticle.h>
43
44#include "AliLog.h"
45#include "AliAnalysisManager.h"
46#include "AliESDEvent.h"
47#include "AliESDkink.h"
48#include "AliMCEvent.h"
49#include "AliESDInputHandler.h"
50#include "AliMCEventHandler.h"
51
52#include "AliESDtrack.h"
53#include "AliMCParticle.h"
54#include "AliPID.h"
55#include "AliStack.h"
56#include "AliTrackReference.h"
57
58#include "AliTRDcheckESD.h"
59
60ClassImp(AliTRDcheckESD)
61
62const Float_t AliTRDcheckESD::fgkxTPC = 290.;
63const Float_t AliTRDcheckESD::fgkxTOF = 365.;
1f6edfbc 64const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
dfd7d48b 658, 4, 2, 20};
9b496b5e 66FILE* AliTRDcheckESD::fgFile = NULL;
1ee39b3a 67
36d21092 68const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;
69const Int_t AliTRDcheckESD::fgkEvVertexN = 1;
70const Float_t AliTRDcheckESD::fgkTrkDCAxy = 40.;
71const Float_t AliTRDcheckESD::fgkTrkDCAz = 15.;
81979445 72const Int_t AliTRDcheckESD::fgkNclTPC = 100;
36d21092 73const Float_t AliTRDcheckESD::fgkPt = 0.2;
74const Float_t AliTRDcheckESD::fgkEta = 0.9;
75
1ee39b3a 76//____________________________________________________________________
77AliTRDcheckESD::AliTRDcheckESD():
f8f46e4d 78 AliAnalysisTaskSE()
1ee39b3a 79 ,fStatus(0)
629ae9b8 80 ,fNRefFigures(0)
9b496b5e 81 ,fESD(NULL)
82 ,fMC(NULL)
83 ,fHistos(NULL)
84 ,fResults(NULL)
1ee39b3a 85{
86 //
87 // Default constructor
88 //
db99a57a 89 SetNameTitle("checkESD", "Check TRD @ ESD level");
1ee39b3a 90 SetMC(kTRUE);
f8f46e4d 91}
92
f2e89a4c 93//____________________________________________________________________
f8f46e4d 94AliTRDcheckESD::AliTRDcheckESD(char* name):
95 AliAnalysisTaskSE(name)
96 ,fStatus(0)
629ae9b8 97 ,fNRefFigures(0)
f8f46e4d 98 ,fESD(NULL)
99 ,fMC(NULL)
100 ,fHistos(NULL)
101 ,fResults(NULL)
102{
103 //
104 // Default constructor
105 //
106 SetMC(kTRUE);
f2e89a4c 107 SetTitle("Check TRD @ ESD level");
f8f46e4d 108 DefineOutput(1, TObjArray::Class());
1ee39b3a 109}
110
111//____________________________________________________________________
112AliTRDcheckESD::~AliTRDcheckESD()
113{
114// Destructor
115 if(fHistos){
116 //fHistos->Delete();
117 delete fHistos;
118 }
119 if(fResults){
120 fResults->Delete();
121 delete fResults;
122 }
123}
124
125//____________________________________________________________________
f8f46e4d 126void AliTRDcheckESD::UserCreateOutputObjects()
1ee39b3a 127{
128 //
129 // Create Output Containers (TObjectArray containing 1D histograms)
130 //
1ee39b3a 131 Histos();
132}
133
629ae9b8 134
1ee39b3a 135//____________________________________________________________________
629ae9b8 136Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
1ee39b3a 137{
629ae9b8 138 if(ifig>=fNRefFigures){
139 AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
140 return kFALSE;
1ee39b3a 141 }
629ae9b8 142 if(!gPad){
143 AliWarning("Please provide a canvas to draw results.");
144 return kFALSE;
94f7dff7 145 } else {
146 gPad->SetLogx(0);gPad->SetLogy(0);
147 gPad->SetMargin(0.125, 0.015, 0.1, 0.015);
1ee39b3a 148 }
149
629ae9b8 150 const Char_t *title[20];
151 TH1 *hF(NULL);
152 if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
153 TLegend *leg(NULL);
735d88d4 154 TList *l(NULL); TVirtualPad *pad(NULL);
629ae9b8 155 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
156 TObjArray *arr(NULL);
157 switch(ifig){
158 case kNCl: // number of clusters/track
d25116b6 159 if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
94f7dff7 160
161 leg = new TLegend(.83, .7, .99, .96);
162 leg->SetHeader("Species");
163 leg->SetBorderSize(0); leg->SetFillStyle(0);
164 for(Int_t ig(0); ig<fgkNgraph[kNCl]; ig++){
165 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
dfd7d48b 166 if(!g->GetN()) continue;
94f7dff7 167 g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");
168 if(ig) continue;
169 hF=g->GetHistogram();
170 hF->SetXTitle("no of clusters");
171 hF->SetYTitle("entries");
172 hF->GetYaxis()->CenterTitle(1);
173 hF->GetYaxis()->SetTitleOffset(1.2);
174 hF->SetMinimum(5);
175 }
176 leg->Draw(); gPad->SetLogy();
629ae9b8 177 break;
178 case kTRDstat: // Efficiency
179 if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
94f7dff7 180 leg = new TLegend(.62, .77, .98, .98);
629ae9b8 181 leg->SetHeader("TRD Efficiency");
94f7dff7 182 leg->SetBorderSize(0); leg->SetFillStyle(0);
629ae9b8 183 title[0] = "Geometrical (TRDin/TPCout)";
184 title[1] = "Tracking (TRDout/TRDin)";
185 title[2] = "PID (TRDpid/TRDin)";
186 title[3] = "Refit (TRDrefit/TRDin)";
187 hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
94f7dff7 188 hF->SetMaximum(1.4);
629ae9b8 189 hF->GetXaxis()->SetMoreLogLabels();
94f7dff7 190 hF->GetYaxis()->CenterTitle(1);
629ae9b8 191 hF->Draw("p");
192 for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
193 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
194 g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
195 //PutTrendValue(name[id], g->GetMean(2));
196 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
1ee39b3a 197 }
629ae9b8 198 leg->Draw(); gPad->SetLogx();
199 break;
200 case kTRDmom: // Energy loss
201 if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
202 leg = new TLegend(.65, .7, .95, .99);
203 leg->SetHeader("Energy Loss");
204 leg->SetBorderSize(1); leg->SetFillColor(0);
205 title[0] = "Max & 90% quantile";
206 title[1] = "Mean & 60% quantile";
207 hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
208 hF->SetMaximum(1.3);hF->SetMinimum(-.3);
209 hF->Draw("p");
210 for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
211 if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
212 ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
213 //PutTrendValue(name[id], g->GetMean(2));
214 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
1ee39b3a 215 }
629ae9b8 216 leg->Draw();gPad->SetLogx(kFALSE);
217 break;
218 case kPtRes: // Pt resolution @ vertex
219 if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
735d88d4 220 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
221 pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();
222 pad->SetMargin(0.1, 0.022, 0.1, 0.023);
223 hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
224 hF->SetMaximum(10.);hF->SetMinimum(-3.);
225 hF->GetXaxis()->SetMoreLogLabels();
226 hF->GetXaxis()->SetTitleOffset(1.2);
227 hF->GetYaxis()->CenterTitle();
228 hF->Draw("p");
81979445 229 //for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){
230 for(Int_t ig(2); ig<6; ig++){
735d88d4 231 if(!(g = (TGraphErrors*)arr->At(ig))) continue;
232 if(!g->GetN()) continue;
233 g->Draw("pl");
234 //PutTrendValue(name[id], g->GetMean(2));
235 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
236 }
237 pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();
238 pad->SetMargin(0.1, 0.22, 0.1, 0.023);
239 hF = (TH1*)hF->Clone("hFcheckESD1");
240 hF->SetTitle("ITS+TPC");
241 hF->SetMaximum(10.);hF->SetMinimum(-3.);
242 hF->Draw("p");
629ae9b8 243 leg = new TLegend(.78, .1, .99, .98);
244 leg->SetHeader("P_{t} @ DCA");
245 leg->SetBorderSize(1); leg->SetFillColor(0);
246 leg->SetTextAlign(22);
247 leg->SetTextFont(12);
248 leg->SetTextSize(0.03813559);
92d6d80c 249 {
250 Int_t nPlots(0);
81979445 251 //for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){
252 for(Int_t ig(12); ig<16; ig++){
735d88d4 253 if(!(g = (TGraphErrors*)arr->At(ig))) continue;
92d6d80c 254 if(!g->GetN()) continue;
255 nPlots++;
256 g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
257 //PutTrendValue(name[id], g->GetMean(2));
258 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
259 }
260 if(nPlots) leg->Draw();
629ae9b8 261 }
735d88d4 262
629ae9b8 263 break;
1ee39b3a 264 }
e2e3cec2 265 return kTRUE;
266}
267
1ee39b3a 268//____________________________________________________________________
f8f46e4d 269void AliTRDcheckESD::UserExec(Option_t *){
1ee39b3a 270 //
271 // Run the Analysis
272 //
f8f46e4d 273 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
274 fMC = MCEvent();
705f8b0a 275
1ee39b3a 276 if(!fESD){
277 AliError("ESD event missing.");
278 return;
279 }
f8f46e4d 280
1ee39b3a 281 // Get MC information if available
9b496b5e 282 AliStack * fStack = NULL;
1ee39b3a 283 if(HasMC()){
284 if(!fMC){
285 AliWarning("MC event missing");
286 SetMC(kFALSE);
287 } else {
288 if(!(fStack = fMC->Stack())){
289 AliWarning("MC stack missing");
290 SetMC(kFALSE);
291 }
292 }
293 }
9b496b5e 294 TH2 *h(NULL);
295
296 AliESDtrack *esdTrack(NULL);
1ee39b3a 297 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
1ee39b3a 298 esdTrack = fESD->GetTrack(itrk);
299
1ee39b3a 300 // track status
fb09668f 301 ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
36d21092 302
303 // track selection
304 Bool_t selected(kTRUE);
305 if(esdTrack->Pt() < fgkPt){
306 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESD->GetEventNumberInFile(), itrk, esdTrack->Pt()));
307 selected = kFALSE;
308 }
309 if(TMath::Abs(esdTrack->Eta()) > fgkEta){
310 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));
311 selected = kFALSE;
312 }
313 if(!Bool_t(status & AliESDtrack::kTPCout)){
314 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] !TPCout", fESD->GetEventNumberInFile(), itrk));
315 selected = kFALSE;
316 }
317 if(esdTrack->GetKinkIndex(0) > 0){
318 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Kink", fESD->GetEventNumberInFile(), itrk));
319 selected = kFALSE;
320 }
321 if(esdTrack->GetTPCNcls() < fgkNclTPC){
322 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESD->GetEventNumberInFile(), itrk, esdTrack->GetTPCNcls()));
323 selected = kFALSE;
324 }
325 Float_t par[2], cov[3];
326 esdTrack->GetImpactParameters(par, cov);
327 if(IsCollision()){ // cuts on DCA
328 if(TMath::Abs(par[0]) > fgkTrkDCAxy){
329 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));
330 selected = kFALSE;
331 }
332 if(TMath::Abs(par[1]) > fgkTrkDCAz){
333 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));
334 selected = kFALSE;
335 }
336 }
337 if(!selected) continue;
1ee39b3a 338
fb09668f 339 //Int_t nTPC(esdTrack->GetNcls(1));
340 Int_t nTRD(esdTrack->GetNcls(2));
341 Double_t pt(esdTrack->Pt());
342 //Double_t eta(esdTrack->Eta());
343 //Double_t phi(esdTrack->Phi());
1ee39b3a 344 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
345 // pid quality
346 //esdTrack->GetTRDntrackletsPID();
fb09668f 347 Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
1ee39b3a 348
349 // look at external track param
350 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
351 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
1ee39b3a 352
fb09668f 353 Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
9b496b5e 354 // read MC info if available
fb09668f 355 Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
356 AliMCParticle *mcParticle(NULL);
9b496b5e 357 if(HasMC()){
358 AliTrackReference *ref(NULL);
359 Int_t fLabel(esdTrack->GetLabel());
fb09668f 360 Int_t fIdx(TMath::Abs(fLabel));
361 if(fIdx > fStack->GetNtrack()) continue;
9b496b5e 362
fb09668f 363 // read MC particle
364 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
9b496b5e 365 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
366 continue;
367 }
fb09668f 368 pt0 = mcParticle->Pt();
369 eta0 = mcParticle->Eta();
370 phi0 = mcParticle->Phi();
371 kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
372
373 // read track references
9b496b5e 374 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
375 if(!nRefs){
fb09668f 376 AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
9b496b5e 377 continue;
378 }
379 Int_t iref = 0;
380 while(iref<nRefs){
381 ref = mcParticle->GetTrackReference(iref);
382 if(ref->LocalX() > fgkxTPC) break;
383 ref=NULL; iref++;
384 }
385 if(ref){
386 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
387 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
388 }
389 } else { // track stopped in TPC
390 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
391 }
fb09668f 392 ptTRD = ref->Pt();kFOUND=kTRUE;
9b496b5e 393 } else { // use reconstructed values
394 if(op){
395 Double_t x(op->GetX());
396 if(x<fgkxTOF && x>fgkxTPC){
fb09668f 397 ptTRD=op->Pt();
9b496b5e 398 kFOUND=kTRUE;
399 }
400 }
1ee39b3a 401
9b496b5e 402 if(!kFOUND && ip){
fb09668f 403 ptTRD=ip->Pt();
9b496b5e 404 kFOUND=kTRUE;
405 }
1ee39b3a 406 }
407
9b496b5e 408 if(kFOUND){
409 h = (TH2I*)fHistos->At(kTRDstat);
fb09668f 410 if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
411 if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
36d21092 412 if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);
fb09668f 413 if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
414 if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
415 }
36d21092 416 Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)
417 ,sgn(esdTrack->Charge()<0?0:1);
418 if(kBarrel && kPhysPrim) {
fb09668f 419 TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
a47a87c5 420 Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
a47a87c5 421 h3->Fill(pt0, 1.e2*(pt/pt0-1.),
36d21092 422 offset + 2*idx + sgn);
1ee39b3a 423 }
36d21092 424 ((TH1*)fHistos->At(kNCl))->Fill(nTRD, 2*idx + sgn);
1ee39b3a 425 if(ip){
426 h = (TH2I*)fHistos->At(kTRDmom);
fb09668f 427 Float_t pTRD(0.);
1ee39b3a 428 for(Int_t ily=6; ily--;){
fb09668f 429 if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
430 h->Fill(ip->GetP()-pTRD, ily);
1ee39b3a 431 }
432 }
433 }
f8f46e4d 434 PostData(1, fHistos);
1ee39b3a 435}
436
437//____________________________________________________________________
438TObjArray* AliTRDcheckESD::Histos()
439{
440// Retrieve histograms array if already build or build it
441
442 if(fHistos) return fHistos;
443
444 fHistos = new TObjArray(kNhistos);
445 //fHistos->SetOwner(kTRUE);
629ae9b8 446
9b496b5e 447 TH1 *h = NULL;
1ee39b3a 448
36d21092 449 // clusters per track
450 const Int_t kNpt(30);
451 Float_t Pt(0.2);
452 Float_t binsPt[kNpt+1];
453 for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;
454 if(!(h = (TH2S*)gROOT->FindObject("hNCl"))){
455 h = new TH2S("hNCl", "Clusters per TRD track;N_{cl}^{TRD};SPECIES;entries", 60, 0., 180., 10, -0.5, 9.5);
456 TAxis *ay(h->GetYaxis());
457 ay->SetLabelOffset(0.015);
458 for(Int_t i(0); i<AliPID::kSPECIES; i++){
459 ay->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
460 ay->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
461 }
1ee39b3a 462 } else h->Reset();
629ae9b8 463 fHistos->AddAt(h, kNCl); fNRefFigures++;
1ee39b3a 464
465 // status bits histogram
36d21092 466 const Int_t kNbits(5);
467 Float_t Bits(.5);
468 Float_t binsBits[kNbits+1];
fb09668f 469 for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
1ee39b3a 470 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
629ae9b8 471 h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
fb09668f 472 TAxis *ay(h->GetYaxis());
473 ay->SetBinLabel(1, "kTPCout");
474 ay->SetBinLabel(2, "kTRDin");
475 ay->SetBinLabel(3, "kTRDout");
476 ay->SetBinLabel(4, "kTRDpid");
477 ay->SetBinLabel(5, "kTRDrefit");
1ee39b3a 478 } else h->Reset();
479 fHistos->AddAt(h, kTRDstat);
480
481 // energy loss
482 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
fb09668f 483 h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
1ee39b3a 484 } else h->Reset();
485 fHistos->AddAt(h, kTRDmom);
629ae9b8 486 if(!HasMC()) return fHistos;
1ee39b3a 487
fb09668f 488 // pt resolution
a47a87c5 489 const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
490 Float_t DPt(-3.), Spec(-0.5);
fb09668f 491 Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
492 for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
493 for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
494 if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
495 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);
496 TAxis *az(h->GetZaxis());
a47a87c5 497 az->SetLabelOffset(0.015);
fb09668f 498 for(Int_t i(0); i<AliPID::kSPECIES; i++){
a47a87c5 499 az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
500 az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
501 az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
502 az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
fb09668f 503 }
504 } else h->Reset();
505 fHistos->AddAt(h, kPtRes);
506
1ee39b3a 507 return fHistos;
508}
509
510//____________________________________________________________________
97d0ccba 511Bool_t AliTRDcheckESD::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
1ee39b3a 512{
513// Load data from performance file
514
97d0ccba 515 if(!TFile::Open(file)){
516 AliWarning(Form("Couldn't open file %s.", file));
1ee39b3a 517 return kFALSE;
518 }
97d0ccba 519 if(dir){
520 if(!gFile->cd(dir)){
521 AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
522 return kFALSE;
523 }
524 }
525 TObjArray *o(NULL);
526 const Char_t *tn=(name ? name : GetName());
527 if(!(o = (TObjArray*)gDirectory->Get(tn))){
528 AliWarning(Form("Missing histogram container %s.", tn));
1ee39b3a 529 return kFALSE;
530 }
531 fHistos = (TObjArray*)o->Clone(GetName());
532 gFile->Close();
533 return kTRUE;
534}
535
536//_______________________________________________________
537Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
538{
539// Dump trending value to default file
540
541 if(!fgFile){
542 fgFile = fopen("TRD.Performance.txt", "at");
543 }
544 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
545 return kTRUE;
546}
547
548//____________________________________________________________________
549void AliTRDcheckESD::Terminate(Option_t *)
550{
551// Steer post-processing
97d0ccba 552 if(!fHistos){
2efff183 553 fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
a96ac33a 554 if(!fHistos){
555 AliError("Histogram container not found in output");
556 return;
557 }
9b496b5e 558 }
1ee39b3a 559
629ae9b8 560 const Char_t *name[kNrefs] = {
561 "Ncl", "Eff", "Eloss", "PtResDCA"
562 };
563 TObjArray *arr(NULL); TGraph *g(NULL);
564 if(!fResults){
565 fResults = new TObjArray(kNrefs);
566 fResults->SetOwner();
567 fResults->SetName("results");
568 for(Int_t iref(0); iref<kNrefs; iref++){
569 fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
570 arr->SetName(name[iref]); arr->SetOwner();
571 switch(iref){
94f7dff7 572 case kNCl:
573 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
574 arr->AddAt(g = new TGraphErrors(), ig);
575 g->SetLineColor(ig+1);
576 g->SetMarkerColor(ig+1);
577 g->SetMarkerStyle(ig+20);
dfd7d48b 578 g->SetName(Form("s%d", ig));
579 switch(ig){
580 case 0: g->SetTitle("ALL"); break;
581 case 1: g->SetTitle("NEG"); break;
582 case 2: g->SetTitle("POS"); break;
583 default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;
584 };
94f7dff7 585 }
586 break;
629ae9b8 587 case kTRDmom:
588 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
589 arr->AddAt(g = new TGraphAsymmErrors(), ig);
590 g->SetLineColor(ig+1);
591 g->SetMarkerColor(ig+1);
592 g->SetMarkerStyle(ig+20);
593 }
594 break;
595 case kPtRes:
735d88d4 596 for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
597 Int_t ig(2*idx);
598 arr->AddAt(g = new TGraphErrors(), ig);
599 g->SetLineColor(kRed-idx);
600 g->SetMarkerColor(kRed-idx);
601 g->SetMarkerStyle(20+idx);
602 g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
603 arr->AddAt(g = new TGraphErrors(), ig+1);
604 g->SetLineColor(kBlue-idx);
605 g->SetMarkerColor(kBlue-idx);
606 g->SetMarkerStyle(20+idx);
607 g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
608
609 ig+=10;
610 arr->AddAt(g = new TGraphErrors(), ig);
629ae9b8 611 g->SetLineColor(kRed-idx);
612 g->SetMarkerColor(kRed-idx);
735d88d4 613 g->SetMarkerStyle(20+idx);
614 g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
615 arr->AddAt(g = new TGraphErrors(), ig+1);
629ae9b8 616 g->SetLineColor(kBlue-idx);
617 g->SetMarkerColor(kBlue-idx);
735d88d4 618 g->SetMarkerStyle(20+idx);
619 g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
629ae9b8 620 }
621 break;
622 default:
623 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
624 arr->AddAt(g = new TGraphErrors(), ig);
625 g->SetLineColor(ig+1);
626 g->SetMarkerColor(ig+1);
627 g->SetMarkerStyle(ig+20);
628 }
629 break;
630 }
631 }
632 }
d25116b6 633 TH1 *h1[2] = {NULL, NULL};
634 TH2I *h2(NULL);
635 TAxis *ax(NULL);
636
637 // No of clusters
94f7dff7 638 if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;
639 ax = h2->GetXaxis();
d25116b6 640 arr = (TObjArray*)fResults->At(kNCl);
dfd7d48b 641 // All tracks
94f7dff7 642 h1[0] = h2->ProjectionX("Ncl_px");
d25116b6 643 TGraphErrors *ge=(TGraphErrors*)arr->At(0);
d25116b6 644 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
645 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
646 }
dfd7d48b 647 // All charged tracks
648 TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};
649 hNclCh[0]->Reset();hNclCh[1]->Reset();
650 for(Int_t is(1); is<=AliPID::kSPECIES; is++){
651 hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg
652 hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is)); // pos
653 }
654 if(Int_t(hNclCh[0]->GetEntries())){
655 ge=(TGraphErrors*)arr->At(1);
656 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
657 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));
658 }
659 }
660 if(Int_t(hNclCh[1]->GetEntries())){
661 ge=(TGraphErrors*)arr->At(2);
662 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
663 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));
664 }
665 }
666 // Species wise
94f7dff7 667 for(Int_t is(1); is<=AliPID::kSPECIES; is++){
668 h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);
dfd7d48b 669 if(!Int_t(h1[0]->GetEntries())) continue;
670 ge=(TGraphErrors*)arr->At(2+is);
94f7dff7 671 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
672 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
673 }
674 }
629ae9b8 675 fNRefFigures = 1;
676
677 // EFFICIENCY
1ee39b3a 678 // geometrical efficiency
629ae9b8 679 if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
680 arr = (TObjArray*)fResults->At(kTRDstat);
1ee39b3a 681 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
682 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
629ae9b8 683 Process(h1, (TGraphErrors*)arr->At(0));
1ee39b3a 684 delete h1[0];delete h1[1];
1ee39b3a 685 // tracking efficiency
686 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
687 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
629ae9b8 688 Process(h1, (TGraphErrors*)arr->At(1));
1ee39b3a 689 delete h1[1];
1ee39b3a 690 // PID efficiency
691 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
629ae9b8 692 Process(h1, (TGraphErrors*)arr->At(2));
1ee39b3a 693 delete h1[1];
1ee39b3a 694 // Refit efficiency
695 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
629ae9b8 696 Process(h1, (TGraphErrors*)arr->At(3));
1ee39b3a 697 delete h1[1];
629ae9b8 698 fNRefFigures++;
699
700 // ENERGY LOSS
1ee39b3a 701 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
629ae9b8 702 arr = (TObjArray*)fResults->At(kTRDmom);
703 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
d25116b6 704 ax=h2->GetXaxis();
1ee39b3a 705 const Int_t nq(4);
706 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
707 Double_t yq[nq];
708 for(Int_t ily=6; ily--;){
709 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
710 h1[0]->GetQuantiles(nq,yq,xq);
711 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
712 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
713 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
714 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
715
716 //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]);
717 delete h1[0];
718 }
629ae9b8 719 fNRefFigures++;
720 if(!HasMC()) return;
721
722 // Pt RESOLUTION @ DCA
723 TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
724 if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
725 arr = (TObjArray*)fResults->At(kPtRes);
726 TAxis *az(h3->GetZaxis());
727 for(Int_t i(0); i<AliPID::kSPECIES; i++){
735d88d4 728 Int_t idx(2*i);
729 az->SetRange(idx+1, idx+2);
730 gg[1] = (TGraphErrors*)arr->At(idx);
731 gg[0] = (TGraphErrors*)arr->At(idx+1);
732 Process2D((TH2*)h3->Project3D("yx"), gg);
733
734 idx+=10;
735 az->SetRange(idx+1, idx+2);
736 gg[1] = (TGraphErrors*)arr->At(idx);
737 gg[0] = (TGraphErrors*)arr->At(idx+1);
629ae9b8 738 Process2D((TH2*)h3->Project3D("yx"), gg);
629ae9b8 739 }
740 fNRefFigures++;
1ee39b3a 741}
742
fb09668f 743//____________________________________________________________________
744Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
745{
746 switch(pdg){
a47a87c5 747 case kElectron: return AliPID::kElectron;
748 case kMuonMinus: return AliPID::kMuon;
749 case kPiPlus: return AliPID::kPion;
750 case kKPlus: return AliPID::kKaon;
751 case kProton: return AliPID::kProton;
fb09668f 752 }
a47a87c5 753 return -1;
fb09668f 754}
755
1ee39b3a 756//____________________________________________________________________
757void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
758{
759// Generic function to process one reference plot
760
761 Int_t n1 = 0, n2 = 0, ip=0;
762 Double_t eff = 0.;
763
764 TAxis *ax = h1[0]->GetXaxis();
765 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
766 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
767 n2 = (Int_t)h1[1]->GetBinContent(ib);
768 eff = n2/Float_t(n1);
769
770 ip=g->GetN();
771 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
772 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
773 }
774}
629ae9b8 775//________________________________________________________
776void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
777{
778 //
779 // Do the processing
780 //
1ee39b3a 781
629ae9b8 782 Int_t n = 0;
783 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
784 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
785 TF1 f("fg", "gaus", -3.,3.);
786 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
787 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
788 TH1D *h = h2->ProjectionY("py", ibin, ibin);
789 if(h->GetEntries()<100) continue;
790 //AdjustF1(h, f);
791
792 h->Fit(&f, "QN");
793 Int_t ip = g[0]->GetN();
794 g[0]->SetPoint(ip, x, f.GetParameter(1));
795 g[0]->SetPointError(ip, 0., f.GetParError(1));
796 g[1]->SetPoint(ip, x, f.GetParameter(2));
797 g[1]->SetPointError(ip, 0., f.GetParError(2));
798 }
799 return;
800}
1ee39b3a 801//____________________________________________________________________
802void AliTRDcheckESD::PrintStatus(ULong_t status)
803{
804// Dump track status to stdout
805
806 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"
807 ,Bool_t(status & AliESDtrack::kITSin)
808 ,Bool_t(status & AliESDtrack::kITSout)
809 ,Bool_t(status & AliESDtrack::kITSrefit)
810 ,Bool_t(status & AliESDtrack::kTPCin)
811 ,Bool_t(status & AliESDtrack::kTPCout)
812 ,Bool_t(status & AliESDtrack::kTPCrefit)
813 ,Bool_t(status & AliESDtrack::kTPCpid)
814 ,Bool_t(status & AliESDtrack::kTRDin)
815 ,Bool_t(status & AliESDtrack::kTRDout)
816 ,Bool_t(status & AliESDtrack::kTRDrefit)
817 ,Bool_t(status & AliESDtrack::kTRDpid)
818 ,Bool_t(status & AliESDtrack::kTRDStop)
819 ,Bool_t(status & AliESDtrack::kHMPIDout)
820 ,Bool_t(status & AliESDtrack::kHMPIDpid)
821 );
822}
823