d0675cba75582c77512fb2377fc24ddd7459e600
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEfficiencyBFPsi.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TH3D.h"
6 #include "TCanvas.h"
7 #include "TParticle.h"
8 #include "TObjArray.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
12
13 #include "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
33 ClassImp(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 //________________________________________________________________________
63 void AliAnalysisTaskEfficiencyBFPsi::UserCreateOutputObjects() {
64   // Create histograms
65   // Called once
66
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
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);
280
281   TH1::AddDirectory(oldStatus);
282
283 }
284
285 //________________________________________________________________________
286 void 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) {
300     AliError("ERROR: Could not retrieve MC event");
301     return;
302   }
303   AliStack* stack = mcEvent->Stack();
304   if (!stack) {
305     AliError("ERROR: Could not retrieve MC stack");
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
354             if(TMath::Abs(vertex->GetXv()) < fVxMax) {
355               if(TMath::Abs(vertex->GetYv()) < fVyMax) {
356                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
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) {
365                       AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
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;
378                     //if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin)) 
379                     //continue;
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;
428                       eta[nMCLabelCounter] = particle->Eta();
429                       pt[nMCLabelCounter]  = particle->Pt();
430                       phi[nMCLabelCounter] = particle->Phi()*TMath::RadToDeg();
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;
462                   
463                     Int_t mcGoods = nMCLabelCounter;
464                     for (Int_t k = 0; k < mcGoods; k++) {
465                       Int_t mcLabel = labelMCArray.At(k);
466                                                               
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;
479                       //if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin)) 
480                       //continue;
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 //________________________________________________________________________
604 void 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 //____________________________________________________________________//
616 Bool_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 }