New TRD QA task
[u/mrichter/AliRoot.git] / ESDCheck / AliTRDQATask.cxx
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 //______________________________________________________________________________
44 AliTRDQATask::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 //______________________________________________________________________________
57 void 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 //______________________________________________________________________________
188 void 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 //______________________________________________________________________________
320 void 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 //______________________________________________________________________________
337 int 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 //______________________________________________________________________________
347 int 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 //______________________________________________________________________________
360 void 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 //______________________________________________________________________________
380 void 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 //______________________________________________________________________________
444 void 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 //______________________________________________________________________________
495 void 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 //______________________________________________________________________________
546 void 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 }