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 <TH3S.h>
  36 #include <TGraphErrors.h>
  37 #include <TGraphAsymmErrors.h>
  38 #include <TFile.h>
  39 #include <TTree.h>
  40 #include <TROOT.h>
  41 #include <TChain.h>
  42 #include <TParticle.h>
  43
  44 #include "AliLog.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDkink.h"
  48 #include "AliMCEvent.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliMCEventHandler.h"
  51
  52 #include "AliESDtrack.h"
  53 #include "AliMCParticle.h"
  54 #include "AliPID.h"
  55 #include "AliStack.h"
  56 #include "AliTrackReference.h"
  57
  58 #include "AliTRDcheckESD.h"
  59
  60 ClassImp(AliTRDcheckESD)
  61
  62 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
  63 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
  64 const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
  65 1, 4, 2, 20};
  66 FILE* AliTRDcheckESD::fgFile = NULL;
  67
  68 //____________________________________________________________________
  69 AliTRDcheckESD::AliTRDcheckESD():
  70   AliAnalysisTaskSE()
  71   ,fStatus(0)
  72   ,fNRefFigures(0)
  73   ,fESD(NULL)
  74   ,fMC(NULL)
  75   ,fHistos(NULL)
  76   ,fResults(NULL)
  77 {
  78   //
  79   // Default constructor
  80   //
  81   SetNameTitle("checkESD", "Check TRD @ ESD level");
  82   SetMC(kTRUE);
  83 }
  84
  85 //____________________________________________________________________
  86 AliTRDcheckESD::AliTRDcheckESD(char* name):
  87   AliAnalysisTaskSE(name)
  88   ,fStatus(0)
  89   ,fNRefFigures(0)
  90   ,fESD(NULL)
  91   ,fMC(NULL)
  92   ,fHistos(NULL)
  93   ,fResults(NULL)
  94 {
  95   //
  96   // Default constructor
  97   //
  98   SetMC(kTRUE);
  99   SetTitle("Check TRD @ ESD level");
 100   DefineOutput(1, TObjArray::Class());
 101 }
 102
 103 //____________________________________________________________________
 104 AliTRDcheckESD::~AliTRDcheckESD()
 105 {
 106 // Destructor
 107   if(fHistos){
 108     //fHistos->Delete();
 109     delete fHistos;
 110   }
 111   if(fResults){
 112     fResults->Delete();
 113     delete fResults;
 114   }
 115 }
 116
 117 //____________________________________________________________________
 118 void AliTRDcheckESD::UserCreateOutputObjects()
 119 {
 120   //
 121   // Create Output Containers (TObjectArray containing 1D histograms)
 122   //
 123   //OpenFile(0, "RECREATE");
 124
 125   Histos();
 126 }
 127
 128
 129 //____________________________________________________________________
 130 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
 131 {
 132   if(ifig>=fNRefFigures){
 133     AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
 134     return kFALSE;
 135   }
 136   if(!gPad){
 137     AliWarning("Please provide a canvas to draw results.");
 138     return kFALSE;
 139   }
 140
 141   const Char_t *title[20];
 142   TH1 *hF(NULL);
 143   if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
 144   TLegend *leg(NULL);
 145   TList *l(NULL); TVirtualPad *pad(NULL);
 146   TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
 147   TObjArray *arr(NULL);
 148   switch(ifig){
 149   case kNCl: // number of clusters/track
 150     if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
 151     g=(TGraphErrors*)arr->At(0);
 152     g->Draw("apc");
 153     hF=g->GetHistogram();
 154     hF->SetXTitle("no of clusters");
 155     hF->SetYTitle("entries");
 156     break;
 157   case kTRDstat: // Efficiency
 158     if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
 159     leg = new TLegend(.65, .7, .95, .99);
 160     leg->SetHeader("TRD Efficiency");
 161     leg->SetBorderSize(1); leg->SetFillColor(0);
 162     title[0] = "Geometrical (TRDin/TPCout)";
 163     title[1] = "Tracking (TRDout/TRDin)";
 164     title[2] = "PID (TRDpid/TRDin)";
 165     title[3] = "Refit (TRDrefit/TRDin)";
 166     hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
 167     hF->SetMaximum(1.3);
 168     hF->GetXaxis()->SetMoreLogLabels();
 169     hF->Draw("p");
 170     for(Int_t ig(0); ig<fgkNgraph[kTRDstat]; ig++){
 171       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
 172       g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
 173       //PutTrendValue(name[id], g->GetMean(2));
 174       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 175     }
 176     leg->Draw(); gPad->SetLogx();
 177     break;
 178   case kTRDmom: // Energy loss
 179     if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
 180     leg = new TLegend(.65, .7, .95, .99);
 181     leg->SetHeader("Energy Loss");
 182     leg->SetBorderSize(1); leg->SetFillColor(0);
 183     title[0] = "Max & 90% quantile";
 184     title[1] = "Mean & 60% quantile";
 185     hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
 186     hF->SetMaximum(1.3);hF->SetMinimum(-.3);
 187     hF->Draw("p");
 188     for(Int_t ig(0); ig<fgkNgraph[kTRDmom]; ig++){
 189       if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
 190       ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
 191       //PutTrendValue(name[id], g->GetMean(2));
 192       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 193     }
 194     leg->Draw();gPad->SetLogx(kFALSE);
 195     break;
 196   case kPtRes: // Pt resolution @ vertex
 197     if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
 198     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 199     pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();
 200     pad->SetMargin(0.1, 0.022, 0.1, 0.023);
 201     hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
 202     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 203     hF->GetXaxis()->SetMoreLogLabels();
 204     hF->GetXaxis()->SetTitleOffset(1.2);
 205     hF->GetYaxis()->CenterTitle();
 206     hF->Draw("p");
 207     for(Int_t ig(0); ig<fgkNgraph[kPtRes]/2; ig++){
 208       if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 209       if(!g->GetN()) continue;
 210       g->Draw("pl");
 211       //PutTrendValue(name[id], g->GetMean(2));
 212       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 213     }
 214     pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();
 215     pad->SetMargin(0.1, 0.22, 0.1, 0.023);
 216     hF = (TH1*)hF->Clone("hFcheckESD1");
 217     hF->SetTitle("ITS+TPC");
 218     hF->SetMaximum(10.);hF->SetMinimum(-3.);
 219     hF->Draw("p");
 220     leg = new TLegend(.78, .1, .99, .98);
 221     leg->SetHeader("P_{t} @ DCA");
 222     leg->SetBorderSize(1); leg->SetFillColor(0);
 223     leg->SetTextAlign(22);
 224     leg->SetTextFont(12);
 225     leg->SetTextSize(0.03813559);
 226     {
 227       Int_t nPlots(0);
 228       for(Int_t ig(fgkNgraph[kPtRes]/2); ig<fgkNgraph[kPtRes]; ig++){
 229         if(!(g = (TGraphErrors*)arr->At(ig))) continue;
 230         if(!g->GetN()) continue;
 231         nPlots++;
 232         g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
 233         //PutTrendValue(name[id], g->GetMean(2));
 234         //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 235       }
 236       if(nPlots) leg->Draw();
 237     }
 238
 239     break;
 240   }
 241   return kTRUE;
 242 }
 243
 244 //____________________________________________________________________
 245 void AliTRDcheckESD::UserExec(Option_t *){
 246   //
 247   // Run the Analysis
 248   //
 249   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 250   fMC = MCEvent();
 251
 252   if(!fESD){
 253     AliError("ESD event missing.");
 254     return;
 255   }
 256
 257   // Get MC information if available
 258   AliStack * fStack = NULL;
 259   if(HasMC()){
 260     if(!fMC){
 261       AliWarning("MC event missing");
 262       SetMC(kFALSE);
 263     } else {
 264       if(!(fStack = fMC->Stack())){
 265         AliWarning("MC stack missing");
 266         SetMC(kFALSE);
 267       }
 268     }
 269   }
 270   TH2 *h(NULL);
 271
 272   AliESDtrack *esdTrack(NULL);
 273   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 274     esdTrack = fESD->GetTrack(itrk);
 275
 276     // track status
 277     ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
 278     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
 279     if(esdTrack->GetKinkIndex(0) > 0) continue;
 280
 281     //Int_t nTPC(esdTrack->GetNcls(1));
 282     Int_t nTRD(esdTrack->GetNcls(2));
 283     Double_t pt(esdTrack->Pt());
 284     //Double_t eta(esdTrack->Eta());
 285     //Double_t phi(esdTrack->Phi());
 286     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 287     // pid quality
 288     //esdTrack->GetTRDntrackletsPID();
 289     Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
 290
 291     // look at external track param
 292     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 293     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
 294
 295     Double_t pt0(0.), eta0(0.), phi0(0.), ptTRD(0.);
 296     // read MC info if available
 297     Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
 298     AliMCParticle *mcParticle(NULL);
 299     if(HasMC()){
 300       AliTrackReference *ref(NULL);
 301       Int_t fLabel(esdTrack->GetLabel());
 302       Int_t fIdx(TMath::Abs(fLabel));
 303       if(fIdx > fStack->GetNtrack()) continue;
 304
 305       // read MC particle
 306       if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
 307         AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 308         continue;
 309       }
 310       pt0  = mcParticle->Pt();
 311       eta0 = mcParticle->Eta();
 312       phi0 = mcParticle->Phi();
 313       kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
 314
 315       // read track references
 316       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 317       if(!nRefs){
 318         AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
 319         continue;
 320       }
 321       Int_t iref = 0;
 322       while(iref<nRefs){
 323         ref = mcParticle->GetTrackReference(iref);
 324         if(ref->LocalX() > fgkxTPC) break;
 325         ref=NULL; iref++;
 326       }
 327       if(ref){
 328         if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 329           ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 330         }
 331       } else { // track stopped in TPC
 332         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 333       }
 334       ptTRD = ref->Pt();kFOUND=kTRUE;
 335     } else { // use reconstructed values
 336       if(op){
 337         Double_t x(op->GetX());
 338         if(x<fgkxTOF && x>fgkxTPC){
 339           ptTRD=op->Pt();
 340           kFOUND=kTRUE;
 341         }
 342       }
 343
 344       if(!kFOUND && ip){
 345         ptTRD=ip->Pt();
 346         kFOUND=kTRUE;
 347       }
 348     }
 349
 350     if(kFOUND){
 351       h = (TH2I*)fHistos->At(kTRDstat);
 352       if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
 353       if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
 354       if(kBarrel && (status & AliESDtrack::kTRDout)){
 355         ((TH1*)fHistos->At(kNCl))->Fill(nTRD);
 356         h->Fill(ptTRD, kTRDout);
 357       }
 358       if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
 359       if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
 360     }
 361     if(HasMC() &&
 362       kBarrel && kPhysPrim &&
 363       TMath::Abs(eta0) < 0.9) {
 364       TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
 365       Int_t sgn = mcParticle->Charge()<0?0:1;
 366       Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10;
 367
 368       h3->Fill(pt0, 1.e2*(pt/pt0-1.),
 369         offset + 2*Pdg2Idx(TMath::Abs(mcParticle->PdgCode())) + sgn);
 370     }
 371     if(ip){
 372       h = (TH2I*)fHistos->At(kTRDmom);
 373       Float_t pTRD(0.);
 374       for(Int_t ily=6; ily--;){
 375         if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
 376         h->Fill(ip->GetP()-pTRD, ily);
 377       }
 378     }
 379   }
 380   PostData(1, fHistos);
 381 }
 382
 383 //____________________________________________________________________
 384 TObjArray* AliTRDcheckESD::Histos()
 385 {
 386 // Retrieve histograms array if already build or build it
 387
 388   if(fHistos) return fHistos;
 389
 390   fHistos = new TObjArray(kNhistos);
 391   //fHistos->SetOwner(kTRUE);
 392
 393   TH1 *h = NULL;
 394
 395   // clusters per tracklet
 396   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
 397     h = new TH1I("hNCl", "Clusters per TRD track;N_{cl}^{TRD};entries", 100, 0., 200.);
 398   } else h->Reset();
 399   fHistos->AddAt(h, kNCl); fNRefFigures++;
 400
 401   // status bits histogram
 402   const Int_t kNpt(30), kNbits(5);
 403   Float_t Pt(0.2), Bits(.5);
 404   Float_t binsPt[kNpt+1], binsBits[kNbits+1];
 405   for(Int_t i=0;i<kNpt+1; i++,Pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=Pt;
 406   for(Int_t i=0; i<kNbits+1; i++,Bits+=1.) binsBits[i]=Bits;
 407   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 408     h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
 409     TAxis *ay(h->GetYaxis());
 410     ay->SetBinLabel(1, "kTPCout");
 411     ay->SetBinLabel(2, "kTRDin");
 412     ay->SetBinLabel(3, "kTRDout");
 413     ay->SetBinLabel(4, "kTRDpid");
 414     ay->SetBinLabel(5, "kTRDrefit");
 415   } else h->Reset();
 416   fHistos->AddAt(h, kTRDstat);
 417
 418   // energy loss
 419   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
 420     h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
 421   } else h->Reset();
 422   fHistos->AddAt(h, kTRDmom);
 423   if(!HasMC()) return fHistos;
 424
 425   // pt resolution
 426   const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
 427   Float_t DPt(-3.), Spec(-0.5);
 428   Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
 429   for(Int_t i=0; i<kNdpt+1; i++,DPt+=6.e-2) binsDPt[i]=DPt;
 430   for(Int_t i=0; i<kNspec+1; i++,Spec+=1.) binsSpec[i]=Spec;
 431   if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
 432     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);
 433     TAxis *az(h->GetZaxis());
 434     az->SetLabelOffset(0.015);
 435     for(Int_t i(0); i<AliPID::kSPECIES; i++){
 436       az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
 437       az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
 438       az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
 439       az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
 440     }
 441   } else h->Reset();
 442   fHistos->AddAt(h, kPtRes);
 443
 444   return fHistos;
 445 }
 446
 447 //____________________________________________________________________
 448 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
 449 {
 450 // Load data from performance file
 451
 452   if(!TFile::Open(filename)){
 453     AliWarning(Form("Couldn't open file %s.", filename));
 454     return kFALSE;
 455   }
 456   TObjArray *o = NULL;
 457   if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
 458     AliWarning("Missing histogram container.");
 459     return kFALSE;
 460   }
 461   fHistos = (TObjArray*)o->Clone(GetName());
 462   gFile->Close();
 463   SETBIT(fStatus, kLoad);
 464   return kTRUE;
 465 }
 466
 467 //_______________________________________________________
 468 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
 469 {
 470 // Dump trending value to default file
 471
 472   if(!fgFile){
 473     fgFile = fopen("TRD.Performance.txt", "at");
 474   }
 475   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
 476   return kTRUE;
 477 }
 478
 479 //____________________________________________________________________
 480 void AliTRDcheckESD::Terminate(Option_t *)
 481 {
 482 // Steer post-processing
 483   if(!IsLoad()){
 484     fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
 485     if(!fHistos){
 486       AliError("Histogram container not found in output");
 487       return;
 488     }
 489   }
 490
 491   const Char_t *name[kNrefs] = {
 492     "Ncl", "Eff", "Eloss", "PtResDCA"
 493   };
 494   TObjArray *arr(NULL); TGraph *g(NULL);
 495   if(!fResults){
 496     fResults = new TObjArray(kNrefs);
 497     fResults->SetOwner();
 498     fResults->SetName("results");
 499     for(Int_t iref(0); iref<kNrefs; iref++){
 500       fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
 501       arr->SetName(name[iref]);  arr->SetOwner();
 502       switch(iref){
 503       case kTRDmom:
 504         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
 505           arr->AddAt(g = new TGraphAsymmErrors(), ig);
 506           g->SetLineColor(ig+1);
 507           g->SetMarkerColor(ig+1);
 508           g->SetMarkerStyle(ig+20);
 509         }
 510         break;
 511       case kPtRes:
 512         for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
 513           Int_t ig(2*idx);
 514           arr->AddAt(g = new TGraphErrors(), ig);
 515           g->SetLineColor(kRed-idx);
 516           g->SetMarkerColor(kRed-idx);
 517           g->SetMarkerStyle(20+idx);
 518           g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
 519           arr->AddAt(g = new TGraphErrors(), ig+1);
 520           g->SetLineColor(kBlue-idx);
 521           g->SetMarkerColor(kBlue-idx);
 522           g->SetMarkerStyle(20+idx);
 523           g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
 524
 525           ig+=10;
 526           arr->AddAt(g = new TGraphErrors(), ig);
 527           g->SetLineColor(kRed-idx);
 528           g->SetMarkerColor(kRed-idx);
 529           g->SetMarkerStyle(20+idx);
 530           g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
 531           arr->AddAt(g = new TGraphErrors(), ig+1);
 532           g->SetLineColor(kBlue-idx);
 533           g->SetMarkerColor(kBlue-idx);
 534           g->SetMarkerStyle(20+idx);
 535           g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
 536         }
 537         break;
 538       default:
 539         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
 540           arr->AddAt(g = new TGraphErrors(), ig);
 541           g->SetLineColor(ig+1);
 542           g->SetMarkerColor(ig+1);
 543           g->SetMarkerStyle(ig+20);
 544         }
 545         break;
 546       }
 547     }
 548   }
 549   TH1 *h1[2] = {NULL, NULL};
 550   TH2I *h2(NULL);
 551   TAxis *ax(NULL);
 552
 553   // No of clusters
 554   if(!(h1[0] = (TH1I*)fHistos->At(kNCl))) return;
 555   arr = (TObjArray*)fResults->At(kNCl);
 556   TGraphErrors *ge=(TGraphErrors*)arr->At(0);
 557   ax = h1[0]->GetXaxis();
 558   for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
 559     ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
 560   }
 561   fNRefFigures = 1;
 562
 563   // EFFICIENCY
 564   // geometrical efficiency
 565   if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
 566   arr = (TObjArray*)fResults->At(kTRDstat);
 567   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
 568   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
 569   Process(h1, (TGraphErrors*)arr->At(0));
 570   delete h1[0];delete h1[1];
 571   // tracking efficiency
 572   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
 573   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
 574   Process(h1, (TGraphErrors*)arr->At(1));
 575   delete h1[1];
 576   // PID efficiency
 577   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
 578   Process(h1, (TGraphErrors*)arr->At(2));
 579   delete h1[1];
 580   // Refit efficiency
 581   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
 582   Process(h1, (TGraphErrors*)arr->At(3));
 583   delete h1[1];
 584   fNRefFigures++;
 585
 586   // ENERGY LOSS
 587   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
 588   arr = (TObjArray*)fResults->At(kTRDmom);
 589   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
 590   ax=h2->GetXaxis();
 591   const Int_t nq(4);
 592   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
 593   Double_t yq[nq];
 594   for(Int_t ily=6; ily--;){
 595     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
 596     h1[0]->GetQuantiles(nq,yq,xq);
 597     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
 598     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
 599     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
 600     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
 601
 602     //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]);
 603     delete h1[0];
 604   }
 605   fNRefFigures++;
 606   if(!HasMC()) return;
 607
 608   // Pt RESOLUTION @ DCA
 609   TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
 610   if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
 611   arr = (TObjArray*)fResults->At(kPtRes);
 612   TAxis *az(h3->GetZaxis());
 613   for(Int_t i(0); i<AliPID::kSPECIES; i++){
 614     Int_t idx(2*i);
 615     az->SetRange(idx+1, idx+2);
 616     gg[1] = (TGraphErrors*)arr->At(idx);
 617     gg[0] = (TGraphErrors*)arr->At(idx+1);
 618     Process2D((TH2*)h3->Project3D("yx"), gg);
 619
 620     idx+=10;
 621     az->SetRange(idx+1, idx+2);
 622     gg[1] = (TGraphErrors*)arr->At(idx);
 623     gg[0] = (TGraphErrors*)arr->At(idx+1);
 624     Process2D((TH2*)h3->Project3D("yx"), gg);
 625   }
 626   fNRefFigures++;
 627 }
 628
 629 //____________________________________________________________________
 630 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg)
 631 {
 632   switch(pdg){
 633   case kElectron: return AliPID::kElectron;
 634   case kMuonMinus: return AliPID::kMuon;
 635   case kPiPlus: return AliPID::kPion;
 636   case kKPlus: return AliPID::kKaon;
 637   case kProton: return AliPID::kProton;
 638   }
 639   return -1;
 640 }
 641
 642 //____________________________________________________________________
 643 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
 644 {
 645 // Generic function to process one reference plot
 646
 647   Int_t n1 = 0, n2 = 0, ip=0;
 648   Double_t eff = 0.;
 649
 650   TAxis *ax = h1[0]->GetXaxis();
 651   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 652     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 653     n2 = (Int_t)h1[1]->GetBinContent(ib);
 654     eff = n2/Float_t(n1);
 655
 656     ip=g->GetN();
 657     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 658     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 659   }
 660 }
 661 //________________________________________________________
 662 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
 663 {
 664   //
 665   // Do the processing
 666   //
 667
 668   Int_t n = 0;
 669   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
 670   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
 671   TF1 f("fg", "gaus", -3.,3.);
 672   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 673     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
 674     TH1D *h = h2->ProjectionY("py", ibin, ibin);
 675     if(h->GetEntries()<100) continue;
 676     //AdjustF1(h, f);
 677
 678     h->Fit(&f, "QN");
 679     Int_t ip = g[0]->GetN();
 680     g[0]->SetPoint(ip, x, f.GetParameter(1));
 681     g[0]->SetPointError(ip, 0., f.GetParError(1));
 682     g[1]->SetPoint(ip, x, f.GetParameter(2));
 683     g[1]->SetPointError(ip, 0., f.GetParError(2));
 684   }
 685   return;
 686 }
 687 //____________________________________________________________________
 688 void AliTRDcheckESD::PrintStatus(ULong_t status)
 689 {
 690 // Dump track status to stdout
 691
 692   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"
 693     ,Bool_t(status & AliESDtrack::kITSin)
 694     ,Bool_t(status & AliESDtrack::kITSout)
 695     ,Bool_t(status & AliESDtrack::kITSrefit)
 696     ,Bool_t(status & AliESDtrack::kTPCin)
 697     ,Bool_t(status & AliESDtrack::kTPCout)
 698     ,Bool_t(status & AliESDtrack::kTPCrefit)
 699     ,Bool_t(status & AliESDtrack::kTPCpid)
 700     ,Bool_t(status & AliESDtrack::kTRDin)
 701     ,Bool_t(status & AliESDtrack::kTRDout)
 702     ,Bool_t(status & AliESDtrack::kTRDrefit)
 703     ,Bool_t(status & AliESDtrack::kTRDpid)
 704     ,Bool_t(status & AliESDtrack::kTRDStop)
 705     ,Bool_t(status & AliESDtrack::kHMPIDout)
 706     ,Bool_t(status & AliESDtrack::kHMPIDpid)
 707   );
 708 }
 709