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 <= dphiBin; 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       return;
 354     }
 355
 356     headers = header->GetNCocktailHeaders();
 357     eventHeader = header->GetCocktailHeader(0);
 358
 359
 360     if (!eventHeader)
 361       {
 362         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 363         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 364         AliError("First event header not found. Skipping this event.");
 365         return;
 366       }
 367
 368     skipParticlesAbove = eventHeader->NProduced();
 369     AliInfo(Form("Injected signals in this event (%d headers). Keeping particles/tracks of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 370   }
 371   // ==============================================================================================
 372
 373
 374   // arrays for 2 particle histograms
 375   Int_t nMCLabelCounter         = 0;
 376   const Int_t maxMCLabelCounter = 20000;
 377
 378   Double_t eta[maxMCLabelCounter];
 379   Double_t pt[maxMCLabelCounter];
 380   Double_t phi[maxMCLabelCounter];
 381   Int_t level[maxMCLabelCounter];
 382   Int_t charge[maxMCLabelCounter];
 383
 384   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fAOD->GetNumberOfTracks()));
 385
 386   fHistEventStats->Fill(1); //all events
 387
 388   //Centrality stuff
 389   Double_t nCentrality = 0;
 390   if(fUseCentrality) {
 391
 392     AliAODHeader *headerAOD = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 393     if (!headerAOD){
 394       AliFatal("AOD header found");
 395       return;
 396     }
 397
 398     AliCentrality *centrality = headerAOD->GetCentralityP();
 399     nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 400
 401
 402     if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
 403                                              fCentralityPercentileMax,
 404                                              fCentralityEstimator.Data()))
 405       return;
 406     else {
 407       fHistEventStats->Fill(2); //triggered + centrality
 408       fHistCentrality->Fill(nCentrality);
 409     }
 410   }
 411   //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
 412
 413   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 414   if(vertex) {
 415     if(vertex->GetNContributors() > 0) {
 416       Double32_t fCov[6];
 417       vertex->GetCovarianceMatrix(fCov);
 418       if(fCov[5] != 0) {
 419         fHistEventStats->Fill(3); //events with a proper vertex
 420         if(TMath::Abs(vertex->GetX()) < fVxMax) {    // antes Xv
 421           //Printf("X Vertex: %lf", vertex->GetX());
 422           //Printf("Y Vertex: %lf", vertex->GetY());
 423           if(TMath::Abs(vertex->GetY()) < fVyMax) {  // antes Yv
 424             if(TMath::Abs(vertex->GetZ()) < fVzMax) {  // antes Zv
 425               //Printf("Z Vertex: %lf", vertex->GetZ());
 426
 427               fHistEventStats->Fill(4); //analyzed events
 428               fHistVz->Fill(vertex->GetZ());
 429
 430               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 431               Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
 432               Int_t nMCParticles = mcEvent->GetNumberOfTracks();
 433               TArrayI labelMCArray(nMCParticles);
 434
 435               for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
 436                 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(jTracks));
 437                 if(!track) AliFatal("Not a standard AOD");
 438                 if(!track) continue;
 439
 440                 if (!track->TestFilterBit(fAODTrackCutBit))
 441                   continue;
 442
 443                 //acceptance
 444                 if(TMath::Abs(track->Eta()) > fMaxEta)
 445                   continue;
 446                 if((track->Pt() > fMaxPt)||(track->Pt() <  fMinPt))
 447                   continue;
 448
 449                 Double_t phiRad = track->Phi();
 450
 451                 Int_t label = TMath::Abs(track->GetLabel());
 452                 if(label > nMCParticles) continue;
 453                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label);
 454                 Short_t gAODmcCharge = AODmcTrack->Charge();////
 455                 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
 456                 //if (!(AODmcTrack->IsPhysicalPrimary())) {
 457                 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
 458                 //}
 459
 460                 // ==============================================================================================
 461                 // Partial copy from AliAnalyseLeadingTrackUE::RemoveInjectedSignals:
 462                 // Skip tracks that come from injected signals
 463                 if (fInjectedSignals)
 464                   {
 465
 466                     AliAODMCParticle* mother = AODmcTrack;
 467
 468                     // find the primary mother (if not already physical primary)
 469                     while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
 470                       {
 471                         if (((AliAODMCParticle*)mother)->GetMother() < 0)
 472                           {
 473                             mother = 0;
 474                             break;
 475                           }
 476
 477                         mother = (AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
 478                         if (!mother)
 479                           break;
 480                       }
 481
 482
 483                     if (!mother)
 484                       {
 485                         AliError(Form("WARNING: No mother found for particle %d:", AODmcTrack->GetLabel()));
 486                         continue;
 487                       }
 488
 489                     if (mother->GetLabel() >= skipParticlesAbove)
 490                       {
 491                         //AliInfo(Form("Remove particle %d (>= %d)",mother->GetLabel(),skipParticlesAbove));
 492                         continue;
 493                       }
 494                   }
 495                 // ==============================================================================================
 496
 497                 if (AODmcTrack->IsPhysicalPrimary()) {
 498                   if(gAODmcCharge > 0){
 499                     fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 500                   }
 501                   if(gAODmcCharge < 0){
 502                     fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 503                   }
 504                 }
 505                 else{
 506                   if(gAODmcCharge > 0){
 507                     fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
 508                   }
 509                   if(gAODmcCharge < 0){
 510                     fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
 511                   }
 512                 }
 513               }//loop over tracks
 514               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
 515
 516               //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
 517               for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 518                 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
 519                 if (!mcTrack) {
 520                   AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
 521                   continue;
 522                 }
 523
 524                 //exclude particles generated out of the acceptance
 525                 Double_t vz = mcTrack->Zv();
 526                 if (TMath::Abs(vz) > 50.) continue;
 527                 //acceptance
 528                 if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
 529                   continue;
 530                 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
 531                   continue;
 532
 533                 if(!mcTrack->IsPhysicalPrimary()) continue;
 534
 535                 Short_t gMCCharge = mcTrack->Charge();
 536                 Double_t phiRad = mcTrack->Phi();
 537
 538                 if(gMCCharge > 0)
 539                   fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
 540                                                    mcTrack->Pt(),
 541                                                    phiRad);
 542                 else if(gMCCharge < 0)
 543                   fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
 544                                                     mcTrack->Pt(),
 545                                                     phiRad);
 546
 547                 Bool_t labelTPC = kTRUE;
 548                 if(labelTPC) {
 549                   labelMCArray.AddAt(iTracks,nMCLabelCounter);
 550                   if(nMCLabelCounter >= maxMCLabelCounter){
 551                     AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
 552                     break;
 553                   }
 554                   //fill the arrays for 2 particle analysis
 555                   eta[nMCLabelCounter]    = mcTrack->Eta();
 556                   pt[nMCLabelCounter]     = mcTrack->Pt();
 557                   phi[nMCLabelCounter]    = mcTrack->Phi();
 558                   charge[nMCLabelCounter] = gMCCharge;
 559
 560                   level[nMCLabelCounter]  = 1;
 561                   nMCLabelCounter += 1;
 562                 }
 563               }//loop over MC particles
 564
 565               fHistNMult->Fill(nMCLabelCounter);
 566
 567               //AOD track loop
 568               Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 569               TArrayI labelArray(nGoodTracks);
 570               Int_t labelCounter = 0;
 571
 572               for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 573                 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 574                 if(!trackAOD) continue;
 575
 576                 //track cuts
 577                 if (!trackAOD->TestFilterBit(fAODTrackCutBit))
 578                   continue;
 579
 580                 Int_t label = TMath::Abs(trackAOD->GetLabel());
 581                 if(IsLabelUsed(labelArray,label)) continue;
 582                 labelArray.AddAt(label,labelCounter);
 583                 labelCounter += 1;
 584
 585                 Int_t mcGoods = nMCLabelCounter;
 586                 for (Int_t k = 0; k < mcGoods; k++) {
 587                   Int_t mcLabel = labelMCArray.At(k);
 588
 589                   if (mcLabel != TMath::Abs(label)) continue;
 590                   if(mcLabel != label) continue;
 591                   if(label > trackAOD->GetLabel()) continue;
 592
 593                   //acceptance
 594                   if(TMath::Abs(trackAOD->Eta()) > fMaxEta)
 595                     continue;
 596                   if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() <  fMinPt))
 597                     continue;
 598
 599                   Short_t gCharge = trackAOD->Charge();
 600                   Double_t phiRad = trackAOD->Phi();
 601
 602                   //===========================PID (so far only for electron rejection)===============================//
 603                   if(fElectronRejection) {
 604
 605                     // get the electron nsigma
 606                     Double_t nSigma = fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kElectron);
 607                     fHistNSigmaTPCvsPtbeforePID->Fill(trackAOD->Pt(),nSigma);
 608
 609                     // check only for given momentum range
 610                     if( trackAOD->Pt() > fElectronRejectionMinPt && trackAOD->Pt() < fElectronRejectionMaxPt ){
 611
 612                       //look only at electron nsigma
 613                       if(!fElectronOnlyRejection){
 614
 615                         //Make the decision based on the n-sigma of electrons
 616                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma) continue;
 617                       }
 618                       else{
 619
 620                         Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kPion));
 621                         Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kKaon));
 622                         Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kProton));
 623
 624                         //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
 625                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma
 626                            && nSigmaPions   > fElectronRejectionNSigma
 627                            && nSigmaKaons   > fElectronRejectionNSigma
 628                            && nSigmaProtons > fElectronRejectionNSigma ) continue;
 629                       }
 630                     }
 631
 632                     fHistNSigmaTPCvsPtafterPID->Fill(trackAOD->Pt(),nSigma);
 633
 634                   }
 635                   //===========================end of PID (so far only for electron rejection)===============================//
 636
 637                   if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){
 638                     level[k]  = 2;
 639
 640                     if(gCharge > 0)
 641                       fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
 642                                                       trackAOD->Pt(),
 643                                                       phiRad);
 644                     else if(gCharge < 0)
 645                       fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
 646                                                        trackAOD->Pt(),
 647                                                        phiRad);
 648                   }//tracks
 649                 }//end of mcGoods
 650               }//AOD track loop
 651
 652               labelMCArray.Reset();
 653               labelArray.Reset();
 654
 655             }//Vz cut
 656           }//Vy cut
 657         }//Vx cut
 658       }//Vz resolution
 659     }//number of contributors
 660   }//valid vertex
 661
 662   // Here comes the 2 particle analysis
 663   // loop over all good MC particles
 664   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
 665     // control 1D histograms (charge might be different?)
 666     if(charge[i] > 0){
 667       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
 668       if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
 669     }
 670     else if(charge[i] < 0){
 671       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
 672       if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
 673     }
 674
 675
 676     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
 677
 678       if(charge[i] > 0 && charge[j] > 0 ){
 679         if(level[i] > 0 && level[j] > 0) {
 680           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 681           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 682             fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 683         }
 684         if(level[i] > 1 && level[j] > 1) {
 685           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 686           if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
 687             fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 688         }
 689       }
 690
 691       else if(charge[i] < 0 && charge[j] < 0 ){
 692         if(level[i] > 0 && level[j] > 0) {
 693           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 694           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 695             fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 696         }
 697         if(level[i] > 2 && level[j] > 1) {
 698           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 699           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 700             fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 701         }
 702       }
 703
 704       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
 705         if(level[i] > 0 && level[j] > 0) {
 706           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 707           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 708             fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 709         }
 710         if(level[i] > 2 && level[j] > 1) {
 711           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 712           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
 713             fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 714         }
 715       }
 716     }
 717   }
 718 }
 719
 720 //________________________________________________________________________
 721 void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
 722   // Draw result to the screen
 723   // Called once at the end of the query
 724
 725   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 726   if (!fOutputList) {
 727     printf("ERROR: Output list not available\n");
 728     return;
 729   }
 730 }
 731
 732 //____________________________________________________________________//
 733 Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
 734   //Checks if the label is used already
 735   Bool_t status = kFALSE;
 736   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
 737     if(labelArray.At(i) == label)
 738       status = kTRUE;
 739   }
 740
 741   return status;
 742 }