PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEffContBF.cxx

   1 #include "TChain.h"
   2 #include "TList.h"
   3 #include "TTree.h"
   4 #include "TH1F.h"
   5 #include "TH2F.h"
   6 #include "TH3D.h"
   7 #include "TCanvas.h"
   8 #include "TParticle.h"
   9 #include "TObjArray.h"
  10
  11 #include "AliAnalysisTask.h"
  12 #include "AliAnalysisManager.h"
  13
  14 #include "AliStack.h"
  15 #include "AliMCEvent.h"
  16 #include "AliAODEvent.h"
  17 #include "AliAODInputHandler.h"
  18 #include "AliAODMCParticle.h"
  19 #include "AliCentrality.h"
  20
  21 #include "AliAnalysisTaskEffContBF.h"
  22
  23 // ---------------------------------------------------------------------
  24 //
  25 // Task for calculating the efficiency of the Balance Function
  26 // for single particles and pairs
  27 //
  28 // ---------------------------------------------------------------------
  29
  30 ClassImp(AliAnalysisTaskEffContBF)
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskEffContBF::AliAnalysisTaskEffContBF(const char *name)
  34   : AliAnalysisTaskSE(name),
  35     fAOD(0),
  36     fArrayMC(0),
  37     fQAList(0),
  38     fOutputList(0),
  39     fHistEventStats(0),
  40     fHistCentrality(0),
  41     fHistNMult(0),
  42     fHistVz(0),
  43     fHistContaminationSecondariesPlus(0),
  44     fHistContaminationSecondariesMinus(0), //
  45     fHistContaminationPrimariesPlus(0),
  46     fHistContaminationPrimariesMinus(0), //
  47     fHistGeneratedEtaPtPhiPlus(0),
  48     fHistSurvivedEtaPtPhiPlus(0),
  49     fHistGeneratedEtaPtPhiMinus(0),
  50     fHistSurvivedEtaPtPhiMinus(0),
  51     fHistGeneratedEtaPtPlusControl(0),
  52     fHistSurvivedEtaPtPlusControl(0),
  53     fHistGeneratedEtaPtMinusControl(0),
  54     fHistSurvivedEtaPtMinusControl(0),
  55     fHistGeneratedEtaPtPlusPlus(0),
  56     fHistSurvivedEtaPtPlusPlus(0),
  57     fHistGeneratedEtaPtMinusMinus(0),
  58     fHistSurvivedEtaPtMinusMinus(0),
  59     fHistGeneratedEtaPtPlusMinus(0),
  60     fHistSurvivedEtaPtPlusMinus(0),
  61     fHistGeneratedPhiEtaPlusPlus(0),
  62     fHistSurvivedPhiEtaPlusPlus(0),
  63     fHistGeneratedPhiEtaMinusMinus(0),
  64     fHistSurvivedPhiEtaMinusMinus(0),
  65     fHistGeneratedPhiEtaPlusMinus(0),
  66     fHistSurvivedPhiEtaPlusMinus(0),
  67     fUseCentrality(kFALSE),
  68     fCentralityEstimator("V0M"),
  69     fCentralityPercentileMin(0.0),
  70     fCentralityPercentileMax(5.0),
  71     fVxMax(3.0),
  72     fVyMax(3.0),
  73     fVzMax(10.),
  74     fAODTrackCutBit(128),
  75     fMinNumberOfTPCClusters(80),
  76     fMaxChi2PerTPCCluster(4.0),
  77     fMaxDCAxy(3.0),
  78     fMaxDCAz(3.0),
  79     fMinPt(0.0),
  80     fMaxPt(20.0),
  81     fMinEta(-0.8),
  82     fMaxEta(0.8),
  83     fEtaRangeMin(0.0),
  84     fEtaRangeMax(1.6),
  85     fPtRangeMin(0.0),
  86     fPtRangeMax(20.0),
  87     fEtaBin(100), //=100 (BF) 16
  88     fdEtaBin(64), //=64 (BF)  16
  89     fPtBin(100), //=100 (BF)  36
  90     fHistSurvived4EtaPtPhiPlus(0),
  91     fHistSurvived8EtaPtPhiPlus(0)
  92
  93 {
  94   // Define input and output slots here
  95   // Input slot #0 works with a TChain
  96   DefineInput(0, TChain::Class());
  97   // Output slot #0 id reserved by the base class for AOD
  98   // Output slot #1 writes into a TH1 container
  99   DefineOutput(1, TList::Class());
 100   DefineOutput(2, TList::Class());
 101 }
 102
 103 //________________________________________________________________________
 104 void AliAnalysisTaskEffContBF::UserCreateOutputObjects() {
 105   // Create histograms
 106   // Called once
 107
 108   fQAList = new TList();
 109   fQAList->SetName("QAList");
 110   fQAList->SetOwner();
 111
 112   fOutputList = new TList();
 113   fOutputList->SetName("OutputList");
 114   fOutputList->SetOwner();
 115
 116   //Event stats.
 117   TString gCutName[4] = {"Total","Offline trigger",
 118                          "Vertex","Analyzed"};
 119   fHistEventStats = new TH1F("fHistEventStats",
 120                              "Event statistics;;N_{events}",
 121                              4,0.5,4.5);
 122   for(Int_t i = 1; i <= 4; i++)
 123     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 124   fQAList->Add(fHistEventStats);
 125
 126   //====================================================//
 127   Int_t ptBin = 36;
 128   Int_t etaBin = 16;
 129   Int_t phiBin = 100;
 130
 131   Double_t nArrayPt[37]={0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0};
 132   Double_t nArrayEta[17]={-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
 133
 134   Double_t nArrayPhi[phiBin+1];
 135   for(Int_t iBin = 0; iBin <= phiBin; iBin++)
 136     nArrayPhi[iBin] = iBin*TMath::TwoPi()/phiBin;
 137
 138   Int_t detaBin = 16;
 139   Int_t dphiBin = 100;
 140   Double_t nArrayDPhi[dphiBin+1];
 141   for(Int_t iBin = 0; iBin <= phiBin; iBin++)
 142     nArrayDPhi[iBin] = iBin*TMath::TwoPi()/dphiBin;
 143   Double_t nArrayDEta[17]={0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6};
 144   //====================================================//
 145
 146   //AOD analysis
 147   fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
 148                              1001,-0.5,100.5);
 149   fQAList->Add(fHistCentrality);
 150
 151   //multiplicity (good MC tracks)
 152   TString histName;
 153   histName = "fHistNMult";
 154   fHistNMult = new TH1F(histName.Data(),
 155                         ";N_{mult.}",
 156                         200,0,20000);
 157   fQAList->Add(fHistNMult);
 158
 159   //Vz addition+++++++++++++++++++++++++++++
 160   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
 161   fQAList->Add(fHistVz);
 162
 163   //Contamination for Secondaries
 164   fHistContaminationSecondariesPlus = new TH3D("fHistContaminationSecondariesPlus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 165   fOutputList->Add(fHistContaminationSecondariesPlus);
 166
 167   fHistContaminationSecondariesMinus = new TH3D("fHistContaminationSecondariesMinus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 168   fOutputList->Add(fHistContaminationSecondariesMinus);
 169
 170   //Contamination for Primaries
 171   fHistContaminationPrimariesPlus = new TH3D("fHistContaminationPrimariesPlus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 172   fOutputList->Add(fHistContaminationPrimariesPlus);
 173
 174   fHistContaminationPrimariesMinus = new TH3D("fHistContaminationPrimariesMinus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 175   fOutputList->Add(fHistContaminationPrimariesMinus);
 176
 177   //eta vs pt for MC positives
 178   fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
 179                                         "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 180                                         etaBin,nArrayEta, ptBin, nArrayPt,phiBin, nArrayPhi);
 181   // fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,nArrayPhi);
 182   fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
 183   fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
 184                                        "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 185                                        etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 186   fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
 187
 188   //eta vs pt for MC negatives
 189   fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
 190                                          "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 191                                          etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 192   fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
 193   fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
 194                                         "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 195                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 196   fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
 197
 198   //eta vs pt for MC positives (control)
 199   fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
 200                                             "Generated positive primaries;#eta;p_{T} (GeV/c)",
 201                                             etaBin,nArrayEta,ptBin,nArrayPt);
 202   fOutputList->Add(fHistGeneratedEtaPtPlusControl);
 203   fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
 204                                            "Survived positive primaries;#eta;p_{T} (GeV/c)",
 205                                            etaBin,nArrayEta,ptBin,nArrayPt);
 206   fOutputList->Add(fHistSurvivedEtaPtPlusControl);
 207
 208   //eta vs pt for MC negatives (control)
 209   fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
 210                                              "Generated positive primaries;#eta;p_{T} (GeV/c)",
 211                                              etaBin,nArrayEta,ptBin,nArrayPt);
 212   fOutputList->Add(fHistGeneratedEtaPtMinusControl);
 213   fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
 214                                             "Survived positive primaries;#eta;p_{T} (GeV/c)",
 215                                             etaBin,nArrayEta,ptBin,nArrayPt);
 216   fOutputList->Add(fHistSurvivedEtaPtMinusControl);
 217
 218   //eta vs pt for MC ++
 219   fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
 220                                          "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 221                                          detaBin,nArrayDEta,ptBin,nArrayPt);
 222   fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
 223   fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
 224                                         "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 225                                         detaBin,nArrayDEta,ptBin,nArrayPt);
 226   fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
 227
 228   //eta vs pt for MC --
 229   fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
 230                                            "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
 231                                            detaBin,nArrayDEta,ptBin,nArrayPt);
 232   fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
 233   fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
 234                                           "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
 235                                           detaBin,nArrayDEta,ptBin,nArrayPt);
 236   fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
 237
 238   //eta vs pt for MC +-
 239   fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
 240                                           "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
 241                                           detaBin,nArrayDEta,ptBin,nArrayPt);
 242   fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
 243   fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
 244                                          "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
 245                                          detaBin,nArrayDEta,ptBin,nArrayPt);
 246   fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
 247
 248   //=============================//
 249   //phi vs eta for MC ++
 250   fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
 251                                           "Generated ++ primaries;#Delta#phi",
 252                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 253   fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
 254   fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
 255                                          "Survived ++ primaries;#Delta#phi;#Delta#eta",
 256                                          dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 257   fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
 258
 259   //phi vs eta for MC --
 260   fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
 261                                             "Generated -- primaries;#Delta#phi;#Delta#eta",
 262                                             dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 263   fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
 264   fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
 265                                            "Survived -- primaries;#Delta#phi;#Delta#eta",
 266                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 267   fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
 268
 269   //phi vs eta for MC +-
 270   fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
 271                                            "Generated +- primaries;#Delta#phi;#Delta#eta",
 272                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 273   fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
 274   fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
 275                                           "Survived +- primaries;#Delta#phi;#Delta#eta",
 276                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 277   fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
 278
 279   fHistSurvived4EtaPtPhiPlus = new TH3F("fHistSurvived4EtaPtPhiPlus",
 280                                         "Survived4 + primaries;#eta;p_{T} (GeV/c);#phi",
 281                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 282   fOutputList->Add(fHistSurvived4EtaPtPhiPlus);
 283   fHistSurvived8EtaPtPhiPlus = new TH3F("fHistSurvived8EtaPtPhiPlus",
 284                                         "Survived8 + primaries;#eta;p_{T} (GeV/c);#phi",
 285                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 286   fOutputList->Add(fHistSurvived8EtaPtPhiPlus);
 287
 288   //fQAList->Print();
 289   //fOutputList->Print();
 290   PostData(1, fQAList);
 291   PostData(2, fOutputList);
 292 }
 293
 294 //________________________________________________________________________
 295 void AliAnalysisTaskEffContBF::UserExec(Option_t *) {
 296   // Main loop
 297   // Called for each event
 298
 299   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 300   if (!fAOD) {
 301     printf("ERROR: fAOD not available\n");
 302     return;
 303   }
 304
 305   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 306   if (!fArrayMC)
 307     AliFatal("No array of MC particles found !!!"); // MW  no AliFatal use return values
 308
 309   AliMCEvent* mcEvent = MCEvent();
 310   if (!mcEvent) {
 311     AliError("ERROR: Could not retrieve MC event");
 312     return;
 313   }
 314
 315   // arrays for 2 particle histograms
 316   Int_t nMCLabelCounter         = 0;
 317   const Int_t maxMCLabelCounter = 20000;
 318
 319   Double_t eta[maxMCLabelCounter];
 320   Double_t pt[maxMCLabelCounter];
 321   Double_t phi[maxMCLabelCounter];
 322   Int_t level[maxMCLabelCounter];
 323   Int_t charge[maxMCLabelCounter];
 324
 325   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fAOD->GetNumberOfTracks()));
 326
 327   fHistEventStats->Fill(1); //all events
 328
 329   //Centrality stuff
 330   AliAODHeader *header = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 331   Double_t nCentrality = 0;
 332
 333   if(fUseCentrality) {
 334     AliCentrality *centrality = header->GetCentralityP();
 335     nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 336
 337
 338     if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
 339                                              fCentralityPercentileMax,
 340                                              fCentralityEstimator.Data()))
 341       return;
 342     else {
 343       fHistEventStats->Fill(2); //triggered + centrality
 344       fHistCentrality->Fill(nCentrality);
 345     }
 346   }
 347   //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
 348
 349   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 350   if(vertex) {
 351     if(vertex->GetNContributors() > 0) {
 352       Double32_t fCov[6];
 353       vertex->GetCovarianceMatrix(fCov);
 354       if(fCov[5] != 0) {
 355         fHistEventStats->Fill(3); //events with a proper vertex
 356         if(TMath::Abs(vertex->GetX()) < fVxMax) {    // antes Xv
 357           //Printf("X Vertex: %lf", vertex->GetX());
 358           //Printf("Y Vertex: %lf", vertex->GetY());
 359           if(TMath::Abs(vertex->GetY()) < fVyMax) {  // antes Yv
 360             if(TMath::Abs(vertex->GetZ()) < fVzMax) {  // antes Zv
 361               //Printf("Z Vertex: %lf", vertex->GetZ());
 362
 363               fHistEventStats->Fill(4); //analyzed events
 364               fHistVz->Fill(vertex->GetZ());
 365
 366               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 367               Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
 368               Int_t nMCParticles = mcEvent->GetNumberOfTracks();
 369               TArrayI labelMCArray(nMCParticles);
 370
 371               for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
 372                 AliAODTrack* track = fAOD->GetTrack(jTracks);
 373                 if(!track) continue;
 374
 375                 if (!track->TestFilterBit(fAODTrackCutBit))
 376                   continue;
 377
 378                 //acceptance
 379                 if(TMath::Abs(track->Eta()) > fMaxEta)
 380                   continue;
 381                 if((track->Pt() > fMaxPt)||(track->Pt() <  fMinPt))
 382                   continue;
 383
 384                 Double_t phiRad = track->Phi();
 385
 386                 Int_t label = TMath::Abs(track->GetLabel());
 387                 if(label > nMCParticles) continue;
 388                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label);
 389                 Short_t gAODmcCharge = AODmcTrack->Charge();////
 390                 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
 391                 //if (!(AODmcTrack->IsPhysicalPrimary())) {
 392                 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
 393                 //}
 394                 if (AODmcTrack->IsPhysicalPrimary()) {
 395                   if(gAODmcCharge > 0){
 396                     fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 397                   }
 398                   if(gAODmcCharge < 0){
 399                     fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 400                   }
 401                 }
 402                 else{
 403                   if(gAODmcCharge > 0){
 404                     fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 405                   }
 406                   if(gAODmcCharge < 0){
 407                     fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 408                   }
 409                 }
 410               }//loop over tracks
 411               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 412
 413               //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
 414               for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 415                 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
 416                 if (!mcTrack) {
 417                   AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
 418                   continue;
 419                 }
 420
 421                 //exclude particles generated out of the acceptance
 422                 Double_t vz = mcTrack->Zv();
 423                 if (TMath::Abs(vz) > 50.) continue;
 424                 //acceptance
 425                 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
 426                   continue;
 427                 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
 428                   continue;
 429
 430                 if(!mcTrack->IsPhysicalPrimary()) continue;
 431
 432                 Short_t gMCCharge = mcTrack->Charge();
 433                 Double_t phiRad = mcTrack->Phi();
 434
 435                 if(gMCCharge > 0)
 436                   fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
 437                                                    mcTrack->Pt(),
 438                                                    phiRad);
 439                 else if(gMCCharge < 0)
 440                   fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
 441                                                     mcTrack->Pt(),
 442                                                     phiRad);
 443
 444                 Bool_t labelTPC = kTRUE;
 445                 if(labelTPC) {
 446                   labelMCArray.AddAt(iTracks,nMCLabelCounter);
 447                   if(nMCLabelCounter >= maxMCLabelCounter){
 448                     AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
 449                     break;
 450                   }
 451                   //fill the arrays for 2 particle analysis
 452                   eta[nMCLabelCounter]    = mcTrack->Eta();
 453                   pt[nMCLabelCounter]     = mcTrack->Pt();
 454                   phi[nMCLabelCounter]    = mcTrack->Phi();
 455                   charge[nMCLabelCounter] = gMCCharge;
 456
 457                   level[nMCLabelCounter]  = 1;
 458                   nMCLabelCounter += 1;
 459                 }
 460               }//loop over MC particles
 461
 462               fHistNMult->Fill(nMCLabelCounter);
 463
 464               //AOD track loop
 465               Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 466               TArrayI labelArray(nGoodTracks);
 467               Int_t labelCounter = 0;
 468
 469               for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 470                 AliAODTrack *trackAOD = static_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 471                 if(!trackAOD) continue;
 472
 473                 //track cuts
 474                 if (!trackAOD->TestFilterBit(fAODTrackCutBit))
 475                   continue;
 476
 477                 Int_t label = TMath::Abs(trackAOD->GetLabel());
 478                 if(IsLabelUsed(labelArray,label)) continue;
 479                 labelArray.AddAt(label,labelCounter);
 480                 labelCounter += 1;
 481
 482                 Bool_t iFound = kFALSE;
 483                 Int_t mcGoods = nMCLabelCounter;
 484                 for (Int_t k = 0; k < mcGoods; k++) {
 485                   Int_t mcLabel = labelMCArray.At(k);
 486                   iFound = kFALSE;
 487
 488                   if (mcLabel != TMath::Abs(label)) continue;
 489                   if(mcLabel != label) continue;
 490                   if(label > trackAOD->GetLabel()) continue;
 491
 492                   //acceptance
 493                   if(TMath::Abs(trackAOD->Eta()) > fMaxEta)
 494                     continue;
 495                   if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() <  fMinPt))
 496                     continue;
 497
 498                   Short_t gCharge = trackAOD->Charge();
 499                   Double_t phiRad = trackAOD->Phi();
 500
 501                   if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){
 502                     level[k]  = 2;
 503
 504                     if(gCharge > 0)
 505                       fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
 506                                                       trackAOD->Pt(),
 507                                                       phiRad);
 508                     else if(gCharge < 0)
 509                       fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
 510                                                        trackAOD->Pt(),
 511                                                        phiRad);
 512                   }//tracks
 513                 }//end of mcGoods
 514               }//AOD track loop
 515
 516               labelMCArray.Reset();
 517               labelArray.Reset();
 518
 519             }//Vz cut
 520           }//Vy cut
 521         }//Vx cut
 522       }//Vz resolution
 523     }//number of contributors
 524   }//valid vertex
 525
 526   // Here comes the 2 particle analysis
 527   // loop over all good MC particles
 528   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
 529     // control 1D histograms (charge might be different?)
 530     if(charge[i] > 0){
 531       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
 532       if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
 533     }
 534     else if(charge[i] < 0){
 535       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
 536       if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
 537     }
 538
 539
 540     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
 541
 542       if(charge[i] > 0 && charge[j] > 0 ){
 543         if(level[i] > 0 && level[j] > 0) {
 544           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 545           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 546             fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 547         }
 548         if(level[i] > 1 && level[j] > 1) {
 549           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 550           if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
 551             fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 552         }
 553       }
 554
 555       else if(charge[i] < 0 && charge[j] < 0 ){
 556         if(level[i] > 0 && level[j] > 0) {
 557           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 558           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 559             fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 560         }
 561         if(level[i] > 2 && level[j] > 1) {
 562           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 563           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 564             fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 565         }
 566       }
 567
 568       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
 569         if(level[i] > 0 && level[j] > 0) {
 570           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 571           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 572             fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 573         }
 574         if(level[i] > 2 && level[j] > 1) {
 575           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 576           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 577             fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 578         }
 579       }
 580     }
 581   }
 582 }
 583
 584 //________________________________________________________________________
 585 void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
 586   // Draw result to the screen
 587   // Called once at the end of the query
 588
 589   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 590   if (!fOutputList) {
 591     printf("ERROR: Output list not available\n");
 592     return;
 593   }
 594 }
 595
 596 //____________________________________________________________________//
 597 Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
 598   //Checks if the label is used already
 599   Bool_t status = kFALSE;
 600   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
 601     if(labelArray.At(i) == label)
 602       status = kTRUE;
 603   }
 604
 605   return status;
 606 }