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