PWG1/TRD/AliTRDcheckESD.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 /////////////////////////////////////////////////////
  17 //
  18 // Check basic detector results at ESD level
  19 //   - Geometrical efficiency
  20 //   - Tracking efficiency
  21 //   - PID efficiency
  22 //   - Refit efficiency
  23 //
  24 // Author
  25 //   Alex Bercuci <A.Bercuci@gsi.de>
  26 //
  27 //////////////////////////////////////////////////////
  28
  29 #include <TClonesArray.h>
  30 #include <TCanvas.h>
  31 #include <TObjArray.h>
  32 #include <TPad.h>
  33 #include <TLegend.h>
  34 #include <TLatex.h>
  35 #include <TLine.h>
  36 #include <TF1.h>
  37 #include <TH2I.h>
  38 #include <TH2F.h>
  39 #include <TH3S.h>
  40 #include <TH3F.h>
  41 #include <TProfile2D.h>
  42 #include <TProfile.h>
  43 #include <TGraphErrors.h>
  44 #include <TGraphAsymmErrors.h>
  45 #include <TFile.h>
  46 #include <TTree.h>
  47 #include <TROOT.h>
  48 #include <TChain.h>
  49 #include <TParticle.h>
  50 #include <TTimeStamp.h>
  51
  52 #include "AliLog.h"
  53 #include "AliAnalysisManager.h"
  54 #include "AliESDEvent.h"
  55 #include "AliESDkink.h"
  56 #include "AliMCEvent.h"
  57 #include "AliESDInputHandler.h"
  58 #include "AliMCEventHandler.h"
  59
  60 #include "AliESDtrack.h"
  61 #include "AliMCParticle.h"
  62 #include "AliPID.h"
  63 #include "AliStack.h"
  64 #include "AliTrackReference.h"
  65
  66 #include "AliTRDcheckESD.h"
  67
  68 ClassImp(AliTRDcheckESD)
  69
  70 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
  71 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
  72 const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
  73 8, 4, 2, 20};
  74 FILE* AliTRDcheckESD::fgFile = NULL;
  75
  76 const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;
  77 const Int_t   AliTRDcheckESD::fgkEvVertexN = 1;
  78 const Float_t AliTRDcheckESD::fgkTrkDCAxy  = 40.;
  79 const Float_t AliTRDcheckESD::fgkTrkDCAz   = 15.;
  80 const Int_t   AliTRDcheckESD::fgkNclTPC    = 100;
  81 const Float_t AliTRDcheckESD::fgkPt        = 0.2;
  82 const Float_t AliTRDcheckESD::fgkEta       = 0.9;
  83 const Float_t AliTRDcheckESD::fgkQs        = 0.002;
  84
  85 //____________________________________________________________________
  86 AliTRDcheckESD::AliTRDcheckESD():
  87   AliAnalysisTaskSE()
  88   ,fStatus(0)
  89   ,fNRefFigures(0)
  90   ,fESD(NULL)
  91   ,fMC(NULL)
  92   ,fHistos(NULL)
  93   ,fResults(NULL)
  94 {
  95   //
  96   // Default constructor
  97   //
  98   SetNameTitle("checkESD", "Check TRD @ ESD level");
  99   SetMC(kTRUE);
 100 }
 101
 102 //____________________________________________________________________
 103 AliTRDcheckESD::AliTRDcheckESD(char* name):
 104   AliAnalysisTaskSE(name)
 105   ,fStatus(0)
 106   ,fNRefFigures(0)
 107   ,fESD(NULL)
 108   ,fMC(NULL)
 109   ,fHistos(NULL)
 110   ,fResults(NULL)
 111 {
 112   //
 113   // Default constructor
 114   //
 115   SetMC(kTRUE);
 116   SetTitle("Check TRD @ ESD level");
 117   DefineOutput(1, TObjArray::Class());
 118 }
 119
 120 //____________________________________________________________________
 121 AliTRDcheckESD::~AliTRDcheckESD()
 122 {
 123 // Destructor
 124   if(fHistos){
 125     //fHistos->Delete();
 126     delete fHistos;
 127   }
 128   if(fResults){
 129     fResults->Delete();
 130     delete fResults;
 131   }
 132 }
 133
 134 //____________________________________________________________________
 135 void AliTRDcheckESD::UserCreateOutputObjects()
 136 {
 137   //
 138   // Create Output Containers (TObjectArray containing 1D histograms)
 139   //
 140   Histos();
 141 }
 142
 143 //____________________________________________________________________
 144 void AliTRDcheckESD::MakeSummary(){
 145   TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);
 146   cOut->cd();
 147   GetRefFigure(4);
 148   cOut->SaveAs(Form("TRDsummary%s1.gif", GetName()));
 149
 150   cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 768);
 151   cOut->cd();
 152   GetRefFigure(5);
 153   cOut->SaveAs(Form("TRDsummary%s2.gif", GetName()));
 154 }
 155
 156 //____________________________________________________________________
 157 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
 158 {
 159   if(ifig>=fNRefFigures){
 160     AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
 161     return kFALSE;
 162   }
 163   if(!gPad){
 164     AliWarning("Please provide a canvas to draw results.");
 165     return kFALSE;
 166   } else {
 167     gPad->SetLogx(0);gPad->SetLogy(0);
 168     gPad->SetMargin(0.125, 0.015, 0.1, 0.015);
 169   }
 170
 171   const Char_t *title[20];
 172   Float_t nada(0.0);
 173   TH1 *hF(NULL);
 174   TH1 *hFeffP(NULL); TH1 *hFeffN(NULL);
 175   TH2 *h2F(NULL); TH2 *h2Feff(NULL);
 176   TH2 *h2FtpcP(NULL); TH2 *h2FtpcN(NULL);
 177   TH2 *h2FtrdP(NULL); TH2 *h2FtrdN(NULL);
 178   TH3 *h3F(NULL);
 179   if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
 180   TLegend *leg(NULL);
 181   TList *l(NULL); TVirtualPad *pad(NULL);
 182   TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
 183   TObjArray *arr(NULL);
 184   TProfile2D *hProf2D(NULL);
 185   TProfile *hProf(NULL);
 186   TLatex *lat=new TLatex();
 187   lat->SetTextSize(0.07);
 188   lat->SetTextColor(2);
 189   TLine line;
 190   TTimeStamp now;
 191   TF1* fitFunc(NULL);
 192   switch(ifig){
 193   case kNCl: // number of clusters/track
 194     if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
 195
 196     leg = new TLegend(.83, .7, .99, .96);
 197     leg->SetHeader("Species");
 198     leg->SetBorderSize(0); leg->SetFillStyle(0);
 199     for(Int_t ig(0); ig<fgkNgraph[kNCl]; ig++){
 200       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
 201       if(!g->GetN()) continue;
 202       g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");
 203       if(ig) continue;
 204       hF=g->GetHistogram();
 205       hF->SetXTitle("no of clusters");
 206       hF->SetYTitle("entries");
 207       hF->GetYaxis()->CenterTitle(1);
 208       hF->GetYaxis()->SetTitleOffset(1.2);
 209       hF->SetMinimum(5);
 210     }
 211     leg->Draw(); gPad->SetLogy();
 212     break;
 213   case kTRDstat: // Efficiency
 214     if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
 215     leg = new TLegend(.62, .77, .98, .98);
 216     leg->SetHeader("TRD Efficiency");
 217     leg->SetBorderSize(0); leg->SetFillStyle(0);
 218     title[0] = "Geometrical (TRDin/TPCout)";
 219     title[1] = "Tracking (TRDout/TRDin)";
 220     title[2] = "PID (TRDpid/TRDin)";
 221     title[3] = "Refit (TRDrefit/TRDin)";
 222     hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
 223     hF->SetMaximum(1.4);
 224     hF->GetXaxis()->SetMoreLogLabels();
 225     hF->GetYaxis()->CenterTitle(1);
 226     hF->Draw("p");
 227     for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
 228       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
 229       g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
 230       //PutTrendValue(name[id], g->GetMean(2));
 231       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 232     }
 233     leg->Draw(); gPad->SetLogx();
 234     break;
 235   case kTRDmom: // Energy loss
 236     if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
 237     leg = new TLegend(.65, .7, .95, .99);
 238     leg->SetHeader("Energy Loss");
 239     leg->SetBorderSize(1); leg->SetFillColor(0);
 240     title[0] = "Max & 90% quantile";
 241     title[1] = "Mean & 60% quantile";
 242     hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
 243     hF->SetMaximum(1.3);hF->SetMinimum(-.3);
 244     hF->Draw("p");
 245     for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
 246       if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
 247       ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
 248       //PutTrendValue(name[id], g->GetMean(2));
 249       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 250     }
 251     leg->Draw();gPad->SetLogx(kFALSE);
 252     break;
 253   case kPtRes: // Pt resolution @ vertex
 254     if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
 255     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 256     pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();
 257     pad->SetMargin(0.1, 0.022, 0.1, 0.023);
 258     hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
 259     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 260     hF->GetXaxis()->SetMoreLogLabels();
 261     hF->GetXaxis()->SetTitleOffset(1.2);
 262     hF->GetYaxis()->CenterTitle();
 263     hF->Draw("p");
 264     //for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){
 265     for(Int_t ig(2); ig<6; ig++){
 266       if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 267       if(!g->GetN()) continue;
 268       g->Draw("pl");
 269       //PutTrendValue(name[id], g->GetMean(2));
 270       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 271     }
 272     pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();
 273     pad->SetMargin(0.1, 0.22, 0.1, 0.023);
 274     hF = (TH1*)hF->Clone("hFcheckESD1");
 275     hF->SetTitle("ITS+TPC");
 276     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 277     hF->Draw("p");
 278     leg = new TLegend(.78, .1, .99, .98);
 279     leg->SetHeader("P_{t} @ DCA");
 280     leg->SetBorderSize(1); leg->SetFillColor(0);
 281     leg->SetTextAlign(22);
 282     leg->SetTextFont(12);
 283     leg->SetTextSize(0.03813559);
 284     {
 285       Int_t nPlots(0);
 286       //for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){
 287       for(Int_t ig(12); ig<16; ig++){
 288         if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 289         if(!g->GetN()) continue;
 290         nPlots++;
 291         g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
 292         //PutTrendValue(name[id], g->GetMean(2));
 293         //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 294       }
 295       if(nPlots) leg->Draw();
 296     }
 297     break;
 298   case 4:            // plot a 3x3 canvas with tracking related histograms
 299     gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
 300     gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
 301     gPad->Divide(3,3,0.,0.);
 302     l=gPad->GetListOfPrimitives();
 303     // eta-phi distr. for positive TPC tracks
 304     pad = ((TVirtualPad*)l->At(0)); pad->cd();
 305     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 306     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 307     h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos));
 308     h2FtpcP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
 309     h2FtpcP->GetXaxis()->SetTitle("#eta");
 310     h2FtpcP->GetXaxis()->CenterTitle();
 311     h2FtpcP->GetXaxis()->SetTitleSize(0.07);
 312     h2FtpcP->GetXaxis()->SetTitleOffset(0.8);
 313     h2FtpcP->GetXaxis()->SetLabelSize(0.05);
 314     h2FtpcP->GetYaxis()->SetTitle("detector #varphi");
 315     h2FtpcP->GetYaxis()->CenterTitle();
 316     h2FtpcP->GetYaxis()->SetTitleSize(0.07);
 317     h2FtpcP->GetYaxis()->SetTitleOffset(0.8);
 318     h2FtpcP->GetYaxis()->SetLabelSize(0.05);
 319     h2FtpcP->SetTitle("");
 320     h2FtpcP->Draw("colz");
 321     lat->DrawLatex(-0.9, 3.6, "TPC positive ref. tracks");
 322     //-----------------
 323     // eta-phi distr. for negative TPC tracks
 324     pad = ((TVirtualPad*)l->At(1)); pad->cd();
 325     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 326     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 327     h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg));
 328     h2FtpcN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
 329     h2FtpcN->GetXaxis()->SetTitle("#eta");
 330     h2FtpcN->GetXaxis()->CenterTitle();
 331     h2FtpcN->GetXaxis()->SetTitleSize(0.07);
 332     h2FtpcN->GetXaxis()->SetTitleOffset(0.8);
 333     h2FtpcN->GetXaxis()->SetLabelSize(0.05);
 334     h2FtpcN->GetYaxis()->SetTitle("detector #varphi");
 335     h2FtpcN->GetYaxis()->CenterTitle();
 336     h2FtpcN->GetYaxis()->SetTitleSize(0.07);
 337     h2FtpcN->GetYaxis()->SetTitleOffset(0.8);
 338     h2FtpcN->GetYaxis()->SetLabelSize(0.05);
 339     h2FtpcN->SetTitle("");
 340     h2FtpcN->Draw("colz");
 341     lat->DrawLatex(-0.9, 3.6, "TPC negative ref. tracks");
 342     // eta-phi distr. for positive TRD tracks
 343     pad = ((TVirtualPad*)l->At(3)); pad->cd();
 344     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 345     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 346     h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos));
 347     h2FtrdP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
 348     h2FtrdP->GetXaxis()->SetTitle("#eta");
 349     h2FtrdP->GetXaxis()->CenterTitle();
 350     h2FtrdP->GetXaxis()->SetTitleSize(0.07);
 351     h2FtrdP->GetXaxis()->SetTitleOffset(0.8);
 352     h2FtrdP->GetXaxis()->SetLabelSize(0.05);
 353     h2FtrdP->GetYaxis()->SetTitle("detector #varphi");
 354     h2FtrdP->GetYaxis()->CenterTitle();
 355     h2FtrdP->GetYaxis()->SetTitleSize(0.07);
 356     h2FtrdP->GetYaxis()->SetTitleOffset(0.8);
 357     h2FtrdP->GetYaxis()->SetLabelSize(0.05);
 358     h2FtrdP->SetMaximum(h2FtpcP->GetMaximum());
 359     h2FtrdP->SetTitle("");
 360     h2FtrdP->Draw("colz");
 361     lat->DrawLatex(-0.9, 3.6, "TRD positive ref. tracks");
 362     //-----------------
 363     // eta-phi distr. for negative TRD tracks
 364     pad = ((TVirtualPad*)l->At(4)); pad->cd();
 365     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 366     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 367     h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg));
 368     h2FtrdN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
 369     h2FtrdN->GetXaxis()->SetTitle("#eta");
 370     h2FtrdN->GetXaxis()->CenterTitle();
 371     h2FtrdN->GetXaxis()->SetTitleSize(0.07);
 372     h2FtrdN->GetXaxis()->SetTitleOffset(0.8);
 373     h2FtrdN->GetXaxis()->SetLabelSize(0.05);
 374     h2FtrdN->GetYaxis()->SetTitle("detector #varphi");
 375     h2FtrdN->GetYaxis()->CenterTitle();
 376     h2FtrdN->GetYaxis()->SetTitleSize(0.07);
 377     h2FtrdN->GetYaxis()->SetTitleOffset(0.8);
 378     h2FtrdN->GetYaxis()->SetLabelSize(0.05);
 379     h2FtrdN->SetMaximum(h2FtpcN->GetMaximum());
 380     h2FtrdN->SetTitle("");
 381     h2FtrdN->Draw("colz");
 382     lat->DrawLatex(-0.9, 3.6, "TRD negative ref. tracks");
 383     // eta-phi efficiency for positive TRD tracks
 384     pad = ((TVirtualPad*)l->At(6)); pad->cd();
 385     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 386     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 387     h2Feff = (TH2F*)h2FtrdP->Clone();
 388     h2Feff->Reset();
 389     h2Feff->Divide(h2FtrdP, h2FtpcP);
 390     h2Feff->GetXaxis()->SetTitle("#eta");
 391     h2Feff->GetXaxis()->CenterTitle();
 392     h2Feff->GetXaxis()->SetTitleSize(0.07);
 393     h2Feff->GetXaxis()->SetTitleOffset(0.8);
 394     h2Feff->GetXaxis()->SetLabelSize(0.05);
 395     h2Feff->GetYaxis()->SetTitle("detector #varphi");
 396     h2Feff->GetYaxis()->CenterTitle();
 397     h2Feff->GetYaxis()->SetTitleSize(0.07);
 398     h2Feff->GetYaxis()->SetTitleOffset(0.8);
 399     h2Feff->GetYaxis()->SetLabelSize(0.05);
 400     h2Feff->SetMaximum(1.0);
 401     h2Feff->SetTitle("");
 402     h2Feff->Draw("colz");
 403     lat->DrawLatex(-0.9, 3.6, "Efficiency positive tracks");
 404     // eta-phi efficiency for negative TRD tracks
 405     pad = ((TVirtualPad*)l->At(7)); pad->cd();
 406     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 407     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 408     h2Feff = (TH2F*)h2FtrdN->Clone();
 409     h2Feff->Reset();
 410     h2Feff->Divide(h2FtrdN, h2FtpcN);
 411     h2Feff->GetXaxis()->SetTitle("#eta");
 412     h2Feff->GetXaxis()->CenterTitle();
 413     h2Feff->GetXaxis()->SetTitleSize(0.07);
 414     h2Feff->GetXaxis()->SetTitleOffset(0.8);
 415     h2Feff->GetXaxis()->SetLabelSize(0.05);
 416     h2Feff->GetYaxis()->SetTitle("detector #varphi");
 417     h2Feff->GetYaxis()->CenterTitle();
 418     h2Feff->GetYaxis()->SetTitleSize(0.07);
 419     h2Feff->GetYaxis()->SetTitleOffset(0.8);
 420     h2Feff->GetYaxis()->SetLabelSize(0.05);
 421     h2Feff->SetMaximum(1.0);
 422     h2Feff->SetTitle("");
 423     h2Feff->Draw("colz");
 424     lat->DrawLatex(-0.9, 3.6, "Efficiency negative tracks");
 425
 426     // <ntracklets> vs (phi,eta)
 427     pad = ((TVirtualPad*)l->At(2)); pad->cd();
 428     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 429     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 430     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvNtrkl));
 431     hProf2D->SetStats(kFALSE);
 432     hProf2D->SetTitle("");
 433     hProf2D->GetXaxis()->SetTitle("#eta");
 434     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 435     hProf2D->GetXaxis()->SetTitleSize(0.07);
 436     hProf2D->GetXaxis()->CenterTitle();
 437     hProf2D->GetXaxis()->SetLabelSize(0.05);
 438     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 439     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 440     hProf2D->GetYaxis()->SetTitleSize(0.07);
 441     hProf2D->GetYaxis()->SetLabelSize(0.05);
 442     hProf2D->GetYaxis()->CenterTitle();
 443     hProf2D->SetMinimum(0.);
 444     hProf2D->SetMaximum(6.);
 445     hProf2D->Draw("colz");
 446     lat->DrawLatex(-0.9, 3.6, "TRD <N_{tracklets}>");
 447     // TPC-TRD matching efficiency vs pt
 448     pad = ((TVirtualPad*)l->At(5)); pad->cd();
 449     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
 450     pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
 451     hFeffP = TRDEfficiency(+1);
 452     hFeffN = TRDEfficiency(-1);
 453     h2F=new TH2F("rangeEffPt", "",10,0.,10.,10,0.,1.1);
 454     h2F->SetStats(kFALSE);
 455     h2F->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 456     h2F->GetXaxis()->SetTitleOffset(0.8);
 457     h2F->GetXaxis()->SetTitleSize(0.07);
 458     h2F->GetXaxis()->CenterTitle();
 459     h2F->GetXaxis()->SetLabelSize(0.05);
 460     h2F->GetYaxis()->SetTitle("TRD-TPC matching efficiency");
 461     h2F->GetYaxis()->SetTitleOffset(0.8);
 462     h2F->GetYaxis()->SetTitleSize(0.07);
 463     h2F->GetYaxis()->SetLabelSize(0.05);
 464     h2F->GetYaxis()->CenterTitle();
 465     h2F->Draw();
 466     line.SetLineStyle(2);
 467     line.SetLineWidth(2);
 468     line.DrawLine(h2F->GetXaxis()->GetXmin(), 0.7, h2F->GetXaxis()->GetXmax(), 0.7);
 469     line.DrawLine(h2F->GetXaxis()->GetXmin(), 0.9, h2F->GetXaxis()->GetXmax(), 0.9);
 470     hFeffP->SetMarkerStyle(20);
 471     hFeffP->SetMarkerColor(2);
 472     hFeffN->SetMarkerStyle(22);
 473     hFeffN->SetMarkerColor(4);
 474     hFeffP->Draw("same");
 475     hFeffN->Draw("same");
 476     leg=new TLegend(0.65, 0.2, 0.95, 0.4);
 477     leg->SetFillColor(0);
 478     leg->AddEntry(hFeffP, "positives", "p");
 479     leg->AddEntry(hFeffN, "negatives", "p");
 480     leg->Draw();
 481     // create trending values for the TPC-TRD matching efficiency
 482     // fit the efficiency histos with a constant in the range [1.0,1.5] GeV/c
 483     fitFunc = new TF1("constantFunc","[0]",1.0,1.5);
 484     hFeffP->Fit(fitFunc,"Q0","",1.0,1.5);
 485     PutTrendValue("TrackingEffPos1GeV", fitFunc->GetParameter(0));
 486     PutTrendValue("TrackingEffPos1GeVErr", fitFunc->GetParError(0));
 487     hFeffN->Fit(fitFunc,"Q0","",1.0,1.5);
 488     PutTrendValue("TrackingEffNeg1GeV", fitFunc->GetParameter(0));
 489     PutTrendValue("TrackingEffNeg1GeVErr", fitFunc->GetParError(0));
 490
 491     // Nclusters per TRD track
 492     pad = ((TVirtualPad*)l->At(8)); pad->cd();
 493     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.12);
 494     pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
 495     pad->SetLogz();
 496     h2F = dynamic_cast<TH2F*>(fHistos->At(kNClsTrackTRD));
 497     h2F->SetStats(kFALSE);
 498     h2F->SetTitle("");
 499     h2F->GetXaxis()->SetTitle("p [GeV/c]");
 500     h2F->GetXaxis()->SetTitleOffset(0.8);
 501     h2F->GetXaxis()->SetTitleSize(0.07);
 502     h2F->GetXaxis()->CenterTitle();
 503     h2F->GetXaxis()->SetLabelSize(0.05);
 504     h2F->GetYaxis()->SetTitle("#clusters per TRD track");
 505     h2F->GetYaxis()->SetTitleOffset(0.8);
 506     h2F->GetYaxis()->SetTitleSize(0.07);
 507     h2F->GetYaxis()->CenterTitle();
 508     h2F->GetYaxis()->SetLabelSize(0.05);
 509     h2F->Draw("colz");
 510     break;
 511   case 5:            // plot a 3x3 canvas with PID related histograms
 512     gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
 513     gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
 514     gPad->Divide(3,3,0.,0.);
 515     l=gPad->GetListOfPrimitives();
 516     // eta-phi distr. for <Qtot> in layer 0
 517     pad = ((TVirtualPad*)l->At(0)); pad->cd();
 518     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 519     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 520     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+0));
 521     hProf2D->SetStats(kFALSE);
 522     hProf2D->SetTitle("");
 523     hProf2D->GetXaxis()->SetTitle("#eta");
 524     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 525     hProf2D->GetXaxis()->SetTitleSize(0.07);
 526     hProf2D->GetXaxis()->CenterTitle();
 527     hProf2D->GetXaxis()->SetLabelSize(0.05);
 528     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 529     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 530     hProf2D->GetYaxis()->SetTitleSize(0.07);
 531     hProf2D->GetYaxis()->SetLabelSize(0.05);
 532     hProf2D->GetYaxis()->CenterTitle();
 533     hProf2D->SetMinimum(0.);
 534     hProf2D->SetMaximum(25.);
 535     hProf2D->Draw("colz");
 536     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 0");
 537     // eta-phi distr. for <Qtot> in layer 1
 538     pad = ((TVirtualPad*)l->At(3)); pad->cd();
 539     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 540     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 541     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+1));
 542     hProf2D->SetStats(kFALSE);
 543     hProf2D->SetTitle("");
 544     hProf2D->GetXaxis()->SetTitle("#eta");
 545     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 546     hProf2D->GetXaxis()->SetTitleSize(0.07);
 547     hProf2D->GetXaxis()->CenterTitle();
 548     hProf2D->GetXaxis()->SetLabelSize(0.05);
 549     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 550     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 551     hProf2D->GetYaxis()->SetTitleSize(0.07);
 552     hProf2D->GetYaxis()->SetLabelSize(0.05);
 553     hProf2D->GetYaxis()->CenterTitle();
 554     hProf2D->SetMinimum(0.);
 555     hProf2D->SetMaximum(25.);
 556     hProf2D->Draw("colz");
 557     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 1");
 558     // eta-phi distr. for <Qtot> in layer 2
 559     pad = ((TVirtualPad*)l->At(6)); pad->cd();
 560     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 561     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 562     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+2));
 563     hProf2D->SetStats(kFALSE);
 564     hProf2D->SetTitle("");
 565     hProf2D->GetXaxis()->SetTitle("#eta");
 566     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 567     hProf2D->GetXaxis()->SetTitleSize(0.07);
 568     hProf2D->GetXaxis()->CenterTitle();
 569     hProf2D->GetXaxis()->SetLabelSize(0.05);
 570     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 571     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 572     hProf2D->GetYaxis()->SetTitleSize(0.07);
 573     hProf2D->GetYaxis()->SetLabelSize(0.05);
 574     hProf2D->GetYaxis()->CenterTitle();
 575     hProf2D->SetMinimum(0.);
 576     hProf2D->SetMaximum(25.);
 577     hProf2D->Draw("colz");
 578     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 2");
 579     // eta-phi distr. for <Qtot> in layer 3
 580     pad = ((TVirtualPad*)l->At(1)); pad->cd();
 581     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 582     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 583     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+3));
 584     hProf2D->SetStats(kFALSE);
 585     hProf2D->SetTitle("");
 586     hProf2D->GetXaxis()->SetTitle("#eta");
 587     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 588     hProf2D->GetXaxis()->SetTitleSize(0.07);
 589     hProf2D->GetXaxis()->CenterTitle();
 590     hProf2D->GetXaxis()->SetLabelSize(0.05);
 591     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 592     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 593     hProf2D->GetYaxis()->SetTitleSize(0.07);
 594     hProf2D->GetYaxis()->SetLabelSize(0.05);
 595     hProf2D->GetYaxis()->CenterTitle();
 596     hProf2D->SetMinimum(0.);
 597     hProf2D->SetMaximum(25.);
 598     hProf2D->Draw("colz");
 599     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 3");
 600     // eta-phi distr. for <Qtot> in layer 4
 601     pad = ((TVirtualPad*)l->At(4)); pad->cd();
 602     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 603     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 604     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+4));
 605     hProf2D->SetStats(kFALSE);
 606     hProf2D->SetTitle("");
 607     hProf2D->GetXaxis()->SetTitle("#eta");
 608     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 609     hProf2D->GetXaxis()->SetTitleSize(0.07);
 610     hProf2D->GetXaxis()->CenterTitle();
 611     hProf2D->GetXaxis()->SetLabelSize(0.05);
 612     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 613     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 614     hProf2D->GetYaxis()->SetTitleSize(0.07);
 615     hProf2D->GetYaxis()->SetLabelSize(0.05);
 616     hProf2D->GetYaxis()->CenterTitle();
 617     hProf2D->SetMinimum(0.);
 618     hProf2D->SetMaximum(25.);
 619     hProf2D->Draw("colz");
 620     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 4");
 621     // eta-phi distr. for <Qtot> in layer 5
 622     pad = ((TVirtualPad*)l->At(7)); pad->cd();
 623     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 624     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 625     hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+5));
 626     hProf2D->SetStats(kFALSE);
 627     hProf2D->SetTitle("");
 628     hProf2D->GetXaxis()->SetTitle("#eta");
 629     hProf2D->GetXaxis()->SetTitleOffset(0.8);
 630     hProf2D->GetXaxis()->SetTitleSize(0.07);
 631     hProf2D->GetXaxis()->CenterTitle();
 632     hProf2D->GetXaxis()->SetLabelSize(0.05);
 633     hProf2D->GetYaxis()->SetTitle("detector #varphi");
 634     hProf2D->GetYaxis()->SetTitleOffset(0.8);
 635     hProf2D->GetYaxis()->SetTitleSize(0.07);
 636     hProf2D->GetYaxis()->SetLabelSize(0.05);
 637     hProf2D->GetYaxis()->CenterTitle();
 638     hProf2D->SetMinimum(0.);
 639     hProf2D->SetMaximum(25.);
 640     hProf2D->Draw("colz");
 641     lat->DrawLatex(-0.9, 3.6, "TRD <Q_{tot}> Layer 5");
 642     // PH versus slice number
 643     pad = ((TVirtualPad*)l->At(2)); pad->cd();
 644     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 645     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 646     h2F = dynamic_cast<TH2F*>(fHistos->At(kPHSlice));
 647     h2F->SetStats(kFALSE);
 648     h2F->SetTitle("");
 649     h2F->GetXaxis()->SetTitle("slice");
 650     h2F->GetXaxis()->SetTitleOffset(0.8);
 651     h2F->GetXaxis()->SetTitleSize(0.07);
 652     h2F->GetXaxis()->CenterTitle();
 653     h2F->GetXaxis()->SetLabelSize(0.05);
 654     h2F->GetYaxis()->SetTitle("PH");
 655     h2F->GetYaxis()->SetTitleOffset(0.8);
 656     h2F->GetYaxis()->SetTitleSize(0.07);
 657     h2F->GetYaxis()->SetLabelSize(0.05);
 658     h2F->GetYaxis()->CenterTitle();
 659     h2F->Draw("colz");
 660     //hProf = h2F->ProfileX("profileX");
 661     //hProf->SetLineWidth(2);
 662     //hProf->Draw("same");
 663     // Qtot vs P
 664     pad = ((TVirtualPad*)l->At(5)); pad->cd();
 665     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
 666     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
 667     pad->SetLogz();
 668     h2F = dynamic_cast<TH2F*>(fHistos->At(kQtotP));
 669     h2F->SetStats(kFALSE);
 670     h2F->SetTitle("");
 671     h2F->GetXaxis()->SetTitle("P [GeV/c]");
 672     h2F->GetXaxis()->SetTitleOffset(0.8);
 673     h2F->GetXaxis()->SetTitleSize(0.07);
 674     h2F->GetXaxis()->CenterTitle();
 675     h2F->GetXaxis()->SetLabelSize(0.05);
 676     h2F->GetYaxis()->SetRangeUser(0.0,100.0);
 677     h2F->GetYaxis()->SetTitle("Q_{tot}");
 678     h2F->GetYaxis()->SetTitleOffset(0.8);
 679     h2F->GetYaxis()->SetTitleSize(0.07);
 680     h2F->GetYaxis()->SetLabelSize(0.05);
 681     h2F->GetYaxis()->CenterTitle();
 682     h2F->Draw("colz");
 683     // create trending value for the average Qtot at 1 GeV/c
 684     hProf = h2F->ProfileX("profileQtot",1,h2F->GetYaxis()->FindBin(40.));
 685     PutTrendValue("AvQtot1GeV", hProf->GetBinContent(hProf->GetXaxis()->FindBin(1.)));
 686     PutTrendValue("AvQtot1GeVErr", hProf->GetBinError(hProf->GetXaxis()->FindBin(1.)));
 687     break;
 688   }
 689   return kTRUE;
 690 }
 691
 692 //____________________________________________________________________
 693 void AliTRDcheckESD::UserExec(Option_t *){
 694   //
 695   // Run the Analysis
 696   //
 697   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 698   fMC = MCEvent();
 699
 700   if(!fESD){
 701     AliError("ESD event missing.");
 702     return;
 703   }
 704
 705   // Get MC information if available
 706   AliStack * fStack = NULL;
 707   if(HasMC()){
 708     if(!fMC){
 709       AliWarning("MC event missing");
 710       SetMC(kFALSE);
 711     } else {
 712       if(!(fStack = fMC->Stack())){
 713         AliWarning("MC stack missing");
 714         SetMC(kFALSE);
 715       }
 716     }
 717   }
 718   TH1 *h(NULL);
 719
 720   // fill event vertex histos
 721   h = (TH1F*)fHistos->At(kTPCVertex);
 722   if(fESD->GetPrimaryVertexTPC()) h->Fill(fESD->GetPrimaryVertexTPC()->GetZv());
 723   h = (TH1F*)fHistos->At(kEventVertex);
 724   if(fESD->GetPrimaryVertex()) h->Fill(fESD->GetPrimaryVertex()->GetZv());
 725   // fill the uncutted number of tracks
 726   h = (TH1I*)fHistos->At(kNTracksAll);
 727   h->Fill(fESD->GetNumberOfTracks());
 728
 729   // counters for number of tracks in acceptance&DCA and for those with a minimum of TPC clusters
 730   Int_t nTracksAcc=0;
 731   Int_t nTracksTPC=0;
 732
 733   AliESDtrack *esdTrack(NULL);
 734   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 735     esdTrack = fESD->GetTrack(itrk);
 736
 737     // track status
 738     ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
 739
 740     // track selection
 741     Bool_t selected(kTRUE);
 742     if(esdTrack->Pt() < fgkPt){
 743       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESD->GetEventNumberInFile(), itrk, esdTrack->Pt()));
 744       selected = kFALSE;
 745     }
 746     if(TMath::Abs(esdTrack->Eta()) > fgkEta){
 747       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));
 748       selected = kFALSE;
 749     }
 750     if(!Bool_t(status & AliESDtrack::kTPCout)){
 751       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] !TPCout", fESD->GetEventNumberInFile(), itrk));
 752       selected = kFALSE;
 753     }
 754     if(esdTrack->GetKinkIndex(0) > 0){
 755       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Kink", fESD->GetEventNumberInFile(), itrk));
 756       selected = kFALSE;
 757     }
 758
 759     Float_t par[2], cov[3];
 760     esdTrack->GetImpactParameters(par, cov);
 761     if(selected && esdTrack->GetTPCNcls()>=10) {
 762       // fill DCA histograms
 763       h = (TH1F*)fHistos->At(kDCAxy); h->Fill(par[0]);
 764       h = (TH1F*)fHistos->At(kDCAz); h->Fill(par[1]);
 765       // fill pt distribution at this stage
 766       h = (TH1F*)fHistos->At(kPt1); h->Fill(esdTrack->Pt());
 767     }
 768     if(IsCollision()){ // cuts on DCA
 769       if(TMath::Abs(par[0]) > fgkTrkDCAxy){
 770         AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));
 771         selected = kFALSE;
 772       }
 773       if(TMath::Abs(par[1]) > fgkTrkDCAz){
 774         AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));
 775         selected = kFALSE;
 776       }
 777     }
 778     Float_t theta=esdTrack->Theta();
 779     Float_t phi=esdTrack->Phi();
 780     Int_t nClustersTPC = esdTrack->GetTPCNcls();
 781     Float_t eta=-TMath::Log(TMath::Tan(theta/2.));
 782     if(selected) {
 783       nTracksAcc++;   // number of tracks in acceptance and DCA cut
 784       // fill pt distribution at this stage
 785       h = (TH1F*)fHistos->At(kPt2); h->Fill(esdTrack->Pt());
 786       // TPC nclusters distribution
 787       h = (TH1I*)fHistos->At(kNTPCCl); h->Fill(nClustersTPC);
 788       if(esdTrack->Pt()>1.0) {
 789         h = (TH1I*)fHistos->At(kNTPCCl2); h->Fill(nClustersTPC);
 790       }
 791       // (eta,nclustersTPC) distrib of TPC ref. tracks
 792       h = (TH2F*)fHistos->At(kEtaNclsTPC); h->Fill(eta, nClustersTPC);
 793       // (phi,nclustersTPC) distrib of TPC ref. tracks
 794       h = (TH2F*)fHistos->At(kPhiNclsTPC); h->Fill(phi, nClustersTPC);
 795
 796     }
 797
 798     if(nClustersTPC < fgkNclTPC){
 799       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESD->GetEventNumberInFile(), itrk, nClustersTPC));
 800       selected = kFALSE;
 801     }
 802     if(!selected) continue;
 803
 804     // number of TPC reference tracks
 805     nTracksTPC++;
 806
 807     Int_t nTRD(esdTrack->GetNcls(2));
 808     Double_t pt(esdTrack->Pt());
 809     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 810     // pid quality
 811     Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
 812
 813     TH3F *hhh;
 814     // find position and momentum of the track at entrance in TRD
 815     Double_t localCoord[3] = {0., 0., 0.};
 816     Bool_t localCoordGood = esdTrack->GetXYZAt(298., fESD->GetMagneticField(), localCoord);
 817     if(localCoordGood) {
 818       hhh = (TH3F*)fHistos->At(kPropagXYvsP); hhh->Fill(localCoord[0], localCoord[1], esdTrack->GetP());
 819       hhh = (TH3F*)fHistos->At(kPropagRZvsP); hhh->Fill(localCoord[2], TMath::Sqrt(localCoord[0]*localCoord[0]+localCoord[1]*localCoord[1]), esdTrack->GetP());
 820     }
 821     Double_t localMom[3] = {0., 0., 0.};
 822     Bool_t localMomGood = esdTrack->GetPxPyPzAt(298., fESD->GetMagneticField(), localMom);
 823     Double_t localPhi = (localMomGood ? TMath::ATan2(localMom[1], localMom[0]) : 0.0);
 824     Double_t localSagitaPhi = (localCoordGood ? TMath::ATan2(localCoord[1], localCoord[0]) : 0.0);
 825
 826     // fill pt distribution at this stage
 827     if(esdTrack->Charge()>0) {
 828       h = (TH1F*)fHistos->At(kPt3pos); h->Fill(pt);
 829       // fill eta-phi map of TPC positive ref. tracks
 830       if(localCoordGood && localMomGood) {
 831         hhh = (TH3F*)fHistos->At(kTPCRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);
 832       }
 833     }
 834     if(esdTrack->Charge()<0) {
 835       h = (TH1F*)fHistos->At(kPt3neg); h->Fill(pt);
 836       // fill eta-phi map of TPC negative ref. tracks
 837       if(localCoordGood && localMomGood) {
 838         hhh = (TH3F*)fHistos->At(kTPCRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);
 839       }
 840     }
 841     // TPC dE/dx vs P
 842     h = (TH2F*)fHistos->At(kTPCDedx); h->Fill(esdTrack->GetP(), esdTrack->GetTPCsignal());
 843     // (eta,phi) distrib of TPC ref. tracks
 844     h = (TH2F*)fHistos->At(kEtaPhi); h->Fill(eta, phi);
 845
 846     Int_t nTRDtrkl = esdTrack->GetTRDntracklets();
 847     // TRD reference tracks
 848     if(nTRDtrkl>=1) {
 849       // fill pt distribution at this stage
 850       if(esdTrack->Charge()>0) {
 851         h = (TH1F*)fHistos->At(kPt4pos); h->Fill(pt);
 852         // fill eta-phi map of TRD positive ref. tracks
 853         if(localCoordGood && localMomGood) {
 854           hhh = (TH3F*)fHistos->At(kTRDRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);
 855         }
 856       }
 857       if(esdTrack->Charge()<0) {
 858         h = (TH1F*)fHistos->At(kPt4neg); h->Fill(pt);
 859         // fill eta-phi map of TRD negative ref. tracks
 860         if(localCoordGood && localMomGood) {
 861           hhh = (TH3F*)fHistos->At(kTRDRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);
 862         }
 863       }
 864       TProfile2D *h2d;
 865       // fill eta-phi map of TRD negative ref. tracks
 866       if(localCoordGood && localMomGood) {
 867         h2d = (TProfile2D*)fHistos->At(kTRDEtaPhiAvNtrkl); h2d->Fill(eta, localSagitaPhi, (Float_t)nTRDtrkl);
 868         h2d = (TProfile2D*)fHistos->At(kTRDEtaDeltaPhiAvNtrkl); h2d->Fill(eta, localPhi-localSagitaPhi, (Float_t)nTRDtrkl);
 869       }
 870       // ntracklets/track vs P
 871       h = (TH2F*)fHistos->At(kNTrackletsTRD); h->Fill(esdTrack->GetP(), nTRDtrkl);
 872       // ntracklets/track vs P
 873       h = (TH2F*)fHistos->At(kNClsTrackTRD); h->Fill(esdTrack->GetP(), esdTrack->GetTRDncls());
 874       // (slicePH,sliceNo) distribution and Qtot from slices
 875       for(Int_t iPlane=0; iPlane<6; iPlane++) {
 876         Float_t Qtot=0;
 877         for(Int_t iSlice=0; iSlice<8; iSlice++) {
 878           if(esdTrack->GetTRDslice(iPlane, iSlice)>20.) {
 879             h = (TH2F*)fHistos->At(kPHSlice); h->Fill(iSlice, esdTrack->GetTRDslice(iPlane, iSlice));
 880             Qtot += esdTrack->GetTRDslice(iPlane, iSlice);
 881           }
 882         }
 883         // Qtot>100 to avoid noise
 884         if(Qtot>100.) {
 885           h = (TH2F*)fHistos->At(kQtotP); h->Fill(esdTrack->GetTRDmomentum(iPlane), fgkQs*Qtot);
 886         }
 887         // Qtot>100 to avoid noise
 888         // fgkQs*Qtot<40. so that the average will give a value close to the peak
 889         if(localCoordGood && localMomGood && Qtot>100. && fgkQs*Qtot<40.) {
 890           h2d = (TProfile2D*)fHistos->At(kTRDEtaPhiAvQtot+iPlane);
 891           h2d->Fill(eta, localSagitaPhi, fgkQs*Qtot);
 892         }
 893       }
 894       // theta distribution
 895       h = (TH1F*)fHistos->At(kTheta); h->Fill(theta);
 896       h = (TH1F*)fHistos->At(kPhi); h->Fill(phi);
 897     }  // end if nTRDtrkl>=1
 898
 899     // look at external track param
 900     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 901     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
 902
 903     Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
 904     // read MC info if available
 905     Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
 906     AliMCParticle *mcParticle(NULL);
 907     if(HasMC()){
 908       AliTrackReference *ref(NULL);
 909       Int_t fLabel(esdTrack->GetLabel());
 910       Int_t fIdx(TMath::Abs(fLabel));
 911       if(fIdx > fStack->GetNtrack()) continue;
 912
 913       // read MC particle
 914       if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
 915         AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 916         continue;
 917       }
 918       pt0  = mcParticle->Pt();
 919       eta0 = mcParticle->Eta();
 920       phi0 = mcParticle->Phi();
 921       kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
 922
 923       // read track references
 924       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 925       if(!nRefs){
 926         AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
 927         continue;
 928       }
 929       Int_t iref = 0;
 930       while(iref<nRefs){
 931         ref = mcParticle->GetTrackReference(iref);
 932         if(ref->LocalX() > fgkxTPC) break;
 933         ref=NULL; iref++;
 934       }
 935       if(ref){
 936         if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 937           ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 938         }
 939       } else { // track stopped in TPC
 940         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 941       }
 942       ptTRD = ref->Pt();kFOUND=kTRUE;
 943     } else { // use reconstructed values
 944       if(op){
 945         Double_t x(op->GetX());
 946         if(x<fgkxTOF && x>fgkxTPC){
 947           ptTRD=op->Pt();
 948           kFOUND=kTRUE;
 949         }
 950       }
 951
 952       if(!kFOUND && ip){
 953         ptTRD=ip->Pt();
 954         kFOUND=kTRUE;
 955       }
 956     }     // end if(HasMC())
 957
 958     if(kFOUND){
 959       h = (TH2I*)fHistos->At(kTRDstat);
 960       if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
 961       if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
 962       if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);
 963       if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
 964       if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
 965     }
 966     Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)
 967          ,sgn(esdTrack->Charge()<0?0:1);
 968     if(kBarrel && kPhysPrim) {
 969       TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
 970       Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
 971       h3->Fill(pt0, 1.e2*(pt/pt0-1.),
 972         offset + 2*idx + sgn);
 973     }
 974     ((TH1*)fHistos->At(kNCl))->Fill(nTRD, 2*idx + sgn);
 975     if(ip){
 976       h = (TH2I*)fHistos->At(kTRDmom);
 977       Float_t pTRD(0.);
 978       for(Int_t ily=6; ily--;){
 979         if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
 980         h->Fill(ip->GetP()-pTRD, ily);
 981       }
 982     }
 983   }  // end loop over tracks
 984
 985   // fill the number of tracks histograms
 986   h = (TH1I*)fHistos->At(kNTracksAcc);
 987   h->Fill(nTracksAcc);
 988   h = (TH1I*)fHistos->At(kNTracksTPC);
 989   h->Fill(nTracksTPC);
 990
 991   PostData(1, fHistos);
 992 }
 993
 994 //____________________________________________________________________
 995 TObjArray* AliTRDcheckESD::Histos()
 996 {
 997 // Retrieve histograms array if already build or build it
 998
 999   if(fHistos) return fHistos;
1000
1001   fHistos = new TObjArray(kNhistos);
1002   //fHistos->SetOwner(kTRUE);
1003
1004   TH1 *h = NULL;
1005
1006   // clusters per track
1007   const Int_t kNpt(30);
1008   Float_t Pt(0.2);
1009   Float_t binsPt[kNpt+1];
1010   for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;
1011   if(!(h = (TH2S*)gROOT->FindObject("hNCl"))){
1012     h = new TH2S("hNCl", "Clusters per TRD track;N_{cl}^{TRD};SPECIES;entries", 60, 0., 180., 10, -0.5, 9.5);
1013     TAxis *ay(h->GetYaxis());
1014     ay->SetLabelOffset(0.015);
1015     for(Int_t i(0); i<AliPID::kSPECIES; i++){
1016       ay->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
1017       ay->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
1018     }
1019   } else h->Reset();
1020   fHistos->AddAt(h, kNCl); fNRefFigures++;
1021
1022   // status bits histogram
1023   const Int_t kNbits(5);
1024   Float_t Bits(.5);
1025   Float_t binsBits[kNbits+1];
1026   for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
1027   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
1028     h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
1029     TAxis *ay(h->GetYaxis());
1030     ay->SetBinLabel(1, "kTPCout");
1031     ay->SetBinLabel(2, "kTRDin");
1032     ay->SetBinLabel(3, "kTRDout");
1033     ay->SetBinLabel(4, "kTRDpid");
1034     ay->SetBinLabel(5, "kTRDrefit");
1035   } else h->Reset();
1036   fHistos->AddAt(h, kTRDstat);
1037
1038   // energy loss
1039   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
1040     h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
1041   } else h->Reset();
1042   fHistos->AddAt(h, kTRDmom);
1043   //if(!HasMC()) return fHistos;
1044
1045   // pt resolution
1046   const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
1047   Float_t DPt(-3.), Spec(-0.5);
1048   Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
1049   for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
1050   for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
1051   if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
1052     h = new TH3S("hPtRes", "P_{t} resolution @ DCA;p_{t}^{MC} [GeV/c];#Delta p_{t}/p_{t}^{MC} [%];SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspec, binsSpec);
1053     TAxis *az(h->GetZaxis());
1054     az->SetLabelOffset(0.015);
1055     for(Int_t i(0); i<AliPID::kSPECIES; i++){
1056       az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
1057       az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
1058       az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
1059       az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
1060     }
1061   } else h->Reset();
1062   fHistos->AddAt(h, kPtRes);
1063
1064   // TPC event vertex distribution
1065   if(!(h = (TH1F*)gROOT->FindObject("hTPCVertex"))){
1066     h = new TH1F("hTPCVertex", "Event vertex Z coord. from TPC tracks", 100, -25., 25.);
1067   } else h->Reset();
1068   fHistos->AddAt(h, kTPCVertex);
1069
1070   // Event vertex
1071   if(!(h = (TH1F*)gROOT->FindObject("hEventVertex"))){
1072     h = new TH1F("hEventVertex", "Event vertex Z coord.", 100, -25., 25.);
1073   } else h->Reset();
1074   fHistos->AddAt(h, kEventVertex);
1075
1076   // Number of all tracks
1077   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAll"))){
1078     h = new TH1I("hNTracksAll", "Number of tracks per event, event vertex cuts", 5000, 0, 5000);
1079   } else h->Reset();
1080   fHistos->AddAt(h, kNTracksAll);
1081
1082   // Number of tracks in acceptance and DCA cut
1083   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAcc"))){
1084     h = new TH1I("hNTracksAcc", Form("Number of tracks per event, |#eta|<%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1085                                      fgkEta, fgkTrkDCAxy, fgkTrkDCAz), 5000, 0, 5000);
1086   } else h->Reset();
1087   fHistos->AddAt(h, kNTracksAcc);
1088
1089   // Number of tracks in TPC (Ncls>10)
1090   if(!(h = (TH1I*)gROOT->FindObject("hNTracksTPC"))){
1091     h = new TH1I("hNTracksTPC", Form("Number of tracks per event, |#eta|<%.1f, pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1092                                      fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 5000, 0, 5000);
1093   } else h->Reset();
1094   fHistos->AddAt(h, kNTracksTPC);
1095
1096   // Distribution of DCA-xy
1097   if(!(h = (TH1F*)gROOT->FindObject("hDCAxy"))){
1098     h = new TH1F("hDCAxy", "Distribution of transverse DCA", 100, -100., 100.);
1099   } else h->Reset();
1100   fHistos->AddAt(h, kDCAxy);
1101
1102   // Distribution of DCA-z
1103   if(!(h = (TH1F*)gROOT->FindObject("hDCAz"))){
1104     h = new TH1F("hDCAz", "Distribution of longitudinal DCA", 100, -100., 100.);
1105   } else h->Reset();
1106   fHistos->AddAt(h, kDCAz);
1107
1108   Float_t binPtLimits[33] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1109                              1.0, 1.1, 1.2, 1.3, 1.4,
1110                              1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,
1111                              3.4, 3.8, 4.2, 4.6, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
1112   // Pt distributions
1113   if(!(h = (TH1F*)gROOT->FindObject("hPt1"))){
1114     h = new TH1F("hPt1", Form("dN/dpt, |#eta|<%.1f and pt>%.1f", fgkEta, fgkPt), 32, binPtLimits);
1115   } else h->Reset();
1116   fHistos->AddAt(h, kPt1);
1117
1118   if(!(h = (TH1F*)gROOT->FindObject("hPt2"))){
1119     h = new TH1F("hPt2", Form("dN/dpt, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1120                               fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 32, binPtLimits);
1121   } else h->Reset();
1122   fHistos->AddAt(h, kPt2);
1123
1124   if(!(h = (TH1F*)gROOT->FindObject("hPt3pos"))){
1125     h = new TH1F("hPt3pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1126                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
1127   } else h->Reset();
1128   fHistos->AddAt(h, kPt3pos);
1129
1130   if(!(h = (TH1F*)gROOT->FindObject("hPt3neg"))){
1131     h = new TH1F("hPt3neg", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1132                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
1133   } else h->Reset();
1134   fHistos->AddAt(h, kPt3neg);
1135
1136   if(!(h = (TH1F*)gROOT->FindObject("hPt4pos"))){
1137     h = new TH1F("hPt4pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
1138                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
1139   } else h->Reset();
1140   fHistos->AddAt(h, kPt4pos);
1141
1142   if(!(h = (TH1F*)gROOT->FindObject("hPt4neg"))){
1143     h = new TH1F("hPt4pos", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
1144                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
1145   } else h->Reset();
1146   fHistos->AddAt(h, kPt4neg);
1147
1148   // theta distribution of TRD tracks
1149   if(!(h = (TH1F*)gROOT->FindObject("hTheta"))){
1150     h = new TH1F("hTheta", Form("dN/d#theta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
1151                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 220,.5,2.7);
1152   } else h->Reset();
1153   fHistos->AddAt(h, kTheta);
1154
1155   // phi distribution of TRD tracks
1156   if(!(h = (TH1F*)gROOT->FindObject("hPhi"))){
1157     h = new TH1F("hPhi", Form("dN/d#varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
1158                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 157,0,6.28);
1159   } else h->Reset();
1160   fHistos->AddAt(h, kPhi);
1161
1162   // TPC cluster distribution
1163   if(!(h = (TH1F*)gROOT->FindObject("hNTPCCl"))){
1164     h = new TH1I("hNTPCCl", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1165                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
1166   } else h->Reset();
1167   fHistos->AddAt(h, kNTPCCl);
1168
1169   if(!(h = (TH1I*)gROOT->FindObject("hNTPCCl2"))){
1170     h = new TH1F("hNTPCCl2", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, pt>1.0 GeV/c",
1171                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
1172   } else h->Reset();
1173   fHistos->AddAt(h, kNTPCCl2);
1174
1175   // dE/dx vs P for TPC reference tracks
1176   if(!(h = (TH2F*)gROOT->FindObject("hTPCDedx"))){
1177     h = new TH2F("hTPCDedx", Form("TPC dE/dx vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1178                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 120, 0,600.);
1179   } else h->Reset();
1180   fHistos->AddAt(h, kTPCDedx);
1181
1182   // eta,phi distribution of TPC reference tracks
1183   if(!(h = (TH2F*)gROOT->FindObject("hEtaPhi"))){
1184     h = new TH2F("hEtaPhi", Form("TPC (#eta,#varphi), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1185                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 50, -1, 1, 157, 0, 6.28);
1186   } else h->Reset();
1187   fHistos->AddAt(h, kEtaPhi);
1188
1189   // Nclusters vs eta distribution for TPC tracks
1190   if(!(h = (TH2F*)gROOT->FindObject("hEtaNclsTPC"))){
1191     h = new TH2F("hEtaNclsTPC", Form("TPC Nclusters vs. #eta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1192                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 50, -1, 1, 160, 0, 160.);
1193   } else h->Reset();
1194   fHistos->AddAt(h, kEtaNclsTPC);
1195
1196   // Nclusters vs phi distribution for TPC reference tracks
1197   if(!(h = (TH2F*)gROOT->FindObject("hPhiNclsTPC"))){
1198     h = new TH2F("hPhiNclsTPC", Form("TPC Nclusters vs. #varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1199                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 157, 0, 6.28, 160, 0, 160.);
1200   } else h->Reset();
1201   fHistos->AddAt(h, kPhiNclsTPC);
1202
1203   // Ntracklets/track vs P for TRD reference tracks
1204   Double_t binsP[19] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.7, 2.0,
1205                         2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 9.0, 12.0};
1206   if(!(h = (TH2F*)gROOT->FindObject("hNTrackletsTRD"))){
1207     h = new TH2F("hNTrackletsTRD", Form("TRD Ntracklets/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1208                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 7, -0.5, 6.5);
1209   } else h->Reset();
1210   fHistos->AddAt(h, kNTrackletsTRD);
1211
1212   // Nclusters/track vs P for TRD reference tracks
1213   if(!(h = (TH2F*)gROOT->FindObject("hNClsTrackTRD"))){
1214     h = new TH2F("hNClsTrackTRD", Form("TRD Nclusters/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1215                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 180, 0., 180.);
1216   } else h->Reset();
1217   fHistos->AddAt(h, kNClsTrackTRD);
1218
1219   // <PH> vs slice number for TRD reference tracklets
1220   if(!(h = (TH2F*)gROOT->FindObject("hPHSlice"))){
1221     h = new TH2F("hPHSlice", Form("<PH> vs sliceNo, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1222                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 2000, 0., 2000.);
1223   } else h->Reset();
1224   fHistos->AddAt(h, kPHSlice);
1225
1226   // Qtot vs P for TRD reference tracklets
1227   if(!(h = (TH2F*)gROOT->FindObject("hQtotP"))){
1228     h = new TH2F("hQtotP", Form("Qtot(from slices) vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1229                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 400, 0., 200);
1230   } else h->Reset();
1231   fHistos->AddAt(h, kQtotP);
1232
1233   // (X,Y,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
1234   if(!(h = (TH3F*)gROOT->FindObject("hPropagXYvsP"))){
1235     h = new TH3F("hPropagXYvsP", Form("(x,y) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1236                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-500,500, 100,-500,500, 10, 0.,10.);
1237   } else h->Reset();
1238   fHistos->AddAt(h, kPropagXYvsP);
1239
1240   // (R,Z,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
1241   if(!(h = (TH3F*)gROOT->FindObject("hPropagRZvsP"))){
1242     h = new TH3F("hPropagRZvsP", Form("(r,z) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1243                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-350., 350., 100,0.,500., 10, 0.,10.);
1244   } else h->Reset();
1245   fHistos->AddAt(h, kPropagRZvsP);
1246
1247   Float_t etaBinLimits[101];
1248   for(Int_t i=0; i<101; i++) etaBinLimits[i] = -1.0 + i*2.0/100.;
1249   Float_t phiBinLimits[151];
1250   for(Int_t i=0; i<151; i++) phiBinLimits[i] = -1.1*TMath::Pi() + i*2.2*TMath::Pi()/150.;
1251   // (eta,detector phi,P) distribution of reference TPC positive tracks
1252   if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksPos"))){
1253     h = new TH3F("hTPCRefTracksPos", Form("(#eta,detector #varphi,p) for TPC positive reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1254                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1255   } else h->Reset();
1256   fHistos->AddAt(h, kTPCRefTracksPos);
1257
1258   // (eta,detector phi,P) distribution of reference TPC negative tracks
1259   if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksNeg"))){
1260     h = new TH3F("hTPCRefTracksNeg", Form("(#eta,detector #varphi,p) for TPC negative reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1261                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1262   } else h->Reset();
1263   fHistos->AddAt(h, kTPCRefTracksNeg);
1264
1265   // (eta,detector phi,P) distribution of reference TRD positive tracks
1266   if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksPos"))){
1267     h = new TH3F("hTRDRefTracksPos", Form("(#eta,detector #varphi,p) for TRD positive reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1268                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1269   } else h->Reset();
1270   fHistos->AddAt(h, kTRDRefTracksPos);
1271
1272   // (eta,detector phi,P) distribution of reference TRD negative tracks
1273   if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksNeg"))){
1274     h = new TH3F("hTRDRefTracksNeg", Form("(#eta,detector #varphi,p) for TRD negative reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1275                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1276   } else h->Reset();
1277   fHistos->AddAt(h, kTRDRefTracksNeg);
1278
1279   // (eta,detector phi) profile of average number of TRD tracklets/track
1280   if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaPhiAvNtrkl"))){
1281     h = new TProfile2D("hTRDEtaPhiAvNtrkl", Form("<Ntracklets/track> vs (#eta,detector #varphi) for TRD reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1282                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 150, -1.1*TMath::Pi(), 1.1*TMath::Pi());
1283   } else h->Reset();
1284   fHistos->AddAt(h, kTRDEtaPhiAvNtrkl);
1285
1286   // (eta,delta phi) profile of average number of TRD tracklets/track
1287   if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaDeltaPhiAvNtrkl"))){
1288     h = new TProfile2D("hTRDEtaDeltaPhiAvNtrkl", Form("<Ntracklets/track> vs (#eta, #Delta#varphi) for TRD reference tracks, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1289                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -0.4*TMath::Pi(), 0.4*TMath::Pi());
1290   } else h->Reset();
1291   fHistos->AddAt(h, kTRDEtaDeltaPhiAvNtrkl);
1292
1293   // (eta, detector phi) profile of average tracklet Qtot from slices
1294   for(Int_t iLayer=0;iLayer<6;iLayer++) {
1295     if(!(h = (TProfile2D*)gROOT->FindObject(Form("hTRDEtaPhiAvQtot_Layer%d",iLayer)))) {
1296       h = new TProfile2D(Form("hTRDEtaPhiAvQtot_Layer%d",iLayer),
1297                          Form("<Q_{tot}> vs (#eta, detector #varphi) for TRD reference tracks (layer %d), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1298                               iLayer, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -1.1*TMath::Pi(), 1.1*TMath::Pi());
1299     } else h->Reset();
1300     fHistos->AddAt(h, kTRDEtaPhiAvQtot+iLayer);
1301   }
1302
1303   return fHistos;
1304 }
1305
1306 //____________________________________________________________________
1307 Bool_t AliTRDcheckESD::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
1308 {
1309 // Load data from performance file
1310
1311   if(!TFile::Open(file)){
1312     AliWarning(Form("Couldn't open file %s.", file));
1313     return kFALSE;
1314   }
1315   if(dir){
1316     if(!gFile->cd(dir)){
1317       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
1318       return kFALSE;
1319     }
1320   }
1321   TObjArray *o(NULL);
1322   const Char_t *tn=(name ? name : GetName());
1323   if(!(o = (TObjArray*)gDirectory->Get(tn))){
1324     AliWarning(Form("Missing histogram container %s.", tn));
1325     return kFALSE;
1326   }
1327   fHistos = (TObjArray*)o->Clone(GetName());
1328   gFile->Close();
1329   return kTRUE;
1330 }
1331
1332 //_______________________________________________________
1333 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
1334 {
1335 // Dump trending value to default file
1336
1337   if(!fgFile){
1338     fgFile = fopen("TRD.Performance.txt", "at");
1339   }
1340   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
1341   return kTRUE;
1342 }
1343
1344 //____________________________________________________________________
1345 void AliTRDcheckESD::Terminate(Option_t *)
1346 {
1347 // Steer post-processing
1348   if(!fHistos){
1349     fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
1350     if(!fHistos){
1351       AliError("Histogram container not found in output");
1352       return;
1353     }
1354   }
1355
1356   const Char_t *name[kNrefs] = {
1357     "Ncl", "Eff", "Eloss", "PtResDCA"
1358   };
1359   TObjArray *arr(NULL); TGraph *g(NULL);
1360   if(!fResults){
1361     fResults = new TObjArray(kNrefs);
1362     fResults->SetOwner();
1363     fResults->SetName("results");
1364     for(Int_t iref(0); iref<kNrefs; iref++){
1365       fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
1366       arr->SetName(name[iref]);  arr->SetOwner();
1367       switch(iref){
1368       case kNCl:
1369         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1370           arr->AddAt(g = new TGraphErrors(), ig);
1371           g->SetLineColor(ig+1);
1372           g->SetMarkerColor(ig+1);
1373           g->SetMarkerStyle(ig+20);
1374           g->SetName(Form("s%d", ig));
1375           switch(ig){
1376           case 0: g->SetTitle("ALL"); break;
1377           case 1: g->SetTitle("NEG"); break;
1378           case 2: g->SetTitle("POS"); break;
1379           default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;
1380           };
1381         }
1382         break;
1383       case kTRDmom:
1384         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1385           arr->AddAt(g = new TGraphAsymmErrors(), ig);
1386           g->SetLineColor(ig+1);
1387           g->SetMarkerColor(ig+1);
1388           g->SetMarkerStyle(ig+20);
1389         }
1390         break;
1391       case kPtRes:
1392         for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
1393           Int_t ig(2*idx);
1394           arr->AddAt(g = new TGraphErrors(), ig);
1395           g->SetLineColor(kRed-idx);
1396           g->SetMarkerColor(kRed-idx);
1397           g->SetMarkerStyle(20+idx);
1398           g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
1399           arr->AddAt(g = new TGraphErrors(), ig+1);
1400           g->SetLineColor(kBlue-idx);
1401           g->SetMarkerColor(kBlue-idx);
1402           g->SetMarkerStyle(20+idx);
1403           g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
1404
1405           ig+=10;
1406           arr->AddAt(g = new TGraphErrors(), ig);
1407           g->SetLineColor(kRed-idx);
1408           g->SetMarkerColor(kRed-idx);
1409           g->SetMarkerStyle(20+idx);
1410           g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
1411           arr->AddAt(g = new TGraphErrors(), ig+1);
1412           g->SetLineColor(kBlue-idx);
1413           g->SetMarkerColor(kBlue-idx);
1414           g->SetMarkerStyle(20+idx);
1415           g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
1416         }
1417         break;
1418       default:
1419         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1420           arr->AddAt(g = new TGraphErrors(), ig);
1421           g->SetLineColor(ig+1);
1422           g->SetMarkerColor(ig+1);
1423           g->SetMarkerStyle(ig+20);
1424         }
1425         break;
1426       }
1427     }
1428   }
1429   TH1 *h1[2] = {NULL, NULL};
1430   TH2I *h2(NULL);
1431   TAxis *ax(NULL);
1432
1433   // No of clusters
1434   if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;
1435   ax = h2->GetXaxis();
1436   arr = (TObjArray*)fResults->At(kNCl);
1437   // All tracks
1438   h1[0] = h2->ProjectionX("Ncl_px");
1439   TGraphErrors *ge=(TGraphErrors*)arr->At(0);
1440   for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1441     ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1442   }
1443   // All charged tracks
1444   TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};
1445   hNclCh[0]->Reset();hNclCh[1]->Reset();
1446   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1447     hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg
1448     hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is));     // pos
1449   }
1450   if(Int_t(hNclCh[0]->GetEntries())){
1451     ge=(TGraphErrors*)arr->At(1);
1452     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1453       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));
1454     }
1455   }
1456   if(Int_t(hNclCh[1]->GetEntries())){
1457     ge=(TGraphErrors*)arr->At(2);
1458     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1459       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));
1460     }
1461   }
1462   // Species wise
1463   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1464     h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);
1465     if(!Int_t(h1[0]->GetEntries())) continue;
1466     ge=(TGraphErrors*)arr->At(2+is);
1467     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1468       ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1469     }
1470   }
1471   fNRefFigures = 1;
1472
1473   // EFFICIENCY
1474   // geometrical efficiency
1475   if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
1476   arr = (TObjArray*)fResults->At(kTRDstat);
1477   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
1478   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
1479   Process(h1, (TGraphErrors*)arr->At(0));
1480   delete h1[0];delete h1[1];
1481   // tracking efficiency
1482   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
1483   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
1484   Process(h1, (TGraphErrors*)arr->At(1));
1485   delete h1[1];
1486   // PID efficiency
1487   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
1488   Process(h1, (TGraphErrors*)arr->At(2));
1489   delete h1[1];
1490   // Refit efficiency
1491   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
1492   Process(h1, (TGraphErrors*)arr->At(3));
1493   delete h1[1];
1494   fNRefFigures++;
1495
1496   // ENERGY LOSS
1497   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
1498   arr = (TObjArray*)fResults->At(kTRDmom);
1499   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
1500   ax=h2->GetXaxis();
1501   const Int_t nq(4);
1502   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
1503   Double_t yq[nq];
1504   for(Int_t ily=6; ily--;){
1505     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
1506     h1[0]->GetQuantiles(nq,yq,xq);
1507     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
1508     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
1509     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
1510     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
1511
1512     //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
1513     delete h1[0];
1514   }
1515   fNRefFigures++;
1516 //  if(!HasMC()) return;
1517
1518   // Pt RESOLUTION @ DCA
1519   TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
1520   if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
1521   arr = (TObjArray*)fResults->At(kPtRes);
1522   TAxis *az(h3->GetZaxis());
1523   for(Int_t i(0); i<AliPID::kSPECIES; i++){
1524     Int_t idx(2*i);
1525     az->SetRange(idx+1, idx+2);
1526     gg[1] = (TGraphErrors*)arr->At(idx);
1527     gg[0] = (TGraphErrors*)arr->At(idx+1);
1528     Process2D((TH2*)h3->Project3D("yx"), gg);
1529
1530     idx+=10;
1531     az->SetRange(idx+1, idx+2);
1532     gg[1] = (TGraphErrors*)arr->At(idx);
1533     gg[0] = (TGraphErrors*)arr->At(idx+1);
1534     Process2D((TH2*)h3->Project3D("yx"), gg);
1535   }
1536   fNRefFigures++;
1537
1538   // 3x3 tracking summary canvas
1539   fNRefFigures++;
1540   // 3x3 PID summary canvas
1541   fNRefFigures++;
1542 }
1543
1544 //____________________________________________________________________
1545 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
1546 {
1547   switch(pdg){
1548   case kElectron:
1549   case kPositron: return AliPID::kElectron;
1550   case kMuonPlus:
1551   case kMuonMinus: return AliPID::kMuon;
1552   case kPiPlus:
1553   case kPiMinus: return AliPID::kPion;
1554   case kKPlus:
1555   case kKMinus: return AliPID::kKaon;
1556   case kProton:
1557   case kProtonBar: return AliPID::kProton;
1558   }
1559   return -1;
1560 }
1561
1562 //____________________________________________________________________
1563 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
1564 {
1565 // Generic function to process one reference plot
1566
1567   Int_t n1 = 0, n2 = 0, ip=0;
1568   Double_t eff = 0.;
1569
1570   TAxis *ax = h1[0]->GetXaxis();
1571   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
1572     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
1573     n2 = (Int_t)h1[1]->GetBinContent(ib);
1574     eff = n2/Float_t(n1);
1575
1576     ip=g->GetN();
1577     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
1578     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
1579   }
1580 }
1581 //________________________________________________________
1582 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
1583 {
1584   //
1585   // Do the processing
1586   //
1587
1588   Int_t n = 0;
1589   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1590   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1591   TF1 f("fg", "gaus", -3.,3.);
1592   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1593     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1594     TH1D *h = h2->ProjectionY("py", ibin, ibin);
1595     if(h->GetEntries()<100) continue;
1596     //AdjustF1(h, f);
1597
1598     h->Fit(&f, "QN");
1599     Int_t ip = g[0]->GetN();
1600     g[0]->SetPoint(ip, x, f.GetParameter(1));
1601     g[0]->SetPointError(ip, 0., f.GetParError(1));
1602     g[1]->SetPoint(ip, x, f.GetParameter(2));
1603     g[1]->SetPointError(ip, 0., f.GetParError(2));
1604   }
1605   return;
1606 }
1607 //____________________________________________________________________
1608 void AliTRDcheckESD::PrintStatus(ULong_t status)
1609 {
1610 // Dump track status to stdout
1611
1612   printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"
1613     ,Bool_t(status & AliESDtrack::kITSin)
1614     ,Bool_t(status & AliESDtrack::kITSout)
1615     ,Bool_t(status & AliESDtrack::kITSrefit)
1616     ,Bool_t(status & AliESDtrack::kTPCin)
1617     ,Bool_t(status & AliESDtrack::kTPCout)
1618     ,Bool_t(status & AliESDtrack::kTPCrefit)
1619     ,Bool_t(status & AliESDtrack::kTPCpid)
1620     ,Bool_t(status & AliESDtrack::kTRDin)
1621     ,Bool_t(status & AliESDtrack::kTRDout)
1622     ,Bool_t(status & AliESDtrack::kTRDrefit)
1623     ,Bool_t(status & AliESDtrack::kTRDpid)
1624     ,Bool_t(status & AliESDtrack::kTRDStop)
1625     ,Bool_t(status & AliESDtrack::kHMPIDout)
1626     ,Bool_t(status & AliESDtrack::kHMPIDpid)
1627   );
1628 }
1629 //____________________________________________________________________
1630 TH2F* AliTRDcheckESD::Proj3D(TH3F* hist, TH2F* accMap, Int_t zbinLow, Int_t zbinHigh, Float_t &entries) {
1631   //
1632   //  Project a 3D histogram to a 2D histogram in the Z axis interval [zbinLow,zbinHigh]
1633   //  Return the 2D histogram and also the number of entries into this projection (entries)
1634
1635   Int_t nBinsX = hist->GetXaxis()->GetNbins();   // X and Y axis bins are assumed to be all equal
1636   Float_t minX = hist->GetXaxis()->GetXmin();
1637   Float_t maxX = hist->GetXaxis()->GetXmax();
1638   Int_t nBinsY = hist->GetYaxis()->GetNbins();
1639   Float_t minY = hist->GetYaxis()->GetXmin();
1640   Float_t maxY = hist->GetYaxis()->GetXmax();
1641   Int_t nBinsZ = hist->GetZaxis()->GetNbins();  // Z axis bins (pt) might have different widths
1642   //Float_t minZ = hist->GetZaxis()->GetXmin();
1643   //Float_t maxZ = hist->GetZaxis()->GetXmax();
1644
1645   TH2F* projHisto = (TH2F*)gROOT->FindObject("projHisto");
1646   if(projHisto)
1647     projHisto->Reset();
1648   else
1649     projHisto = new TH2F("projHisto", "projection", nBinsX, minX, maxX, nBinsY, minY, maxY);
1650
1651   const Double_t kMinAccFraction = 0.1;
1652   entries = 0.0;
1653   Double_t maxAcc = (accMap ? accMap->GetMaximum() : -1);
1654
1655   for(Int_t iZ=1; iZ<=nBinsZ; iZ++) {
1656     if(iZ<zbinLow) continue;
1657     if(iZ>zbinHigh) continue;
1658     for(Int_t iX=1; iX<=nBinsX; iX++) {
1659       for(Int_t iY=1; iY<=nBinsY; iY++) {
1660         if(accMap && maxAcc>0) {
1661           if(accMap->GetBinContent(iX,iY)/maxAcc>kMinAccFraction)
1662             projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));
1663         }
1664         else    // no acc. cut
1665           projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));
1666         // count only the entries which are inside the acceptance map
1667         if(accMap && maxAcc>0) {
1668           if(accMap->GetBinContent(iX,iY)/maxAcc>kMinAccFraction)
1669             entries+=hist->GetBinContent(iX,iY,iZ);
1670         }
1671         else    // no acc. cut
1672           entries+=hist->GetBinContent(iX,iY,iZ);
1673       }
1674     }
1675   }
1676   return projHisto;
1677 }
1678 //____________________________________________________________________
1679 TH1F* AliTRDcheckESD::TRDEfficiency(Short_t positives) {
1680   //
1681   // Calculate the TRD-TPC matching efficiency as function of pt
1682   //
1683   TH3F* tpc3D(NULL); TH3F* trd3D(NULL);
1684   if(positives>0) {      // get the histos for positive tracks
1685     if(!(tpc3D = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos)))) return 0x0;
1686     if(!(trd3D = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos)))) return 0x0;
1687   }
1688   else {      // get the histos for positive tracks
1689     if(!(tpc3D = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg)))) return 0x0;
1690     if(!(trd3D = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg)))) return 0x0;
1691   }
1692
1693   Int_t nBinsZ = trd3D->GetZaxis()->GetNbins();
1694   // project everything on the eta-phi map to obtain an acceptance map (make sure there is enough statistics)
1695   Float_t nada = 0.;
1696   TH2F *trdAcc = (TH2F*)Proj3D(trd3D, 0x0, 1, nBinsZ, nada)->Clone();
1697   // get the bin limits from the Z axis of 3D histos
1698   Float_t *ptBinLimits = new Float_t[nBinsZ+1];
1699   for(Int_t i=1; i<=nBinsZ; i++) {
1700     ptBinLimits[i-1] = trd3D->GetZaxis()->GetBinLowEdge(i);
1701   }
1702   ptBinLimits[nBinsZ] = trd3D->GetZaxis()->GetBinUpEdge(nBinsZ);
1703   TH1F *efficiency = new TH1F("eff", "TRD-TPC matching efficiency", nBinsZ, ptBinLimits);
1704   // loop over Z bins
1705   for(Int_t i=1; i<=nBinsZ; i++) {
1706     Float_t tpcEntries = 0.0; Float_t trdEntries = 0.0;
1707     Proj3D(tpc3D, trdAcc, i, i, tpcEntries);
1708     Proj3D(trd3D, trdAcc, i, i, trdEntries);
1709     Float_t ratio = 0;
1710     if(tpcEntries>0) ratio = trdEntries/tpcEntries;
1711     Float_t error = 0;
1712     if(tpcEntries>0 && trdEntries>0)
1713       error = ratio*TMath::Sqrt((1.0/tpcEntries)+(1.0/trdEntries));
1714     if(ratio>0.0) {
1715       efficiency->SetBinContent(i,ratio);
1716       efficiency->SetBinError(i,error);
1717     }
1718   }     // end loop over Z bins
1719
1720   efficiency->SetLineColor((positives>0 ? 2 : 4));
1721   return efficiency;
1722 }