PWG1/TRD/AliTRDcheckESD.cxx

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