From Rongrong: minor changes (binning etc)

[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Correlations / AliAnalysisTaskV0ChCorrelations.cxx

/**************************************************************************
 * Copyright(c) 1998-2013, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: Marek Bombara                                                  *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* The task selects candidates for K0s, Lambdas and AntiLambdas (trigger particles)
 * and calculates correlations with charged unidentified particles (associated particles) in phi and eta. 
 * The task works with AOD (with or without MC info) events only and containes also mixing for acceptance corrections.
 * Last update edited by Marek Bombara, August 2013, Marek.Bombara@cern.ch
 */

#include <TCanvas.h>
#include <TList.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <THnSparse.h>
#include <TObjArray.h>
#include <TObject.h>

#include "AliLog.h"
#include "AliAnalysisManager.h"
#include "AliAODTrack.h"
#include "AliAODEvent.h"
#include "AliAODv0.h"
#include "AliAODcascade.h"
#include "AliAODVertex.h"

#include "AliMCEvent.h"
#include "AliMCVertex.h" 
#include "AliGenEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAnalyseLeadingTrackUE.h"

#include "AliAODPid.h"
#include "AliPIDResponse.h"
#include "AliEventPoolManager.h"
#include "AliCentrality.h"

#include "AliAnalysisTaskV0ChCorrelations.h"
#include "AliPhysicsSelectionTask.h"

#include <AliMultiInputEventHandler.h>
#include <AliMixInputEventHandler.h>

ClassImp(AliAnalysisTaskV0ChCorrelations)
ClassImp(AliV0ChBasicParticle)

//________________________________________________________________________
AliAnalysisTaskV0ChCorrelations::AliAnalysisTaskV0ChCorrelations(const char *name) // All data members should be initialised here
   : AliAnalysisTaskSE(name),
     fAnalysisMC(kFALSE),
	 fFillMixed(kTRUE),
	 fMixingTracks(50000),
	 fPoolMgr(0x0),
     
     fDcaDToPV(0.5),
	 fDcaV0D(0.1),
	 fWithChCh(kFALSE),
	 fOStatus(1),

	 fOutput(0),
	 fPIDResponse(0),
     fHistCentVtx(0),
     fHistMultiMain(0),
	 fHistMassK0(0),
	 fHistMassLambda(0),
	 fHistMassAntiLambda(0),
 
	 fHistdPhidEtaSib(0),
	 fHistdPhidEtaMix(0),
	 fHistNTrigSib(0),

   	 fHistMCPtCentTrig(0),
     fHistRCPtCentTrig(0),
     fHistMCPtCentAs(0),
     fHistRCPtCentAs(0),
     fHistRCPtCentAll(0),
	 
	 fHistTemp(0),
	 fHistTemp2(0)// The last in the above list should not have a comma after it
{
   // Constructor
   // Define input and output slots here (never in the dummy constructor)
   // Input slot #0 works with a TChain - it is connected to the default input container
   // Output slot #1 writes into a TH1 container
   DefineOutput(1, TList::Class());                                            // for output list
}

//________________________________________________________________________
AliAnalysisTaskV0ChCorrelations::~AliAnalysisTaskV0ChCorrelations()
{
   // Destructor. Clean-up the output list, but not the histograms that are put inside
   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
      delete fOutput;
   }
}

//________________________________________________________________________
void AliAnalysisTaskV0ChCorrelations::UserCreateOutputObjects()
{
   // Create histograms

   fOutput = new TList();
   fOutput->SetOwner();  // IMPORTANT!
   // defining bins for centrality
   Int_t nCentralityBins  = 9;
   Double_t centBins[] = {0.,10.,20.,30.,40.,50.,60.,70.,80.,90.};
   const Double_t* centralityBins = centBins;
   // defining bins for Z vertex
   Int_t nZvtxBins  = 7;
   //Double_t vertexBins[] = {-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.};
   Double_t vertexBins[] = {-7.,-5.,-3.,-1.,1.,3.,5.,7.};
   const Double_t* zvtxBins = vertexBins;
   // pt bins of associated particles for the analysis
   Int_t nPtBins = 7;
   const Double_t PtBins[8] = {2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0}; 
   //Int_t nPtBins = 1;
   //const Double_t PtBins[2] = {3.0,15.0}; 
   // pt bins of trigger particles for the analysis
   Int_t nPtBinsV0 = 11;
   //const Double_t PtBinsV0[2] = {6.0,15.0}; 
   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}; 
   // V0 candidate: 1 - K0sig, 2 - Lamsig, 3 - Alamsig, 4 - K0bg, 5 - Lambg, 6 - Alambg
   Int_t nTrigC = 7;
   const Double_t TrigC[8] = {0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5};
   
   // Create general histograms 
   fHistCentVtx = new TH2F("fHistCentVtx", "Centrality vs. Z vertex", nCentralityBins, centralityBins, nZvtxBins, zvtxBins);
   fHistMultiMain = new TH1F("fHistMultiMain", "Multiplicity of main events", 2000, 0, 2000);

   const Int_t mrBins[3] = {nPtBinsV0, nCentralityBins, nTrigC};
   const Double_t mrMin[3] = {PtBinsV0[0], centralityBins[0], TrigC[0]};
   const Double_t mrMax[3] = {PtBinsV0[11], centralityBins[9], TrigC[6]};

   // Create histograms for reconstruction track and V0 efficiency
   fHistMCPtCentTrig = new THnSparseF("fHistMCPtCentTrig", "MC Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);
   fHistRCPtCentTrig = new THnSparseF("fHistRCPtCentTrig", "Rec Pt vs. Cent. Trig", 3, mrBins, mrMin, mrMax);

   // pt bins of associated particles for the efficiency
   Int_t nPtBinsAs = 13;
   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};

   fHistMCPtCentAs = new TH2D("fHistMCPtCentAs", "MC Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
   fHistRCPtCentAs = new TH2D("fHistRCPtCentAs", "Rec Pt vs. Cent. Assoc", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);
   fHistRCPtCentAll = new TH2D("fHistRCPtCentAll", "Rec Pt vs. Cent. All Prim+Sec", nPtBinsAs, PtBinsAs[0], PtBinsAs[13], nCentralityBins, centBins[0], centBins[9]);

   // defining bins for mass distributions
   Int_t nBins = 200;
   Double_t mk0Min = 0.40;
   Double_t mk0Max = 0.60;
   Double_t mlaMin = 1.07;
   Double_t mlaMax = 1.17;
   Double_t malMin = 1.07;
   Double_t malMax = 1.17;
   
   const Int_t spBins[3] = {nBins, nPtBinsV0, nCentralityBins};
   const Double_t spMinK0[3] = {mk0Min, PtBinsV0[0], centralityBins[0]};
   const Double_t spMaxK0[3] = {mk0Max, PtBinsV0[11], centralityBins[9]};
   const Double_t spMinLa[3] = {mlaMin, PtBinsV0[0], centralityBins[0]};
   const Double_t spMaxLa[3] = {mlaMax, PtBinsV0[11], centralityBins[9]};
   const Double_t spMinAl[3] = {malMin, PtBinsV0[0], centralityBins[0]};
   const Double_t spMaxAl[3] = {malMax, PtBinsV0[11], centralityBins[9]};
   // Create mass histograms
   fHistMassK0 = new THnSparseF("fHistMassK0","V0 mass for K0 hypothesis", 3, spBins, spMinK0, spMaxK0);
   fHistMassLambda = new THnSparseF("fHistMassLambda","V0 mass for Lambda hypothesis", 3, spBins, spMinLa, spMaxLa);
   fHistMassAntiLambda = new THnSparseF("fHistMassAntiLambda","V0 mass for AntiLambda hypothesis", 3, spBins, spMinAl, spMaxAl);
   
   // defining bins for dPhi distributions
   const Int_t nbPhiBins = 72;
   //const Double_t kPi = TMath::Pi();
   //Double_t PhiMin = -kPi/2.;
   //Double_t PhiMax = -kPi/2. + 2*kPi;
   Double_t PhiMin = -1.57;
   Double_t PhiMax = 4.71;
   Double_t PhiBins[nbPhiBins+1] = {0.};
   PhiBins[0] = PhiMin;
   for (Int_t i=0; i<nbPhiBins; i++) { PhiBins[i+1] = PhiBins[i] + (PhiMax - PhiMin)/nbPhiBins; }

   // defining bins for dEta distributions
   const Int_t nbEtaBins = 40;
   Double_t EtaMin = -2.0;
   Double_t EtaMax = 2.0;
   Double_t EtaBins[nbEtaBins+1] = {0.};
   EtaBins[0] = EtaMin;
   for (Int_t i=0; i<nbEtaBins; i++) { EtaBins[i+1] = EtaBins[i] + (EtaMax - EtaMin)/nbEtaBins; }

   const Int_t corBins[7] = {nbPhiBins, nbEtaBins, nPtBinsV0, nPtBins, nCentralityBins, nZvtxBins, nTrigC};
   const Double_t corMin[7] = {PhiBins[0], EtaBins[0], PtBinsV0[0], PtBins[0], centralityBins[0], zvtxBins[0], TrigC[0]};
   const Double_t corMax[7] = {PhiBins[72], EtaBins[40], PtBinsV0[11], PtBins[7], centralityBins[9], zvtxBins[7], TrigC[7]};
   // Create correlation histograms
   fHistdPhidEtaSib = new THnSparseF("fHistdPhidEtaSib","dPhi vs. dEta siblings", 7, corBins, corMin, corMax); 
   fHistdPhidEtaMix = new THnSparseF("fHistdPhidEtaMix","dPhi vs. dEta mixed", 7, corBins, corMin, corMax);

   // Create histograms for counting the numbers of trigger particles
   const Int_t corNTrigBins[5] = {nPtBinsV0, nCentralityBins, nZvtxBins, nTrigC};
   const Double_t corNTrigMin[5] = {PtBinsV0[0], centBins[0], vertexBins[0], TrigC[0]};
   const Double_t corNTrigMax[5] = {PtBinsV0[11], centBins[9], vertexBins[7], TrigC[7]};
   fHistNTrigSib = new THnSparseF("fHistNTrigSib","Number trigger sib", 4, corNTrigBins, corNTrigMin, corNTrigMax);

   // Histograms for debugging
   fHistTemp = new TH1D("fHistTemp", "Temporary", 100, -10., 10.);
   fHistTemp2 = new TH1D("fHistTemp2", "Temporary2", 100, -10., 10.);

   fOutput->Add(fHistCentVtx);

   fOutput->Add(fHistMultiMain);
   fOutput->Add(fHistMassK0);
   fOutput->Add(fHistMassLambda);
   fOutput->Add(fHistMassAntiLambda);
   
   fOutput->Add(fHistdPhidEtaSib);
   fOutput->Add(fHistdPhidEtaMix);
   fOutput->Add(fHistNTrigSib);

   fOutput->Add(fHistMCPtCentTrig);
   fOutput->Add(fHistRCPtCentTrig);
   fOutput->Add(fHistMCPtCentAs);
   fOutput->Add(fHistRCPtCentAs);
   fOutput->Add(fHistRCPtCentAll);
   
   fOutput->Add(fHistTemp);
   fOutput->Add(fHistTemp2);

   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram

   // Settings for event mixing -------------------------------------
   Int_t trackDepth = fMixingTracks;
   Int_t poolSize   = 200;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager

   fPoolMgr = new AliEventPoolManager(poolSize, trackDepth, nCentralityBins, centBins, nZvtxBins, vertexBins);
   //----------------------------------------------
}

//________________________________________________________________________
void AliAnalysisTaskV0ChCorrelations::Terminate(Option_t *)
{
   // Draw result to screen, or perform fitting, normalizations
   // Called once at the end of the query

   fOutput = dynamic_cast<TList*>(GetOutputData(1));
   if (!fOutput) { AliError("Could not retrieve TList fOutput"); return; }

   // NEW HISTO should be retrieved from the TList container in the above way,
   // so it is available to draw on a canvas such as below
}

//_________________________________________________________________________
void AliAnalysisTaskV0ChCorrelations::UserExec(Option_t *)
{
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
    AliInputEventHandler *inEvMain = (AliInputEventHandler*)mgr->GetInputEventHandler();

    // physics selection
    UInt_t maskIsSelected = inEvMain->IsEventSelected();
   // 2010 data trigger selection
    //Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
   // 2011 data trigger selection
    Bool_t isSelected = ((maskIsSelected & AliVEvent::kMB) || (maskIsSelected & AliVEvent::kCentral) || (maskIsSelected & AliVEvent::kSemiCentral));
    if (!isSelected) return;
	// calculating the event types
    if (maskIsSelected & AliVEvent::kMB) fHistTemp->Fill(2);
    if (maskIsSelected & AliVEvent::kCentral) fHistTemp->Fill(4);
    if (maskIsSelected & AliVEvent::kSemiCentral) fHistTemp->Fill(6);

    AliAODEvent* aod = (AliAODEvent*)inEvMain->GetEvent();
    fPIDResponse = inEvMain->GetPIDResponse();

  // pt intervals for trigger particles  
  const Double_t kPi = TMath::Pi();
  Double_t PtTrigMin = 4.0;
  Double_t PtTrigMax = 15.0;
  // pt interval for associated particles
  Double_t PtAssocMin = 2.0;
  fHistMultiMain->Fill(aod->GetNumberOfTracks());

  // Vertex cut
  Double_t cutPrimVertex = 7.0;
  AliAODVertex *myPrimVertex = aod->GetPrimaryVertex();
  if (!myPrimVertex) return;
  if ( ( TMath::Abs(myPrimVertex->GetZ()) ) >= cutPrimVertex) return ;
  Double_t lPVx = myPrimVertex->GetX();
  Double_t lPVy = myPrimVertex->GetY();
  Double_t lPVz = myPrimVertex->GetZ();

  if (TMath::Abs(lPVx)<10e-5 && TMath::Abs(lPVy)<10e-5 && TMath::Abs(lPVz)<10e-5) return;

  // Centrality definition
  Double_t lCent = 0.0;
  AliCentrality *centralityObj = 0;
  centralityObj = aod->GetHeader()->GetCentralityP();
  lCent = centralityObj->GetCentralityPercentile("V0M");
  if ((lCent < 0.)||(lCent > 90.)) return;
  fHistCentVtx->Fill(lCent,lPVz);

//=========== MC loop ===============================
  if (fAnalysisMC)
  {
    AliAODMCHeader *aodMCheader = (AliAODMCHeader*)aod->FindListObject(AliAODMCHeader::StdBranchName());
    Float_t vzMC = aodMCheader->GetVtxZ();
    if (TMath::Abs(vzMC) >= 7.) return;
    //retrieve MC particles from event
    TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
    if(!mcArray){
        Printf("No MC particle branch found");
        return;
    }
    
    Int_t nMCAllTracks = mcArray->GetEntriesFast();
    // new tracks array - without injected signal
    TObjArray * mcTracks = new TObjArray;
    //selectedMCTracks->SetOwner(kTRUE);
  
    for (Int_t i = 0; i < nMCAllTracks; i++)
    { 
        AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcArray->At(i);
        if (!mcTrack) {
            Error("ReadEventAODMC", "Could not receive particle %d", i);
            continue;
        }
        mcTracks->Add(mcTrack);
    }

    // get labels for injected particles -------
    TObject* mc = mcArray;
    //TObjArray * mcTracks = mcArray;
    Int_t skipParticlesAbove = 0;
    //AliGenEventHeader* eventHeader = 0;
    //Int_t headers = 0;
    //headers = aodMCheader->GetNCocktailHeaders();
    //eventHeader = aodMCheader->GetCocktailHeader(0);
    AliGenEventHeader* eventHeader = aodMCheader->GetCocktailHeader(0);
    Int_t headers = aodMCheader->GetNCocktailHeaders();
    if (!eventHeader)
    {
      // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
      // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
      AliError("First event header not found. Skipping this event.");
      return;
    }
    skipParticlesAbove = eventHeader->NProduced();
    //cout << "ski label = " << skipParticlesAbove << endl;
    AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s.Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
    //cout << "before MC compressing = " << mcTracks->GetEntriesFast() << endl;
    RemovingInjectedSignal(mcTracks,mc,skipParticlesAbove);
    //cout << "after MC compressing = " << mcTracks->GetEntriesFast() << endl;
    // -------------

    Int_t nMCTracks = mcTracks->GetEntriesFast();
    //cout << "number of MC tracks = " << nMCTracks << endl;
    for (Int_t iMC = 0; iMC<nMCTracks; iMC++)
    {
      AliAODMCParticle *mcTrack = (AliAODMCParticle*)mcTracks->At(iMC);
        if (!mcTrack) {
            Error("ReadEventAODMC", "Could not receive particle %d", iMC);
            continue;
      }
      // track part
      Double_t mcTrackEta = mcTrack->Eta();
      Double_t mcTrackPt = mcTrack->Pt();
      Bool_t TrIsPrim = mcTrack->IsPhysicalPrimary();
      Bool_t TrEtaMax = TMath::Abs(mcTrackEta)<0.9;
      Bool_t TrPtMin = mcTrackPt>PtAssocMin;
      Bool_t TrCharge = (mcTrack->Charge())!=0;
      if (TrIsPrim && TrEtaMax && TrPtMin && TrCharge) fHistMCPtCentAs->Fill(mcTrackPt,lCent);
      // V0 part
      Int_t mcMotherPdg = 0;
	  Int_t mcPartPdg = mcTrack->GetPdgCode();

      // Keep only K0s, Lambda and AntiLambda
	  if ((mcPartPdg != 310) && (mcPartPdg != 3122) && (mcPartPdg != (-3122))) continue;

      //cout << " mc pdg is " << mcPartPdg << endl;

      Bool_t IsK0 = mcPartPdg==310;
      Bool_t IsLambda = mcPartPdg==3122;
      Bool_t IsAntiLambda = mcPartPdg==-3122;
	  Bool_t IsSigma = kFALSE;
	  Int_t mcMotherLabel = mcTrack->GetMother();
      AliAODMCParticle *mcMother = (AliAODMCParticle*)mcArray->At(mcMotherLabel);
	  if (mcMotherLabel < 0) {mcMotherPdg = 0;} else {mcMotherPdg = mcMother->GetPdgCode();}
	  //if ((mcMotherLabel >= 0) && mcMother) 
	  //if ((mcMotherLabel >= 0) && mcMother) 
	  //{
      	//  Bool_t IsFromCascade = ((mcMotherPdg==3312)||(mcMotherPdg==-3312)||(mcMotherPdg==3334)||(mcMotherPdg==-3334));
          Bool_t IsFromSigma = ((mcMotherPdg==3212)||(mcMotherPdg==-3212));
          IsFromSigma = IsFromSigma || ((mcMotherPdg==3224)||(mcMotherPdg==-3224));
          IsFromSigma = IsFromSigma || ((mcMotherPdg==3214)||(mcMotherPdg==-3214));
          IsFromSigma = IsFromSigma || ((mcMotherPdg==3114)||(mcMotherPdg==-3114));
		  if ((IsFromSigma) && (mcMother->IsPhysicalPrimary()) && (IsLambda || IsAntiLambda)) IsSigma = kTRUE;
      	  Double_t mcRapidity = mcTrack->Y();
      	  Bool_t V0RapMax = TMath::Abs(mcRapidity)<0.75;
          Bool_t PtInterval = ((mcTrackPt>PtTrigMin)&&(mcTrackPt<PtTrigMax));
          //Bool_t PtInterval = kTRUE;
          IsK0 = IsK0 && (mcTrack->IsPhysicalPrimary());
          IsLambda = IsLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
          IsAntiLambda = IsAntiLambda && (mcTrack->IsPhysicalPrimary() || IsSigma);
      	  Double_t mcK0[3] = {mcTrackPt, lCent, 1};
      	  Double_t mcLa[3] = {mcTrackPt, lCent, 2};
      	  Double_t mcAl[3] = {mcTrackPt, lCent, 3};
      	  if (IsK0 && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcK0); 
      	  if (IsLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcLa);
      	  if (IsAntiLambda && V0RapMax && PtInterval) fHistMCPtCentTrig->Fill(mcAl);
	   //}
    }
    // ------- access the real data -----------
    Int_t nTracks = aod->GetNumberOfTracks();
    // new tracks array - without injected signal
    TObjArray * selectedMCTracks = new TObjArray;
    //selectedMCTracks->SetOwner(kTRUE);
  
    for (Int_t i = 0; i < nTracks; i++)
    {
        AliAODTrack* tr = aod->GetTrack(i); 
        selectedMCTracks->Add(tr);
    }
   // cout << "before compressing = " << selectedMCTracks->GetEntriesFast() << endl;
    RemovingInjectedSignal(selectedMCTracks,mc,skipParticlesAbove);
    //AliAnalyseLeadingTrackUE::RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
    //RemoveInjectedSignals(selectedMCTracks,mc,skipParticlesAbove);
   // cout << "after compressing = " << selectedMCTracks->GetEntriesFast() << endl;
    //-----------------------
    //cout << "before getting label -4" << endl;
    Int_t nRecTracks = selectedMCTracks->GetEntriesFast();
    //cout << "before getting label -3" << endl;
    for (Int_t i = 0; i < nRecTracks; i++)
    {
      //AliAODTrack* tras = aod->GetTrack(i);
      AliAODTrack* tras = (AliAODTrack*)selectedMCTracks->At(i);
      //cout << "before getting label -2" << endl;
      //cout << " and i = " << i << endl;
      //cout << "pt = " << tras->Pt() << " and i = " << i << endl;
      if ((tras->Pt())<PtAssocMin) continue; 
      //cout << "before getting label -1" << endl;
      if (!(IsMyGoodPrimaryTrack(tras))) continue;
      //cout << "before getting label 0" << endl;
      Int_t AssocLabel = tras->GetLabel();
      if (AssocLabel<=0) continue;
      //cout << "before getting label 1" << endl;
      Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
      //Bool_t isPhyPrim = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary();
      //if(!(static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->IsPhysicalPrimary())) continue;
	     //cout << "before getting label 2" << endl;
    Double_t mcPt = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Pt();
	     //cout << "before getting label 3" << endl;
    //Double_t mcEta = static_cast<AliAODMCParticle*>(mcArray->At(tras->GetLabel()))->Eta();
	  //if (mcPt<PtAssocMin) continue;
	  //if (TMath::Abs(mcEta)>0.8) continue;
      //Double_t trPt = tras->Pt();
      fHistRCPtCentAll->Fill(mcPt,lCent);
      if (isPhyPrim) fHistRCPtCentAs->Fill(mcPt,lCent);
    }
    // ------- end of real data access, for V0s see the main V0 loop -------- 
  }
//============= End of MC loop ======================
	
	// Track selection loop
	//--------------------------------
	Int_t nTracks = aod->GetNumberOfTracks();
	// new tracks array
	TObjArray * selectedTracks = new TObjArray;
	selectedTracks->SetOwner(kTRUE);
	
	Bool_t ChChWith = GetWithChCh();
    TObjArray * selectedChargedTriggers = new TObjArray;
    selectedChargedTriggers->SetOwner(kTRUE);
    for (Int_t i = 0; i < nTracks; i++)
    {
        AliAODTrack* tr = aod->GetTrack(i);
        if ((tr->Pt())<PtAssocMin) continue;
        if (!(IsMyGoodPrimaryTrack(tr))) continue;
        selectedTracks->Add(tr);
        // saving the Charged trigger particles
        if ((tr->Pt()>4.)&&(tr->Pt()<15.))
        {
           selectedChargedTriggers->Add(new AliV0ChBasicParticle(tr->Eta(), tr->Phi(), tr->Pt(), 7));
        }
    }

	//---------------------------------
	
	// V0 loop for reconstructed event
    TObjArray * selectedV0s = new TObjArray;
	selectedV0s->SetOwner(kTRUE);

	Int_t nV0s = aod->GetNumberOfV0s();

	for (Int_t i = 0; i < nV0s; i++)
	{ // start of V0 slection loop
		AliAODv0* aodV0 = dynamic_cast<AliAODv0 *>(aod->GetV0(i));
		if (!aodV0) {
         AliError(Form("ERROR: Could not retrieve aodv0 %d", i));
         continue;
        }
		if (aodV0->GetOnFlyStatus()) fHistTemp->Fill(6.);
		if (!aodV0->GetOnFlyStatus()) fHistTemp->Fill(8.);
		//cout << "pt of v0: " << aodV0->Pt() << endl;
		if (((aodV0->Pt())<PtTrigMin)||((aodV0->Pt())>PtTrigMax)) continue;
        // get daughters
        const AliAODTrack *myTrackPos;
        const AliAODTrack *myTrackNeg;
   	    AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
        AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));

        if (!myTrackPosTest || !myTrackNegTest) {
           Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
           continue;
        }

        if( myTrackPosTest->Charge() ==1){
           myTrackPos = myTrackPosTest;
           myTrackNeg = myTrackNegTest;
        }

        if( myTrackPosTest->Charge() ==-1){
            myTrackPos = myTrackNegTest;
            myTrackNeg = myTrackPosTest;
        }

        // effective mass calculations for each hypothesis
		Double_t lInvMassK0 = aodV0->MassK0Short();
        Double_t lInvMassAntiLambda = aodV0->MassAntiLambda();
        Double_t lInvMassLambda = aodV0->MassLambda();

		// calculations for c*tau cut--------------------------------------
	//	Double_t lDVx = aodV0->GetSecVtxX();
    //    Double_t lDVy = aodV0->GetSecVtxY();
    //    Double_t lDVz = aodV0->GetSecVtxZ();
    //    const Double_t kLambdaMass = 1.115683;
    //    const Double_t kK0Mass = 0.497648;
    //    Double_t cutcTauLam = 3*7.89;
    //    Double_t cutcTauK0 = 3*2.68;
     //   Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lDVx - lPVx,2) + TMath::Power(lDVy- lPVy,2) + TMath::Power(lDVz - lPVz,2 ));
      //  Double_t lPV0 = TMath::Sqrt((aodV0->Pt())*(aodV0->Pt())+(aodV0->Pz())*(aodV0->Pz()));
     //   Double_t lcTauLam = (lV0DecayLength*kLambdaMass)/lPV0;
     //   Double_t lcTauK0 = (lV0DecayLength*kK0Mass)/lPV0;
        // sc - standard cuts
		//Bool_t cutK0sc = (lcTauK0<cutcTauK0);
        //Bool_t cutLambdasc = (lcTauLam<cutcTauLam);
        //Bool_t cutAntiLambdasc = (lcTauLam<cutcTauLam);
		Bool_t cutK0sc = kTRUE;
        Bool_t cutLambdasc = kTRUE;
        Bool_t cutAntiLambdasc = kTRUE;

		//------------------------------------------------

		// preparations for daughter's PID cut------------------------------
        Float_t nSigmaPosPion   = 0.;
        Float_t nSigmaNegPion   = 0.;
        Float_t nSigmaPosProton = 0.;
        Float_t nSigmaNegProton = 0.;

        const AliAODPid *pPid = myTrackPos->GetDetPid();
        const AliAODPid *nPid = myTrackNeg->GetDetPid();

        if (pPid)
        {
            Double_t pdMom = pPid->GetTPCmomentum();
            if (pdMom<1.)
            {
                nSigmaPosPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
                nSigmaPosProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
            }
        }

        if (nPid)
        {
            Double_t ndMom = nPid->GetTPCmomentum();
            if (ndMom<1.)
            {
                nSigmaNegPion = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion);
                nSigmaNegProton = fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton);
            }
        }
		Bool_t bpPion = TMath::Abs(nSigmaPosPion) <= 3.;
        Bool_t bpProton = TMath::Abs(nSigmaPosProton) <= 3.;
        Bool_t bnPion = TMath::Abs(nSigmaNegPion) <= 3.;
        Bool_t bnProton = TMath::Abs(nSigmaNegProton) <= 3.;

        Bool_t cutK0Pid = (bpPion && bnPion);
        Bool_t cutLambdaPid = (bpProton && bnPion);
        Bool_t cutAntiLambdaPid = (bpPion && bnProton);
        //--------------------------------------------------
        // mass cuts
        Bool_t cutMassLambda = ((lInvMassLambda>1.10) && (lInvMassLambda<1.13));
        Bool_t cutMassAntiLambda = ((lInvMassAntiLambda>1.10) && (lInvMassAntiLambda<1.13));
        Bool_t cutMassK0 = (lInvMassK0>0.47) && (lInvMassK0<0.53);

        cutK0sc = cutK0Pid && (!cutMassLambda) && (!cutMassAntiLambda);
		cutLambdasc = cutLambdaPid && (!cutMassK0); 
		cutAntiLambdasc = cutAntiLambdaPid && (!cutMassK0);
        // special cut related to AP diagram for K0s
		Bool_t k0APcut = (aodV0->PtArmV0()>(TMath::Abs(0.2*aodV0->AlphaV0())));
		cutK0sc = cutK0sc && k0APcut;
        // fill the mass histograms

        Int_t oStatus = GetOStatus();

        if (!IsMyGoodV0(aod,aodV0,myTrackPos,myTrackNeg,oStatus)) continue;
		Double_t spK0[3] = {lInvMassK0, aodV0->Pt(), lCent};
		Double_t spLa[3] = {lInvMassLambda, aodV0->Pt(), lCent};
		Double_t spAl[3] = {lInvMassAntiLambda, aodV0->Pt(), lCent};
        if (cutK0sc) fHistMassK0->Fill(spK0);
		if (cutLambdasc) fHistMassLambda->Fill(spLa);
		if (cutAntiLambdasc) fHistMassAntiLambda->Fill(spAl);
        // select final V0s for correlation, selected background to study its contribution to correlation function
        // the values for signal might change in the future ...
        Bool_t K0Signal = (lInvMassK0>0.48)&&(lInvMassK0<0.52);
        Bool_t K0Bckg = ((lInvMassK0>0.40)&&(lInvMassK0<0.44)) || ((lInvMassK0>0.56)&&(lInvMassK0<0.60));

        Bool_t LamSignal = (lInvMassLambda>1.108)&&(lInvMassLambda<1.125);
        Bool_t LamBckg = ((lInvMassLambda>1.090)&&(lInvMassLambda<1.100)) || ((lInvMassLambda>1.135)&&(lInvMassLambda<1.145));

        Bool_t ALamSignal = (lInvMassAntiLambda>1.108)&&(lInvMassAntiLambda<1.125);
        Bool_t ALamBckg = ((lInvMassAntiLambda>1.090)&&(lInvMassAntiLambda<1.100)) || ((lInvMassAntiLambda>1.135)&&(lInvMassAntiLambda<1.145));
 
		// Fill selected V0 particle array 
        if ((cutK0sc)&&(K0Signal))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 1.));
            Double_t nTrigK0Sig[4] = {aodV0->Pt(), lCent, lPVz, 1.};
            fHistNTrigSib->Fill(nTrigK0Sig);
        }
        if ((cutK0sc)&&(K0Bckg))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 4.));
            Double_t nTrigK0Bkg[4] = {aodV0->Pt(), lCent, lPVz, 4.};
            fHistNTrigSib->Fill(nTrigK0Bkg);
        }
        if ((cutLambdasc)&&(LamSignal))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 2.));
            Double_t nTrigLaSig[4] = {aodV0->Pt(), lCent, lPVz, 2.};
            fHistNTrigSib->Fill(nTrigLaSig);
        }
        if ((cutLambdasc)&&(LamBckg))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 5.));
            Double_t nTrigLaBkg[4] = {aodV0->Pt(), lCent, lPVz, 5.};
            fHistNTrigSib->Fill(nTrigLaBkg);
        }
        if ((cutAntiLambdasc)&&(ALamSignal))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 3.));
            Double_t nTrigAlSig[4] = {aodV0->Pt(), lCent, lPVz, 3.};
            fHistNTrigSib->Fill(nTrigAlSig);
        }
        if ((cutAntiLambdasc)&&(ALamBckg))
        {
            selectedV0s->Add(new AliV0ChBasicParticle(aodV0->Eta(), aodV0->Phi(), aodV0->Pt(), 6.));
            Double_t nTrigAlBkg[4] = {aodV0->Pt(), lCent, lPVz, 6.};
            fHistNTrigSib->Fill(nTrigAlBkg);
        }


    	//===== MC part for V0 reconstruction efficiency ==============
    	if (fAnalysisMC)
    	{
      		TClonesArray *mcArray = (TClonesArray*)aod->FindListObject(AliAODMCParticle::StdBranchName());
      		if(!mcArray){
        	Printf("No MC particle branch found");
        	return;
      		}

			Int_t MotherOfMotherPdg = 0;

      		Int_t myTrackPosLabel = TMath::Abs(myTrackPos->GetLabel());
      		Int_t myTrackNegLabel = TMath::Abs(myTrackNeg->GetLabel());
      		AliAODMCParticle *mcPosTrack = (AliAODMCParticle*)mcArray->At(myTrackPosLabel);
			if (!mcPosTrack) continue;
      		AliAODMCParticle *mcNegTrack = (AliAODMCParticle*)mcArray->At(myTrackNegLabel);
			if (!mcNegTrack) continue;

      		Int_t PosTrackPdg = mcPosTrack->GetPdgCode();
      		Int_t NegTrackPdg = mcNegTrack->GetPdgCode();
        	//if (!mcPosTrack->IsPrimary()) continue;
        	//if (!mcNegTrack->IsPrimary()) continue;

      		Int_t myTrackPosMotherLabel = mcPosTrack->GetMother();
      		Int_t myTrackNegMotherLabel = mcNegTrack->GetMother();

      		if ((myTrackPosMotherLabel==-1)||(myTrackNegMotherLabel==-1)) continue;

      		AliAODMCParticle *mcPosMother = (AliAODMCParticle*)mcArray->At(myTrackPosMotherLabel);
			if (!mcPosMother) continue;	
      		Int_t MotherPdg = mcPosMother->GetPdgCode();
			Int_t MotherOfMother = mcPosMother->GetMother();
			//if (MotherOfMother == -1) MotherOfMotherPdg = 0;

			if (myTrackPosMotherLabel!=myTrackNegMotherLabel) continue;
        	//if (!mcPosMother->IsPrimary()) continue;
			Bool_t IsK0FromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==310)&&(PosTrackPdg==211)&&(NegTrackPdg==-211);
			Bool_t IsLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
			Bool_t IsAntiLambdaFromMC = (mcPosMother->IsPhysicalPrimary())&&(MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);

			Bool_t ComeFromSigma = kFALSE;
			Bool_t ComeFromSigmaLa = kFALSE;
			Bool_t ComeFromSigmaAl = kFALSE;
		
		    if (MotherOfMother != -1)
			{
      			AliAODMCParticle *mcPosMotherOfMother = (AliAODMCParticle*)mcArray->At(MotherOfMother);
				MotherOfMotherPdg = mcPosMotherOfMother->GetPdgCode();
            	Int_t MoMPdg = TMath::Abs(MotherOfMotherPdg); 
				ComeFromSigma = (mcPosMotherOfMother->IsPhysicalPrimary())&&((MoMPdg==3212)||(MoMPdg==3224)||(MoMPdg==3214)||(MoMPdg==3114));
				ComeFromSigmaLa = ComeFromSigma && (MotherPdg==3122)&&(PosTrackPdg==2212)&&(NegTrackPdg==-211);
				ComeFromSigmaAl = ComeFromSigma && (MotherPdg==-3122)&&(PosTrackPdg==211)&&(NegTrackPdg==-2212);
			}

            IsLambdaFromMC = IsLambdaFromMC || ComeFromSigmaLa;
            IsAntiLambdaFromMC = IsAntiLambdaFromMC || ComeFromSigmaAl;
			
      		Double_t RecMotherPt = aodV0->Pt();
			//cout << "Pt of rec v0 = " << RecMotherPt << " nMC = " << mcArray->GetEntries() << endl;
			//cout << "Pos Label = " << myTrackPosLabel << " Neg Label = " << myTrackNegLabel << endl;
      		Double_t rcK0[3] = {RecMotherPt, lCent, 1};
      		Double_t rcLa[3] = {RecMotherPt, lCent, 2};
      		Double_t rcAl[3] = {RecMotherPt, lCent, 3};
      		if ((cutK0sc)&&(K0Signal)&&(IsK0FromMC)) fHistRCPtCentTrig->Fill(rcK0);
      		if ((cutLambdasc)&&(LamSignal)&&(IsLambdaFromMC)) fHistRCPtCentTrig->Fill(rcLa);
      		if ((cutAntiLambdasc)&&(ALamSignal)&&(IsAntiLambdaFromMC)) fHistRCPtCentTrig->Fill(rcAl);

    	}

    	//===== End of the MC part for V0 reconstruction efficiency ===

		Int_t nSelectedTracks = selectedTracks->GetEntries();
		// Correlation part
		//===================================
		for (Int_t j = 0; j < nSelectedTracks; j++)
		{
			AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
			if ((atr->Pt())>=(aodV0->Pt())) continue;
			Double_t dEta = atr->Eta() - aodV0->Eta();
			Double_t dPhi = atr->Phi() - aodV0->Phi();
			if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
			if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;
			// removing autocorrelations 
			//----------------------------------
			Int_t negID = myTrackNeg->GetID();
			Int_t posID = myTrackPos->GetID();
			Int_t atrID = atr->GetID();

			if ((TMath::Abs(negID)+1)==(TMath::Abs(atrID))) continue;
			if ((TMath::Abs(posID)+1)==(TMath::Abs(atrID))) continue;
			//----------------------------------

			fHistNTrigSib->Sumw2();
			fHistdPhidEtaSib->Sumw2();
			// Filling correlation histograms and histograms for triggers counting
			//----------------- K0 ---------------------
			if ((cutK0sc)&&(K0Signal)) 
			{
				Double_t spK0Sig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 1.};
				fHistdPhidEtaSib->Fill(spK0Sig);
			}

			if ((cutK0sc)&&(K0Bckg)) 
			{
				Double_t spK0Bkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 4.};
				fHistdPhidEtaSib->Fill(spK0Bkg);
			}
			//---------------- Lambda -------------------
			if ((cutLambdasc)&&(LamSignal)) 
			{
				Double_t spLaSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 2.};
				fHistdPhidEtaSib->Fill(spLaSig);
			}

			if ((cutLambdasc)&&(LamBckg)) 
			{
				Double_t spLaBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 5.};
				fHistdPhidEtaSib->Fill(spLaBkg);
			}
			//------------- AntiLambda -------------------
			if ((cutAntiLambdasc)&&(ALamSignal)) 
			{
				Double_t spAlSig[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 3.};
				fHistdPhidEtaSib->Fill(spAlSig);
			}

			if ((cutAntiLambdasc)&&(ALamBckg)) 
			{
				Double_t spAlBkg[7] = {dPhi, dEta, aodV0->Pt(), atr->Pt(), lCent, lPVz, 6.};
				fHistdPhidEtaSib->Fill(spAlBkg);
			}
	
		} // end of correlation loop
		//===================================


	} // end of V0 selection loop

	// ===========================================
    // Ch-Ch correlation part

    if (ChChWith)
    {
        Int_t nSelectedChargedTriggers = selectedChargedTriggers->GetEntries();
        for (Int_t i = 0; i < nSelectedChargedTriggers; i++)
        {
            AliV0ChBasicParticle* chTrig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
            Double_t chTrigPhi = chTrig->Phi();
            Double_t chTrigEta = chTrig->Eta();
            Double_t chTrigPt = chTrig->Pt();

            Double_t nTrigChSig[4] = {chTrigPt, lCent, lPVz, 7.};
            fHistNTrigSib->Fill(nTrigChSig);

            Int_t nSelectedTracks = selectedTracks->GetEntries();
            for (Int_t j = 0; j < nSelectedTracks; j++)
            {
                AliAODTrack* atr = (AliAODTrack*) selectedTracks->At(j);
                if ((atr->Pt())>=chTrigPt) continue;
                Double_t dEta = atr->Eta() - chTrigEta;
                Double_t dPhi = atr->Phi() - chTrigPhi;
                if (dPhi > (1.5*kPi)) dPhi -= 2.0*kPi;
                if (dPhi < (-0.5*kPi)) dPhi += 2.0*kPi;

                Double_t spCh[7] = {dPhi, dEta, chTrigPt, atr->Pt(), lCent, lPVz, 7.};
                fHistdPhidEtaSib->Fill(spCh);
            }
        }
    }
    // end of Ch-Ch correlation part

	// Mixing ==============================================
	
	fHistdPhidEtaMix->Sumw2();
    AliEventPool* pool = fPoolMgr->GetEventPool(lCent, lPVz);
	if (!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", lCent, lPVz));
	//pool->SetDebug(1);
    if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
	{
		Int_t nMix = pool->GetCurrentNEvents();
		for (Int_t jMix=0; jMix<nMix; jMix++)
		{// loop through mixing events
			TObjArray* bgTracks = pool->GetEvent(jMix);
			for (Int_t i=0; i<selectedV0s->GetEntriesFast(); i++)
			{// loop through selected V0 particles
				AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedV0s->At(i);
				Double_t trigPhi = trig->Phi();
				Double_t trigEta = trig->Eta();
				Double_t trigPt = trig->Pt();
				Short_t trigC = trig->WhichCandidate();
				for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
				{ // mixing tracks loop 
					AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
					// be careful tracks may have bigger pt than v0s.
					if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;
					Double_t dEtaMix = assoc->Eta() - trigEta;
					Double_t dPhiMix = assoc->Phi() - trigPhi;
					if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
					if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;
				
					Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
				    fHistdPhidEtaMix->Fill(spMix);
					
				} // end of mixing track loop
			}// end of loop through selected V0 particles
			if (ChChWith)
            {
                for (Int_t i=0; i<selectedChargedTriggers->GetEntriesFast(); i++)
                {// loop through selected charged trigger particles
                    AliV0ChBasicParticle* trig = (AliV0ChBasicParticle*) selectedChargedTriggers->At(i);
                    Double_t trigPhi = trig->Phi();
                    Double_t trigEta = trig->Eta();
                    Double_t trigPt = trig->Pt();
                    Short_t trigC = trig->WhichCandidate();
                    for (Int_t j=0; j<bgTracks->GetEntriesFast(); j++)
                    { // mixing tracks loop 
                        AliVParticle* assoc = (AliVParticle*) bgTracks->At(j);
                        // be careful tracks may have bigger pt than v0s.
                        if ( ( (assoc->Pt())>=trigPt )||( (assoc->Pt())<PtAssocMin ) ) continue;

                        Double_t dEtaMix = assoc->Eta() - trigEta;
                        Double_t dPhiMix = assoc->Phi() - trigPhi;
                        if (dPhiMix > (1.5*kPi)) dPhiMix -= 2.0*kPi;
                        if (dPhiMix < (-0.5*kPi)) dPhiMix += 2.0*kPi;

                        Double_t spMix[7] = {dPhiMix, dEtaMix, trigPt, assoc->Pt(), lCent, lPVz, trigC};
                        fHistdPhidEtaMix->Fill(spMix);

                    } // end of mixing track loop
                }// end of loop through selected Ch particles
            } // end of ChCh mixing

		}// end of loop of mixing events
	}

	TObjArray* tracksClone = (TObjArray*) selectedTracks->Clone();
	tracksClone->SetOwner(kTRUE);
	pool->UpdatePool(tracksClone);
	//delete selectedtracks;

    PostData(1, fOutput);

}
//___________________________________________
Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodPrimaryTrack(const AliAODTrack *t)
{
	// Pseudorapidity cut
	if (TMath::Abs(t->Eta())>0.9) return kFALSE;
	// Should correspond to set of cuts suitable for correlation analysis, 768 - hybrid tracks for 2011 data
    if (!t->TestFilterBit(768)) return kFALSE;

	return kTRUE;
}
//_____________________________________________
Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodDaughterTrack(const AliAODTrack *t)
{
	// TPC refit
	if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
	// Minimum number of clusters
	Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
	if (nCrossedRowsTPC < 70) return kFALSE;
	Int_t findable=t->GetTPCNclsF();
	if (findable <= 0) return kFALSE;
	if (nCrossedRowsTPC/findable < 0.8) return kFALSE;

	return kTRUE;
}
//______________________________________________
Bool_t AliAnalysisTaskV0ChCorrelations::IsMyGoodV0(const AliAODEvent* aod, const AliAODv0* aodV0, const AliAODTrack* myTrackPos, const AliAODTrack* myTrackNeg, Int_t oSta)
{
	if (!aodV0) {
       AliError(Form("ERROR: Could not retrieve aodV0"));
       return kFALSE;
    }

	// Offline reconstructed V0 only
    if (oSta==1) {if (aodV0->GetOnFlyStatus()) return kFALSE;}
    if (oSta==3) {if (!aodV0->GetOnFlyStatus()) return kFALSE;}
   
	if (oSta==2) 
	{
		if (aodV0->GetOnFlyStatus()) 
		{
			return kTRUE;
		} else {
			return kFALSE;
		}
	}
    if (oSta==4) 
	{
		if (!aodV0->GetOnFlyStatus()) 
		{
			return kTRUE;
		} else {
			return kFALSE;
		}
	}
    
    // Get daughters and check them
	AliAODTrack *myTrackNegTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(1));
	AliAODTrack *myTrackPosTest=dynamic_cast<AliAODTrack *>(aodV0->GetDaughter(0));
	
	if (!myTrackPosTest || !myTrackNegTest) {
		Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
		return kFALSE;
	}

    if( myTrackPosTest->Charge() ==1){
            myTrackPos = myTrackPosTest;
            myTrackNeg = myTrackNegTest;
    }

    if( myTrackPosTest->Charge() ==-1){
            myTrackPos = myTrackNegTest;
            myTrackNeg = myTrackPosTest;
    }

	// Track cuts for daughter tracks
    if ( !(IsMyGoodDaughterTrack(myTrackPos)) || !(IsMyGoodDaughterTrack(myTrackNeg)) ) return kFALSE;

	// Unlike signs of daughters
    if (myTrackNegTest->Charge() == myTrackPosTest->Charge()) return kFALSE;

	// Rapidity cut
	//Double_t lCutRap = 0.75;
	//Double_t lRapK0s = aodV0->Y(310);
	//Double_t lRapLambda = aodV0->Y(3122);
	//Double_t lRapAntiLambda = aodV0->Y(-3122);

	// Pseudorapidity cut - there are high pt V0
	Double_t lCutEta = 0.75;
	Double_t lEtaV0 = aodV0->Eta();
	if (TMath::Abs(lEtaV0)>=lCutEta) return kFALSE;
	//if (TMath::Abs(lRapK0s)>=lCutRap) return kFALSE;
	//if (TMath::Abs(lRapLambda)>=lCutRap) return kFALSE;
	//if (TMath::Abs(lRapAntiLambda)>=lCutRap) return kFALSE;
    
    // getting global variables
	Float_t dcaDaughtersToPrimVtx = GetDcaDToPV();
	Float_t dcaBetweenDaughters = GetDcaV0D();

    // DCA of daughter track to Primary Vertex
    Float_t xyn=aodV0->DcaNegToPrimVertex();
    if (TMath::Abs(xyn)<dcaDaughtersToPrimVtx) return kFALSE;
    Float_t xyp=aodV0->DcaPosToPrimVertex();
    if (TMath::Abs(xyp)<dcaDaughtersToPrimVtx) return kFALSE;

	// DCA of daughter tracks 
    Double_t dca=aodV0->DcaV0Daughters();
    if (dca>dcaBetweenDaughters) return kFALSE;

	// Cosinus of pointing angle
    Double_t cpa=aodV0->CosPointingAngle(aod->GetPrimaryVertex());
    if (cpa<0.998) return kFALSE;

	// Fiducial volume cut
    Double_t xyz[3]; aodV0->GetSecondaryVtx(xyz);
    Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
    if (r2<0.9*0.9) return kFALSE;
    if (r2>100*100) return kFALSE;

	// c*tau cut - in main V0 loop - depends on particle hypothesis

	// Minimum pt of daughters
    Double_t  lMomPos[3] = {999,999,999};
	Double_t  lMomNeg[3] = {999,999,999};

	lMomPos[0] = aodV0->MomPosX();
    lMomPos[1] = aodV0->MomPosY();
    lMomPos[2] = aodV0->MomPosZ();

    lMomNeg[0] = aodV0->MomNegX();
    lMomNeg[1] = aodV0->MomNegY();
    lMomNeg[2] = aodV0->MomNegZ();

    Double_t lPtPos = TMath::Sqrt(lMomPos[0]*lMomPos[0] + lMomPos[1]*lMomPos[1]);
    Double_t lPtNeg = TMath::Sqrt(lMomNeg[0]*lMomNeg[0] + lMomNeg[1]*lMomNeg[1]);
	
	Double_t cutMinPtDaughter = 0.160;
	if (lPtPos<cutMinPtDaughter || lPtNeg<cutMinPtDaughter) return kFALSE;

	// Daughter PID cut - in main V0 loop - depends on particle hypothesis

	return kTRUE;
}

void AliAnalysisTaskV0ChCorrelations::RemovingInjectedSignal(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
{
  // remove injected signals (primaries above <maxLabel>)
  // <tracks> can be the following cases:
  // a. tracks: in this case the label is taken and then case b.
  // b. particles: the first stable mother is searched and checked if it is <= <maxLabel>
  // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)

  TClonesArray* arrayMC = 0;
  AliMCEvent* mcEvent = 0;
  if (mcObj->InheritsFrom("AliMCEvent"))
    mcEvent = static_cast<AliMCEvent*>(mcObj);
  else if (mcObj->InheritsFrom("TClonesArray"))
    arrayMC = static_cast<TClonesArray*>(mcObj);
  else
  {
    mcObj->Dump();
    AliFatal("Invalid object passed");
  }

  Int_t before = tracks->GetEntriesFast();

  for (Int_t i=0; i<before; ++i)
  {
    AliVParticle* part = (AliVParticle*) tracks->At(i);

    if (part->InheritsFrom("AliESDtrack")) part = mcEvent->GetTrack(TMath::Abs(part->GetLabel()));
    if (part->InheritsFrom("AliAODTrack"))
    {
      part =(AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel()));
      //cout << "toto musi byt len pri reco trackoch" << endl;
    }

    AliVParticle* mother = part;
    if (mcEvent)
    {
      while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
      {
    if (((AliMCParticle*)mother)->GetMother() < 0)
    {
      mother = 0;
      break;
    }

    mother = (AliMCParticle*) mcEvent->GetTrack(((AliMCParticle*)mother)->GetMother());
    if (!mother)
      break;
      }
    }
    else
    {
      // find the primary mother
     while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
      {
    if (((AliAODMCParticle*)mother)->GetMother() < 0)
    {
      mother = 0;
      break;
    }

    mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
    if (!mother)
      break;
      }
    }

    if (!mother)
    {
      //Printf("WARNING: No mother found for particle %d:", part->GetLabel());
      continue;
    }

    if (mother->GetLabel() > maxLabel)
    {
//      Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
      TObject* object = tracks->RemoveAt(i);
      if (tracks->IsOwner())
    delete object;
    }
  }

  tracks->Compress();

  AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));

}