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