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