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