10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
14 #include "AliMCEvent.h"
15 #include "AliESDEvent.h"
16 #include "AliESDInputHandler.h"
17 #include "AliESDtrackCuts.h"
18 #include "AliCentrality.h"
21 #include "AliAnalysisTaskEfficiencyBF.h"
23 // ---------------------------------------------------------------------
25 // Task for calculating the efficiency of the Balance Function
26 // for single particles and pairs
28 // Authors: Panos Christakoglou, Michael Weber
30 // ---------------------------------------------------------------------
32 ClassImp(AliAnalysisTaskEfficiencyBF)
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());
72 //________________________________________________________________________
73 void AliAnalysisTaskEfficiencyBF::UserCreateOutputObjects() {
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);
82 fQAList = new TList();
83 fQAList->SetName("QAList");
86 fOutputList = new TList();
87 fOutputList->SetName("OutputList");
88 fOutputList->SetOwner();
91 TString gCutName[4] = {"Total","Offline trigger",
93 fHistEventStats = new TH1F("fHistEventStats",
94 "Event statistics;;N_{events}",
96 for(Int_t i = 1; i <= 4; i++)
97 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
98 fQAList->Add(fHistEventStats);
101 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
103 fQAList->Add(fHistCentrality);
105 //multiplicity (good MC tracks)
107 histName = "fHistNMult";
108 fHistNMult = new TH1F(histName.Data(),
111 fQAList->Add(fHistNMult);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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 //=============================//
296 fOutputList->Print();
298 PostData(1, fQAList);
299 PostData(2, fOutputList);
301 TH1::AddDirectory(oldStatus);
305 //________________________________________________________________________
306 void AliAnalysisTaskEfficiencyBF::UserExec(Option_t *) {
308 // Called for each event
312 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
314 printf("ERROR: fESD not available\n");
318 AliMCEvent* mcEvent = MCEvent();
320 AliError("ERROR: Could not retrieve MC event");
323 AliStack* stack = mcEvent->Stack();
325 AliError("ERROR: Could not retrieve MC stack");
329 // arrays for 2 particle histograms
330 Int_t nMCLabelCounter = 0;
331 const Int_t maxMCLabelCounter = 20000;
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];
339 //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks()));
340 fHistEventStats->Fill(1); //all events
343 AliCentrality *centrality = fESD->GetCentrality();
344 Int_t nCentrality = 0;
345 nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.);
347 //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
349 if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
350 fCentralityPercentileMax,
351 fCentralityEstimator.Data())) {
352 fHistEventStats->Fill(2); //triggered + centrality
353 fHistCentrality->Fill(nCentrality+1);
355 //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
357 if(fAnalysisMode.CompareTo("TPC") == 0 ) {
358 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
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
368 Int_t nMCParticles = mcEvent->GetNumberOfTracks();
369 TArrayI labelMCArray(nMCParticles);
371 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
372 AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks);
374 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
378 //exclude particles generated out of the acceptance
379 Double_t vz = mcTrack->Zv();
380 if (TMath::Abs(vz) > 50.) continue;
383 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
385 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
387 if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin))
390 TParticle* particle = mcTrack->Particle();
391 if(!particle) continue;
392 if(!stack->IsPhysicalPrimary(iTracks)) continue;
394 if(iTracks <= stack->GetNprimary()) {
395 Short_t gMCCharge = mcTrack->Charge();
396 Double_t phiRad = particle->Phi();
397 Double_t phiDeg = phiRad*TMath::RadToDeg();
400 fHistGeneratedEtaPtPhiPlus->Fill(particle->Eta(),
403 else if(gMCCharge < 0)
404 fHistGeneratedEtaPtPhiMinus->Fill(particle->Eta(),
409 // findable tracks --> DOES NOT WORK????
410 // Loop over Track References
411 Bool_t labelTPC = kTRUE;
412 AliTrackReference* trackRef = 0;
414 for (Int_t iTrackRef = 0; iTrackRef < mcTrack->GetNumberOfTrackReferences(); iTrackRef++) {
415 trackRef = mcTrack->GetTrackReference(iTrackRef);
417 Int_t detectorId = trackRef->DetectorId();
418 if (detectorId == AliTrackReference::kTPC) {
423 }//loop over track references
426 labelMCArray.AddAt(iTracks,nMCLabelCounter);
428 if(nMCLabelCounter >= maxMCLabelCounter){
429 AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
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;
439 // findable = generated in this case!
441 level[nMCLabelCounter] = 1;
442 nMCLabelCounter += 1;
446 fHistFindableEtaPtPhiPlus->Fill(particle->Eta(),
449 else if(gMCCharge < 0)
450 fHistFindableEtaPtPhiMinus->Fill(particle->Eta(),
455 }//loop over MC particles
457 fHistNMult->Fill(nMCLabelCounter);
460 //Float_t dcaXY = 0.0, dcaZ = 0.0;
463 Int_t nGoodTracks = fESD->GetNumberOfTracks();
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);
472 AliESDtrack *tpcOnlyTrack = new AliESDtrack();
474 if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
479 Int_t label = TMath::Abs(track->GetTPCLabel());
480 if(IsLabelUsed(labelArray,label)) continue;
481 labelArray.AddAt(label,labelCounter);
484 Bool_t iFound = kFALSE;
485 Int_t mcGoods = nMCLabelCounter;
486 for (Int_t k = 0; k < mcGoods; k++) {
487 Int_t mcLabel = labelMCArray.At(k);
490 if (mcLabel != TMath::Abs(label)) continue;
491 if(mcLabel != label) continue;
492 if(label > stack->GetNtrack()) continue;
494 TParticle *particle = stack->Particle(label);
495 if(!particle) continue;
498 if(TMath::Abs(particle->Eta()) > fMaxEta)
500 if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt))
502 if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin))
505 if(!stack->IsPhysicalPrimary(label)) continue;
507 if(label <= stack->GetNprimary()) {
512 Short_t gCharge = track->Charge();
513 Double_t phiRad = particle->Phi();
514 Double_t phiDeg = phiRad*TMath::RadToDeg();
517 fHistReconstructedEtaPtPhiPlus->Fill(particle->Eta(),
521 fHistReconstructedEtaPtPhiMinus->Fill(particle->Eta(),
525 // track cuts + analysis kinematic cuts
526 if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
532 fHistSurvivedEtaPtPhiPlus->Fill(particle->Eta(),
536 fHistSurvivedEtaPtPhiMinus->Fill(particle->Eta(),
542 }//findable track loop
545 labelMCArray.Reset();
553 }//number of contributors
560 // Here comes the 2 particle analysis
561 // loop over all good MC particles
562 for (Int_t i = 0; i < nMCLabelCounter ; i++) {
564 // control 1D histograms (charge might be different?)
566 if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
567 if(level[i] > 1) fHistReconstructedEtaPtPlusControl->Fill(eta[i],pt[i]);
568 if(level[i] > 2) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
570 else if(charge[i] < 0){
571 if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
572 if(level[i] > 1) fHistReconstructedEtaPtMinusControl->Fill(eta[i],pt[i]);
573 if(level[i] > 2) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
577 for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
579 if(charge[i] > 0 && charge[j] > 0 ){
580 if(level[i] > 0 && level[j] > 0) {
581 fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
582 if (TMath::Abs(phi[i]-phi[j]) < 180)
583 fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
585 if(level[i] > 1 && level[j] > 1) {
586 fHistReconstructedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
587 if (TMath::Abs(phi[i]-phi[j]) < 180)
588 fHistReconstructedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
590 if(level[i] > 2 && level[j] > 2) {
591 fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
592 if (TMath::Abs(phi[i]-phi[j]) < 180)
593 fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
597 else if(charge[i] < 0 && charge[j] < 0 ){
598 if(level[i] > 0 && level[j] > 0) {
599 fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
600 if (TMath::Abs(phi[i]-phi[j]) < 180)
601 fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
603 if(level[i] > 1 && level[j] > 1) {
604 fHistReconstructedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
605 fHistReconstructedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
607 if(level[i] > 2 && level[j] > 2) {
608 fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
609 if (TMath::Abs(phi[i]-phi[j]) < 180)
610 fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
614 else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
615 if(level[i] > 0 && level[j] > 0) {
616 fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
617 if (TMath::Abs(phi[i]-phi[j]) < 180)
618 fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
620 if(level[i] > 1 && level[j] > 1) {
621 fHistReconstructedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
622 fHistReconstructedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
624 if(level[i] > 2 && level[j] > 2) {
625 fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
626 if (TMath::Abs(phi[i]-phi[j]) < 180)
627 fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
633 //Checking Entries for generated and survived
634 /*TH1F* GeneratedEtaPt = (TH1F*)fHistGeneratedEtaPtPlusMinus->ProjectionX("GeneratedEtaPt",5,15);
635 Int_t xGeneratedPt = fHistGeneratedEtaPtPlusMinus->GetNbinsX();
636 for (Int_t h=1; h < xGeneratedPt+1; h++){
637 Double_t binEntriesGenerated = GeneratedEtaPt->GetBinContent(h);
638 Printf("binEntriesGenerated: %lf - xGeneratedPt: %d",binEntriesGenerated,h);
640 TH1F* GeneratedPhiEta = (TH1F*)fHistGeneratedPhiEtaPlusMinus->ProjectionY("GeneratedPhiEta",5,15);
641 Int_t yGeneratedPhi = fHistGeneratedPhiEtaPlusMinus->GetNbinsY();
642 for (Int_t h=1; h < yGeneratedPhi+1; h++){
643 Double_t binEntriesGenerated = GeneratedPhiEta->GetBinContent(h);
644 Printf("binEntriesGenerated: %lf - yGeneratedPhi: %d",binEntriesGenerated,h);
647 /*TH1F* SurvivedEtaPt = (TH1F*)fHistSurvivedEtaPtPlusMinus->ProjectionX("SurvivedEtaPt",5,15);
648 Int_t xSurvivedPt = fHistSurvivedEtaPtPlusMinus->GetNbinsX();
649 for (Int_t h=1; h < xSurvivedPt+1; h++){
650 Double_t binEntriesSurvived = SurvivedEtaPt->GetBinContent(h);
651 Printf("binEntriesSurvived: %lf - xSurvivedPt: %d",binEntriesSurvived,h);
653 TH1F* SurvivedPhiEta = (TH1F*)fHistSurvivedPhiEtaPlusMinus->ProjectionY("SurvivedPhiEta",5,15);
654 Int_t ySurvivedPhi = fHistSurvivedPhiEtaPlusMinus->GetNbinsY();
655 for (Int_t h=1; h < ySurvivedPhi+1; h++){
656 Double_t binEntriesSurvived = SurvivedPhiEta->GetBinContent(h);
657 Printf("binEntriesSurvived: %lf - ySurvivedPhi: %d",binEntriesSurvived,h);
660 //________________________________________________________________________
661 void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) {
662 // Draw result to the screen
663 // Called once at the end of the query
665 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
667 printf("ERROR: Output list not available\n");
672 //____________________________________________________________________//
673 Bool_t AliAnalysisTaskEfficiencyBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
674 //Checks if the label is used already
675 Bool_t status = kFALSE;
676 for(Int_t i = 0; i < labelArray.GetSize(); i++) {
677 if(labelArray.At(i) == label)