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