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