1)AliAODParticleCorrelation.h: Data member fIsolation changed from float to bool
[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
16
17#include <TClonesArray.h>
18#include <TObjArray.h>
19#include <TObject.h>
20#include <TH2I.h>
48797466 21#include <TGraphErrors.h>
9a281e83 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
44ClassImp(AliTRDcheckESD)
45
46const Float_t AliTRDcheckESD::xTPC = 290.;
47const Float_t AliTRDcheckESD::xTOF = 365.;
48
49//____________________________________________________________________
50AliTRDcheckESD::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//____________________________________________________________________
66AliTRDcheckESD::~AliTRDcheckESD()
67{
68 if(fHistos){
69 //fHistos->Delete();
70 delete fHistos;
71 }
72}
73
74//____________________________________________________________________
75void 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//____________________________________________________________________
92void 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);
48797466 119
120 TObjArray *res = new TObjArray();
121 res->SetName("Results");
122 fHistos->AddAt(res, kResults);
9a281e83 123}
124
125//____________________________________________________________________
126void AliTRDcheckESD::Exec(Option_t *){
127 //
128 // Run the Analysis
129 //
130 if(!fESD){
48797466 131 AliError("ESD event missing.");
9a281e83 132 return;
133 }
134
135 // Get MC information if available
136 AliStack * fStack = 0x0;
137 if(HasMC() && !fMC){
48797466 138 AliWarning("MC event missing");
9a281e83 139 SetMC(kFALSE);
140 } else {
141 if(!(fStack = fMC->Stack())){
48797466 142 AliWarning("MC stack missing");
9a281e83 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
48797466 159 //ULong_t status = esdTrack->GetStatus();
9a281e83 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)))){
48797466 188 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 189 continue;
190 }
191
192 AliTrackReference *ref = 0x0;
193 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 194 if(!nRefs){
195 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
196 continue;
197 }
9a281e83 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();
48797466 207 //Int_t fPdg = tParticle->GetPdgCode();
9a281e83 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 [");
48797466 216 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 217 } else {
218 //printf("%7.2f [", ref->LocalX());
48797466 219 if(ref->LocalX() < 300. && tParticle->IsPrimary()){
220 TRDin=1;
221 if(esdTrack->GetNcls(2)) TRDout=1;
222 if(esdTrack->GetTRDpidQuality()) TRDpid=1;
223 }
9a281e83 224 }
225 } else {
226 //printf(" TPC [");
48797466 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");
9a281e83 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//____________________________________________________________________
252void AliTRDcheckESD::Terminate(Option_t *)
253{
48797466 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 }
9a281e83 308}