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