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 <TFile.h>
  35 #include <TTree.h>
  36 #include <TROOT.h>
  37 #include <TChain.h>
  38 #include <TParticle.h>
  39
  40 #include "AliLog.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDkink.h"
  44 #include "AliMCEvent.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliMCEventHandler.h"
  47
  48 #include "AliESDtrack.h"
  49 #include "AliMCParticle.h"
  50 #include "AliPID.h"
  51 #include "AliStack.h"
  52 #include "AliTrackReference.h"
  53
  54 #include "AliTRDcheckESD.h"
  55
  56 ClassImp(AliTRDcheckESD)
  57
  58 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
  59 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
  60 FILE* AliTRDcheckESD::fgFile = 0x0;
  61
  62 //____________________________________________________________________
  63 AliTRDcheckESD::AliTRDcheckESD():
  64   AliAnalysisTask("checkESD", "ESD checker for TRD info")
  65   ,fStatus(0)
  66   ,fESD(0x0)
  67   ,fMC(0x0)
  68   ,fHistos(0x0)
  69 {
  70   //
  71   // Default constructor
  72   //
  73   SetMC(kTRUE);
  74   DefineInput(0, TChain::Class());
  75   DefineOutput(0, TObjArray::Class());
  76 }
  77
  78 //____________________________________________________________________
  79 AliTRDcheckESD::~AliTRDcheckESD()
  80 {
  81 // Destructor
  82   if(fHistos){
  83     //fHistos->Delete();
  84     delete fHistos;
  85   }
  86 }
  87
  88 //____________________________________________________________________
  89 void AliTRDcheckESD::ConnectInputData(Option_t *)
  90 {
  91   //
  92   // Link the Input Data
  93   //
  94   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
  95   if(tree) tree->SetBranchStatus("Tracks", 1);
  96
  97   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  98   fESD = esdH ? esdH->GetEvent() : 0x0;
  99
 100   if(!HasMC()) return;
 101   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 102   fMC = mcH ? mcH->MCEvent() : 0x0;
 103 }
 104
 105 //____________________________________________________________________
 106 void AliTRDcheckESD::CreateOutputObjects()
 107 {
 108   //
 109   // Create Output Containers (TObjectArray containing 1D histograms)
 110   //
 111   OpenFile(0, "RECREATE");
 112   Histos();
 113 }
 114
 115 //____________________________________________________________________
 116 TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
 117 {
 118 // Retrieve graph with "id"
 119 // Possible options are :
 120 //   "b" - build graph if none found
 121 //   "c" - clear existing graph
 122
 123   Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
 124          kCLEAR = strstr(opt, "c"); // clear existing graph
 125
 126   const Char_t *name[] = {
 127     "Geo", "Trk", "Pid", "Ref"
 128   };
 129   const Char_t *title[] = {
 130     "TRD geometrical efficiency (TRDin/TPCout)"
 131     ,"TRD tracking efficiency (TRDout/TRDin)"
 132     ,"TRD PID efficiency (TRDpid/TRDin)"
 133     ,"TRD refit efficiency (TRDrefit/TRDin)"
 134   };
 135   const Int_t ngr = sizeof(name)/sizeof(Char_t*);
 136   if(ngr != kNgraphs){
 137     AliWarning("No of graphs defined different from definition");
 138     return 0x0;
 139   }
 140
 141   TObjArray *res = 0x0;
 142   if(!(res = (TObjArray*)fHistos->At(kResults)) ||
 143       (id < 0 || id >= ngr)){
 144     AliWarning("Graph array missing.");
 145     return 0x0;
 146   }
 147
 148   TGraphErrors *g = 0x0;
 149   if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
 150     if(kCLEAR){
 151       for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
 152     } else {
 153       PutTrendValue(name[id], g->GetMean(2), g->GetRMS(2));
 154     }
 155   } else {
 156     if(kBUILD){
 157       g = new TGraphErrors();
 158       g->SetNameTitle(name[id], title[id]);
 159       res->AddAt(g, id);
 160     }
 161   }
 162   return g;
 163 }
 164
 165 //____________________________________________________________________
 166 void AliTRDcheckESD::Exec(Option_t *){
 167   //
 168   // Run the Analysis
 169   //
 170   if(!fESD){
 171     AliError("ESD event missing.");
 172     return;
 173   }
 174
 175   // Get MC information if available
 176   AliStack * fStack = 0x0;
 177   if(HasMC()){
 178     if(!fMC){
 179       AliWarning("MC event missing");
 180       SetMC(kFALSE);
 181     } else {
 182       if(!(fStack = fMC->Stack())){
 183         AliWarning("MC stack missing");
 184         SetMC(kFALSE);
 185       }
 186     }
 187   }
 188   Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
 189
 190   AliESDtrack *esdTrack = 0x0;
 191   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 192     bTRDin=0;bTRDout=0;bTRDpid=0;
 193     esdTrack = fESD->GetTrack(itrk);
 194
 195 //     if(esdTrack->GetNcls(1)) nTPC++;
 196 //     if(esdTrack->GetNcls(2)) nTRD++;
 197
 198     // track status
 199     ULong_t status = esdTrack->GetStatus();
 200     //PrintStatus(status);
 201
 202     // define TPC out tracks
 203     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
 204     if(esdTrack->GetKinkIndex(0) > 0) continue;
 205
 206     // TRD PID
 207     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 208     // pid quality
 209     //esdTrack->GetTRDntrackletsPID();
 210
 211     // look at external track param
 212     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 213     Double_t xyz[3];
 214     if(op){
 215       op->GetXYZ(xyz);
 216       op->Global2LocalPosition(xyz, op->GetAlpha());
 217       //printf("op @ X[%7.3f]\n", xyz[0]);
 218     }
 219
 220     // read MC info
 221     if(!HasMC()) continue;
 222
 223     Int_t fLabel = esdTrack->GetLabel();
 224     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
 225
 226     // read MC particle
 227     AliMCParticle *mcParticle = 0x0;
 228     if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
 229       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 230       continue;
 231     }
 232
 233     AliTrackReference *ref = 0x0;
 234     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 235     if(!nRefs){
 236       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
 237       continue;
 238     }
 239     Int_t iref = 0;
 240     while(iref<nRefs){
 241       ref = mcParticle->GetTrackReference(iref);
 242       if(ref->LocalX() > fgkxTPC) break;
 243       ref=0x0; iref++;
 244     }
 245
 246     // read TParticle
 247     //TParticle *tParticle = mcParticle->Particle();
 248     //Int_t fPdg = tParticle->GetPdgCode();
 249     // reject secondaries
 250     //if(!tParticle->IsPrimary()) continue;
 251
 252     if(ref){
 253       if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
 254         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 255       } else {
 256         bTRDin=1;
 257         if(esdTrack->GetNcls(2)) bTRDout=1;
 258         if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
 259       }
 260     } else { // track stopped in TPC
 261       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 262     }
 263     // get the MC pt !!
 264     Float_t pt = ref->Pt();
 265
 266     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
 267     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
 268     if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
 269     if(/*status & AliESDtrack::k*/bTRDout){
 270       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
 271       h->Fill(pt, kTRDout);
 272     }
 273     if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
 274     if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
 275   }
 276   PostData(0, fHistos);
 277 }
 278
 279 //____________________________________________________________________
 280 TObjArray* AliTRDcheckESD::Histos()
 281 {
 282 // Retrieve histograms array if already build or build it
 283
 284   if(fHistos) return fHistos;
 285
 286   fHistos = new TObjArray(kNhistos);
 287   //fHistos->SetOwner(kTRUE);
 288
 289   TH1 *h = 0x0;
 290
 291   // clusters per tracklet
 292   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
 293     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
 294     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
 295     h->GetYaxis()->SetTitle("entries");
 296   } else h->Reset();
 297   fHistos->AddAt(h, kNCl);
 298
 299   // status bits histogram
 300   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 301     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
 302     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 303     h->GetYaxis()->SetTitle("status bits");
 304     h->GetZaxis()->SetTitle("entries");
 305   } else h->Reset();
 306   fHistos->AddAt(h, kTRDstat);
 307
 308   // results array
 309   TObjArray *res = new TObjArray();
 310   res->SetName("Results");
 311   fHistos->AddAt(res, kResults);
 312   return fHistos;
 313 }
 314
 315 //____________________________________________________________________
 316 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
 317 {
 318 // Load data from performance file
 319
 320   if(!TFile::Open(filename)){
 321     AliWarning(Form("Couldn't open file %s.", filename));
 322     return kFALSE;
 323   }
 324   TObjArray *o = 0x0;
 325   if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
 326     AliWarning("Missing histogram container.");
 327     return kFALSE;
 328   }
 329   fHistos = (TObjArray*)o->Clone(GetName());
 330   gFile->Close();
 331   return kTRUE;
 332 }
 333
 334 //_______________________________________________________
 335 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val, Double_t err)
 336 {
 337 // Dump trending value to default file
 338
 339   if(!fgFile){
 340     fgFile = fopen("TRD.Performance.txt", "at");
 341   }
 342   fprintf(fgFile, "%s_%s %f %f\n", GetName(), name, val, err);
 343   return kTRUE;
 344 }
 345
 346 //____________________________________________________________________
 347 void AliTRDcheckESD::Terminate(Option_t *)
 348 {
 349 // Steer post-processing
 350
 351   TObjArray *res = 0x0;
 352   if(!(res = (TObjArray*)fHistos->At(kResults))){
 353     AliWarning("Graph container missing.");
 354     return;
 355   }
 356   if(!res->GetEntriesFast()) res->Expand(kNgraphs);
 357
 358   // geometrical efficiency
 359   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
 360   TH1 *h1[2] = {0x0, 0x0};
 361   h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
 362   h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
 363   Process(h1, GetGraph(0));
 364   delete h1[0];delete h1[1];
 365
 366   // tracking efficiency
 367   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
 368   h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
 369   Process(h1, GetGraph(1));
 370   delete h1[1];
 371
 372   // PID efficiency
 373   //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
 374   h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
 375   Process(h1, GetGraph(2));
 376   delete h1[1];
 377
 378   // Refit efficiency
 379   //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
 380   h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
 381   Process(h1, GetGraph(3));
 382   delete h1[1];
 383 }
 384
 385 //____________________________________________________________________
 386 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
 387 {
 388 // Generic function to process one reference plot
 389
 390   Int_t n1 = 0, n2 = 0, ip=0;
 391   Double_t eff = 0.;
 392
 393   TAxis *ax = h1[0]->GetXaxis();
 394   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 395     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 396     n2 = (Int_t)h1[1]->GetBinContent(ib);
 397     eff = n2/Float_t(n1);
 398
 399     ip=g->GetN();
 400     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 401     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 402   }
 403 }
 404
 405 //____________________________________________________________________
 406 void AliTRDcheckESD::PrintStatus(ULong_t status)
 407 {
 408 // Dump track status to stdout
 409
 410   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"
 411     ,Bool_t(status & AliESDtrack::kITSin)
 412     ,Bool_t(status & AliESDtrack::kITSout)
 413     ,Bool_t(status & AliESDtrack::kITSrefit)
 414     ,Bool_t(status & AliESDtrack::kTPCin)
 415     ,Bool_t(status & AliESDtrack::kTPCout)
 416     ,Bool_t(status & AliESDtrack::kTPCrefit)
 417     ,Bool_t(status & AliESDtrack::kTPCpid)
 418     ,Bool_t(status & AliESDtrack::kTRDin)
 419     ,Bool_t(status & AliESDtrack::kTRDout)
 420     ,Bool_t(status & AliESDtrack::kTRDrefit)
 421     ,Bool_t(status & AliESDtrack::kTRDpid)
 422     ,Bool_t(status & AliESDtrack::kTRDStop)
 423     ,Bool_t(status & AliESDtrack::kHMPIDout)
 424     ,Bool_t(status & AliESDtrack::kHMPIDpid)
 425   );
 426 }
 427