- Clean-up of data members and methods in base classes and derived classes
[u/mrichter/AliRoot.git] / ESDCheck / AliTRDQATask.cxx
CommitLineData
8d14dc14 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// An analysis task to check the TRD data in simulated data
17//
18//*-- Sylwester Radomski
19//////////////////////////////////////////////////////////////////////////////
20// track type codding
21//
22// tpci = kTPCin
23// tpco = kTPCout
24// tpcz = kTPCout && !kTRDout
25// trdo = kTRDout
26// trdr = kTRDref
27// trdz = kTRDout && !kTRDref
28//
29
30#include "AliTRDQATask.h"
31
32#include "TChain.h"
33#include "TH1D.h"
34#include "TH2D.h"
35#include "TFile.h"
36#include "TStyle.h"
37#include "TGaxis.h"
38#include "TCanvas.h"
39
40#include "AliESD.h"
41#include "AliLog.h"
42
43//______________________________________________________________________________
44AliTRDQATask::AliTRDQATask(const char *name) :
45 AliAnalysisTask(name,""),
46 fChain(0),
47 fESD(0)
48{
49 // Constructor.
50 // Input slot #0 works with an Ntuple
51 DefineInput(0, TChain::Class());
52 // Output slot #0 writes into a TH1 container
53 DefineOutput(0, TObjArray::Class()) ;
54}
55
56//______________________________________________________________________________
57void AliTRDQATask::Init(const Option_t *)
58{
59 // Initialisation of branch container and histograms
60
61 AliInfo(Form("*** Initialization of %s", GetName())) ;
62
63 // Get input data
64 fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
65 if (!fChain) {
66 AliError(Form("Input 0 for %s not found\n", GetName()));
67 return ;
68 }
69
70 if (!fESD) {
71 // One should first check if the branch address was taken by some other task
72 char ** address = (char **)GetBranchAddress(0, "ESD");
73 if (address) {
74 fESD = (AliESD *)(*address);
75 AliInfo("Old ESD found.");
76 }
77 if (!fESD) {
78 fESD = new AliESD();
79 SetBranchAddress(0, "ESD", &fESD);
80 if (fESD) AliInfo("ESD connected.");
81 }
82 }
83 // The output objects will be written to
84 TDirectory * cdir = gDirectory ;
85 OpenFile(0, Form("%s.root", GetName()), "RECREATE");
86 if (cdir)
87 cdir->cd() ;
88
89 // create histograms
90
91 fNTracks = new TH1D("ntracks", ";number of all tracks", 500, -0.5, 499.5);
92 fEventSize = new TH1D("evSize", ";event size (MB)", 100, 0, 5);
93
94 fTrackStatus = new TH1D("trackStatus", ";status bit", 32, -0.5, 31.5);
95 fKinkIndex = new TH1D("kinkIndex", ";kink index", 41, -20.5, 20.5);
96
97 fParIn = new TH1D("parIn", "Inner Plane", 2, -0.5, 1.5);
98 fParOut = new TH1D("parOut", "Outer Plane", 2, -0.5, 1.5);
99
100 fXIn = new TH1D("xIn", ";X at the inner plane (cm)", 200, 50, 250);
101 fXOut = new TH1D("xOut", ";X at the outer plane (cm)", 300, 50, 400);
102
103 const int nNameAlpha = 4;
104 const char *namesAlpha[nNameAlpha] = {"alphaTPCi", "alphaTPCo", "alphaTRDo", "alphaTRDr"};
105 //TH1D *fAlpha[4];
106 for(int i=0; i<nNameAlpha; i++) {
107 fAlpha[i] = new TH1D(namesAlpha[i], "alpha", 100, -4, 4);
108 }
109 fSectorTRD = new TH1D ("sectorTRD", ";sector TRD", 20, -0.5, 19.5);
110
111
112 // track parameters
113 const int nbits = 6;
114 const char *suf[nbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
115 for(int i=0; i<nbits; i++) {
116 fPt[i] = new TH1D(Form("pt%s",suf[i]), ";p_{T} (GeV/c);entries TPC", 50, 0, 10);
117 fTheta[i] = new TH1D(Form("theta%s", suf[i]), "theta (rad)", 100, -4, 4);
118 fSigmaY[i] = new TH1D(Form("sigmaY%s",suf[i]), ";sigma Y (cm)", 200, 0, 1);
119 fChi2[i] = new TH1D(Form("Chi2%s", suf[i]), ";#chi2 / ndf", 100, 0, 10);
120 fPlaneYZ[i] = new TH2D(Form("planeYZ%s", suf[i]), Form("%sy (cm);z (cm)", suf[i]),
121 100, -60, 60, 500, -500, 500);
122 }
123
124 // efficiency
125 fEffPt[0] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[0], suf[1]));
126 fEffPt[1] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[1], suf[3]));
127 fEffPt[2] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[3], suf[4]));
128 fEffPt[3] = (TH1D*) fPt[0]->Clone(Form("eff_%s_%s", suf[1], suf[4]));
129
130 for(int i=0; i<4; i++) {
131 fEffPt[i]->Sumw2();
132 fEffPt[i]->SetMarkerStyle(20);
133 fEffPt[i]->SetMinimum(0);
134 fEffPt[i]->SetMaximum(1.1);
135 }
136
137 // track features
138 fClustersTRD[0] = new TH1D("clsTRDo", "TRDo;number of clusters", 130, -0.5, 129.5);;
139 fClustersTRD[1] = new TH1D("clsTRDr", "TRDr;number of clusters", 130, -0.5, 129.5);;
140 fClustersTRD[2] = new TH1D("clsTRDz", "TRDz;number of clusters", 130, -0.5, 129.5);;
141
142 // for good refitted tracks only
143 fTime = new TH1D("time", ";time bin", 25, -0.5, 24.5);
144 fBudget = new TH1D("budget", ";material budget", 100, 0, 100);
145 fQuality = new TH1D("quality", ";track quality", 100, 0, 1.1);
146 fSignal = new TH1D("signal", ";signal", 100, 0, 1e3);
147
148 // dEdX and PID
149 fTrdSigMom = new TH2D("trdSigMom", ";momentum (GeV/c);signal", 100, 0, 3, 100, 0, 1e3);
150 fTpcSigMom = new TH2D("tpcSigMom", ";momentum (GeV/c);signal", 100, 0, 3, 100, 0, 200);
151
152 const char *pidName[6] = {"El", "Mu", "Pi", "K", "P", "Ch"};
153 for(int i=0; i<6; i++) {
154
155 // TPC
156 fTpcPID[i] = new TH1D(Form("tpcPid%s",pidName[i]), pidName[i], 100, 0, 1.5);
157 fTpcPID[i]->GetXaxis()->SetTitle("probability");
158
159 fTpcSigMomPID[i] = new TH2D(Form("tpcSigMom%s",pidName[i]), "", 100, 0, 3, 100, 0, 200);
160 fTpcSigMomPID[i]->SetTitle(Form("%s;momentum (GeV/c);signal",pidName[i]));
161
162 // TRD
163 fTrdPID[i] = new TH1D(Form("trdPid%s",pidName[i]), pidName[i], 100, 0, 1.5);
164 fTrdPID[i]->GetXaxis()->SetTitle("probability");
165
166 fTrdSigMomPID[i] = new TH2D(Form("trdSigMom%s",pidName[i]), "", 100, 0, 3, 100, 0, 1e3);
167 fTrdSigMomPID[i]->SetTitle(Form("%s;momentum (GeV/c);signal",pidName[i]));
168 }
169
170
171 // create output container
172 fOutputContainer = new TObjArray(80);
173
174 // register histograms to the container
175 TIter next(gDirectory->GetList());
176 TObject *obj;
177 int counter = 0;
178
179 while (obj = next.Next()) {
180 if (obj->InheritsFrom("TH1")) fOutputContainer->AddAt(obj, counter++);
181 }
182
183 AliInfo(Form("Number of histograms = %d", counter));
184
185 }
186
187//______________________________________________________________________________
188void AliTRDQATask::Exec(Option_t *)
189{
190 // Process one event
191
192 Long64_t entry = fChain->GetReadEntry() ;
193
194 // Processing of one event
195
196 if (!fESD) {
197 AliError("fESD is not connected to the input!") ;
198 return ;
199 }
200
201 if ( !((entry-1)%100) )
202 AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
203
204 int nTracks = fESD->GetNumberOfTracks();
205 fNTracks->Fill(nTracks);
206
207 // track loop
208 for(int i=0; i<nTracks; i++) {
209
210 AliESDtrack *track = fESD->GetTrack(i);
211 const AliExternalTrackParam *paramOut = track->GetOuterParam();
212 const AliExternalTrackParam *paramIn = track->GetInnerParam();
213
214 fParIn->Fill(!!paramIn);
215 if (!paramIn) continue;
216 fXIn->Fill(paramIn->GetX());
217
218 fParOut->Fill(!!paramOut);
219 if (!paramOut) continue;
220 fXOut->Fill(paramOut->GetX());
221
222 int sector = GetSector(paramOut->GetAlpha());
223 if (!CheckSector(sector)) continue;
224 fSectorTRD->Fill(sector);
225
226 fKinkIndex->Fill(track->GetKinkIndex(0));
227 if (track->GetKinkIndex(0)) continue;
228
229 UInt_t u = 1;
230 UInt_t status = track->GetStatus();
231 for(int bit=0; bit<32; bit++)
232 if (u<<bit & status) fTrackStatus->Fill(bit);
233
234 const int nbits = 6;
235 int bit[6] = {0,0,0,0,0,0};
236 bit[0] = status & AliESDtrack::kTPCin;
237 bit[1] = status & AliESDtrack::kTPCout;
238 bit[2] = (status & AliESDtrack::kTPCout) && !(status & AliESDtrack::kTRDout);
239 bit[3] = status & AliESDtrack::kTRDout;
240 bit[4] = status & AliESDtrack::kTRDrefit;
241 bit[5] = (status & AliESDtrack::kTRDout) && !(status & AliESDtrack::kTRDrefit);
242
243
244 // transverse momentum
245 const double *val = track->GetParameter(); // parameters at the vertex
246 double pt = 1./TMath::Abs(val[4]);
247
248 for(int b=0; b<nbits; b++) {
249 if (bit[b]) {
250 fPt[b]->Fill(pt);
251 fTheta[b]->Fill(val[3]);
252 fSigmaY[b]->Fill(TMath::Sqrt(paramOut->GetSigmaY2()));
253 fChi2[b]->Fill(track->GetTRDchi2()/track->GetTRDncls());
254 fPlaneYZ[b]->Fill(paramOut->GetY(), paramOut->GetZ());
255 }
256 }
257
258 // sectors
259 if (bit[1]) {
260 fAlpha[0]->Fill(paramIn->GetAlpha());
261 fAlpha[1]->Fill(paramOut->GetAlpha());
262 }
263
264 if (bit[3]) fAlpha[2]->Fill(paramOut->GetAlpha());
265 if (bit[4]) fAlpha[3]->Fill(paramOut->GetAlpha());
266
267 // clusters
268 for(int b=0; b<3; b++)
269 if (bit[3+b]) fClustersTRD[b]->Fill(track->GetTRDncls());
270
271 // refitted only
272 if (!bit[4]) continue;
273
274 fQuality->Fill(track->GetTRDQuality());
275 fBudget->Fill(track->GetTRDBudget());
276 fSignal->Fill(track->GetTRDsignal());
277
278 fTrdSigMom->Fill(track->GetP(), track->GetTRDsignal());
279 fTpcSigMom->Fill(track->GetP(), track->GetTPCsignal());
280
281 // PID only
282 if (status & AliESDtrack::kTRDpid) {
283
284 for(int l=0; l<6; l++) fTime->Fill(track->GetTRDTimBin(l));
285
286 // fill pid histograms
287 double trdr0 = 0, tpcr0 = 0;
288 int trdBestPid = 5, tpcBestPid = 5; // charged
289 const double minPidValue = 0.9;
290
291 double pp[5];
292 track->GetTPCpid(pp); // ESD inconsequence
293
294 for(int pid=0; pid<5; pid++) {
295
296 trdr0 += track->GetTRDpid(pid);
297 tpcr0 += pp[pid];
298
299 fTrdPID[pid]->Fill(track->GetTRDpid(pid));
300 fTpcPID[pid]->Fill(pp[pid]);
301
302 if (track->GetTRDpid(pid) > minPidValue) trdBestPid = pid;
303 if (pp[pid] > minPidValue) tpcBestPid = pid;
304 }
305
306 fTrdPID[5]->Fill(trdr0); // check unitarity
307 fTrdSigMomPID[trdBestPid]->Fill(track->GetP(), track->GetTRDsignal());
308
309 fTpcPID[5]->Fill(tpcr0); // check unitarity
310 fTpcSigMomPID[tpcBestPid]->Fill(track->GetP(), track->GetTPCsignal());
311 }
312
313 }
314
315 CalculateEff();
316 PostData(0, fOutputContainer);
317}
318
319//______________________________________________________________________________
320void AliTRDQATask::Terminate(Option_t *)
321{
322 // Processing when the event loop is ended
323 AliInfo("TRD QA module");
324
325 // create efficiency histograms
326
327 CalculateEff();
328 PostData(0, fOutputContainer);
329
330 DrawESD() ;
331 DrawGeoESD() ;
332 //DrawConvESD() ;
333 DrawPidESD() ;
334}
335
336//______________________________________________________________________________
337int AliTRDQATask::GetSector(double alpha)
338{
339 // Gets the sector number
340
341 double size = TMath::DegToRad() * 20.;
342 int sector = (int)((alpha + TMath::Pi())/size);
343 return sector;
344}
345
346//______________________________________________________________________________
347int AliTRDQATask::CheckSector(int sector)
348{
349 // Checks the sector number
350 const int nSec = 8;
351 int sec[] = {2,3,5,6,11,12,13,15};
352
353 for(int i=0; i<nSec; i++)
354 if (sector == sec[i]) return 1;
355
356 return 0;
357}
358
359//______________________________________________________________________________
360void AliTRDQATask::CalculateEff()
361{
362 // calculates the efficiency
363
364 for(int i=0; i<4; i++) fEffPt[i]->Reset();
365
366 fEffPt[0]->Add(fPt[1]);
367 fEffPt[0]->Divide(fPt[0]);
368
369 fEffPt[1]->Add(fPt[3]);
370 fEffPt[1]->Divide(fPt[1]);
371
372 fEffPt[2]->Add(fPt[4]);
373 fEffPt[2]->Divide(fPt[3]);
374
375 fEffPt[3]->Add(fPt[4]);
376 fEffPt[3]->Divide(fPt[1]);
377}
378
379//______________________________________________________________________________
380void AliTRDQATask::DrawESD()
381{
382 // Makes a few plots
383
384 AliInfo("Plotting....") ;
385
386 gROOT->SetStyle("Plain");
387 gStyle->SetPalette(1);
388 gStyle->SetOptStat(0);
389
390 TGaxis::SetMaxDigits(3);
391
392 gStyle->SetLabelFont(52, "XYZ");
393 gStyle->SetTitleFont(62, "XYZ");
394 gStyle->SetPadRightMargin(0.02);
395
396 // draw all
397
398 const int nplots = 18;
399 const int nover[nplots] = {1,1,1,4,1,1,1,1,1,1,2,1,1,3,1,1,1,1};
400 const int nnames = 24;
401 const char *names[nnames] = {
402 "ntracks", "kinkIndex", "trackStatus",
403 "ptTPCi", "ptTPCo", "ptTRDo", "ptTRDr", "ptTPCz", "ptTRDz",
404 "eff_TPCi_TPCo", "eff_TPCo_TRDo", "eff_TRDo_TRDr", "eff_TPCo_TRDr",
405 "clsTRDo", "clsTRDr", "clsTRDz",
406 "alphaTPCi", "alphaTPCo", "alphaTRDo", "alphaTRDr", "sectorTRD",
407 "time", "budget", "signal"
408 };
409
410 const int logy[nnames] = {
411 1,1,1,
412 1,1,1,
413 0,0,0,0,
414 1,1,
415 0,0,0,0,0,
416 0,1,1
417 };
418
419 int nhist=0;
420 for(int i=0; i<nplots; i++) {
421
422 new TCanvas(names[i], names[nhist], 500, 300);
423 gPad->SetLogy(logy[i]);
424
425 for(int j=0; j<nover[i]; j++) {
426 TH1D *hist = dynamic_cast<TH1D*>(gDirectory->FindObject(names[nhist++]));
427 if (!hist) continue;
428
429 if (strstr(hist->GetName(), "eff")) {
430 hist->SetMarkerStyle(20);
431 hist->SetMinimum(0);
432 hist->SetMaximum(1.2);
433 }
434
435 if (!j) hist->Draw();
436 else hist->Draw("SAME");
437 }
438
439 gPad->Print(Form("%s.gif", names[i]));
440 }
441}
442
443//______________________________________________________________________________
444void AliTRDQATask::DrawGeoESD()
445{
446 // Makes a few plots
447
448 AliInfo("Plotting....") ;
449
450 gStyle->SetOptStat(0);
451 TGaxis::SetMaxDigits(3);
452
453 gStyle->SetLabelFont(52, "XYZ");
454 gStyle->SetTitleFont(62, "XYZ");
455
456 gStyle->SetPadTopMargin(0.06);
457 gStyle->SetTitleBorderSize(0);
458
459 // draw all
460 const int nnames = 7;
461 const char *names[nnames] = {
462 "xIn", "xOut",
463 "planeYZTPCo", "planeYZTPCz", "planeYZTRDo", "planeYZTRDr", "planeYZTRDz",
464 };
465
466 const char *opt[nnames] = {
467 "", "",
468 "colz","colz", "colz", "colz", "colz"
469 };
470
471 const int logy[nnames] = {
472 1,1,
473 0,0,0,0,0
474 };
475
476 for(int i=0; i<nnames; i++) {
477
478 TH1D *hist = dynamic_cast<TH1D*>(gDirectory->FindObject(names[i]));
479 if (!hist) continue;
480
481 if (i<2) new TCanvas(names[i], names[i], 500, 300);
482 else new TCanvas(names[i], names[i], 300, 900);
483
484 gPad->SetLogy(logy[i]);
485 if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1);
486
487 hist->Draw(opt[i]);
488 AliInfo(Form("%s\t%d", names[i], hist->GetEntries()));
489
490 gPad->Print(Form("%s.gif", names[i]));
491 }
492}
493
494//______________________________________________________________________________
495void AliTRDQATask::DrawConvESD()
496{
497 // Makes a few plots
498
499 AliInfo("Plotting....") ;
500
501 gROOT->SetStyle("Plain");
502 gROOT->ForceStyle();
503 gStyle->SetPalette(1);
504
505 TGaxis::SetMaxDigits(3);
506
507 gStyle->SetLabelFont(52, "XYZ");
508 gStyle->SetTitleFont(62, "XYZ");
509 gStyle->SetPadRightMargin(0.02);
510
511 const int nnames = 9;
512 const int nplots = 5;
513 const int nover[nplots] = {3,1,1,3,1};
514
515 const char *names[nnames] = {
516 "sigmaYTPCo","sigmaYTRDo", "sigmaYTRDr", "sigmaYTPCz", "sigmaYTRDz",
517 "Chi2TPCo", "Chi2TRDo", "Chi2TRDr", "Chi2TRDz"
518 };
519
520 const char *opt[nplots] = {
521 "", "", "","","",
522 };
523
524 const int logy[nplots] = {
525 0,0,0,1,1
526 };
527
528 int nhist = 0;
529 for(int i=0; i<nplots; i++) {
530
531 new TCanvas(names[i], names[i], 500, 300);
532 gPad->SetLogy(logy[i]);
533 if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1);
534
535 for(int j=0; j<nover[i]; j++) {
536 TH1D *hist = dynamic_cast<TH1D*>(gDirectory->FindObject(names[nhist++]));
537 if (!j) hist->Draw(opt[i]);
538 else hist->Draw("same");
539 }
540
541 gPad->Print(Form("%s.eps", names[i]));
542 }
543}
544
545//______________________________________________________________________________
546void AliTRDQATask::DrawPidESD()
547{
548 // Makes a few plots
549
550 AliInfo("Plotting....") ;
551
552 gROOT->SetStyle("Plain");
553 gROOT->ForceStyle();
554 gStyle->SetPalette(1);
555 gStyle->SetOptStat(0);
556
557 TGaxis::SetMaxDigits(3);
558
559 gStyle->SetLabelFont(52, "XYZ");
560 gStyle->SetTitleFont(62, "XYZ");
561
562 gStyle->SetPadTopMargin(0.05);
563 gStyle->SetPadRightMargin(0.02);
564
565 // draw all
566
567 const int nnames = 27;
568 const char *names[nnames] = {
569
570 "signal", "trdSigMom", "tpcSigMom",
571
572 "trdPidEl", "trdPidMu", "trdPidPi", "trdPidK", "trdPidP", "trdPidCh",
573 "trdSigMomEl", "trdSigMomMu", "trdSigMomPi", "trdSigMomK", "trdSigMomP", "trdSigMomCh",
574
575 "tpcPidEl", "tpcPidMu", "tpcPidPi", "tpcPidK", "tpcPidP", "tpcPidCh",
576 "tpcSigMomEl", "tpcSigMomMu", "tpcSigMomPi", "tpcSigMomK", "tpcSigMomP", "tpcSigMomCh"
577
578 };
579
580 const char *opt[nnames] = {
581
582 "", "colz", "colz",
583
584 "", "", "", "", "", "" ,
585 "colz", "colz", "colz", "colz", "colz", "colz",
586
587 "", "", "", "", "", "" ,
588 "colz", "colz", "colz", "colz", "colz", "colz"
589 };
590
591 const int logy[nnames] = {
592
593 0,0,0,
594
595 1,1,1,1,1,1,
596 0,0,0,0,0,0,
597
598 1,1,1,1,1,1,
599 0,0,0,0,0,0
600 };
601
602 for(int i=0; i<nnames; i++) {
603
604 TH1D *hist = dynamic_cast<TH1D*>(gDirectory->FindObject(names[i]));
605 if (!hist) continue;
606
607 new TCanvas(names[i], names[i], 500, 300);
608 gPad->SetLogy(logy[i]);
609 if (strstr(opt[i],"colz")) gPad->SetRightMargin(0.1);
610
611 if (strstr(names[i],"sigMom")) gPad->SetLogz(1);
612 if (strstr(names[i],"SigMom")) gPad->SetLogz(1);
613
614 hist->Draw(opt[i]);
615 AliInfo(Form("%s\t%d", names[i], hist->GetEntries()));
616
617 gPad->Print(Form("%s.gif", names[i]));
618 }
619
620}