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 /* $Id: AliTRDefficiency.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Reconstruction QA //
23 // Markus Fasel <M.Fasel@gsi.de> //
25 ////////////////////////////////////////////////////////////////////////////
30 #include <TClonesArray.h>
31 #include <TObjArray.h>
36 #include <THnSparse.h>
40 #include "TTreeStream.h"
43 #include "AliESDtrack.h"
44 #include "AliTrackReference.h"
45 #include "AliExternalTrackParam.h"
46 #include "AliTracker.h"
47 #include "AliAnalysisManager.h"
49 #include "AliTRDgeometry.h"
50 #include "AliTRDtrackV1.h"
51 #include "Cal/AliTRDCalPID.h"
52 #include "AliTRDefficiency.h"
53 #include "info/AliTRDtrackInfo.h"
55 ClassImp(AliTRDefficiency)
57 //____________________________________________________________________
58 AliTRDefficiency::AliTRDefficiency()
64 // Default constructor
66 SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker");
69 //____________________________________________________________________
70 AliTRDefficiency::AliTRDefficiency(char* name)
71 :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker")
76 // Default constructor
80 //____________________________________________________________________
81 AliTRDefficiency::~AliTRDefficiency()
90 // //____________________________________________________________________
91 // void AliTRDefficiency::UserCreateOutputObjects()
94 // // Create output objects
97 // const Int_t nbins = AliTRDCalPID::kNMom;
98 // Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
101 // fContainer = new TObjArray(); fContainer->SetOwner();
102 // for(Int_t is=0; is<AliPID::kSPECIES; is++){
103 // fContainer->Add(h = new TProfile(Form("h%s", AliTRDCalPID::GetPartSymb(is)), AliPID::ParticleShortName(is), nbins, xbins));
104 // h->SetLineColor(AliTRDCalPID::GetPartColor(is));
105 // h->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
106 // h->SetMarkerStyle(24);
108 // fContainer->Add(h = new TProfile("h", "", nbins, xbins));
109 // h->SetMarkerStyle(7);
110 // PostData(1, fContainer);
113 //____________________________________________________________________
114 TH1* AliTRDefficiency::PlotBasicEff(const AliTRDtrackV1 *track)
116 // plot TRD efficiency based on ESD info
119 AliDebug(4, "No ESD info.");
124 if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hEFF")))){
125 AliWarning("No output container defined.");
128 if(track) fkTrack = track;
130 Double_t val[7]; memset(val, 0, 7*sizeof(Double_t));
131 ULong_t status(fkESD->GetStatus());
132 val[0] =((status&AliESDtrack::kTRDin)?1:0) +
133 ((status&AliESDtrack::kTRDStop)?1:0) +
134 ((status&AliESDtrack::kTRDout)?2:0);
135 val[1] = fPhi;//fkESD->Phi();
136 val[2] = fEta;//fkESD->Eta();
137 val[3] = GetPtBin(fPt/*fkESD->Pt()*/);
140 if(fkMC->GetLabel() == fkMC->GetTRDlabel()) val[4] = 0.;
141 else if(fkMC->GetLabel() == -fkMC->GetTRDlabel()) val[4] = 1.;
144 val[5] = fkTrack?fkTrack->GetNumberOfTracklets():0;
145 // down scale PID resolution
146 Int_t spc(fSpecies); if(spc==-6) spc=0; if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
149 if(fkTrack) AliDebug(2, Form("%3d[%s] tracklets[%d] label[%2d %+2d] species[%d %d] in[%c] out[%c] stop[%c]",
150 fkESD->GetId(), fkTrack->IsTPCseeded()?"TPC":"ITS", Int_t(val[5]), fkMC->GetLabel(), fkMC->GetTRDlabel(), fSpecies, spc,
151 (status&AliESDtrack::kTRDin)?'x':'o', (status&AliESDtrack::kTRDout)?'x':'o', (status&AliESDtrack::kTRDStop)?'x':'o'));
156 // //____________________________________________________________________
157 // TH1* AliTRDefficiency::PlotMC(const AliTRDtrackV1 *track)
159 // // plot TRD efficiency based on MC info
161 // if(!HasMC()) return NULL;
163 // AliDebug(4, "No ESD info.");
167 // AliDebug(4, "No MC info.");
171 // THnSparse *H(NULL);
172 // if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hMC")))){
173 // AliWarning("No output container defined.");
176 // if(track) fkTrack = track;
177 // Double_t val[11]; memset(val, 0, 11*sizeof(Double_t));
178 // ULong_t status(fkESD->GetStatus());
179 // val[0] =((status&AliESDtrack::kTRDin)?1:0) +
180 // ((status&AliESDtrack::kTRDStop)?1:0) +
181 // ((status&AliESDtrack::kTRDout)?2:0);
182 // val[1] = fkESD->Phi();
183 // val[2] = fkESD->Eta();
184 // val[3] = DebugLevel()>=1?GetPtBin(fkESD->Pt()):GetPtBinSignificant(fkESD->Pt());
185 // if(fkTrack){ // read track status in debug mode with friends
186 // val[4] = fkTrack->GetStatusTRD(-1);
187 // for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++) val[5+ily]=fkTrack->GetStatusTRD(ily);
193 //____________________________________________________________________
194 void AliTRDefficiency::LocalUserExec(Option_t *)
200 Int_t labelsacc[10000];
201 memset(labelsacc, 0, sizeof(Int_t) * 10000);
203 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
205 if(!fTracks->GetEntriesFast()) return;
206 else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));
208 fMissed = new TClonesArray("AliTRDtrackInfo", 10);
213 Int_t selection[10000], nselect = 0;
214 ULong_t status; Int_t pidx;
215 Int_t nTRD = 0, nTPC = 0, nMiss = 0;
216 AliTRDtrackInfo *track = NULL;
217 AliTrackReference *ref = NULL;
218 AliExternalTrackParam *esd = NULL;
219 for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
220 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
222 if(!track->HasESDtrack()) continue;
223 status = track->GetStatus();
225 // missing TPC propagation - interesting for SA
226 if(!(status&AliESDtrack::kTPCout)) continue;
229 if(HasMCdata() && track->GetNTrackRefs() <= 1) continue;
232 selection[nselect++]=itrk;
233 ref = track->GetTrackRef(0);
234 esd = track->GetESDinfo()->GetOuterParam();
235 mom = ref ? ref->P(): esd->P();
236 pidx = AliTRDCalPID::GetPartIndex(track->GetPDG());
237 pidx = TMath::Max(pidx, 0);
238 AliDebug(4, Form("PID: %d", pidx));
240 //Int_t n = track->GetNumberOfClusters();
241 // where are this tracklets ???
242 //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n);
243 if(track->GetNumberOfClustersRefit()){
244 ((TProfile*)fContainer->At(pidx))->Fill(mom, 1.);
245 labelsacc[nTRD] = track->GetLabel();
253 Int_t iref = 1; Bool_t found = kFALSE;
254 while((ref = track->GetTrackRef(iref))){
255 xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX());
256 xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX());
257 if(TMath::Abs(xmed - 298.5) < .5 &&
258 TMath::Abs(xleng - 3.7) < .1){
266 // track missing first layer. Maybe interesting for SA.
270 new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track);
273 AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
276 // Find double tracks
277 Float_t threshold = 10.;
278 AliTrackReference *refMiss = NULL;
279 AliExternalTrackParam *op = NULL;
280 AliTRDtrackInfo *tt = NULL;
281 for(Int_t imiss=0; imiss<nMiss; imiss++){
282 //printf("Searching missing %d ...\n", imiss);
284 // get outer param of missed
285 tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss);
286 op = tt->GetESDinfo()->GetOuterParam();
287 Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha);
289 Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d;
291 Bool_t bFOUND = kFALSE;
292 for(Int_t iselect=0; iselect<nselect; iselect++){
293 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]);
295 // check first MC ... if available
297 for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){
298 if(!(ref = track->GetTrackRef(iref))) continue;
299 if((refMiss = tt->GetTrackRef(iref))){
300 dy = ref->LocalY() - refMiss->LocalY();
301 dz = ref->Z() - refMiss->Z();
303 // compare missOP with refTrackRef in LTC
305 op->GetYAt(x0, AliTracker::GetBz(), y0);
306 op->GetZAt(x0, AliTracker::GetBz(), z0);
307 dy = y0 - ref->LocalY();
310 d += (dy*dy + dz*dz);
312 //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs());
313 if((track->GetNTrackRefs())){
314 d /= track->GetNTrackRefs();
316 //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs());
317 bFOUND = kTRUE; break;
321 // process outer param ... always available
322 // compare missOP with OP in GTC
323 esd = track->GetESDinfo()->GetOuterParam();
326 op->GetYAt(x0, AliTracker::GetBz(), y0);
327 op->GetZAt(x0, AliTracker::GetBz(), z0);
328 x = x0*cosa - y0*sina;
329 y = x0*sina + y0*cosa;
334 d = dx*dx+dy*dy+dz*dz;
335 //printf("\td[%d] = %f op\n", selection[iselect], d);
337 //printf("\t\tFound %2d in op[%3d] : d[%f] dx[%5.2f] dy[%5.2f] dz[%5.2f]\n", imiss, selection[iselect], d, dx, dy, dz);
338 bFOUND = kTRUE; break;
343 ref = tt->GetTrackRef(0);
344 mom = ref ? ref->P(): op->P();
345 pidx = AliTRDCalPID::GetPartIndex(tt->GetPDG());
346 pidx = TMath::Max(pidx, 0);
347 ((TProfile*)fContainer->At(pidx))->Fill(mom, 0.);
348 AliDebug(2, Form(" NOT bFOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom));
352 AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
355 // check for double countings
356 Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000);
357 TMath::Sort(nTRD, labelsacc, indices);
358 if(DebugLevel() > 2){
359 for(Int_t itk = 0; itk < nTRD - 1; itk++)
360 if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]);
364 //____________________________________________________________________
365 Int_t AliTRDefficiency::GetPtBin(Float_t pt)
367 // Get logaritmic pt bin
369 Float_t pt0(0.5), dpt(0.002);
373 ipt++; pt0+=(TMath::Exp(ipt*ipt*dpt)-1.);
378 //____________________________________________________________________
379 Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig)
381 // Steer reference figures
384 AliWarning("Please provide a canvas to draw results.");
390 Bool_t bFIRST(kTRUE);
394 h = (TProfile*)fContainer->At(AliPID::kSPECIES);
395 for(Int_t is=0; is<AliPID::kSPECIES; is++){
396 h->Add((TProfile*)fContainer->At(is));
398 h->SetMarkerStyle(24);
399 h->SetMarkerColor(kBlack);
400 h->SetLineColor(kBlack);
401 h->SetTitle("TRD Efficiency integrated");
402 h->SetXTitle("p [GeV/c]");
403 h->GetXaxis()->SetMoreLogLabels();
404 h->SetYTitle("Efficiency");
405 h->GetYaxis()->CenterTitle();
410 for(Int_t is=0; is<AliPID::kSPECIES; is++){
411 if(!(h = (TProfile*)fContainer->At(is))) continue;
412 h->SetMarkerStyle(24);
415 h->SetXTitle("p [GeV/c]");
416 h->GetXaxis()->SetMoreLogLabels();
417 h->SetYTitle("Efficiency");
418 h->GetYaxis()->CenterTitle();
419 h->GetYaxis()->SetRangeUser(0.8, 1.05);
420 leg=new TLegend(.7, .2, .98, .6);
421 leg->SetHeader("Species");
422 leg->SetBorderSize(0);
423 leg->SetFillStyle(0);
424 leg->AddEntry(h, h->GetTitle(), "pl");
426 leg->AddEntry(h, h->GetTitle(), "pl");
437 //________________________________________________________
438 TObjArray* AliTRDefficiency::Histos()
444 if(fContainer) return fContainer;
446 fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE);
450 //++++++++++++++++++++++
451 // cluster to detector
452 if(!(H = (THnSparseI*)gROOT->FindObject("hEFF"))){
455 const Char_t *eTitle[mdim] = {"status", "#phi [rad]", "eta", "p_{t} [bin]", "label", "N_{trklt}", "chg*spec*rc"};
456 const Int_t eNbins[mdim] = { 5, 180, 50, fNpt, nlabel, AliTRDgeometry::kNlayer-2, 5};
457 const Double_t eMin[mdim] = { -0.5, -TMath::Pi(), -1., -0.5, -0.5, 1.5, -2.5},
458 eMax[mdim] = { 4.5, TMath::Pi(), 1., fNpt-.5, nlabel-0.5, 5.5, 2.5};
459 st = "basic efficiency;";
460 // define minimum info to be saved in non debug mode
461 Int_t ndim=DebugLevel()>=1?mdim:(HasMCdata()?5:4);
462 for(Int_t idim(0); idim<ndim; idim++){ st += eTitle[idim]; st+=";";}
463 H = new THnSparseI("hEFF", st.Data(), ndim, eNbins, eMin, eMax);
464 /* TAxis *ax(H->GetAxis(0)); const Char_t *lTRDflag[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"};
465 for(Int_t ibin(1); ibin<=ax->GetNbins(); ibin++) ax->SetBinLabel(ibin, lTRDflag[ibin-1]);*/
467 fContainer->AddAt(H, 0);
472 //____________________________________________________________________
473 Bool_t AliTRDefficiency::PostProcess()
475 // Fit, Project, Combine, Extract values from the containers filled during execution
478 AliError("ERROR: list not available");
482 AliInfo("Building array of projections ...");
483 fProj = new TObjArray(200); fProj->SetOwner(kTRUE);
485 // set pt/p segmentation. guess from data
487 if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){
488 AliError("Missing/Wrong data @ hEFF.");
491 fNpt=H->GetAxis(3)->GetNbins()+1;
492 if(!MakeMomSegmentation()) return kFALSE;
494 if(!MakeProjectionBasicEff()) return kFALSE;
498 //____________________________________________________________________
499 Bool_t AliTRDefficiency::MakeProjectionBasicEff()
501 // Make basic efficiency plots
503 if(!fContainer || !fProj){
504 AliError("Missing data container.");
508 if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){
509 AliError("Missing/Wrong data @ hEFF.");
512 Int_t ndim(H->GetNdimensions()); //Bool_t debug(ndim>Int_t(kNdimCl));
513 TAxis *aa[11], *al(NULL); memset(aa, 0, sizeof(TAxis*) * 11);
514 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
515 if(H->GetNdimensions() > 4) al = H->GetAxis(4);
518 // define rebinning strategy
519 //const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
520 AliTRDrecoProjection hp[15]; TObjArray php(15);
521 const Char_t *stat[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"};
522 const Char_t *lab[] = {"MC-Bad", "MC-Good", "MC-Accept"};
524 for(Int_t ilab(0); ilab<nlab; ilab++){
525 for(Int_t istat(0); istat<5; istat++){
526 hp[ih].Build(Form("HEff%d%d", ilab, istat),
527 Form("Efficiency :: Stat[#bf{%s}] %s", stat[istat], nlab>1?lab[ilab]:""), 2, 1, 3, aa);
528 //hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
529 php.AddLast(&hp[ih++]); //np[isel]++;
532 AliInfo(Form("Build %3d 3D projections.", ih));
534 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
535 Int_t istatus, ilab(0), coord[11]; memset(coord, 0, sizeof(Int_t) * 11); Double_t v = 0.;
536 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
537 v = H->GetBinContent(ib, coord); if(v<1.) continue;
538 istatus = coord[0]-1;
540 Int_t isel = ilab*5+istatus;
542 AliError(Form("Wrong selection %d [%3d]", isel, ih));
545 if(!(pr0=(AliTRDrecoProjection*)php.At(isel))) {
546 AliError(Form("Missing projection %d", isel));
549 if(strcmp(pr0->H()->GetName(), Form("HEff%d%d", ilab, istatus))!=0){
550 AliError(Form("Projection mismatch :: request[HEff%d%d] found[%s]", ilab, istatus, pr0->H()->GetName()));
553 for(Int_t jh(0); jh<1/*np[isel]*/; jh++) ((AliTRDrecoProjection*)php.At(isel+jh))->Increment(coord, v);
556 TDirectory *cwd = gDirectory;
557 TFile::Open(Form("EffDump_%s.root", H->GetName()), "RECREATE");
558 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
559 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
560 if(!pr0->H()) continue;
561 TH3 *h3=(TH3*)pr0->H()->Clone();
570 if(!hp[ih].H()) continue;
571 for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(Projection2D(hp[ih].H(), ipt), jh++);
574 /* AliTRDrecoProjection prLab; TH2 *hLab[3] = {0}; TH1 *hpLab[3] = {0};
575 for(ilab=0; ilab<nlab; ilab++){
576 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 3)))) continue;
578 prLab.SetNameTitle(Form("HEffLb%d", ilab), "Sum over status");
579 prLab.H()->SetNameTitle(Form("HEffLb%d", ilab), Form("Efficiency :: #bf{%s} Propagated Tracks", lab[ilab]));
580 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 4)))) continue;
582 prLab.Projection2D(1, 10, -1, kFALSE);
583 if((hLab[ilab] = (TH2*)gDirectory->Get(Form("%sEn", prLab.H()->GetName())))) fProj->AddAt(hLab[ilab], jh++);
584 if((hpLab[ilab] = prLab.H()->Project3D("z"))) fProj->AddAt(hpLab[ilab], jh++);
587 for(Int_t istat(0); istat<5; istat++) {
588 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, istat)))) {
589 // sum over MC labels if available
590 for(ilab=1; ilab<nlab; ilab++){
591 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, istat)))) continue;
594 pr0->H()->SetNameTitle(Form("HEff%d", istat), Form("Efficiency :: Stat[#bf{%s}]", stat[istat]));
595 for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(Projection2D(pr0->H(), ipt), jh++);
597 if(istat>1 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff01"))) (*pr1)+=(*pr0);
598 if(istat>2 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff02"))) (*pr1)+=(*pr0);
599 if(istat>3 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff03"))) (*pr1)+=(*pr0);
603 const char suffix[] = {'A', 'T', 'P'};
604 const char *sname[] = {"All", "Trk", "Prp"};
605 for(Int_t istat(0); istat<3; istat++){
606 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, istat+1)))) continue;
607 pr0->H()->SetNameTitle(Form("HEff%c", suffix[istat]), Form("Efficiency :: %s Tracks", sname[istat]));
608 for(Int_t ipt(-1); ipt<=fNpt; ipt++) fProj->AddAt(Projection2D(pr0->H(), ipt), jh++);
612 TH2F *h2T(NULL), *h2P(NULL);
613 for(Int_t ipt(-1); ipt<=fNpt; ipt++){
614 if(!(h2T = (TH2F*)fProj->FindObject(Form("HEffT%d_2D", ipt)))) continue;
615 if(!(h2P = (TH2F*)fProj->FindObject(Form("HEffP%d_2D", ipt)))) continue;
617 PutTrendValue(ipt<0?"pt":(Form("pt%d", ipt)), GetMeanStat(h2P, 0.01, ">"));
619 /* // process MC label
621 for(ilab=0; ilab<nlab; ilab++){
622 if(!hLab[ilab]) continue;
624 TH2 *hEff1 = (TH2*)hLab[ilab]->Clone(Form("%sN", hLab[ilab]->GetName()));
625 for(Int_t ix(1); ix<=hLab[ilab]->GetNbinsX(); ix++){
626 for(Int_t iy(1); iy<=hLab[ilab]->GetNbinsY(); iy++){
627 if(hLab[ilab]->GetBinContent(ix, iy)<5) hEff1->SetBinContent(ix, iy, 0.);
629 hEff1->Divide(hEff[2]);
630 fProj->AddAt(hEff1, jh++);
634 for(ilab=0; ilab<nlab; ilab++){
635 if(!hpLab[ilab]) continue;
636 TH1 *hpEff1 = (TH1*)hpLab[ilab]->Clone(Form("%sN", hpLab[ilab]->GetName()));
637 hpEff1->Divide(hpEff[2]);
638 fProj->AddAt(hpEff1, jh++);
641 AliInfo(Form("Done %3d 2D projections.", jh));
645 //____________________________________________________________________
646 void AliTRDefficiency::MakeSummary()
648 // Build summary plots
650 AliError("Missing results");
653 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
655 gStyle->SetPalette(1);
657 Int_t nbins(DebugLevel()==0?3:20);
658 //calculate true pt bin
659 Float_t ptBins[23]; ptBins[0] = 0.;
660 if(nbins==3){ // significant bins
666 } else if(nbins==20){ // logarithmic bins
669 for(Int_t ipt(1); ipt<21; ipt++) ptBins[ipt+1] = ptBins[ipt]+(TMath::Exp(ipt*ipt*dpt)-1.);
672 AliError(Form("Unknown no.[%d] of bins in the p_t spectrum", nbins));
676 const Char_t cid[]={'T','P'};
677 cOut = new TCanvas(Form("%s_Eff", GetName()), "TRD Efficiency", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
678 // tracking eff :: eta-phi dependence
679 for(Int_t it(0); it<2; it++){
680 if(!(h2 = (TH2*)fProj->FindObject(Form("HEff%cEnN", cid[it])))) {
681 AliError(Form("Missing \"HEff%cEnN\".", cid[it]));
684 h2->SetContour(10); h2->Scale(1.e2); SetRangeZ(h2, 80, 100, 5);
685 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
686 p=cOut->cd(it+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
690 if(!(h2 = (TH2*)fProj->FindObject("HEff0En"))) {
691 AliError("Missing \"HEff0En\".");
694 p=cOut->cd(3);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
695 h2->Draw("colz"); MakeDetectorPlot();
696 // tracking eff :: pt dependence
698 if(!(h[0] = (TH1*)fProj->FindObject("HEffP_zN"))){
699 AliError("Missing \"HEffP_zN\".");
702 if(!(h[1] = (TH1*)fProj->FindObject("HEffT_zN"))){
703 AliError("Missing \"HEffT_zN\".");
707 Color_t color[] = {kGreen, kBlue, kRed};
708 for(Int_t il=0;il<3;il++){
709 h1[il]=new TH1F(Form("h1Eff%d", il), "", nbins+2, ptBins);
710 h1[il]->SetFillColor(color[il]);
711 h1[il]->SetFillStyle(il==2?3002:1001);
712 h1[il]->SetLineColor(color[il]);
713 h1[il]->SetMarkerStyle(4);
714 h1[il]->SetMarkerColor(color[il]);
715 h1[il]->SetLineWidth(1);
717 for(Int_t ip(0);ip<=(nbins+1);ip++){
718 h1[0]->SetBinContent(ip+1, 1.e2*h[0]->GetBinContent(ip)); // propagated
719 h1[1]->SetBinContent(ip+1, 1.e2*(h[1]->GetBinContent(ip) - h[0]->GetBinContent(ip))); // stopped
720 h1[2]->SetBinContent(ip+1, 1.e2*(1 - h[1]->GetBinContent(ip))); // missed
722 const Char_t *labEff[] = {"Propagated", "Stopped", "Missed"};
723 THStack *hs = new THStack("hEff","Tracking Efficiency;p_{t} [GeV/c];Efficiency [%]");
724 TLegend *leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
725 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
726 for(Int_t ic(0); ic<3;ic++){ hs->Add(h1[ic]);leg->AddEntry(h1[ic], labEff[ic], "f");}
727 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932);p->SetLogx();
728 hs->Draw(); //hs->Draw("same nostack,e1p");
730 hs->GetXaxis()->SetRangeUser(0.6,10.);
731 hs->GetXaxis()->SetMoreLogLabels();
732 hs->GetXaxis()->SetTitleOffset(1.2);
733 hs->GetYaxis()->SetNdivisions(513);
735 hs->GetYaxis()->CenterTitle();
736 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
738 cOut = new TCanvas(Form("%s_MC", GetName()), "TRD Label", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
739 for(Int_t ipad(0); ipad<3; ipad++){
740 p=cOut->cd(ipad+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
741 if(!(h2 = (TH2*)fProj->FindObject(Form("HEffLb%dEnN", ipad)))) continue;
743 h2->Scale(1.e2); SetRangeZ(h2, ipad==1?50:0., ipad==1?90.:50., ipad==1?0.01:0.01);
744 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
748 color[0] = kRed; color[1] = kGreen; color[2] = kBlue;
749 for(Int_t il=0;il<3;il++){
750 if(!(h[il] = (TH1D*)fProj->FindObject(Form("HEffLb%d_zN", il)))) continue;
751 h1[il]=new TH1F(Form("h1Lab%d", il), "", nbins+2, ptBins);
752 for(Int_t ip(0);ip<=(nbins+1);ip++) h1[il]->SetBinContent(ip+1, 1.e2*h[il]->GetBinContent(ip));
753 h1[il]->SetFillColor(color[il]);
754 h1[il]->SetFillStyle(il==2?3002:1001);
755 h1[il]->SetLineColor(color[il]);
756 h1[il]->SetLineWidth(1);
758 const Char_t *labMC[] = {"TRD != ESD [bad]", "TRD == ESD [good]", "TRD == -ESD [accept]"};
759 leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
760 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
761 hs = new THStack("hLab","TRD Label;p_{t} [GeV/c];Efficiency [%]");
762 hs->Add(h1[1]);leg->AddEntry(h1[1], labMC[1], "f"); // good
763 hs->Add(h1[2]);leg->AddEntry(h1[2], labMC[2], "f"); // accept
764 hs->Add(h1[0]);leg->AddEntry(h1[0], labMC[0], "f"); // bad
765 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932); p->SetLogx();
766 hs->Draw(/*"nostack,e1p"*/); leg->Draw();
767 cOut->Modified();cOut->Update();
768 hs->GetXaxis()->SetRangeUser(0.6,10.);
769 hs->GetXaxis()->SetMoreLogLabels();
770 hs->GetXaxis()->SetTitleOffset(1.2);
771 hs->GetYaxis()->SetNdivisions(513);
773 hs->GetYaxis()->CenterTitle();
774 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
778 //____________________________________________________________________
779 TH2* AliTRDefficiency::Projection2D(TH3 *h3, Int_t ipt)
781 TAxis *ax(h3->GetXaxis()), *ay(h3->GetYaxis());
784 h2 =(TH2F*)h3->Project3D("yx");
785 h2->SetNameTitle(Form("%s%d_2D", h3->GetName(), ipt), h3->GetTitle());
787 h2 = new TH2F(Form("%s%d_2D", h3->GetName(), ipt),
788 Form("%s | #it{%4.2f<=p_{t}[GeV/c]<%4.2f};%s;%s;Entries", h3->GetTitle(),
789 ipt?fgPt[ipt-1]:0., ipt==fNpt?9.99:fgPt[ipt], ax->GetTitle(), ay->GetTitle()),
790 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
791 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
792 for(Int_t ix(1); ix<=ax->GetNbins(); ix++){
793 for(Int_t iy(1); iy<=ay->GetNbins(); iy++){
794 h2->SetBinContent(ix, iy, h3->GetBinContent(ix, iy, ipt));