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->GetX()) < fVxMax) {
 355               if(TMath::Abs(vertex->GetY()) < fVyMax) {
 356                 if(TMath::Abs(vertex->GetZ()) < 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 }