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 "AliAODMCHeader.h"
  20 #include "AliCentrality.h"
  21 #include "AliGenEventHeader.h"
  22
  23 #include "AliLog.h"
  24 #include "AliAnalysisTaskEffContBF.h"
  25
  26 // ---------------------------------------------------------------------
  27 //
  28 // Task for calculating the efficiency of the Balance Function
  29 // for single particles and pairs
  30 //
  31 // ---------------------------------------------------------------------
  32
  33 ClassImp(AliAnalysisTaskEffContBF)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskEffContBF::AliAnalysisTaskEffContBF(const char *name)
  37   : AliAnalysisTaskSE(name),
  38     fAOD(0),
  39     fArrayMC(0),
  40     fQAList(0),
  41     fOutputList(0),
  42     fHistEventStats(0),
  43     fHistCentrality(0),
  44     fHistNMult(0),
  45     fHistVz(0),
  46     fHistNSigmaTPCvsPtbeforePID(0),
  47     fHistNSigmaTPCvsPtafterPID(0),
  48     fHistContaminationSecondariesPlus(0),
  49     fHistContaminationSecondariesMinus(0), //
  50     fHistContaminationPrimariesPlus(0),
  51     fHistContaminationPrimariesMinus(0), //
  52     fHistGeneratedEtaPtPhiPlus(0),
  53     fHistSurvivedEtaPtPhiPlus(0),
  54     fHistGeneratedEtaPtPhiMinus(0),
  55     fHistSurvivedEtaPtPhiMinus(0),
  56     fHistGeneratedEtaPtPlusControl(0),
  57     fHistSurvivedEtaPtPlusControl(0),
  58     fHistGeneratedEtaPtMinusControl(0),
  59     fHistSurvivedEtaPtMinusControl(0),
  60     fHistGeneratedEtaPtPlusPlus(0),
  61     fHistSurvivedEtaPtPlusPlus(0),
  62     fHistGeneratedEtaPtMinusMinus(0),
  63     fHistSurvivedEtaPtMinusMinus(0),
  64     fHistGeneratedEtaPtPlusMinus(0),
  65     fHistSurvivedEtaPtPlusMinus(0),
  66     fHistGeneratedPhiEtaPlusPlus(0),
  67     fHistSurvivedPhiEtaPlusPlus(0),
  68     fHistGeneratedPhiEtaMinusMinus(0),
  69     fHistSurvivedPhiEtaMinusMinus(0),
  70     fHistGeneratedPhiEtaPlusMinus(0),
  71     fHistSurvivedPhiEtaPlusMinus(0),
  72     fUseCentrality(kFALSE),
  73     fCentralityEstimator("V0M"),
  74     fCentralityPercentileMin(0.0),
  75     fCentralityPercentileMax(5.0),
  76     fInjectedSignals(kFALSE),
  77     fPIDResponse(0),
  78     fElectronRejection(kFALSE),
  79     fElectronOnlyRejection(kFALSE),
  80     fElectronRejectionNSigma(-1.),
  81     fElectronRejectionMinPt(0.),
  82     fElectronRejectionMaxPt(1000.),
  83     fVxMax(3.0),
  84     fVyMax(3.0),
  85     fVzMax(10.),
  86     fAODTrackCutBit(128),
  87     fMinNumberOfTPCClusters(80),
  88     fMaxChi2PerTPCCluster(4.0),
  89     fMaxDCAxy(3.0),
  90     fMaxDCAz(3.0),
  91     fMinPt(0.0),
  92     fMaxPt(20.0),
  93     fMinEta(-0.8),
  94     fMaxEta(0.8),
  95     fEtaRangeMin(0.0),
  96     fEtaRangeMax(1.6),
  97     fPtRangeMin(0.0),
  98     fPtRangeMax(20.0),
  99     fEtaBin(100), //=100 (BF) 16
 100     fdEtaBin(64), //=64 (BF)  16
 101     fPtBin(100), //=100 (BF)  36
 102     fHistSurvived4EtaPtPhiPlus(0),
 103     fHistSurvived8EtaPtPhiPlus(0)
 104
 105 {
 106   // Define input and output slots here
 107   // Input slot #0 works with a TChain
 108   DefineInput(0, TChain::Class());
 109   // Output slot #0 id reserved by the base class for AOD
 110   // Output slot #1 writes into a TH1 container
 111   DefineOutput(1, TList::Class());
 112   DefineOutput(2, TList::Class());
 113 }
 114
 115 //________________________________________________________________________
 116 void AliAnalysisTaskEffContBF::UserCreateOutputObjects() {
 117   // Create histograms
 118   // Called once
 119
 120   fQAList = new TList();
 121   fQAList->SetName("QAList");
 122   fQAList->SetOwner();
 123
 124   fOutputList = new TList();
 125   fOutputList->SetName("OutputList");
 126   fOutputList->SetOwner();
 127
 128   //Event stats.
 129   TString gCutName[4] = {"Total","Offline trigger",
 130                          "Vertex","Analyzed"};
 131   fHistEventStats = new TH1F("fHistEventStats",
 132                              "Event statistics;;N_{events}",
 133                              4,0.5,4.5);
 134   for(Int_t i = 1; i <= 4; i++)
 135     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 136   fQAList->Add(fHistEventStats);
 137
 138   //====================================================//
 139   Int_t ptBin = 36;
 140   Int_t etaBin = 16;
 141   Int_t phiBin = 100;
 142
 143   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};
 144   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};
 145
 146   Double_t nArrayPhi[phiBin+1];
 147   for(Int_t iBin = 0; iBin <= phiBin; iBin++)
 148     nArrayPhi[iBin] = iBin*TMath::TwoPi()/phiBin;
 149
 150   Int_t detaBin = 16;
 151   Int_t dphiBin = 100;
 152   Double_t nArrayDPhi[dphiBin+1];
 153   for(Int_t iBin = 0; iBin <= phiBin; iBin++)
 154     nArrayDPhi[iBin] = iBin*TMath::TwoPi()/dphiBin;
 155   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};
 156   //====================================================//
 157
 158   //AOD analysis
 159   fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
 160                              1001,-0.5,100.5);
 161   fQAList->Add(fHistCentrality);
 162
 163   //multiplicity (good MC tracks)
 164   TString histName;
 165   histName = "fHistNMult";
 166   fHistNMult = new TH1F(histName.Data(),
 167                         ";N_{mult.}",
 168                         200,0,20000);
 169   fQAList->Add(fHistNMult);
 170
 171   //Vz addition+++++++++++++++++++++++++++++
 172   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
 173   fQAList->Add(fHistVz);
 174
 175   //Electron cuts -> PID QA
 176   fHistNSigmaTPCvsPtbeforePID = new TH2F ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore",200, 0, 20, 200, -10, 10);
 177   fQAList->Add(fHistNSigmaTPCvsPtbeforePID);
 178
 179   fHistNSigmaTPCvsPtafterPID = new TH2F ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter",200, 0, 20, 200, -10, 10);
 180   fQAList->Add(fHistNSigmaTPCvsPtafterPID);
 181
 182   //Contamination for Secondaries
 183   fHistContaminationSecondariesPlus = new TH3D("fHistContaminationSecondariesPlus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 184   fOutputList->Add(fHistContaminationSecondariesPlus);
 185
 186   fHistContaminationSecondariesMinus = new TH3D("fHistContaminationSecondariesMinus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 187   fOutputList->Add(fHistContaminationSecondariesMinus);
 188
 189   //Contamination for Primaries
 190   fHistContaminationPrimariesPlus = new TH3D("fHistContaminationPrimariesPlus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 191   fOutputList->Add(fHistContaminationPrimariesPlus);
 192
 193   fHistContaminationPrimariesMinus = new TH3D("fHistContaminationPrimariesMinus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 194   fOutputList->Add(fHistContaminationPrimariesMinus);
 195
 196   //eta vs pt for MC positives
 197   fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
 198                                         "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 199                                         etaBin,nArrayEta, ptBin, nArrayPt,phiBin, nArrayPhi);
 200   // fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,nArrayPhi);
 201   fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
 202   fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
 203                                        "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 204                                        etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 205   fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
 206
 207   //eta vs pt for MC negatives
 208   fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
 209                                          "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 210                                          etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 211   fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
 212   fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
 213                                         "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 214                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 215   fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
 216
 217   //eta vs pt for MC positives (control)
 218   fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
 219                                             "Generated positive primaries;#eta;p_{T} (GeV/c)",
 220                                             etaBin,nArrayEta,ptBin,nArrayPt);
 221   fOutputList->Add(fHistGeneratedEtaPtPlusControl);
 222   fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
 223                                            "Survived positive primaries;#eta;p_{T} (GeV/c)",
 224                                            etaBin,nArrayEta,ptBin,nArrayPt);
 225   fOutputList->Add(fHistSurvivedEtaPtPlusControl);
 226
 227   //eta vs pt for MC negatives (control)
 228   fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
 229                                              "Generated positive primaries;#eta;p_{T} (GeV/c)",
 230                                              etaBin,nArrayEta,ptBin,nArrayPt);
 231   fOutputList->Add(fHistGeneratedEtaPtMinusControl);
 232   fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
 233                                             "Survived positive primaries;#eta;p_{T} (GeV/c)",
 234                                             etaBin,nArrayEta,ptBin,nArrayPt);
 235   fOutputList->Add(fHistSurvivedEtaPtMinusControl);
 236
 237   //eta vs pt for MC ++
 238   fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
 239                                          "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 240                                          detaBin,nArrayDEta,ptBin,nArrayPt);
 241   fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
 242   fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
 243                                         "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 244                                         detaBin,nArrayDEta,ptBin,nArrayPt);
 245   fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
 246
 247   //eta vs pt for MC --
 248   fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
 249                                            "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
 250                                            detaBin,nArrayDEta,ptBin,nArrayPt);
 251   fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
 252   fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
 253                                           "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
 254                                           detaBin,nArrayDEta,ptBin,nArrayPt);
 255   fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
 256
 257   //eta vs pt for MC +-
 258   fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
 259                                           "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
 260                                           detaBin,nArrayDEta,ptBin,nArrayPt);
 261   fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
 262   fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
 263                                          "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
 264                                          detaBin,nArrayDEta,ptBin,nArrayPt);
 265   fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
 266
 267   //=============================//
 268   //phi vs eta for MC ++
 269   fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
 270                                           "Generated ++ primaries;#Delta#phi",
 271                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 272   fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
 273   fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
 274                                          "Survived ++ primaries;#Delta#phi;#Delta#eta",
 275                                          dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 276   fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
 277
 278   //phi vs eta for MC --
 279   fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
 280                                             "Generated -- primaries;#Delta#phi;#Delta#eta",
 281                                             dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 282   fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
 283   fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
 284                                            "Survived -- primaries;#Delta#phi;#Delta#eta",
 285                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 286   fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
 287
 288   //phi vs eta for MC +-
 289   fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
 290                                            "Generated +- primaries;#Delta#phi;#Delta#eta",
 291                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 292   fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
 293   fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
 294                                           "Survived +- primaries;#Delta#phi;#Delta#eta",
 295                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
 296   fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
 297
 298   fHistSurvived4EtaPtPhiPlus = new TH3F("fHistSurvived4EtaPtPhiPlus",
 299                                         "Survived4 + primaries;#eta;p_{T} (GeV/c);#phi",
 300                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 301   fOutputList->Add(fHistSurvived4EtaPtPhiPlus);
 302   fHistSurvived8EtaPtPhiPlus = new TH3F("fHistSurvived8EtaPtPhiPlus",
 303                                         "Survived8 + primaries;#eta;p_{T} (GeV/c);#phi",
 304                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
 305   fOutputList->Add(fHistSurvived8EtaPtPhiPlus);
 306
 307   //fQAList->Print();
 308   //fOutputList->Print();
 309   PostData(1, fQAList);
 310   PostData(2, fOutputList);
 311 }
 312
 313 //________________________________________________________________________
 314 void AliAnalysisTaskEffContBF::UserExec(Option_t *) {
 315   // Main loop
 316   // Called for each event
 317
 318   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 319   if (!fAOD) {
 320     printf("ERROR: fAOD not available\n");
 321     return;
 322   }
 323
 324   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 325   if (!fArrayMC)
 326     AliFatal("No array of MC particles found !!!"); // MW  no AliFatal use return values
 327
 328   AliMCEvent* mcEvent = MCEvent();
 329   if (!mcEvent) {
 330     AliError("ERROR: Could not retrieve MC event");
 331     return;
 332   }
 333
 334   // PID Response task active?
 335   if(fElectronRejection) {
 336     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
 337     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
 338   }
 339
 340   // ==============================================================================================
 341   // Copy from AliAnalysisTaskPhiCorrelations:
 342   // For productions with injected signals, figure out above which label to skip particles/tracks
 343   Int_t skipParticlesAbove = 0;
 344   if (fInjectedSignals)
 345   {
 346     AliGenEventHeader* eventHeader = 0;
 347     Int_t headers = 0;
 348
 349     // AOD only
 350     AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 351     if (!header)
 352       AliFatal("fInjectedSignals set but no MC header found");
 353
 354     headers = header->GetNCocktailHeaders();
 355     eventHeader = header->GetCocktailHeader(0);
 356
 357
 358     if (!eventHeader)
 359       {
 360         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 361         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 362         AliError("First event header not found. Skipping this event.");
 363         return;
 364       }
 365
 366     skipParticlesAbove = eventHeader->NProduced();
 367     AliInfo(Form("Injected signals in this event (%d headers). Keeping particles/tracks of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 368   }
 369   // ==============================================================================================
 370
 371
 372   // arrays for 2 particle histograms
 373   Int_t nMCLabelCounter         = 0;
 374   const Int_t maxMCLabelCounter = 20000;
 375
 376   Double_t eta[maxMCLabelCounter];
 377   Double_t pt[maxMCLabelCounter];
 378   Double_t phi[maxMCLabelCounter];
 379   Int_t level[maxMCLabelCounter];
 380   Int_t charge[maxMCLabelCounter];
 381
 382   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fAOD->GetNumberOfTracks()));
 383
 384   fHistEventStats->Fill(1); //all events
 385
 386   //Centrality stuff
 387   AliAODHeader *header = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 388   Double_t nCentrality = 0;
 389
 390   if(fUseCentrality) {
 391     AliCentrality *centrality = header->GetCentralityP();
 392     nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 393
 394
 395     if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
 396                                              fCentralityPercentileMax,
 397                                              fCentralityEstimator.Data()))
 398       return;
 399     else {
 400       fHistEventStats->Fill(2); //triggered + centrality
 401       fHistCentrality->Fill(nCentrality);
 402     }
 403   }
 404   //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
 405
 406   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 407   if(vertex) {
 408     if(vertex->GetNContributors() > 0) {
 409       Double32_t fCov[6];
 410       vertex->GetCovarianceMatrix(fCov);
 411       if(fCov[5] != 0) {
 412         fHistEventStats->Fill(3); //events with a proper vertex
 413         if(TMath::Abs(vertex->GetX()) < fVxMax) {    // antes Xv
 414           //Printf("X Vertex: %lf", vertex->GetX());
 415           //Printf("Y Vertex: %lf", vertex->GetY());
 416           if(TMath::Abs(vertex->GetY()) < fVyMax) {  // antes Yv
 417             if(TMath::Abs(vertex->GetZ()) < fVzMax) {  // antes Zv
 418               //Printf("Z Vertex: %lf", vertex->GetZ());
 419
 420               fHistEventStats->Fill(4); //analyzed events
 421               fHistVz->Fill(vertex->GetZ());
 422
 423               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 424               Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
 425               Int_t nMCParticles = mcEvent->GetNumberOfTracks();
 426               TArrayI labelMCArray(nMCParticles);
 427
 428               for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
 429                 AliAODTrack* track = fAOD->GetTrack(jTracks);
 430                 if(!track) continue;
 431
 432                 if (!track->TestFilterBit(fAODTrackCutBit))
 433                   continue;
 434
 435                 //acceptance
 436                 if(TMath::Abs(track->Eta()) > fMaxEta)
 437                   continue;
 438                 if((track->Pt() > fMaxPt)||(track->Pt() <  fMinPt))
 439                   continue;
 440
 441                 Double_t phiRad = track->Phi();
 442
 443                 Int_t label = TMath::Abs(track->GetLabel());
 444                 if(label > nMCParticles) continue;
 445                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label);
 446                 Short_t gAODmcCharge = AODmcTrack->Charge();////
 447                 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
 448                 //if (!(AODmcTrack->IsPhysicalPrimary())) {
 449                 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
 450                 //}
 451
 452                 // ==============================================================================================
 453                 // Partial copy from AliAnalyseLeadingTrackUE::RemoveInjectedSignals:
 454                 // Skip tracks that come from injected signals
 455                 if (fInjectedSignals)
 456                   {
 457
 458                     AliAODMCParticle* mother = AODmcTrack;
 459
 460                     // find the primary mother (if not already physical primary)
 461                     while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 462                       {
 463                         if (((AliAODMCParticle*)mother)->GetMother() < 0)
 464                           {
 465                             mother = 0;
 466                             break;
 467                           }
 468
 469                         mother = (AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 470                         if (!mother)
 471                           break;
 472                       }
 473
 474
 475                     if (!mother)
 476                       {
 477                         AliError(Form("WARNING: No mother found for particle %d:", AODmcTrack->GetLabel()));
 478                         continue;
 479                       }
 480
 481                     if (mother->GetLabel() >= skipParticlesAbove)
 482                       {
 483                         //AliInfo(Form("Remove particle %d (>= %d)",mother->GetLabel(),skipParticlesAbove));
 484                         continue;
 485                       }
 486                   }
 487                 // ==============================================================================================
 488
 489                 if (AODmcTrack->IsPhysicalPrimary()) {
 490                   if(gAODmcCharge > 0){
 491                     fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 492                   }
 493                   if(gAODmcCharge < 0){
 494                     fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 495                   }
 496                 }
 497                 else{
 498                   if(gAODmcCharge > 0){
 499                     fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 500                   }
 501                   if(gAODmcCharge < 0){
 502                     fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 503                   }
 504                 }
 505               }//loop over tracks
 506               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 507
 508               //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
 509               for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 510                 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
 511                 if (!mcTrack) {
 512                   AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
 513                   continue;
 514                 }
 515
 516                 //exclude particles generated out of the acceptance
 517                 Double_t vz = mcTrack->Zv();
 518                 if (TMath::Abs(vz) > 50.) continue;
 519                 //acceptance
 520                 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
 521                   continue;
 522                 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
 523                   continue;
 524
 525                 if(!mcTrack->IsPhysicalPrimary()) continue;
 526
 527                 Short_t gMCCharge = mcTrack->Charge();
 528                 Double_t phiRad = mcTrack->Phi();
 529
 530                 if(gMCCharge > 0)
 531                   fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
 532                                                    mcTrack->Pt(),
 533                                                    phiRad);
 534                 else if(gMCCharge < 0)
 535                   fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
 536                                                     mcTrack->Pt(),
 537                                                     phiRad);
 538
 539                 Bool_t labelTPC = kTRUE;
 540                 if(labelTPC) {
 541                   labelMCArray.AddAt(iTracks,nMCLabelCounter);
 542                   if(nMCLabelCounter >= maxMCLabelCounter){
 543                     AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
 544                     break;
 545                   }
 546                   //fill the arrays for 2 particle analysis
 547                   eta[nMCLabelCounter]    = mcTrack->Eta();
 548                   pt[nMCLabelCounter]     = mcTrack->Pt();
 549                   phi[nMCLabelCounter]    = mcTrack->Phi();
 550                   charge[nMCLabelCounter] = gMCCharge;
 551
 552                   level[nMCLabelCounter]  = 1;
 553                   nMCLabelCounter += 1;
 554                 }
 555               }//loop over MC particles
 556
 557               fHistNMult->Fill(nMCLabelCounter);
 558
 559               //AOD track loop
 560               Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 561               TArrayI labelArray(nGoodTracks);
 562               Int_t labelCounter = 0;
 563
 564               for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 565                 AliAODTrack *trackAOD = static_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 566                 if(!trackAOD) continue;
 567
 568                 //track cuts
 569                 if (!trackAOD->TestFilterBit(fAODTrackCutBit))
 570                   continue;
 571
 572                 Int_t label = TMath::Abs(trackAOD->GetLabel());
 573                 if(IsLabelUsed(labelArray,label)) continue;
 574                 labelArray.AddAt(label,labelCounter);
 575                 labelCounter += 1;
 576
 577                 Bool_t iFound = kFALSE;
 578                 Int_t mcGoods = nMCLabelCounter;
 579                 for (Int_t k = 0; k < mcGoods; k++) {
 580                   Int_t mcLabel = labelMCArray.At(k);
 581                   iFound = kFALSE;
 582
 583                   if (mcLabel != TMath::Abs(label)) continue;
 584                   if(mcLabel != label) continue;
 585                   if(label > trackAOD->GetLabel()) continue;
 586
 587                   //acceptance
 588                   if(TMath::Abs(trackAOD->Eta()) > fMaxEta)
 589                     continue;
 590                   if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() <  fMinPt))
 591                     continue;
 592
 593                   Short_t gCharge = trackAOD->Charge();
 594                   Double_t phiRad = trackAOD->Phi();
 595
 596                   //===========================PID (so far only for electron rejection)===============================//
 597                   if(fElectronRejection) {
 598
 599                     // get the electron nsigma
 600                     Double_t nSigma = fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kElectron);
 601                     fHistNSigmaTPCvsPtbeforePID->Fill(trackAOD->Pt(),nSigma);
 602
 603                     // check only for given momentum range
 604                     if( trackAOD->Pt() > fElectronRejectionMinPt && trackAOD->Pt() < fElectronRejectionMaxPt ){
 605
 606                       //look only at electron nsigma
 607                       if(!fElectronOnlyRejection){
 608
 609                         //Make the decision based on the n-sigma of electrons
 610                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma) continue;
 611                       }
 612                       else{
 613
 614                         Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kPion));
 615                         Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kKaon));
 616                         Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kProton));
 617
 618                         //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
 619                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma
 620                            && nSigmaPions   > fElectronRejectionNSigma
 621                            && nSigmaKaons   > fElectronRejectionNSigma
 622                            && nSigmaProtons > fElectronRejectionNSigma ) continue;
 623                       }
 624                     }
 625
 626                     fHistNSigmaTPCvsPtafterPID->Fill(trackAOD->Pt(),nSigma);
 627
 628                   }
 629                   //===========================end of PID (so far only for electron rejection)===============================//
 630
 631                   if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){
 632                     level[k]  = 2;
 633
 634                     if(gCharge > 0)
 635                       fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
 636                                                       trackAOD->Pt(),
 637                                                       phiRad);
 638                     else if(gCharge < 0)
 639                       fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
 640                                                        trackAOD->Pt(),
 641                                                        phiRad);
 642                   }//tracks
 643                 }//end of mcGoods
 644               }//AOD track loop
 645
 646               labelMCArray.Reset();
 647               labelArray.Reset();
 648
 649             }//Vz cut
 650           }//Vy cut
 651         }//Vx cut
 652       }//Vz resolution
 653     }//number of contributors
 654   }//valid vertex
 655
 656   // Here comes the 2 particle analysis
 657   // loop over all good MC particles
 658   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
 659     // control 1D histograms (charge might be different?)
 660     if(charge[i] > 0){
 661       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
 662       if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
 663     }
 664     else if(charge[i] < 0){
 665       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
 666       if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
 667     }
 668
 669
 670     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
 671
 672       if(charge[i] > 0 && charge[j] > 0 ){
 673         if(level[i] > 0 && level[j] > 0) {
 674           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 675           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 676             fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 677         }
 678         if(level[i] > 1 && level[j] > 1) {
 679           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 680           if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
 681             fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 682         }
 683       }
 684
 685       else if(charge[i] < 0 && charge[j] < 0 ){
 686         if(level[i] > 0 && level[j] > 0) {
 687           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 688           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 689             fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 690         }
 691         if(level[i] > 2 && level[j] > 1) {
 692           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 693           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 694             fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 695         }
 696       }
 697
 698       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
 699         if(level[i] > 0 && level[j] > 0) {
 700           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 701           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 702             fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 703         }
 704         if(level[i] > 2 && level[j] > 1) {
 705           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 706           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 707             fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 708         }
 709       }
 710     }
 711   }
 712 }
 713
 714 //________________________________________________________________________
 715 void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
 716   // Draw result to the screen
 717   // Called once at the end of the query
 718
 719   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 720   if (!fOutputList) {
 721     printf("ERROR: Output list not available\n");
 722     return;
 723   }
 724 }
 725
 726 //____________________________________________________________________//
 727 Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
 728   //Checks if the label is used already
 729   Bool_t status = kFALSE;
 730   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
 731     if(labelArray.At(i) == label)
 732       status = kTRUE;
 733   }
 734
 735   return status;
 736 }