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