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