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