last modification of the task before going into production
[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 #include "AliTRDgeometry.h"
42
43 #include "AliTRDcheckESD.h"
44
45 ClassImp(AliTRDcheckESD)
46
47 const Float_t AliTRDcheckESD::xTPC = 290.;
48 const Float_t AliTRDcheckESD::xTOF = 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
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   fHistos = new TObjArray(5);
100   //fHistos->SetOwner(kTRUE);
101   
102   TH1 *h = 0x0;
103
104   // clusters per tracklet
105   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
106     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
107     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
108     h->GetYaxis()->SetTitle("entries");
109   } else h->Reset();
110   fHistos->AddAt(h, kNCl);
111
112   // status bits histogram
113   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
114     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., kNbits, .5, kNbits+.5);
115     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
116     h->GetYaxis()->SetTitle("status bits");
117     h->GetZaxis()->SetTitle("entries");
118   } else h->Reset();
119   fHistos->AddAt(h, kTRDstat);
120
121   // results array
122   TObjArray *res = new TObjArray();
123   res->SetName("Results");
124   fHistos->AddAt(res, kResults);
125 }
126
127 //____________________________________________________________________
128 void AliTRDcheckESD::Exec(Option_t *){
129   //
130   // Run the Analysis
131   //
132   if(!fESD){
133     AliError("ESD event missing.");
134     return;
135   }
136
137   // Get MC information if available
138   AliStack * fStack = 0x0;
139   if(HasMC() && !fMC){ 
140     AliWarning("MC event missing");
141     SetMC(kFALSE);
142   } else {
143     if(!(fStack = fMC->Stack())){
144       AliWarning("MC stack missing");
145       SetMC(kFALSE);
146     }
147   }
148   Bool_t TRDin(0), TRDout(0), TRDpid(0);
149
150   AliESDtrack *esdTrack = 0x0;
151   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
152     TRDin=0;TRDout=0;TRDpid=0;
153     esdTrack = fESD->GetTrack(itrk);
154
155 //     if(esdTrack->GetNcls(1)) nTPC++;
156 //     if(esdTrack->GetNcls(2)) nTRD++;
157
158     // track status
159     ULong_t status = esdTrack->GetStatus();
160
161     // define TPC out tracks
162     if(!Bool_t(status & AliESDtrack::kTPCout)) continue;
163     if(esdTrack->GetKinkIndex(0) > 0) continue;
164
165     // TRD PID
166     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
167     // pid quality
168     esdTrack->GetTRDpidQuality();
169
170     // look at external track param
171     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
172     Double_t xyz[3];
173     if(op){
174       op->GetXYZ(xyz);
175       op->Global2LocalPosition(xyz, op->GetAlpha());
176       //printf("op @ X[%7.3f]\n", xyz[0]);
177     }
178
179     // read MC info
180     if(!HasMC()) continue;
181
182     Int_t fLabel = esdTrack->GetLabel();
183     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue; 
184     
185     // read MC particle
186     AliMCParticle *mcParticle = 0x0; 
187     if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
188       AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
189       continue;
190     }
191
192     AliTrackReference *ref = 0x0; 
193     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
194     if(!nRefs){
195       AliWarning(Form("Track refs missing. Label[%d].", fLabel));
196       continue;
197     }
198     Int_t iref = 0;
199     while(iref<nRefs){
200       ref = mcParticle->GetTrackReference(iref);
201       if(ref->LocalX() > xTPC) break;
202       ref=0x0; iref++;
203     }
204
205     // read TParticle
206     TParticle *tParticle = mcParticle->Particle(); 
207     //Int_t fPdg = tParticle->GetPdgCode();
208     // reject secondaries
209     //if(!tParticle->IsPrimary()) continue;
210
211     if(ref){ 
212       if(ref->LocalX() > xTOF){ // track skipping TRD fiducial volume
213         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
214       } else {
215         TRDin=1;
216         if(esdTrack->GetNcls(2)) TRDout=1;
217         if(esdTrack->GetTRDpidQuality()) TRDpid=1;
218       }
219     } else { // track stopped in TPC 
220       ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
221     }
222     // get the MC pt !!
223     Float_t pt = ref->Pt();
224
225     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
226     if(status & AliESDtrack::kTPCout) h->Fill(pt, kTPCout);
227     if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, kTRDin);
228     if(/*status & AliESDtrack::k*/TRDout){ 
229       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
230       h->Fill(pt, kTRDout);
231     }
232     if(status & AliESDtrack::kTRDpid) h->Fill(pt, kTRDpid);
233   }  
234   PostData(0, fHistos);
235 }
236
237
238 //____________________________________________________________________
239 void AliTRDcheckESD::Terminate(Option_t *)
240 {
241   TH2I *h2 = (TH2I*)fHistos->At(kTRDstat);
242   TAxis *ax = h2->GetXaxis();
243   TObjArray *res = (TObjArray*)fHistos->At(kResults);
244   
245   TH1 *h1[2] = {0x0, 0x0};
246   TGraphErrors *g = 0x0;
247   res->Expand(3);
248   Int_t n1 = 0, n2 = 0, ip=0;
249   Double_t eff = 0.;
250   
251   // geometrical efficiency
252   h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
253   h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
254   res->Add(g = new TGraphErrors());
255   g->SetNameTitle("geom", "TRD geometrical efficiency (TRDin/TPCout)");
256   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
257     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
258     n2 = (Int_t)h1[1]->GetBinContent(ib);
259     eff = n2/Float_t(n1);
260
261     ip=g->GetN();
262     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
263     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
264   }
265
266   // tracking efficiency
267   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
268   h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
269   res->Add(g = new TGraphErrors());
270   g->SetNameTitle("tracking", "TRD tracking efficiency (TRDout/TRDin)");
271   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
272     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
273     n2 = (Int_t)h1[1]->GetBinContent(ib);
274     eff = n2/Float_t(n1);
275
276     ip=g->GetN();
277     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
278     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
279   }
280
281   // PID efficiency
282   h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
283   h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
284   res->Add(g = new TGraphErrors());
285   g->SetNameTitle("PID", "TRD PID efficiency (TRDpid/TRDin)");
286   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
287     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
288     n2 = (Int_t)h1[1]->GetBinContent(ib);
289     eff = n2/Float_t(n1);
290
291     ip=g->GetN();
292     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
293     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
294   }
295 }
296
297 // if(TRDin && !TRDout){
298 //   printf("[%c] x[%2d]=%7.2f pt[%7.3f] id[%d] kink[%d] label[%d] in[%d] out[%d] pdg[%d]\n", tParticle->IsPrimary() ? 'P' : 'S', iref, ref->LocalX(), pt, esdTrack->GetID(), esdTrack->GetKinkIndex(0), fLabel, TRDin, TRDout, tParticle->GetPdgCode());
299 //   printf("\tstatus ITS[%d %d] TPC[%d %d %d] \n"
300 //     ,Bool_t(status & AliESDtrack::kITSin)
301 //     ,Bool_t(status & AliESDtrack::kITSout)
302 //     ,Bool_t(status & AliESDtrack::kTPCin)
303 //     ,Bool_t(status & AliESDtrack::kTPCout)
304 //     ,Bool_t(status & AliESDtrack::kTPCrefit)
305 //   );
306 //   printf("clusters ITS[%d] TPC[%d] TRD[%d]\n", esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2));
307 //   if(esdTrack->GetKinkIndex(0) != 0){
308 //     AliESDkink *kink = fESD->GetKink(TMath::Abs(esdTrack->GetKinkIndex(0)));
309 //     if(kink) printf("index[%d %d] label[%d %d]\n", kink->GetIndex(0), kink->GetIndex(1), kink->GetLabel(0), kink->GetLabel(1));
310 //   }
311 // }