warnings and few violations corrected
[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"
d12237d6 31#include "AliESDkink.h"
9a281e83 32#include "AliMCEvent.h"
33#include "AliESDInputHandler.h"
34#include "AliMCEventHandler.h"
35
36#include "AliESDtrack.h"
37#include "AliMCParticle.h"
38#include "AliPID.h"
39#include "AliStack.h"
40#include "AliTrackReference.h"
9a281e83 41
42#include "AliTRDcheckESD.h"
43
44ClassImp(AliTRDcheckESD)
45
965e895b 46const Int_t AliTRDcheckESD::fgkNgraphs = 4;
47const Float_t AliTRDcheckESD::fgkxTPC = 290.;
48const Float_t AliTRDcheckESD::fgkxTOF = 365.;
59589fc9 49FILE* AliTRDcheckESD::fgFile = 0x0;
9a281e83 50
51//____________________________________________________________________
52AliTRDcheckESD::AliTRDcheckESD():
6da3eee3 53 AliAnalysisTask("checkESD", "ESD checker for TRD info")
9a281e83 54 ,fStatus(0)
55 ,fESD(0x0)
56 ,fMC(0x0)
57 ,fHistos(0x0)
58{
59 //
60 // Default constructor
61 //
49c7f0a5 62 SetMC(kTRUE);
9a281e83 63 DefineInput(0, TChain::Class());
64 DefineOutput(0, TObjArray::Class());
65}
66
67//____________________________________________________________________
68AliTRDcheckESD::~AliTRDcheckESD()
69{
70 if(fHistos){
71 //fHistos->Delete();
72 delete fHistos;
73 }
74}
75
76//____________________________________________________________________
77void AliTRDcheckESD::ConnectInputData(Option_t *)
78{
79 //
80 // Link the Input Data
81 //
82 TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
83 if(tree) tree->SetBranchStatus("Tracks", 1);
84
85 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
86 fESD = esdH ? esdH->GetEvent() : 0x0;
87
88 if(!HasMC()) return;
89 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
90 fMC = mcH ? mcH->MCEvent() : 0x0;
91}
92
93//____________________________________________________________________
94void AliTRDcheckESD::CreateOutputObjects()
95{
96 //
97 // Create Output Containers (TObjectArray containing 1D histograms)
98 //
99 OpenFile(0, "RECREATE");
c2fd5953 100 Histos();
9a281e83 101}
102
965e895b 103//____________________________________________________________________
ea4502f6 104TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
965e895b 105{
ea4502f6 106 Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
107 kCLEAR = strstr(opt, "c"); // clear existing graph
965e895b 108
109 const Char_t *name[] = {
110 "Geo", "Trk", "Pid", "Ref"
111 };
112 const Char_t *title[] = {
113 "TRD geometrical efficiency (TRDin/TPCout)"
114 ,"TRD tracking efficiency (TRDout/TRDin)"
115 ,"TRD PID efficiency (TRDpid/TRDin)"
116 ,"TRD refit efficiency (TRDrefit/TRDin)"
117 };
118 const Int_t ngr = sizeof(name)/sizeof(Char_t*);
119 if(ngr != fgkNgraphs){
120 AliWarning("No of graphs defined different from definition");
121 return 0x0;
122 }
123
124 TObjArray *res = 0x0;
125 if(!(res = (TObjArray*)fHistos->At(kResults)) ||
126 (id < 0 || id >= ngr)){
ea4502f6 127 AliWarning("Graph array missing.");
965e895b 128 return 0x0;
129 }
130
131 TGraphErrors *g = 0x0;
132 if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
59589fc9 133 if(kCLEAR){
134 for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
135 } else {
136 PutTrendValue(name[id], g->GetMean(2), g->GetRMS(2));
137 }
965e895b 138 } else {
139 if(kBUILD){
140 g = new TGraphErrors();
141 g->SetNameTitle(name[id], title[id]);
142 res->AddAt(g, id);
143 }
144 }
145 return g;
146}
147
9a281e83 148//____________________________________________________________________
149void AliTRDcheckESD::Exec(Option_t *){
150 //
151 // Run the Analysis
152 //
153 if(!fESD){
48797466 154 AliError("ESD event missing.");
9a281e83 155 return;
156 }
157
158 // Get MC information if available
159 AliStack * fStack = 0x0;
2b104c26 160 if(HasMC()){
161 if(!fMC){
162 AliWarning("MC event missing");
9a281e83 163 SetMC(kFALSE);
2b104c26 164 } else {
165 if(!(fStack = fMC->Stack())){
166 AliWarning("MC stack missing");
167 SetMC(kFALSE);
168 }
9a281e83 169 }
170 }
d12237d6 171 Bool_t TRDin(0), TRDout(0), TRDpid(0);
9a281e83 172
9a281e83 173 AliESDtrack *esdTrack = 0x0;
174 for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
d12237d6 175 TRDin=0;TRDout=0;TRDpid=0;
9a281e83 176 esdTrack = fESD->GetTrack(itrk);
d12237d6 177
178// if(esdTrack->GetNcls(1)) nTPC++;
179// if(esdTrack->GetNcls(2)) nTRD++;
9a281e83 180
181 // track status
d12237d6 182 ULong_t status = esdTrack->GetStatus();
2b104c26 183 //PrintStatus(status);
d12237d6 184
185 // define TPC out tracks
186 if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
187 if(esdTrack->GetKinkIndex(0) > 0) continue;
9a281e83 188
189 // TRD PID
190 Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
191 // pid quality
ed15ef4f 192 //esdTrack->GetTRDntrackletsPID();
9a281e83 193
194 // look at external track param
195 const AliExternalTrackParam *op = esdTrack->GetOuterParam();
196 Double_t xyz[3];
197 if(op){
198 op->GetXYZ(xyz);
199 op->Global2LocalPosition(xyz, op->GetAlpha());
200 //printf("op @ X[%7.3f]\n", xyz[0]);
201 }
202
203 // read MC info
204 if(!HasMC()) continue;
205
206 Int_t fLabel = esdTrack->GetLabel();
207 if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue;
208
209 // read MC particle
210 AliMCParticle *mcParticle = 0x0;
7aad0c47 211 if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
48797466 212 AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
9a281e83 213 continue;
214 }
215
216 AliTrackReference *ref = 0x0;
217 Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
48797466 218 if(!nRefs){
219 AliWarning(Form("Track refs missing. Label[%d].", fLabel));
220 continue;
221 }
9a281e83 222 Int_t iref = 0;
223 while(iref<nRefs){
224 ref = mcParticle->GetTrackReference(iref);
965e895b 225 if(ref->LocalX() > fgkxTPC) break;
9a281e83 226 ref=0x0; iref++;
227 }
228
229 // read TParticle
965e895b 230 //TParticle *tParticle = mcParticle->Particle();
48797466 231 //Int_t fPdg = tParticle->GetPdgCode();
d12237d6 232 // reject secondaries
233 //if(!tParticle->IsPrimary()) continue;
9a281e83 234
d12237d6 235 if(ref){
965e895b 236 if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
48797466 237 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
9a281e83 238 } else {
d12237d6 239 TRDin=1;
240 if(esdTrack->GetNcls(2)) TRDout=1;
ed15ef4f 241 if(esdTrack->GetTRDntrackletsPID()>=4) TRDpid=1;
9a281e83 242 }
d12237d6 243 } else { // track stopped in TPC
48797466 244 ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
245 }
d12237d6 246 // get the MC pt !!
9a281e83 247 Float_t pt = ref->Pt();
d12237d6 248
9a281e83 249 TH2 *h = (TH2I*)fHistos->At(kTRDstat);
d12237d6 250 if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
251 if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
9a281e83 252 if(/*status & AliESDtrack::k*/TRDout){
253 ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
d12237d6 254 h->Fill(pt, kTRDout);
9a281e83 255 }
965e895b 256 if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
257 if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
9a281e83 258 }
9a281e83 259 PostData(0, fHistos);
260}
261
c2fd5953 262//____________________________________________________________________
263TObjArray* AliTRDcheckESD::Histos()
264{
265 if(fHistos) return fHistos;
266
267 fHistos = new TObjArray(kNhistos);
268 //fHistos->SetOwner(kTRUE);
269
270 TH1 *h = 0x0;
271
272 // clusters per tracklet
273 if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
274 h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
275 h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
276 h->GetYaxis()->SetTitle("entries");
277 } else h->Reset();
278 fHistos->AddAt(h, kNCl);
279
280 // status bits histogram
281 if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
282 h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
283 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
284 h->GetYaxis()->SetTitle("status bits");
285 h->GetZaxis()->SetTitle("entries");
286 } else h->Reset();
287 fHistos->AddAt(h, kTRDstat);
288
289 // results array
290 TObjArray *res = new TObjArray();
291 res->SetName("Results");
292 fHistos->AddAt(res, kResults);
293 return fHistos;
294}
295
965e895b 296//____________________________________________________________________
297Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
298{
299 if(!TFile::Open(filename)){
300 AliWarning(Form("Couldn't open file %s.", filename));
301 return kFALSE;
302 }
303 TObjArray *o = 0x0;
304 if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
305 AliWarning("Missing histogram container.");
306 return kFALSE;
307 }
308 fHistos = (TObjArray*)o->Clone(GetName());
309 gFile->Close();
310 return kTRUE;
311}
9a281e83 312
59589fc9 313//_______________________________________________________
314Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val, Double_t err)
315{
316 if(!fgFile){
317 fgFile = fopen("TRD.Performance.txt", "at");
318 }
319 fprintf(fgFile, "%s_%s %f %f\n", GetName(), name, val, err);
320 return kTRUE;
321}
322
9a281e83 323//____________________________________________________________________
324void AliTRDcheckESD::Terminate(Option_t *)
325{
965e895b 326 TObjArray *res = 0x0;
327 if(!(res = (TObjArray*)fHistos->At(kResults))){
328 AliWarning("Graph container missing.");
329 return;
330 }
331 if(!res->GetEntriesFast()) res->Expand(fgkNgraphs);
48797466 332
333 // geometrical efficiency
965e895b 334 TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
335 TH1 *h1[2] = {0x0, 0x0};
d12237d6 336 h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
337 h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
965e895b 338 Process(h1, GetGraph(0));
ea4502f6 339 delete h1[0];delete h1[1];
48797466 340
341 // tracking efficiency
d12237d6 342 h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
343 h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
965e895b 344 Process(h1, GetGraph(1));
ea4502f6 345 delete h1[1];
48797466 346
347 // PID efficiency
965e895b 348 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
d12237d6 349 h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
965e895b 350 Process(h1, GetGraph(2));
ea4502f6 351 delete h1[1];
965e895b 352
353 // Refit efficiency
354 //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
355 h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
356 Process(h1, GetGraph(3));
ea4502f6 357 delete h1[1];
965e895b 358}
359
360//____________________________________________________________________
361void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
362{
363 Int_t n1 = 0, n2 = 0, ip=0;
364 Double_t eff = 0.;
365
366 TAxis *ax = h1[0]->GetXaxis();
48797466 367 for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
368 if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
369 n2 = (Int_t)h1[1]->GetBinContent(ib);
370 eff = n2/Float_t(n1);
371
372 ip=g->GetN();
373 g->SetPoint(ip, ax->GetBinCenter(ib), eff);
374 g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
375 }
2b104c26 376}
377
378//____________________________________________________________________
379void AliTRDcheckESD::PrintStatus(ULong_t status)
380{
381 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"
382 ,Bool_t(status & AliESDtrack::kITSin)
383 ,Bool_t(status & AliESDtrack::kITSout)
384 ,Bool_t(status & AliESDtrack::kITSrefit)
385 ,Bool_t(status & AliESDtrack::kTPCin)
386 ,Bool_t(status & AliESDtrack::kTPCout)
387 ,Bool_t(status & AliESDtrack::kTPCrefit)
388 ,Bool_t(status & AliESDtrack::kTPCpid)
389 ,Bool_t(status & AliESDtrack::kTRDin)
390 ,Bool_t(status & AliESDtrack::kTRDout)
391 ,Bool_t(status & AliESDtrack::kTRDrefit)
392 ,Bool_t(status & AliESDtrack::kTRDpid)
393 ,Bool_t(status & AliESDtrack::kTRDStop)
394 ,Bool_t(status & AliESDtrack::kHMPIDout)
395 ,Bool_t(status & AliESDtrack::kHMPIDpid)
396 );
9a281e83 397}
2b104c26 398