changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / test / vmctest / scripts / efficiency / AliAnalysisTaskEfficiency.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TProfile.h"
6
7 #include "TCanvas.h"
8 #include "TList.h"
9 #include "TParticle.h"
10 #include "TParticlePDG.h"
11 #include "TProfile.h"
12
13 #include "AliAnalysisTask.h"
14 #include "AliAnalysisManager.h"
15
16 #include "AliESDEvent.h"
17 #include "AliStack.h"
18 #include "AliESDVertex.h"
19 #include "AliESDInputHandler.h"
20 #include "AliESDtrackCuts.h"
21 #include "AliMultiplicity.h"
22
23 #include "AliMCParticle.h"
24 #include "AliMCEvent.h"
25 #include "AliAnalysisTaskEfficiency.h"
26 #include "AliExternalTrackParam.h"
27 #include "AliTrackReference.h"
28 #include "AliStack.h"
29 #include "AliHeader.h"
30 #include "AliGenEventHeader.h"
31 #include "AliGenCocktailEventHeader.h"
32 #include "AliGenPythiaEventHeader.h"
33 #include "AliGenDPMjetEventHeader.h"
34
35 // Analysis Task for basic QA on the ESD
36 // Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing, Veronica Canoa, AM, Eva Sicking
37
38 ClassImp(AliAnalysisTaskEfficiency)
39
40 //________________________________________________________________________
41   AliAnalysisTaskEfficiency::AliAnalysisTaskEfficiency(const char *name) 
42     : AliAnalysisTaskSE(name) 
43     ,fFieldOn(kTRUE)
44     ,fHists(0)
45     ,fHistRECpt(0)
46     ,fHistMCpt(0)
47     ,fHistMCNRpt(0)
48     ,fHistFAKEpt(0)
49     ,fHistMULTpt(0)
50     ,fHistRECeta(0)
51     ,fHistMCeta(0)
52     ,fHistMCNReta(0)
53     ,fHistFAKEeta(0)
54     ,fHistMULTeta(0)
55     ,fHistRECphi(0)
56     ,fHistMCphi(0)
57     ,fHistMCNRphi(0)
58     ,fHistFAKEphi(0)
59     ,fHistMULTphi(0)
60     ,fHistRECHPTeta(0)
61     ,fHistMCHPTeta(0)
62     ,fHistMCNRHPTeta(0)
63     ,fHistFAKEHPTeta(0)
64     ,fHistRECHPTphi(0)
65     ,fHistMCHPTphi(0)
66     ,fHistMCNRHPTphi(0)
67     ,fHistFAKEHPTphi(0)
68     ,fHistRecMult(0)
69     ,fHistNCluster(0)  
70     ,fh1VtxEff(0) 
71     ,fh2VtxTpcSpd(0)
72     ,fh2EtaCorrelation(0)      
73     ,fh2EtaCorrelationShift(0)  
74     ,fh2PhiCorrelation(0)      
75     ,fh2PhiCorrelationShift(0)  
76     ,fh2PtCorrelation(0)      
77     ,fh2PtCorrelationShift(0)
78     ,fHistDeltaphiprimaries(0)
79     ,fHistDeltaphisecondaries(0)           
80     ,fHistDeltaphireject(0) 
81     ,fHistTPCITSdeltax(0)
82     ,fHistTPCITSdeltay(0)
83     ,fHistTPCITSdeltaz(0)
84     ,fHistTPCITSdeltar(0)
85     ,fHistTPCTRDdeltax(0)
86     ,fHistTPCTRDdeltay(0)
87     ,fHistTPCTRDdeltaz(0)
88     ,fHistTPCTRDdeltar(0)
89     ,fHistTPCTRDdeltaphi(0)
90     ,fPtPullsPos(0)
91     ,fPtPullsNeg(0)
92     ,fPtShiftPos(0)
93     ,fPtShiftNeg(0)
94     ,fTrackType(0)
95     ,fCuts(0)
96     ,fVertexRvsZ(0)
97     ,fVertexRvsZC(0)
98     ,fEtaMultiFluc(0)
99     ,fEtaMulti(0)
100     ,fEtaMultiH(0)
101     ,fPhiMultiFluc(0)
102     ,fPhiMulti(0)
103     ,fPhiMultiH(0)
104 {
105   //
106   // Constructor
107   //
108
109   for(int i = 0;i< 2;++i){
110     fh2MultSpdChips[i] = 0;
111   }
112
113   for(int i = 0;i < kTotalAC;++i){
114     fh2PhiPadRow[i] = fh2PhiLayer[i] = 0;
115   }
116
117   for(int i = 0; i < kTotalVtx; ++i){ 
118     fh2VertexCorrelation[i] = 
119       fh2VertexCorrelationShift[i] = 0;
120     fh1VertexShift[i] =         
121       fh1VertexShiftNorm[i] = 0;      
122   }
123
124   for(int i = 0;i< 8;++i){
125     fHistRECptCharge[i]=0;
126     fHistMCptCharge[i]=0;
127     fHistMCNRptCharge[i]=0;
128     fHistFAKEptCharge[i]=0;
129
130     fHistRECetaCharge[i]=0;
131     fHistMCetaCharge[i]=0;
132     fHistMCNRetaCharge[i]=0;
133     fHistFAKEetaCharge[i]=0;
134
135     fHistRECphiCharge[i]=0;
136     fHistMCphiCharge[i]=0;
137     fHistMCNRphiCharge[i]=0;
138     fHistFAKEphiCharge[i]=0;
139
140   }
141
142   DefineOutput(1,  TList::Class());
143 }
144
145
146 //________________________________________________________________________
147 void AliAnalysisTaskEfficiency::UserCreateOutputObjects()
148 {
149   // Create histograms
150
151
152   fHists = new TList();
153   fHistRECpt   = new TH1F("fHistRECpt",  " p_{T} distribution: all reconstructed",        100, 0., 20.);
154   fHistMCpt    = new TH1F("fHistMCpt",   " p_{T} distribution: all MC",                   100, 0., 20.);
155   fHistMCNRpt  = new TH1F("fHistMCNRpt", " p_{T} distribution: all not-reconstructed MC", 100, 0., 20.);
156   fHistFAKEpt  = new TH1F("fHistFAKEpt", " p_{T} distribution: all Fake",                 100, 0., 20.);
157   fHistMULTpt  = new TH1F("fHistMULTpt", " p_{T} distribution: multiply rec.",            100, 0., 20.);
158     
159   fHistRECeta   = new TH1F("fHistRECeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);
160   fHistMCeta    = new TH1F("fHistMCeta",   " #eta distribution: all MC",                   100,-1.0,1.0);
161   fHistMCNReta  = new TH1F("fHistMCNReta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);
162   fHistFAKEeta  = new TH1F("fHistFAKEeta", " #eta distribution: all Fake",                 100,-1.0,1.0);
163   fHistMULTeta  = new TH1F("fHistMULTeta", " #eta distribution: multiply rec.",            100,-1.0,1.0);
164
165   fHistRECphi   = new TH1F("fHistRECphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);
166   fHistMCphi    = new TH1F("fHistMCphi",   " #phi distribution: all MC",                   314, 0., 6.28);
167   fHistMCNRphi  = new TH1F("fHistMCNRphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);
168   fHistFAKEphi  = new TH1F("fHistFAKEphi", " #phi distribution: all Fake",                 314, 0., 6.28);
169   fHistMULTphi  = new TH1F("fHistMULTphi", " #phi distribution: multipli rec.",            314, 0., 6.28);
170
171   fHistRECHPTeta   = new TH1F("fHistRECHPTeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);
172   fHistMCHPTeta    = new TH1F("fHistMCHPTeta",   " #eta distribution: all MC",                   100,-1.0,1.0);
173   fHistMCNRHPTeta  = new TH1F("fHistMCNRHPTeta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);
174   fHistFAKEHPTeta  = new TH1F("fHistFAKEHPTeta", " #eta distribution: all Fake",                 100,-1.0,1.0);
175
176   fHistRECHPTphi   = new TH1F("fHistRECHPTphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);
177   fHistMCHPTphi    = new TH1F("fHistMCHPTphi",   " #phi distribution: all MC",                   314, 0., 6.28);
178   fHistMCNRHPTphi  = new TH1F("fHistMCNRHPTphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);
179   fHistFAKEHPTphi  = new TH1F("fHistFAKEHPTphi", " #phi distribution: all Fake",                 314, 0., 6.28);
180
181   fHistRecMult     = new TH1F("fHistRecMult",  "Multiple reconstructed tracks", 50, 0., 50.);
182   fHistNCluster    = new TH1F("fHistNCluster", "Number of clusters for suspicious tracks", 300, 0., 300.);
183
184   // CKB
185
186   fh1VtxEff = new TH1F("fh1VtxEff","VtxEff TPC",4,-0.5,3.5);
187   fh1VtxEff->GetXaxis()->SetBinLabel(1,"NO TPC VTX");
188   fh1VtxEff->GetXaxis()->SetBinLabel(2,"TPC VTX");
189   fh1VtxEff->GetXaxis()->SetBinLabel(3,"NO SPD VTX");
190   fh1VtxEff->GetXaxis()->SetBinLabel(4,"SPD VTX");
191
192   TString labels[kTotalAC];
193   labels[kNegA]="NegA";
194   labels[kPosA]="PosA";
195   labels[kNegC]="NegC";
196   labels[kPosC]="PosC";
197
198   for(int i = 0;i< kTotalAC;++i){
199     fh2PhiPadRow[i] = new TH2F(Form("fh2PhiPadRow_%d",i),Form("Padrow vs phi AC %s;phi;padrow",labels[i].Data()),360,0.,360.,159,-0.5,158.5);
200     fh2PhiLayer[i] = new TH2F(Form("fh2PhiLayer_%d",i),Form("layer vs phi AC %s;phi;layer",labels[i].Data()),360,0.,360.,6,-0.5,5.5);
201   }
202   for(int i = 0;i<2;++i){
203     fh2MultSpdChips[i] = new TH2F(Form("fh2ChipsSpdMult_%d",i),"mult SPD vs chips;chips;Mult",201,-1.5,199.5,201,-1.5,199.5);
204   }
205   fh2VtxTpcSpd  = new TH2F("fh2VtxTpcSpd","SPD vtx vs TPC;tpc-z;spd-z",200,-20.,20.,200,-20,20.);
206
207   TString labelsVtx[kTotalVtx];
208   labelsVtx[kTPC] = "TPC";
209   labelsVtx[kSPD] = "SPD";
210   for(int i = 0; i < kTotalVtx; ++i){ 
211     fh2VertexCorrelation[i] = new TH2F(Form("fh2VertexCorrelation_%d",i),Form("VertexCorrelation %s;MC z-vtx;ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 120, -30, 30);
212     fh2VertexCorrelationShift[i] = new TH2F(Form("fh2VertexCorrelationShift_%d",i), Form("VertexCorrelationShift %s;MC z-vtx;MC z-vtx - ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 100, -1, 1); 
213     fh1VertexShift[i] =  new TH1F(Form("fh1VertexShift_%d",i), Form("VertexShift %s;(MC z-vtx - ESD z-vtx);Entries",labelsVtx[i].Data()), 201, -2, 2);       
214     fh1VertexShiftNorm[i] = new TH1F(Form("fh1VertexShiftNorm_%d",i), Form("VertexShiftNorm %s;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries",labelsVtx[i].Data()), 200, -100, 100);      
215   }
216
217     
218   fh2EtaCorrelation = new TH2F("fh2EtaCorrelation", "EtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
219   fh2EtaCorrelationShift = new TH2F("fh2EtaCorrelationShift", "EtaCorrelationShift;ESD #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
220   fh2PhiCorrelation = new TH2F("fh2PhiCorrelation", "PhiCorrelation;MC #phi;ESD #phi", 120, 0, 360, 120, 0, 360);
221   fh2PhiCorrelationShift = new TH2F("fh2PhiCorrelationShift", "PhiCorrelationShift;ESD #phi;MC #phi - ESD #phi", 120, 0, 360, 100, -10, 10);
222   fh2PtCorrelation = new TH2F("fh2PtCorrelation", "PtCorrelation;MC p_T;ESD p_T", 200, 0., 200., 200, 0, 200);
223   fh2PtCorrelationShift = new TH2F("fh2PtCorrelationShift", "PtCorrelationShift;ESD p_T;MC p_T - ESD #p_T", 120, 0, 10, 100, -2, 2);
224   
225   fHistDeltaphiprimaries  =new TH1F("fHistDeltaphiprimaries", " #Delta#phi distribution: primaries",     314, -0.2,0.2);               
226   fHistDeltaphisecondaries=new TH1F("fHistDeltaphisecondaries", "#Delta#phi distribution: secondaries", 314, -0.2,0.2);              
227   fHistDeltaphireject     =new TH1F("fHistDeltaphireject", " #Delta#phi distribution: reject trackled",  314, -0.2,0.2); 
228   
229   fHistTPCITSdeltax =new TH1F("fHistTPCITSdeltax", "TPC-ITS matching:#Delta x", 100,-20,20); 
230   fHistTPCITSdeltay =new TH1F("fHistTPCITSdeltay", "TPC-ITS matching:#Delta y", 100,-20,20); 
231   fHistTPCITSdeltaz =new TH1F("fHistTPCITSdeltaz", "TPC-ITS matching:#Delta z", 100,-20,20); 
232   fHistTPCITSdeltar =new TH1F("fHistTPCITSdeltar", "TPC-ITS matching:#Delta r", 100,0,20); 
233   
234   fHistTPCTRDdeltax   = new TH1F("fHistTPCTRDdeltax", "TPC-TRD matching:#Delta x", 400,-400, 400); 
235   fHistTPCTRDdeltay   = new TH1F("fHistTPCTRDdeltay", "TPC-TRD matching:#Delta y", 400,-400, 400); 
236   fHistTPCTRDdeltaz   = new TH1F("fHistTPCTRDdeltaz", "TPC-TRD matching:#Delta z", 400,-400, 400); 
237   fHistTPCTRDdeltar   = new TH1F("fHistTPCTRDdeltar", "TPC-TRD matching:#Delta r", 100,0,20); 
238   fHistTPCTRDdeltaphi = new TH1F("fHistTPCTRDdeltarphi", "TPC-TRD matching:#Delta #phi", 128,-3.14 , 3.14); 
239   // Pulls
240   fPtPullsPos = new TProfile("fPtPullsPos", "Pt Pulls for pos. primaries", 50 , 0., 10., -5., 5.,"S");
241   fPtPullsNeg = new TProfile("fPtPullsNeg", "Pt Pulls for neg. primaries", 50 , 0., 10., -5., 5.,"S");
242   fPtShiftPos = new TProfile("fPtShiftPos", "Pt Shift for pos. primaries", 50 , 0., 10., -0.2, 0.2,"S");
243   fPtShiftNeg = new TProfile("fPtShiftNeg", "Pt Shift for neg. primaries", 50 , 0., 10., -0.2, 0.2,"S");
244
245   fEtaMultiFluc = new TProfile("fEtaMultiFluc", "eta multiplicity fluctuations", 10, -2., 2., 0., 600., "S");
246   fEtaMulti     = new TH1F("fEtaMulti",     "eta multiplicity fluctuations", 10, -2., 2.);
247   fEtaMultiH    = new TH1F("fEtaMultiH",    "eta multiplicity fluctuations", 600,  0., 600.);
248
249   fPhiMultiFluc = new TProfile("fPhiMultiFluc", "phi multiplicity fluctuations", 64, 0., 6.4, 0., 100., "S");
250   fPhiMulti     = new TH1F("fPhiMulti",         "phi multiplicity fluctuations", 64, 0., 6.4);
251   fPhiMultiH    = new TH1F("fPhiMultiH",        "phi multiplicity fluctuations", 100, 0., 100.);
252
253   // SPD Vertex
254   fVertexRvsZ  = new TH2F("fVertexRvsZ", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);
255   fVertexRvsZC = new TH2F("fVertexRvsZC", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);
256   
257   TString charge[8]={"Deu", "AntiDeu", "Tri", "AntiTri", "He3", "AntiHe3", "He4", "AntiHe4"};
258   for(Int_t i=0;i<8;i++){
259     fHistRECptCharge[i]   = new TH1F(Form("fHistRECptCharge%s",charge[i].Data()), 
260                                      "p_{T} distribution: all reconstructed",
261                                      100, 0., 20.);
262     fHistMCptCharge[i]    = new TH1F(Form("fHistMCptCharge%s",charge[i].Data()), 
263                                      "p_{T} distribution: all MC",
264                                      100, 0., 20.);
265     fHistMCNRptCharge[i]  = new TH1F(Form("fHistMCNRptCharge%s",charge[i].Data()), 
266                                      "p_{T} distribution: all not-reconstructed MC", 
267                                      100, 0., 20.);
268     fHistFAKEptCharge[i]  = new TH1F( Form("fHistFAKEptCharge%s",charge[i].Data()), 
269                                      "p_{T} distribution: all Fake",                 
270                                      100, 0., 20.);
271
272     fHistRECetaCharge[i]   = new TH1F(Form("fHistRECetaCharge%s",charge[i].Data()), 
273                                      "p_{T} distribution: all reconstructed",
274                                      100, 0., 20.);
275     fHistMCetaCharge[i]    = new TH1F(Form("fHistMCetaCharge%s",charge[i].Data()), 
276                                      "p_{T} distribution: all MC",
277                                      100, 0., 20.);
278     fHistMCNRetaCharge[i]  = new TH1F(Form("fHistMCNRetaCharge%s",charge[i].Data()), 
279                                      "p_{T} distribution: all not-reconstructed MC", 
280                                      100, 0., 20.);
281     fHistFAKEetaCharge[i]  = new TH1F( Form("fHistFAKEetaCharge%s",charge[i].Data()), 
282                                      "p_{T} distribution: all Fake",                 
283                                      100, 0., 20.);
284
285     fHistRECphiCharge[i]   = new TH1F(Form("fHistRECphiCharge%s",charge[i].Data()), 
286                                      "p_{T} distribution: all reconstructed",
287                                      100, 0., 20.);
288     fHistMCphiCharge[i]    = new TH1F(Form("fHistMCphiCharge%s",charge[i].Data()), 
289                                      "p_{T} distribution: all MC",
290                                      100, 0., 20.);
291     fHistMCNRphiCharge[i]  = new TH1F(Form("fHistMCNRphiCharge%s",charge[i].Data()), 
292                                      "p_{T} distribution: all not-reconstructed MC", 
293                                      100, 0., 20.);
294     fHistFAKEphiCharge[i]  = new TH1F( Form("fHistFAKEphiCharge%s",charge[i].Data()), 
295                                      "p_{T} distribution: all Fake",                 
296                                      100, 0., 20.);
297
298   }
299   // BKC
300
301
302   fHists->SetOwner();
303
304   fHists->Add(fHistRECpt);
305   fHists->Add(fHistMCpt);
306   fHists->Add(fHistMCNRpt);
307   fHists->Add(fHistFAKEpt);
308   fHists->Add(fHistMULTpt);
309
310   fHists->Add(fHistRECeta);
311   fHists->Add(fHistMCeta);
312   fHists->Add(fHistMCNReta);
313   fHists->Add(fHistFAKEeta);
314   fHists->Add(fHistMULTeta);
315
316   fHists->Add(fHistRECphi);
317   fHists->Add(fHistMCphi);
318   fHists->Add(fHistMCNRphi);
319   fHists->Add(fHistFAKEphi);
320   fHists->Add(fHistMULTphi);
321
322   fHists->Add(fHistRECHPTeta);
323   fHists->Add(fHistMCHPTeta);
324   fHists->Add(fHistMCNRHPTeta);
325   fHists->Add(fHistFAKEHPTeta);
326
327   fHists->Add(fHistRECHPTphi);
328   fHists->Add(fHistMCHPTphi);
329   fHists->Add(fHistMCNRHPTphi);
330   fHists->Add(fHistFAKEHPTphi);
331
332   // CKB
333
334   fHists->Add(fh1VtxEff);
335   for(int i = 0;i < kTotalAC;++i){
336     fHists->Add(fh2PhiPadRow[i]);
337     fHists->Add(fh2PhiLayer[i]);
338   }
339   for(int i = 0;i<2;++i){
340     fHists->Add(fh2MultSpdChips[i]);
341   }
342   fHists->Add(fh2VtxTpcSpd);
343     
344   for(int i = 0; i < kTotalVtx; ++i){ 
345     fHists->Add(fh2VertexCorrelation[i]);
346     fHists->Add(fh2VertexCorrelationShift[i]);
347     fHists->Add(fh1VertexShift[i]);      
348     fHists->Add(fh1VertexShiftNorm[i]);
349   }
350
351
352   fHists->Add(fh2EtaCorrelation);
353   fHists->Add(fh2EtaCorrelationShift);
354   fHists->Add(fh2PhiCorrelation);
355   fHists->Add(fh2PhiCorrelationShift);
356   fHists->Add(fh2PtCorrelation);
357   fHists->Add(fh2PtCorrelationShift);
358
359   fHists->Add(fHistDeltaphiprimaries);              
360   fHists->Add(fHistDeltaphisecondaries);
361   fHists->Add(fHistDeltaphireject);  
362
363   fHists->Add(fHistTPCITSdeltax);  
364   fHists->Add(fHistTPCITSdeltay);  
365   fHists->Add(fHistTPCITSdeltaz);  
366   fHists->Add(fHistTPCITSdeltar);  
367     
368   fHists->Add(fHistTPCTRDdeltax);  
369   fHists->Add(fHistTPCTRDdeltay);  
370   fHists->Add(fHistTPCTRDdeltaz);  
371   fHists->Add(fHistTPCTRDdeltar);  
372   fHists->Add(fHistTPCTRDdeltaphi);  
373
374   fHists->Add(fHistRecMult);
375   fHists->Add(fHistNCluster);
376  
377
378   fHists->Add(fPtPullsPos);
379   fHists->Add(fPtPullsNeg);
380   fHists->Add(fPtShiftPos);
381   fHists->Add(fPtShiftNeg);
382
383   fHists->Add(fVertexRvsZ);
384   fHists->Add(fVertexRvsZC);
385   fHists->Add(fEtaMultiFluc);
386   fHists->Add(fEtaMulti);
387   fHists->Add(fEtaMultiH);
388   fHists->Add(fPhiMultiFluc);
389   fHists->Add(fPhiMulti);
390   fHists->Add(fPhiMultiH);
391
392   for(Int_t i=0;i<8;i++){
393     fHists->Add(fHistRECptCharge[i]);
394     fHists->Add(fHistMCptCharge[i]);
395     fHists->Add(fHistMCNRptCharge[i]);
396     fHists->Add(fHistFAKEptCharge[i]);
397
398     fHists->Add(fHistRECetaCharge[i]);
399     fHists->Add(fHistMCetaCharge[i]);
400     fHists->Add(fHistMCNRetaCharge[i]);
401     fHists->Add(fHistFAKEetaCharge[i]);
402
403     fHists->Add(fHistRECphiCharge[i]);
404     fHists->Add(fHistMCphiCharge[i]);
405     fHists->Add(fHistMCNRphiCharge[i]);
406     fHists->Add(fHistFAKEphiCharge[i]);
407   }
408
409   for (Int_t i=0; i<fHists->GetEntries(); ++i) {
410     TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));
411     if (h1){
412       h1->Sumw2();
413     }
414   }
415   // BKC
416
417
418   // Post output data.
419   PostData(1, fHists);
420
421 }
422
423 //________________________________________________________________________
424 void AliAnalysisTaskEfficiency::UserExec(Option_t *) 
425 {
426   // Event loop
427   // Analysis of MC true particles and reconstructed tracks
428   // Different track types (Global, TPC, ITS) can be selected
429
430   if (!fInputEvent) {
431     Printf("ERROR: fESD not available");
432     return;
433   }
434   static Int_t nfc = 0;
435
436   //AliESDInputHandler* esdI = (AliESDInputHandler*) 
437   //(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler();
438   //AliESDEvent* hltEvent = esdI->GetHLTEvent();
439
440   AliStack* stack = MCEvent()->Stack();
441   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
442     
443   // CKB
444   // Fetch the SPD vertex:
445   Double_t vtxSPD[3];
446   const AliESDVertex* spdVertex = esdE->GetPrimaryVertexSPD();
447   if (spdVertex) {
448     if(spdVertex->GetNContributors()<=0){
449       spdVertex = 0;
450     }
451     else{
452       spdVertex->GetXYZ(vtxSPD);
453     }
454   }
455   
456   // Fetch the TPC vertex
457   Double_t vtxTPC[3];
458   const AliESDVertex* tpcVertex = esdE->GetPrimaryVertexTPC();
459   if (tpcVertex) {
460     if(tpcVertex->GetNContributors()<=0){
461       tpcVertex = 0;
462     }
463     else{
464       tpcVertex->GetXYZ(vtxTPC);
465     }
466   }
467   // SPD Vertex
468   if (spdVertex) {
469     Double_t x,y,z;
470     x = spdVertex->GetX() + 0.07;
471     y = spdVertex->GetY() - 0.25;
472     z = spdVertex->GetZ();
473     if (TMath::Abs(z) < 10.) {
474       fVertexRvsZ->Fill(x,y);
475       if (TMath::Sqrt(x*x + y*y) > 5. * 0.028) 
476         fVertexRvsZC->Fill(x,y);
477     }
478   }
479   // BKC
480   //Printf("%s:%d %5d",(char*)__FILE__,__LINE__, esdE->GetNumberOfTracks());
481   TArrayI labels(esdE->GetNumberOfTracks());
482   Int_t igood = 0;
483   // Track loop to fill a pT spectrum
484   //printf("ESD track loop \n");
485
486   AliESDtrack *tpcP = 0x0;
487
488   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
489
490     //prevent mem leak for TPConly track
491     if(fTrackType==2&&tpcP){
492       delete tpcP;
493       tpcP = 0;
494     }
495
496     AliVParticle *track = esdE->GetTrack(iTracks);
497     AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
498     esdtrack->PropagateToDCA(esdE->GetPrimaryVertex(),
499                              esdE->GetMagneticField(), 10000.);
500
501     if (!track) {
502       Printf("ERROR: Could not receive track %d", iTracks);
503       continue;
504     }
505     //__________
506     // run Task for global tracks or ITS tracks or TPC tracks
507     const AliExternalTrackParam *tpcPin = 0x0;
508     //Double_t phiIn=0.;
509     Float_t phi = 0.;
510     Float_t eta = 0.;
511     Float_t pt  = 0.;
512     Int_t indexAC = GetIndexAC(esdtrack);  
513
514     //select in the steering macro, which track type should be analysed. 
515     // Four track types are selectable:
516     // 0. Global tracks
517     // 1. ITS tracks (SA or Pure SA)
518     // 2. TPC only tracks
519     // Track selection copied from PWGPP/AliAnalysisTaskQAsym
520
521     if(fTrackType==0){
522       //Fill all histograms with global tracks
523       tpcP = esdtrack;
524       if (!tpcP) continue;
525       if (!fCuts->AcceptTrack(tpcP)) continue;
526       phi = tpcP->Phi();
527       eta = tpcP->Eta();
528       pt  = tpcP->Pt();
529     }
530     else if(fTrackType==1){
531       //Fill all histograms with ITS tracks
532       tpcP = esdtrack;
533       if (!tpcP) continue;
534       if (!fCuts->AcceptTrack(tpcP)) continue;
535       phi = tpcP->Phi();
536       eta = tpcP->Eta();
537       pt  = tpcP->Pt();
538     }
539     else if(fTrackType==2){     
540       //Fill all histograms with TPC track information
541       tpcPin = esdtrack->GetInnerParam();
542       if (!tpcPin) continue;
543       tpcP = AliESDtrackCuts::
544         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdE),
545                         esdtrack->GetID());
546       if (!tpcP) continue;
547       if (!fCuts->AcceptTrack(tpcP)) continue;
548       if(tpcP->GetNcls(1)>160)continue;//jacek's track cut
549       if(tpcP->GetConstrainedChi2TPC()<0)continue; // jacek's track cut
550       phi=tpcPin->Phi();
551       eta=tpcPin->Eta();
552       pt=tpcPin->Pt();
553     }
554     else{
555       Printf("ERROR: wrong track type \n");
556       continue;
557     }
558     //___________
559     //
560
561
562     labels.AddAt(-1, iTracks);
563       
564
565       // Selected
566       if (TMath::Abs(eta) > 0.9) {
567       } else {
568
569         //Int_t charge=track->Charge();
570         // Within acceptance
571         Int_t ind = TMath::Abs(esdtrack->GetLabel());
572         AliMCParticle* mcPtest = (AliMCParticle*) MCEvent()->GetTrack(ind);
573         if (stack->IsPhysicalPrimary(mcPtest->Label())){
574           if(TMath::Abs(mcPtest->PdgCode())>=1000020030){//all helium (+-1000020030,+-1000020040)
575             pt*=2;// reconstruction takes charge=1 -> for particles with charge=2, pT,rec need to be corrected
576           }
577           
578           Int_t index=ConvertePDG(mcPtest->PdgCode());
579           if(fDebug>1)Printf("PDG=%d, index=%d", mcPtest->PdgCode(), index);
580           //fill tracks comming from d,t,3He, 4He
581           if(index<8){
582             fHistRECptCharge [index] ->Fill(pt);
583             fHistRECetaCharge[index]->Fill(eta);
584             fHistRECphiCharge[index]->Fill(phi);
585           }
586
587         }
588         fHistRECpt ->Fill(pt);
589         fHistRECeta->Fill(eta);
590         fHistRECphi->Fill(phi);
591
592         if (pt > 2.) {
593           fHistRECHPTeta->Fill(eta);
594           fHistRECHPTphi->Fill(phi);
595         }
596           
597         //      Int_t ind = TMath::Abs(esdtrack->GetLabel());
598         AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(ind);
599         if (!(stack->IsPhysicalPrimary(mcP->Label()))) {
600           fHistFAKEpt ->Fill(pt);
601           fHistFAKEeta->Fill(eta);                
602           fHistFAKEphi->Fill(phi);
603           
604           if (pt > 2.) {
605             fHistFAKEHPTeta->Fill(eta);
606             fHistFAKEHPTphi->Fill(phi);
607           }
608             
609         }
610           
611         labels.AddAt(TMath::Abs(esdtrack->GetLabel()), iTracks);
612         igood++;
613           
614           
615         Float_t phiDeg = TMath::RadToDeg()*track->Phi();
616           
617         //TPC-ITS matching
618         //TPC-TRD matching  
619               
620         if (tpcP){
621           Int_t labeltpcits = esdtrack->GetLabel();
622           AliExternalTrackParam * tpcPCopy = new AliExternalTrackParam(*tpcP);
623           Double_t xk = 43.6;  // cm
624           Double_t bz = 5.0;   // kG
625           if(tpcPCopy->PropagateTo(xk,bz)){
626             Double_t alpha=tpcPCopy->GetAlpha();
627             // if(tpcPCopy->Rotate(0.)){ 
628             Float_t xtpc=tpcPCopy->GetX();
629             Float_t ytpc=tpcPCopy->GetY();
630             Float_t ztpc=tpcPCopy->GetZ();
631             Float_t xpr,ypr ; 
632             xpr=xtpc*TMath::Cos(alpha)-ytpc*TMath::Sin(alpha);
633             ypr=xtpc*TMath::Sin(alpha)+ytpc*TMath::Cos(alpha);
634                 
635             AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpcits));
636             Int_t ntref = mcPart->GetNumberOfTrackReferences();
637                  
638             for (Int_t k = 0; k < ntref; k++) {
639               AliTrackReference* ref = mcPart->GetTrackReference(k);
640               Float_t xhits = ref->X();
641               Float_t yhits = ref->Y();
642               Float_t radio = TMath::Sqrt(xhits*xhits+yhits*yhits);
643                 
644               if(ref->DetectorId() == 0 && (radio > 42.8 && radio < 43.7)) {
645                   
646                 Float_t xits = ref->X();
647                 Float_t yits = ref->Y();
648                 Float_t zits = ref->Z();
649                   
650                 Float_t difx=(xits-xpr);
651                 fHistTPCITSdeltax->Fill(difx);
652                 Float_t dify=(yits-ypr);
653                 fHistTPCITSdeltay->Fill(dify);
654                 Float_t difz=(zits-ztpc);
655                 fHistTPCITSdeltaz->Fill(difz);
656                 Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);
657                 fHistTPCITSdeltar->Fill(deltar);
658                 break;
659               }
660               // }
661             } // trackrefs
662           } 
663         } //ITS-TPC maching
664           //TPC-TRD maching 
665                
666         const AliExternalTrackParam *trd = esdtrack->GetOuterParam();
667           
668         if (trd){Int_t labeltpctrd = track->GetLabel();
669           
670           AliExternalTrackParam * trdCopy = new AliExternalTrackParam(*trd);
671           Double_t xktrd=296.0;
672           Double_t bztrd=5.0;
673           if(trdCopy->PropagateTo(xktrd,bztrd)){
674             Float_t xtpc2=trdCopy->GetX();
675             Float_t ytpc2=trdCopy->GetY();
676             Float_t ztpc2=trdCopy->GetZ();
677             Double_t alpha=trdCopy->GetAlpha();
678             Float_t xpr,ypr ; 
679             xpr=xtpc2*TMath::Cos(alpha)-ytpc2*TMath::Sin(alpha);
680             ypr=xtpc2*TMath::Sin(alpha)+ytpc2*TMath::Cos(alpha);
681             AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpctrd));
682             Int_t ntreftrd = mcPart->GetNumberOfTrackReferences(); 
683             
684             for (Int_t k = 0; k < ntreftrd; k++) {
685               AliTrackReference* ref = mcPart->GetTrackReference(k);
686               if(ref->DetectorId() == 3 && ref->Label() == labeltpctrd){
687                 
688                 Float_t xtrd = ref->X();
689                 Float_t ytrd = ref->Y();
690                 Float_t ztrd = ref->Z();
691                 
692                 Float_t difx=(xtrd-xpr);
693                 fHistTPCTRDdeltax->Fill(difx);
694                 Float_t dify=(ytrd-ypr);
695                 fHistTPCTRDdeltay->Fill(dify);
696                 Float_t difz=(ztrd-ztpc2);
697                 fHistTPCTRDdeltaz->Fill(difz);
698                 Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);
699                 fHistTPCTRDdeltar->Fill(deltar);
700                 Float_t phi_tpc = TMath::ATan2(ypr, xpr);
701                 Float_t phi_trd = TMath::ATan2(ytrd, xtrd);
702                 Float_t dphi = phi_trd - phi_tpc;
703                 if (dphi >   TMath::Pi()) dphi -= 2. * TMath::Pi();
704                 if (dphi < - TMath::Pi()) dphi += 2. * TMath::Pi();
705
706                 fHistTPCTRDdeltaphi->Fill(dphi);
707                 break;
708               }
709               
710             }
711           }
712
713         }//TRD-TPC maching
714
715         // CKB
716         if(pt>2.&&fFieldOn){
717           TBits bits(esdtrack->GetTPCClusterMap());
718           for(unsigned int ip = 0;ip < bits.GetNbits();++ip){
719             if(bits[ip]){
720               fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);
721             }
722           }
723           for(int i = 0;i < 6;++i){
724             if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);
725           }
726         }
727         else if(!fFieldOn){ // field not on all momenta are set to MPV 
728           TBits bits(esdtrack->GetTPCClusterMap());
729           for(unsigned int ip = 0;ip < bits.GetNbits();++ip){
730             if(bits[ip]){
731               fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);
732             }
733             for(int i = 0;i < 6;++i){
734               if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);
735             }
736           }      
737
738         }
739         // Fill the correlation
740         if(MCEvent()){
741           TParticle *part = MCEvent()->Stack()->Particle(TMath::Abs(esdtrack->GetLabel()));
742           if(part){
743             Float_t mcEta = part->Eta();
744             Float_t mcPhi = TMath::RadToDeg()*part->Phi();
745             Float_t mcPt  = part->Pt();
746             //             if (pt - mcPt > 20.) {
747             fh2EtaCorrelation->Fill(mcEta,eta);
748             fh2EtaCorrelationShift->Fill(eta,mcEta-eta);
749             fh2PhiCorrelation->Fill(mcPhi,phiDeg);
750             fh2PhiCorrelationShift->Fill(phiDeg,mcPhi-phiDeg);
751             fh2PtCorrelation->Fill(mcPt,pt);
752             fh2PtCorrelationShift->Fill(pt,mcPt-pt);
753             //}
754             Double_t sigmaPt = TMath::Sqrt(esdtrack->GetSigma1Pt2());
755             if (track->Charge() > 0.) {
756               fPtPullsPos->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt); 
757               fPtShiftPos->Fill(mcPt, (1./pt - 1./mcPt));
758             } else {
759               fPtPullsNeg->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt); 
760               fPtShiftNeg->Fill(mcPt, (1./pt - 1./mcPt));  
761             }
762           }
763           // BKC
764         }
765       }
766    
767   } //track loop 
768
769   //prevent mem leak for TPConly track
770   if(fTrackType==2&&tpcP){
771     delete tpcP;
772     tpcP = 0;
773   }
774
775   //Printf("%s:%d",(char*)__FILE__,__LINE__);
776   // CKB
777   // Vertex eff
778   if(tpcVertex){
779     fh1VtxEff->Fill("TPC VTX",1);
780   }
781   else{
782     fh1VtxEff->Fill("NO TPC VTX",1);
783   }
784
785   if(spdVertex){
786     fh1VtxEff->Fill("SPD VTX",1);
787   }
788   else{
789     fh1VtxEff->Fill("NO SPD VTX",1);
790   }
791
792
793   // Vertex correlation SPD vs. TPC
794   if(tpcVertex&&spdVertex){
795     fh2VtxTpcSpd->Fill(vtxTPC[2],vtxSPD[2]);
796   }
797   //  Printf("%s:%d",(char*)__FILE__,__LINE__);
798   // Multiplicity checks in the SPD
799   const AliMultiplicity *spdMult = esdE->GetMultiplicity();
800   // Multiplicity Analysis
801   Int_t nt = spdMult->GetNumberOfTracklets();
802   nfc++;
803   if (nfc == 520) {
804     nfc = 0;
805     for (Int_t ib = 1; ib <= 10; ib++) {
806       Int_t ic = Int_t(fEtaMulti->GetBinContent(ib));
807       Float_t eta = -1.8 + Float_t(ib-1) * 0.4;
808       fEtaMultiFluc->Fill(eta, Float_t(ic));
809       if (ib == 5) fEtaMultiH->Fill(Float_t(ic));
810     }
811
812     for (Int_t ib = 1; ib <= 64; ib++) {
813       Int_t ic = Int_t(fPhiMulti->GetBinContent(ib));
814       Float_t phi = 0.05 + Float_t(ib-1) * 0.1;
815       fPhiMultiFluc->Fill(phi, Float_t(ic));
816       if (ib  == 2) fPhiMultiH->Fill(Float_t(ic));
817     }
818
819     fEtaMulti->Reset();
820     fPhiMulti->Reset();    
821   }
822         
823   for (Int_t j = 0; j < nt; j++) {
824     fEtaMulti->Fill(spdMult->GetEta(j));
825     fPhiMulti->Fill(spdMult->GetPhi(j));
826   }
827
828   Short_t chips0 = spdMult->GetNumberOfFiredChips(0);
829   Short_t chips1 = spdMult->GetNumberOfFiredChips(1);
830   //Printf("%s:%d",(char*)__FILE__,__LINE__);
831   Int_t inputCount = 0;
832   
833   if(spdVertex){
834     // get multiplicity from ITS tracklets
835     for (Int_t i=0; i<spdMult->GetNumberOfTracklets(); ++i){
836       Float_t deltaPhi = spdMult->GetDeltaPhi(i);
837       // prevent values to be shifted by 2 Pi()
838       if (deltaPhi < -TMath::Pi())
839         deltaPhi += TMath::Pi() * 2;
840       if (deltaPhi > TMath::Pi())
841         deltaPhi -= TMath::Pi() * 2;      
842       if (TMath::Abs(deltaPhi) > 1)
843         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, spdMult->GetTheta(i), spdMult->GetPhi(i), deltaPhi);
844       //trackled Deltaphi
845       Int_t label1=spdMult->GetLabel(i,0);
846       Int_t label2=spdMult->GetLabel(i,1);
847       if ( label1==label2){
848         if(stack->IsPhysicalPrimary(label1) == kTRUE)
849           fHistDeltaphiprimaries->Fill(deltaPhi);
850         else
851           fHistDeltaphisecondaries->Fill(deltaPhi);
852       }
853       else{
854         fHistDeltaphireject->Fill(deltaPhi);
855       }
856       ++inputCount;
857     }
858     // Printf("%s:%d",(char*)__FILE__,__LINE__);
859     fh2MultSpdChips[0]->Fill(chips0,inputCount);
860     fh2MultSpdChips[1]->Fill(chips1,inputCount);
861     //Printf("%s:%d",(char*)__FILE__,__LINE__);
862   }
863   // BKC
864
865   // MC analysis
866   Int_t nmc = MCEvent()->GetNumberOfTracks();
867   AliGenEventHeader*       header = MCEvent()->GenEventHeader();
868   // some test
869   AliGenCocktailEventHeader* cH = dynamic_cast<AliGenCocktailEventHeader*> (header);
870   AliGenPythiaEventHeader* pH;
871   if (cH == 0) 
872     {
873       header->Dump();
874       pH = dynamic_cast<AliGenPythiaEventHeader*>(header);
875     } else {
876     TObject* entry = cH->GetHeaders()->FindObject("Pythia");
877     pH = dynamic_cast<AliGenPythiaEventHeader*>(entry);
878   }
879   Int_t iproc = -1;
880
881   if (pH) iproc = pH->ProcessType();
882
883   TArrayF mcV;
884   header->PrimaryVertex(mcV);
885   Float_t vzMC = mcV[2];
886
887   // CKB
888   // Fill the correlation with MC
889   if(spdVertex&&MCEvent()){
890     Float_t diff = vzMC-vtxSPD[2];
891     fh2VertexCorrelation[kSPD]->Fill(vzMC,vtxSPD[2]);
892     fh2VertexCorrelationShift[kSPD]->Fill(vzMC,diff);
893     fh1VertexShift[kSPD]->Fill(diff);
894     if(spdVertex->GetZRes()>0)fh1VertexShiftNorm[kSPD]->Fill(diff/spdVertex->GetZRes());
895   }
896   if(tpcVertex&&MCEvent()){
897     Float_t diff = vzMC-vtxTPC[2];
898     fh2VertexCorrelation[kTPC]->Fill(vzMC,vtxTPC[2]);
899     fh2VertexCorrelationShift[kTPC]->Fill(vzMC,diff);
900     fh1VertexShift[kTPC]->Fill(diff);
901     if(tpcVertex->GetZRes()>0)fh1VertexShiftNorm[kTPC]->Fill(diff/tpcVertex->GetZRes());
902   }
903   // BKC
904
905
906   // MC event loop
907   //printf("MC particle loop \n");
908   for (Int_t i = 0; i < nmc; i++) {
909     AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(i);
910     //printf("MC particle loop %5d \n", i);
911     // Primaries only
912     if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
913     if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
914     // Int_t charge=Int_t(mcP->Particle()->GetPDG()->Charge());
915     Float_t phi = mcP->Phi();
916     Float_t eta = mcP->Eta();
917     Float_t pt  = mcP->Pt();
918     if (TMath::Abs(eta) > 0.9) continue;
919     Int_t ntr = mcP->GetNumberOfTrackReferences();
920     Int_t nITS = 0;
921     Int_t nTPC = 0;
922     Int_t nFRA = 0;
923     Float_t  x = 0.;
924     Float_t  y = 0.;
925       
926     for (Int_t j = 0; j < ntr; j++) {
927          
928       AliTrackReference* ref = mcP->GetTrackReference(j);
929         
930       if (ref->DetectorId() == 0) nITS++;
931       if (ref->DetectorId() == 1) nTPC++;
932       if (ref->DetectorId() == 2) nFRA++;
933       if (nTPC == 1) {
934         x = ref->X();
935         y = ref->Y();
936         break;
937       }
938     }
939
940     fHistMCpt ->Fill(pt);
941     fHistMCeta->Fill(eta);
942     fHistMCphi->Fill(phi);
943     Int_t index=ConvertePDG(mcP->PdgCode());
944     if(index<8){
945       fHistMCptCharge [index] ->Fill(pt);
946       fHistMCetaCharge[index]->Fill(eta);
947       fHistMCphiCharge[index]->Fill(phi);
948     }
949
950
951     if (pt > 2.) {
952       fHistMCHPTeta->Fill(eta);
953       fHistMCHPTphi->Fill(phi);
954     }
955
956     Bool_t reco = kFALSE;
957     Int_t multR =  0;
958     Int_t jold  = -1;
959       
960     for (Int_t j = 0; j < esdE->GetNumberOfTracks(); j++) {
961       if (i == labels.At(j)) {
962         reco = kTRUE;
963         multR++;
964         //AliESDtrack* test = esdE->GetTrack(j);
965         if (multR > 1) {
966           Int_t nclus = 0;
967           AliESDtrack* track = esdE->GetTrack(jold);
968           nclus = track->GetTPCclusters(0);              
969           fHistNCluster->Fill(nclus);
970           
971                   
972           track = esdE->GetTrack(j);
973           nclus = track->GetTPCclusters(0);
974           fHistNCluster->Fill(nclus);
975           
976
977         }
978         jold = j;
979       }
980     }
981
982     fHistRecMult->Fill(Float_t(multR), 1.);
983     if (multR == 0) {
984       fHistMCNRpt ->Fill(pt);
985       fHistMCNReta->Fill(eta);
986       fHistMCNRphi->Fill(phi);
987       if(index<8){
988         fHistMCNRptCharge [index] ->Fill(pt);
989         fHistMCNRetaCharge[index]->Fill(eta);
990         fHistMCNRphiCharge[index]->Fill(phi);
991       }
992      
993
994       if (pt > 2.) {
995         fHistMCNRHPTeta->Fill(eta);
996         fHistMCNRHPTphi->Fill(phi);
997       }
998     }
999     else if (multR > 1)
1000       {
1001         fHistMULTpt ->Fill(pt);
1002         fHistMULTeta->Fill(eta);
1003         fHistMULTphi->Fill(phi);
1004       } 
1005   } // MC particle loop
1006   //printf("End of MC particle loop \n");
1007   
1008   
1009 }      
1010
1011 //________________________________________________________________________
1012 void AliAnalysisTaskEfficiency::Terminate(Option_t *) 
1013 {
1014
1015 }  
1016
1017
1018 Bool_t AliAnalysisTaskEfficiency::SelectJouri(AliESDtrack* track) 
1019 {
1020
1021   Bool_t selected = kTRUE;
1022   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;    
1023   // > 50 Clusters
1024   if (track->GetTPCclusters(0) < 50) selected = kFALSE;
1025
1026   const AliESDVertex *vtx=esdE->GetPrimaryVertexSPD();
1027   if (!vtx->GetStatus()) selected = kFALSE;
1028     
1029   Double_t zv=vtx->GetZ();
1030
1031     
1032   const AliExternalTrackParam *ct=track->GetTPCInnerParam();
1033   if (!ct)  selected = kFALSE;
1034     
1035   Double_t xyz[3];
1036   ct->GetXYZ(xyz);
1037   Float_t rv = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
1038   if (rv > 3.0)                     selected = kFALSE;
1039   if (TMath::Abs(xyz[2] - zv)>2.5) selected = kFALSE;
1040   return (selected);
1041     
1042 }
1043
1044 Int_t AliAnalysisTaskEfficiency::GetIndexAC(AliESDtrack *track){
1045   if(!track)return -1;
1046
1047   // Crude selection for A and C side
1048   // just with eta
1049   if (track->Eta() > 0) { 
1050     if (track->Charge() > 0) 
1051       return kPosA;
1052     else
1053       return kNegA;     
1054   }
1055   else { // C side
1056     if (track->Charge() > 0) 
1057       return kPosC;
1058     else
1059       return kNegC;     
1060   }
1061   
1062   return -1;
1063
1064 }
1065
1066
1067 Int_t AliAnalysisTaskEfficiency::ConvertePDG(Int_t pdg)
1068 {
1069
1070   // function converts the pdg values for nuclei (d, t, 3He, 
1071   // 4He and their antiparticles) into indices for histo-array
1072   // filled in UserExec
1073   
1074   if      ( pdg ==  1000010020 )return 0;
1075   else if ( pdg == -1000010020 )return 1;
1076   else if ( pdg ==  1000010030 )return 2;
1077   else if ( pdg == -1000010030 )return 3;
1078   else if ( pdg ==  1000020030 )return 4;
1079   else if ( pdg == -1000020030 )return 5;
1080   else if ( pdg ==  1000020040 )return 6;
1081   else if ( pdg == -1000020040 )return 7;
1082   else return 9;
1083   
1084   
1085 }