changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEfficiencyBFPsi.cxx
CommitLineData
74a4d448 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 "AliTHn.h"
14#include "AliStack.h"
15#include "AliMCEvent.h"
16#include "AliESDEvent.h"
17#include "AliESDInputHandler.h"
18#include "AliESDtrackCuts.h"
19#include "AliCentrality.h"
20#include "AliGenHijingEventHeader.h"
21
22#include "AliAnalysisTaskEfficiencyBFPsi.h"
23
24// ---------------------------------------------------------------------
25//
26// Task for calculating the efficiency of the Balance Function
27// for single particles and pairs
28//
29// Authors: Panos Christakoglou, Michael Weber
30//
31// ---------------------------------------------------------------------
32
33ClassImp(AliAnalysisTaskEfficiencyBFPsi)
34
35//________________________________________________________________________
36 AliAnalysisTaskEfficiencyBFPsi::AliAnalysisTaskEfficiencyBFPsi(const char *name):
37 AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0),
38 fHistEventStats(0), fHistCentrality(0), fHistNMult(0),
39 fHistGeneratedPlus(0), fHistSurvivedPlus(0),
40 fHistGeneratedMinus(0), fHistSurvivedMinus(0),
41 fHistGeneratedPlusPlus(0), fHistSurvivedPlusPlus(0),
42 fHistGeneratedMinusMinus(0), fHistSurvivedMinusMinus(0),
43 fHistGeneratedPlusMinus(0), fHistSurvivedPlusMinus(0),
44 fHistGeneratedMinusPlus(0), fHistSurvivedMinusPlus(0),
45 fESDtrackCuts(0), fAnalysisMode(0),
46 fCentralityEstimator("V0M"),
47 fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0),
48 fVxMax(3.0), fVyMax(3.0), fVzMax(10.),
49 fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0),
50 fMaxDCAxy(3.0), fMaxDCAz(3.0),
51 fMinPt(0.3), fMaxPt(1.5), fMaxEta(0.8), fEtaRangeMax(0.8),
52 fPtRangeMin(0.1), fPtRangeMax(5.0), fPhiRangeMin(0.0),fPhiRangeMax(6.28) {
53 // Define input and output slots here
54 // Input slot #0 works with a TChain
55 DefineInput(0, TChain::Class());
56 // Output slot #0 id reserved by the base class for AOD
57 // Output slot #1 writes into a TH1 container
58 DefineOutput(1, TList::Class());
59 DefineOutput(2, TList::Class());
60}
61
62//________________________________________________________________________
63void AliAnalysisTaskEfficiencyBFPsi::UserCreateOutputObjects() {
64 // Create histograms
65 // Called once
66
211b716d 67 // global switch disabling the reference
68 // (to avoid "Replacing existing TH1" if several wagons are created in train)
69 Bool_t oldStatus = TH1::AddDirectoryStatus();
70 TH1::AddDirectory(kFALSE);
71
74a4d448 72 fQAList = new TList();
73 fQAList->SetName("QAList");
74 fQAList->SetOwner();
75
76 fOutputList = new TList();
77 fOutputList->SetName("OutputList");
78 fOutputList->SetOwner();
79
80 //Event stats.
81 TString gCutName[4] = {"Total","Offline trigger",
82 "Vertex","Analyzed"};
83 fHistEventStats = new TH1F("fHistEventStats",
84 "Event statistics;;N_{events}",
85 4,0.5,4.5);
86 for(Int_t i = 1; i <= 4; i++)
87 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
88 fQAList->Add(fHistEventStats);
89
90 //ESD analysis
91 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
92 20,0.5,20.5);
93 fQAList->Add(fHistCentrality);
94
95 //multiplicity (good MC tracks)
96 TString histName;
97 histName = "fHistNMult";
98 fHistNMult = new TH1F(histName.Data(),
99 ";N_{mult.}",
100 200,0,20000);
101 fQAList->Add(fHistNMult);
102
103 //Setting up the AliTHn
104 Int_t anaSteps = 1; // analysis steps
105 Int_t iBinSingle[kVariablesSingle]; // binning for track variables
106 Double_t* dBinsSingle[kVariablesSingle]; // bins for track variables
107 TString axisTitleSingle[kVariablesSingle]; // axis titles for track variables
108
109 // two particle histograms
110 Int_t iBinPair[kVariablesPair]; // binning for track variables
111 Double_t* dBinsPair[kVariablesPair]; // bins for track variables
112 TString axisTitlePair[kVariablesPair]; // axis titles for track variables
113
114 //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (all)
115 const Int_t kNPsi2Bins = 4;
116 Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
117 iBinSingle[0] = kNPsi2Bins;
118 dBinsSingle[0] = psi2Bins;
119 axisTitleSingle[0] = "#phi - #Psi_{2} (a.u.)";
120 iBinPair[0] = kNPsi2Bins;
121 dBinsPair[0] = psi2Bins;
122 axisTitlePair[0] = "#phi - #Psi_{2} (a.u.)";
123
124 // delta eta
125 const Int_t kNDeltaEtaBins = 80;
126 Double_t deltaEtaBins[kNDeltaEtaBins+1];
127 for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
128 deltaEtaBins[i] = -2.0 + i * 0.05;
129 iBinPair[1] = kNDeltaEtaBins;
130 dBinsPair[1] = deltaEtaBins;
131 axisTitlePair[1] = "#Delta #eta";
132
133 // delta phi
134 const Int_t kNDeltaPhiBins = 72;
135 Double_t deltaPhiBins[kNDeltaPhiBins+1];
136 for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
137 deltaPhiBins[i] = -180.0 + i * 5.;
138 }
139 iBinPair[2] = kNDeltaPhiBins;
140 dBinsPair[2] = deltaPhiBins;
141 axisTitlePair[2] = "#Delta #phi (#circ)";
142
143 // pt(trigger-associated)
144 const Int_t kNPtBins = 16;
145 Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
146 iBinSingle[1] = kNPtBins;
147 dBinsSingle[1] = ptBins;
148 axisTitleSingle[1] = "p_{t}^{trig.} (GeV/c)";
149
150 iBinPair[3] = kNPtBins;
151 dBinsPair[3] = ptBins;
152 axisTitlePair[3] = "p_{t}^{trig.} (GeV/c)";
153
154 iBinPair[4] = kNPtBins;
155 dBinsPair[4] = ptBins;
156 axisTitlePair[4] = "p_{t}^{assoc.} (GeV/c)";
157
158 //=============================//
159 //Generated: Single particle distributions
160 fHistGeneratedPlus = new AliTHn("fHistGeneratedPlus",
161 "Generated positive primaries",
162 anaSteps,kVariablesSingle,iBinSingle);
163 for (Int_t j = 0; j < kVariablesSingle; j++) {
164 fHistGeneratedPlus->SetBinLimits(j, dBinsSingle[j]);
165 fHistGeneratedPlus->SetVarTitle(j, axisTitleSingle[j]);
166 }
167 fOutputList->Add(fHistGeneratedPlus);
168
169 fHistGeneratedMinus = new AliTHn("fHistGeneratedMinus",
170 "Generated positive primaries",
171 anaSteps,kVariablesSingle,iBinSingle);
172 for (Int_t j = 0; j < kVariablesSingle; j++) {
173 fHistGeneratedMinus->SetBinLimits(j, dBinsSingle[j]);
174 fHistGeneratedMinus->SetVarTitle(j, axisTitleSingle[j]);
175 }
176 fOutputList->Add(fHistGeneratedMinus);
177
178 //Survived: Single particle distributions
179 fHistSurvivedPlus = new AliTHn("fHistSurvivedPlus",
180 "Survived positive primaries",
181 anaSteps,kVariablesSingle,iBinSingle);
182 for (Int_t j = 0; j < kVariablesSingle; j++) {
183 fHistSurvivedPlus->SetBinLimits(j, dBinsSingle[j]);
184 fHistSurvivedPlus->SetVarTitle(j, axisTitleSingle[j]);
185 }
186 fOutputList->Add(fHistSurvivedPlus);
187
188 fHistSurvivedMinus = new AliTHn("fHistSurvivedMinus",
189 "Survived positive primaries",
190 anaSteps,kVariablesSingle,iBinSingle);
191 for (Int_t j = 0; j < kVariablesSingle; j++) {
192 fHistSurvivedMinus->SetBinLimits(j, dBinsSingle[j]);
193 fHistSurvivedMinus->SetVarTitle(j, axisTitleSingle[j]);
194 }
195 fOutputList->Add(fHistSurvivedMinus);
196
197 //=============================//
198 //Generated-Survived: Particle pairs +-
199 fHistGeneratedPlusMinus = new AliTHn("fHistGeneratedPlusMinus",
200 "Generated +- primaries",
201 anaSteps,kVariablesPair,iBinPair);
202 for (Int_t j = 0; j < kVariablesPair; j++) {
203 fHistGeneratedPlusMinus->SetBinLimits(j, dBinsPair[j]);
204 fHistGeneratedPlusMinus->SetVarTitle(j, axisTitlePair[j]);
205 }
206 fOutputList->Add(fHistGeneratedPlusMinus);
207
208 fHistSurvivedPlusMinus = new AliTHn("fHistSurvivedPlusMinus",
209 "Survived +- primaries",
210 anaSteps,kVariablesPair,iBinPair);
211 for (Int_t j = 0; j < kVariablesPair; j++) {
212 fHistSurvivedPlusMinus->SetBinLimits(j, dBinsPair[j]);
213 fHistSurvivedPlusMinus->SetVarTitle(j, axisTitlePair[j]);
214 }
215 fOutputList->Add(fHistSurvivedPlusMinus);
216
217 //Generated-Survived: Particle pairs -+
218 fHistGeneratedMinusPlus = new AliTHn("fHistGeneratedMinusPlus",
219 "Generated -+ primaries",
220 anaSteps,kVariablesPair,iBinPair);
221 for (Int_t j = 0; j < kVariablesPair; j++) {
222 fHistGeneratedMinusPlus->SetBinLimits(j, dBinsPair[j]);
223 fHistGeneratedMinusPlus->SetVarTitle(j, axisTitlePair[j]);
224 }
225 fOutputList->Add(fHistGeneratedMinusPlus);
226
227 fHistSurvivedMinusPlus = new AliTHn("fHistSurvivedMinusPlus",
228 "Survived -+ primaries",
229 anaSteps,kVariablesPair,iBinPair);
230 for (Int_t j = 0; j < kVariablesPair; j++) {
231 fHistSurvivedMinusPlus->SetBinLimits(j, dBinsPair[j]);
232 fHistSurvivedMinusPlus->SetVarTitle(j, axisTitlePair[j]);
233 }
234 fOutputList->Add(fHistSurvivedMinusPlus);
235
236 //Generated-Survived: Particle pairs ++
237 fHistGeneratedPlusPlus = new AliTHn("fHistGeneratedPlusPlus",
238 "Generated ++ primaries",
239 anaSteps,kVariablesPair,iBinPair);
240 for (Int_t j = 0; j < kVariablesPair; j++) {
241 fHistGeneratedPlusPlus->SetBinLimits(j, dBinsPair[j]);
242 fHistGeneratedPlusPlus->SetVarTitle(j, axisTitlePair[j]);
243 }
244 fOutputList->Add(fHistGeneratedPlusPlus);
245
246 fHistSurvivedPlusPlus = new AliTHn("fHistSurvivedPlusPlus",
247 "Survived ++ primaries",
248 anaSteps,kVariablesPair,iBinPair);
249 for (Int_t j = 0; j < kVariablesPair; j++) {
250 fHistSurvivedPlusPlus->SetBinLimits(j, dBinsPair[j]);
251 fHistSurvivedPlusPlus->SetVarTitle(j, axisTitlePair[j]);
252 }
253 fOutputList->Add(fHistSurvivedPlusPlus);
254
255 //Generated-Survived: Particle pairs --
256 fHistGeneratedMinusMinus = new AliTHn("fHistGeneratedMinusMinus",
257 "Generated -- primaries",
258 anaSteps,kVariablesPair,iBinPair);
259 for (Int_t j = 0; j < kVariablesPair; j++) {
260 fHistGeneratedMinusMinus->SetBinLimits(j, dBinsPair[j]);
261 fHistGeneratedMinusMinus->SetVarTitle(j, axisTitlePair[j]);
262 }
263 fOutputList->Add(fHistGeneratedMinusMinus);
264
265 fHistSurvivedMinusMinus = new AliTHn("fHistSurvivedMinusMinus",
266 "Survived -- primaries",
267 anaSteps,kVariablesPair,iBinPair);
268 for (Int_t j = 0; j < kVariablesPair; j++) {
269 fHistSurvivedMinusMinus->SetBinLimits(j, dBinsPair[j]);
270 fHistSurvivedMinusMinus->SetVarTitle(j, axisTitlePair[j]);
271 }
272 fOutputList->Add(fHistSurvivedMinusMinus);
273 //=============================//
274
275 fQAList->Print();
276 fOutputList->Print();
277
278 PostData(1, fQAList);
279 PostData(2, fOutputList);
211b716d 280
281 TH1::AddDirectory(oldStatus);
282
74a4d448 283}
284
285//________________________________________________________________________
286void AliAnalysisTaskEfficiencyBFPsi::UserExec(Option_t *) {
287 // Main loop
288 // Called for each event
289
290 // Post output data.
291 //ESD analysis
292 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
293 if (!fESD) {
294 printf("ERROR: fESD not available\n");
295 return;
296 }
297
298 AliMCEvent* mcEvent = MCEvent();
299 if (!mcEvent) {
6d6eb2df 300 AliError("ERROR: Could not retrieve MC event");
74a4d448 301 return;
302 }
303 AliStack* stack = mcEvent->Stack();
304 if (!stack) {
6d6eb2df 305 AliError("ERROR: Could not retrieve MC stack");
74a4d448 306 return;
307 }
308
309 // arrays for 2 particle histograms
310 Int_t nMCLabelCounter = 0;
311 const Int_t maxMCLabelCounter = 20000;
312
313 Double_t phiMinusPsi[maxMCLabelCounter];
314 Double_t eta[maxMCLabelCounter];
315 Double_t pt[maxMCLabelCounter];
316 Double_t phi[maxMCLabelCounter];
317 Int_t level[maxMCLabelCounter];
318 Int_t charge[maxMCLabelCounter];
319
320 Double_t trackVariablesSingle[kVariablesSingle];
321 Double_t trackVariablesPair[kVariablesPair];
322
323 Double_t gReactionPlane = 0.;
324 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(mcEvent)->GenEventHeader());
325 if (headerH) {
326 gReactionPlane = headerH->ReactionPlaneAngle();
327 gReactionPlane *= TMath::RadToDeg();
328 }
329
330 //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks()));
331 fHistEventStats->Fill(1); //all events
332
333 //Centrality stuff
334 AliCentrality *centrality = fESD->GetCentrality();
335 Int_t nCentrality = 0;
336 nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.);
337
338 //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
339
340 if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
341 fCentralityPercentileMax,
342 fCentralityEstimator.Data())) {
343 fHistEventStats->Fill(2); //triggered + centrality
344 fHistCentrality->Fill(nCentrality+1);
345
346 //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
347
348 if(fAnalysisMode.CompareTo("TPC") == 0 ) {
349 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
350 if(vertex) {
351 if(vertex->GetNContributors() > 0) {
352 if(vertex->GetZRes() != 0) {
353 fHistEventStats->Fill(3); //events with a proper vertex
e690d4d0 354 if(TMath::Abs(vertex->GetX()) < fVxMax) {
355 if(TMath::Abs(vertex->GetY()) < fVyMax) {
356 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
74a4d448 357 fHistEventStats->Fill(4); //analyzed events
358
359 Int_t nMCParticles = mcEvent->GetNumberOfTracks();
360 TArrayI labelMCArray(nMCParticles);
361
362 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
363 AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks);
364 if (!mcTrack) {
6d6eb2df 365 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
74a4d448 366 continue;
367 }
368
369 //exclude particles generated out of the acceptance
370 Double_t vz = mcTrack->Zv();
371 if (TMath::Abs(vz) > 50.) continue;
372
373 //acceptance
374 if(TMath::Abs(mcTrack->Eta()) > fEtaRangeMax)
375 continue;
376 if((mcTrack->Pt() > fPtRangeMax)||(mcTrack->Pt() < fPtRangeMin))
377 continue;
d130918a 378 //if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin))
379 //continue;
74a4d448 380
381 TParticle* particle = mcTrack->Particle();
382 if(!particle) continue;
383 if(!stack->IsPhysicalPrimary(iTracks)) continue;
384
385 if(iTracks <= stack->GetNprimary()) {
386 Short_t gMCCharge = mcTrack->Charge();
387 Float_t firstPhi = mcTrack->Phi()*TMath::RadToDeg();
388 Float_t firstPt = mcTrack->Pt();
389
390 // Event plane (determine psi bin)
391 Double_t gPsiMinusPhi = 0.;
392 Double_t gPsiMinusPhiBin = -10.;
393 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
394 //in-plane
395 if((gPsiMinusPhi <= 7.5)||
396 ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
397 gPsiMinusPhiBin = 0.0;
398 //intermediate
399 else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
400 ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
401 ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
402 ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
403 gPsiMinusPhiBin = 1.0;
404 //out of plane
405 else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
406 ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
407 gPsiMinusPhiBin = 2.0;
408 //everything else
409 else
410 gPsiMinusPhiBin = 3.0;
411
412 trackVariablesSingle[0] = gPsiMinusPhiBin;
413 trackVariablesSingle[1] = firstPt;
414
415 if(gMCCharge > 0)
416 fHistGeneratedPlus->Fill(trackVariablesSingle,0,1.);
417 else if(gMCCharge < 0)
418 fHistGeneratedMinus->Fill(trackVariablesSingle,0,1.);
419
420 labelMCArray.AddAt(iTracks,nMCLabelCounter);
421 if(nMCLabelCounter >= maxMCLabelCounter){
422 AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
423 break;
424 }
425
426 //fill the arrays for 2 particle analysis
427 phiMinusPsi[nMCLabelCounter] = gPsiMinusPhiBin;
d130918a 428 eta[nMCLabelCounter] = particle->Eta();
429 pt[nMCLabelCounter] = particle->Pt();
430 phi[nMCLabelCounter] = particle->Phi()*TMath::RadToDeg();
74a4d448 431 charge[nMCLabelCounter] = gMCCharge;
432 // findable = generated in this case!
433
434 level[nMCLabelCounter] = 1;
435 nMCLabelCounter += 1;
436 }//primaries
437 }//loop over MC particles
438
439 fHistNMult->Fill(nMCLabelCounter);
440
441 //ESD track loop
442 Int_t nGoodTracks = fESD->GetNumberOfTracks();
443
444 TArrayI labelArray(nGoodTracks);
445 Int_t labelCounter = 0;
446 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
447 AliESDtrack* track = fESD->GetTrack(iTracks);
448 //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks);
449 if(!track) continue;
450
451 AliESDtrack *tpcOnlyTrack = new AliESDtrack();
452
453 if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
454 delete tpcOnlyTrack;
455 continue;
456 }
457
458 Int_t label = TMath::Abs(track->GetTPCLabel());
459 if(IsLabelUsed(labelArray,label)) continue;
460 labelArray.AddAt(label,labelCounter);
461 labelCounter += 1;
fed402ec 462
74a4d448 463 Int_t mcGoods = nMCLabelCounter;
464 for (Int_t k = 0; k < mcGoods; k++) {
465 Int_t mcLabel = labelMCArray.At(k);
fed402ec 466
74a4d448 467 if (mcLabel != TMath::Abs(label)) continue;
468 if(mcLabel != label) continue;
469 if(label > stack->GetNtrack()) continue;
470
471 TParticle *particle = stack->Particle(label);
472 if(!particle) continue;
473
474 //acceptance
475 if(TMath::Abs(particle->Eta()) > fEtaRangeMax)
476 continue;
477 if((particle->Pt() > fPtRangeMax)||(particle->Pt() < fPtRangeMin))
478 continue;
d130918a 479 //if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin))
480 //continue;
74a4d448 481
482 if(!stack->IsPhysicalPrimary(label)) continue;
483
484 if(label <= stack->GetNprimary()) {
485 Short_t gCharge = track->Charge();
486 // track cuts + analysis kinematic cuts
487 if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
488 // survived
489 level[k] = 2;
490
491 Float_t firstPhi = particle->Phi()*TMath::RadToDeg();
492 Float_t firstPt = particle->Pt();
493
494 // Event plane (determine psi bin)
495 Double_t gPsiMinusPhi = 0.;
496 Double_t gPsiMinusPhiBin = -10.;
497 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
498 //in-plane
499 if((gPsiMinusPhi <= 7.5)||
500 ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
501 gPsiMinusPhiBin = 0.0;
502 //intermediate
503 else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
504 ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
505 ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
506 ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
507 gPsiMinusPhiBin = 1.0;
508 //out of plane
509 else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
510 ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
511 gPsiMinusPhiBin = 2.0;
512 //everything else
513 else
514 gPsiMinusPhiBin = 3.0;
515
516 trackVariablesSingle[0] = gPsiMinusPhiBin;
517 trackVariablesSingle[1] = firstPt;
518
519 if(gCharge > 0)
520 fHistSurvivedPlus->Fill(trackVariablesSingle,0,1.);
521 else if(gCharge < 0)
522 fHistSurvivedMinus->Fill(trackVariablesSingle,0,1.);
523 }//track cuts
524 }//primary particles
525 }//loop over the generated
526 }//ESD track loop
527
528 labelMCArray.Reset();
529 labelArray.Reset();
530 }//Vz cut
531 }//Vy cut
532 }//Vx cut
533 }//Vz resolution
534 }//number of contributors
535 }//valid vertex
536 }//TPC analysis mode
537 }//centrality
538
539 // Here comes the 2 particle analysis
540 // loop over all good MC particles
541 for (Int_t i = 0; i < nMCLabelCounter ; i++) {
542 Float_t firstEta = eta[i];
543 Float_t firstPhi = phi[i];
544 Float_t firstPt = pt[i];
545 Float_t gPhisMinusPsiBin = phiMinusPsi[i];
546 for (Int_t j = 0; j < nMCLabelCounter ; j++) {
547 if(i == j) continue;
548
549 Float_t secondEta = eta[j];
550 Float_t secondPhi = phi[j];
551 Float_t secondPt = pt[j];
552
553 trackVariablesPair[0] = gPhisMinusPsiBin;
554 trackVariablesPair[1] = firstEta - secondEta; // delta eta
555 trackVariablesPair[2] = firstPhi - secondPhi; // delta phi
556 if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
557 trackVariablesPair[2] -= 360.;
558 if (trackVariablesPair[2] < - 180.)
559 trackVariablesPair[2] += 360.;
560 trackVariablesPair[3] = firstPt; // pt trigger
561 trackVariablesPair[4] = secondPt; // pt
562
563 //++ pairs
564 if(charge[i] > 0 && charge[j] > 0 ) {
565 if(level[i] > 0 && level[j] > 0)
566 fHistGeneratedPlusPlus->Fill(trackVariablesPair,0,1.);
567
568 if(level[i] > 1 && level[j] > 1)
569 fHistSurvivedPlusPlus->Fill(trackVariablesPair,0,1.);
570 }
571
572 //-- pairs
573 else if(charge[i] < 0 && charge[j] < 0 ) {
574 if(level[i] > 0 && level[j] > 0)
575 fHistGeneratedMinusMinus->Fill(trackVariablesPair,0,1.);
576
577 if(level[i] > 1 && level[j] > 1)
578 fHistSurvivedMinusMinus->Fill(trackVariablesPair,0,1.);
579 }
580
581 //+- pairs
582 else if(charge[i] > 0 && charge[j] < 0 ) {
583 if(level[i] > 0 && level[j] > 0)
584 fHistGeneratedPlusMinus->Fill(trackVariablesPair,0,1.);
585
586 if(level[i] > 1 && level[j] > 1)
587 fHistSurvivedPlusMinus->Fill(trackVariablesPair,0,1.);
588 }
589
590 //-+ pairs
591 else if(charge[i] < 0 && charge[j] > 0 ) {
592 if(level[i] > 0 && level[j] > 0)
593 fHistGeneratedMinusPlus->Fill(trackVariablesPair,0,1.);
594
595 if(level[i] > 1 && level[j] > 1)
596 fHistSurvivedMinusPlus->Fill(trackVariablesPair,0,1.);
597 }
598 }//second particle loop
599 }//first particle loop
600
601}
602
603//________________________________________________________________________
604void AliAnalysisTaskEfficiencyBFPsi::Terminate(Option_t *) {
605 // Draw result to the screen
606 // Called once at the end of the query
607
608 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
609 if (!fOutputList) {
610 printf("ERROR: Output list not available\n");
611 return;
612 }
613}
614
615//____________________________________________________________________//
616Bool_t AliAnalysisTaskEfficiencyBFPsi::IsLabelUsed(TArrayI labelArray, Int_t label) {
617 //Checks if the label is used already
618 Bool_t status = kFALSE;
619 for(Int_t i = 0; i < labelArray.GetSize(); i++) {
620 if(labelArray.At(i) == label)
621 status = kTRUE;
622 }
623
624 return status;
625}