PWG1/TRD/AliTRDcheckESD.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /////////////////////////////////////////////////////
  17 //
  18 // Check basic detector results at ESD level
  19 //   - Geometrical efficiency
  20 //   - Tracking efficiency
  21 //   - PID efficiency
  22 //   - Refit efficiency
  23 //
  24 // Author
  25 //   Alex Bercuci <A.Bercuci@gsi.de>
  26 //
  27 //////////////////////////////////////////////////////
  28
  29 #include <TClonesArray.h>
  30 #include <TObjArray.h>
  31 #include <TPad.h>
  32 #include <TLegend.h>
  33 #include <TF1.h>
  34 #include <TH2I.h>
  35 #include <TH2F.h>
  36 #include <TH3S.h>
  37 #include <TH3F.h>
  38 #include <TProfile2D.h>
  39 #include <TGraphErrors.h>
  40 #include <TGraphAsymmErrors.h>
  41 #include <TFile.h>
  42 #include <TTree.h>
  43 #include <TROOT.h>
  44 #include <TChain.h>
  45 #include <TParticle.h>
  46
  47 #include "AliLog.h"
  48 #include "AliAnalysisManager.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESDkink.h"
  51 #include "AliMCEvent.h"
  52 #include "AliESDInputHandler.h"
  53 #include "AliMCEventHandler.h"
  54
  55 #include "AliESDtrack.h"
  56 #include "AliMCParticle.h"
  57 #include "AliPID.h"
  58 #include "AliStack.h"
  59 #include "AliTrackReference.h"
  60
  61 #include "AliTRDcheckESD.h"
  62
  63 ClassImp(AliTRDcheckESD)
  64
  65 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
  66 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
  67 const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
  68 8, 4, 2, 20};
  69 FILE* AliTRDcheckESD::fgFile = NULL;
  70
  71 const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;
  72 const Int_t   AliTRDcheckESD::fgkEvVertexN = 1;
  73 const Float_t AliTRDcheckESD::fgkTrkDCAxy  = 40.;
  74 const Float_t AliTRDcheckESD::fgkTrkDCAz   = 15.;
  75 const Int_t   AliTRDcheckESD::fgkNclTPC    = 100;
  76 const Float_t AliTRDcheckESD::fgkPt        = 0.2;
  77 const Float_t AliTRDcheckESD::fgkEta       = 0.9;
  78 const Float_t AliTRDcheckESD::fgkQs        = 0.002;
  79
  80 //____________________________________________________________________
  81 AliTRDcheckESD::AliTRDcheckESD():
  82   AliAnalysisTaskSE()
  83   ,fStatus(0)
  84   ,fNRefFigures(0)
  85   ,fESD(NULL)
  86   ,fMC(NULL)
  87   ,fHistos(NULL)
  88   ,fResults(NULL)
  89 {
  90   //
  91   // Default constructor
  92   //
  93   SetNameTitle("checkESD", "Check TRD @ ESD level");
  94   SetMC(kTRUE);
  95 }
  96
  97 //____________________________________________________________________
  98 AliTRDcheckESD::AliTRDcheckESD(char* name):
  99   AliAnalysisTaskSE(name)
 100   ,fStatus(0)
 101   ,fNRefFigures(0)
 102   ,fESD(NULL)
 103   ,fMC(NULL)
 104   ,fHistos(NULL)
 105   ,fResults(NULL)
 106 {
 107   //
 108   // Default constructor
 109   //
 110   SetMC(kTRUE);
 111   SetTitle("Check TRD @ ESD level");
 112   DefineOutput(1, TObjArray::Class());
 113 }
 114
 115 //____________________________________________________________________
 116 AliTRDcheckESD::~AliTRDcheckESD()
 117 {
 118 // Destructor
 119   if(fHistos){
 120     //fHistos->Delete();
 121     delete fHistos;
 122   }
 123   if(fResults){
 124     fResults->Delete();
 125     delete fResults;
 126   }
 127 }
 128
 129 //____________________________________________________________________
 130 void AliTRDcheckESD::UserCreateOutputObjects()
 131 {
 132   //
 133   // Create Output Containers (TObjectArray containing 1D histograms)
 134   //
 135   Histos();
 136 }
 137
 138
 139 //____________________________________________________________________
 140 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
 141 {
 142   if(ifig>=fNRefFigures){
 143     AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
 144     return kFALSE;
 145   }
 146   if(!gPad){
 147     AliWarning("Please provide a canvas to draw results.");
 148     return kFALSE;
 149   } else {
 150     gPad->SetLogx(0);gPad->SetLogy(0);
 151     gPad->SetMargin(0.125, 0.015, 0.1, 0.015);
 152   }
 153
 154   const Char_t *title[20];
 155   TH1 *hF(NULL);
 156   if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
 157   TLegend *leg(NULL);
 158   TList *l(NULL); TVirtualPad *pad(NULL);
 159   TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
 160   TObjArray *arr(NULL);
 161   switch(ifig){
 162   case kNCl: // number of clusters/track
 163     if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
 164
 165     leg = new TLegend(.83, .7, .99, .96);
 166     leg->SetHeader("Species");
 167     leg->SetBorderSize(0); leg->SetFillStyle(0);
 168     for(Int_t ig(0); ig<fgkNgraph[kNCl]; ig++){
 169       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
 170       if(!g->GetN()) continue;
 171       g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");
 172       if(ig) continue;
 173       hF=g->GetHistogram();
 174       hF->SetXTitle("no of clusters");
 175       hF->SetYTitle("entries");
 176       hF->GetYaxis()->CenterTitle(1);
 177       hF->GetYaxis()->SetTitleOffset(1.2);
 178       hF->SetMinimum(5);
 179     }
 180     leg->Draw(); gPad->SetLogy();
 181     break;
 182   case kTRDstat: // Efficiency
 183     if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
 184     leg = new TLegend(.62, .77, .98, .98);
 185     leg->SetHeader("TRD Efficiency");
 186     leg->SetBorderSize(0); leg->SetFillStyle(0);
 187     title[0] = "Geometrical (TRDin/TPCout)";
 188     title[1] = "Tracking (TRDout/TRDin)";
 189     title[2] = "PID (TRDpid/TRDin)";
 190     title[3] = "Refit (TRDrefit/TRDin)";
 191     hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
 192     hF->SetMaximum(1.4);
 193     hF->GetXaxis()->SetMoreLogLabels();
 194     hF->GetYaxis()->CenterTitle(1);
 195     hF->Draw("p");
 196     for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
 197       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
 198       g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
 199       //PutTrendValue(name[id], g->GetMean(2));
 200       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 201     }
 202     leg->Draw(); gPad->SetLogx();
 203     break;
 204   case kTRDmom: // Energy loss
 205     if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
 206     leg = new TLegend(.65, .7, .95, .99);
 207     leg->SetHeader("Energy Loss");
 208     leg->SetBorderSize(1); leg->SetFillColor(0);
 209     title[0] = "Max & 90% quantile";
 210     title[1] = "Mean & 60% quantile";
 211     hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
 212     hF->SetMaximum(1.3);hF->SetMinimum(-.3);
 213     hF->Draw("p");
 214     for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
 215       if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
 216       ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
 217       //PutTrendValue(name[id], g->GetMean(2));
 218       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 219     }
 220     leg->Draw();gPad->SetLogx(kFALSE);
 221     break;
 222   case kPtRes: // Pt resolution @ vertex
 223     if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
 224     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 225     pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();
 226     pad->SetMargin(0.1, 0.022, 0.1, 0.023);
 227     hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
 228     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 229     hF->GetXaxis()->SetMoreLogLabels();
 230     hF->GetXaxis()->SetTitleOffset(1.2);
 231     hF->GetYaxis()->CenterTitle();
 232     hF->Draw("p");
 233     //for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){
 234     for(Int_t ig(2); ig<6; ig++){
 235       if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 236       if(!g->GetN()) continue;
 237       g->Draw("pl");
 238       //PutTrendValue(name[id], g->GetMean(2));
 239       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 240     }
 241     pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();
 242     pad->SetMargin(0.1, 0.22, 0.1, 0.023);
 243     hF = (TH1*)hF->Clone("hFcheckESD1");
 244     hF->SetTitle("ITS+TPC");
 245     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 246     hF->Draw("p");
 247     leg = new TLegend(.78, .1, .99, .98);
 248     leg->SetHeader("P_{t} @ DCA");
 249     leg->SetBorderSize(1); leg->SetFillColor(0);
 250     leg->SetTextAlign(22);
 251     leg->SetTextFont(12);
 252     leg->SetTextSize(0.03813559);
 253     {
 254       Int_t nPlots(0);
 255       //for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){
 256       for(Int_t ig(12); ig<16; ig++){
 257         if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 258         if(!g->GetN()) continue;
 259         nPlots++;
 260         g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
 261         //PutTrendValue(name[id], g->GetMean(2));
 262         //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 263       }
 264       if(nPlots) leg->Draw();
 265     }
 266
 267     break;
 268   }
 269   return kTRUE;
 270 }
 271
 272 //____________________________________________________________________
 273 void AliTRDcheckESD::UserExec(Option_t *){
 274   //
 275   // Run the Analysis
 276   //
 277   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 278   fMC = MCEvent();
 279
 280   if(!fESD){
 281     AliError("ESD event missing.");
 282     return;
 283   }
 284
 285   // Get MC information if available
 286   AliStack * fStack = NULL;
 287   if(HasMC()){
 288     if(!fMC){
 289       AliWarning("MC event missing");
 290       SetMC(kFALSE);
 291     } else {
 292       if(!(fStack = fMC->Stack())){
 293         AliWarning("MC stack missing");
 294         SetMC(kFALSE);
 295       }
 296     }
 297   }
 298   TH1 *h(NULL);
 299
 300   // fill event vertex histos
 301   h = (TH1F*)fHistos->At(kTPCVertex);
 302   if(fESD->GetPrimaryVertexTPC()) h->Fill(fESD->GetPrimaryVertexTPC()->GetZv());
 303   h = (TH1F*)fHistos->At(kEventVertex);
 304   if(fESD->GetPrimaryVertex()) h->Fill(fESD->GetPrimaryVertex()->GetZv());
 305   // fill the uncutted number of tracks
 306   h = (TH1I*)fHistos->At(kNTracksAll);
 307   h->Fill(fESD->GetNumberOfTracks());
 308
 309   // counters for number of tracks in acceptance&DCA and for those with a minimum of TPC clusters
 310   Int_t nTracksAcc=0;
 311   Int_t nTracksTPC=0;
 312
 313   AliESDtrack *esdTrack(NULL);
 314   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 315     esdTrack = fESD->GetTrack(itrk);
 316
 317     // track status
 318     ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
 319
 320     // track selection
 321     Bool_t selected(kTRUE);
 322     if(esdTrack->Pt() < fgkPt){
 323       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESD->GetEventNumberInFile(), itrk, esdTrack->Pt()));
 324       selected = kFALSE;
 325     }
 326     if(TMath::Abs(esdTrack->Eta()) > fgkEta){
 327       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));
 328       selected = kFALSE;
 329     }
 330     if(!Bool_t(status & AliESDtrack::kTPCout)){
 331       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] !TPCout", fESD->GetEventNumberInFile(), itrk));
 332       selected = kFALSE;
 333     }
 334     if(esdTrack->GetKinkIndex(0) > 0){
 335       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Kink", fESD->GetEventNumberInFile(), itrk));
 336       selected = kFALSE;
 337     }
 338
 339     Float_t par[2], cov[3];
 340     esdTrack->GetImpactParameters(par, cov);
 341     if(selected && esdTrack->GetTPCNcls()>=10) {
 342       // fill DCA histograms
 343       h = (TH1F*)fHistos->At(kDCAxy); h->Fill(par[0]);
 344       h = (TH1F*)fHistos->At(kDCAz); h->Fill(par[1]);
 345       // fill pt distribution at this stage
 346       h = (TH1F*)fHistos->At(kPt1); h->Fill(esdTrack->Pt());
 347     }
 348     if(IsCollision()){ // cuts on DCA
 349       if(TMath::Abs(par[0]) > fgkTrkDCAxy){
 350         AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));
 351         selected = kFALSE;
 352       }
 353       if(TMath::Abs(par[1]) > fgkTrkDCAz){
 354         AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESD->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));
 355         selected = kFALSE;
 356       }
 357     }
 358     Float_t theta=esdTrack->Theta();
 359     Float_t phi=esdTrack->Phi();
 360     Int_t nClustersTPC = esdTrack->GetTPCNcls();
 361     Float_t eta=-TMath::Log(TMath::Tan(theta/2.));
 362     if(selected) {
 363       nTracksAcc++;   // number of tracks in acceptance and DCA cut
 364       // fill pt distribution at this stage
 365       h = (TH1F*)fHistos->At(kPt2); h->Fill(esdTrack->Pt());
 366       // TPC nclusters distribution
 367       h = (TH1I*)fHistos->At(kNTPCCl); h->Fill(nClustersTPC);
 368       if(esdTrack->Pt()>1.0) {
 369         h = (TH1I*)fHistos->At(kNTPCCl2); h->Fill(nClustersTPC);
 370       }
 371       // (eta,nclustersTPC) distrib of TPC ref. tracks
 372       h = (TH2F*)fHistos->At(kEtaNclsTPC); h->Fill(eta, nClustersTPC);
 373       // (phi,nclustersTPC) distrib of TPC ref. tracks
 374       h = (TH2F*)fHistos->At(kPhiNclsTPC); h->Fill(phi, nClustersTPC);
 375
 376     }
 377
 378     if(nClustersTPC < fgkNclTPC){
 379       AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESD->GetEventNumberInFile(), itrk, nClustersTPC));
 380       selected = kFALSE;
 381     }
 382     if(!selected) continue;
 383
 384     // number of TPC reference tracks
 385     nTracksTPC++;
 386
 387     Int_t nTRD(esdTrack->GetNcls(2));
 388     Double_t pt(esdTrack->Pt());
 389     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 390     // pid quality
 391     Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
 392
 393     TH3F *hhh;
 394     // find position and momentum of the track at entrance in TRD
 395     Double_t localCoord[3] = {0., 0., 0.};
 396     Bool_t localCoordGood = esdTrack->GetXYZAt(298., fESD->GetMagneticField(), localCoord);
 397     if(localCoordGood) {
 398       hhh = (TH3F*)fHistos->At(kPropagXYvsP); hhh->Fill(localCoord[0], localCoord[1], esdTrack->GetP());
 399       hhh = (TH3F*)fHistos->At(kPropagRZvsP); hhh->Fill(localCoord[2], TMath::Sqrt(localCoord[0]*localCoord[0]+localCoord[1]*localCoord[1]), esdTrack->GetP());
 400     }
 401     Double_t localMom[3] = {0., 0., 0.};
 402     Bool_t localMomGood = esdTrack->GetPxPyPzAt(298., fESD->GetMagneticField(), localMom);
 403     Double_t localPhi = (localMomGood ? TMath::ATan2(localMom[1], localMom[0]) : 0.0);
 404     Double_t localSagitaPhi = (localCoordGood ? TMath::ATan2(localCoord[1], localCoord[0]) : 0.0);
 405
 406     // fill pt distribution at this stage
 407     if(esdTrack->Charge()>0) {
 408       h = (TH1F*)fHistos->At(kPt3pos); h->Fill(pt);
 409       // fill eta-phi map of TPC positive ref. tracks
 410       if(localCoordGood && localMomGood) {
 411         hhh = (TH3F*)fHistos->At(kTPCRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);
 412       }
 413     }
 414     if(esdTrack->Charge()<0) {
 415       h = (TH1F*)fHistos->At(kPt3neg); h->Fill(pt);
 416       // fill eta-phi map of TPC negative ref. tracks
 417       if(localCoordGood && localMomGood) {
 418         hhh = (TH3F*)fHistos->At(kTPCRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);
 419       }
 420     }
 421     // TPC dE/dx vs P
 422     h = (TH2F*)fHistos->At(kTPCDedx); h->Fill(esdTrack->GetP(), esdTrack->GetTPCsignal());
 423     // (eta,phi) distrib of TPC ref. tracks
 424     h = (TH2F*)fHistos->At(kEtaPhi); h->Fill(eta, phi);
 425
 426     Int_t nTRDtrkl = esdTrack->GetTRDntracklets();
 427     // TRD reference tracks
 428     if(nTRDtrkl>=1) {
 429       // fill pt distribution at this stage
 430       if(esdTrack->Charge()>0) {
 431         h = (TH1F*)fHistos->At(kPt4pos); h->Fill(pt);
 432         // fill eta-phi map of TRD positive ref. tracks
 433         if(localCoordGood && localMomGood) {
 434           hhh = (TH3F*)fHistos->At(kTRDRefTracksPos); hhh->Fill(eta, localSagitaPhi, pt);
 435         }
 436       }
 437       if(esdTrack->Charge()<0) {
 438         h = (TH1F*)fHistos->At(kPt4neg); h->Fill(pt);
 439         // fill eta-phi map of TRD negative ref. tracks
 440         if(localCoordGood && localMomGood) {
 441           hhh = (TH3F*)fHistos->At(kTRDRefTracksNeg); hhh->Fill(eta, localSagitaPhi, pt);
 442         }
 443       }
 444       TProfile2D *h2d;
 445       // fill eta-phi map of TRD negative ref. tracks
 446       if(localCoordGood && localMomGood) {
 447         h2d = (TProfile2D*)fHistos->At(kTRDEtaPhiAvNtrkl); h2d->Fill(eta, localSagitaPhi, (Float_t)nTRDtrkl);
 448         h2d = (TProfile2D*)fHistos->At(kTRDEtaDeltaPhiAvNtrkl); h2d->Fill(eta, localPhi-localSagitaPhi, (Float_t)nTRDtrkl);
 449       }
 450       // ntracklets/track vs P
 451       h = (TH2F*)fHistos->At(kNTrackletsTRD); h->Fill(esdTrack->GetP(), nTRDtrkl);
 452       // ntracklets/track vs P
 453       h = (TH2F*)fHistos->At(kNClsTrackTRD); h->Fill(esdTrack->GetP(), esdTrack->GetTRDncls());
 454       // (slicePH,sliceNo) distribution and Qtot from slices
 455       for(Int_t iPlane=0; iPlane<6; iPlane++) {
 456         Float_t Qtot=0;
 457         for(Int_t iSlice=0; iSlice<8; iSlice++) {
 458           if(esdTrack->GetTRDslice(iPlane, iSlice)>20.) {
 459             h = (TH2F*)fHistos->At(kPHSlice); h->Fill(iSlice, esdTrack->GetTRDslice(iPlane, iSlice));
 460             Qtot += esdTrack->GetTRDslice(iPlane, iSlice);
 461           }
 462         }
 463         h = (TH2F*)fHistos->At(kQtotP); h->Fill(esdTrack->GetTRDmomentum(iPlane), fgkQs*Qtot);
 464       }
 465       // theta distribution
 466       h = (TH1F*)fHistos->At(kTheta); h->Fill(theta);
 467       h = (TH1F*)fHistos->At(kPhi); h->Fill(phi);
 468     }  // end if nTRDtrkl>=1
 469
 470     // look at external track param
 471     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 472     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
 473
 474     Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
 475     // read MC info if available
 476     Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
 477     AliMCParticle *mcParticle(NULL);
 478     if(HasMC()){
 479       AliTrackReference *ref(NULL);
 480       Int_t fLabel(esdTrack->GetLabel());
 481       Int_t fIdx(TMath::Abs(fLabel));
 482       if(fIdx > fStack->GetNtrack()) continue;
 483
 484       // read MC particle
 485       if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
 486         AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 487         continue;
 488       }
 489       pt0  = mcParticle->Pt();
 490       eta0 = mcParticle->Eta();
 491       phi0 = mcParticle->Phi();
 492       kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
 493
 494       // read track references
 495       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 496       if(!nRefs){
 497         AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
 498         continue;
 499       }
 500       Int_t iref = 0;
 501       while(iref<nRefs){
 502         ref = mcParticle->GetTrackReference(iref);
 503         if(ref->LocalX() > fgkxTPC) break;
 504         ref=NULL; iref++;
 505       }
 506       if(ref){
 507         if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 508           ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 509         }
 510       } else { // track stopped in TPC
 511         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 512       }
 513       ptTRD = ref->Pt();kFOUND=kTRUE;
 514     } else { // use reconstructed values
 515       if(op){
 516         Double_t x(op->GetX());
 517         if(x<fgkxTOF && x>fgkxTPC){
 518           ptTRD=op->Pt();
 519           kFOUND=kTRUE;
 520         }
 521       }
 522
 523       if(!kFOUND && ip){
 524         ptTRD=ip->Pt();
 525         kFOUND=kTRUE;
 526       }
 527     }     // end if(HasMC())
 528
 529     if(kFOUND){
 530       h = (TH2I*)fHistos->At(kTRDstat);
 531       if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
 532       if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
 533       if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);
 534       if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
 535       if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
 536     }
 537     Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)
 538          ,sgn(esdTrack->Charge()<0?0:1);
 539     if(kBarrel && kPhysPrim) {
 540       TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
 541       Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
 542       h3->Fill(pt0, 1.e2*(pt/pt0-1.),
 543         offset + 2*idx + sgn);
 544     }
 545     ((TH1*)fHistos->At(kNCl))->Fill(nTRD, 2*idx + sgn);
 546     if(ip){
 547       h = (TH2I*)fHistos->At(kTRDmom);
 548       Float_t pTRD(0.);
 549       for(Int_t ily=6; ily--;){
 550         if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
 551         h->Fill(ip->GetP()-pTRD, ily);
 552       }
 553     }
 554   }  // end loop over tracks
 555
 556   // fill the number of tracks histograms
 557   h = (TH1I*)fHistos->At(kNTracksAcc);
 558   h->Fill(nTracksAcc);
 559   h = (TH1I*)fHistos->At(kNTracksTPC);
 560   h->Fill(nTracksTPC);
 561
 562   PostData(1, fHistos);
 563 }
 564
 565 //____________________________________________________________________
 566 TObjArray* AliTRDcheckESD::Histos()
 567 {
 568 // Retrieve histograms array if already build or build it
 569
 570   if(fHistos) return fHistos;
 571
 572   fHistos = new TObjArray(kNhistos);
 573   //fHistos->SetOwner(kTRUE);
 574
 575   TH1 *h = NULL;
 576
 577   // clusters per track
 578   const Int_t kNpt(30);
 579   Float_t Pt(0.2);
 580   Float_t binsPt[kNpt+1];
 581   for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;
 582   if(!(h = (TH2S*)gROOT->FindObject("hNCl"))){
 583     h = new TH2S("hNCl", "Clusters per TRD track;N_{cl}^{TRD};SPECIES;entries", 60, 0., 180., 10, -0.5, 9.5);
 584     TAxis *ay(h->GetYaxis());
 585     ay->SetLabelOffset(0.015);
 586     for(Int_t i(0); i<AliPID::kSPECIES; i++){
 587       ay->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
 588       ay->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
 589     }
 590   } else h->Reset();
 591   fHistos->AddAt(h, kNCl); fNRefFigures++;
 592
 593   // status bits histogram
 594   const Int_t kNbits(5);
 595   Float_t Bits(.5);
 596   Float_t binsBits[kNbits+1];
 597   for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
 598   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 599     h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
 600     TAxis *ay(h->GetYaxis());
 601     ay->SetBinLabel(1, "kTPCout");
 602     ay->SetBinLabel(2, "kTRDin");
 603     ay->SetBinLabel(3, "kTRDout");
 604     ay->SetBinLabel(4, "kTRDpid");
 605     ay->SetBinLabel(5, "kTRDrefit");
 606   } else h->Reset();
 607   fHistos->AddAt(h, kTRDstat);
 608
 609   // energy loss
 610   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
 611     h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
 612   } else h->Reset();
 613   fHistos->AddAt(h, kTRDmom);
 614   //if(!HasMC()) return fHistos;
 615
 616   // pt resolution
 617   const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
 618   Float_t DPt(-3.), Spec(-0.5);
 619   Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
 620   for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
 621   for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
 622   if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
 623     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);
 624     TAxis *az(h->GetZaxis());
 625     az->SetLabelOffset(0.015);
 626     for(Int_t i(0); i<AliPID::kSPECIES; i++){
 627       az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
 628       az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
 629       az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
 630       az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
 631     }
 632   } else h->Reset();
 633   fHistos->AddAt(h, kPtRes);
 634
 635   // TPC event vertex distribution
 636   if(!(h = (TH1F*)gROOT->FindObject("hTPCVertex"))){
 637     h = new TH1F("hTPCVertex", "Event vertex Z coord. from TPC tracks", 100, -25., 25.);
 638   } else h->Reset();
 639   fHistos->AddAt(h, kTPCVertex);
 640
 641   // Event vertex
 642   if(!(h = (TH1F*)gROOT->FindObject("hEventVertex"))){
 643     h = new TH1F("hEventVertex", "Event vertex Z coord.", 100, -25., 25.);
 644   } else h->Reset();
 645   fHistos->AddAt(h, kEventVertex);
 646
 647   // Number of all tracks
 648   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAll"))){
 649     h = new TH1I("hNTracksAll", "Number of tracks per event, event vertex cuts", 5000, 0, 5000);
 650   } else h->Reset();
 651   fHistos->AddAt(h, kNTracksAll);
 652
 653   // Number of tracks in acceptance and DCA cut
 654   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAcc"))){
 655     h = new TH1I("hNTracksAcc", Form("Number of tracks per event, |#eta|<%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
 656                                      fgkEta, fgkTrkDCAxy, fgkTrkDCAz), 5000, 0, 5000);
 657   } else h->Reset();
 658   fHistos->AddAt(h, kNTracksAcc);
 659
 660   // Number of tracks in TPC (Ncls>10)
 661   if(!(h = (TH1I*)gROOT->FindObject("hNTracksTPC"))){
 662     h = new TH1I("hNTracksTPC", Form("Number of tracks per event, |#eta|<%.1f, pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 663                                      fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 5000, 0, 5000);
 664   } else h->Reset();
 665   fHistos->AddAt(h, kNTracksTPC);
 666
 667   // Distribution of DCA-xy
 668   if(!(h = (TH1F*)gROOT->FindObject("hDCAxy"))){
 669     h = new TH1F("hDCAxy", "Distribution of transverse DCA", 100, -100., 100.);
 670   } else h->Reset();
 671   fHistos->AddAt(h, kDCAxy);
 672
 673   // Distribution of DCA-z
 674   if(!(h = (TH1F*)gROOT->FindObject("hDCAz"))){
 675     h = new TH1F("hDCAz", "Distribution of longitudinal DCA", 100, -100., 100.);
 676   } else h->Reset();
 677   fHistos->AddAt(h, kDCAz);
 678
 679   Float_t binPtLimits[33] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
 680                              1.0, 1.1, 1.2, 1.3, 1.4,
 681                              1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,
 682                              3.4, 3.8, 4.2, 4.6, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
 683   // Pt distributions
 684   if(!(h = (TH1F*)gROOT->FindObject("hPt1"))){
 685     h = new TH1F("hPt1", Form("dN/dpt, |#eta|<%.1f and pt>%.1f", fgkEta, fgkPt), 32, binPtLimits);
 686   } else h->Reset();
 687   fHistos->AddAt(h, kPt1);
 688
 689   if(!(h = (TH1F*)gROOT->FindObject("hPt2"))){
 690     h = new TH1F("hPt2", Form("dN/dpt, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
 691                               fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 32, binPtLimits);
 692   } else h->Reset();
 693   fHistos->AddAt(h, kPt2);
 694
 695   if(!(h = (TH1F*)gROOT->FindObject("hPt3pos"))){
 696     h = new TH1F("hPt3pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 697                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
 698   } else h->Reset();
 699   fHistos->AddAt(h, kPt3pos);
 700
 701   if(!(h = (TH1F*)gROOT->FindObject("hPt3neg"))){
 702     h = new TH1F("hPt3neg", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 703                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
 704   } else h->Reset();
 705   fHistos->AddAt(h, kPt3neg);
 706
 707   if(!(h = (TH1F*)gROOT->FindObject("hPt4pos"))){
 708     h = new TH1F("hPt4pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
 709                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
 710   } else h->Reset();
 711   fHistos->AddAt(h, kPt4pos);
 712
 713   if(!(h = (TH1F*)gROOT->FindObject("hPt4neg"))){
 714     h = new TH1F("hPt4pos", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
 715                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
 716   } else h->Reset();
 717   fHistos->AddAt(h, kPt4neg);
 718
 719   // theta distribution of TRD tracks
 720   if(!(h = (TH1F*)gROOT->FindObject("hTheta"))){
 721     h = new TH1F("hTheta", Form("dN/d#theta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
 722                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 220,.5,2.7);
 723   } else h->Reset();
 724   fHistos->AddAt(h, kTheta);
 725
 726   // phi distribution of TRD tracks
 727   if(!(h = (TH1F*)gROOT->FindObject("hPhi"))){
 728     h = new TH1F("hPhi", Form("dN/d#varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
 729                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 157,0,6.28);
 730   } else h->Reset();
 731   fHistos->AddAt(h, kPhi);
 732
 733   // TPC cluster distribution
 734   if(!(h = (TH1F*)gROOT->FindObject("hNTPCCl"))){
 735     h = new TH1I("hNTPCCl", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
 736                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
 737   } else h->Reset();
 738   fHistos->AddAt(h, kNTPCCl);
 739
 740   if(!(h = (TH1I*)gROOT->FindObject("hNTPCCl2"))){
 741     h = new TH1F("hNTPCCl2", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, pt>1.0 GeV/c",
 742                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
 743   } else h->Reset();
 744   fHistos->AddAt(h, kNTPCCl2);
 745
 746   // dE/dx vs P for TPC reference tracks
 747   if(!(h = (TH2F*)gROOT->FindObject("hTPCDedx"))){
 748     h = new TH2F("hTPCDedx", Form("TPC dE/dx vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 749                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 120, 0,600.);
 750   } else h->Reset();
 751   fHistos->AddAt(h, kTPCDedx);
 752
 753   // eta,phi distribution of TPC reference tracks
 754   if(!(h = (TH2F*)gROOT->FindObject("hEtaPhi"))){
 755     h = new TH2F("hEtaPhi", Form("TPC (#eta,#varphi), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 756                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 50, -1, 1, 157, 0, 6.28);
 757   } else h->Reset();
 758   fHistos->AddAt(h, kEtaPhi);
 759
 760   // Nclusters vs eta distribution for TPC tracks
 761   if(!(h = (TH2F*)gROOT->FindObject("hEtaNclsTPC"))){
 762     h = new TH2F("hEtaNclsTPC", Form("TPC Nclusters vs. #eta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
 763                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 50, -1, 1, 160, 0, 160.);
 764   } else h->Reset();
 765   fHistos->AddAt(h, kEtaNclsTPC);
 766
 767   // Nclusters vs phi distribution for TPC reference tracks
 768   if(!(h = (TH2F*)gROOT->FindObject("hPhiNclsTPC"))){
 769     h = new TH2F("hPhiNclsTPC", Form("TPC Nclusters vs. #varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
 770                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 157, 0, 6.28, 160, 0, 160.);
 771   } else h->Reset();
 772   fHistos->AddAt(h, kPhiNclsTPC);
 773
 774   // Ntracklets/track vs P for TRD reference tracks
 775   Double_t binsP[19] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.7, 2.0,
 776                         2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 9.0, 12.0};
 777   if(!(h = (TH2F*)gROOT->FindObject("hNTrackletsTRD"))){
 778     h = new TH2F("hNTrackletsTRD", Form("TRD Ntracklets/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 779                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 7, -0.5, 6.5);
 780   } else h->Reset();
 781   fHistos->AddAt(h, kNTrackletsTRD);
 782
 783   // Nclusters/track vs P for TRD reference tracks
 784   if(!(h = (TH2F*)gROOT->FindObject("hNClsTrackTRD"))){
 785     h = new TH2F("hNClsTrackTRD", Form("TRD Nclusters/track vs. P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 786                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 180, 0., 180.);
 787   } else h->Reset();
 788   fHistos->AddAt(h, kNClsTrackTRD);
 789
 790   // <PH> vs slice number for TRD reference tracklets
 791   if(!(h = (TH2F*)gROOT->FindObject("hPHSlice"))){
 792     h = new TH2F("hPHSlice", Form("<PH> vs sliceNo, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 793                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 2000, 0., 2000.);
 794   } else h->Reset();
 795   fHistos->AddAt(h, kPHSlice);
 796
 797   // Qtot vs P for TRD reference tracklets
 798   if(!(h = (TH2F*)gROOT->FindObject("hQtotP"))){
 799     h = new TH2F("hQtotP", Form("Qtot(from slices) vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
 800                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 400, 0., 200);
 801   } else h->Reset();
 802   fHistos->AddAt(h, kQtotP);
 803
 804   // (X,Y,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
 805   if(!(h = (TH3F*)gROOT->FindObject("hPropagXYvsP"))){
 806     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",
 807                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-500,500, 100,-500,500, 10, 0.,10.);
 808   } else h->Reset();
 809   fHistos->AddAt(h, kPropagXYvsP);
 810
 811   // (R,Z,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
 812   if(!(h = (TH3F*)gROOT->FindObject("hPropagRZvsP"))){
 813     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",
 814                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-350., 350., 100,0.,500., 10, 0.,10.);
 815   } else h->Reset();
 816   fHistos->AddAt(h, kPropagRZvsP);
 817
 818   Float_t etaBinLimits[101];
 819   for(Int_t i=0; i<101; i++) etaBinLimits[i] = -1.0 + i*2.0/100.;
 820   Float_t phiBinLimits[151];
 821   for(Int_t i=0; i<151; i++) phiBinLimits[i] = -1.1*TMath::Pi() + i*2.2*TMath::Pi()/150.;
 822   // (eta,detector phi,P) distribution of reference TPC positive tracks
 823   if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksPos"))){
 824     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",
 825                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
 826   } else h->Reset();
 827   fHistos->AddAt(h, kTPCRefTracksPos);
 828
 829   // (eta,detector phi,P) distribution of reference TPC negative tracks
 830   if(!(h = (TH3F*)gROOT->FindObject("hTPCRefTracksNeg"))){
 831     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",
 832                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
 833   } else h->Reset();
 834   fHistos->AddAt(h, kTPCRefTracksNeg);
 835
 836   // (eta,detector phi,P) distribution of reference TRD positive tracks
 837   if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksPos"))){
 838     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",
 839                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
 840   } else h->Reset();
 841   fHistos->AddAt(h, kTRDRefTracksPos);
 842
 843   // (eta,detector phi,P) distribution of reference TRD negative tracks
 844   if(!(h = (TH3F*)gROOT->FindObject("hTRDRefTracksNeg"))){
 845     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",
 846                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
 847   } else h->Reset();
 848   fHistos->AddAt(h, kTRDRefTracksNeg);
 849
 850   // (eta,detector phi) profile of average number of TRD tracklets/track
 851   if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaPhiAvNtrkl"))){
 852     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",
 853                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 150, -1.1*TMath::Pi(), 1.1*TMath::Pi());
 854   } else h->Reset();
 855   fHistos->AddAt(h, kTRDEtaPhiAvNtrkl);
 856
 857   // (eta,delta phi) profile of average number of TRD tracklets/track
 858   if(!(h = (TProfile2D*)gROOT->FindObject("hTRDEtaDeltaPhiAvNtrkl"))){
 859     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",
 860                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -0.4*TMath::Pi(), 0.4*TMath::Pi());
 861   } else h->Reset();
 862   fHistos->AddAt(h, kTRDEtaDeltaPhiAvNtrkl);
 863
 864   return fHistos;
 865 }
 866
 867 //____________________________________________________________________
 868 Bool_t AliTRDcheckESD::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
 869 {
 870 // Load data from performance file
 871
 872   if(!TFile::Open(file)){
 873     AliWarning(Form("Couldn't open file %s.", file));
 874     return kFALSE;
 875   }
 876   if(dir){
 877     if(!gFile->cd(dir)){
 878       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
 879       return kFALSE;
 880     }
 881   }
 882   TObjArray *o(NULL);
 883   const Char_t *tn=(name ? name : GetName());
 884   if(!(o = (TObjArray*)gDirectory->Get(tn))){
 885     AliWarning(Form("Missing histogram container %s.", tn));
 886     return kFALSE;
 887   }
 888   fHistos = (TObjArray*)o->Clone(GetName());
 889   gFile->Close();
 890   return kTRUE;
 891 }
 892
 893 //_______________________________________________________
 894 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
 895 {
 896 // Dump trending value to default file
 897
 898   if(!fgFile){
 899     fgFile = fopen("TRD.Performance.txt", "at");
 900   }
 901   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
 902   return kTRUE;
 903 }
 904
 905 //____________________________________________________________________
 906 void AliTRDcheckESD::Terminate(Option_t *)
 907 {
 908 // Steer post-processing
 909   if(!fHistos){
 910     fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
 911     if(!fHistos){
 912       AliError("Histogram container not found in output");
 913       return;
 914     }
 915   }
 916
 917   const Char_t *name[kNrefs] = {
 918     "Ncl", "Eff", "Eloss", "PtResDCA"
 919   };
 920   TObjArray *arr(NULL); TGraph *g(NULL);
 921   if(!fResults){
 922     fResults = new TObjArray(kNrefs);
 923     fResults->SetOwner();
 924     fResults->SetName("results");
 925     for(Int_t iref(0); iref<kNrefs; iref++){
 926       fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
 927       arr->SetName(name[iref]);  arr->SetOwner();
 928       switch(iref){
 929       case kNCl:
 930         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
 931           arr->AddAt(g = new TGraphErrors(), ig);
 932           g->SetLineColor(ig+1);
 933           g->SetMarkerColor(ig+1);
 934           g->SetMarkerStyle(ig+20);
 935           g->SetName(Form("s%d", ig));
 936           switch(ig){
 937           case 0: g->SetTitle("ALL"); break;
 938           case 1: g->SetTitle("NEG"); break;
 939           case 2: g->SetTitle("POS"); break;
 940           default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;
 941           };
 942         }
 943         break;
 944       case kTRDmom:
 945         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
 946           arr->AddAt(g = new TGraphAsymmErrors(), ig);
 947           g->SetLineColor(ig+1);
 948           g->SetMarkerColor(ig+1);
 949           g->SetMarkerStyle(ig+20);
 950         }
 951         break;
 952       case kPtRes:
 953         for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
 954           Int_t ig(2*idx);
 955           arr->AddAt(g = new TGraphErrors(), ig);
 956           g->SetLineColor(kRed-idx);
 957           g->SetMarkerColor(kRed-idx);
 958           g->SetMarkerStyle(20+idx);
 959           g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
 960           arr->AddAt(g = new TGraphErrors(), ig+1);
 961           g->SetLineColor(kBlue-idx);
 962           g->SetMarkerColor(kBlue-idx);
 963           g->SetMarkerStyle(20+idx);
 964           g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
 965
 966           ig+=10;
 967           arr->AddAt(g = new TGraphErrors(), ig);
 968           g->SetLineColor(kRed-idx);
 969           g->SetMarkerColor(kRed-idx);
 970           g->SetMarkerStyle(20+idx);
 971           g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
 972           arr->AddAt(g = new TGraphErrors(), ig+1);
 973           g->SetLineColor(kBlue-idx);
 974           g->SetMarkerColor(kBlue-idx);
 975           g->SetMarkerStyle(20+idx);
 976           g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
 977         }
 978         break;
 979       default:
 980         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
 981           arr->AddAt(g = new TGraphErrors(), ig);
 982           g->SetLineColor(ig+1);
 983           g->SetMarkerColor(ig+1);
 984           g->SetMarkerStyle(ig+20);
 985         }
 986         break;
 987       }
 988     }
 989   }
 990   TH1 *h1[2] = {NULL, NULL};
 991   TH2I *h2(NULL);
 992   TAxis *ax(NULL);
 993
 994   // No of clusters
 995   if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;
 996   ax = h2->GetXaxis();
 997   arr = (TObjArray*)fResults->At(kNCl);
 998   // All tracks
 999   h1[0] = h2->ProjectionX("Ncl_px");
1000   TGraphErrors *ge=(TGraphErrors*)arr->At(0);
1001   for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1002     ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1003   }
1004   // All charged tracks
1005   TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};
1006   hNclCh[0]->Reset();hNclCh[1]->Reset();
1007   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1008     hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg
1009     hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is));     // pos
1010   }
1011   if(Int_t(hNclCh[0]->GetEntries())){
1012     ge=(TGraphErrors*)arr->At(1);
1013     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1014       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));
1015     }
1016   }
1017   if(Int_t(hNclCh[1]->GetEntries())){
1018     ge=(TGraphErrors*)arr->At(2);
1019     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1020       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));
1021     }
1022   }
1023   // Species wise
1024   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1025     h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);
1026     if(!Int_t(h1[0]->GetEntries())) continue;
1027     ge=(TGraphErrors*)arr->At(2+is);
1028     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1029       ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1030     }
1031   }
1032   fNRefFigures = 1;
1033
1034   // EFFICIENCY
1035   // geometrical efficiency
1036   if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
1037   arr = (TObjArray*)fResults->At(kTRDstat);
1038   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
1039   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
1040   Process(h1, (TGraphErrors*)arr->At(0));
1041   delete h1[0];delete h1[1];
1042   // tracking efficiency
1043   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
1044   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
1045   Process(h1, (TGraphErrors*)arr->At(1));
1046   delete h1[1];
1047   // PID efficiency
1048   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
1049   Process(h1, (TGraphErrors*)arr->At(2));
1050   delete h1[1];
1051   // Refit efficiency
1052   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
1053   Process(h1, (TGraphErrors*)arr->At(3));
1054   delete h1[1];
1055   fNRefFigures++;
1056
1057   // ENERGY LOSS
1058   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
1059   arr = (TObjArray*)fResults->At(kTRDmom);
1060   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
1061   ax=h2->GetXaxis();
1062   const Int_t nq(4);
1063   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
1064   Double_t yq[nq];
1065   for(Int_t ily=6; ily--;){
1066     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
1067     h1[0]->GetQuantiles(nq,yq,xq);
1068     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
1069     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
1070     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
1071     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
1072
1073     //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]);
1074     delete h1[0];
1075   }
1076   fNRefFigures++;
1077   if(!HasMC()) return;
1078
1079   // Pt RESOLUTION @ DCA
1080   TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
1081   if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
1082   arr = (TObjArray*)fResults->At(kPtRes);
1083   TAxis *az(h3->GetZaxis());
1084   for(Int_t i(0); i<AliPID::kSPECIES; i++){
1085     Int_t idx(2*i);
1086     az->SetRange(idx+1, idx+2);
1087     gg[1] = (TGraphErrors*)arr->At(idx);
1088     gg[0] = (TGraphErrors*)arr->At(idx+1);
1089     Process2D((TH2*)h3->Project3D("yx"), gg);
1090
1091     idx+=10;
1092     az->SetRange(idx+1, idx+2);
1093     gg[1] = (TGraphErrors*)arr->At(idx);
1094     gg[0] = (TGraphErrors*)arr->At(idx+1);
1095     Process2D((TH2*)h3->Project3D("yx"), gg);
1096   }
1097   fNRefFigures++;
1098 }
1099
1100 //____________________________________________________________________
1101 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
1102 {
1103   switch(pdg){
1104   case kElectron:
1105   case kPositron: return AliPID::kElectron;
1106   case kMuonPlus:
1107   case kMuonMinus: return AliPID::kMuon;
1108   case kPiPlus:
1109   case kPiMinus: return AliPID::kPion;
1110   case kKPlus:
1111   case kKMinus: return AliPID::kKaon;
1112   case kProton:
1113   case kProtonBar: return AliPID::kProton;
1114   }
1115   return -1;
1116 }
1117
1118 //____________________________________________________________________
1119 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
1120 {
1121 // Generic function to process one reference plot
1122
1123   Int_t n1 = 0, n2 = 0, ip=0;
1124   Double_t eff = 0.;
1125
1126   TAxis *ax = h1[0]->GetXaxis();
1127   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
1128     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
1129     n2 = (Int_t)h1[1]->GetBinContent(ib);
1130     eff = n2/Float_t(n1);
1131
1132     ip=g->GetN();
1133     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
1134     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
1135   }
1136 }
1137 //________________________________________________________
1138 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
1139 {
1140   //
1141   // Do the processing
1142   //
1143
1144   Int_t n = 0;
1145   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1146   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1147   TF1 f("fg", "gaus", -3.,3.);
1148   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1149     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1150     TH1D *h = h2->ProjectionY("py", ibin, ibin);
1151     if(h->GetEntries()<100) continue;
1152     //AdjustF1(h, f);
1153
1154     h->Fit(&f, "QN");
1155     Int_t ip = g[0]->GetN();
1156     g[0]->SetPoint(ip, x, f.GetParameter(1));
1157     g[0]->SetPointError(ip, 0., f.GetParError(1));
1158     g[1]->SetPoint(ip, x, f.GetParameter(2));
1159     g[1]->SetPointError(ip, 0., f.GetParError(2));
1160   }
1161   return;
1162 }
1163 //____________________________________________________________________
1164 void AliTRDcheckESD::PrintStatus(ULong_t status)
1165 {
1166 // Dump track status to stdout
1167
1168   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"
1169     ,Bool_t(status & AliESDtrack::kITSin)
1170     ,Bool_t(status & AliESDtrack::kITSout)
1171     ,Bool_t(status & AliESDtrack::kITSrefit)
1172     ,Bool_t(status & AliESDtrack::kTPCin)
1173     ,Bool_t(status & AliESDtrack::kTPCout)
1174     ,Bool_t(status & AliESDtrack::kTPCrefit)
1175     ,Bool_t(status & AliESDtrack::kTPCpid)
1176     ,Bool_t(status & AliESDtrack::kTRDin)
1177     ,Bool_t(status & AliESDtrack::kTRDout)
1178     ,Bool_t(status & AliESDtrack::kTRDrefit)
1179     ,Bool_t(status & AliESDtrack::kTRDpid)
1180     ,Bool_t(status & AliESDtrack::kTRDStop)
1181     ,Bool_t(status & AliESDtrack::kHMPIDout)
1182     ,Bool_t(status & AliESDtrack::kHMPIDpid)
1183   );
1184 }
1185