PWGLF/STRANGENESS/Correlations/AliAnalysisTaskV0ChCorrelations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2013, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Marek Bombara                                                  *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* The task selects candidates for K0s, Lambdas and AntiLambdas (trigger particles)
  17  * and calculates correlations with charged unidentified particles (associated particles) in phi and eta.
  18  * The task works with AOD (with or without MC info) events only and containes also mixing for acceptance corrections.
  19  * Last update edited by Marek Bombara, August 2013, Marek.Bombara@cern.ch
  20  */
  21
  22 #include <TCanvas.h>
  23 #include <TList.h>
  24 #include <TH1.h>
  25 #include <TH2.h>
  26 #include <TH3.h>
  27 #include <THnSparse.h>
  28 #include <TObjArray.h>
  29 #include <TObject.h>
  30
  31 #include "AliLog.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliAODTrack.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODv0.h"
  36 #include "AliAODcascade.h"
  37 #include "AliAODVertex.h"
  38
  39 #include "AliMCEvent.h"
  40 #include "AliMCVertex.h"
  41 #include "AliGenEventHeader.h"
  42 #include "AliAODMCHeader.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliAnalyseLeadingTrackUE.h"
  45
  46 #include "AliAODPid.h"
  47 #include "AliPIDResponse.h"
  48 #include "AliEventPoolManager.h"
  49 #include "AliCentrality.h"
  50
  51 #include "AliAnalysisTaskV0ChCorrelations.h"
  52 #include "AliPhysicsSelectionTask.h"
  53
  54 #include <AliMultiInputEventHandler.h>
  55 #include <AliMixInputEventHandler.h>
  56
  57 ClassImp(AliAnalysisTaskV0ChCorrelations)
  58 ClassImp(AliV0ChBasicParticle)
  59
  60 //________________________________________________________________________
  61 AliAnalysisTaskV0ChCorrelations::AliAnalysisTaskV0ChCorrelations(const char *name) // All data members should be initialised here
  62    : AliAnalysisTaskSE(name),
  63      fAnalysisMC(kFALSE),
  64          fFillMixed(kTRUE),
  65          fMixingTracks(50000),
  66          fPoolMgr(0x0),
  67
  68      fDcaDToPV(0.5),
  69          fDcaV0D(0.1),
  70          fWithChCh(kFALSE),
  71          fOStatus(1),
  72
  73          fOutput(0),
  74          fPIDResponse(0),
  75      fHistCentVtx(0),
  76      fHistMultiMain(0),
  77          fHistMassK0(0),
  78          fHistMassLambda(0),
  79          fHistMassAntiLambda(0),
  80
  81          fHistdPhidEtaSib(0),
  82          fHistdPhidEtaMix(0),
  83          fHistNTrigSib(0),
  84
  85          fHistMCPtCentTrig(0),
  86      fHistRCPtCentTrig(0),
  87      fHistMCPtCentAs(0),
  88      fHistRCPtCentAs(0),
  89      fHistRCPtCentAll(0),
  90
  91          fHistTemp(0),
  92          fHistTemp2(0)// The last in the above list should not have a comma after it
  93 {
  94    // Constructor
  95    // Define input and output slots here (never in the dummy constructor)
  96    // Input slot #0 works with a TChain - it is connected to the default input container
  97    // Output slot #1 writes into a TH1 container
  98    DefineOutput(1, TList::Class());                                            // for output list
  99 }
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskV0ChCorrelations::~AliAnalysisTaskV0ChCorrelations()
 103 {
 104    // Destructor. Clean-up the output list, but not the histograms that are put inside
 105    // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 106    if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 107       delete fOutput;
 108    }
 109 }
 110
 111 //________________________________________________________________________
 112 void AliAnalysisTaskV0ChCorrelations::UserCreateOutputObjects()
 113 {
 114    // Create histograms
 115
 116    fOutput = new TList();
 117    fOutput->SetOwner();  // IMPORTANT!
 118    // defining bins for centrality
 119    Int_t nCentralityBins  = 9;
 120    Double_t centBins[] = {0.,10.,20.,30.,40.,50.,60.,70.,80.,90.};
 121    const Double_t* centralityBins = centBins;
 122    // defining bins for Z vertex
 123    Int_t nZvtxBins  = 7;
 124    //Double_t vertexBins[] = {-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.};
 125    Double_t vertexBins[] = {-7.,-5.,-3.,-1.,1.,3.,5.,7.};
 126    const Double_t* zvtxBins = vertexBins;
 127    // pt bins of associated particles for the analysis
 128    Int_t nPtBins = 7;
 129    const Double_t PtBins[8] = {2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0};
 130    //Int_t nPtBins = 1;
 131    //const Double_t PtBins[2] = {3.0,15.0};
 132    // pt bins of trigger particles for the analysis
 133    Int_t nPtBinsV0 = 11;
 134    //const Double_t PtBinsV0[2] = {6.0,15.0};
 135    const Double_t PtBinsV0[12] = {4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0};
 136    // V0 candidate: 1 - K0sig, 2 - Lamsig, 3 - Alamsig, 4 - K0bg, 5 - Lambg, 6 - Alambg
 137    Int_t nTrigC = 7;
 138    const Double_t TrigC[8] = {0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5};
 139
 140    // Create general histograms
 141    fHistCentVtx = new TH2F("fHistCentVtx", "Centrality vs. Z vertex", nCentralityBins, centralityBins, nZvtxBins, zvtxBins);
 142    fHistMultiMain = new TH1F("fHistMultiMain", "Multiplicity of main events", 2000, 0, 2000);
 143
 144    const Int_t mrBins[3] = {nPtBinsV0, nCentralityBins, nTrigC};
 145    const Double_t mrMin[3] = {PtBinsV0[0], centralityBins[0], TrigC[0]};
 146    const Double_t mrMax[3] = {PtBinsV0[11], centralityBins[9], TrigC[6]};
 147
 148    // Create histograms for reconstruction track and V0 efficiency
 149    fHistMCPtCentTrig = new THnSparseF("fHistMCPtCentTrig", "MC Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);
 150    fHistRCPtCentTrig = new THnSparseF("fHistRCPtCentTrig", "Rec Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);
 151
 152    // pt bins of associated particles for the efficiency
 153    Int_t nPtBinsAs = 13;
 154    const Double_t PtBinsAs[14] = {2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0};
 155
 156    fHistMCPtCentAs = new TH2D("fHistMCPtCentAs", "MC Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
 157    fHistRCPtCentAs = new TH2D("fHistRCPtCentAs", "Rec Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
 158    fHistRCPtCentAll = new TH2D("fHistRCPtCentAll", "Rec Pt vs. Cent. All Prim+Sec", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
 159
 160    // defining bins for mass distributions
 161    Int_t nBins = 200;
 162    Double_t mk0Min = 0.40;
 163    Double_t mk0Max = 0.60;
 164    Double_t mlaMin = 1.07;
 165    Double_t mlaMax = 1.17;
 166    Double_t malMin = 1.07;
 167    Double_t malMax = 1.17;
 168
 169    const Int_t spBins[3] = {nBins, nPtBinsV0, nCentralityBins};
 170    const Double_t spMinK0[3] = {mk0Min, PtBinsV0[0], centralityBins[0]};
 171    const Double_t spMaxK0[3] = {mk0Max, PtBinsV0[11], centralityBins[9]};
 172    const Double_t spMinLa[3] = {mlaMin, PtBinsV0[0], centralityBins[0]};
 173    const Double_t spMaxLa[3] = {mlaMax, PtBinsV0[11], centralityBins[9]};
 174    const Double_t spMinAl[3] = {malMin, PtBinsV0[0], centralityBins[0]};
 175    const Double_t spMaxAl[3] = {malMax, PtBinsV0[11], centralityBins[9]};
 176    // Create mass histograms
 177    fHistMassK0 = new THnSparseF("fHistMassK0","V0 mass for K0 hypothesis", 3, spBins, spMinK0, spMaxK0);
 178    fHistMassLambda = new THnSparseF("fHistMassLambda","V0 mass for Lambda hypothesis", 3, spBins, spMinLa, spMaxLa);
 179    fHistMassAntiLambda = new THnSparseF("fHistMassAntiLambda","V0 mass for AntiLambda hypothesis", 3, spBins, spMinAl, spMaxAl);
 180
 181    // defining bins for dPhi distributions
 182    const Int_t nbPhiBins = 72;
 183    //const Double_t kPi = TMath::Pi();
 184    //Double_t PhiMin = -kPi/2.;
 185    //Double_t PhiMax = -kPi/2. + 2*kPi;
 186    Double_t PhiMin = -1.57;
 187    Double_t PhiMax = 4.71;
 188    Double_t PhiBins[nbPhiBins+1] = {0.};
 189    PhiBins[0] = PhiMin;
 190    for (Int_t i=0; i<nbPhiBins; i++) { PhiBins[i+1] = PhiBins[i] + (PhiMax - PhiMin)/nbPhiBins; }
 191
 192    // defining bins for dEta distributions
 193    const Int_t nbEtaBins = 40;
 194    Double_t EtaMin = -2.0;
 195    Double_t EtaMax = 2.0;
 196    Double_t EtaBins[nbEtaBins+1] = {0.};
 197    EtaBins[0] = EtaMin;
 198    for (Int_t i=0; i<nbEtaBins; i++) { EtaBins[i+1] = EtaBins[i] + (EtaMax - EtaMin)/nbEtaBins; }
 199
 200    const Int_t corBins[7] = {nbPhiBins, nbEtaBins, nPtBinsV0, nPtBins, nCentralityBins, nZvtxBins, nTrigC};
 201    const Double_t corMin[7] = {PhiBins[0], EtaBins[0], PtBinsV0[0], PtBins[0], centralityBins[0], zvtxBins[0], TrigC[0]};
 202    const Double_t corMax[7] = {PhiBins[72], EtaBins[40], PtBinsV0[11], PtBins[7], centralityBins[9], zvtxBins[7], TrigC[7]};
 203    // Create correlation histograms
 204    fHistdPhidEtaSib = new THnSparseF("fHistdPhidEtaSib","dPhi vs. dEta siblings", 7, corBins, corMin, corMax);
 205    fHistdPhidEtaMix = new THnSparseF("fHistdPhidEtaMix","dPhi vs. dEta mixed", 7, corBins, corMin, corMax);
 206
 207    // Create histograms for counting the numbers of trigger particles
 208    const Int_t corNTrigBins[5] = {nPtBinsV0, nCentralityBins, nZvtxBins, nTrigC};
 209    const Double_t corNTrigMin[5] = {PtBinsV0[0], centBins[0], vertexBins[0], TrigC[0]};
 210    const Double_t corNTrigMax[5] = {PtBinsV0[11], centBins[9], vertexBins[7], TrigC[7]};
 211    fHistNTrigSib = new THnSparseF("fHistNTrigSib","Number trigger sib", 4, corNTrigBins, corNTrigMin, corNTrigMax);
 212
 213    // Histograms for debugging
 214    fHistTemp = new TH1D("fHistTemp", "Temporary", 100, -10., 10.);
 215    fHistTemp2 = new TH1D("fHistTemp2", "Temporary2", 100, -10., 10.);
 216
 217    fOutput->Add(fHistCentVtx);
 218
 219    fOutput->Add(fHistMultiMain);
 220    fOutput->Add(fHistMassK0);
 221    fOutput->Add(fHistMassLambda);
 222    fOutput->Add(fHistMassAntiLambda);
 223
 224    fOutput->Add(fHistdPhidEtaSib);
 225    fOutput->Add(fHistdPhidEtaMix);
 226    fOutput->Add(fHistNTrigSib);
 227
 228    fOutput->Add(fHistMCPtCentTrig);
 229    fOutput->Add(fHistRCPtCentTrig);
 230    fOutput->Add(fHistMCPtCentAs);
 231    fOutput->Add(fHistRCPtCentAs);
 232    fOutput->Add(fHistRCPtCentAll);
 233
 234    fOutput->Add(fHistTemp);
 235    fOutput->Add(fHistTemp2);
 236
 237    PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 238
 239    // Settings for event mixing -------------------------------------
 240    Int_t trackDepth = fMixingTracks;
 241    Int_t poolSize   = 200;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 242
 243    fPoolMgr = new AliEventPoolManager(poolSize, trackDepth, nCentralityBins, centBins, nZvtxBins, vertexBins);
 244    //----------------------------------------------
 245 }
 246
 247 //________________________________________________________________________
 248 void AliAnalysisTaskV0ChCorrelations::Terminate(Option_t *)
 249 {
 250    // Draw result to screen, or perform fitting, normalizations
 251    // Called once at the end of the query
 252
 253    fOutput = dynamic_cast<TList*>(GetOutputData(1));
 254    if (!fOutput) { AliError("Could not retrieve TList fOutput"); return; }
 255
 256    // NEW HISTO should be retrieved from the TList container in the above way,
 257    // so it is available to draw on a canvas such as below
 258 }
 259
 260 //_________________________________________________________________________
 261 void AliAnalysisTaskV0ChCorrelations::UserExec(Option_t *)
 262 {
 263     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 264     AliInputEventHandler *inEvMain = (AliInputEventHandler*)mgr->GetInputEventHandler();
 265
 266     // physics selection
 267     UInt_t maskIsSelected = inEvMain->IsEventSelected();
 268    // 2010 data trigger selection
 269     //Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
 270    // 2011 data trigger selection
 271     Bool_t isSelected = ((maskIsSelected & AliVEvent::kMB) || (maskIsSelected & AliVEvent::kCentral) || (maskIsSelected & AliVEvent::kSemiCentral));
 272     if (!isSelected) return;
 273         // calculating the event types
 274     if (maskIsSelected & AliVEvent::kMB) fHistTemp->Fill(2);
 275     if (maskIsSelected & AliVEvent::kCentral) fHistTemp->Fill(4);
 276     if (maskIsSelected & AliVEvent::kSemiCentral) fHistTemp->Fill(6);
 277
 278     AliAODEvent* aod = (AliAODEvent*)inEvMain->GetEvent();
 279     fPIDResponse = inEvMain->GetPIDResponse();
 280
 281   // pt intervals for trigger particles
 282   const Double_t kPi = TMath::Pi();
 283   Double_t PtTrigMin = 4.0;
 284   Double_t PtTrigMax = 15.0;
 285   // pt interval for associated particles
 286   Double_t PtAssocMin = 2.0;
 287   fHistMultiMain->Fill(aod->GetNumberOfTracks());
 288
 289   // Vertex cut
 290   Double_t cutPrimVertex = 7.0;
 291   AliAODVertex *myPrimVertex = aod->GetPrimaryVertex();
 292   if (!myPrimVertex) return;
 293   if ( ( TMath::Abs(myPrimVertex->GetZ()) ) >= cutPrimVertex) return ;
 294   Double_t lPVx = myPrimVertex->GetX();
 295   Double_t lPVy = myPrimVertex->GetY();
 296   Double_t lPVz = myPrimVertex->GetZ();
 297
 298   if (TMath::Abs(lPVx)<10e-5 && TMath::Abs(lPVy)<10e-5 && TMath::Abs(lPVz)<10e-5) return;
 299
 300   // Centrality definition
 301   Double_t lCent = 0.0;
 302   AliCentrality *centralityObj = 0;
 303   centralityObj = aod->GetHeader()->GetCentralityP();
 304   lCent = centralityObj->GetCentralityPercentile("V0M");
 305   if ((lCent < 0.)||(lCent > 90.)) return;
 306   fHistCentVtx->Fill(lCent,lPVz);
 307
 308 //=========== MC loop ===============================
 309   if (fAnalysisMC)
 310   {
 311     AliAODMCHeader *aodMCheader = (AliAODMCHeader*)aod->FindListObject(AliAODMCHeader::StdBranchName());
 312     Float_t vzMC = aodMCheader->GetVtxZ();
 313     if (TMath::Abs(vzMC) >= 7.) return;
 314     //retrieve MC particles from event
 315     TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
 316     if(!mcArray){
 317         Printf("No MC particle branch found");
 318         return;
 319     }
 320
 321     Int_t nMCAllTracks = mcArray->GetEntriesFast();
 322     // new tracks array - without injected signal
 323     TObjArray * mcTracks = new TObjArray;
 324     //selectedMCTracks->SetOwner(kTRUE);
 325
 326     for (Int_t i = 0; i < nMCAllTracks; i++)
 327     {
 328         AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcArray->At(i);
 329         if (!mcTrack) {
 330             Error("ReadEventAODMC", "Could not receive particle %d", i);
 331             continue;
 332         }
 333         mcTracks->Add(mcTrack);
 334     }
 335
 336     // get labels for injected particles -------
 337     TObject* mc = mcArray;
 338     //TObjArray * mcTracks = mcArray;
 339     Int_t skipParticlesAbove = 0;
 340     //AliGenEventHeader* eventHeader = 0;
 341     //Int_t headers = 0;
 342     //headers = aodMCheader->GetNCocktailHeaders();
 343     //eventHeader = aodMCheader->GetCocktailHeader(0);
 344     AliGenEventHeader* eventHeader = aodMCheader->GetCocktailHeader(0);
 345     Int_t headers = aodMCheader->GetNCocktailHeaders();
 346     if (!eventHeader)
 347     {
 348       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 349       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 350       AliError("First event header not found. Skipping this event.");
 351       return;
 352     }
 353     skipParticlesAbove = eventHeader->NProduced();
 354     //cout << "ski label = " << skipParticlesAbove << endl;
 355     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s.Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 356     //cout << "before MC compressing = " << mcTracks->GetEntriesFast() << endl;
 357     RemovingInjectedSignal(mcTracks,mc,skipParticlesAbove);
 358     //cout << "after MC compressing = " << mcTracks->GetEntriesFast() << endl;
 359     // -------------
 360
 361     Int_t nMCTracks = mcTracks->GetEntriesFast();
 362     //cout << "number of MC tracks = " << nMCTracks << endl;
 363     for (Int_t iMC = 0; iMC<nMCTracks; iMC++)
 364     {
 365       AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcTracks->At(iMC);
 366         if (!mcTrack) {
 367             Error("ReadEventAODMC", "Could not receive particle %d", iMC);
 368             continue;
 369       }
 370       // track part
 371       Double_t mcTrackEta = mcTrack->Eta();
 372       Double_t mcTrackPt = mcTrack->Pt();
 373       Bool_t TrIsPrim = mcTrack->IsPhysicalPrimary();
 374       Bool_t TrEtaMax = TMath::Abs(mcTrackEta)<0.9;
 375       Bool_t TrPtMin = mcTrackPt>PtAssocMin;
 376       Bool_t TrCharge = (mcTrack->Charge())!=0;
 377       if (TrIsPrim && TrEtaMax && TrPtMin && TrCharge) fHistMCPtCentAs->Fill(mcTrackPt,lCent);
 378       // V0 part
 379       Int_t mcMotherPdg = 0;
 380           Int_t mcPartPdg = mcTrack->GetPdgCode();
 381
 382       // Keep only K0s, Lambda and AntiLambda
 383           if ((mcPartPdg != 310) && (mcPartPdg != 3122) && (mcPartPdg != (-3122))) continue;
 384
 385       //cout << " mc pdg is " << mcPartPdg << endl;
 386
 387       Bool_t IsK0 = mcPartPdg==310;
 388       Bool_t IsLambda = mcPartPdg==3122;
 389       Bool_t IsAntiLambda = mcPartPdg==-3122;
 390           Bool_t IsSigma = kFALSE;
 391           Int_t mcMotherLabel = mcTrack->GetMother();
 392       AliAODMCParticle *mcMother = (AliAODMCParticle*)mcArray->At(mcMotherLabel);
 393           if (mcMotherLabel < 0) {mcMotherPdg = 0;} else {mcMotherPdg = mcMother->GetPdgCode();}
 394           //if ((mcMotherLabel >= 0) && mcMother)
 395           //if ((mcMotherLabel >= 0) && mcMother)
 396           //{
 397         //  Bool_t IsFromCascade = ((mcMotherPdg==3312)||(mcMotherPdg==-3312)||(mcMotherPdg==3334)||(mcMotherPdg==-3334));
 398           Bool_t IsFromSigma = ((mcMotherPdg==3212)||(mcMotherPdg==-3212));
 399           IsFromSigma = IsFromSigma || ((mcMotherPdg==3224)||(mcMotherPdg==-3224));
 400           IsFromSigma = IsFromSigma || ((mcMotherPdg==3214)||(mcMotherPdg==-3214));
 401           IsFromSigma = IsFromSigma || ((mcMotherPdg==3114)||(mcMotherPdg==-3114));
 402                   if ((IsFromSigma) && (mcMother->IsPhysicalPrimary()) && (IsLambda || IsAntiLambda)) IsSigma = kTRUE;
 403           Double_t mcRapidity = mcTrack->Y();
 404           Bool_t V0RapMax = TMath::Abs(mcRapidity)<0.75;
 405           Bool_t PtInterval = ((mcTrackPt>PtTrigMin)&&(mcTrackPt<PtTrigMax));
 406           //Bool_t PtInterval = kTRUE;
 407           IsK0 = IsK0 && (mcTrack->IsPhysicalPrimary());
 408           IsLambda = IsLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
 409           IsAntiLambda = IsAntiLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
 410           Double_t mcK0[3] = {mcTrackPt, lCent, 1};
 411           Double_t mcLa[3] = {mcTrackPt, lCent, 2};
 412           Double_t mcAl[3] = {mcTrackPt, lCent, 3};
 413           if (IsK0 && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcK0);
 414           if (IsLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcLa);
 415           if (IsAntiLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcAl);
 416            //}
 417     }
 418     // ------- access the real data -----------
 419     Int_t nTracks = aod->GetNumberOfTracks();
 420     // new tracks array - without injected signal
 421     TObjArray * selectedMCTracks = new TObjArray;
 422     //selectedMCTracks->SetOwner(kTRUE);
 423
 424     for (Int_t i = 0; i < nTracks; i++)
 425     {
 426         AliAODTrack* tr = aod->GetTrack(i);
 427         selectedMCTracks->Add(tr);
 428     }
 429    // cout << "before compressing = " << selectedMCTracks->GetEntriesFast() << endl;
 430     RemovingInjectedSignal(selectedMCTracks,mc,skipParticlesAbove);
 431     //AliAnalyseLeadingTrackUE::RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
 432     //RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
 433    // cout << "after compressing = " << selectedMCTracks->GetEntriesFast() << endl;
 434     //-----------------------
 435     //cout << "before getting label -4" << endl;
 436     Int_t nRecTracks = selectedMCTracks->GetEntriesFast();
 437     //cout << "before getting label -3" << endl;
 438     for (Int_t i = 0; i < nRecTracks; i++)
 439     {
 440       //AliAODTrack* tras = aod->GetTrack(i);
 441       AliAODTrack* tras = (AliAODTrack*)selectedMCTracks->At(i);
 442       //cout << "before getting label -2" << endl;
 443       //cout << " and i = " << i << endl;
 444       //cout << "pt = " << tras->Pt() << " and i = " << i << endl;
 445       if ((tras->Pt())<PtAssocMin) continue;
 446       //cout << "before getting label -1" << endl;
 447       if (!(IsMyGoodPrimaryTrack(tras))) continue;
 448       //cout << "before getting label 0" << endl;
 449       Int_t AssocLabel = tras->GetLabel();
 450       if (AssocLabel<=0) continue;
 451       //cout << "before getting label 1" << endl;
 452       Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
 453       //Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
 454       //if(!(static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary())) continue;
 455              //cout << "before getting label 2" << endl;
 456     Double_t mcPt = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Pt();
 457              //cout << "before getting label 3" << endl;
 458     //Double_t mcEta = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Eta();
 459           //if (mcPt<PtAssocMin) continue;
 460           //if (TMath::Abs(mcEta)>0.8) continue;
 461       //Double_t trPt = tras->Pt();
 462       fHistRCPtCentAll->Fill(mcPt,lCent);
 463       if (isPhyPrim) fHistRCPtCentAs->Fill(mcPt,lCent);
 464     }
 465     // ------- end of real data access, for V0s see the main V0 loop --------
 466   }
 467 //============= End of MC loop ======================
 468
 469         // Track selection loop
 470         //--------------------------------
 471         Int_t nTracks = aod->GetNumberOfTracks();
 472         // new tracks array
 473         TObjArray * selectedTracks = new TObjArray;
 474         selectedTracks->SetOwner(kTRUE);
 475
 476         Bool_t ChChWith = GetWithChCh();
 477     TObjArray * selectedChargedTriggers = new TObjArray;
 478     selectedChargedTriggers->SetOwner(kTRUE);
 479     for (Int_t i = 0; i < nTracks; i++)
 480     {
 481         AliAODTrack* tr = aod->GetTrack(i);
 482         if ((tr->Pt())<PtAssocMin) continue;
 483         if (!(IsMyGoodPrimaryTrack(tr))) continue;
 484         selectedTracks->Add(tr);
 485         // saving the Charged trigger particles
 486         if ((tr->Pt()>4.)&&(tr->Pt()<15.))
 487         {
 488            selectedChargedTriggers->Add(new AliV0ChBasicParticle(tr->Eta(), tr->Phi(), tr->Pt(), 7));
 489         }
 490     }
 491
 492         //---------------------------------
 493
 494         // V0 loop for reconstructed event
 495     TObjArray * selectedV0s = new TObjArray;
 496         selectedV0s->SetOwner(kTRUE);
 497
 498         Int_t nV0s = aod->GetNumberOfV0s();
 499
 500         for (Int_t i = 0; i < nV0s; i++)
 501         { // start of V0 slection loop
 502                 AliAODv0* aodV0 = dynamic_cast<AliAODv0 *>(aod->GetV0(i));
 503                 if (!aodV0) {
 504          AliError(Form("ERROR: Could not retrieve aodv0 %d", i));
 505          continue;
 506         }
 507                 if (aodV0->GetOnFlyStatus()) fHistTemp->Fill(6.);
 508                 if (!aodV0->GetOnFlyStatus()) fHistTemp->Fill(8.);
 509                 //cout << "pt of v0: " << aodV0->Pt() << endl;
 510                 if (((aodV0->Pt())<PtTrigMin)||((aodV0->Pt())>PtTrigMax)) continue;
 511         // get daughters
 512         const AliAODTrack *myTrackPos;
 513         const AliAODTrack *myTrackNeg;
 514             AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
 515         AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));
 516
 517         if (!myTrackPosTest || !myTrackNegTest) {
 518            Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
 519            continue;
 520         }
 521
 522         if( myTrackPosTest->Charge() ==1){
 523            myTrackPos = myTrackPosTest;
 524            myTrackNeg = myTrackNegTest;
 525         }
 526
 527         if( myTrackPosTest->Charge() ==-1){
 528             myTrackPos = myTrackNegTest;
 529             myTrackNeg = myTrackPosTest;
 530         }
 531
 532         // effective mass calculations for each hypothesis
 533                 Double_t lInvMassK0 = aodV0->MassK0Short();
 534         Double_t lInvMassAntiLambda = aodV0->MassAntiLambda();
 535         Double_t lInvMassLambda = aodV0->MassLambda();
 536
 537                 // calculations for c*tau cut--------------------------------------
 538         //      Double_t lDVx = aodV0->GetSecVtxX();
 539     //    Double_t lDVy = aodV0->GetSecVtxY();
 540     //    Double_t lDVz = aodV0->GetSecVtxZ();
 541     //    const Double_t kLambdaMass = 1.115683;
 542     //    const Double_t kK0Mass = 0.497648;
 543     //    Double_t cutcTauLam = 3*7.89;
 544     //    Double_t cutcTauK0 = 3*2.68;
 545      //   Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lDVx - lPVx,2) + TMath::Power(lDVy- lPVy,2) + TMath::Power(lDVz - lPVz,2 ));
 546       //  Double_t lPV0 = TMath::Sqrt((aodV0->Pt())*(aodV0->Pt())+(aodV0->Pz())*(aodV0->Pz()));
 547      //   Double_t lcTauLam = (lV0DecayLength*kLambdaMass)/lPV0;
 548      //   Double_t lcTauK0 = (lV0DecayLength*kK0Mass)/lPV0;
 549         // sc - standard cuts
 550                 //Bool_t cutK0sc = (lcTauK0<cutcTauK0);
 551         //Bool_t cutLambdasc = (lcTauLam<cutcTauLam);
 552         //Bool_t cutAntiLambdasc = (lcTauLam<cutcTauLam);
 553                 Bool_t cutK0sc = kTRUE;
 554         Bool_t cutLambdasc = kTRUE;
 555         Bool_t cutAntiLambdasc = kTRUE;
 556
 557                 //------------------------------------------------
 558
 559                 // preparations for daughter's PID cut------------------------------
 560         Float_t nSigmaPosPion   = 0.;
 561         Float_t nSigmaNegPion   = 0.;
 562         Float_t nSigmaPosProton = 0.;
 563         Float_t nSigmaNegProton = 0.;
 564
 565         const AliAODPid *pPid = myTrackPos->GetDetPid();
 566         const AliAODPid *nPid = myTrackNeg->GetDetPid();
 567
 568         if (pPid)
 569         {
 570             Double_t pdMom = pPid->GetTPCmomentum();
 571             if (pdMom<1.)
 572             {
 573                 nSigmaPosPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
 574                 nSigmaPosProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
 575             }
 576         }
 577
 578         if (nPid)
 579         {
 580             Double_t ndMom = nPid->GetTPCmomentum();
 581             if (ndMom<1.)
 582             {
 583                 nSigmaNegPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
 584                 nSigmaNegProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
 585             }
 586         }
 587                 Bool_t bpPion = TMath::Abs(nSigmaPosPion) <= 3.;
 588         Bool_t bpProton = TMath::Abs(nSigmaPosProton) <= 3.;
 589         Bool_t bnPion = TMath::Abs(nSigmaNegPion) <= 3.;
 590         Bool_t bnProton = TMath::Abs(nSigmaNegProton) <= 3.;
 591
 592         Bool_t cutK0Pid = (bpPion && bnPion);
 593         Bool_t cutLambdaPid = (bpProton && bnPion);
 594         Bool_t cutAntiLambdaPid = (bpPion && bnProton);
 595         //--------------------------------------------------
 596         // mass cuts
 597         Bool_t cutMassLambda = ((lInvMassLambda>1.10) && (lInvMassLambda<1.13));
 598         Bool_t cutMassAntiLambda = ((lInvMassAntiLambda>1.10) && (lInvMassAntiLambda<1.13));
 599         Bool_t cutMassK0 = (lInvMassK0>0.47) && (lInvMassK0<0.53);
 600
 601         cutK0sc = cutK0Pid && (!cutMassLambda) && (!cutMassAntiLambda);
 602                 cutLambdasc = cutLambdaPid && (!cutMassK0);
 603                 cutAntiLambdasc = cutAntiLambdaPid && (!cutMassK0);
 604         // special cut related to AP diagram for K0s
 605                 Bool_t k0APcut = (aodV0->PtArmV0()>(TMath::Abs(0.2*aodV0->AlphaV0())));
 606                 cutK0sc = cutK0sc && k0APcut;
 607         // fill the mass histograms
 608
 609         Int_t oStatus = GetOStatus();
 610
 611         if (!IsMyGoodV0(aod,aodV0,myTrackPos,myTrackNeg,oStatus)) continue;
 612                 Double_t spK0[3] = {lInvMassK0, aodV0->Pt(), lCent};
 613                 Double_t spLa[3] = {lInvMassLambda, aodV0->Pt(), lCent};
 614                 Double_t spAl[3] = {lInvMassAntiLambda, aodV0->Pt(), lCent};
 615         if (cutK0sc) fHistMassK0->Fill(spK0);
 616                 if (cutLambdasc) fHistMassLambda->Fill(spLa);
 617                 if (cutAntiLambdasc) fHistMassAntiLambda->Fill(spAl);
 618         // select final V0s for correlation, selected background to study its contribution to correlation function
 619         // the values for signal might change in the future ...
 620         Bool_t K0Signal = (lInvMassK0>0.48)&&(lInvMassK0<0.52);
 621         Bool_t K0Bckg = ((lInvMassK0>0.40)&&(lInvMassK0<0.44)) || ((lInvMassK0>0.56)&&(lInvMassK0<0.60));
 622
 623         Bool_t LamSignal = (lInvMassLambda>1.108)&&(lInvMassLambda<1.125);
 624         Bool_t LamBckg = ((lInvMassLambda>1.090)&&(lInvMassLambda<1.100)) || ((lInvMassLambda>1.135)&&(lInvMassLambda<1.145));
 625
 626         Bool_t ALamSignal = (lInvMassAntiLambda>1.108)&&(lInvMassAntiLambda<1.125);
 627         Bool_t ALamBckg = ((lInvMassAntiLambda>1.090)&&(lInvMassAntiLambda<1.100)) || ((lInvMassAntiLambda>1.135)&&(lInvMassAntiLambda<1.145));
 628
 629                 // Fill selected V0 particle array
 630         if ((cutK0sc)&&(K0Signal))
 631         {
 632             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 1.));
 633             Double_t nTrigK0Sig[4] = {aodV0->Pt(), lCent, lPVz, 1.};
 634             fHistNTrigSib->Fill(nTrigK0Sig);
 635         }
 636         if ((cutK0sc)&&(K0Bckg))
 637         {
 638             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 4.));
 639             Double_t nTrigK0Bkg[4] = {aodV0->Pt(), lCent, lPVz, 4.};
 640             fHistNTrigSib->Fill(nTrigK0Bkg);
 641         }
 642         if ((cutLambdasc)&&(LamSignal))
 643         {
 644             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 2.));
 645             Double_t nTrigLaSig[4] = {aodV0->Pt(), lCent, lPVz, 2.};
 646             fHistNTrigSib->Fill(nTrigLaSig);
 647         }
 648         if ((cutLambdasc)&&(LamBckg))
 649         {
 650             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 5.));
 651             Double_t nTrigLaBkg[4] = {aodV0->Pt(), lCent, lPVz, 5.};
 652             fHistNTrigSib->Fill(nTrigLaBkg);
 653         }
 654         if ((cutAntiLambdasc)&&(ALamSignal))
 655         {
 656             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 3.));
 657             Double_t nTrigAlSig[4] = {aodV0->Pt(), lCent, lPVz, 3.};
 658             fHistNTrigSib->Fill(nTrigAlSig);
 659         }
 660         if ((cutAntiLambdasc)&&(ALamBckg))
 661         {
 662             selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 6.));
 663             Double_t nTrigAlBkg[4] = {aodV0->Pt(), lCent, lPVz, 6.};
 664             fHistNTrigSib->Fill(nTrigAlBkg);
 665         }
 666
 667
 668         //===== MC part for V0 reconstruction efficiency ==============
 669         if (fAnalysisMC)
 670         {
 671                 TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
 672                 if(!mcArray){
 673                 Printf("No MC particle branch found");
 674                 return;
 675                 }
 676
 677                         Int_t MotherOfMotherPdg = 0;
 678
 679                 Int_t myTrackPosLabel = TMath::Abs(myTrackPos->GetLabel());
 680                 Int_t myTrackNegLabel = TMath::Abs(myTrackNeg->GetLabel());
 681                 AliAODMCParticle *mcPosTrack = (AliAODMCParticle*)mcArray->At(myTrackPosLabel);
 682                         if (!mcPosTrack) continue;
 683                 AliAODMCParticle *mcNegTrack = (AliAODMCParticle*)mcArray->At(myTrackNegLabel);
 684                         if (!mcNegTrack) continue;
 685
 686                 Int_t PosTrackPdg = mcPosTrack->GetPdgCode();
 687                 Int_t NegTrackPdg = mcNegTrack->GetPdgCode();
 688                 //if (!mcPosTrack->IsPrimary()) continue;
 689                 //if (!mcNegTrack->IsPrimary()) continue;
 690
 691                 Int_t myTrackPosMotherLabel = mcPosTrack->GetMother();
 692                 Int_t myTrackNegMotherLabel = mcNegTrack->GetMother();
 693
 694                 if ((myTrackPosMotherLabel==-1)||(myTrackNegMotherLabel==-1)) continue;
 695
 696                 AliAODMCParticle *mcPosMother = (AliAODMCParticle*)mcArray->At(myTrackPosMotherLabel);
 697                         if (!mcPosMother) continue;
 698                 Int_t MotherPdg = mcPosMother->GetPdgCode();
 699                         Int_t MotherOfMother = mcPosMother->GetMother();
 700                         //if (MotherOfMother == -1) MotherOfMotherPdg = 0;
 701
 702                         if (myTrackPosMotherLabel!=myTrackNegMotherLabel) continue;
 703                 //if (!mcPosMother->IsPrimary()) continue;
 704                         Bool_t IsK0FromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==310)&&(PosTrackPdg==211)&&(NegTrackPdg==-211);
 705                         Bool_t IsLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
 706                         Bool_t IsAntiLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);
 707
 708                         Bool_t ComeFromSigma = kFALSE;
 709                         Bool_t ComeFromSigmaLa = kFALSE;
 710                         Bool_t ComeFromSigmaAl = kFALSE;
 711
 712                     if (MotherOfMother != -1)
 713                         {
 714                         AliAODMCParticle *mcPosMotherOfMother = (AliAODMCParticle*)mcArray->At(MotherOfMother);
 715                                 MotherOfMotherPdg = mcPosMotherOfMother->GetPdgCode();
 716                 Int_t MoMPdg = TMath::Abs(MotherOfMotherPdg);
 717                                 ComeFromSigma = (mcPosMotherOfMother->IsPhysicalPrimary())&&((MoMPdg==3212)||(MoMPdg==3224)||(MoMPdg==3214)||(MoMPdg==3114));
 718                                 ComeFromSigmaLa = ComeFromSigma && (MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
 719                                 ComeFromSigmaAl = ComeFromSigma && (MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);
 720                         }
 721
 722             IsLambdaFromMC = IsLambdaFromMC || ComeFromSigmaLa;
 723             IsAntiLambdaFromMC = IsAntiLambdaFromMC || ComeFromSigmaAl;
 724
 725                 Double_t RecMotherPt = aodV0->Pt();
 726                         //cout << "Pt of rec v0 = " << RecMotherPt << " nMC = " << mcArray->GetEntries() << endl;
 727                         //cout << "Pos Label = " << myTrackPosLabel << " Neg Label = " << myTrackNegLabel << endl;
 728                 Double_t rcK0[3] = {RecMotherPt, lCent, 1};
 729                 Double_t rcLa[3] = {RecMotherPt, lCent, 2};
 730                 Double_t rcAl[3] = {RecMotherPt, lCent, 3};
 731                 if ((cutK0sc)&&(K0Signal)&&(IsK0FromMC)) fHistRCPtCentTrig->Fill(rcK0);
 732                 if ((cutLambdasc)&&(LamSignal)&&(IsLambdaFromMC)) fHistRCPtCentTrig->Fill(rcLa);
 733                 if ((cutAntiLambdasc)&&(ALamSignal)&&(IsAntiLambdaFromMC)) fHistRCPtCentTrig->Fill(rcAl);
 734
 735         }
 736
 737         //===== End of the MC part for V0 reconstruction efficiency ===
 738
 739                 Int_t nSelectedTracks = selectedTracks->GetEntries();
 740                 // Correlation part
 741                 //===================================
 742                 for (Int_t j = 0; j < nSelectedTracks; j++)
 743                 {
 744                         AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
 745                         if ((atr->Pt())>=(aodV0->Pt())) continue;
 746                         Double_t dEta = atr->Eta() - aodV0->Eta();
 747                         Double_t dPhi = atr->Phi() - aodV0->Phi();
 748                         if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
 749                         if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;
 750                         // removing autocorrelations
 751                         //----------------------------------
 752                         Int_t negID = myTrackNeg->GetID();
 753                         Int_t posID = myTrackPos->GetID();
 754                         Int_t atrID = atr->GetID();
 755
 756                         if ((TMath::Abs(negID)+1)==(TMath::Abs(atrID))) continue;
 757                         if ((TMath::Abs(posID)+1)==(TMath::Abs(atrID))) continue;
 758                         //----------------------------------
 759
 760                         fHistNTrigSib->Sumw2();
 761                         fHistdPhidEtaSib->Sumw2();
 762                         // Filling correlation histograms and histograms for triggers counting
 763                         //----------------- K0 ---------------------
 764                         if ((cutK0sc)&&(K0Signal))
 765                         {
 766                                 Double_t spK0Sig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 1.};
 767                                 fHistdPhidEtaSib->Fill(spK0Sig);
 768                         }
 769
 770                         if ((cutK0sc)&&(K0Bckg))
 771                         {
 772                                 Double_t spK0Bkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 4.};
 773                                 fHistdPhidEtaSib->Fill(spK0Bkg);
 774                         }
 775                         //---------------- Lambda -------------------
 776                         if ((cutLambdasc)&&(LamSignal))
 777                         {
 778                                 Double_t spLaSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 2.};
 779                                 fHistdPhidEtaSib->Fill(spLaSig);
 780                         }
 781
 782                         if ((cutLambdasc)&&(LamBckg))
 783                         {
 784                                 Double_t spLaBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 5.};
 785                                 fHistdPhidEtaSib->Fill(spLaBkg);
 786                         }
 787                         //------------- AntiLambda -------------------
 788                         if ((cutAntiLambdasc)&&(ALamSignal))
 789                         {
 790                                 Double_t spAlSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 3.};
 791                                 fHistdPhidEtaSib->Fill(spAlSig);
 792                         }
 793
 794                         if ((cutAntiLambdasc)&&(ALamBckg))
 795                         {
 796                                 Double_t spAlBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 6.};
 797                                 fHistdPhidEtaSib->Fill(spAlBkg);
 798                         }
 799
 800                 } // end of correlation loop
 801                 //===================================
 802
 803
 804         } // end of V0 selection loop
 805
 806         // ===========================================
 807     // Ch-Ch correlation part
 808
 809     if (ChChWith)
 810     {
 811         Int_t nSelectedChargedTriggers = selectedChargedTriggers->GetEntries();
 812         for (Int_t i = 0; i < nSelectedChargedTriggers; i++)
 813         {
 814             AliV0ChBasicParticle* chTrig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
 815             Double_t chTrigPhi = chTrig->Phi();
 816             Double_t chTrigEta = chTrig->Eta();
 817             Double_t chTrigPt = chTrig->Pt();
 818
 819             Double_t nTrigChSig[4] = {chTrigPt, lCent, lPVz, 7.};
 820             fHistNTrigSib->Fill(nTrigChSig);
 821
 822             Int_t nSelectedTracks = selectedTracks->GetEntries();
 823             for (Int_t j = 0; j < nSelectedTracks; j++)
 824             {
 825                 AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
 826                 if ((atr->Pt())>=chTrigPt) continue;
 827                 Double_t dEta = atr->Eta() - chTrigEta;
 828                 Double_t dPhi = atr->Phi() - chTrigPhi;
 829                 if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
 830                 if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;
 831
 832                 Double_t spCh[7] = {dPhi, dEta, chTrigPt, atr->Pt(), lCent, lPVz, 7.};
 833                 fHistdPhidEtaSib->Fill(spCh);
 834             }
 835         }
 836     }
 837     // end of Ch-Ch correlation part
 838
 839         // Mixing ==============================================
 840
 841         fHistdPhidEtaMix->Sumw2();
 842     AliEventPool* pool = fPoolMgr->GetEventPool(lCent, lPVz);
 843         if (!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", lCent, lPVz));
 844         //pool->SetDebug(1);
 845     if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
 846         {
 847                 Int_t nMix = pool->GetCurrentNEvents();
 848                 for (Int_t jMix=0; jMix<nMix; jMix++)
 849                 {// loop through mixing events
 850                         TObjArray* bgTracks = pool->GetEvent(jMix);
 851                         for (Int_t i=0; i<selectedV0s->GetEntriesFast(); i++)
 852                         {// loop through selected V0 particles
 853                                 AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedV0s->At(i);
 854                                 Double_t trigPhi = trig->Phi();
 855                                 Double_t trigEta = trig->Eta();
 856                                 Double_t trigPt = trig->Pt();
 857                                 Short_t trigC = trig->WhichCandidate();
 858                                 for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
 859                                 { // mixing tracks loop
 860                                         AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
 861                                         // be careful tracks may have bigger pt than v0s.
 862                                         if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;
 863                                         Double_t dEtaMix = assoc->Eta() - trigEta;
 864                                         Double_t dPhiMix = assoc->Phi() - trigPhi;
 865                                         if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
 866                                         if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;
 867
 868                                         Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
 869                                     fHistdPhidEtaMix->Fill(spMix);
 870
 871                                 } // end of mixing track loop
 872                         }// end of loop through selected V0 particles
 873                         if (ChChWith)
 874             {
 875                 for (Int_t i=0; i<selectedChargedTriggers->GetEntriesFast(); i++)
 876                 {// loop through selected charged trigger particles
 877                     AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
 878                     Double_t trigPhi = trig->Phi();
 879                     Double_t trigEta = trig->Eta();
 880                     Double_t trigPt = trig->Pt();
 881                     Short_t trigC = trig->WhichCandidate();
 882                     for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
 883                     { // mixing tracks loop
 884                         AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
 885                         // be careful tracks may have bigger pt than v0s.
 886                         if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;
 887
 888                         Double_t dEtaMix = assoc->Eta() - trigEta;
 889                         Double_t dPhiMix = assoc->Phi() - trigPhi;
 890                         if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
 891                         if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;
 892
 893                         Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
 894                         fHistdPhidEtaMix->Fill(spMix);
 895
 896                     } // end of mixing track loop
 897                 }// end of loop through selected Ch particles
 898             } // end of ChCh mixing
 899
 900                 }// end of loop of mixing events
 901         }
 902
 903         TObjArray* tracksClone = (TObjArray*) selectedTracks->Clone();
 904         tracksClone->SetOwner(kTRUE);
 905         pool->UpdatePool(tracksClone);
 906         //delete selectedtracks;
 907
 908     PostData(1, fOutput);
 909
 910 }
 911 //___________________________________________
 912 Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodPrimaryTrack(const AliAODTrack *t)
 913 {
 914         // Pseudorapidity cut
 915         if (TMath::Abs(t->Eta())>0.9) return kFALSE;
 916         // Should correspond to set of cuts suitable for correlation analysis, 768 - hybrid tracks for 2011 data
 917     if (!t->TestFilterBit(768)) return kFALSE;
 918
 919         return kTRUE;
 920 }
 921 //_____________________________________________
 922 Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodDaughterTrack(const AliAODTrack *t)
 923 {
 924         // TPC refit
 925         if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
 926         // Minimum number of clusters
 927         Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
 928         if (nCrossedRowsTPC < 70) return kFALSE;
 929         Int_t findable=t->GetTPCNclsF();
 930         if (findable <= 0) return kFALSE;
 931         if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
 932
 933         return kTRUE;
 934 }
 935 //______________________________________________
 936 Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodV0(const AliAODEvent* aod, const AliAODv0* aodV0, const AliAODTrack* myTrackPos, const AliAODTrack* myTrackNeg, Int_t oSta)
 937 {
 938         if (!aodV0) {
 939        AliError(Form("ERROR: Could not retrieve aodV0"));
 940        return kFALSE;
 941     }
 942
 943         // Offline reconstructed V0 only
 944     if (oSta==1) {if (aodV0->GetOnFlyStatus()) return kFALSE;}
 945     if (oSta==3) {if (!aodV0->GetOnFlyStatus()) return kFALSE;}
 946
 947         if (oSta==2)
 948         {
 949                 if (aodV0->GetOnFlyStatus())
 950                 {
 951                         return kTRUE;
 952                 } else {
 953                         return kFALSE;
 954                 }
 955         }
 956     if (oSta==4)
 957         {
 958                 if (!aodV0->GetOnFlyStatus())
 959                 {
 960                         return kTRUE;
 961                 } else {
 962                         return kFALSE;
 963                 }
 964         }
 965
 966     // Get daughters and check them
 967         AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
 968         AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));
 969
 970         if (!myTrackPosTest || !myTrackNegTest) {
 971                 Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
 972                 return kFALSE;
 973         }
 974
 975     if( myTrackPosTest->Charge() ==1){
 976             myTrackPos = myTrackPosTest;
 977             myTrackNeg = myTrackNegTest;
 978     }
 979
 980     if( myTrackPosTest->Charge() ==-1){
 981             myTrackPos = myTrackNegTest;
 982             myTrackNeg = myTrackPosTest;
 983     }
 984
 985         // Track cuts for daughter tracks
 986     if ( !(IsMyGoodDaughterTrack(myTrackPos)) || !(IsMyGoodDaughterTrack(myTrackNeg)) ) return kFALSE;
 987
 988         // Unlike signs of daughters
 989     if (myTrackNegTest->Charge() == myTrackPosTest->Charge()) return kFALSE;
 990
 991         // Rapidity cut
 992         //Double_t lCutRap = 0.75;
 993         //Double_t lRapK0s = aodV0->Y(310);
 994         //Double_t lRapLambda = aodV0->Y(3122);
 995         //Double_t lRapAntiLambda = aodV0->Y(-3122);
 996
 997         // Pseudorapidity cut - there are high pt V0
 998         Double_t lCutEta = 0.75;
 999         Double_t lEtaV0 = aodV0->Eta();
1000         if (TMath::Abs(lEtaV0)>=lCutEta) return kFALSE;
1001         //if (TMath::Abs(lRapK0s)>=lCutRap) return kFALSE;
1002         //if (TMath::Abs(lRapLambda)>=lCutRap) return kFALSE;
1003         //if (TMath::Abs(lRapAntiLambda)>=lCutRap) return kFALSE;
1004
1005     // getting global variables
1006         Float_t dcaDaughtersToPrimVtx = GetDcaDToPV();
1007         Float_t dcaBetweenDaughters = GetDcaV0D();
1008
1009     // DCA of daughter track to Primary Vertex
1010     Float_t xyn=aodV0->DcaNegToPrimVertex();
1011     if (TMath::Abs(xyn)<dcaDaughtersToPrimVtx) return kFALSE;
1012     Float_t xyp=aodV0->DcaPosToPrimVertex();
1013     if (TMath::Abs(xyp)<dcaDaughtersToPrimVtx) return kFALSE;
1014
1015         // DCA of daughter tracks
1016     Double_t dca=aodV0->DcaV0Daughters();
1017     if (dca>dcaBetweenDaughters) return kFALSE;
1018
1019         // Cosinus of pointing angle
1020     Double_t cpa=aodV0->CosPointingAngle(aod->GetPrimaryVertex());
1021     if (cpa<0.998) return kFALSE;
1022
1023         // Fiducial volume cut
1024     Double_t xyz[3]; aodV0->GetSecondaryVtx(xyz);
1025     Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
1026     if (r2<0.9*0.9) return kFALSE;
1027     if (r2>100*100) return kFALSE;
1028
1029         // c*tau cut - in main V0 loop - depends on particle hypothesis
1030
1031         // Minimum pt of daughters
1032     Double_t  lMomPos[3] = {999,999,999};
1033         Double_t  lMomNeg[3] = {999,999,999};
1034
1035         lMomPos[0] = aodV0->MomPosX();
1036     lMomPos[1] = aodV0->MomPosY();
1037     lMomPos[2] = aodV0->MomPosZ();
1038
1039     lMomNeg[0] = aodV0->MomNegX();
1040     lMomNeg[1] = aodV0->MomNegY();
1041     lMomNeg[2] = aodV0->MomNegZ();
1042
1043     Double_t lPtPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1]);
1044     Double_t lPtNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1]);
1045
1046         Double_t cutMinPtDaughter = 0.160;
1047         if (lPtPos<cutMinPtDaughter || lPtNeg<cutMinPtDaughter) return kFALSE;
1048
1049         // Daughter PID cut - in main V0 loop - depends on particle hypothesis
1050
1051         return kTRUE;
1052 }
1053
1054 void AliAnalysisTaskV0ChCorrelations::RemovingInjectedSignal(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
1055 {
1056   // remove injected signals (primaries above <maxLabel>)
1057   // <tracks> can be the following cases:
1058   // a. tracks: in this case the label is taken and then case b.
1059   // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
1060   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
1061
1062   TClonesArray* arrayMC = 0;
1063   AliMCEvent* mcEvent = 0;
1064   if (mcObj->InheritsFrom("AliMCEvent"))
1065     mcEvent = static_cast<AliMCEvent*>(mcObj);
1066   else if (mcObj->InheritsFrom("TClonesArray"))
1067     arrayMC = static_cast<TClonesArray*>(mcObj);
1068   else
1069   {
1070     mcObj->Dump();
1071     AliFatal("Invalid object passed");
1072   }
1073
1074   Int_t before = tracks->GetEntriesFast();
1075
1076   for (Int_t i=0; i<before; ++i)
1077   {
1078     AliVParticle* part = (AliVParticle*) tracks->At(i);
1079
1080     if (part->InheritsFrom("AliESDtrack")) part = mcEvent->GetTrack(TMath::Abs(part->GetLabel()));
1081     if (part->InheritsFrom("AliAODTrack"))
1082     {
1083       part =(AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel()));
1084       //cout << "toto musi byt len pri reco trackoch" << endl;
1085     }
1086
1087     AliVParticle* mother = part;
1088     if (mcEvent)
1089     {
1090       while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
1091       {
1092     if (((AliMCParticle*)mother)->GetMother() < 0)
1093     {
1094       mother = 0;
1095       break;
1096     }
1097
1098     mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
1099     if (!mother)
1100       break;
1101       }
1102     }
1103     else
1104     {
1105       // find the primary mother
1106      while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
1107       {
1108     if (((AliAODMCParticle*)mother)->GetMother() < 0)
1109     {
1110       mother = 0;
1111       break;
1112     }
1113
1114     mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1115     if (!mother)
1116       break;
1117       }
1118     }
1119
1120     if (!mother)
1121     {
1122       //Printf("WARNING: No mother found for particle %d:", part->GetLabel());
1123       continue;
1124     }
1125
1126     if (mother->GetLabel() > maxLabel)
1127     {
1128 //      Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
1129       TObject* object = tracks->RemoveAt(i);
1130       if (tracks->IsOwner())
1131     delete object;
1132     }
1133   }
1134
1135   tracks->Compress();
1136
1137   AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
1138
1139 }
1140