]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TRD/AliTRDefficiency.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
expand trending functionality by generic PutTrendValue/GetTrendValue
[u/mrichter/AliRoot.git] / PWGPP / TRD / AliTRDefficiency.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/* $Id: AliTRDefficiency.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Reconstruction QA //
21// //
22// Authors: //
23// Markus Fasel <M.Fasel@gsi.de> //
24// //
25////////////////////////////////////////////////////////////////////////////
26
eb05d549 27#include <TROOT.h>
668a0654 28#include <TFile.h>
eb05d549 29#include <TStyle.h>
1ee39b3a 30#include <TClonesArray.h>
31#include <TObjArray.h>
32#include <TProfile.h>
e2e3cec2 33#include <TPad.h>
eb05d549 34#include <TCanvas.h>
a310e49b 35#include <TLegend.h>
eb05d549 36#include <THnSparse.h>
37#include <TH2.h>
38#include <TH3.h>
39#include <THStack.h>
1ee39b3a 40#include "TTreeStream.h"
41
42#include "AliPID.h"
43#include "AliESDtrack.h"
44#include "AliTrackReference.h"
45#include "AliExternalTrackParam.h"
46#include "AliTracker.h"
47#include "AliAnalysisManager.h"
48
eb05d549 49#include "AliTRDgeometry.h"
50#include "AliTRDtrackV1.h"
1ee39b3a 51#include "Cal/AliTRDCalPID.h"
52#include "AliTRDefficiency.h"
53#include "info/AliTRDtrackInfo.h"
54
55ClassImp(AliTRDefficiency)
56
57//____________________________________________________________________
58AliTRDefficiency::AliTRDefficiency()
f2e89a4c 59 :AliTRDrecoTask()
db99a57a 60 ,fMissed(NULL)
eb05d549 61 ,fProj(NULL)
f8f46e4d 62{
63 //
64 // Default constructor
65 //
f2e89a4c 66 SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker");
f8f46e4d 67}
68
f2e89a4c 69//____________________________________________________________________
f8f46e4d 70AliTRDefficiency::AliTRDefficiency(char* name)
71 :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker")
db99a57a 72 ,fMissed(NULL)
eb05d549 73 ,fProj(NULL)
1ee39b3a 74{
75 //
76 // Default constructor
77 //
78}
79
80//____________________________________________________________________
81AliTRDefficiency::~AliTRDefficiency()
82{
83 // Destructor
84 if(fMissed){
85 fMissed->Delete();
86 delete fMissed;
87 }
88}
89
eb05d549 90// //____________________________________________________________________
91// void AliTRDefficiency::UserCreateOutputObjects()
92// {
93// //
94// // Create output objects
95// //
96//
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.};
99//
100// TH1 *h = NULL;
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);
107// }
108// fContainer->Add(h = new TProfile("h", "", nbins, xbins));
109// h->SetMarkerStyle(7);
110// PostData(1, fContainer);
111// }
112
1ee39b3a 113//____________________________________________________________________
eb05d549 114TH1* AliTRDefficiency::PlotBasicEff(const AliTRDtrackV1 *track)
1ee39b3a 115{
eb05d549 116// plot TRD efficiency based on ESD info
1ee39b3a 117
eb05d549 118 if(!fkESD){
119 AliDebug(4, "No ESD info.");
120 return NULL;
121 }
1ee39b3a 122
eb05d549 123 THnSparse *H(NULL);
124 if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hEFF")))){
125 AliWarning("No output container defined.");
126 return NULL;
1ee39b3a 127 }
eb05d549 128 if(track) fkTrack = track;
129
668a0654 130 Double_t val[7]; memset(val, 0, 7*sizeof(Double_t));
eb05d549 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] = fkESD->Phi();
136 val[2] = fkESD->Eta();
137 val[3] = DebugLevel()>=1?GetPtBin(fkESD->Pt()):GetPtBinSignificant(fkESD->Pt());
138 val[4] = 0.;
139 if(fkMC){
140 if(fkMC->GetLabel() == fkMC->GetTRDlabel()) val[4] = 0.;
141 else if(fkMC->GetLabel() == -fkMC->GetTRDlabel()) val[4] = 1.;
142 else val[4] = -1.;
143 }
668a0654 144 val[6] = 0;//fkTrack->GetNumberOfTracklets();
145 // down scale PID resolution
146 Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2;
147 val[7] = spc;
eb05d549 148 H->Fill(val);
149 return NULL;
150}
151
152// //____________________________________________________________________
153// TH1* AliTRDefficiency::PlotMC(const AliTRDtrackV1 *track)
154// {
155// // plot TRD efficiency based on MC info
156//
157// if(!HasMC()) return NULL;
158// if(!fkESD){
159// AliDebug(4, "No ESD info.");
160// return NULL;
161// }
162// if(!fkMC){
163// AliDebug(4, "No MC info.");
164// return NULL;
165// }
166//
167// THnSparse *H(NULL);
168// if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hMC")))){
169// AliWarning("No output container defined.");
170// return NULL;
171// }
172// if(track) fkTrack = track;
173// Double_t val[11]; memset(val, 0, 11*sizeof(Double_t));
174// ULong_t status(fkESD->GetStatus());
175// val[0] =((status&AliESDtrack::kTRDin)?1:0) +
176// ((status&AliESDtrack::kTRDStop)?1:0) +
177// ((status&AliESDtrack::kTRDout)?2:0);
178// val[1] = fkESD->Phi();
179// val[2] = fkESD->Eta();
180// val[3] = DebugLevel()>=1?GetPtBin(fkESD->Pt()):GetPtBinSignificant(fkESD->Pt());
181// if(fkTrack){ // read track status in debug mode with friends
182// val[4] = fkTrack->GetStatusTRD(-1);
183// for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++) val[5+ily]=fkTrack->GetStatusTRD(ily);
184// }
185// H->Fill(val);
186//
187// }
1ee39b3a 188
189//____________________________________________________________________
eb05d549 190void AliTRDefficiency::LocalUserExec(Option_t *)
1ee39b3a 191{
192 //
eb05d549 193 // Do it obsolete
1ee39b3a 194 //
195
eb05d549 196 Int_t labelsacc[10000];
1ee39b3a 197 memset(labelsacc, 0, sizeof(Int_t) * 10000);
198
5935a6da 199 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
200 if(!fTracks) return;
201 if(!fTracks->GetEntriesFast()) return;
202 else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));
1ee39b3a 203 if(!fMissed){
204 fMissed = new TClonesArray("AliTRDtrackInfo", 10);
205 fMissed->SetOwner();
206 }
207
208 Float_t mom;
209 Int_t selection[10000], nselect = 0;
210 ULong_t status; Int_t pidx;
211 Int_t nTRD = 0, nTPC = 0, nMiss = 0;
db99a57a 212 AliTRDtrackInfo *track = NULL;
213 AliTrackReference *ref = NULL;
214 AliExternalTrackParam *esd = NULL;
1ee39b3a 215 for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
216 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
217
218 if(!track->HasESDtrack()) continue;
219 status = track->GetStatus();
220
221 // missing TPC propagation - interesting for SA
222 if(!(status&AliESDtrack::kTPCout)) continue;
223
224 // missing MC info.
225 if(HasMCdata() && track->GetNTrackRefs() <= 1) continue;
226
227 nTPC++;
228 selection[nselect++]=itrk;
229 ref = track->GetTrackRef(0);
230 esd = track->GetESDinfo()->GetOuterParam();
231 mom = ref ? ref->P(): esd->P();
232 pidx = AliTRDCalPID::GetPartIndex(track->GetPDG());
233 pidx = TMath::Max(pidx, 0);
8f7b1226 234 AliDebug(4, Form("PID: %d", pidx));
1ee39b3a 235
236 //Int_t n = track->GetNumberOfClusters();
237 // where are this tracklets ???
238 //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n);
239 if(track->GetNumberOfClustersRefit()){
240 ((TProfile*)fContainer->At(pidx))->Fill(mom, 1.);
241 labelsacc[nTRD] = track->GetLabel();
242 nTRD++;
243 continue;
244 }
245
246
247
248 Float_t xmed, xleng;
249 Int_t iref = 1; Bool_t found = kFALSE;
250 while((ref = track->GetTrackRef(iref))){
251 xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX());
252 xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX());
253 if(TMath::Abs(xmed - 298.5) < .5 &&
254 TMath::Abs(xleng - 3.7) < .1){
255 found = kTRUE;
256 break;
257 }
258 iref++;
259 }
260 if(!found){
261 nTPC--;
262 // track missing first layer. Maybe interesting for SA.
263 continue;
264 }
265 nselect--;
266 new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track);
267 nMiss++;
268 }
269 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()));
270
271
272 // Find double tracks
273 Float_t threshold = 10.;
db99a57a 274 AliTrackReference *refMiss = NULL;
275 AliExternalTrackParam *op = NULL;
276 AliTRDtrackInfo *tt = NULL;
1ee39b3a 277 for(Int_t imiss=0; imiss<nMiss; imiss++){
278 //printf("Searching missing %d ...\n", imiss);
279
280 // get outer param of missed
281 tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss);
282 op = tt->GetESDinfo()->GetOuterParam();
283 Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha);
284
285 Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d;
286
287 Bool_t bFOUND = kFALSE;
288 for(Int_t iselect=0; iselect<nselect; iselect++){
289 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]);
290
291 // check first MC ... if available
292 d = 0;
293 for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){
294 if(!(ref = track->GetTrackRef(iref))) continue;
295 if((refMiss = tt->GetTrackRef(iref))){
296 dy = ref->LocalY() - refMiss->LocalY();
297 dz = ref->Z() - refMiss->Z();
298 } else {
299 // compare missOP with refTrackRef in LTC
300 x0 = ref->LocalX();
301 op->GetYAt(x0, AliTracker::GetBz(), y0);
302 op->GetZAt(x0, AliTracker::GetBz(), z0);
303 dy = y0 - ref->LocalY();
304 dz = z0 - ref->Z();
305 }
306 d += (dy*dy + dz*dz);
307 }
308 //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs());
309 if((track->GetNTrackRefs())){
310 d /= track->GetNTrackRefs();
311 if(d < threshold){
312 //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs());
313 bFOUND = kTRUE; break;
314 }
315 }
316
317 // process outer param ... always available
318 // compare missOP with OP in GTC
319 esd = track->GetESDinfo()->GetOuterParam();
320 esd->GetXYZ(xyz);
321 x0 = esd->GetX();
322 op->GetYAt(x0, AliTracker::GetBz(), y0);
323 op->GetZAt(x0, AliTracker::GetBz(), z0);
324 x = x0*cosa - y0*sina;
325 y = x0*sina + y0*cosa;
326 z = z0;
327 dx=xyz[0]-x;
328 dy=xyz[1]-y;
329 dz=xyz[2]-z;
330 d = dx*dx+dy*dy+dz*dz;
331 //printf("\td[%d] = %f op\n", selection[iselect], d);
332 if(d < threshold){
333 //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);
334 bFOUND = kTRUE; break;
335 }
336 }
337 if(bFOUND) nTPC--;
338 else{
339 ref = tt->GetTrackRef(0);
340 mom = ref ? ref->P(): op->P();
341 pidx = AliTRDCalPID::GetPartIndex(tt->GetPDG());
342 pidx = TMath::Max(pidx, 0);
343 ((TProfile*)fContainer->At(pidx))->Fill(mom, 0.);
344 AliDebug(2, Form(" NOT bFOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom));
345 }
346 }
347
348 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()));
349
350 //fMissed->Delete();
351 // check for double countings
352 Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000);
353 TMath::Sort(nTRD, labelsacc, indices);
354 if(DebugLevel() > 2){
355 for(Int_t itk = 0; itk < nTRD - 1; itk++)
356 if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]);
357 }
1ee39b3a 358}
359
eb05d549 360//____________________________________________________________________
361Int_t AliTRDefficiency::GetPtBin(Float_t pt)
362{
363// Get logaritmic pt bin
364
365 Float_t pt0(0.5), dpt(0.002);
366 Int_t ipt(0);
367 while(ipt<30){
368 if(pt<pt0) break;
369 ipt++; pt0+=(TMath::Exp(ipt*ipt*dpt)-1.);
370 }
371 return ipt-1;
372}
373
1ee39b3a 374//____________________________________________________________________
375Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig)
376{
377// Steer reference figures
378
e2e3cec2 379 if(!gPad){
380 AliWarning("Please provide a canvas to draw results.");
381 return kFALSE;
382 }
383 gPad->SetLogx();
384
a310e49b 385 TLegend *leg(NULL);
386 Bool_t bFIRST(kTRUE);
387 TProfile *h(NULL);
1ee39b3a 388 switch(ifig){
389 case 0:
390 h = (TProfile*)fContainer->At(AliPID::kSPECIES);
391 for(Int_t is=0; is<AliPID::kSPECIES; is++){
392 h->Add((TProfile*)fContainer->At(is));
393 }
a310e49b 394 h->SetMarkerStyle(24);
1ee39b3a 395 h->SetMarkerColor(kBlack);
396 h->SetLineColor(kBlack);
e2e3cec2 397 h->SetTitle("TRD Efficiency integrated");
a310e49b 398 h->SetXTitle("p [GeV/c]");
e2e3cec2 399 h->GetXaxis()->SetMoreLogLabels();
a310e49b 400 h->SetYTitle("Efficiency");
401 h->GetYaxis()->CenterTitle();
1ee39b3a 402 h->Draw("e1");
403 break;
404 case 1:
405 bFIRST = kTRUE;
406 for(Int_t is=0; is<AliPID::kSPECIES; is++){
407 if(!(h = (TProfile*)fContainer->At(is))) continue;
a310e49b 408 h->SetMarkerStyle(24);
1ee39b3a 409 if(bFIRST){
410 h->Draw("e1");
a310e49b 411 h->SetXTitle("p [GeV/c]");
e2e3cec2 412 h->GetXaxis()->SetMoreLogLabels();
a310e49b 413 h->SetYTitle("Efficiency");
414 h->GetYaxis()->CenterTitle();
415 h->GetYaxis()->SetRangeUser(0.8, 1.05);
416 leg=new TLegend(.7, .2, .98, .6);
417 leg->SetHeader("Species");
418 leg->SetBorderSize(0);
419 leg->SetFillStyle(0);
420 leg->AddEntry(h, h->GetTitle(), "pl");
421 } else {
422 leg->AddEntry(h, h->GetTitle(), "pl");
423 h->Draw("same e1");
424 }
1ee39b3a 425 bFIRST = kFALSE;
426 }
7fe4e88b 427 if(leg) leg->Draw();
1ee39b3a 428 break;
429 }
430 return kTRUE;
431}
432
eb05d549 433//________________________________________________________
434TObjArray* AliTRDefficiency::Histos()
435{
436 //
437 // Define histograms
438 //
439
440 if(fContainer) return fContainer;
441
442 fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE);
443 THnSparse *H(NULL);
444 TString st;
445
446 //++++++++++++++++++++++
447 // cluster to detector
448 if(!(H = (THnSparseI*)gROOT->FindObject("hEFF"))){
668a0654 449 const Int_t mdim(7);
eb05d549 450 Int_t npt=DebugLevel()>=1?20:3;
451 Int_t nlabel(1);
668a0654 452 const Char_t *eTitle[mdim] = {"label", "#phi [rad]", "eta", "p_{t} [bin]", "MClabel", "n_trklt", "chg*spec*rc"};
453 const Int_t eNbins[mdim] = { 5, 180, 50, npt, nlabel, AliTRDgeometry::kNlayer-2, 5};
454 const Double_t eMin[mdim] = { -0.5, -TMath::Pi(), -1., -0.5, -0.5, 0.5, -2.5},
455 eMax[mdim] = { 4.5, TMath::Pi(), 1., npt-.5, nlabel-0.5, 4.5, 2.5};
eb05d549 456 st = "basic efficiency;";
457 // define minimum info to be saved in non debug mode
458 Int_t ndim=DebugLevel()>=1?mdim:(HasMCdata()?5:4);
459 for(Int_t idim(0); idim<ndim; idim++){ st += eTitle[idim]; st+=";";}
460 H = new THnSparseI("hEFF", st.Data(), ndim, eNbins, eMin, eMax);
461/* TAxis *ax(H->GetAxis(0)); const Char_t *lTRDflag[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"};
462 for(Int_t ibin(1); ibin<=ax->GetNbins(); ibin++) ax->SetBinLabel(ibin, lTRDflag[ibin-1]);*/
463 } else H->Reset();
464 fContainer->AddAt(H, 0);
465
466 return fContainer;
467}
1ee39b3a 468
469//____________________________________________________________________
470Bool_t AliTRDefficiency::PostProcess()
471{
eb05d549 472// Fit, Project, Combine, Extract values from the containers filled during execution
473
474 if (!fContainer) {
475 AliError("ERROR: list not available");
476 return kFALSE;
477 }
478 if(!fProj){
479 AliInfo("Building array of projections ...");
480 fProj = new TObjArray(50); fProj->SetOwner(kTRUE);
481 }
482 if(!MakeProjectionBasicEff()) return kFALSE;
483 return kTRUE;
484}
485
486//____________________________________________________________________
487Bool_t AliTRDefficiency::MakeProjectionBasicEff()
488{
489// Make basic efficiency plots
490
491 if(!fContainer || !fProj){
492 AliError("Missing data container.");
493 return kFALSE;
494 }
495 THnSparse *H(NULL);
496 if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){
497 AliError("Missing/Wrong data @ hEFF.");
498 return kFALSE;
499 }
500 Int_t ndim(H->GetNdimensions()); //Bool_t debug(ndim>Int_t(kNdimCl));
501 TAxis *aa[11], *al(NULL); memset(aa, 0, sizeof(TAxis*) * 11);
502 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
4b92200c 503 if(aa[3]->GetNbins()>3) SetDebugLevel(1);
eb05d549 504 if(H->GetNdimensions() > 4) al = H->GetAxis(4);
505 Int_t nlab=al?3:1;
506
507 // define rebinning strategy
508 //const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
509 AliTRDrecoProjection hp[15]; TObjArray php(15);
510 const Char_t *stat[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"};
511 const Char_t *lab[] = {"Bad", "Good", "Accept"};
512 Int_t ih(0);
513 for(Int_t ilab(0); ilab<nlab; ilab++){
514 for(Int_t istat(0); istat<5; istat++){
515// isel++; // new selection
516 hp[ih].Build(Form("HEff%d%d", ilab, istat),
517 Form("Efficiency :: Lab[%s] Stat[#bf{%s}]", lab[ilab], stat[istat]),
518 2, 1, 3, aa);
519 //hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
520 php.AddLast(&hp[ih++]); //np[isel]++;
521 }
522 }
523 AliInfo(Form("Build %3d 3D projections.", ih));
524
668a0654 525 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
eb05d549 526 Int_t istatus, ilab(0), coord[11]; memset(coord, 0, sizeof(Int_t) * 11); Double_t v = 0.;
527 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
528 v = H->GetBinContent(ib, coord); if(v<1.) continue;
529 istatus = coord[0]-1;
530 if(al) ilab = coord[4];
531 Int_t isel = ilab*5+istatus;
668a0654 532 if(isel>=ih){
533 AliError(Form("Wrong selection %d [%3d]", isel, ih));
534 return kFALSE;
535 }
536 if(!(pr0=(AliTRDrecoProjection*)php.At(isel))) {
537 AliError(Form("Missing projection %d", isel));
538 return kFALSE;
539 }
540 if(strcmp(pr0->H()->GetName(), Form("HEff%d%d", ilab, istatus))!=0){
541 AliError(Form("Projection mismatch :: request[HEff%d%d] found[%s]", ilab, istatus, pr0->H()->GetName()));
542 return kFALSE;
543 }
eb05d549 544 for(Int_t jh(0); jh<1/*np[isel]*/; jh++) ((AliTRDrecoProjection*)php.At(isel+jh))->Increment(coord, v);
545 }
668a0654 546 if(HasDump3D()){
547 TDirectory *cwd = gDirectory;
548 TFile::Open(Form("EffDump_%s.root", H->GetName()), "RECREATE");
549 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
550 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
551 if(!pr0->H()) continue;
552 TH3 *h3=(TH3*)pr0->H()->Clone();
553 h3->Write();
554 }
555 gFile->Close();
556 cwd->cd();
557 }
558
eb05d549 559 TH2 *h2(NULL); Int_t jh(0);
560 for(; ih--; ){
561 if(!hp[ih].H()) continue;
562 hp[ih].Projection2D(1, 10, -1, kFALSE);
563 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) fProj->AddAt(h2, jh++);
564 }
565
eb05d549 566 AliTRDrecoProjection prLab; TH2 *hLab[3] = {0}; TH1 *hpLab[3] = {0};
567 for(ilab=0; ilab<nlab; ilab++){
568 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 3)))) continue;
569 prLab=(*pr0);
570 prLab.SetNameTitle(Form("HEffLb%d", ilab), "Sum over status");
571 prLab.H()->SetNameTitle(Form("HEffLb%d", ilab), Form("Efficiency :: #bf{%s} Propagated Tracks", lab[ilab]));
572 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, 4)))) continue;
573 prLab+=(*pr1);
4b92200c 574 prLab.Projection2D(1, 10, -1, kFALSE);
eb05d549 575 if((hLab[ilab] = (TH2*)gDirectory->Get(Form("%sEn", prLab.H()->GetName())))) fProj->AddAt(hLab[ilab], jh++);
576 if((hpLab[ilab] = prLab.H()->Project3D("z"))) fProj->AddAt(hpLab[ilab], jh++);
577 }
578
579 for(Int_t istat(0); istat<5; istat++) {
580 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, istat)))) {
581 for(ilab=1; ilab<nlab; ilab++){
582 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", ilab, istat)))) continue;
583 (*pr0)+=(*pr1);
584 }
585 pr0->H()->SetNameTitle(Form("HEff%d", istat), Form("Efficiency :: Stat[#bf{%s}]", stat[istat]));
4b92200c 586 pr0->Projection2D(1, 10, -1, kFALSE);
eb05d549 587 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) fProj->AddAt(h2, jh++);
588
589 if(istat>1 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff01"))) (*pr1)+=(*pr0);
590 if(istat>2 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff02"))) (*pr1)+=(*pr0);
591 if(istat>3 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff03"))) (*pr1)+=(*pr0);
592 }
593 }
594 // All tracks
595 TH2 *hEff[3] = {0};TH1 *hpEff[3] = {0};
596 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, 1)))) {
597 pr0->H()->SetNameTitle("HEff", "Efficiency :: All Tracks");
4b92200c 598 pr0->Projection2D(1, 10, -1, kFALSE);
eb05d549 599 hEff[0] = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName()));
600 hpEff[0]= pr0->H()->Project3D("z");
601 }
602 // Tracked tracks
603 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, 2)))) {
668a0654 604 pr0->H()->SetNameTitle("HEffT", "Efficiency :: Tracked Tracks");
4b92200c 605 pr0->Projection2D(1, 10, -1, kFALSE);
eb05d549 606 hEff[1] = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName()));
668a0654 607 // remove fakes
608 for(Int_t ix(1); ix<=hEff[1]->GetNbinsX(); ix++){
609 for(Int_t iy(1); iy<=hEff[1]->GetNbinsY(); iy++){
610 if(hEff[1]->GetBinContent(ix, iy)<5) hEff[1]->SetBinContent(ix, iy, 0.);
611 }}
eb05d549 612 hpEff[1]= pr0->H()->Project3D("z");
613 }
614 // Propagated tracks
615 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, 3)))) {
668a0654 616 pr0->H()->SetNameTitle("HEffP", "Efficiency :: Propagated Tracks");
4b92200c 617 pr0->Projection2D(1, 10, -1, kFALSE);
eb05d549 618 hEff[2] = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName()));
668a0654 619 // remove fakes
620 for(Int_t ix(1); ix<=hEff[2]->GetNbinsX(); ix++){
621 for(Int_t iy(1); iy<=hEff[2]->GetNbinsY(); iy++){
622 if(hEff[2]->GetBinContent(ix, iy)<5) hEff[2]->SetBinContent(ix, iy, 0.);
623 }}
eb05d549 624 hpEff[2]= pr0->H()->Project3D("z");
625 }
626 if(hEff[0]){
627 if(hEff[1]){
668a0654 628 TH2 *hEff1 = (TH2*)hEff[1]->Clone(Form("%sN", hEff[1]->GetName()));
629 hEff1->Divide(hEff[0]);
630 fProj->AddAt(hEff1, jh++);
eb05d549 631 }
632 if(hEff[2]){
668a0654 633 TH2 *hEff1 = (TH2*)hEff[2]->Clone(Form("%sN", hEff[2]->GetName()));
eb05d549 634 hEff1->Divide(hEff[0]);
635 fProj->AddAt(hEff1, jh++);
636 }
637 }
638 if(hpEff[0]){
639 if(hpEff[1]){
668a0654 640 TH1 *hpEff1 = (TH1*)hpEff[1]->Clone(Form("%sN", hpEff[1]->GetName()));
641 hpEff1->Divide(hpEff[0]);
642 fProj->AddAt(hpEff1, jh++);
eb05d549 643 }
644 if(hEff[2]){
668a0654 645 TH1 *hpEff1 = (TH1*)hpEff[2]->Clone(Form("%sN", hpEff[2]->GetName()));
eb05d549 646 hpEff1->Divide(hpEff[0]);
647 fProj->AddAt(hpEff1, jh++);
648 }
649 }
650 // process MC label
651 if(hEff[2]){
652 for(ilab=0; ilab<nlab; ilab++){
653 if(!hLab[ilab]) continue;
668a0654 654 // remove fakes
655 TH2 *hEff1 = (TH2*)hLab[ilab]->Clone(Form("%sN", hLab[ilab]->GetName()));
656 for(Int_t ix(1); ix<=hLab[ilab]->GetNbinsX(); ix++){
657 for(Int_t iy(1); iy<=hLab[ilab]->GetNbinsY(); iy++){
658 if(hLab[ilab]->GetBinContent(ix, iy)<5) hEff1->SetBinContent(ix, iy, 0.);
659 }}
660 hEff1->Divide(hEff[2]);
661 fProj->AddAt(hEff1, jh++);
eb05d549 662 }
663 }
664 if(hpEff[2]){
665 for(ilab=0; ilab<nlab; ilab++){
666 if(!hpLab[ilab]) continue;
668a0654 667 TH1 *hpEff1 = (TH1*)hpLab[ilab]->Clone(Form("%sN", hpLab[ilab]->GetName()));
668 hpEff1->Divide(hpEff[2]);
669 fProj->AddAt(hpEff1, jh++);
eb05d549 670 }
671 }
672 AliInfo(Form("Done %3d 2D projections.", jh));
1ee39b3a 673 return kTRUE;
674}
eb05d549 675
676//____________________________________________________________________
677void AliTRDefficiency::MakeSummary()
678{
679// Build summary plots
680 if(!fProj){
681 AliError("Missing results");
682 return;
683 }
684 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
685 TH2 *h2(NULL);
686 gStyle->SetPalette(1);
687
4b92200c 688 Int_t nbins(DebugLevel()==0?3:20);
eb05d549 689 //calculate true pt bin
690 Float_t ptBins[23]; ptBins[0] = 0.;
691 if(nbins==3){ // significant bins
692 ptBins[1] = 0.5;
693 ptBins[2] = 0.8;
694 ptBins[3] = 1.5;
695 ptBins[4] = 5.;
696 ptBins[5] = 10.;
697 } else if(nbins==20){ // logarithmic bins
698 ptBins[1] = 0.5;
699 Float_t dpt(0.002);
700 for(Int_t ipt(1); ipt<21; ipt++) ptBins[ipt+1] = ptBins[ipt]+(TMath::Exp(ipt*ipt*dpt)-1.);
701 ptBins[22] = 10.;
702 } else {
703 AliError(Form("Unknown no.[%d] of bins in the p_t spectrum", nbins));
704 return;// kFALSE;
705 }
706
33056e04 707 const Char_t cid[]={'T','P'};
eb05d549 708 cOut = new TCanvas(Form("%s_Eff", GetName()), "TRD Efficiency", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
709 // tracking eff :: eta-phi dependence
710 for(Int_t it(0); it<2; it++){
668a0654 711 if(!(h2 = (TH2*)fProj->FindObject(Form("HEff%cEnN", cid[it])))) {
712 AliError(Form("Missing \"HEff%cEnN\".", cid[it]));
eb05d549 713 continue;
714 }
33056e04 715 h2->SetContour(10); h2->Scale(1.e2); SetRangeZ(h2, 80, 100, 5);
eb05d549 716 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
717 p=cOut->cd(it+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
718 h2->Draw("colz");
33056e04 719 MakeDetectorPlot();
eb05d549 720 }
721 if(!(h2 = (TH2*)fProj->FindObject("HEff0En"))) {
722 AliError("Missing \"HEff0En\".");
723 return;
724 }
725 p=cOut->cd(3);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
33056e04 726 h2->Draw("colz"); MakeDetectorPlot();
eb05d549 727 // tracking eff :: pt dependence
728 TH1 *h[2] = {0};
668a0654 729 if(!(h[0] = (TH1*)fProj->FindObject("HEffP_zN"))){
730 AliError("Missing \"HEffP_zN\".");
eb05d549 731 return;
732 }
668a0654 733 if(!(h[1] = (TH1*)fProj->FindObject("HEffT_zN"))){
734 AliError("Missing \"HEffT_zN\".");
eb05d549 735 return;
736 }
737 TH1 *h1[3] = {0};
738 Color_t color[] = {kGreen, kBlue, kRed};
739 for(Int_t il=0;il<3;il++){
740 h1[il]=new TH1F(Form("h1Eff%d", il), "", nbins+2, ptBins);
741 h1[il]->SetFillColor(color[il]);
742 h1[il]->SetFillStyle(il==2?3002:1001);
743 h1[il]->SetLineColor(color[il]);
668a0654 744 h1[il]->SetMarkerStyle(4);
745 h1[il]->SetMarkerColor(color[il]);
eb05d549 746 h1[il]->SetLineWidth(1);
747 }
748 for(Int_t ip(0);ip<=(nbins+1);ip++){
749 h1[0]->SetBinContent(ip+1, 1.e2*h[0]->GetBinContent(ip)); // propagated
750 h1[1]->SetBinContent(ip+1, 1.e2*(h[1]->GetBinContent(ip) - h[0]->GetBinContent(ip))); // stopped
751 h1[2]->SetBinContent(ip+1, 1.e2*(1 - h[1]->GetBinContent(ip))); // missed
752 }
33056e04 753 const Char_t *labEff[] = {"Propagated", "Stopped", "Missed"};
eb05d549 754 THStack *hs = new THStack("hEff","Tracking Efficiency;p_{t} [GeV/c];Efficiency [%]");
755 TLegend *leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
756 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
757 for(Int_t ic(0); ic<3;ic++){ hs->Add(h1[ic]);leg->AddEntry(h1[ic], labEff[ic], "f");}
758 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932);p->SetLogx();
668a0654 759 hs->Draw(); //hs->Draw("same nostack,e1p");
760 leg->Draw();
eb05d549 761 hs->GetXaxis()->SetRangeUser(0.6,10.);
762 hs->GetXaxis()->SetMoreLogLabels();
763 hs->GetXaxis()->SetTitleOffset(1.2);
764 hs->GetYaxis()->SetNdivisions(513);
668a0654 765 hs->SetMinimum(75.);
eb05d549 766 hs->GetYaxis()->CenterTitle();
767 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
768
769 cOut = new TCanvas(Form("%s_MC", GetName()), "TRD Label", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
770 for(Int_t ipad(0); ipad<3; ipad++){
771 p=cOut->cd(ipad+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
668a0654 772 if(!(h2 = (TH2*)fProj->FindObject(Form("HEffLb%dEnN", ipad)))) continue;
eb05d549 773 h2->SetContour(10);
668a0654 774 h2->Scale(1.e2); SetRangeZ(h2, ipad==1?50:0., ipad==1?90.:50., ipad==1?0.01:0.01);
eb05d549 775 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
776 h2->Draw("colz");
33056e04 777 MakeDetectorPlot();
eb05d549 778 }
33056e04 779 color[0] = kRed; color[1] = kGreen; color[2] = kBlue;
eb05d549 780 for(Int_t il=0;il<3;il++){
668a0654 781 if(!(h[il] = (TH1D*)fProj->FindObject(Form("HEffLb%d_zN", il)))) continue;
eb05d549 782 h1[il]=new TH1F(Form("h1Lab%d", il), "", nbins+2, ptBins);
783 for(Int_t ip(0);ip<=(nbins+1);ip++) h1[il]->SetBinContent(ip+1, 1.e2*h[il]->GetBinContent(ip));
33056e04 784 h1[il]->SetFillColor(color[il]);
eb05d549 785 h1[il]->SetFillStyle(il==2?3002:1001);
33056e04 786 h1[il]->SetLineColor(color[il]);
eb05d549 787 h1[il]->SetLineWidth(1);
788 }
33056e04 789 const Char_t *labMC[] = {"TRD != ESD [bad]", "TRD == ESD [good]", "TRD == -ESD [accept]"};
eb05d549 790 leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
791 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
792 hs = new THStack("hLab","TRD Label;p_{t} [GeV/c];Efficiency [%]");
793 hs->Add(h1[1]);leg->AddEntry(h1[1], labMC[1], "f"); // good
794 hs->Add(h1[2]);leg->AddEntry(h1[2], labMC[2], "f"); // accept
795 hs->Add(h1[0]);leg->AddEntry(h1[0], labMC[0], "f"); // bad
796 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932); p->SetLogx();
668a0654 797 hs->Draw(/*"nostack,e1p"*/); leg->Draw();
eb05d549 798 cOut->Modified();cOut->Update();
799 hs->GetXaxis()->SetRangeUser(0.6,10.);
800 hs->GetXaxis()->SetMoreLogLabels();
801 hs->GetXaxis()->SetTitleOffset(1.2);
802 hs->GetYaxis()->SetNdivisions(513);
668a0654 803 hs->SetMinimum(50.);
eb05d549 804 hs->GetYaxis()->CenterTitle();
805 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
806}