TRD/qaRec/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 <TObject.h>
  32 #include <TH2I.h>
  33 #include <TGraphErrors.h>
  34 #include <TGraphAsymmErrors.h>
  35 #include <TFile.h>
  36 #include <TTree.h>
  37 #include <TROOT.h>
  38 #include <TChain.h>
  39 #include <TParticle.h>
  40
  41 #include "AliLog.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDkink.h"
  45 #include "AliMCEvent.h"
  46 #include "AliESDInputHandler.h"
  47 #include "AliMCEventHandler.h"
  48
  49 #include "AliESDtrack.h"
  50 #include "AliMCParticle.h"
  51 #include "AliPID.h"
  52 #include "AliStack.h"
  53 #include "AliTrackReference.h"
  54
  55 #include "AliTRDcheckESD.h"
  56
  57 ClassImp(AliTRDcheckESD)
  58
  59 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
  60 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
  61 FILE* AliTRDcheckESD::fgFile = 0x0;
  62
  63 //____________________________________________________________________
  64 AliTRDcheckESD::AliTRDcheckESD():
  65   AliAnalysisTask("checkESD", "ESD checker for TRD info")
  66   ,fStatus(0)
  67   ,fESD(0x0)
  68   ,fMC(0x0)
  69   ,fHistos(0x0)
  70   ,fResults(0x0)
  71 {
  72   //
  73   // Default constructor
  74   //
  75   SetMC(kTRUE);
  76   DefineInput(0, TChain::Class());
  77   DefineOutput(0, TObjArray::Class());
  78 }
  79
  80 //____________________________________________________________________
  81 AliTRDcheckESD::~AliTRDcheckESD()
  82 {
  83 // Destructor
  84   if(fHistos){
  85     //fHistos->Delete();
  86     delete fHistos;
  87   }
  88   if(fResults){
  89     fResults->Delete();
  90     delete fResults;
  91   }
  92 }
  93
  94 //____________________________________________________________________
  95 void AliTRDcheckESD::ConnectInputData(Option_t *)
  96 {
  97   //
  98   // Link the Input Data
  99   //
 100   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
 101   if(tree) tree->SetBranchStatus("Tracks", 1);
 102
 103   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 104   fESD = esdH ? esdH->GetEvent() : 0x0;
 105
 106   if(!HasMC()) return;
 107   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 108   fMC = mcH ? mcH->MCEvent() : 0x0;
 109 }
 110
 111 //____________________________________________________________________
 112 void AliTRDcheckESD::CreateOutputObjects()
 113 {
 114   //
 115   // Create Output Containers (TObjectArray containing 1D histograms)
 116   //
 117   OpenFile(0, "RECREATE");
 118   Histos();
 119 }
 120
 121 //____________________________________________________________________
 122 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
 123 {
 124   if(!gPad){
 125     AliWarning("Please provide a canvas to draw results.");
 126     return kFALSE;
 127   }
 128   TGraph *g(0x0); TH1 *h(0x0); TLegend *leg(0x0);
 129   switch(ifig){
 130   case kNCl:
 131     if(!fHistos || !(h=(TH1I*)fHistos->At(kNCl))) break;
 132     h->Draw("c");
 133     return kTRUE;
 134   case kTRDstat:
 135     leg=new TLegend(.5, .75, .98, .98);
 136     leg->SetBorderSize(1); leg->SetFillColor(0);
 137     leg->SetHeader("Tracking Efficiency");
 138     for(Int_t ieff=0; ieff<4; ieff++){
 139       if(!(g=GetGraph(ieff, ""))) break;
 140       if(!ieff){
 141         if((h=(TH1S*)gROOT->FindObject("frame"))) delete h;
 142         h=new TH1S("h","",100, 0., 15.);
 143         h->SetLineColor(1);h->SetLineWidth(1);
 144         h->SetMinimum(0.2);h->SetMaximum(1.2);
 145         h->SetXTitle("p_{T} [GeV/c]");
 146         h->SetYTitle("Efficiency");
 147         h->Draw();
 148       }
 149       g->Draw("pc"); leg->AddEntry(g, g->GetTitle(), "pl");
 150     }
 151     leg->Draw();
 152     return kTRUE;
 153   case kTRDmom:
 154     leg=new TLegend(.6, .75, .98, .98);
 155     leg->SetBorderSize(1); leg->SetFillColor(0);
 156     leg->SetHeader("Energy loss");
 157     for(Int_t ieff=0; ieff<2; ieff++){
 158       if(!(g=GetGraph(4+ieff, ""))) break;
 159       if(!ieff){
 160         if((h=(TH1S*)gROOT->FindObject("frame"))) delete h;
 161         h=new TH1S("h","",100, -0.5, 5.5);
 162         h->SetLineColor(0);h->SetLineWidth(1);
 163         h->SetMinimum(-0.5);h->SetMaximum(1.2);
 164         h->SetXTitle("layer");
 165         h->SetYTitle("p_{T}^{Inner TPC} - p_{T}^{layer} [GeV/c]");
 166         h->Draw();
 167       }
 168       g->Draw("p"); leg->AddEntry(g, g->GetTitle(), "pl");
 169     }
 170     leg->Draw();
 171     return kTRUE;
 172   default:
 173     break;
 174   }
 175   AliInfo(Form("Reference plot [%d] missing result", ifig));
 176   return kFALSE;
 177 }
 178
 179 //____________________________________________________________________
 180 TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
 181 {
 182 // Retrieve graph with "id"
 183 // Possible options are :
 184 //   "b" - build graph if none found
 185 //   "c" - clear existing graph
 186
 187   Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
 188          kCLEAR = strstr(opt, "c"); // clear existing graph
 189
 190   const Char_t *name[] = {
 191     "Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
 192   };
 193   const Char_t *title[] = {
 194     "TRD geometrical efficiency (TRDin/TPCout)"
 195     ,"TRD tracking efficiency (TRDout/TRDin)"
 196     ,"TRD PID efficiency (TRDpid/TRDin)"
 197     ,"TRD refit efficiency (TRDrefit/TRDin)"
 198     ,"TRD Eloss (Max/90% quantile)"
 199     ,"TRD Eloss (Mean/60% quantile)"
 200   };
 201   const Int_t ngr = sizeof(name)/sizeof(Char_t*);
 202   if(ngr != kNgraphs){
 203     AliWarning("No of graphs defined different from definition");
 204     return 0x0;
 205   }
 206
 207   if(!fResults){
 208     fResults = new TObjArray(kNgraphs);
 209     fResults->SetOwner();
 210     fResults->SetName("results");
 211   }
 212
 213   TGraph *g = 0x0;
 214   if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
 215     if(kCLEAR){
 216       for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
 217     } else {
 218       PutTrendValue(name[id], g->GetMean(2));
 219       PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 220     }
 221   } else {
 222     if(kBUILD){
 223       switch(id){
 224       case 0:
 225         g = new TGraphErrors();
 226         g->SetMarkerStyle(7);g->SetMarkerColor(kBlack);
 227         g->SetLineColor(kBlack);
 228         break;
 229       case 1:
 230         g = new TGraphErrors();
 231         g->SetMarkerStyle(7);g->SetMarkerColor(kRed);
 232         g->SetLineColor(kRed);
 233         break;
 234       case 2:
 235         g = new TGraphErrors();
 236         g->SetMarkerStyle(7);g->SetMarkerColor(kBlue);
 237         g->SetLineColor(kBlue);
 238         break;
 239       case 3:
 240         g = new TGraphErrors();
 241         g->SetMarkerStyle(7);g->SetMarkerColor(kGreen);
 242         g->SetLineColor(kGreen);
 243         break;
 244       case 4:
 245         g = new TGraphAsymmErrors(6);
 246         g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
 247         g->SetLineColor(kBlack);g->SetLineWidth(2);
 248         break;
 249       case 5:
 250         g = new TGraphAsymmErrors(6);
 251         g->SetMarkerStyle(21);
 252         g->SetLineColor(kRed);g->SetLineWidth(2);
 253         break;
 254       default:
 255         AliWarning(Form("Graph index[%d] missing/not defined.", id));
 256         return 0x0;
 257       }
 258       g->SetNameTitle(name[id], title[id]);
 259       fResults->AddAt(g, id);
 260     }
 261   }
 262   return g;
 263 }
 264
 265 //____________________________________________________________________
 266 void AliTRDcheckESD::Exec(Option_t *){
 267   //
 268   // Run the Analysis
 269   //
 270   if(!fESD){
 271     AliError("ESD event missing.");
 272     return;
 273   }
 274
 275   // Get MC information if available
 276   AliStack * fStack = 0x0;
 277   if(HasMC()){
 278     if(!fMC){
 279       AliWarning("MC event missing");
 280       SetMC(kFALSE);
 281     } else {
 282       if(!(fStack = fMC->Stack())){
 283         AliWarning("MC stack missing");
 284         SetMC(kFALSE);
 285       }
 286     }
 287   }
 288   Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
 289
 290   AliESDtrack *esdTrack = 0x0;
 291   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 292     bTRDin=0;bTRDout=0;bTRDpid=0;
 293     esdTrack = fESD->GetTrack(itrk);
 294
 295 //     if(esdTrack->GetNcls(1)) nTPC++;
 296 //     if(esdTrack->GetNcls(2)) nTRD++;
 297
 298     // track status
 299     ULong_t status = esdTrack->GetStatus();
 300     //PrintStatus(status);
 301
 302     // define TPC out tracks
 303     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
 304     if(esdTrack->GetKinkIndex(0) > 0) continue;
 305
 306     // TRD PID
 307     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 308     // pid quality
 309     //esdTrack->GetTRDntrackletsPID();
 310
 311     // look at external track param
 312     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 313     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
 314     Double_t xyz[3];
 315     if(op){
 316       op->GetXYZ(xyz);
 317       op->Global2LocalPosition(xyz, op->GetAlpha());
 318       //printf("op @ X[%7.3f]\n", xyz[0]);
 319     }
 320
 321     // read MC info
 322     if(!HasMC()) continue;
 323
 324     Int_t fLabel = esdTrack->GetLabel();
 325     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
 326
 327     // read MC particle
 328     AliMCParticle *mcParticle = 0x0;
 329     if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
 330       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 331       continue;
 332     }
 333
 334     AliTrackReference *ref = 0x0;
 335     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 336     if(!nRefs){
 337       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
 338       continue;
 339     }
 340     Int_t iref = 0;
 341     while(iref<nRefs){
 342       ref = mcParticle->GetTrackReference(iref);
 343       if(ref->LocalX() > fgkxTPC) break;
 344       ref=0x0; iref++;
 345     }
 346
 347     // read TParticle
 348     //TParticle *tParticle = mcParticle->Particle();
 349     //Int_t fPdg = tParticle->GetPdgCode();
 350     // reject secondaries
 351     //if(!tParticle->IsPrimary()) continue;
 352
 353     if(ref){
 354       if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 355         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 356       } else {
 357         bTRDin=1;
 358         if(esdTrack->GetNcls(2)) bTRDout=1;
 359         if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
 360       }
 361     } else { // track stopped in TPC
 362       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 363     }
 364     // get the MC pt !!
 365     Float_t pt = ref->Pt();
 366
 367     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
 368     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
 369     if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
 370     if(/*status & AliESDtrack::k*/bTRDout){
 371       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
 372       h->Fill(pt, kTRDout);
 373     }
 374     if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
 375     if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
 376
 377     if(ip){
 378       h = (TH2I*)fHistos->At(kTRDmom);
 379       Float_t pp(0.);
 380       for(Int_t ily=6; ily--;){
 381         if((pp=esdTrack->GetTRDmomentum(ily))<0.) continue;
 382         h->Fill(ip->GetP()-pp, ily);
 383       }
 384     }
 385   }
 386   PostData(0, fHistos);
 387 }
 388
 389 //____________________________________________________________________
 390 TObjArray* AliTRDcheckESD::Histos()
 391 {
 392 // Retrieve histograms array if already build or build it
 393
 394   if(fHistos) return fHistos;
 395
 396   fHistos = new TObjArray(kNhistos);
 397   //fHistos->SetOwner(kTRUE);
 398
 399   TH1 *h = 0x0;
 400
 401   // clusters per tracklet
 402   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
 403     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
 404     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
 405     h->GetYaxis()->SetTitle("entries");
 406   } else h->Reset();
 407   fHistos->AddAt(h, kNCl);
 408
 409   // status bits histogram
 410   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 411     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
 412     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 413     h->GetYaxis()->SetTitle("status bits");
 414     h->GetZaxis()->SetTitle("entries");
 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", 100, -1., 2., 6, -0.5, 5.5);
 421     h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
 422     h->GetYaxis()->SetTitle("layer");
 423     h->GetZaxis()->SetTitle("entries");
 424   } else h->Reset();
 425   fHistos->AddAt(h, kTRDmom);
 426
 427   return fHistos;
 428 }
 429
 430 //____________________________________________________________________
 431 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
 432 {
 433 // Load data from performance file
 434
 435   if(!TFile::Open(filename)){
 436     AliWarning(Form("Couldn't open file %s.", filename));
 437     return kFALSE;
 438   }
 439   TObjArray *o = 0x0;
 440   if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
 441     AliWarning("Missing histogram container.");
 442     return kFALSE;
 443   }
 444   fHistos = (TObjArray*)o->Clone(GetName());
 445   gFile->Close();
 446   return kTRUE;
 447 }
 448
 449 //_______________________________________________________
 450 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
 451 {
 452 // Dump trending value to default file
 453
 454   if(!fgFile){
 455     fgFile = fopen("TRD.Performance.txt", "at");
 456   }
 457   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
 458   return kTRUE;
 459 }
 460
 461 //____________________________________________________________________
 462 void AliTRDcheckESD::Terminate(Option_t *)
 463 {
 464 // Steer post-processing
 465
 466
 467   // geometrical efficiency
 468   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
 469   TH1 *h1[2] = {0x0, 0x0};
 470   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
 471   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
 472   Process(h1, (TGraphErrors*)GetGraph(0));
 473   delete h1[0];delete h1[1];
 474
 475   // tracking efficiency
 476   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
 477   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
 478   Process(h1, (TGraphErrors*)GetGraph(1));
 479   delete h1[1];
 480
 481   // PID efficiency
 482   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
 483   Process(h1, (TGraphErrors*)GetGraph(2));
 484   delete h1[1];
 485
 486   // Refit efficiency
 487   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
 488   Process(h1, (TGraphErrors*)GetGraph(3));
 489   delete h1[1];
 490   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
 491
 492   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)GetGraph(4), *g09 = (TGraphAsymmErrors*)GetGraph(5);
 493   TAxis *ax=h2->GetXaxis();
 494   const Int_t nq(4);
 495   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
 496   Double_t yq[nq];
 497   for(Int_t ily=6; ily--;){
 498     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
 499     h1[0]->GetQuantiles(nq,yq,xq);
 500     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
 501     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
 502     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
 503     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
 504
 505     //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]);
 506     delete h1[0];
 507   }
 508 }
 509
 510 //____________________________________________________________________
 511 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
 512 {
 513 // Generic function to process one reference plot
 514
 515   Int_t n1 = 0, n2 = 0, ip=0;
 516   Double_t eff = 0.;
 517
 518   TAxis *ax = h1[0]->GetXaxis();
 519   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 520     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 521     n2 = (Int_t)h1[1]->GetBinContent(ib);
 522     eff = n2/Float_t(n1);
 523
 524     ip=g->GetN();
 525     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 526     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 527   }
 528 }
 529
 530 //____________________________________________________________________
 531 void AliTRDcheckESD::PrintStatus(ULong_t status)
 532 {
 533 // Dump track status to stdout
 534
 535   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"
 536     ,Bool_t(status & AliESDtrack::kITSin)
 537     ,Bool_t(status & AliESDtrack::kITSout)
 538     ,Bool_t(status & AliESDtrack::kITSrefit)
 539     ,Bool_t(status & AliESDtrack::kTPCin)
 540     ,Bool_t(status & AliESDtrack::kTPCout)
 541     ,Bool_t(status & AliESDtrack::kTPCrefit)
 542     ,Bool_t(status & AliESDtrack::kTPCpid)
 543     ,Bool_t(status & AliESDtrack::kTRDin)
 544     ,Bool_t(status & AliESDtrack::kTRDout)
 545     ,Bool_t(status & AliESDtrack::kTRDrefit)
 546     ,Bool_t(status & AliESDtrack::kTRDpid)
 547     ,Bool_t(status & AliESDtrack::kTRDStop)
 548     ,Bool_t(status & AliESDtrack::kHMPIDout)
 549     ,Bool_t(status & AliESDtrack::kHMPIDpid)
 550   );
 551 }
 552