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