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 //____________________________________________________________________
69 AliTRDcheckESD::AliTRDcheckESD():
79 // Default constructor
81 SetNameTitle("checkESD", "Check TRD @ ESD level");
85 //____________________________________________________________________
86 AliTRDcheckESD::AliTRDcheckESD(char* name):
87 AliAnalysisTaskSE(name)
96 // Default constructor
99 SetTitle("Check TRD @ ESD level");
100 DefineOutput(1, TObjArray::Class());
103 //____________________________________________________________________
104 AliTRDcheckESD::~AliTRDcheckESD()
117 //____________________________________________________________________
118 void AliTRDcheckESD::UserCreateOutputObjects()
121 // Create Output Containers (TObjectArray containing 1D histograms)
123 //OpenFile(0, "RECREATE");
129 //____________________________________________________________________
130 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
132 if(ifig>=fNRefFigures){
133 AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
137 AliWarning("Please provide a canvas to draw results.");
141 const Char_t *title[20];
143 if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
145 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
146 TObjArray *arr(NULL);
148 case kNCl: // number of clusters/track
149 ((TH1I*)fResults->At(kNCl))->Draw("c");
151 case kTRDstat: // Efficiency
152 if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
153 leg = new TLegend(.65, .7, .95, .99);
154 leg->SetHeader("TRD Efficiency");
155 leg->SetBorderSize(1); leg->SetFillColor(0);
156 title[0] = "Geometrical (TRDin/TPCout)";
157 title[1] = "Tracking (TRDout/TRDin)";
158 title[2] = "PID (TRDpid/TRDin)";
159 title[3] = "Refit (TRDrefit/TRDin)";
160 hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
162 hF->GetXaxis()->SetMoreLogLabels();
164 for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
165 if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
166 g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
167 //PutTrendValue(name[id], g->GetMean(2));
168 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
170 leg->Draw(); gPad->SetLogx();
172 case kTRDmom: // Energy loss
173 if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
174 leg = new TLegend(.65, .7, .95, .99);
175 leg->SetHeader("Energy Loss");
176 leg->SetBorderSize(1); leg->SetFillColor(0);
177 title[0] = "Max & 90% quantile";
178 title[1] = "Mean & 60% quantile";
179 hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
180 hF->SetMaximum(1.3);hF->SetMinimum(-.3);
182 for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
183 if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
184 ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
185 //PutTrendValue(name[id], g->GetMean(2));
186 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
188 leg->Draw();gPad->SetLogx(kFALSE);
190 case kPtRes: // Pt resolution @ vertex
191 if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
192 gPad->SetMargin(0.1, 0.22, 0.1, 0.023);
193 leg = new TLegend(.78, .1, .99, .98);
194 leg->SetHeader("P_{t} @ DCA");
195 leg->SetBorderSize(1); leg->SetFillColor(0);
196 leg->SetTextAlign(22);
197 leg->SetTextFont(12);
198 leg->SetTextSize(0.03813559);
199 hF = new TH1S("hFcheckESD", ";p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
200 hF->SetMaximum(10.);hF->SetMinimum(-3.);
201 hF->GetXaxis()->SetMoreLogLabels();
202 hF->GetXaxis()->SetTitleOffset(1.2);
203 hF->GetYaxis()->CenterTitle();
207 for(Int_t ig(0); ig<fgkNgraph[kPtRes]; ig++){
208 if(!(g = (TGraphErrors*)arr->At(ig))) continue;//return kFALSE;
209 if(!g->GetN()) continue;
211 g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
212 //PutTrendValue(name[id], g->GetMean(2));
213 //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
215 if(nPlots) leg->Draw();
223 //____________________________________________________________________
224 void AliTRDcheckESD::UserExec(Option_t *){
228 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
232 AliError("ESD event missing.");
236 // Get MC information if available
237 AliStack * fStack = NULL;
240 AliWarning("MC event missing");
243 if(!(fStack = fMC->Stack())){
244 AliWarning("MC stack missing");
251 AliESDtrack *esdTrack(NULL);
252 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
253 esdTrack = fESD->GetTrack(itrk);
256 ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
257 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
258 if(esdTrack->GetKinkIndex(0) > 0) continue;
260 //Int_t nTPC(esdTrack->GetNcls(1));
261 Int_t nTRD(esdTrack->GetNcls(2));
262 Double_t pt(esdTrack->Pt());
263 //Double_t eta(esdTrack->Eta());
264 //Double_t phi(esdTrack->Phi());
265 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
267 //esdTrack->GetTRDntrackletsPID();
268 Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
270 // look at external track param
271 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
272 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
274 Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
275 // read MC info if available
276 Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
277 AliMCParticle *mcParticle(NULL);
279 AliTrackReference *ref(NULL);
280 Int_t fLabel(esdTrack->GetLabel());
281 Int_t fIdx(TMath::Abs(fLabel));
282 if(fIdx > fStack->GetNtrack()) continue;
285 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
286 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
289 pt0 = mcParticle->Pt();
290 eta0 = mcParticle->Eta();
291 phi0 = mcParticle->Phi();
292 kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
294 // read track references
295 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
297 AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
302 ref = mcParticle->GetTrackReference(iref);
303 if(ref->LocalX() > fgkxTPC) break;
307 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
308 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
310 } else { // track stopped in TPC
311 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
313 ptTRD = ref->Pt();kFOUND=kTRUE;
314 } else { // use reconstructed values
316 Double_t x(op->GetX());
317 if(x<fgkxTOF && x>fgkxTPC){
330 h = (TH2I*)fHistos->At(kTRDstat);
331 if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
332 if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
333 if(kBarrel && (status & AliESDtrack::kTRDout)){
334 ((TH1*)fHistos->At(kNCl))->Fill(nTRD);
335 h->Fill(ptTRD, kTRDout);
337 if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
338 if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
341 kBarrel && kPhysPrim &&
342 TMath::Abs(eta0) < 0.9) {
343 TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
344 Int_t sgn = mcParticle->Charge()<0?0:1;
345 Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
347 h3->Fill(pt0, 1.e2*(pt/pt0-1.),
348 offset + 2*Pdg2Idx(TMath::Abs(mcParticle->PdgCode())) + sgn);
351 h = (TH2I*)fHistos->At(kTRDmom);
353 for(Int_t ily=6; ily--;){
354 if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
355 h->Fill(ip->GetP()-pTRD, ily);
359 PostData(1, fHistos);
362 //____________________________________________________________________
363 TObjArray* AliTRDcheckESD::Histos()
365 // Retrieve histograms array if already build or build it
367 if(fHistos) return fHistos;
369 fHistos = new TObjArray(kNhistos);
370 //fHistos->SetOwner(kTRUE);
374 // clusters per tracklet
375 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
376 h = new TH1I("hNCl", "Clusters per TRD track;N_{cl}^{TRD};entries", 100, 0., 200.);
378 fHistos->AddAt(h, kNCl); fNRefFigures++;
380 // status bits histogram
381 const Int_t kNpt(30), kNbits(5);
382 Float_t Pt(0.2), Bits(.5);
383 Float_t binsPt[kNpt+1], binsBits[kNbits+1];
384 for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;
385 for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
386 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
387 h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
388 TAxis *ay(h->GetYaxis());
389 ay->SetBinLabel(1, "kTPCout");
390 ay->SetBinLabel(2, "kTRDin");
391 ay->SetBinLabel(3, "kTRDout");
392 ay->SetBinLabel(4, "kTRDpid");
393 ay->SetBinLabel(5, "kTRDrefit");
395 fHistos->AddAt(h, kTRDstat);
398 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
399 h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
401 fHistos->AddAt(h, kTRDmom);
402 if(!HasMC()) return fHistos;
405 const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
406 Float_t DPt(-3.), Spec(-0.5);
407 Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
408 for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
409 for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
410 if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
411 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);
412 TAxis *az(h->GetZaxis());
413 az->SetLabelOffset(0.015);
414 for(Int_t i(0); i<AliPID::kSPECIES; i++){
415 az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
416 az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
417 az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
418 az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
421 fHistos->AddAt(h, kPtRes);
426 //____________________________________________________________________
427 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
429 // Load data from performance file
431 if(!TFile::Open(filename)){
432 AliWarning(Form("Couldn't open file %s.", filename));
436 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
437 AliWarning("Missing histogram container.");
440 fHistos = (TObjArray*)o->Clone(GetName());
442 SETBIT(fStatus, kLoad);
446 //_______________________________________________________
447 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
449 // Dump trending value to default file
452 fgFile = fopen("TRD.Performance.txt", "at");
454 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
458 //____________________________________________________________________
459 void AliTRDcheckESD::Terminate(Option_t *)
461 // Steer post-processing
463 fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
465 AliError("Histogram container not found in output");
470 const Char_t *name[kNrefs] = {
471 "Ncl", "Eff", "Eloss", "PtResDCA"
473 TObjArray *arr(NULL); TGraph *g(NULL);
475 fResults = new TObjArray(kNrefs);
476 fResults->SetOwner();
477 fResults->SetName("results");
478 for(Int_t iref(0); iref<kNrefs; iref++){
479 fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
480 arr->SetName(name[iref]); arr->SetOwner();
483 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
484 arr->AddAt(g = new TGraphAsymmErrors(), ig);
485 g->SetLineColor(ig+1);
486 g->SetMarkerColor(ig+1);
487 g->SetMarkerStyle(ig+20);
491 for(Int_t im(fgkNgraph[iref]/2), ig(im), idx(0); ig--; idx+=ig){
492 arr->AddAt(g = new TGraphErrors(), ig);
493 g->SetLineColor(kRed-idx);
494 g->SetMarkerColor(kRed-idx);
495 g->SetMarkerStyle(20+ig);
496 g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(ig)));
497 arr->AddAt(g = new TGraphErrors(), im+ig);
498 g->SetLineColor(kBlue-idx);
499 g->SetMarkerColor(kBlue-idx);
500 g->SetMarkerStyle(20+ig);
501 g->SetNameTitle(Form("m%d", ig), Form("sys %s", AliPID::ParticleLatexName(ig)));
505 for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
506 arr->AddAt(g = new TGraphErrors(), ig);
507 g->SetLineColor(ig+1);
508 g->SetMarkerColor(ig+1);
509 g->SetMarkerStyle(ig+20);
518 // geometrical efficiency
520 if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
521 arr = (TObjArray*)fResults->At(kTRDstat);
522 TH1 *h1[2] = {NULL, NULL};
523 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
524 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
525 Process(h1, (TGraphErrors*)arr->At(0));
526 delete h1[0];delete h1[1];
527 // tracking efficiency
528 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
529 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
530 Process(h1, (TGraphErrors*)arr->At(1));
533 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
534 Process(h1, (TGraphErrors*)arr->At(2));
537 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
538 Process(h1, (TGraphErrors*)arr->At(3));
543 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
544 arr = (TObjArray*)fResults->At(kTRDmom);
545 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
546 TAxis *ax=h2->GetXaxis();
548 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
550 for(Int_t ily=6; ily--;){
551 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
552 h1[0]->GetQuantiles(nq,yq,xq);
553 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
554 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
555 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
556 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
558 //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]);
564 // Pt RESOLUTION @ DCA
565 TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
566 if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
567 arr = (TObjArray*)fResults->At(kPtRes);
568 TAxis *az(h3->GetZaxis());
569 for(Int_t i(0); i<AliPID::kSPECIES; i++){
570 az->SetRange(5-i, 5-i);
571 gg[1] = (TGraphErrors*)arr->At(4-i);
572 gg[0] = (TGraphErrors*)arr->At(9-i);
573 Process2D((TH2*)h3->Project3D("yx"), gg);
574 //az->SetRange(7+i, 7+i);
575 //Process2D((TH2*)h3->Project3D("yx"), gg);
580 //____________________________________________________________________
581 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
584 case kElectron: return AliPID::kElectron;
585 case kMuonMinus: return AliPID::kMuon;
586 case kPiPlus: return AliPID::kPion;
587 case kKPlus: return AliPID::kKaon;
588 case kProton: return AliPID::kProton;
593 //____________________________________________________________________
594 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
596 // Generic function to process one reference plot
598 Int_t n1 = 0, n2 = 0, ip=0;
601 TAxis *ax = h1[0]->GetXaxis();
602 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
603 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
604 n2 = (Int_t)h1[1]->GetBinContent(ib);
605 eff = n2/Float_t(n1);
608 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
609 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
612 //________________________________________________________
613 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
620 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
621 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
622 TF1 f("fg", "gaus", -3.,3.);
623 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
624 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
625 TH1D *h = h2->ProjectionY("py", ibin, ibin);
626 if(h->GetEntries()<100) continue;
630 Int_t ip = g[0]->GetN();
631 g[0]->SetPoint(ip, x, f.GetParameter(1));
632 g[0]->SetPointError(ip, 0., f.GetParError(1));
633 g[1]->SetPoint(ip, x, f.GetParameter(2));
634 g[1]->SetPointError(ip, 0., f.GetParError(2));
638 //____________________________________________________________________
639 void AliTRDcheckESD::PrintStatus(ULong_t status)
641 // Dump track status to stdout
643 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"
644 ,Bool_t(status & AliESDtrack::kITSin)
645 ,Bool_t(status & AliESDtrack::kITSout)
646 ,Bool_t(status & AliESDtrack::kITSrefit)
647 ,Bool_t(status & AliESDtrack::kTPCin)
648 ,Bool_t(status & AliESDtrack::kTPCout)
649 ,Bool_t(status & AliESDtrack::kTPCrefit)
650 ,Bool_t(status & AliESDtrack::kTPCpid)
651 ,Bool_t(status & AliESDtrack::kTRDin)
652 ,Bool_t(status & AliESDtrack::kTRDout)
653 ,Bool_t(status & AliESDtrack::kTRDrefit)
654 ,Bool_t(status & AliESDtrack::kTRDpid)
655 ,Bool_t(status & AliESDtrack::kTRDStop)
656 ,Bool_t(status & AliESDtrack::kHMPIDout)
657 ,Bool_t(status & AliESDtrack::kHMPIDpid)