fix compilation mismatches in the new framework
[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>
eede30c8 34#include <TGraphAsymmErrors.h>
9a281e83 35#include <TFile.h>
36#include <TTree.h>
37#include <TROOT.h>
38#include <TChain.h>
39#include <TParticle.h>
40
41#include "AliLog.h"
42#include "AliAnalysisManager.h"
43#include "AliESDEvent.h"
d12237d6 44#include "AliESDkink.h"
9a281e83 45#include "AliMCEvent.h"
46#include "AliESDInputHandler.h"
47#include "AliMCEventHandler.h"
48
49#include "AliESDtrack.h"
50#include "AliMCParticle.h"
51#include "AliPID.h"
52#include "AliStack.h"
53#include "AliTrackReference.h"
9a281e83 54
55#include "AliTRDcheckESD.h"
56
57ClassImp(AliTRDcheckESD)
58
965e895b 59const Float_t AliTRDcheckESD::fgkxTPC = 290.;
60const Float_t AliTRDcheckESD::fgkxTOF = 365.;
59589fc9 61FILE* AliTRDcheckESD::fgFile = 0x0;
9a281e83 62
63//____________________________________________________________________
64AliTRDcheckESD::AliTRDcheckESD():
6da3eee3 65 AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83 66 ,fStatus(0)
67 ,fESD(0x0)
68 ,fMC(0x0)
69 ,fHistos(0x0)
eede30c8 70 ,fResults(0x0)
9a281e83 71{
72 //
73 // Default constructor
74 //
49c7f0a5 75 SetMC(kTRUE);
9a281e83 76 DefineInput(0, TChain::Class());
77 DefineOutput(0, TObjArray::Class());
78}
79
80//____________________________________________________________________
81AliTRDcheckESD::~AliTRDcheckESD()
82{
ddd1a033 83// Destructor
9a281e83 84 if(fHistos){
85 //fHistos->Delete();
86 delete fHistos;
87 }
eede30c8 88 if(fResults){
89 fResults->Delete();
90 delete fResults;
91 }
9a281e83 92}
93
94//____________________________________________________________________
95void AliTRDcheckESD::ConnectInputData(Option_t *)
96{
97 //
98 // Link the Input Data
99 //
100 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
101 if(tree) tree->SetBranchStatus("Tracks", 1);
102
103 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104 fESD = esdH ? esdH->GetEvent() : 0x0;
105
106 if(!HasMC()) return;
107 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
108 fMC = mcH ? mcH->MCEvent() : 0x0;
109}
110
111//____________________________________________________________________
112void AliTRDcheckESD::CreateOutputObjects()
113{
114 //
115 // Create Output Containers (TObjectArray containing 1D histograms)
116 //
117 OpenFile(0, "RECREATE");
c2fd5953 118 Histos();
9a281e83 119}
120
965e895b 121//____________________________________________________________________
eede30c8 122TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b 123{
ddd1a033 124// Retrieve graph with "id"
125// Possible options are :
126// "b" - build graph if none found
127// "c" - clear existing graph
128
ea4502f6 129 Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
130 kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b 131
132 const Char_t *name[] = {
eede30c8 133 "Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
965e895b 134 };
135 const Char_t *title[] = {
136 "TRD geometrical efficiency (TRDin/TPCout)"
137 ,"TRD tracking efficiency (TRDout/TRDin)"
138 ,"TRD PID efficiency (TRDpid/TRDin)"
139 ,"TRD refit efficiency (TRDrefit/TRDin)"
eede30c8 140 ,"TRD Eloss (Max/90% quantile)"
141 ,"TRD Eloss (Mean/60% quantile)"
965e895b 142 };
143 const Int_t ngr = sizeof(name)/sizeof(Char_t*);
ddd1a033 144 if(ngr != kNgraphs){
965e895b 145 AliWarning("No of graphs defined different from definition");
146 return 0x0;
147 }
148
eede30c8 149 if(!fResults){
150 fResults = new TObjArray(kNgraphs);
151 fResults->SetOwner();
152 fResults->SetName("results");
965e895b 153 }
154
eede30c8 155 TGraph *g = 0x0;
156 if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
59589fc9 157 if(kCLEAR){
158 for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
159 } else {
5e17c62a 160 PutTrendValue(name[id], g->GetMean(2));
161 PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
59589fc9 162 }
965e895b 163 } else {
164 if(kBUILD){
eede30c8 165 switch(id){
166 case 0:
167 g = new TGraphErrors();
168 break;
169 case 1:
170 g = new TGraphErrors();
171 break;
172 case 2:
173 g = new TGraphErrors();
174 break;
175 case 3:
176 g = new TGraphErrors();
177 break;
178 case 4:
179 g = new TGraphAsymmErrors(6);
180 g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
181 g->SetLineColor(kBlack);g->SetLineWidth(2);
182 break;
183 case 5:
184 g = new TGraphAsymmErrors(6);
185 g->SetMarkerStyle(21);
186 g->SetLineColor(kRed);g->SetLineWidth(2);
187 break;
188 default:
189 AliWarning(Form("Graph index[%d] missing/not defined.", id));
190 return 0x0;
191 }
965e895b 192 g->SetNameTitle(name[id], title[id]);
eede30c8 193 fResults->AddAt(g, id);
965e895b 194 }
195 }
196 return g;
197}
198
9a281e83 199//____________________________________________________________________
200void AliTRDcheckESD::Exec(Option_t *){
201 //
202 // Run the Analysis
203 //
204 if(!fESD){
48797466 205 AliError("ESD event missing.");
9a281e83 206 return;
207 }
208
209 // Get MC information if available
210 AliStack * fStack = 0x0;
2b104c26 211 if(HasMC()){
212 if(!fMC){
213 AliWarning("MC event missing");
9a281e83 214 SetMC(kFALSE);
2b104c26 215 } else {
216 if(!(fStack = fMC->Stack())){
217 AliWarning("MC stack missing");
218 SetMC(kFALSE);
219 }
9a281e83 220 }
221 }
ddd1a033 222 Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
9a281e83 223
9a281e83 224 AliESDtrack *esdTrack = 0x0;
225 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
ddd1a033 226 bTRDin=0;bTRDout=0;bTRDpid=0;
9a281e83 227 esdTrack = fESD->GetTrack(itrk);
d12237d6 228
229// if(esdTrack->GetNcls(1)) nTPC++;
230// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 231
232 // track status
d12237d6 233 ULong_t status = esdTrack->GetStatus();
2b104c26 234 //PrintStatus(status);
d12237d6 235
236 // define TPC out tracks
237 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
238 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 239
240 // TRD PID
241 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
242 // pid quality
ed15ef4f 243 //esdTrack->GetTRDntrackletsPID();
9a281e83 244
245 // look at external track param
246 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
eede30c8 247 const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
9a281e83 248 Double_t xyz[3];
249 if(op){
250 op->GetXYZ(xyz);
251 op->Global2LocalPosition(xyz, op->GetAlpha());
252 //printf("op @ X[%7.3f]\n", xyz[0]);
253 }
254
255 // read MC info
256 if(!HasMC()) continue;
257
258 Int_t fLabel = esdTrack->GetLabel();
259 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
260
261 // read MC particle
262 AliMCParticle *mcParticle = 0x0;
7aad0c47 263 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 264 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 265 continue;
266 }
267
268 AliTrackReference *ref = 0x0;
269 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 270 if(!nRefs){
271 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
272 continue;
273 }
9a281e83 274 Int_t iref = 0;
275 while(iref<nRefs){
276 ref = mcParticle->GetTrackReference(iref);
965e895b 277 if(ref->LocalX() > fgkxTPC) break;
9a281e83 278 ref=0x0; iref++;
279 }
280
281 // read TParticle
965e895b 282 //TParticle *tParticle = mcParticle->Particle();
48797466 283 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 284 // reject secondaries
285 //if(!tParticle->IsPrimary()) continue;
9a281e83 286
d12237d6 287 if(ref){
965e895b 288 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466 289 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 290 } else {
ddd1a033 291 bTRDin=1;
292 if(esdTrack->GetNcls(2)) bTRDout=1;
293 if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
9a281e83 294 }
d12237d6 295 } else { // track stopped in TPC
48797466 296 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
297 }
d12237d6 298 // get the MC pt !!
9a281e83 299 Float_t pt = ref->Pt();
d12237d6 300
9a281e83 301 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 302 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
ddd1a033 303 if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
304 if(/*status & AliESDtrack::k*/bTRDout){
9a281e83 305 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 306 h->Fill(pt, kTRDout);
9a281e83 307 }
ddd1a033 308 if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
965e895b 309 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
eede30c8 310
311 if(ip){
312 h = (TH2I*)fHistos->At(kTRDmom);
33711b35 313 Float_t pp(0.);
eede30c8 314 for(Int_t ily=6; ily--;){
33711b35 315 if((pp=esdTrack->GetTRDmomentum(ily))<0.) continue;
316 h->Fill(ip->GetP()-pp, ily);
eede30c8 317 }
318 }
9a281e83 319 }
9a281e83 320 PostData(0, fHistos);
321}
322
c2fd5953 323//____________________________________________________________________
324TObjArray* AliTRDcheckESD::Histos()
325{
ddd1a033 326// Retrieve histograms array if already build or build it
327
c2fd5953 328 if(fHistos) return fHistos;
329
330 fHistos = new TObjArray(kNhistos);
331 //fHistos->SetOwner(kTRUE);
332
333 TH1 *h = 0x0;
334
335 // clusters per tracklet
336 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
337 h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
338 h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
339 h->GetYaxis()->SetTitle("entries");
340 } else h->Reset();
341 fHistos->AddAt(h, kNCl);
342
343 // status bits histogram
344 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
345 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
346 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
347 h->GetYaxis()->SetTitle("status bits");
348 h->GetZaxis()->SetTitle("entries");
349 } else h->Reset();
350 fHistos->AddAt(h, kTRDstat);
351
eede30c8 352 // energy loss
353 if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
354 h = new TH2I("hTRDmom", "TRD energy loss", 100, -1., 2., 6, -0.5, 5.5);
355 h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
356 h->GetYaxis()->SetTitle("layer");
357 h->GetZaxis()->SetTitle("entries");
358 } else h->Reset();
359 fHistos->AddAt(h, kTRDmom);
360
c2fd5953 361 return fHistos;
362}
363
965e895b 364//____________________________________________________________________
365Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
366{
ddd1a033 367// Load data from performance file
368
965e895b 369 if(!TFile::Open(filename)){
370 AliWarning(Form("Couldn't open file %s.", filename));
371 return kFALSE;
372 }
373 TObjArray *o = 0x0;
374 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
375 AliWarning("Missing histogram container.");
376 return kFALSE;
377 }
378 fHistos = (TObjArray*)o->Clone(GetName());
379 gFile->Close();
380 return kTRUE;
381}
9a281e83 382
59589fc9 383//_______________________________________________________
5e17c62a 384Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
59589fc9 385{
ddd1a033 386// Dump trending value to default file
387
59589fc9 388 if(!fgFile){
389 fgFile = fopen("TRD.Performance.txt", "at");
390 }
5e17c62a 391 fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
59589fc9 392 return kTRUE;
393}
394
9a281e83 395//____________________________________________________________________
396void AliTRDcheckESD::Terminate(Option_t *)
397{
ddd1a033 398// Steer post-processing
399
eede30c8 400
48797466 401 // geometrical efficiency
965e895b 402 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
403 TH1 *h1[2] = {0x0, 0x0};
dad3aa6b 404 h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
405 h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
eede30c8 406 Process(h1, (TGraphErrors*)GetGraph(0));
ea4502f6 407 delete h1[0];delete h1[1];
48797466 408
409 // tracking efficiency
dad3aa6b 410 h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
411 h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
eede30c8 412 Process(h1, (TGraphErrors*)GetGraph(1));
ea4502f6 413 delete h1[1];
48797466 414
415 // PID efficiency
dad3aa6b 416 h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
eede30c8 417 Process(h1, (TGraphErrors*)GetGraph(2));
ea4502f6 418 delete h1[1];
965e895b 419
420 // Refit efficiency
dad3aa6b 421 h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
eede30c8 422 Process(h1, (TGraphErrors*)GetGraph(3));
ea4502f6 423 delete h1[1];
eede30c8 424 if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
425
426 TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)GetGraph(4), *g09 = (TGraphAsymmErrors*)GetGraph(5);
427 TAxis *ax=h2->GetXaxis();
428 const Int_t nq(4);
429 const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
430 Double_t yq[nq];
431 for(Int_t ily=6; ily--;){
432 h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
433 h1[0]->GetQuantiles(nq,yq,xq);
434 g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
435 g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
436 g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
437 g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
438
439 //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
440 delete h1[0];
441 }
965e895b 442}
443
444//____________________________________________________________________
445void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
446{
ddd1a033 447// Generic function to process one reference plot
448
965e895b 449 Int_t n1 = 0, n2 = 0, ip=0;
450 Double_t eff = 0.;
451
452 TAxis *ax = h1[0]->GetXaxis();
48797466 453 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
454 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
455 n2 = (Int_t)h1[1]->GetBinContent(ib);
456 eff = n2/Float_t(n1);
457
458 ip=g->GetN();
459 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
460 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
461 }
2b104c26 462}
463
464//____________________________________________________________________
465void AliTRDcheckESD::PrintStatus(ULong_t status)
466{
ddd1a033 467// Dump track status to stdout
468
2b104c26 469 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"
470 ,Bool_t(status & AliESDtrack::kITSin)
471 ,Bool_t(status & AliESDtrack::kITSout)
472 ,Bool_t(status & AliESDtrack::kITSrefit)
473 ,Bool_t(status & AliESDtrack::kTPCin)
474 ,Bool_t(status & AliESDtrack::kTPCout)
475 ,Bool_t(status & AliESDtrack::kTPCrefit)
476 ,Bool_t(status & AliESDtrack::kTPCpid)
477 ,Bool_t(status & AliESDtrack::kTRDin)
478 ,Bool_t(status & AliESDtrack::kTRDout)
479 ,Bool_t(status & AliESDtrack::kTRDrefit)
480 ,Bool_t(status & AliESDtrack::kTRDpid)
481 ,Bool_t(status & AliESDtrack::kTRDStop)
482 ,Bool_t(status & AliESDtrack::kHMPIDout)
483 ,Bool_t(status & AliESDtrack::kHMPIDpid)
484 );
9a281e83 485}
2b104c26 486