fix coding convention errors
[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
ddd1a033 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//////////////////////////////////////////////////////
9a281e83 28
29#include <TClonesArray.h>
30#include <TObjArray.h>
31#include <TObject.h>
32#include <TH2I.h>
48797466 33#include <TGraphErrors.h>
9a281e83 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"
d12237d6 43#include "AliESDkink.h"
9a281e83 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"
9a281e83 53
54#include "AliTRDcheckESD.h"
55
56ClassImp(AliTRDcheckESD)
57
965e895b 58const Float_t AliTRDcheckESD::fgkxTPC = 290.;
59const Float_t AliTRDcheckESD::fgkxTOF = 365.;
59589fc9 60FILE* AliTRDcheckESD::fgFile = 0x0;
9a281e83 61
62//____________________________________________________________________
63AliTRDcheckESD::AliTRDcheckESD():
6da3eee3 64 AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83 65 ,fStatus(0)
66 ,fESD(0x0)
67 ,fMC(0x0)
68 ,fHistos(0x0)
69{
70 //
71 // Default constructor
72 //
49c7f0a5 73 SetMC(kTRUE);
9a281e83 74 DefineInput(0, TChain::Class());
75 DefineOutput(0, TObjArray::Class());
76}
77
78//____________________________________________________________________
79AliTRDcheckESD::~AliTRDcheckESD()
80{
ddd1a033 81// Destructor
9a281e83 82 if(fHistos){
83 //fHistos->Delete();
84 delete fHistos;
85 }
86}
87
88//____________________________________________________________________
89void 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//____________________________________________________________________
106void AliTRDcheckESD::CreateOutputObjects()
107{
108 //
109 // Create Output Containers (TObjectArray containing 1D histograms)
110 //
111 OpenFile(0, "RECREATE");
c2fd5953 112 Histos();
9a281e83 113}
114
965e895b 115//____________________________________________________________________
ea4502f6 116TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b 117{
ddd1a033 118// Retrieve graph with "id"
119// Possible options are :
120// "b" - build graph if none found
121// "c" - clear existing graph
122
ea4502f6 123 Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
124 kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b 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*);
ddd1a033 136 if(ngr != kNgraphs){
965e895b 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)){
ea4502f6 144 AliWarning("Graph array missing.");
965e895b 145 return 0x0;
146 }
147
148 TGraphErrors *g = 0x0;
149 if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
59589fc9 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 }
965e895b 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
9a281e83 165//____________________________________________________________________
166void AliTRDcheckESD::Exec(Option_t *){
167 //
168 // Run the Analysis
169 //
170 if(!fESD){
48797466 171 AliError("ESD event missing.");
9a281e83 172 return;
173 }
174
175 // Get MC information if available
176 AliStack * fStack = 0x0;
2b104c26 177 if(HasMC()){
178 if(!fMC){
179 AliWarning("MC event missing");
9a281e83 180 SetMC(kFALSE);
2b104c26 181 } else {
182 if(!(fStack = fMC->Stack())){
183 AliWarning("MC stack missing");
184 SetMC(kFALSE);
185 }
9a281e83 186 }
187 }
ddd1a033 188 Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
9a281e83 189
9a281e83 190 AliESDtrack *esdTrack = 0x0;
191 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
ddd1a033 192 bTRDin=0;bTRDout=0;bTRDpid=0;
9a281e83 193 esdTrack = fESD->GetTrack(itrk);
d12237d6 194
195// if(esdTrack->GetNcls(1)) nTPC++;
196// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 197
198 // track status
d12237d6 199 ULong_t status = esdTrack->GetStatus();
2b104c26 200 //PrintStatus(status);
d12237d6 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
ed15ef4f 209 //esdTrack->GetTRDntrackletsPID();
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;
7aad0c47 228 if(!(mcParticle = (AliMCParticle*) 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 {
ddd1a033 256 bTRDin=1;
257 if(esdTrack->GetNcls(2)) bTRDout=1;
258 if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=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);
ddd1a033 268 if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
269 if(/*status & AliESDtrack::k*/bTRDout){
9a281e83 270 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 271 h->Fill(pt, kTRDout);
9a281e83 272 }
ddd1a033 273 if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
965e895b 274 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83 275 }
9a281e83 276 PostData(0, fHistos);
277}
278
c2fd5953 279//____________________________________________________________________
280TObjArray* AliTRDcheckESD::Histos()
281{
ddd1a033 282// Retrieve histograms array if already build or build it
283
c2fd5953 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
965e895b 315//____________________________________________________________________
316Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
317{
ddd1a033 318// Load data from performance file
319
965e895b 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}
9a281e83 333
59589fc9 334//_______________________________________________________
335Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val, Double_t err)
336{
ddd1a033 337// Dump trending value to default file
338
59589fc9 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
9a281e83 346//____________________________________________________________________
347void AliTRDcheckESD::Terminate(Option_t *)
348{
ddd1a033 349// Steer post-processing
350
965e895b 351 TObjArray *res = 0x0;
352 if(!(res = (TObjArray*)fHistos->At(kResults))){
353 AliWarning("Graph container missing.");
354 return;
355 }
ddd1a033 356 if(!res->GetEntriesFast()) res->Expand(kNgraphs);
48797466 357
358 // geometrical efficiency
965e895b 359 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
360 TH1 *h1[2] = {0x0, 0x0};
d12237d6 361 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
362 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
965e895b 363 Process(h1, GetGraph(0));
ea4502f6 364 delete h1[0];delete h1[1];
48797466 365
366 // tracking efficiency
d12237d6 367 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
368 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
965e895b 369 Process(h1, GetGraph(1));
ea4502f6 370 delete h1[1];
48797466 371
372 // PID efficiency
965e895b 373 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6 374 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
965e895b 375 Process(h1, GetGraph(2));
ea4502f6 376 delete h1[1];
965e895b 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));
ea4502f6 382 delete h1[1];
965e895b 383}
384
385//____________________________________________________________________
386void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
387{
ddd1a033 388// Generic function to process one reference plot
389
965e895b 390 Int_t n1 = 0, n2 = 0, ip=0;
391 Double_t eff = 0.;
392
393 TAxis *ax = h1[0]->GetXaxis();
48797466 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 }
2b104c26 403}
404
405//____________________________________________________________________
406void AliTRDcheckESD::PrintStatus(ULong_t status)
407{
ddd1a033 408// Dump track status to stdout
409
2b104c26 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 );
9a281e83 426}
2b104c26 427