fix compilation
[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 <TPad.h>
34 #include <TLegend.h>
35 #include <TGraphErrors.h>
36 #include <TGraphAsymmErrors.h>
37 #include <TFile.h>
38 #include <TTree.h>
39 #include <TROOT.h>
40 #include <TChain.h>
41 #include <TParticle.h>
42
43 #include "AliLog.h"
44 #include "AliAnalysisManager.h"
45 #include "AliESDEvent.h"
46 #include "AliESDkink.h"
47 #include "AliMCEvent.h"
48 #include "AliESDInputHandler.h"
49 #include "AliMCEventHandler.h"
50
51 #include "AliESDtrack.h"
52 #include "AliMCParticle.h"
53 #include "AliPID.h"
54 #include "AliStack.h"
55 #include "AliTrackReference.h"
56
57 #include "AliTRDcheckESD.h"
58
59 ClassImp(AliTRDcheckESD)
60
61 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
62 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
63 FILE* AliTRDcheckESD::fgFile = 0x0;
64
65 //____________________________________________________________________
66 AliTRDcheckESD::AliTRDcheckESD():
67   AliAnalysisTask("checkESD", "ESD checker for TRD info")
68   ,fStatus(0)
69   ,fESD(0x0)
70   ,fMC(0x0)
71   ,fHistos(0x0)
72   ,fResults(0x0)
73 {
74   //
75   // Default constructor
76   //
77   SetMC(kTRUE);
78   DefineInput(0, TChain::Class());
79   DefineOutput(0, TObjArray::Class());
80 }
81
82 //____________________________________________________________________
83 AliTRDcheckESD::~AliTRDcheckESD()
84 {
85 // Destructor
86   if(fHistos){
87     //fHistos->Delete();
88     delete fHistos;
89   }
90   if(fResults){
91     fResults->Delete();
92     delete fResults;
93   }
94 }
95
96 //____________________________________________________________________
97 void AliTRDcheckESD::ConnectInputData(Option_t *)
98 {
99   //
100   // Link the Input Data
101   //
102   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
103   if(tree) tree->SetBranchStatus("Tracks", 1);
104
105   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
106   fESD = esdH ? esdH->GetEvent() : 0x0;
107
108   if(!HasMC()) return;
109   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
110   fMC = mcH ? mcH->MCEvent() : 0x0;
111 }
112
113 //____________________________________________________________________
114 void AliTRDcheckESD::CreateOutputObjects()
115 {       
116   //
117   // Create Output Containers (TObjectArray containing 1D histograms)
118   //
119   OpenFile(0, "RECREATE");  
120   Histos();
121 }
122
123 //____________________________________________________________________
124 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
125 {
126   if(!gPad){
127     AliWarning("Please provide a canvas to draw results.");
128     return kFALSE;
129   }
130   TGraph *g(0x0); TH1 *h(0x0); TLegend *leg(0x0);
131   switch(ifig){
132   case kNCl:
133     if(!fHistos || !(h=(TH1I*)fHistos->At(kNCl))) break;
134     h->Draw("c");
135     return kTRUE;
136   case kTRDstat:
137     leg=new TLegend(.5, .75, .98, .98);
138     leg->SetBorderSize(1); leg->SetFillColor(0);
139     leg->SetHeader("Tracking Efficiency");
140     for(Int_t ieff=0; ieff<4; ieff++){
141       if(!(g=GetGraph(ieff, ""))) break;
142       if(!ieff){
143         if((h=(TH1S*)gROOT->FindObject("frame"))) delete h;
144         h=new TH1S("h","",100, 0., 15.);
145         h->SetLineColor(1);h->SetLineWidth(1);
146         h->SetMinimum(0.2);h->SetMaximum(1.2);
147         h->SetXTitle("p_{T} [GeV/c]");
148         h->SetYTitle("Efficiency");
149         h->Draw();
150       }
151       g->Draw("pc"); leg->AddEntry(g, g->GetTitle(), "pl");
152     }
153     leg->Draw();
154     return kTRUE;
155   case kTRDmom:
156     leg=new TLegend(.6, .75, .98, .98);
157     leg->SetBorderSize(1); leg->SetFillColor(0);
158     leg->SetHeader("Energy loss");
159     for(Int_t ieff=0; ieff<2; ieff++){
160       if(!(g=GetGraph(4+ieff, ""))) break;
161       if(!ieff){
162         if((h=(TH1S*)gROOT->FindObject("frame"))) delete h;
163         h=new TH1S("h","",100, -0.5, 5.5);
164         h->SetLineColor(0);h->SetLineWidth(1);
165         h->SetMinimum(-0.5);h->SetMaximum(1.2);
166         h->SetXTitle("layer");
167         h->SetYTitle("p_{T}^{Inner TPC} - p_{T}^{layer} [GeV/c]");
168         h->Draw();
169       }
170       g->Draw("p"); leg->AddEntry(g, g->GetTitle(), "pl");
171     }
172     leg->Draw();
173     return kTRUE;
174   default:
175     break;
176   }
177   AliInfo(Form("Reference plot [%d] missing result", ifig));
178   return kFALSE;
179 }
180
181 //____________________________________________________________________
182 TGraph* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
183 {
184 // Retrieve graph with "id"
185 // Possible options are :
186 //   "b" - build graph if none found
187 //   "c" - clear existing graph
188
189   Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
190          kCLEAR = strstr(opt, "c"); // clear existing graph
191
192   const Char_t *name[] = {
193     "Geo", "Trk", "Pid", "Ref", "Max06", "Mean09"
194   };
195   const Char_t *title[] = {
196     "TRD geometrical efficiency (TRDin/TPCout)"
197     ,"TRD tracking efficiency (TRDout/TRDin)"
198     ,"TRD PID efficiency (TRDpid/TRDin)"
199     ,"TRD refit efficiency (TRDrefit/TRDin)"
200     ,"TRD Eloss (Max/90% quantile)"
201     ,"TRD Eloss (Mean/60% quantile)"
202   };
203   const Int_t ngr = sizeof(name)/sizeof(Char_t*);
204   if(ngr != kNgraphs){
205     AliWarning("No of graphs defined different from definition");
206     return 0x0;
207   }
208
209   if(!fResults){
210     fResults = new TObjArray(kNgraphs);
211     fResults->SetOwner();
212     fResults->SetName("results");
213   }
214
215   TGraph *g = 0x0;
216   if((g = dynamic_cast<TGraph*>(fResults->At(id)))){
217     if(kCLEAR){ 
218       for(Int_t ip=g->GetN(); ip--;) g->RemovePoint(ip);
219     } else {
220       PutTrendValue(name[id], g->GetMean(2));
221       PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
222     }
223   } else {
224     if(kBUILD){
225       switch(id){
226       case 0:
227         g = new TGraphErrors();
228         g->SetMarkerStyle(7);g->SetMarkerColor(kBlack);
229         g->SetLineColor(kBlack);
230         break;
231       case 1:
232         g = new TGraphErrors();
233         g->SetMarkerStyle(7);g->SetMarkerColor(kRed);
234         g->SetLineColor(kRed);
235         break;
236       case 2:
237         g = new TGraphErrors();
238         g->SetMarkerStyle(7);g->SetMarkerColor(kBlue);
239         g->SetLineColor(kBlue);
240         break;
241       case 3:
242         g = new TGraphErrors();
243         g->SetMarkerStyle(7);g->SetMarkerColor(kGreen);
244         g->SetLineColor(kGreen);
245         break;
246       case 4:
247         g = new TGraphAsymmErrors(6);
248         g->SetMarkerStyle(22);g->SetMarkerColor(kRed);
249         g->SetLineColor(kBlack);g->SetLineWidth(2);
250         break;
251       case 5:
252         g = new TGraphAsymmErrors(6);
253         g->SetMarkerStyle(21);
254         g->SetLineColor(kRed);g->SetLineWidth(2);
255         break;
256       default:
257         AliWarning(Form("Graph index[%d] missing/not defined.", id));
258         return 0x0;
259       }
260       g->SetNameTitle(name[id], title[id]);
261       fResults->AddAt(g, id);
262     }
263   }
264   return g;
265 }
266
267 //____________________________________________________________________
268 void AliTRDcheckESD::Exec(Option_t *){
269   //
270   // Run the Analysis
271   //
272   if(!fESD){
273     AliError("ESD event missing.");
274     return;
275   }
276
277   // Get MC information if available
278   AliStack * fStack = 0x0;
279   if(HasMC()){
280     if(!fMC){ 
281       AliWarning("MC event missing");
282       SetMC(kFALSE);
283     } else {
284       if(!(fStack = fMC->Stack())){
285         AliWarning("MC stack missing");
286         SetMC(kFALSE);
287       }
288     }
289   }
290   Bool_t bTRDin(0), bTRDout(0), bTRDpid(0);
291
292   AliESDtrack *esdTrack = 0x0;
293   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
294     bTRDin=0;bTRDout=0;bTRDpid=0;
295     esdTrack = fESD->GetTrack(itrk);
296
297 //     if(esdTrack->GetNcls(1)) nTPC++;
298 //     if(esdTrack->GetNcls(2)) nTRD++;
299
300     // track status
301     ULong_t status = esdTrack->GetStatus();
302     //PrintStatus(status);
303
304     // define TPC out tracks
305     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
306     if(esdTrack->GetKinkIndex(0) > 0) continue;
307
308     // TRD PID
309     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
310     // pid quality
311     //esdTrack->GetTRDntrackletsPID();
312
313     // look at external track param
314     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
315     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
316     Double_t xyz[3];
317     if(op){
318       op->GetXYZ(xyz);
319       op->Global2LocalPosition(xyz, op->GetAlpha());
320       //printf("op @ X[%7.3f]\n", xyz[0]);
321     }
322
323     // read MC info
324     if(!HasMC()) continue;
325
326     Int_t fLabel = esdTrack->GetLabel();
327     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue; 
328     
329     // read MC particle
330     AliMCParticle *mcParticle = 0x0; 
331     if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(TMath::Abs(fLabel)))){
332       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
333       continue;
334     }
335
336     AliTrackReference *ref = 0x0; 
337     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
338     if(!nRefs){
339       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
340       continue;
341     }
342     Int_t iref = 0;
343     while(iref<nRefs){
344       ref = mcParticle->GetTrackReference(iref);
345       if(ref->LocalX() > fgkxTPC) break;
346       ref=0x0; iref++;
347     }
348
349     // read TParticle
350     //TParticle *tParticle = mcParticle->Particle(); 
351     //Int_t fPdg = tParticle->GetPdgCode();
352     // reject secondaries
353     //if(!tParticle->IsPrimary()) continue;
354
355     if(ref){ 
356       if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
357         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
358       } else {
359         bTRDin=1;
360         if(esdTrack->GetNcls(2)) bTRDout=1;
361         if(esdTrack->GetTRDntrackletsPID()>=4) bTRDpid=1;
362       }
363     } else { // track stopped in TPC 
364       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
365     }
366     // get the MC pt !!
367     Float_t pt = ref->Pt();
368
369     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
370     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
371     if(/*status & AliESDtrack::k*/bTRDin) h->Fill(pt, kTRDin);
372     if(/*status & AliESDtrack::k*/bTRDout){ 
373       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
374       h->Fill(pt, kTRDout);
375     }
376     if(/*status & AliESDtrack::k*/bTRDpid) h->Fill(pt, kTRDpid);
377     if(status & AliESDtrack::kTRDrefit) h->Fill(pt, kTRDref);
378
379     if(ip){
380       h = (TH2I*)fHistos->At(kTRDmom);
381       Float_t pp(0.);
382       for(Int_t ily=6; ily--;){
383         if((pp=esdTrack->GetTRDmomentum(ily))<0.) continue;
384         h->Fill(ip->GetP()-pp, ily);
385       }
386     }
387   }  
388   PostData(0, fHistos);
389 }
390
391 //____________________________________________________________________
392 TObjArray* AliTRDcheckESD::Histos()
393 {
394 // Retrieve histograms array if already build or build it
395
396   if(fHistos) return fHistos;
397
398   fHistos = new TObjArray(kNhistos);
399   //fHistos->SetOwner(kTRUE);
400   
401   TH1 *h = 0x0;
402
403   // clusters per tracklet
404   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
405     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
406     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
407     h->GetYaxis()->SetTitle("entries");
408   } else h->Reset();
409   fHistos->AddAt(h, kNCl);
410
411   // status bits histogram
412   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
413     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
414     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
415     h->GetYaxis()->SetTitle("status bits");
416     h->GetZaxis()->SetTitle("entries");
417   } else h->Reset();
418   fHistos->AddAt(h, kTRDstat);
419
420   // energy loss
421   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
422     h = new TH2I("hTRDmom", "TRD energy loss", 100, -1., 2., 6, -0.5, 5.5);
423     h->GetXaxis()->SetTitle("p_{inner} - p_{ly} [GeV/c]");
424     h->GetYaxis()->SetTitle("layer");
425     h->GetZaxis()->SetTitle("entries");
426   } else h->Reset();
427   fHistos->AddAt(h, kTRDmom);
428
429   return fHistos;
430 }
431
432 //____________________________________________________________________
433 Bool_t AliTRDcheckESD::Load(const Char_t *filename, const Char_t *name)
434 {
435 // Load data from performance file
436
437   if(!TFile::Open(filename)){
438     AliWarning(Form("Couldn't open file %s.", filename));
439     return kFALSE;
440   }
441   TObjArray *o = 0x0;
442   if(!(o = (TObjArray*)gFile->Get(name ? name : GetName()))){
443     AliWarning("Missing histogram container.");
444     return kFALSE;
445   }
446   fHistos = (TObjArray*)o->Clone(GetName());
447   gFile->Close();
448   return kTRUE;
449 }
450
451 //_______________________________________________________
452 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
453 {
454 // Dump trending value to default file
455
456   if(!fgFile){
457     fgFile = fopen("TRD.Performance.txt", "at");
458   }
459   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
460   return kTRUE;
461 }
462
463 //____________________________________________________________________
464 void AliTRDcheckESD::Terminate(Option_t *)
465 {
466 // Steer post-processing 
467
468
469   // geometrical efficiency
470   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
471   TH1 *h1[2] = {0x0, 0x0};
472   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
473   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
474   Process(h1, (TGraphErrors*)GetGraph(0));
475   delete h1[0];delete h1[1];
476
477   // tracking efficiency
478   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
479   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
480   Process(h1, (TGraphErrors*)GetGraph(1));
481   delete h1[1];
482
483   // PID efficiency
484   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
485   Process(h1, (TGraphErrors*)GetGraph(2));
486   delete h1[1];
487
488   // Refit efficiency
489   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
490   Process(h1, (TGraphErrors*)GetGraph(3));
491   delete h1[1];
492   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
493  
494   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)GetGraph(4), *g09 = (TGraphAsymmErrors*)GetGraph(5);
495   TAxis *ax=h2->GetXaxis();
496   const Int_t nq(4);
497   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
498   Double_t yq[nq];
499   for(Int_t ily=6; ily--;){
500     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
501     h1[0]->GetQuantiles(nq,yq,xq);
502     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
503     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
504     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
505     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
506
507     //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
508     delete h1[0];
509   }
510 }
511
512 //____________________________________________________________________
513 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
514 {
515 // Generic function to process one reference plot
516
517   Int_t n1 = 0, n2 = 0, ip=0;
518   Double_t eff = 0.;
519
520   TAxis *ax = h1[0]->GetXaxis();
521   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
522     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
523     n2 = (Int_t)h1[1]->GetBinContent(ib);
524     eff = n2/Float_t(n1);
525
526     ip=g->GetN();
527     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
528     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
529   }
530 }  
531
532 //____________________________________________________________________
533 void AliTRDcheckESD::PrintStatus(ULong_t status)
534 {
535 // Dump track status to stdout
536
537   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"
538     ,Bool_t(status & AliESDtrack::kITSin)
539     ,Bool_t(status & AliESDtrack::kITSout)
540     ,Bool_t(status & AliESDtrack::kITSrefit)
541     ,Bool_t(status & AliESDtrack::kTPCin)
542     ,Bool_t(status & AliESDtrack::kTPCout)
543     ,Bool_t(status & AliESDtrack::kTPCrefit)
544     ,Bool_t(status & AliESDtrack::kTPCpid)
545     ,Bool_t(status & AliESDtrack::kTRDin)
546     ,Bool_t(status & AliESDtrack::kTRDout)
547     ,Bool_t(status & AliESDtrack::kTRDrefit)
548     ,Bool_t(status & AliESDtrack::kTRDpid)
549     ,Bool_t(status & AliESDtrack::kTRDStop)
550     ,Bool_t(status & AliESDtrack::kHMPIDout)
551     ,Bool_t(status & AliESDtrack::kHMPIDpid)
552   );
553 }
554