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