a few extra null pointer checks
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckESD.cxx
CommitLineData
9a281e83 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>
48797466 21#include <TGraphErrors.h>
9a281e83 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"
d12237d6 31#include "AliESDkink.h"
9a281e83 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"
9a281e83 41
42#include "AliTRDcheckESD.h"
43
44ClassImp(AliTRDcheckESD)
45
965e895b 46const Int_t AliTRDcheckESD::fgkNgraphs = 4;
47const Float_t AliTRDcheckESD::fgkxTPC = 290.;
48const Float_t AliTRDcheckESD::fgkxTOF = 365.;
9a281e83 49
50//____________________________________________________________________
51AliTRDcheckESD::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//____________________________________________________________________
67AliTRDcheckESD::~AliTRDcheckESD()
68{
69 if(fHistos){
70 //fHistos->Delete();
71 delete fHistos;
72 }
73}
74
75//____________________________________________________________________
76void 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//____________________________________________________________________
93void 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
d12237d6 112 // status bits histogram
9a281e83 113 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
d12237d6 114 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
9a281e83 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);
48797466 120
d12237d6 121 // results array
48797466 122 TObjArray *res = new TObjArray();
123 res->SetName("Results");
124 fHistos->AddAt(res, kResults);
9a281e83 125}
126
965e895b 127//____________________________________________________________________
128TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t*)
129{
130 Bool_t kBUILD = 1, // build graph if none found
131 kCLEAR = 1; // clear existing graph
132
133 const Char_t *name[] = {
134 "Geo", "Trk", "Pid", "Ref"
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 };
142 const Int_t ngr = sizeof(name)/sizeof(Char_t*);
143 if(ngr != fgkNgraphs){
144 AliWarning("No of graphs defined different from definition");
145 return 0x0;
146 }
147
148 TObjArray *res = 0x0;
149 if(!(res = (TObjArray*)fHistos->At(kResults)) ||
150 (id < 0 || id >= ngr)){
151 AliWarning("Graph missing.");
152 return 0x0;
153 }
154
155 TGraphErrors *g = 0x0;
156 if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
157 if(kCLEAR) for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
158 } else {
159 if(kBUILD){
160 g = new TGraphErrors();
161 g->SetNameTitle(name[id], title[id]);
162 res->AddAt(g, id);
163 }
164 }
165 return g;
166}
167
9a281e83 168//____________________________________________________________________
169void AliTRDcheckESD::Exec(Option_t *){
170 //
171 // Run the Analysis
172 //
173 if(!fESD){
48797466 174 AliError("ESD event missing.");
9a281e83 175 return;
176 }
177
178 // Get MC information if available
179 AliStack * fStack = 0x0;
180 if(HasMC() && !fMC){
48797466 181 AliWarning("MC event missing");
9a281e83 182 SetMC(kFALSE);
183 } else {
184 if(!(fStack = fMC->Stack())){
48797466 185 AliWarning("MC stack missing");
9a281e83 186 SetMC(kFALSE);
187 }
188 }
d12237d6 189 Bool_t TRDin(0), TRDout(0), TRDpid(0);
9a281e83 190
9a281e83 191 AliESDtrack *esdTrack = 0x0;
192 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
d12237d6 193 TRDin=0;TRDout=0;TRDpid=0;
9a281e83 194 esdTrack = fESD->GetTrack(itrk);
d12237d6 195
196// if(esdTrack->GetNcls(1)) nTPC++;
197// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 198
199 // track status
d12237d6 200 ULong_t status = esdTrack->GetStatus();
201
202 // define TPC out tracks
203 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
204 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 205
206 // TRD PID
207 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
208 // pid quality
209 esdTrack->GetTRDpidQuality();
9a281e83 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 = fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 229 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 230 continue;
231 }
232
233 AliTrackReference *ref = 0x0;
234 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 235 if(!nRefs){
236 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
237 continue;
238 }
9a281e83 239 Int_t iref = 0;
240 while(iref<nRefs){
241 ref = mcParticle->GetTrackReference(iref);
965e895b 242 if(ref->LocalX() > fgkxTPC) break;
9a281e83 243 ref=0x0; iref++;
244 }
245
246 // read TParticle
965e895b 247 //TParticle *tParticle = mcParticle->Particle();
48797466 248 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 249 // reject secondaries
250 //if(!tParticle->IsPrimary()) continue;
9a281e83 251
d12237d6 252 if(ref){
965e895b 253 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466 254 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 255 } else {
d12237d6 256 TRDin=1;
257 if(esdTrack->GetNcls(2)) TRDout=1;
965e895b 258 if(esdTrack->GetTRDpidQuality()>=4) TRDpid=1;
9a281e83 259 }
d12237d6 260 } else { // track stopped in TPC
48797466 261 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
262 }
d12237d6 263 // get the MC pt !!
9a281e83 264 Float_t pt = ref->Pt();
d12237d6 265
9a281e83 266 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 267 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
268 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
9a281e83 269 if(/*status & AliESDtrack::k*/TRDout){
270 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 271 h->Fill(pt, kTRDout);
9a281e83 272 }
965e895b 273 if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
274 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83 275 }
9a281e83 276 PostData(0, fHistos);
277}
278
965e895b 279//____________________________________________________________________
280Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
281{
282 if(!TFile::Open(filename)){
283 AliWarning(Form("Couldn't open file %s.", filename));
284 return kFALSE;
285 }
286 TObjArray *o = 0x0;
287 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
288 AliWarning("Missing histogram container.");
289 return kFALSE;
290 }
291 fHistos = (TObjArray*)o->Clone(GetName());
292 gFile->Close();
293 return kTRUE;
294}
9a281e83 295
296//____________________________________________________________________
297void AliTRDcheckESD::Terminate(Option_t *)
298{
965e895b 299 TObjArray *res = 0x0;
300 if(!(res = (TObjArray*)fHistos->At(kResults))){
301 AliWarning("Graph container missing.");
302 return;
303 }
304 if(!res->GetEntriesFast()) res->Expand(fgkNgraphs);
48797466 305
306 // geometrical efficiency
965e895b 307 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
308 TH1 *h1[2] = {0x0, 0x0};
d12237d6 309 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
310 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
965e895b 311 Process(h1, GetGraph(0));
48797466 312
313 // tracking efficiency
d12237d6 314 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
315 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
965e895b 316 Process(h1, GetGraph(1));
48797466 317
318 // PID efficiency
965e895b 319 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6 320 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
965e895b 321 Process(h1, GetGraph(2));
322
323 // Refit efficiency
324 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
325 h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
326 Process(h1, GetGraph(3));
327}
328
329//____________________________________________________________________
330void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
331{
332 Int_t n1 = 0, n2 = 0, ip=0;
333 Double_t eff = 0.;
334
335 TAxis *ax = h1[0]->GetXaxis();
48797466 336 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
337 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
338 n2 = (Int_t)h1[1]->GetBinContent(ib);
339 eff = n2/Float_t(n1);
340
341 ip=g->GetN();
342 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
343 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
344 }
9a281e83 345}