trend plots for DET and ESD performance tasks
[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 "AliESDkink.h"
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"
41
42 #include "AliTRDcheckESD.h"
43
44 ClassImp(AliTRDcheckESD)
45
46 const Int_t   AliTRDcheckESD::fgkNgraphs = 4;
47 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
48 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
49 FILE* AliTRDcheckESD::fgFile = 0x0;
50
51 //____________________________________________________________________
52 AliTRDcheckESD::AliTRDcheckESD():
53   AliAnalysisTask("checkESD", "ESD checker for TRD info")
54   ,fStatus(0)
55   ,fESD(0x0)
56   ,fMC(0x0)
57   ,fHistos(0x0)
58 {
59   //
60   // Default constructor
61   //
62   SetMC(kTRUE);
63   DefineInput(0, TChain::Class());
64   DefineOutput(0, TObjArray::Class());
65 }
66
67 //____________________________________________________________________
68 AliTRDcheckESD::~AliTRDcheckESD()
69 {
70   if(fHistos){
71     //fHistos->Delete();
72     delete fHistos;
73   }
74 }
75
76 //____________________________________________________________________
77 void 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 //____________________________________________________________________
94 void AliTRDcheckESD::CreateOutputObjects()
95 {       
96   //
97   // Create Output Containers (TObjectArray containing 1D histograms)
98   //
99   OpenFile(0, "RECREATE");  
100   Histos();
101 }
102
103 //____________________________________________________________________
104 TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
105 {
106   Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
107          kCLEAR = strstr(opt, "c"); // clear existing graph
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)){
127     AliWarning("Graph array missing.");
128     return 0x0;
129   }
130
131   TGraphErrors *g = 0x0;
132   if((g = dynamic_cast<TGraphErrors*>(res->At(id)))){
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     }
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
148 //____________________________________________________________________
149 void AliTRDcheckESD::Exec(Option_t *){
150   //
151   // Run the Analysis
152   //
153   if(!fESD){
154     AliError("ESD event missing.");
155     return;
156   }
157
158   // Get MC information if available
159   AliStack * fStack = 0x0;
160   if(HasMC()){
161     if(!fMC){ 
162       AliWarning("MC event missing");
163       SetMC(kFALSE);
164     } else {
165       if(!(fStack = fMC->Stack())){
166         AliWarning("MC stack missing");
167         SetMC(kFALSE);
168       }
169     }
170   }
171   Bool_t TRDin(0), TRDout(0), TRDpid(0);
172
173   AliESDtrack *esdTrack = 0x0;
174   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
175     TRDin=0;TRDout=0;TRDpid=0;
176     esdTrack = fESD->GetTrack(itrk);
177
178 //     if(esdTrack->GetNcls(1)) nTPC++;
179 //     if(esdTrack->GetNcls(2)) nTRD++;
180
181     // track status
182     ULong_t status = esdTrack->GetStatus();
183     //PrintStatus(status);
184
185     // define TPC out tracks
186     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
187     if(esdTrack->GetKinkIndex(0) > 0) continue;
188
189     // TRD PID
190     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
191     // pid quality
192     //esdTrack->GetTRDntrackletsPID();
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; 
211     if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
212       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
213       continue;
214     }
215
216     AliTrackReference *ref = 0x0; 
217     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
218     if(!nRefs){
219       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
220       continue;
221     }
222     Int_t iref = 0;
223     while(iref<nRefs){
224       ref = mcParticle->GetTrackReference(iref);
225       if(ref->LocalX() > fgkxTPC) break;
226       ref=0x0; iref++;
227     }
228
229     // read TParticle
230     //TParticle *tParticle = mcParticle->Particle(); 
231     //Int_t fPdg = tParticle->GetPdgCode();
232     // reject secondaries
233     //if(!tParticle->IsPrimary()) continue;
234
235     if(ref){ 
236       if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
237         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
238       } else {
239         TRDin=1;
240         if(esdTrack->GetNcls(2)) TRDout=1;
241         if(esdTrack->GetTRDntrackletsPID()>=4) TRDpid=1;
242       }
243     } else { // track stopped in TPC 
244       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
245     }
246     // get the MC pt !!
247     Float_t pt = ref->Pt();
248
249     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
250     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
251     if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
252     if(/*status & AliESDtrack::k*/TRDout){ 
253       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
254       h->Fill(pt, kTRDout);
255     }
256     if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, kTRDpid);
257     if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
258   }  
259   PostData(0, fHistos);
260 }
261
262 //____________________________________________________________________
263 TObjArray* 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
296 //____________________________________________________________________
297 Bool_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 }
312
313 //_______________________________________________________
314 Bool_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
323 //____________________________________________________________________
324 void AliTRDcheckESD::Terminate(Option_t *)
325 {
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);
332   
333   // geometrical efficiency
334   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
335   TH1 *h1[2] = {0x0, 0x0};
336   h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
337   h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
338   Process(h1, GetGraph(0));
339   delete h1[0];delete h1[1];
340
341   // tracking efficiency
342   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
343   h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
344   Process(h1, GetGraph(1));
345   delete h1[1];
346
347   // PID efficiency
348   //h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
349   h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
350   Process(h1, GetGraph(2));
351   delete h1[1];
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));
357   delete h1[1];
358 }
359
360 //____________________________________________________________________
361 void 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();
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   }
376 }  
377
378 //____________________________________________________________________
379 void 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   );
397 }
398