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