- add output graphs for efficiency
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckESD.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
17 #include <TClonesArray.h>
18 #include <TObjArray.h>
19 #include <TObject.h>
20 #include <TH2I.h>
21 #include <TGraphErrors.h>
22 #include <TFile.h>
23 #include <TTree.h>
24 #include <TROOT.h>
25 #include <TChain.h>
26 #include <TParticle.h>
27
28 #include "AliLog.h"
29 #include "AliAnalysisManager.h"
30 #include "AliESDEvent.h"
31 #include "AliMCEvent.h"
32 #include "AliESDInputHandler.h"
33 #include "AliMCEventHandler.h"
34
35 #include "AliESDtrack.h"
36 #include "AliMCParticle.h"
37 #include "AliPID.h"
38 #include "AliStack.h"
39 #include "AliTrackReference.h"
40 #include "AliTRDgeometry.h"
41
42 #include "AliTRDcheckESD.h"
43
44 ClassImp(AliTRDcheckESD)
45
46 const Float_t AliTRDcheckESD::xTPC = 290.;
47 const Float_t AliTRDcheckESD::xTOF = 365.;
48
49 //____________________________________________________________________
50 AliTRDcheckESD::AliTRDcheckESD():
51   AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
52   ,fStatus(0)
53   ,fESD(0x0)
54   ,fMC(0x0)
55   ,fHistos(0x0)
56 {
57   //
58   // Default constructor
59   //
60
61   DefineInput(0, TChain::Class());
62   DefineOutput(0, TObjArray::Class());
63 }
64
65 //____________________________________________________________________
66 AliTRDcheckESD::~AliTRDcheckESD()
67 {
68   if(fHistos){
69     //fHistos->Delete();
70     delete fHistos;
71   }
72 }
73
74 //____________________________________________________________________
75 void AliTRDcheckESD::ConnectInputData(Option_t *)
76 {
77   //
78   // Link the Input Data
79   //
80   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
81   if(tree) tree->SetBranchStatus("Tracks", 1);
82
83   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
84   fESD = esdH ? esdH->GetEvent() : 0x0;
85
86   if(!HasMC()) return;
87   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
88   fMC = mcH ? mcH->MCEvent() : 0x0;
89 }
90
91 //____________________________________________________________________
92 void AliTRDcheckESD::CreateOutputObjects()
93 {       
94   //
95   // Create Output Containers (TObjectArray containing 1D histograms)
96   //
97   OpenFile(0, "RECREATE");  
98   fHistos = new TObjArray(5);
99   //fHistos->SetOwner(kTRUE);
100   
101   TH1 *h = 0x0;
102
103   // clusters per tracklet
104   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
105     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
106     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
107     h->GetYaxis()->SetTitle("entries");
108   } else h->Reset();
109   fHistos->AddAt(h, kNCl);
110
111   // TPC out
112   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
113     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., 4, -.5, 3.5);
114     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
115     h->GetYaxis()->SetTitle("status bits");
116     h->GetZaxis()->SetTitle("entries");
117   } else h->Reset();
118   fHistos->AddAt(h, kTRDstat);
119
120   TObjArray *res = new TObjArray();
121   res->SetName("Results");
122   fHistos->AddAt(res, kResults);
123 }
124
125 //____________________________________________________________________
126 void AliTRDcheckESD::Exec(Option_t *){
127   //
128   // Run the Analysis
129   //
130   if(!fESD){
131     AliError("ESD event missing.");
132     return;
133   }
134
135   // Get MC information if available
136   AliStack * fStack = 0x0;
137   if(HasMC() && !fMC){ 
138     AliWarning("MC event missing");
139     SetMC(kFALSE);
140   } else {
141     if(!(fStack = fMC->Stack())){
142       AliWarning("MC stack missing");
143       SetMC(kFALSE);
144     }
145   }
146   Bool_t TPCout(0), TRDin(0), TRDout(0), TRDpid(0);
147
148
149   Int_t nTRD = 0, nTPC = 0;
150   //Int_t nTracks = fESD->GetNumberOfTracks();
151   AliESDtrack *esdTrack = 0x0;
152   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
153     TPCout=0;TRDin=0;TRDout=0;TRDpid=0;
154     esdTrack = fESD->GetTrack(itrk);
155     if(esdTrack->GetNcls(1)) nTPC++;
156     if(esdTrack->GetNcls(2)) nTRD++;
157
158     // track status
159     //ULong_t status = esdTrack->GetStatus();
160
161     // TRD PID
162     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
163     // pid quality
164     esdTrack->GetTRDpidQuality();
165     // kink index
166     esdTrack->GetKinkIndex(0);
167     // TPC clusters
168     esdTrack->GetNcls(1);
169
170     // look at external track param
171     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
172     Double_t xyz[3];
173     if(op){
174       op->GetXYZ(xyz);
175       op->Global2LocalPosition(xyz, op->GetAlpha());
176       //printf("op @ X[%7.3f]\n", xyz[0]);
177     }
178
179     // read MC info
180     if(!HasMC()) continue;
181
182     Int_t fLabel = esdTrack->GetLabel();
183     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue; 
184     
185     // read MC particle
186     AliMCParticle *mcParticle = 0x0; 
187     if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
188       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
189       continue;
190     }
191
192     AliTrackReference *ref = 0x0; 
193     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
194     if(!nRefs){
195       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
196       continue;
197     }
198     Int_t iref = 0;
199     while(iref<nRefs){
200       ref = mcParticle->GetTrackReference(iref);
201       if(ref->LocalX() > xTPC) break;
202       ref=0x0; iref++;
203     }
204
205     // read TParticle
206     TParticle *tParticle = mcParticle->Particle(); 
207     //Int_t fPdg = tParticle->GetPdgCode();
208     //tParticle->IsPrimary();
209
210     //printf("[%c] ref[%2d]=", tParticle->IsPrimary() ? 'P' : 'S', iref);
211
212     TPCout=1;
213     if(ref){
214       if(ref->LocalX() > xTOF){ 
215         //printf("  TOF   [");
216         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
217       } else {
218         //printf("%7.2f [", ref->LocalX());
219         if(ref->LocalX() < 300. && tParticle->IsPrimary()){
220           TRDin=1;
221           if(esdTrack->GetNcls(2)) TRDout=1;
222           if(esdTrack->GetTRDpidQuality()) TRDpid=1;
223         }
224       }
225     } else { 
226       //printf("  TPC   [");
227       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
228     }
229     if(!ref){
230       printf("[%c] ref[%2d] [", tParticle->IsPrimary() ? 'P' : 'S', iref);
231       for(Int_t ir=0; ir<nRefs; ir++) printf(" %6.2f", mcParticle->GetTrackReference(ir)->LocalX());
232       printf("]\n");
233     }
234     Float_t pt = ref->Pt();
235     //printf("%f]\n", pt);
236     
237     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
238     if(/*status & AliESDtrack::k*/TPCout) h->Fill(pt, 0);
239     if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, 1);
240     if(/*status & AliESDtrack::k*/TRDout){ 
241       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
242       h->Fill(pt, 2);
243     }
244     if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, 3);
245   }  
246
247   PostData(0, fHistos);
248 }
249
250
251 //____________________________________________________________________
252 void AliTRDcheckESD::Terminate(Option_t *)
253 {
254   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
255   TAxis *ax = h2->GetXaxis();
256   TObjArray *res = (TObjArray*)fHistos->At(kResults);
257   
258   TH1 *h1[2] = {0x0, 0x0};
259   TGraphErrors *g = 0x0;
260   res->Expand(1);
261   Int_t n1 = 0, n2 = 0, ip=0;
262   Double_t eff = 0.;
263   
264   // geometrical efficiency
265   h1[0] = h2->ProjectionX("px0", 1, 1);
266   h1[1] = h2->ProjectionX("px1", 2, 2);
267   res->Add(g = new TGraphErrors());
268   g->SetNameTitle("geom", "TRD geometrical efficiency (TRDin/TPCout)");
269   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
270     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
271     n2 = (Int_t)h1[1]->GetBinContent(ib);
272     eff = n2/Float_t(n1);
273
274     ip=g->GetN();
275     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
276     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
277   }
278
279   // tracking efficiency
280   h1[0] = h2->ProjectionX("px0", 2, 2);
281   h1[1] = h2->ProjectionX("px1", 3, 3);
282   res->Add(g = new TGraphErrors());
283   g->SetNameTitle("tracking", "TRD tracking efficiency (TRDout/TRDin)");
284   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
285     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
286     n2 = (Int_t)h1[1]->GetBinContent(ib);
287     eff = n2/Float_t(n1);
288
289     ip=g->GetN();
290     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
291     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
292   }
293
294   // PID efficiency
295   h1[0] = h2->ProjectionX("px0", 2, 2);
296   h1[1] = h2->ProjectionX("px1", 4, 4);
297   res->Add(g = new TGraphErrors());
298   g->SetNameTitle("PID", "TRD PID efficiency (TRDpid/TRDin)");
299   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
300     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
301     n2 = (Int_t)h1[1]->GetBinContent(ib);
302     eff = n2/Float_t(n1);
303
304     ip=g->GetN();
305     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
306     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
307   }
308 }