9ab6861fc2ec3441860bd7cc751b867745b825a7
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEfficiencyBF.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TH3D.h"
6 #include "TCanvas.h"
7 #include "TParticle.h"
8 #include "TObjArray.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
12
13 #include "AliStack.h"
14 #include "AliMCEvent.h"
15 #include "AliESDEvent.h"
16 #include "AliESDInputHandler.h"
17 #include "AliESDtrackCuts.h"
18 #include "AliCentrality.h"
19 #include "AliLog.h"
20
21 #include "AliAnalysisTaskEfficiencyBF.h"
22
23 // ---------------------------------------------------------------------
24 //
25 // Task for calculating the efficiency of the Balance Function 
26 // for single particles and pairs
27 //
28 // Authors: Panos Christakoglou, Michael Weber
29 // 
30 // ---------------------------------------------------------------------
31
32 ClassImp(AliAnalysisTaskEfficiencyBF)
33
34 //________________________________________________________________________
35 AliAnalysisTaskEfficiencyBF::AliAnalysisTaskEfficiencyBF(const char *name) 
36   : AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0), 
37   fHistEventStats(0), fHistCentrality(0), fHistNMult(0), 
38   fHistGeneratedEtaPtPhiPlus(0), fHistFindableEtaPtPhiPlus(0), 
39   fHistReconstructedEtaPtPhiPlus(0), fHistSurvivedEtaPtPhiPlus(0),
40   fHistGeneratedEtaPtPhiMinus(0), fHistFindableEtaPtPhiMinus(0), 
41   fHistReconstructedEtaPtPhiMinus(0), fHistSurvivedEtaPtPhiMinus(0),
42   fHistGeneratedEtaPtPlusControl(0), fHistFindableEtaPtPlusControl(0), 
43   fHistReconstructedEtaPtPlusControl(0), fHistSurvivedEtaPtPlusControl(0),
44   fHistGeneratedEtaPtMinusControl(0), fHistFindableEtaPtMinusControl(0), 
45   fHistReconstructedEtaPtMinusControl(0), fHistSurvivedEtaPtMinusControl(0),
46   fHistGeneratedEtaPtPlusPlus(0), fHistFindableEtaPtPlusPlus(0), 
47   fHistReconstructedEtaPtPlusPlus(0), fHistSurvivedEtaPtPlusPlus(0),
48   fHistGeneratedEtaPtMinusMinus(0), fHistFindableEtaPtMinusMinus(0), 
49   fHistReconstructedEtaPtMinusMinus(0), fHistSurvivedEtaPtMinusMinus(0),
50   fHistGeneratedEtaPtPlusMinus(0), fHistFindableEtaPtPlusMinus(0), 
51   fHistReconstructedEtaPtPlusMinus(0), fHistSurvivedEtaPtPlusMinus(0),
52   fHistGeneratedPhiEtaPlusPlus(0), fHistFindablePhiEtaPlusPlus(0), 
53   fHistReconstructedPhiEtaPlusPlus(0), fHistSurvivedPhiEtaPlusPlus(0),
54   fHistGeneratedPhiEtaMinusMinus(0), fHistFindablePhiEtaMinusMinus(0), 
55   fHistReconstructedPhiEtaMinusMinus(0), fHistSurvivedPhiEtaMinusMinus(0),
56   fHistGeneratedPhiEtaPlusMinus(0), fHistFindablePhiEtaPlusMinus(0), 
57   fHistReconstructedPhiEtaPlusMinus(0), fHistSurvivedPhiEtaPlusMinus(0),
58   fESDtrackCuts(0), fAnalysisMode(0), 
59   fCentralityEstimator("V0M"), fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0), 
60   fVxMax(3.0), fVyMax(3.0), fVzMax(10.), 
61   fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), fMaxDCAxy(3.0), fMaxDCAz(3.0),
62     fMinPt(0.3), fMaxPt(1.5), fMinEta(-0.8), fMaxEta(0.8),fEtaRangeMin(0.0), fEtaRangeMax(1.6), fPtRangeMin(0.3), fPtRangeMax(1.5), fPhiRangeMin(0.0),fPhiRangeMax(360.),fdPhiRangeMax(180.), fEtaBin(100),fdEtaBin(64),fPtBin(49),fPhiBin(100),fdPhiBin(90)   {  
63   // Define input and output slots here
64   // Input slot #0 works with a TChain
65   DefineInput(0, TChain::Class());
66   // Output slot #0 id reserved by the base class for AOD
67   // Output slot #1 writes into a TH1 container
68   DefineOutput(1, TList::Class());
69   DefineOutput(2, TList::Class());
70 }
71
72 //________________________________________________________________________
73 void AliAnalysisTaskEfficiencyBF::UserCreateOutputObjects() {
74   // Create histograms
75   // Called once
76
77   // global switch disabling the reference 
78   // (to avoid "Replacing existing TH1" if several wagons are created in train)
79   Bool_t oldStatus = TH1::AddDirectoryStatus();
80   TH1::AddDirectory(kFALSE);
81
82   fQAList = new TList();
83   fQAList->SetName("QAList");
84   fQAList->SetOwner();
85
86   fOutputList = new TList();
87   fOutputList->SetName("OutputList");
88   fOutputList->SetOwner();
89
90   //Event stats.
91   TString gCutName[4] = {"Total","Offline trigger",
92                          "Vertex","Analyzed"};
93   fHistEventStats = new TH1F("fHistEventStats",
94                              "Event statistics;;N_{events}",
95                              4,0.5,4.5);
96   for(Int_t i = 1; i <= 4; i++)
97     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
98   fQAList->Add(fHistEventStats);
99
100   //ESD analysis
101   fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
102                              20,0.5,20.5);
103   fQAList->Add(fHistCentrality);
104   
105   //multiplicity (good MC tracks)
106   TString histName;
107   histName = "fHistNMult";
108   fHistNMult = new TH1F(histName.Data(), 
109                         ";N_{mult.}",
110                         200,0,20000);
111   fQAList->Add(fHistNMult);
112   
113   //eta vs pt for MC positives
114   fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
115                                      "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
116                                         fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
117   fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
118   fHistFindableEtaPtPhiPlus = new TH3D("fHistFindableEtaPtPhiPlus",
119                                      "Findable positive primaries;#eta;p_{T} (GeV/c);#phi",
120                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
121   fOutputList->Add(fHistFindableEtaPtPhiPlus);
122   fHistReconstructedEtaPtPhiPlus = new TH3D("fHistReconstructedEtaPtPhiPlus",
123                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi",
124                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
125   fOutputList->Add(fHistReconstructedEtaPtPhiPlus);
126   fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
127                                      "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
128                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
129   fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
130
131   //eta vs pt for MC negatives
132   fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
133                                      "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
134                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
135   fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
136   fHistFindableEtaPtPhiMinus = new TH3D("fHistFindableEtaPtPhiMinus",
137                                      "Findable positive primaries;#eta;p_{T} (GeV/c);#phi",
138                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
139   fOutputList->Add(fHistFindableEtaPtPhiMinus);
140   fHistReconstructedEtaPtPhiMinus = new TH3D("fHistReconstructedEtaPtPhiMinus",
141                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi",
142                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
143   fOutputList->Add(fHistReconstructedEtaPtPhiMinus);
144   fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
145                                      "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
146                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
147   fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
148
149   //eta vs pt for MC positives (control)
150   fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
151                                      "Generated positive primaries;#eta;p_{T} (GeV/c)",
152                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
153   fOutputList->Add(fHistGeneratedEtaPtPlusControl);
154   fHistFindableEtaPtPlusControl = new TH2F("fHistFindableEtaPtPlusControl",
155                                      "Findable positive primaries;#eta;p_{T} (GeV/c)",
156                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
157   fOutputList->Add(fHistFindableEtaPtPlusControl);
158   fHistReconstructedEtaPtPlusControl = new TH2F("fHistReconstructedEtaPtPlusControl",
159                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
160                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
161   fOutputList->Add(fHistReconstructedEtaPtPlusControl);
162   fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
163                                      "Survived positive primaries;#eta;p_{T} (GeV/c)",
164                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
165   fOutputList->Add(fHistSurvivedEtaPtPlusControl);
166
167   //eta vs pt for MC negatives (control)
168   fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
169                                      "Generated positive primaries;#eta;p_{T} (GeV/c)",
170                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
171   fOutputList->Add(fHistGeneratedEtaPtMinusControl);
172   fHistFindableEtaPtMinusControl = new TH2F("fHistFindableEtaPtMinusControl",
173                                      "Findable positive primaries;#eta;p_{T} (GeV/c)",
174                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
175   fOutputList->Add(fHistFindableEtaPtMinusControl);
176   fHistReconstructedEtaPtMinusControl = new TH2F("fHistReconstructedEtaPtMinusControl",
177                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
178                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
179   fOutputList->Add(fHistReconstructedEtaPtMinusControl);
180   fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
181                                      "Survived positive primaries;#eta;p_{T} (GeV/c)",
182                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
183   fOutputList->Add(fHistSurvivedEtaPtMinusControl);
184
185   //eta vs pt for MC ++
186   fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
187                                      "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
188                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
189   fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
190   fHistFindableEtaPtPlusPlus = new TH2F("fHistFindableEtaPtPlusPlus",
191                                      "Findable ++ primaries;#Delta#eta;p_{T} (GeV/c)",
192                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
193   fOutputList->Add(fHistFindableEtaPtPlusPlus);
194   fHistReconstructedEtaPtPlusPlus = new TH2F("fHistReconstructedEtaPtPlusPlus",
195                                      "Reconstructed ++ primaries;#Delta#eta;p_{T} (GeV/c)",
196                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
197   fOutputList->Add(fHistReconstructedEtaPtPlusPlus);
198   fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
199                                      "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
200                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
201   fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
202
203   //eta vs pt for MC --
204   fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
205                                      "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
206                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
207   fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
208   fHistFindableEtaPtMinusMinus = new TH2F("fHistFindableEtaPtMinusMinus",
209                                      "Findable -- primaries;#Delta#eta;p_{T} (GeV/c)",
210                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
211   fOutputList->Add(fHistFindableEtaPtMinusMinus);
212   fHistReconstructedEtaPtMinusMinus = new TH2F("fHistReconstructedEtaPtMinusMinus",
213                                      "Reconstructed -- primaries;#Delta#eta;p_{T} (GeV/c)",
214                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
215   fOutputList->Add(fHistReconstructedEtaPtMinusMinus);
216   fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
217                                      "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
218                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
219   fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
220
221   //eta vs pt for MC +-
222   fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
223                                      "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
224                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
225   fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
226   fHistFindableEtaPtPlusMinus = new TH2F("fHistFindableEtaPtPlusMinus",
227                                      "Findable +- primaries;#Delta#eta;p_{T} (GeV/c)",
228                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
229   fOutputList->Add(fHistFindableEtaPtPlusMinus);
230   fHistReconstructedEtaPtPlusMinus = new TH2F("fHistReconstructedEtaPtPlusMinus",
231                                      "Reconstructed +- primaries;#Delta#eta;p_{T} (GeV/c)",
232                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
233   fOutputList->Add(fHistReconstructedEtaPtPlusMinus);
234   fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
235                                      "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
236                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
237   fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
238
239   //=============================//
240    //phi vs eta for MC ++
241   fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
242                                      "Generated ++ primaries;#Delta#phi",
243                                           fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
244   fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
245   fHistFindablePhiEtaPlusPlus = new TH2F("fHistFindablePhiEtaPlusPlus",
246                                      "Findable ++ primaries;#Delta#phi;#Delta#eta",
247                                          fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
248   fOutputList->Add(fHistFindablePhiEtaPlusPlus);
249   fHistReconstructedPhiEtaPlusPlus = new TH2F("fHistReconstructedPhiEtaPlusPlus",
250                                      "Reconstructed ++ primaries;#Delta#phi;#Delta#eta",
251                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
252   fOutputList->Add(fHistReconstructedPhiEtaPlusPlus);
253   fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
254                                      "Survived ++ primaries;#Delta#phi;#Delta#eta",
255                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
256   fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
257
258   //phi vs eta for MC --
259   fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
260                                      "Generated -- primaries;#Delta#phi;#Delta#eta",
261                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
262   fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
263   fHistFindablePhiEtaMinusMinus = new TH2F("fHistFindablePhiEtaMinusMinus",
264                                      "Findable -- primaries;#Delta#phi;#Delta#eta",
265                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
266   fOutputList->Add(fHistFindablePhiEtaMinusMinus);
267   fHistReconstructedPhiEtaMinusMinus = new TH2F("fHistReconstructedPhiEtaMinusMinus",
268                                      "Reconstructed -- primaries;#Delta#phi;#Delta#eta",
269                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
270   fOutputList->Add(fHistReconstructedPhiEtaMinusMinus);
271   fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
272                                      "Survived -- primaries;#Delta#phi;#Delta#eta",
273                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
274   fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
275
276   //phi vs eta for MC +-
277   fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
278                                      "Generated +- primaries;#Delta#phi;#Delta#eta",
279                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
280   fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
281   fHistFindablePhiEtaPlusMinus = new TH2F("fHistFindablePhiEtaPlusMinus",
282                                      "Findable +- primaries;#Delta#phi;#Delta#eta",
283                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
284   fOutputList->Add(fHistFindablePhiEtaPlusMinus);
285   fHistReconstructedPhiEtaPlusMinus = new TH2F("fHistReconstructedPhiEtaPlusMinus",
286                                      "Reconstructed +- primaries;#Delta#phi;#Delta#eta",
287                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
288   fOutputList->Add(fHistReconstructedPhiEtaPlusMinus);
289   fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
290                                      "Survived +- primaries;#Delta#phi;#Delta#eta",
291                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
292   fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
293   //=============================//
294
295   fQAList->Print();
296   fOutputList->Print();
297
298   PostData(1, fQAList);
299   PostData(2, fOutputList);
300
301   TH1::AddDirectory(oldStatus);
302
303 }
304
305 //________________________________________________________________________
306 void AliAnalysisTaskEfficiencyBF::UserExec(Option_t *) {
307   // Main loop
308   // Called for each event
309
310   // Post output data.
311   //ESD analysis
312   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
313   if (!fESD) {
314     printf("ERROR: fESD not available\n");
315     return;
316   }
317   
318   AliMCEvent* mcEvent = MCEvent();
319   if (!mcEvent) {
320     AliError("ERROR: Could not retrieve MC event");
321     return;
322   }
323   AliStack* stack = mcEvent->Stack();
324   if (!stack) {
325     AliError("ERROR: Could not retrieve MC stack");
326     return;
327   }
328
329   // arrays for 2 particle histograms
330   Int_t nMCLabelCounter         = 0;
331   const Int_t maxMCLabelCounter = 20000;
332
333   Double_t eta[maxMCLabelCounter];
334   Double_t pt[maxMCLabelCounter];
335   Double_t phi[maxMCLabelCounter];
336   Int_t level[maxMCLabelCounter];
337   Int_t charge[maxMCLabelCounter];
338
339   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks()));
340   fHistEventStats->Fill(1); //all events
341     
342   //Centrality stuff
343   AliCentrality *centrality = fESD->GetCentrality();
344   Int_t nCentrality = 0;
345   nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.);
346
347   //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
348
349   if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
350                                           fCentralityPercentileMax,
351                                           fCentralityEstimator.Data())) {
352     fHistEventStats->Fill(2); //triggered + centrality
353     fHistCentrality->Fill(nCentrality+1);
354
355     //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
356   
357     if(fAnalysisMode.CompareTo("TPC") == 0 ) {
358       const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
359       if(vertex) {
360         if(vertex->GetNContributors() > 0) {
361           if(vertex->GetZRes() != 0) {
362             fHistEventStats->Fill(3); //events with a proper vertex
363             if(TMath::Abs(vertex->GetXv()) < fVxMax) {
364               if(TMath::Abs(vertex->GetYv()) < fVyMax) {
365                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
366                   fHistEventStats->Fill(4); //analyzed events
367                   
368                   Int_t nMCParticles = mcEvent->GetNumberOfTracks();
369                   TArrayI labelMCArray(nMCParticles);
370                   
371                   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
372                     AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks);
373                     if (!mcTrack) {
374                       AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
375                       continue;
376                     }
377                     
378                     //exclude particles generated out of the acceptance
379                     Double_t vz = mcTrack->Zv();
380                     if (TMath::Abs(vz) > 50.) continue;
381                    
382                     //acceptance
383                     if(TMath::Abs(mcTrack->Eta()) > fMaxEta) 
384                       continue;
385                     if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt)) 
386                       continue;
387                     if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin)) 
388                       continue;
389                     
390                     TParticle* particle = mcTrack->Particle();
391                     if(!particle) continue;
392                     if(!stack->IsPhysicalPrimary(iTracks)) continue;
393
394                     if(iTracks <= stack->GetNprimary()) {                     
395                       Short_t gMCCharge = mcTrack->Charge();
396                       Double_t phiRad = particle->Phi();
397                       Double_t phiDeg = phiRad*TMath::RadToDeg();       
398                         
399                       if(gMCCharge > 0)
400                         fHistGeneratedEtaPtPhiPlus->Fill(particle->Eta(),
401                                                          particle->Pt(),
402                                                          phiDeg);
403                       else if(gMCCharge < 0)
404                         fHistGeneratedEtaPtPhiMinus->Fill(particle->Eta(),
405                                                           particle->Pt(),
406                                                           phiDeg);
407
408                       
409                       // findable tracks --> DOES NOT WORK????
410                       // Loop over Track References
411                       Bool_t labelTPC = kTRUE;
412                       AliTrackReference* trackRef = 0;
413
414                       for (Int_t iTrackRef = 0; iTrackRef < mcTrack->GetNumberOfTrackReferences(); iTrackRef++) {
415                         trackRef = mcTrack->GetTrackReference(iTrackRef);
416                         if(trackRef) {
417                           Int_t detectorId = trackRef->DetectorId();
418                           if (detectorId == AliTrackReference::kTPC) {
419                             labelTPC = kTRUE;
420                             break;
421                           }
422                         }
423                       }//loop over track references
424
425                       if(labelTPC) {
426                         labelMCArray.AddAt(iTracks,nMCLabelCounter);
427
428                         if(nMCLabelCounter >= maxMCLabelCounter){
429                           AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
430                           break;
431                         }
432
433                         //fill the arrays for 2 particle analysis
434                         eta[nMCLabelCounter]    = particle->Eta();
435                         pt[nMCLabelCounter]     = particle->Pt();
436                         phi[nMCLabelCounter]     = particle->Phi()*TMath::RadToDeg();
437                         charge[nMCLabelCounter] = gMCCharge;
438                 
439                         // findable = generated in this case!
440
441                         level[nMCLabelCounter]  = 1;
442                         nMCLabelCounter += 1;
443                                 
444
445                         if(gMCCharge > 0)
446                           fHistFindableEtaPtPhiPlus->Fill(particle->Eta(),
447                                                           particle->Pt(),
448                                                           phiDeg);
449                         else if(gMCCharge < 0)
450                           fHistFindableEtaPtPhiMinus->Fill(particle->Eta(),
451                                                            particle->Pt(),
452                                                            phiDeg);
453                       }
454                     }//primaries
455                   }//loop over MC particles
456                 
457                   fHistNMult->Fill(nMCLabelCounter);
458                   
459                   // not used so far
460                   //Float_t dcaXY = 0.0, dcaZ = 0.0;
461
462                   //ESD track loop
463                   Int_t nGoodTracks = fESD->GetNumberOfTracks();
464                   
465                   TArrayI labelArray(nGoodTracks);
466                   Int_t labelCounter = 0;
467                   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
468                     AliESDtrack* track = fESD->GetTrack(iTracks);
469                     //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks);
470                     if(!track) continue;
471
472                     AliESDtrack *tpcOnlyTrack = new AliESDtrack();
473                     
474                     if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
475                       delete tpcOnlyTrack;
476                       continue;
477                     }
478
479                     Int_t label = TMath::Abs(track->GetTPCLabel());
480                     if(IsLabelUsed(labelArray,label)) continue;
481                     labelArray.AddAt(label,labelCounter);
482                     labelCounter += 1;
483
484                     Int_t mcGoods = nMCLabelCounter;
485                     for (Int_t k = 0; k < mcGoods; k++) {
486                       Int_t mcLabel = labelMCArray.At(k);
487                                               
488                       if (mcLabel != TMath::Abs(label)) continue;
489                       if(mcLabel != label) continue;
490                       if(label > stack->GetNtrack()) continue;
491                       
492                       TParticle *particle = stack->Particle(label);
493                       if(!particle) continue;
494                       
495                       //acceptance
496                       if(TMath::Abs(particle->Eta()) > fMaxEta) 
497                         continue;
498                       if((particle->Pt() > fMaxPt)||(particle->Pt() <  fMinPt)) 
499                         continue;
500                       if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin)) 
501                         continue;
502
503                       if(!stack->IsPhysicalPrimary(label)) continue;
504                       
505                       if(label <= stack->GetNprimary()) {
506                         
507                         // reconstructed
508                         level[k]  = 2;
509                         
510                         Short_t gCharge = track->Charge();
511                         Double_t phiRad = particle->Phi();
512                         Double_t phiDeg = phiRad*TMath::RadToDeg();
513
514                         if(gCharge > 0)
515                           fHistReconstructedEtaPtPhiPlus->Fill(particle->Eta(),
516                                                                particle->Pt(),
517                                                                phiDeg);
518                         else if(gCharge < 0)
519                           fHistReconstructedEtaPtPhiMinus->Fill(particle->Eta(),
520                                                                 particle->Pt(),
521                                                                 phiDeg);
522                         
523                         // track cuts + analysis kinematic cuts
524                         if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
525
526                           // survived
527                           level[k]  = 3;
528
529                           if(gCharge > 0)
530                             fHistSurvivedEtaPtPhiPlus->Fill(particle->Eta(),
531                                                             particle->Pt(),
532                                                             phiDeg);
533                           else if(gCharge < 0)
534                             fHistSurvivedEtaPtPhiMinus->Fill(particle->Eta(),
535                                                              particle->Pt(),
536                                                              phiDeg);
537                           
538                         }//track cuts
539                       }//primary particles
540                     }//findable track loop
541                   }//ESD track loop
542
543                   labelMCArray.Reset();
544                   labelArray.Reset();
545                   
546                   
547                 }//Vz cut
548               }//Vy cut
549             }//Vx cut
550           }//Vz resolution
551         }//number of contributors
552       }//valid vertex
553     }//TPC analysis mode
554   }//centrality  
555       
556
557
558   // Here comes the 2 particle analysis
559   // loop over all good MC particles
560   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
561     
562     // control 1D histograms (charge might be different?)
563     if(charge[i] > 0){
564       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
565       if(level[i] > 1) fHistReconstructedEtaPtPlusControl->Fill(eta[i],pt[i]);
566       if(level[i] > 2) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
567     }
568     else if(charge[i] < 0){
569       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
570       if(level[i] > 1) fHistReconstructedEtaPtMinusControl->Fill(eta[i],pt[i]);
571       if(level[i] > 2) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
572     }
573     
574     
575     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
576       
577       if(charge[i] > 0 && charge[j] > 0 ){
578         if(level[i] > 0 && level[j] > 0) {  
579           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
580           if (TMath::Abs(phi[i]-phi[j]) < 180)
581           fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
582         }
583         if(level[i] > 1 && level[j] > 1) {
584           fHistReconstructedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
585           if (TMath::Abs(phi[i]-phi[j]) < 180)
586           fHistReconstructedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
587         }
588         if(level[i] > 2 && level[j] > 2) {
589           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
590           if (TMath::Abs(phi[i]-phi[j]) < 180)
591           fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
592         }
593       }
594       
595       else if(charge[i] < 0 && charge[j] < 0 ){
596         if(level[i] > 0 && level[j] > 0) {
597           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);     
598           if (TMath::Abs(phi[i]-phi[j]) < 180)
599           fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));        
600         }
601         if(level[i] > 1 && level[j] > 1) {
602           fHistReconstructedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);   
603           fHistReconstructedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));           
604         }
605         if(level[i] > 2 && level[j] > 2) {
606           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
607           if (TMath::Abs(phi[i]-phi[j]) < 180)
608           fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
609         }
610       }
611       
612       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
613         if(level[i] > 0 && level[j] > 0) {
614           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
615           if (TMath::Abs(phi[i]-phi[j]) < 180)
616           fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));     
617         }
618         if(level[i] > 1 && level[j] > 1) {
619           fHistReconstructedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
620           fHistReconstructedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); 
621         }
622         if(level[i] > 2 && level[j] > 2) {
623           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
624           if (TMath::Abs(phi[i]-phi[j]) < 180)
625           fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
626         }             
627       }
628     }
629   }
630   
631   //Checking Entries for generated and survived
632   /*TH1F* GeneratedEtaPt = (TH1F*)fHistGeneratedEtaPtPlusMinus->ProjectionX("GeneratedEtaPt",5,15);
633     Int_t xGeneratedPt = fHistGeneratedEtaPtPlusMinus->GetNbinsX();
634     for (Int_t h=1; h < xGeneratedPt+1; h++){
635     Double_t binEntriesGenerated = GeneratedEtaPt->GetBinContent(h);
636     Printf("binEntriesGenerated: %lf - xGeneratedPt: %d",binEntriesGenerated,h);  
637     }
638     TH1F* GeneratedPhiEta = (TH1F*)fHistGeneratedPhiEtaPlusMinus->ProjectionY("GeneratedPhiEta",5,15);
639     Int_t yGeneratedPhi = fHistGeneratedPhiEtaPlusMinus->GetNbinsY();
640     for (Int_t h=1; h < yGeneratedPhi+1; h++){
641     Double_t binEntriesGenerated = GeneratedPhiEta->GetBinContent(h);
642     Printf("binEntriesGenerated: %lf - yGeneratedPhi: %d",binEntriesGenerated,h);  
643     }*/
644   
645   /*TH1F* SurvivedEtaPt = (TH1F*)fHistSurvivedEtaPtPlusMinus->ProjectionX("SurvivedEtaPt",5,15);
646     Int_t xSurvivedPt = fHistSurvivedEtaPtPlusMinus->GetNbinsX();
647     for (Int_t h=1; h < xSurvivedPt+1; h++){
648     Double_t binEntriesSurvived = SurvivedEtaPt->GetBinContent(h);
649     Printf("binEntriesSurvived: %lf - xSurvivedPt: %d",binEntriesSurvived,h);
650     }
651     TH1F* SurvivedPhiEta = (TH1F*)fHistSurvivedPhiEtaPlusMinus->ProjectionY("SurvivedPhiEta",5,15);
652     Int_t ySurvivedPhi = fHistSurvivedPhiEtaPlusMinus->GetNbinsY();
653     for (Int_t h=1; h < ySurvivedPhi+1; h++){
654     Double_t binEntriesSurvived = SurvivedPhiEta->GetBinContent(h);
655     Printf("binEntriesSurvived: %lf - ySurvivedPhi: %d",binEntriesSurvived,h);
656     }*/  
657 }
658   //________________________________________________________________________
659 void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) {
660   // Draw result to the screen
661   // Called once at the end of the query
662
663   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
664   if (!fOutputList) {
665     printf("ERROR: Output list not available\n");
666     return;
667   }
668 }
669
670 //____________________________________________________________________//
671 Bool_t AliAnalysisTaskEfficiencyBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
672   //Checks if the label is used already
673   Bool_t status = kFALSE;
674   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
675     if(labelArray.At(i) == label)
676       status = kTRUE;
677   }
678
679   return status;
680 }