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 <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 = NULL;
  62
  63 //____________________________________________________________________
  64 AliTRDcheckESD::AliTRDcheckESD():
  65   AliAnalysisTask("checkESD", "ESD checker for TRD info")
  66   ,fStatus(0)
  67   ,fESD(NULL)
  68   ,fMC(NULL)
  69   ,fHistos(NULL)
  70   ,fResults(NULL)
  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() : NULL;
 105
 106   if(!HasMC()) return;
 107   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 108   fMC = mcH ? mcH->MCEvent() : NULL;
 109 }
 110
 111 //____________________________________________________________________
 112 void AliTRDcheckESD::CreateOutputObjects()
 113 {
 114   //
 115   // Create Output Containers (TObjectArray containing 1D histograms)
 116   //
 117   //OpenFile(0, "RECREATE");
 118
 119   Histos();
 120 }
 121
 122 //____________________________________________________________________
 123 TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
 124 {
 125 // Retrieve graph with "id"
 126 // Possible options are :
 127 //   "b" - build graph if none found
 128 //   "c" - clear existing graph
 129
 130   Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
 131          kCLEAR = strstr(opt, "c"); // clear existing graph
 132
 133   const Char_t *name[] = {
 134     "Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
 135   };
 136   const Char_t *title[] = {
 137     "TRD geometrical efficiency (TRDin/TPCout)"
 138     ,"TRD tracking efficiency (TRDout/TRDin)"
 139     ,"TRD PID efficiency (TRDpid/TRDin)"
 140     ,"TRD refit efficiency (TRDrefit/TRDin)"
 141     ,"TRD Eloss (Max/90% quantile)"
 142     ,"TRD Eloss (Mean/60% quantile)"
 143   };
 144   const Int_t ngr = sizeof(name)/sizeof(Char_t*);
 145   if(ngr != kNgraphs){
 146     AliWarning("No of graphs defined different from definition");
 147     return NULL;
 148   }
 149
 150   if(!fResults){
 151     fResults = new TObjArray(kNgraphs);
 152     fResults->SetOwner();
 153     fResults->SetName("results");
 154   }
 155
 156   TGraph *g = NULL;
 157   if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
 158     if(kCLEAR){
 159       for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
 160     } else {
 161       PutTrendValue(name[id], g->GetMean(2));
 162       PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
 163     }
 164   } else {
 165     if(kBUILD){
 166       switch(id){
 167       case 0:
 168         g = new TGraphErrors();
 169         break;
 170       case 1:
 171         g = new TGraphErrors();
 172         break;
 173       case 2:
 174         g = new TGraphErrors();
 175         break;
 176       case 3:
 177         g = new TGraphErrors();
 178         break;
 179       case 4:
 180         g = new TGraphAsymmErrors(6);
 181         g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
 182         g->SetLineColor(kBlack);g->SetLineWidth(2);
 183         break;
 184       case 5:
 185         g = new TGraphAsymmErrors(6);
 186         g->SetMarkerStyle(21);
 187         g->SetLineColor(kRed);g->SetLineWidth(2);
 188         break;
 189       default:
 190         AliWarning(Form("Graph index[%d] missing/not defined.", id));
 191         return NULL;
 192       }
 193       g->SetNameTitle(name[id], title[id]);
 194       fResults->AddAt(g, id);
 195     }
 196   }
 197   return g;
 198 }
 199
 200 //____________________________________________________________________
 201 void AliTRDcheckESD::Exec(Option_t *){
 202   //
 203   // Run the Analysis
 204   //
 205   if(!fESD){
 206     AliError("ESD event missing.");
 207     return;
 208   }
 209
 210   // Get MC information if available
 211   AliStack * fStack = NULL;
 212   if(HasMC()){
 213     if(!fMC){
 214       AliWarning("MC event missing");
 215       SetMC(kFALSE);
 216     } else {
 217       if(!(fStack = fMC->Stack())){
 218         AliWarning("MC stack missing");
 219         SetMC(kFALSE);
 220       }
 221     }
 222   }
 223   TH2 *h(NULL);
 224
 225   AliESDtrack *esdTrack(NULL);
 226   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 227     esdTrack = fESD->GetTrack(itrk);
 228
 229 //     if(esdTrack->GetNcls(1)) nTPC++;
 230 //     if(esdTrack->GetNcls(2)) nTRD++;
 231
 232     // track status
 233     ULong_t status = esdTrack->GetStatus();
 234     //PrintStatus(status);
 235
 236     // define TPC out tracks
 237     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
 238     if(esdTrack->GetKinkIndex(0) > 0) continue;
 239
 240     // TRD PID
 241     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 242     // pid quality
 243     //esdTrack->GetTRDntrackletsPID();
 244
 245     // look at external track param
 246     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 247     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
 248
 249     Float_t pt(0.); Bool_t kFOUND(kFALSE);
 250
 251     // read MC info if available
 252     if(HasMC()){
 253       AliTrackReference *ref(NULL);
 254       Int_t fLabel(esdTrack->GetLabel());
 255       if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
 256
 257       // read MC particle
 258       AliMCParticle *mcParticle(NULL);
 259       if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
 260         AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 261         continue;
 262       }
 263
 264       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 265       if(!nRefs){
 266         AliWarning(Form("Track refs missing. Label[%d].", fLabel));
 267         continue;
 268       }
 269       Int_t iref = 0;
 270       while(iref<nRefs){
 271         ref = mcParticle->GetTrackReference(iref);
 272         if(ref->LocalX() > fgkxTPC) break;
 273         ref=NULL; iref++;
 274       }
 275       if(ref){
 276         if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 277           ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 278         }
 279       } else { // track stopped in TPC
 280         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 281       }
 282       pt = ref->Pt();kFOUND=kTRUE;
 283     } else { // use reconstructed values
 284       if(op){
 285         Double_t x(op->GetX());
 286         if(x<fgkxTOF && x>fgkxTPC){
 287           pt=op->Pt();
 288           kFOUND=kTRUE;
 289         }
 290       }
 291
 292       if(!kFOUND && ip){
 293         pt=ip->Pt();
 294         kFOUND=kTRUE;
 295       }
 296     }
 297
 298     if(kFOUND){
 299       h = (TH2I*)fHistos->At(kTRDstat);
 300       if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
 301       if(status & AliESDtrack::kTRDin) h->Fill(pt, kTRDin);
 302       if(status & AliESDtrack::kTRDout){
 303         ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
 304         h->Fill(pt, kTRDout);
 305       }
 306       if(status & AliESDtrack::kTRDpid) h->Fill(pt, kTRDpid);
 307       if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
 308     }
 309
 310     if(ip){
 311       h = (TH2I*)fHistos->At(kTRDmom);
 312       Float_t p(0.);
 313       for(Int_t ily=6; ily--;){
 314         if((p=esdTrack->GetTRDmomentum(ily))<0.) continue;
 315         h->Fill(ip->GetP()-p, ily);
 316       }
 317     }
 318   }
 319   PostData(0, fHistos);
 320 }
 321
 322 //____________________________________________________________________
 323 TObjArray* AliTRDcheckESD::Histos()
 324 {
 325 // Retrieve histograms array if already build or build it
 326
 327   if(fHistos) return fHistos;
 328
 329   fHistos = new TObjArray(kNhistos);
 330   //fHistos->SetOwner(kTRUE);
 331
 332   TH1 *h = NULL;
 333
 334   // clusters per tracklet
 335   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
 336     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
 337     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
 338     h->GetYaxis()->SetTitle("entries");
 339   } else h->Reset();
 340   fHistos->AddAt(h, kNCl);
 341
 342   // status bits histogram
 343   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 344     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
 345     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 346     h->GetYaxis()->SetTitle("status bits");
 347     h->GetZaxis()->SetTitle("entries");
 348   } else h->Reset();
 349   fHistos->AddAt(h, kTRDstat);
 350
 351   // energy loss
 352   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
 353     h = new TH2I("hTRDmom", "TRD energy loss", 100, -1., 2., 6, -0.5, 5.5);
 354     h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
 355     h->GetYaxis()->SetTitle("layer");
 356     h->GetZaxis()->SetTitle("entries");
 357   } else h->Reset();
 358   fHistos->AddAt(h, kTRDmom);
 359
 360   return fHistos;
 361 }
 362
 363 //____________________________________________________________________
 364 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
 365 {
 366 // Load data from performance file
 367
 368   if(!TFile::Open(filename)){
 369     AliWarning(Form("Couldn't open file %s.", filename));
 370     return kFALSE;
 371   }
 372   TObjArray *o = NULL;
 373   if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
 374     AliWarning("Missing histogram container.");
 375     return kFALSE;
 376   }
 377   fHistos = (TObjArray*)o->Clone(GetName());
 378   gFile->Close();
 379   SETBIT(fStatus, kLoad);
 380   return kTRUE;
 381 }
 382
 383 //_______________________________________________________
 384 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
 385 {
 386 // Dump trending value to default file
 387
 388   if(!fgFile){
 389     fgFile = fopen("TRD.Performance.txt", "at");
 390   }
 391   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
 392   return kTRUE;
 393 }
 394
 395 //____________________________________________________________________
 396 void AliTRDcheckESD::Terminate(Option_t *)
 397 {
 398 // Steer post-processing
 399   if(!IsLoad()){
 400     fHistos = dynamic_cast<TObjArray *>(GetOutputData(0));
 401     if(!fHistos){
 402       AliError("Histogram container not found in output");
 403       return;
 404     }
 405   }
 406
 407   // geometrical efficiency
 408   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
 409   TH1 *h1[2] = {NULL, NULL};
 410   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
 411   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
 412   Process(h1, (TGraphErrors*)GetGraph(0));
 413   delete h1[0];delete h1[1];
 414
 415   // tracking efficiency
 416   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
 417   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
 418   Process(h1, (TGraphErrors*)GetGraph(1));
 419   delete h1[1];
 420
 421   // PID efficiency
 422   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
 423   Process(h1, (TGraphErrors*)GetGraph(2));
 424   delete h1[1];
 425
 426   // Refit efficiency
 427   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
 428   Process(h1, (TGraphErrors*)GetGraph(3));
 429   delete h1[1];
 430   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
 431
 432   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)GetGraph(4), *g09 = (TGraphAsymmErrors*)GetGraph(5);
 433   TAxis *ax=h2->GetXaxis();
 434   const Int_t nq(4);
 435   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
 436   Double_t yq[nq];
 437   for(Int_t ily=6; ily--;){
 438     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
 439     h1[0]->GetQuantiles(nq,yq,xq);
 440     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
 441     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
 442     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
 443     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
 444
 445     //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]);
 446     delete h1[0];
 447   }
 448 }
 449
 450 //____________________________________________________________________
 451 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
 452 {
 453 // Generic function to process one reference plot
 454
 455   Int_t n1 = 0, n2 = 0, ip=0;
 456   Double_t eff = 0.;
 457
 458   TAxis *ax = h1[0]->GetXaxis();
 459   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 460     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 461     n2 = (Int_t)h1[1]->GetBinContent(ib);
 462     eff = n2/Float_t(n1);
 463
 464     ip=g->GetN();
 465     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 466     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 467   }
 468 }
 469
 470 //____________________________________________________________________
 471 void AliTRDcheckESD::PrintStatus(ULong_t status)
 472 {
 473 // Dump track status to stdout
 474
 475   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"
 476     ,Bool_t(status & AliESDtrack::kITSin)
 477     ,Bool_t(status & AliESDtrack::kITSout)
 478     ,Bool_t(status & AliESDtrack::kITSrefit)
 479     ,Bool_t(status & AliESDtrack::kTPCin)
 480     ,Bool_t(status & AliESDtrack::kTPCout)
 481     ,Bool_t(status & AliESDtrack::kTPCrefit)
 482     ,Bool_t(status & AliESDtrack::kTPCpid)
 483     ,Bool_t(status & AliESDtrack::kTRDin)
 484     ,Bool_t(status & AliESDtrack::kTRDout)
 485     ,Bool_t(status & AliESDtrack::kTRDrefit)
 486     ,Bool_t(status & AliESDtrack::kTRDpid)
 487     ,Bool_t(status & AliESDtrack::kTRDStop)
 488     ,Bool_t(status & AliESDtrack::kHMPIDout)
 489     ,Bool_t(status & AliESDtrack::kHMPIDpid)
 490   );
 491 }
 492