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 <TFile.h>
  29 #include <TStyle.h>
  30 #include <TClonesArray.h>
  31 #include <TObjArray.h>
  32 #include <TProfile.h>
  33 #include <TPad.h>
  34 #include <TCanvas.h>
  35 #include <TLegend.h>
  36 #include <THnSparse.h>
  37 #include <TH2.h>
  38 #include <TH3.h>
  39 #include <THStack.h>
  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
  49 #include "AliTRDgeometry.h"
  50 #include "AliTRDtrackV1.h"
  51 #include "Cal/AliTRDCalPID.h"
  52 #include "AliTRDefficiency.h"
  53 #include "info/AliTRDtrackInfo.h"
  54
  55 ClassImp(AliTRDefficiency)
  56
  57 //____________________________________________________________________
  58 AliTRDefficiency::AliTRDefficiency()
  59   :AliTRDrecoTask()
  60   ,fMissed(NULL)
  61   ,fProj(NULL)
  62 {
  63   //
  64   // Default constructor
  65   //
  66   SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker");
  67 }
  68
  69 //____________________________________________________________________
  70 AliTRDefficiency::AliTRDefficiency(char* name)
  71   :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker")
  72   ,fMissed(NULL)
  73   ,fProj(NULL)
  74 {
  75   //
  76   // Default constructor
  77   //
  78 }
  79
  80 //____________________________________________________________________
  81 AliTRDefficiency::~AliTRDefficiency()
  82 {
  83   // Destructor
  84   if(fMissed){
  85     fMissed->Delete();
  86     delete fMissed;
  87   }
  88 }
  89
  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
 113 //____________________________________________________________________
 114 TH1* AliTRDefficiency::PlotBasicEff(const AliTRDtrackV1 *track)
 115 {
 116 // plot TRD efficiency based on ESD info
 117
 118   if(!fkESD){
 119     AliDebug(4, "No ESD info.");
 120     return NULL;
 121   }
 122
 123   THnSparse *H(NULL);
 124   if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hEFF")))){
 125     AliWarning("No output container defined.");
 126     return NULL;
 127   }
 128   if(track) fkTrack = track;
 129
 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] = 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   }
 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;
 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 // }
 188
 189 //____________________________________________________________________
 190 void AliTRDefficiency::LocalUserExec(Option_t *)
 191 {
 192   //
 193   // Do it obsolete
 194   //
 195
 196   Int_t labelsacc[10000];
 197   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 198
 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()));
 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;
 212   AliTRDtrackInfo     *track = NULL;
 213   AliTrackReference     *ref = NULL;
 214   AliExternalTrackParam *esd = NULL;
 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);
 234     AliDebug(4, Form("PID: %d", pidx));
 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.;
 274   AliTrackReference *refMiss = NULL;
 275   AliExternalTrackParam *op = NULL;
 276   AliTRDtrackInfo       *tt = NULL;
 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         }
 358 }
 359
 360 //____________________________________________________________________
 361 Int_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
 374 //____________________________________________________________________
 375 Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig)
 376 {
 377 // Steer reference figures
 378
 379   if(!gPad){
 380     AliWarning("Please provide a canvas to draw results.");
 381     return kFALSE;
 382   }
 383   gPad->SetLogx();
 384
 385   TLegend *leg(NULL);
 386   Bool_t bFIRST(kTRUE);
 387   TProfile *h(NULL);
 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     }
 394     h->SetMarkerStyle(24);
 395     h->SetMarkerColor(kBlack);
 396     h->SetLineColor(kBlack);
 397     h->SetTitle("TRD Efficiency integrated");
 398     h->SetXTitle("p [GeV/c]");
 399     h->GetXaxis()->SetMoreLogLabels();
 400     h->SetYTitle("Efficiency");
 401     h->GetYaxis()->CenterTitle();
 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;
 408       h->SetMarkerStyle(24);
 409       if(bFIRST){
 410         h->Draw("e1");
 411         h->SetXTitle("p [GeV/c]");
 412         h->GetXaxis()->SetMoreLogLabels();
 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       }
 425       bFIRST = kFALSE;
 426     }
 427     if(leg) leg->Draw();
 428     break;
 429   }
 430   return kTRUE;
 431 }
 432
 433 //________________________________________________________
 434 TObjArray* 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"))){
 449     const Int_t mdim(7);
 450     Int_t npt=DebugLevel()>=1?20:3;
 451     Int_t nlabel(1);
 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};
 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 }
 468
 469 //____________________________________________________________________
 470 Bool_t AliTRDefficiency::PostProcess()
 471 {
 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 //____________________________________________________________________
 487 Bool_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);
 503   if(aa[3]->GetNbins()>3) SetDebugLevel(1);
 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
 525   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
 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;
 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     }
 544     for(Int_t jh(0); jh<1/*np[isel]*/; jh++) ((AliTRDrecoProjection*)php.At(isel+jh))->Increment(coord, v);
 545   }
 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
 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
 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);
 574     prLab.Projection2D(1, 10, -1, kFALSE);
 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]));
 586       pr0->Projection2D(1, 10, -1, kFALSE);
 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");
 598     pr0->Projection2D(1, 10, -1, kFALSE);
 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)))) {
 604     pr0->H()->SetNameTitle("HEffT", "Efficiency :: Tracked Tracks");
 605     pr0->Projection2D(1, 10, -1, kFALSE);
 606     hEff[1] = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName()));
 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     }}
 612     hpEff[1]= pr0->H()->Project3D("z");
 613   }
 614   // Propagated tracks
 615   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d", 0, 3)))) {
 616     pr0->H()->SetNameTitle("HEffP", "Efficiency :: Propagated Tracks");
 617     pr0->Projection2D(1, 10, -1, kFALSE);
 618     hEff[2] = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName()));
 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     }}
 624     hpEff[2]= pr0->H()->Project3D("z");
 625   }
 626   if(hEff[0]){
 627     if(hEff[1]){
 628       TH2 *hEff1 = (TH2*)hEff[1]->Clone(Form("%sN", hEff[1]->GetName()));
 629       hEff1->Divide(hEff[0]);
 630       fProj->AddAt(hEff1, jh++);
 631     }
 632     if(hEff[2]){
 633       TH2 *hEff1 = (TH2*)hEff[2]->Clone(Form("%sN", hEff[2]->GetName()));
 634       hEff1->Divide(hEff[0]);
 635       fProj->AddAt(hEff1, jh++);
 636     }
 637   }
 638   if(hpEff[0]){
 639     if(hpEff[1]){
 640       TH1 *hpEff1 = (TH1*)hpEff[1]->Clone(Form("%sN", hpEff[1]->GetName()));
 641       hpEff1->Divide(hpEff[0]);
 642       fProj->AddAt(hpEff1, jh++);
 643     }
 644     if(hEff[2]){
 645       TH1 *hpEff1 = (TH1*)hpEff[2]->Clone(Form("%sN", hpEff[2]->GetName()));
 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;
 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++);
 662     }
 663   }
 664   if(hpEff[2]){
 665     for(ilab=0; ilab<nlab; ilab++){
 666       if(!hpLab[ilab]) continue;
 667       TH1 *hpEff1 = (TH1*)hpLab[ilab]->Clone(Form("%sN", hpLab[ilab]->GetName()));
 668       hpEff1->Divide(hpEff[2]);
 669       fProj->AddAt(hpEff1, jh++);
 670     }
 671   }
 672   AliInfo(Form("Done %3d 2D projections.", jh));
 673   return kTRUE;
 674 }
 675
 676 //____________________________________________________________________
 677 void 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
 688   Int_t nbins(DebugLevel()==0?3:20);
 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
 707   const Char_t cid[]={'T','P'};
 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++){
 711     if(!(h2 = (TH2*)fProj->FindObject(Form("HEff%cEnN", cid[it])))) {
 712       AliError(Form("Missing \"HEff%cEnN\".", cid[it]));
 713       continue;
 714     }
 715     h2->SetContour(10); h2->Scale(1.e2); SetRangeZ(h2, 80, 100, 5);
 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");
 719     MakeDetectorPlot();
 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);
 726   h2->Draw("colz"); MakeDetectorPlot();
 727   // tracking eff :: pt dependence
 728   TH1 *h[2] = {0};
 729   if(!(h[0] = (TH1*)fProj->FindObject("HEffP_zN"))){
 730     AliError("Missing \"HEffP_zN\".");
 731     return;
 732   }
 733   if(!(h[1] = (TH1*)fProj->FindObject("HEffT_zN"))){
 734     AliError("Missing \"HEffT_zN\".");
 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]);
 744     h1[il]->SetMarkerStyle(4);
 745     h1[il]->SetMarkerColor(color[il]);
 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   }
 753   const Char_t *labEff[] = {"Propagated", "Stopped", "Missed"};
 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();
 759   hs->Draw(); //hs->Draw("same nostack,e1p");
 760   leg->Draw();
 761   hs->GetXaxis()->SetRangeUser(0.6,10.);
 762   hs->GetXaxis()->SetMoreLogLabels();
 763   hs->GetXaxis()->SetTitleOffset(1.2);
 764   hs->GetYaxis()->SetNdivisions(513);
 765   hs->SetMinimum(75.);
 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);
 772     if(!(h2 = (TH2*)fProj->FindObject(Form("HEffLb%dEnN", ipad)))) continue;
 773     h2->SetContour(10);
 774     h2->Scale(1.e2); SetRangeZ(h2, ipad==1?50:0., ipad==1?90.:50., ipad==1?0.01:0.01);
 775     h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
 776     h2->Draw("colz");
 777     MakeDetectorPlot();
 778   }
 779   color[0] = kRed; color[1] = kGreen; color[2] = kBlue;
 780   for(Int_t il=0;il<3;il++){
 781     if(!(h[il] = (TH1D*)fProj->FindObject(Form("HEffLb%d_zN", il)))) continue;
 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));
 784     h1[il]->SetFillColor(color[il]);
 785     h1[il]->SetFillStyle(il==2?3002:1001);
 786     h1[il]->SetLineColor(color[il]);
 787     h1[il]->SetLineWidth(1);
 788   }
 789   const Char_t *labMC[] = {"TRD != ESD [bad]", "TRD == ESD [good]", "TRD == -ESD [accept]"};
 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();
 797   hs->Draw(/*"nostack,e1p"*/); leg->Draw();
 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);
 803   hs->SetMinimum(50.);
 804   hs->GetYaxis()->CenterTitle();
 805   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 806 }