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