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