Adding scripts for filtering
[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():
6da3eee3 52 AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83 53 ,fStatus(0)
54 ,fESD(0x0)
55 ,fMC(0x0)
56 ,fHistos(0x0)
57{
58 //
59 // Default constructor
60 //
49c7f0a5 61 SetMC(kTRUE);
9a281e83 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");
c2fd5953 99 Histos();
9a281e83 100}
101
965e895b 102//____________________________________________________________________
ea4502f6 103TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b 104{
ea4502f6 105 Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
106 kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b 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)){
ea4502f6 126 AliWarning("Graph array missing.");
965e895b 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
9a281e83 143//____________________________________________________________________
144void AliTRDcheckESD::Exec(Option_t *){
145 //
146 // Run the Analysis
147 //
148 if(!fESD){
48797466 149 AliError("ESD event missing.");
9a281e83 150 return;
151 }
152
153 // Get MC information if available
154 AliStack * fStack = 0x0;
2b104c26 155 if(HasMC()){
156 if(!fMC){
157 AliWarning("MC event missing");
9a281e83 158 SetMC(kFALSE);
2b104c26 159 } else {
160 if(!(fStack = fMC->Stack())){
161 AliWarning("MC stack missing");
162 SetMC(kFALSE);
163 }
9a281e83 164 }
165 }
d12237d6 166 Bool_t TRDin(0), TRDout(0), TRDpid(0);
9a281e83 167
9a281e83 168 AliESDtrack *esdTrack = 0x0;
169 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
d12237d6 170 TRDin=0;TRDout=0;TRDpid=0;
9a281e83 171 esdTrack = fESD->GetTrack(itrk);
d12237d6 172
173// if(esdTrack->GetNcls(1)) nTPC++;
174// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 175
176 // track status
d12237d6 177 ULong_t status = esdTrack->GetStatus();
2b104c26 178 //PrintStatus(status);
d12237d6 179
180 // define TPC out tracks
181 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
182 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 183
184 // TRD PID
185 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
186 // pid quality
ed15ef4f 187 //esdTrack->GetTRDntrackletsPID();
9a281e83 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;
7aad0c47 206 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 207 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 208 continue;
209 }
210
211 AliTrackReference *ref = 0x0;
212 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 213 if(!nRefs){
214 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
215 continue;
216 }
9a281e83 217 Int_t iref = 0;
218 while(iref<nRefs){
219 ref = mcParticle->GetTrackReference(iref);
965e895b 220 if(ref->LocalX() > fgkxTPC) break;
9a281e83 221 ref=0x0; iref++;
222 }
223
224 // read TParticle
965e895b 225 //TParticle *tParticle = mcParticle->Particle();
48797466 226 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 227 // reject secondaries
228 //if(!tParticle->IsPrimary()) continue;
9a281e83 229
d12237d6 230 if(ref){
965e895b 231 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466 232 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 233 } else {
d12237d6 234 TRDin=1;
235 if(esdTrack->GetNcls(2)) TRDout=1;
ed15ef4f 236 if(esdTrack->GetTRDntrackletsPID()>=4) TRDpid=1;
9a281e83 237 }
d12237d6 238 } else { // track stopped in TPC
48797466 239 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
240 }
d12237d6 241 // get the MC pt !!
9a281e83 242 Float_t pt = ref->Pt();
d12237d6 243
9a281e83 244 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 245 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
246 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
9a281e83 247 if(/*status & AliESDtrack::k*/TRDout){
248 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 249 h->Fill(pt, kTRDout);
9a281e83 250 }
965e895b 251 if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
252 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83 253 }
9a281e83 254 PostData(0, fHistos);
255}
256
c2fd5953 257//____________________________________________________________________
258TObjArray* 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
965e895b 291//____________________________________________________________________
292Bool_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}
9a281e83 307
308//____________________________________________________________________
309void AliTRDcheckESD::Terminate(Option_t *)
310{
965e895b 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);
48797466 317
318 // geometrical efficiency
965e895b 319 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
320 TH1 *h1[2] = {0x0, 0x0};
d12237d6 321 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
322 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
965e895b 323 Process(h1, GetGraph(0));
ea4502f6 324 delete h1[0];delete h1[1];
48797466 325
326 // tracking efficiency
d12237d6 327 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
328 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
965e895b 329 Process(h1, GetGraph(1));
ea4502f6 330 delete h1[1];
48797466 331
332 // PID efficiency
965e895b 333 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6 334 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
965e895b 335 Process(h1, GetGraph(2));
ea4502f6 336 delete h1[1];
965e895b 337
338 // Refit efficiency
339 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
340 h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
341 Process(h1, GetGraph(3));
ea4502f6 342 delete h1[1];
965e895b 343}
344
345//____________________________________________________________________
346void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
347{
348 Int_t n1 = 0, n2 = 0, ip=0;
349 Double_t eff = 0.;
350
351 TAxis *ax = h1[0]->GetXaxis();
48797466 352 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
353 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
354 n2 = (Int_t)h1[1]->GetBinContent(ib);
355 eff = n2/Float_t(n1);
356
357 ip=g->GetN();
358 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
359 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
360 }
2b104c26 361}
362
363//____________________________________________________________________
364void AliTRDcheckESD::PrintStatus(ULong_t status)
365{
366 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"
367 ,Bool_t(status & AliESDtrack::kITSin)
368 ,Bool_t(status & AliESDtrack::kITSout)
369 ,Bool_t(status & AliESDtrack::kITSrefit)
370 ,Bool_t(status & AliESDtrack::kTPCin)
371 ,Bool_t(status & AliESDtrack::kTPCout)
372 ,Bool_t(status & AliESDtrack::kTPCrefit)
373 ,Bool_t(status & AliESDtrack::kTPCpid)
374 ,Bool_t(status & AliESDtrack::kTRDin)
375 ,Bool_t(status & AliESDtrack::kTRDout)
376 ,Bool_t(status & AliESDtrack::kTRDrefit)
377 ,Bool_t(status & AliESDtrack::kTRDpid)
378 ,Bool_t(status & AliESDtrack::kTRDStop)
379 ,Bool_t(status & AliESDtrack::kHMPIDout)
380 ,Bool_t(status & AliESDtrack::kHMPIDpid)
381 );
9a281e83 382}
2b104c26 383