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