1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////
18 // Check basic detector results at ESD level
19 // - Geometrical efficiency
20 // - Tracking efficiency
25 // Alex Bercuci <A.Bercuci@gsi.de>
27 //////////////////////////////////////////////////////
29 #include <TClonesArray.h>
30 #include <TObjArray.h>
36 #include <TGraphErrors.h>
37 #include <TGraphAsymmErrors.h>
42 #include <TParticle.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"
52 #include "AliESDtrack.h"
53 #include "AliMCParticle.h"
56 #include "AliTrackReference.h"
58 #include "AliTRDcheckESD.h"
60 ClassImp(AliTRDcheckESD)
62 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
63 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
64 const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
66 FILE* AliTRDcheckESD::fgFile = NULL;
68 const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;
69 const Int_t AliTRDcheckESD::fgkEvVertexN = 1;
70 const Float_t AliTRDcheckESD::fgkTrkDCAxy = 40.;
71 const Float_t AliTRDcheckESD::fgkTrkDCAz = 15.;
72 const Int_t AliTRDcheckESD::fgkNclTPC = 100;
73 const Float_t AliTRDcheckESD::fgkPt = 0.2;
74 const Float_t AliTRDcheckESD::fgkEta = 0.9;
76 //____________________________________________________________________
77 AliTRDcheckESD::AliTRDcheckESD():
87 // Default constructor
89 SetNameTitle("checkESD", "Check TRD @ ESD level");
93 //____________________________________________________________________
94 AliTRDcheckESD::AliTRDcheckESD(char* name):
95 AliAnalysisTaskSE(name)
104 // Default constructor
107 SetTitle("Check TRD @ ESD level");
108 DefineOutput(1, TObjArray::Class());
111 //____________________________________________________________________
112 AliTRDcheckESD::~AliTRDcheckESD()
125 //____________________________________________________________________
126 void AliTRDcheckESD::UserCreateOutputObjects()
129 // Create Output Containers (TObjectArray containing 1D histograms)
135 //____________________________________________________________________
136 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
138 if(ifig>=fNRefFigures){
139 AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
143 AliWarning("Please provide a canvas to draw results.");
146 gPad->SetLogx(0);gPad->SetLogy(0);
147 gPad->SetMargin(0.125, 0.015, 0.1, 0.015);
150 const Char_t *title[20];
152 if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
154 TList *l(NULL); TVirtualPad *pad(NULL);
155 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
156 TObjArray *arr(NULL);
158 case kNCl: // number of clusters/track
159 if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
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;
166 if(!g->GetN()) continue;
167 g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");
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);
176 leg->Draw(); gPad->SetLogy();
178 case kTRDstat: // Efficiency
179 if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
180 leg = new TLegend(.62, .77, .98, .98);
181 leg->SetHeader("TRD Efficiency");
182 leg->SetBorderSize(0); leg->SetFillStyle(0);
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.);
189 hF->GetXaxis()->SetMoreLogLabels();
190 hF->GetYaxis()->CenterTitle(1);
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));
198 leg->Draw(); gPad->SetLogx();
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);
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));
216 leg->Draw();gPad->SetLogx(kFALSE);
218 case kPtRes: // Pt resolution @ vertex
219 if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
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();
229 //for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){
230 for(Int_t ig(2); ig<6; ig++){
231 if(!(g = (TGraphErrors*)arr->At(ig))) continue;
232 if(!g->GetN()) continue;
234 //PutTrendValue(name[id], g->GetMean(2));
235 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
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.);
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);
251 //for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){
252 for(Int_t ig(12); ig<16; ig++){
253 if(!(g = (TGraphErrors*)arr->At(ig))) continue;
254 if(!g->GetN()) continue;
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));
260 if(nPlots) leg->Draw();
268 //____________________________________________________________________
269 void AliTRDcheckESD::UserExec(Option_t *){
273 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
277 AliError("ESD event missing.");
281 // Get MC information if available
282 AliStack * fStack = NULL;
285 AliWarning("MC event missing");
288 if(!(fStack = fMC->Stack())){
289 AliWarning("MC stack missing");
296 AliESDtrack *esdTrack(NULL);
297 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
298 esdTrack = fESD->GetTrack(itrk);
301 ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
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()));
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())));
313 if(!Bool_t(status & AliESDtrack::kTPCout)){
314 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] !TPCout", fESD->GetEventNumberInFile(), itrk));
317 if(esdTrack->GetKinkIndex(0) > 0){
318 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Kink", fESD->GetEventNumberInFile(), itrk));
321 if(esdTrack->GetTPCNcls() < fgkNclTPC){
322 AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESD->GetEventNumberInFile(), itrk, esdTrack->GetTPCNcls()));
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])));
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])));
337 if(!selected) continue;
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());
344 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
346 //esdTrack->GetTRDntrackletsPID();
347 Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
349 // look at external track param
350 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
351 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
353 Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
354 // read MC info if available
355 Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
356 AliMCParticle *mcParticle(NULL);
358 AliTrackReference *ref(NULL);
359 Int_t fLabel(esdTrack->GetLabel());
360 Int_t fIdx(TMath::Abs(fLabel));
361 if(fIdx > fStack->GetNtrack()) continue;
364 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
365 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
368 pt0 = mcParticle->Pt();
369 eta0 = mcParticle->Eta();
370 phi0 = mcParticle->Phi();
371 kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
373 // read track references
374 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
376 AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
381 ref = mcParticle->GetTrackReference(iref);
382 if(ref->LocalX() > fgkxTPC) break;
386 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
387 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
389 } else { // track stopped in TPC
390 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
392 ptTRD = ref->Pt();kFOUND=kTRUE;
393 } else { // use reconstructed values
395 Double_t x(op->GetX());
396 if(x<fgkxTOF && x>fgkxTPC){
409 h = (TH2I*)fHistos->At(kTRDstat);
410 if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
411 if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
412 if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);
413 if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
414 if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
416 Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)
417 ,sgn(esdTrack->Charge()<0?0:1);
418 if(kBarrel && kPhysPrim) {
419 TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
420 Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
421 h3->Fill(pt0, 1.e2*(pt/pt0-1.),
422 offset + 2*idx + sgn);
424 ((TH1*)fHistos->At(kNCl))->Fill(nTRD, 2*idx + sgn);
426 h = (TH2I*)fHistos->At(kTRDmom);
428 for(Int_t ily=6; ily--;){
429 if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
430 h->Fill(ip->GetP()-pTRD, ily);
434 PostData(1, fHistos);
437 //____________________________________________________________________
438 TObjArray* AliTRDcheckESD::Histos()
440 // Retrieve histograms array if already build or build it
442 if(fHistos) return fHistos;
444 fHistos = new TObjArray(kNhistos);
445 //fHistos->SetOwner(kTRUE);
449 // clusters per track
450 const Int_t kNpt(30);
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)));
463 fHistos->AddAt(h, kNCl); fNRefFigures++;
465 // status bits histogram
466 const Int_t kNbits(5);
468 Float_t binsBits[kNbits+1];
469 for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
470 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
471 h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
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");
479 fHistos->AddAt(h, kTRDstat);
482 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
483 h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
485 fHistos->AddAt(h, kTRDmom);
486 if(!HasMC()) return fHistos;
489 const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
490 Float_t DPt(-3.), Spec(-0.5);
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());
497 az->SetLabelOffset(0.015);
498 for(Int_t i(0); i<AliPID::kSPECIES; i++){
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)));
505 fHistos->AddAt(h, kPtRes);
510 //____________________________________________________________________
511 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
513 // Load data from performance file
515 if(!TFile::Open(filename)){
516 AliWarning(Form("Couldn't open file %s.", filename));
520 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
521 AliWarning("Missing histogram container.");
524 fHistos = (TObjArray*)o->Clone(GetName());
526 SETBIT(fStatus, kLoad);
530 //_______________________________________________________
531 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
533 // Dump trending value to default file
536 fgFile = fopen("TRD.Performance.txt", "at");
538 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
542 //____________________________________________________________________
543 void AliTRDcheckESD::Terminate(Option_t *)
545 // Steer post-processing
547 fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
549 AliError("Histogram container not found in output");
554 const Char_t *name[kNrefs] = {
555 "Ncl", "Eff", "Eloss", "PtResDCA"
557 TObjArray *arr(NULL); TGraph *g(NULL);
559 fResults = new TObjArray(kNrefs);
560 fResults->SetOwner();
561 fResults->SetName("results");
562 for(Int_t iref(0); iref<kNrefs; iref++){
563 fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
564 arr->SetName(name[iref]); arr->SetOwner();
567 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
568 arr->AddAt(g = new TGraphErrors(), ig);
569 g->SetLineColor(ig+1);
570 g->SetMarkerColor(ig+1);
571 g->SetMarkerStyle(ig+20);
572 g->SetName(Form("s%d", ig));
574 case 0: g->SetTitle("ALL"); break;
575 case 1: g->SetTitle("NEG"); break;
576 case 2: g->SetTitle("POS"); break;
577 default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;
582 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
583 arr->AddAt(g = new TGraphAsymmErrors(), ig);
584 g->SetLineColor(ig+1);
585 g->SetMarkerColor(ig+1);
586 g->SetMarkerStyle(ig+20);
590 for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
592 arr->AddAt(g = new TGraphErrors(), ig);
593 g->SetLineColor(kRed-idx);
594 g->SetMarkerColor(kRed-idx);
595 g->SetMarkerStyle(20+idx);
596 g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
597 arr->AddAt(g = new TGraphErrors(), ig+1);
598 g->SetLineColor(kBlue-idx);
599 g->SetMarkerColor(kBlue-idx);
600 g->SetMarkerStyle(20+idx);
601 g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
604 arr->AddAt(g = new TGraphErrors(), ig);
605 g->SetLineColor(kRed-idx);
606 g->SetMarkerColor(kRed-idx);
607 g->SetMarkerStyle(20+idx);
608 g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
609 arr->AddAt(g = new TGraphErrors(), ig+1);
610 g->SetLineColor(kBlue-idx);
611 g->SetMarkerColor(kBlue-idx);
612 g->SetMarkerStyle(20+idx);
613 g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
617 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
618 arr->AddAt(g = new TGraphErrors(), ig);
619 g->SetLineColor(ig+1);
620 g->SetMarkerColor(ig+1);
621 g->SetMarkerStyle(ig+20);
627 TH1 *h1[2] = {NULL, NULL};
632 if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;
634 arr = (TObjArray*)fResults->At(kNCl);
636 h1[0] = h2->ProjectionX("Ncl_px");
637 TGraphErrors *ge=(TGraphErrors*)arr->At(0);
638 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
639 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
641 // All charged tracks
642 TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};
643 hNclCh[0]->Reset();hNclCh[1]->Reset();
644 for(Int_t is(1); is<=AliPID::kSPECIES; is++){
645 hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg
646 hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is)); // pos
648 if(Int_t(hNclCh[0]->GetEntries())){
649 ge=(TGraphErrors*)arr->At(1);
650 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
651 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));
654 if(Int_t(hNclCh[1]->GetEntries())){
655 ge=(TGraphErrors*)arr->At(2);
656 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
657 ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));
661 for(Int_t is(1); is<=AliPID::kSPECIES; is++){
662 h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);
663 if(!Int_t(h1[0]->GetEntries())) continue;
664 ge=(TGraphErrors*)arr->At(2+is);
665 for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
666 ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
672 // geometrical efficiency
673 if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
674 arr = (TObjArray*)fResults->At(kTRDstat);
675 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
676 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
677 Process(h1, (TGraphErrors*)arr->At(0));
678 delete h1[0];delete h1[1];
679 // tracking efficiency
680 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
681 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
682 Process(h1, (TGraphErrors*)arr->At(1));
685 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
686 Process(h1, (TGraphErrors*)arr->At(2));
689 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
690 Process(h1, (TGraphErrors*)arr->At(3));
695 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
696 arr = (TObjArray*)fResults->At(kTRDmom);
697 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
700 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
702 for(Int_t ily=6; ily--;){
703 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
704 h1[0]->GetQuantiles(nq,yq,xq);
705 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
706 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
707 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
708 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
710 //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]);
716 // Pt RESOLUTION @ DCA
717 TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
718 if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
719 arr = (TObjArray*)fResults->At(kPtRes);
720 TAxis *az(h3->GetZaxis());
721 for(Int_t i(0); i<AliPID::kSPECIES; i++){
723 az->SetRange(idx+1, idx+2);
724 gg[1] = (TGraphErrors*)arr->At(idx);
725 gg[0] = (TGraphErrors*)arr->At(idx+1);
726 Process2D((TH2*)h3->Project3D("yx"), gg);
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);
737 //____________________________________________________________________
738 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
741 case kElectron: return AliPID::kElectron;
742 case kMuonMinus: return AliPID::kMuon;
743 case kPiPlus: return AliPID::kPion;
744 case kKPlus: return AliPID::kKaon;
745 case kProton: return AliPID::kProton;
750 //____________________________________________________________________
751 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
753 // Generic function to process one reference plot
755 Int_t n1 = 0, n2 = 0, ip=0;
758 TAxis *ax = h1[0]->GetXaxis();
759 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
760 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
761 n2 = (Int_t)h1[1]->GetBinContent(ib);
762 eff = n2/Float_t(n1);
765 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
766 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
769 //________________________________________________________
770 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
777 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
778 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
779 TF1 f("fg", "gaus", -3.,3.);
780 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
781 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
782 TH1D *h = h2->ProjectionY("py", ibin, ibin);
783 if(h->GetEntries()<100) continue;
787 Int_t ip = g[0]->GetN();
788 g[0]->SetPoint(ip, x, f.GetParameter(1));
789 g[0]->SetPointError(ip, 0., f.GetParError(1));
790 g[1]->SetPoint(ip, x, f.GetParameter(2));
791 g[1]->SetPointError(ip, 0., f.GetParError(2));
795 //____________________________________________________________________
796 void AliTRDcheckESD::PrintStatus(ULong_t status)
798 // Dump track status to stdout
800 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"
801 ,Bool_t(status & AliESDtrack::kITSin)
802 ,Bool_t(status & AliESDtrack::kITSout)
803 ,Bool_t(status & AliESDtrack::kITSrefit)
804 ,Bool_t(status & AliESDtrack::kTPCin)
805 ,Bool_t(status & AliESDtrack::kTPCout)
806 ,Bool_t(status & AliESDtrack::kTPCrefit)
807 ,Bool_t(status & AliESDtrack::kTPCpid)
808 ,Bool_t(status & AliESDtrack::kTRDin)
809 ,Bool_t(status & AliESDtrack::kTRDout)
810 ,Bool_t(status & AliESDtrack::kTRDrefit)
811 ,Bool_t(status & AliESDtrack::kTRDpid)
812 ,Bool_t(status & AliESDtrack::kTRDStop)
813 ,Bool_t(status & AliESDtrack::kHMPIDout)
814 ,Bool_t(status & AliESDtrack::kHMPIDpid)