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 #include "AliTRDgeometry.h"
  42
  43 #include "AliTRDcheckESD.h"
  44
  45 ClassImp(AliTRDcheckESD)
  46
  47 const Float_t AliTRDcheckESD::xTPC = 290.;
  48 const Float_t AliTRDcheckESD::xTOF = 365.;
  49
  50 //____________________________________________________________________
  51 AliTRDcheckESD::AliTRDcheckESD():
  52   AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
  53   ,fStatus(0)
  54   ,fESD(0x0)
  55   ,fMC(0x0)
  56   ,fHistos(0x0)
  57 {
  58   //
  59   // Default constructor
  60   //
  61
  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   fHistos = new TObjArray(5);
 100   //fHistos->SetOwner(kTRUE);
 101
 102   TH1 *h = 0x0;
 103
 104   // clusters per tracklet
 105   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
 106     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
 107     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
 108     h->GetYaxis()->SetTitle("entries");
 109   } else h->Reset();
 110   fHistos->AddAt(h, kNCl);
 111
 112   // status bits histogram
 113   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
 114     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
 115     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 116     h->GetYaxis()->SetTitle("status bits");
 117     h->GetZaxis()->SetTitle("entries");
 118   } else h->Reset();
 119   fHistos->AddAt(h, kTRDstat);
 120
 121   // results array
 122   TObjArray *res = new TObjArray();
 123   res->SetName("Results");
 124   fHistos->AddAt(res, kResults);
 125 }
 126
 127 //____________________________________________________________________
 128 void AliTRDcheckESD::Exec(Option_t *){
 129   //
 130   // Run the Analysis
 131   //
 132   if(!fESD){
 133     AliError("ESD event missing.");
 134     return;
 135   }
 136
 137   // Get MC information if available
 138   AliStack * fStack = 0x0;
 139   if(HasMC() && !fMC){
 140     AliWarning("MC event missing");
 141     SetMC(kFALSE);
 142   } else {
 143     if(!(fStack = fMC->Stack())){
 144       AliWarning("MC stack missing");
 145       SetMC(kFALSE);
 146     }
 147   }
 148   Bool_t TRDin(0), TRDout(0), TRDpid(0);
 149
 150   AliESDtrack *esdTrack = 0x0;
 151   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
 152     TRDin=0;TRDout=0;TRDpid=0;
 153     esdTrack = fESD->GetTrack(itrk);
 154
 155 //     if(esdTrack->GetNcls(1)) nTPC++;
 156 //     if(esdTrack->GetNcls(2)) nTRD++;
 157
 158     // track status
 159     ULong_t status = esdTrack->GetStatus();
 160
 161     // define TPC out tracks
 162     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
 163     if(esdTrack->GetKinkIndex(0) > 0) continue;
 164
 165     // TRD PID
 166     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 167     // pid quality
 168     esdTrack->GetTRDpidQuality();
 169
 170     // look at external track param
 171     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 172     Double_t xyz[3];
 173     if(op){
 174       op->GetXYZ(xyz);
 175       op->Global2LocalPosition(xyz, op->GetAlpha());
 176       //printf("op @ X[%7.3f]\n", xyz[0]);
 177     }
 178
 179     // read MC info
 180     if(!HasMC()) continue;
 181
 182     Int_t fLabel = esdTrack->GetLabel();
 183     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
 184
 185     // read MC particle
 186     AliMCParticle *mcParticle = 0x0;
 187     if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
 188       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
 189       continue;
 190     }
 191
 192     AliTrackReference *ref = 0x0;
 193     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 194     if(!nRefs){
 195       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
 196       continue;
 197     }
 198     Int_t iref = 0;
 199     while(iref<nRefs){
 200       ref = mcParticle->GetTrackReference(iref);
 201       if(ref->LocalX() > xTPC) break;
 202       ref=0x0; iref++;
 203     }
 204
 205     // read TParticle
 206     TParticle *tParticle = mcParticle->Particle();
 207     //Int_t fPdg = tParticle->GetPdgCode();
 208     // reject secondaries
 209     //if(!tParticle->IsPrimary()) continue;
 210
 211     if(ref){
 212       if(ref->LocalX() > xTOF){ // track skipping TRD fiducial volume
 213         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 214       } else {
 215         TRDin=1;
 216         if(esdTrack->GetNcls(2)) TRDout=1;
 217         if(esdTrack->GetTRDpidQuality()) TRDpid=1;
 218       }
 219     } else { // track stopped in TPC
 220       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
 221     }
 222     // get the MC pt !!
 223     Float_t pt = ref->Pt();
 224
 225     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
 226     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
 227     if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
 228     if(/*status & AliESDtrack::k*/TRDout){
 229       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
 230       h->Fill(pt, kTRDout);
 231     }
 232     if(status & AliESDtrack::kTRDpid) h->Fill(pt, kTRDpid);
 233   }
 234   PostData(0, fHistos);
 235 }
 236
 237
 238 //____________________________________________________________________
 239 void AliTRDcheckESD::Terminate(Option_t *)
 240 {
 241   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
 242   TAxis *ax = h2->GetXaxis();
 243   TObjArray *res = (TObjArray*)fHistos->At(kResults);
 244
 245   TH1 *h1[2] = {0x0, 0x0};
 246   TGraphErrors *g = 0x0;
 247   res->Expand(3);
 248   Int_t n1 = 0, n2 = 0, ip=0;
 249   Double_t eff = 0.;
 250
 251   // geometrical efficiency
 252   h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
 253   h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
 254   res->Add(g = new TGraphErrors());
 255   g->SetNameTitle("geom", "TRD geometrical efficiency (TRDin/TPCout)");
 256   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 257     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 258     n2 = (Int_t)h1[1]->GetBinContent(ib);
 259     eff = n2/Float_t(n1);
 260
 261     ip=g->GetN();
 262     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 263     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 264   }
 265
 266   // tracking efficiency
 267   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
 268   h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
 269   res->Add(g = new TGraphErrors());
 270   g->SetNameTitle("tracking", "TRD tracking efficiency (TRDout/TRDin)");
 271   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 272     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 273     n2 = (Int_t)h1[1]->GetBinContent(ib);
 274     eff = n2/Float_t(n1);
 275
 276     ip=g->GetN();
 277     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 278     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 279   }
 280
 281   // PID efficiency
 282   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
 283   h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
 284   res->Add(g = new TGraphErrors());
 285   g->SetNameTitle("PID", "TRD PID efficiency (TRDpid/TRDin)");
 286   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
 287     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
 288     n2 = (Int_t)h1[1]->GetBinContent(ib);
 289     eff = n2/Float_t(n1);
 290
 291     ip=g->GetN();
 292     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
 293     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
 294   }
 295 }
 296
 297 // if(TRDin && !TRDout){
 298 //   printf("[%c] x[%2d]=%7.2f pt[%7.3f] id[%d] kink[%d] label[%d] in[%d] out[%d] pdg[%d]\n", tParticle->IsPrimary() ? 'P' : 'S', iref, ref->LocalX(), pt, esdTrack->GetID(), esdTrack->GetKinkIndex(0), fLabel, TRDin, TRDout, tParticle->GetPdgCode());
 299 //   printf("\tstatus ITS[%d %d] TPC[%d %d %d] \n"
 300 //     ,Bool_t(status & AliESDtrack::kITSin)
 301 //     ,Bool_t(status & AliESDtrack::kITSout)
 302 //     ,Bool_t(status & AliESDtrack::kTPCin)
 303 //     ,Bool_t(status & AliESDtrack::kTPCout)
 304 //     ,Bool_t(status & AliESDtrack::kTPCrefit)
 305 //   );
 306 //   printf("clusters ITS[%d] TPC[%d] TRD[%d]\n", esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2));
 307 //   if(esdTrack->GetKinkIndex(0) != 0){
 308 //     AliESDkink *kink = fESD->GetKink(TMath::Abs(esdTrack->GetKinkIndex(0)));
 309 //     if(kink) printf("index[%d %d] label[%d %d]\n", kink->GetIndex(0), kink->GetIndex(1), kink->GetLabel(0), kink->GetLabel(1));
 310 //   }
 311 // }