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