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 Int_t mcGoods = nMCLabelCounter;
485 for (Int_t k = 0; k < mcGoods; k++) {
486 Int_t mcLabel = labelMCArray.At(k);
488 if (mcLabel != TMath::Abs(label)) continue;
489 if(mcLabel != label) continue;
490 if(label > stack->GetNtrack()) continue;
492 TParticle *particle = stack->Particle(label);
493 if(!particle) continue;
496 if(TMath::Abs(particle->Eta()) > fMaxEta)
498 if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt))
500 if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin))
503 if(!stack->IsPhysicalPrimary(label)) continue;
505 if(label <= stack->GetNprimary()) {
510 Short_t gCharge = track->Charge();
511 Double_t phiRad = particle->Phi();
512 Double_t phiDeg = phiRad*TMath::RadToDeg();
515 fHistReconstructedEtaPtPhiPlus->Fill(particle->Eta(),
519 fHistReconstructedEtaPtPhiMinus->Fill(particle->Eta(),
523 // track cuts + analysis kinematic cuts
524 if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
530 fHistSurvivedEtaPtPhiPlus->Fill(particle->Eta(),
534 fHistSurvivedEtaPtPhiMinus->Fill(particle->Eta(),
540 }//findable track loop
543 labelMCArray.Reset();
551 }//number of contributors
558 // Here comes the 2 particle analysis
559 // loop over all good MC particles
560 for (Int_t i = 0; i < nMCLabelCounter ; i++) {
562 // control 1D histograms (charge might be different?)
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]);
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]);
575 for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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);
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);
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);
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);
658 //________________________________________________________________________
659 void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) {
660 // Draw result to the screen
661 // Called once at the end of the query
663 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
665 printf("ERROR: Output list not available\n");
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)