PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskV0ForRAA.cxx

   1 /***************************************************************          *
   2  * Authors : Simone Schuchmann
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //-----------------------------------------------------------------
  15 // AliAnalysisTaskV0ForRAA class
  16 // This task is for analysing Lambda and K0s pt spectra in PbPb and
  17 // pp as well as with MC. The flag for pp and MC  must be set
  18 // accordingly, default is PbPb data.
  19 // It works with ESD files only.
  20 //-----------------------------------------------------------------
  21
  22
  23 #define AliAnalysisTaskV0ForRAA_cxx
  24
  25 #include "AliAnalysisTaskV0ForRAA.h"
  26
  27 #include "Riostream.h"
  28
  29 #include "TROOT.h"
  30 #include "TH1.h"
  31 #include "TH2.h"
  32 #include "TLorentzVector.h"
  33
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliMCEventHandler.h"
  38
  39 #include "AliESDVZERO.h"
  40 #include "AliMultiplicity.h"
  41 #include "AliCentrality.h"
  42
  43 #include "AliKFParticle.h"
  44 #include "AliKFVertex.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliESDpid.h"
  47 #include "AliESDv0.h"
  48 #include "AliESDEvent.h"
  49 #include "AliMCEvent.h"
  50 #include "AliStack.h"
  51 #include "AliGenEventHeader.h"
  52
  53
  54 ClassImp(AliAnalysisTaskV0ForRAA)
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskV0ForRAA::AliAnalysisTaskV0ForRAA(const char *name)
  58 :AliAnalysisTaskSE(name),
  59   fESD(0),
  60   fMCev(0),
  61 //----------other objects------------
  62   fESDpid(0),
  63   fESDTrackCuts(0),
  64   fESDTrackCutsCharged(0),
  65   fESDTrackCutsLowPt(0),
  66   fOutputContainer(0),
  67 //-----------event histos-----------
  68   fHistITSLayerHits(0),
  69   fHistOneHitWithSDD(0),
  70   fHistNEvents(0),
  71   fHistPrimVtxZESDVSNContributors(0),
  72   fHistPrimVtxZESDTPCVSNContributors(0),
  73   fHistPrimVtxZESDSPDVSNContributors(0),
  74   fHistPrimVtxZESDVSNContributorsMC(0),
  75   fHistPrimVtxZESDTPCVSNContributorsMC(0),
  76   fHistPrimVtxZESDSPDVSNContributorsMC(0),
  77   fHistPrimVtxZESD(0),
  78   fHistPrimVtxZESDTPC(0),
  79   fHistPrimVtxZESDSPD(0),
  80   fHistESDVertexZ(0),
  81   fHistMCVertexZ(0),
  82   fHistMuliplicity(0),
  83   fHistMuliplicityRaw(0),
  84   fHistCentBinRaw(0),
  85   fHistCentBin(0),
  86   fHistMultiplicityPrimary(0),
  87   fHistNPrim(0),
  88 //--------MC pdg code histos-------------
  89 //K0s
  90   fHistPiPiPDGCode(0),
  91 //Lambda
  92   fHistPiPPDGCode(0),
  93   fHistPiPCosPointAngXiVsPt(0),
  94   fHistPiPMassVSPtSecXiMCTruth(0),
  95   fHistPiPMassVSPtSecOmegaMCTruth(0),
  96 //AntiLambda
  97   fHistPiAPPDGCode(0),
  98   fHistPiAPCosPointAngXiVsPt(0),
  99   fHistPiAPMassVSPtSecXiMCTruth(0),
 100   fHistPiAPMassVSPtSecOmegaMCTruth(0),
 101 // fHistUserPtShift(0),
 102
 103 //----------selection booleans and values---------
 104   fMCMode(0),
 105   fMCTruthMode(0),
 106   fSelectInjected(0),
 107   fSelectMBMotherMC(0),
 108   fCheckNegLabelReco(0),
 109   fOnlyFoundRecoV0(0),
 110   fUseCentrality(0),
 111   fUseCentralityBin(0),
 112   fUseCentralityRange(0),
 113   fAnapp(0),
 114   fSelSDD(0),
 115   fSelNoSDD(0),
 116   fOntheFly(0),
 117   fVertexZCut(0),
 118   fVtxStatus(0),
 119   fUsePID(0),
 120   fNSigma(0),
 121   fPPIDcut(0),
 122   fPtTPCCut(0),
 123   fMoreNclsThanRows(0),
 124   fMoreNclsThanFindable(0),
 125   fMoreNclsThanFindableMax(0),
 126   fRatioFoundOverFindable(0),
 127   fRatioMaxCRowsOverFindable(0),
 128   fChi2PerClusterITS(0),
 129   fRapCutV0(0),
 130   fRap(0),
 131   fEtaCutMCDaughters(0),
 132   fEtaCutMCDaughtersVal(0),
 133   fMinPt(0),
 134   fAlfaCut(0),
 135   fQtCut(0),
 136   fQtCutPtLow(0),
 137   fArmCutK0(0),
 138   fArmCutL(0),
 139   fArmQtSlope(0),
 140   fExcludeLambdaFromK0s(0),
 141   fExcludeK0sFromLambda(0),
 142   fDCAToVertexK0(0),
 143   fDCAToVertexL(0),
 144   fDCAXK(0),
 145   fDCAYK(0),
 146   fDCAXL(0),
 147   fDCAYL(0),
 148   fDCAZ(0),
 149   fDCADaughtersL(0),
 150   fDCADaughtersAL(0),
 151   fDCADaughtersK0(0),
 152   fDCADaughtersToVtxLarge(0),
 153   fDCADaughtersToVtxSmall(0),
 154   fDecayRadXYMin(0),
 155   fDecayRadXYMax(0),
 156   fCosPointAngL(0),
 157   fCosPointAngK(0),
 158   fCPAPtCutK0(0),
 159   fCPAPtCutL(0),
 160   fOpengAngleDaughters(0),
 161   fOpAngPtCut(0),
 162   fDecayLengthMax(0),
 163   fDecayLengthMin(0),
 164   fCtauK0s(0),
 165   fCtauL(0),
 166   fCtauPtCutK0(0),
 167   fCtauPtCutL(0),
 168   fChiCutKf(0)
 169   //  fShift(0),
 170   // fDeltaInvP(0)
 171 {  // Constructor
 172
 173   DefineOutput(1,TList::Class());
 174
 175   // define defaults for globals
 176   /*
 177     fShift = kFALSE;                       // shift in charge/pt yes/no
 178     fDeltaInvP = 0.00;                     // shift value
 179   */
 180
 181   //MC general
 182   fMCMode      = kFALSE; // run on MC reco or MC as data
 183   fMCTruthMode = kFALSE; // run on MC truth (needed kTRUE for MC reco as well)
 184   //MC mode
 185   fSelectInjected    = kFALSE; // reject injected particles in MC
 186   fSelectMBMotherMC  = kFALSE; // reject injected mothers of particles in MC
 187   fCheckNegLabelReco = kFALSE; // check if particle with negative reco label is bad in MC truth
 188   fOnlyFoundRecoV0   = kFALSE; // select only found reco V0 in MC truth
 189
 190   //centrality
 191   fUseCentrality      = 0; // flag for centrality usage
 192   fUseCentralityBin   = 0; // select centrality bin
 193   fUseCentralityRange = 0; // set centrality bin range
 194
 195   //pp analysis
 196   fAnapp    = kFALSE; // analyse pp collisions
 197   fSelSDD   = kFALSE; // select events with SDD
 198   fSelNoSDD = kFALSE; // select events without SDD
 199
 200   //vertex
 201   fVertexZCut = 100000.0; // select z vertex window
 202   fVtxStatus  = kFALSE;   // check vertex status
 203
 204   //V0 finder
 205   fOntheFly   = kTRUE;    // on-the-fly V0 finder default
 206
 207
 208   //----- define defaults for V0 and track cuts ----//
 209
 210   fMoreNclsThanRows          = kFALSE; // cut on more TPC clusters than crossed rows
 211   fMoreNclsThanFindable      = kFALSE; // cut on more TPC clusters than findable clusters
 212   fMoreNclsThanFindableMax   = kFALSE; // cut on more TPC clusters than findalbe max
 213   fRatioFoundOverFindable    = -1.0;   // cut value on ratio TPC clusters found over findable
 214   fRatioMaxCRowsOverFindable = 1000.0; // cut value on max ratio TPC crossed rows over findable clusters
 215
 216   fChi2PerClusterITS = 100000.0; // cut value on chi2 per ITS cluster
 217
 218   //PID
 219   fUsePID                    = kFALSE; // use PID dEdx TPC
 220   fNSigma            = 100000.0; // number of sigmas in TPC dEdx
 221   fPPIDcut           = 100.0;    // select momentum for PID cut
 222   fPtTPCCut          = -1.0;     // select pt value for low and high pt esd track cuts
 223
 224
 225   fRapCutV0 = kFALSE; // cut on rapidity
 226   fRap      = 1000.0; // cut value on rapidity
 227
 228   fAlfaCut    = -100.0; // cut value on alpha of armenteros-podolanski diagram
 229   fQtCut      =  100.0; // cut value on qt of armenteros-podolanski diagram
 230   fQtCutPtLow = -1.0;   // select pt max for cut on qt in armenteros-podolanski diagram
 231   fArmCutK0   = kFALSE; // cut in armenteros-podolanski diagram  for K0s
 232   fArmCutL    = kFALSE; // cut in armenteros-podolanski diagram  for Lambda
 233   fArmQtSlope = 0.2;    // slope of diagonal cut in armenteros-podolanski diagram  for K0s
 234
 235   fExcludeLambdaFromK0s = -1.0; //cut value on mass exclusion of Lambdas from K0s
 236   fExcludeK0sFromLambda = -1.0; //cut value on mass exclusion of K0s from Lambdas
 237
 238   fEtaCutMCDaughters    = kFALSE; // cut on eta in MC truth
 239   fEtaCutMCDaughtersVal = 50.0;   // cut value on eta for data, MC reco and truth
 240
 241   fMinPt= -1.0; // cut value on min pt
 242
 243   //DCA to vertex in 2D and 1D cut value
 244   fDCAToVertexK0 = 10000.0;// 2D K0s
 245   fDCAToVertexL  = 10000.0;// 2D Lambda
 246   fDCAXK         = 10000.0;// 1D x K0s
 247   fDCAYK         = 10000.0;// 1D y K0s
 248   fDCAXL         = 10000.0;// 1D x Lambda
 249   fDCAYL         = 10000.0;// 1D y Lambda
 250   fDCAZ          = 10000.0;// 1D z all
 251
 252   // DCA between daughters cut value
 253   fDCADaughtersL  = 10000.0;
 254   fDCADaughtersAL = 10000.0;
 255   fDCADaughtersK0 = 10000.0;
 256
 257   // DCA of daughters to vertex cut value
 258   fDCADaughtersToVtxLarge = -1.0; //large value ((a)proton)
 259   fDCADaughtersToVtxSmall = -1.0; //small value (pion)
 260
 261   //decay radii 2D cut value
 262   fDecayRadXYMin  = -100000.0;
 263   fDecayRadXYMax  =  100000.0;
 264
 265   //decay length 3D cut value
 266   fDecayLengthMax =  100000.0;
 267   fDecayLengthMin = -100000.0;
 268
 269   //cosine of pointing angle cut value
 270   fCosPointAngL = -1.0;
 271   fCosPointAngK = -1.0;
 272   fCPAPtCutK0   = 1000.0; // pt max for cut
 273   fCPAPtCutL    = 1000.0; // pt max for cut
 274   //opening angle between daughters cut value
 275   fOpengAngleDaughters = -1.0;
 276   fOpAngPtCut          = -1.0; // pt max for cut
 277
 278   //ctau cut values
 279   fCtauK0s     = 1000.0;
 280   fCtauL       = 1000.0;
 281   fCtauPtCutK0 = 1000.0; // pt max for cut
 282   fCtauPtCutL  = 1000.0; // pt max for cut
 283
 284   fChiCutKf = 1000000.0; // cut value on chi2 of KF particle mass fit
 285
 286
 287   //---- histograms ----//
 288   for(Int_t j=0;j<2;j++){
 289     fHistV0RadiusZ[j]     =NULL;
 290     fHistV0RadiusZVSPt[j] =NULL;
 291     fHistV0RadiusXY[j]    =NULL;
 292     fHistV0RadiusXYVSY[j] =NULL;
 293     fHistArmenteros[j]    =NULL;
 294
 295     if(j<1){ //K0
 296       fHistPiPiMass[j]                         =NULL;
 297       fHistPiPiMassVSPt[j]                     =NULL;
 298       fHistPiPiMassVSPtMCTruth[j]              =NULL;
 299       // fHistPiPiMassVSAlpha[j]=NULL;
 300       fHistPiPiRadiusXY[j]                     =NULL;
 301       fHistPiPiCosPointAng[j]                  =NULL;
 302       fHistPiPiDecayLengthVsPt[j]              =NULL;
 303       fHistPiPiDecayLengthVsMass[j]            =NULL;
 304       fHistPiPiDecayLengthVsCtau[j]            =NULL;
 305       fHistPiPiDCADaughterPosToPrimVtxVSMass[j]=NULL;
 306       fHistPiPiDCADaughters[j]                 =NULL;
 307       //    fHistPiPiPtDaughters[j]=NULL;
 308       fHistPiPiPtVSY[j]                        =NULL;
 309       fHistPiPiDCAVSMass[j]                    =NULL;
 310       fHistPiPiMonitorCuts[j]                  =NULL;
 311       fHistPiPiMonitorMCCuts[j]                =NULL;
 312       fHistPiPiDecayLengthResolution[j]        =NULL;
 313
 314       fHistNclsITSPosK0[j]                     =NULL;
 315       fHistNclsITSNegK0[j]                     =NULL;
 316       fHistNclsTPCPosK0[j]                     =NULL;
 317       fHistNclsTPCNegK0[j]                     =NULL;
 318       fHistChi2PerNclsITSPosK0[j]              =NULL;
 319       fHistChi2PerNclsITSNegK0[j]              =NULL;
 320     }
 321
 322     //Lambda
 323     fHistPiPMass[j]                            =NULL;
 324     fHistPiPMassVSPt[j]                        =NULL;
 325     fHistPiPMassVSPtMCTruth[j]                 =NULL;
 326     fHistPiPRadiusXY[j]                        =NULL;
 327     fHistPiPCosPointAng[j]                     =NULL;
 328     fHistPiPDecayLengthVsPt[j]                 =NULL;
 329     fHistPiPDecayLengthVsMass[j]               =NULL;
 330     fHistPiPDecayLengthVsCtau[j]               =NULL;
 331     fHistPiPDCADaughterPosToPrimVtxVSMass[j]   =NULL;
 332     fHistPiPMassVSPtSecSigma[j]                =NULL;
 333     fHistPiPMassVSPtSecXi[j]                   =NULL;
 334     fHistPiPMassVSPtSecOmega[j]                =NULL;
 335     fHistPiPMassVSYSecXi[j]                    =NULL;
 336     fHistPiPXi0PtVSLambdaPt[j]                 =NULL;
 337     fHistPiPXiMinusPtVSLambdaPt[j]             =NULL;
 338     fHistPiPOmegaPtVSLambdaPt[j]               =NULL;
 339     fHistPiPDCADaughters[j]                    =NULL;
 340     //  fHistPiPPtDaughters[j]=NULL;
 341     fHistPiPPtVSY[j]                           =NULL;
 342     fHistPiPDCAVSMass[j]                       =NULL;
 343     fHistPiPMonitorCuts[j]                     =NULL;
 344     fHistPiPMonitorMCCuts[j]                   =NULL;
 345     fHistPiPDecayLengthResolution[j]           =NULL;
 346
 347     //AntiLambda
 348     fHistPiAPMass[j]                           =NULL;
 349     fHistPiAPMassVSPt[j]                       =NULL;
 350     fHistPiAPMassVSPtMCTruth[j]                =NULL;
 351     fHistPiAPRadiusXY[j]                       =NULL;
 352     fHistPiAPCosPointAng[j]                    =NULL;
 353     fHistPiAPDecayLengthVsPt[j]                =NULL;
 354     fHistPiAPDecayLengthVsMass[j]              =NULL;
 355     fHistPiAPDecayLengthVsCtau[j]              =NULL;
 356     fHistPiAPDCADaughterPosToPrimVtxVSMass[j]  =NULL;
 357     fHistPiAPMassVSPtSecSigma[j]               =NULL;
 358     fHistPiAPMassVSPtSecXi[j]                  =NULL;
 359     fHistPiAPMassVSPtSecOmega[j]               =NULL;
 360     fHistPiAPMassVSYSecXi[j]                   =NULL;
 361     fHistPiAPXi0PtVSLambdaPt[j]                =NULL;
 362     fHistPiAPXiMinusPtVSLambdaPt[j]            =NULL;
 363     fHistPiAPOmegaPtVSLambdaPt[j]              =NULL;
 364     fHistPiAPDCADaughters[j]                   =NULL;
 365     // fHistPiAPPtDaughters[j]=NULL;
 366     fHistPiAPPtVSY[j]                          =NULL;
 367     fHistPiAPDCAVSMass[j]                      =NULL;
 368     fHistPiAPMonitorCuts[j]                    =NULL;
 369     fHistPiAPMonitorMCCuts[j]                  =NULL;
 370     fHistPiAPDecayLengthResolution[j]          =NULL;
 371
 372     //----- other ------------------//
 373     //dedx
 374     fHistDedxSecProt[j]               =NULL;
 375     fHistDedxSecAProt[j]              =NULL;
 376     fHistDedxSecPiMinus[j]            =NULL;
 377     fHistDedxSecPiPlus[j]             =NULL;
 378     fHistDedxProt[j]                  =NULL;
 379     fHistDedxAProt[j]                 =NULL;
 380     fHistDedxPiMinus[j]               =NULL;
 381     fHistDedxPiPlus[j]                =NULL;
 382     //TPC Lambda
 383     fHistNclsITSPosL[j]               =NULL;
 384     fHistNclsITSNegL[j]               =NULL;
 385     fHistNclsTPCPosL[j]               =NULL;
 386     fHistNclsTPCNegL[j]               =NULL;
 387     fHistChi2PerNclsITSPosL[j]        =NULL;
 388     fHistChi2PerNclsITSNegL[j]        =NULL;
 389     //TPC all
 390     fHistNclsITSPos[j]                =NULL;
 391     fHistNclsITSNeg[j]                =NULL;
 392     fHistNclsTPCPos[j]                =NULL;
 393     fHistNclsTPCNeg[j]                =NULL;
 394     fHistChi2PerNclsITSPos[j]         =NULL;
 395     fHistChi2PerNclsITSNeg[j]         =NULL;
 396     fHistNclsITS[j]                   =NULL;
 397     fHistNclsTPC[j]                   =NULL;
 398     fHistNCRowsTPCPos[j]              =NULL;
 399     fHistNCRowsTPCNeg[j]              =NULL;
 400     fHistRatioFoundOverFinableTPCK0[j]=NULL;
 401     fHistRatioFoundOverFinableTPCL[j] =NULL;
 402     //eta
 403     fHistPiPiEtaDMC[j]                =NULL;
 404     fHistPiPiEtaDReco[j]              =NULL;
 405     fHistPiPEtaDMC[j]                 =NULL;
 406     fHistPiPEtaDReco[j]               =NULL;
 407   }
 408
 409 }
 410
 411 //_____________________________________________________
 412 AliAnalysisTaskV0ForRAA::~AliAnalysisTaskV0ForRAA()
 413 {
 414   //---- Remove all pointers ----//
 415   if(fOutputContainer)     delete fOutputContainer;    fOutputContainer     = 0;
 416   if(fESDTrackCuts)        delete fESDTrackCuts;       fESDTrackCuts        = 0;
 417   if(fESDTrackCutsCharged) delete fESDTrackCutsCharged;fESDTrackCutsCharged = 0;
 418   if(fESDTrackCutsLowPt)   delete fESDTrackCutsLowPt;  fESDTrackCutsLowPt   = 0;
 419 }
 420
 421 //________________________________________________________________________
 422 void AliAnalysisTaskV0ForRAA::UserCreateOutputObjects(){
 423   //create output objects
 424
 425   Int_t nbPt   = 800;
 426   Int_t nbMass = 500;
 427
 428
 429   //-----------------  create output container -----------------//
 430
 431   fOutputContainer = new TList() ;
 432   fOutputContainer->SetName(GetName()) ;
 433   fOutputContainer->SetOwner();
 434
 435   Int_t mchist = 1;   //for data or MC as data
 436   if((fMCMode && fMCTruthMode) || fMCTruthMode) mchist = 2; //for MC reco or truth
 437
 438
 439   //------------ create allways -----------------------//
 440   fHistNEvents = new TH1F("fHistNEvents","no of events before cuts =0, after cuts=1, after process =2",5,0.0,5.0);
 441   fOutputContainer->Add(fHistNEvents);
 442
 443   fHistMuliplicity =  new TH1F("fHistMuliplicity","V0 multiplicity",3000,0.0,30000);
 444   fOutputContainer->Add(fHistMuliplicity);
 445
 446   fHistMuliplicityRaw =  new TH1F("fHistMuliplicityRaw","V0 multiplicity before process",3000,0.0,30000);
 447   fOutputContainer->Add(fHistMuliplicityRaw);
 448
 449   fHistMultiplicityPrimary = new TH1F("fHistMultiplicityPrimary","number of charged tracks",5000,0.0,20000);
 450   fOutputContainer->Add(fHistMultiplicityPrimary);
 451
 452   fHistESDVertexZ  = new TH1F("fHistESDVertexZ"," z vertex distr in cm",500,-50,50);
 453   fOutputContainer->Add(fHistESDVertexZ);
 454
 455   fHistPrimVtxZESD = new TH1F("fHistPrimVtxZESD","z vertex pos ESD",250,-50,50);
 456   fOutputContainer->Add(fHistPrimVtxZESD);
 457
 458   fHistPrimVtxZESDVSNContributors = new TH2F("fHistPrimVtxZESDVSNContributors","prim vtx pos z ESD vs no. of contributers TPC",250,-50,50,500,0.0,500.0);
 459   fOutputContainer->Add(fHistPrimVtxZESDVSNContributors);
 460
 461   fHistNPrim       = new TH1F("fHistNPrim","Number of contributers to vertex",2500,0.0,5000);
 462   fOutputContainer->Add(fHistNPrim);
 463
 464   //------------------------ pp analysis only -------------------------//
 465   if(fAnapp){
 466     fHistITSLayerHits   = new TH1F("fHistITSLayerHits","SDD layer -1=0,1=1,2=2 ... 5=5,0=nothing",7,-1.5,5.5);
 467     fOutputContainer->Add(fHistITSLayerHits);
 468     fHistOneHitWithSDD  = new TH1F("fHistOneHitWithSDD","min one hit in SDD",2,-0.5,1.5);
 469     fOutputContainer->Add(fHistOneHitWithSDD);
 470     fHistPrimVtxZESDTPC = new TH1F("fHistPrimVtxZESDTPC","z vertex pos TPC",250,-50,50);
 471     fOutputContainer->Add(fHistPrimVtxZESDTPC);
 472     fHistPrimVtxZESDSPD = new TH1F("fHistPrimVtxZESDSPD","z vertex pos SPD",250,-50,50);
 473     fOutputContainer->Add(fHistPrimVtxZESDSPD);
 474     fHistPrimVtxZESDTPCVSNContributors = new TH2F("fHistPrimVtxZESDTPCVSNContributors","prim vtx pos z TPC vs no. of contributers TPC",250,-50,50,500,0.0,500.0);
 475     fOutputContainer->Add(fHistPrimVtxZESDTPCVSNContributors);
 476     fHistPrimVtxZESDSPDVSNContributors = new TH2F("fHistPrimVtxZESDSPDVSNContributors","prim vtx pos z SPD vs no. of contributers TPC",250,-50,50,500,0.0,500.0);
 477     fOutputContainer->Add(fHistPrimVtxZESDSPDVSNContributors);
 478
 479   }
 480   else {
 481     Double_t binsCent[12]={0.0,5.0,10.0,20.0,30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0};
 482     fHistCentBinRaw = new TH1F("fHistCentBinRaw","centrality bin before cent selection",11,binsCent);
 483     fOutputContainer->Add(fHistCentBinRaw);
 484     fHistCentBin    = new TH1F("fHistCentBin","centrality bin",11,binsCent);
 485     fOutputContainer->Add(fHistCentBin);
 486
 487   }
 488
 489   // ------------------- add always ---------------------------//
 490   fHistV0RadiusZ[0]      = new TH2F("fHistV0RadiusZ","z of decay radius vs 2D radius",100,0.0,100.0,250,-125.0,125.0);
 491   fHistV0RadiusZVSPt[0]  = new TH2F("fHistV0RadiusZVSPt","z of decay radius vs pt radius",200,0.0,20.0,125,0.0,125.0);
 492   fHistV0RadiusXY[0]     = new TH2F("fHistV0RadiusXY","y vs x decay radius",250,-125.0,125.0,250,-125.0,125.0);
 493   fHistV0RadiusXYVSY[0]  = new TH2F("fHistV0RadiusXYVSY","2D decay radius vs rap",100,-1,1,100,0.0,100.0);
 494   fHistArmenteros[0]     = new TH2F("fHistArmenteros"," pi+pi- armenteros",nbMass,-1.,1.,500,0.,0.5);
 495
 496
 497   //***************************************** K0s **********************************//
 498   fHistPiPiMass[0]              = new TH1F("fHistPiPiMass"," pi+pi- InvMass distribution",2*nbMass,0.,2.);
 499   fHistPiPiMassVSPt[0]          = new TH2F("fHistPiPiMassVSPt","pi+pi- InvMass distribution",nbMass,0.25,0.75,300,0.0,30.0);
 500   fHistPiPiMassVSPtMCTruth[0]   = new TH2F("fHistPiPiMassVSPtMCTruth","pi+pi- InvMass distribution vs pt MCTruth",nbMass,0.25,0.75,300,0.0,30.0);
 501   fHistPiPiPtVSY[0]             = new TH2F("fHistPiPiPtVSY","phi vs mass",100,-1,1,100,0.0,20);
 502   fHistPiPiDecayLengthVsPt[0]   = new TH2F("fHistPiPiDecayLengthVsPt","K0 decay length vs pt",200,0.0,20.0,220,0.0,110.0);
 503   fHistPiPiDecayLengthVsCtau[0] = new TH2F("fHistPiPiDecayLengthVsCtau","K0s ctau vs mass",nbMass,0.25,0.75,250,0.0,50.0);
 504   fHistPiPiDecayLengthVsMass[0] = new TH2F("fHistPiPiDecayLengthVsMass","K0s decay length vs mass",nbMass,0.25,0.75,220,0.0,110.0);
 505   fHistPiPiMonitorCuts[0]       = new TH1F("fHistPiPiMonitorCuts","K0 cut monitor",25,0.5,25.5);
 506   fHistPiPiMonitorMCCuts[0]     = new TH1F("fHistPiPiMonitorMCCuts","K0 cut monitor mc",25,0.5,25.5);
 507
 508   //***************************************** Lambda **********************************//
 509   fHistPiPMass[0]               = new TH1F("fHistPiPMass"," p+pi- InvMass distribution",2*nbMass,0.,2.);
 510   fHistPiPMassVSPt[0]           = new TH2F("fHistPiPMassVSPt","p+pi- InvMass distribution",nbMass,1.05,1.25,300,0.0,30.0);
 511   fHistPiPMassVSPtMCTruth[0]    = new TH2F("fHistPiPMassVSPtMCTruth","p+pi- InvMass distribution vs pt MCTruth",nbMass,1.05,1.25,300,0.0,30.0);
 512   fHistPiPPtVSY[0]              = new TH2F("fHistPiPPtVSY","p{t} vs y",100,-1,1,100,0.0,20);
 513   fHistPiPDecayLengthVsPt[0]    = new TH2F("fHistPiPDecayLengthVsPt","#Lambda decay length vs pt",200,0.0,20.0,220,0.0,110.0);
 514   fHistPiPDecayLengthVsCtau[0]  = new TH2F("fHistPiPDecayLengthVsCtau","L ctau vs mass",nbMass,1.05,1.25,250,0.0,50.0);
 515   fHistPiPDecayLengthVsMass[0]  = new TH2F("fHistPiPDecayLengthVsMass","#Lambda decay length vs mass",nbMass,1.05,1.25,220,0.0,110.0);
 516   fHistPiPMonitorCuts[0]        = new TH1F("fHistPiPMonitorCuts","#Lambda cut monitor",25,0.5,25.5);
 517   fHistPiPMonitorMCCuts[0]      = new TH1F("fHistPiPMonitorMCCuts","#Lambda cut monitor mc ",25,0.5,25.5);
 518
 519   //***************************************** AntiLambda **********************************//
 520   fHistPiAPMass[0]              = new TH1F("fHistPiAPMass"," ap-pi+ InvMass distribution",2*nbMass,0.,2.);
 521   fHistPiAPMassVSPt[0]          = new TH2F("fHistPiAPMassVSPt","p-pi+ InvMass distribution",nbMass,1.05,1.25,300,0.0,30.0);
 522   fHistPiAPMassVSPtMCTruth[0]   = new TH2F("fHistPiAPMassVSPtMCTruth","p-pi+ InvMass distribution vs pt MCTruth",nbMass,1.05,1.25,300,0.0,30.0);
 523   fHistPiAPPtVSY[0]             = new TH2F("fHistPiAPPtVSY","p{t} vs y",100,-1,1,100,0.0,20);
 524   fHistPiAPDecayLengthVsPt[0]   = new TH2F("fHistPiAPDecayLengthVsPt","#bar{#Lambda} decay length vs pt",200,0.0,20.0,220,0.0,110.0);
 525   fHistPiAPDecayLengthVsCtau[0] = new TH2F("fHistPiAPDecayLengthVsCtau","AL ctau vs mass",nbMass,1.05,1.25,200,0.0,100.0);
 526   fHistPiAPDecayLengthVsMass[0] = new TH2F("fHistPiAPDecayLengthVsMass","#bar{#Lambda} decay length vs mass",nbMass,1.05,1.25,220,0.0,110.0);
 527   fHistPiAPMonitorCuts[0]       = new TH1F("fHistPiAPMonitorCuts","#bar{#Lambda} cut monitor",25,0.5,25.5);
 528   fHistPiAPMonitorMCCuts[0]     = new TH1F("fHistPiAPMonitorMCCuts","#bar{#Lambda} cut monitor mc",25,0.5,25.5);
 529
 530   if(mchist==2){// for MC reco
 531     fHistV0RadiusZ[1]      = new TH2F("fHistV0RadiusZSec","z of decay radius vs 2D radius",100,0.0,100.0,250,-125.0,125.0);
 532     fHistV0RadiusZVSPt[1]  = new TH2F("fHistV0RadiusZVSPtSec","z of decay radius vs pt radius",200,0.0,20.0,125,0.0,125.0);
 533     fHistV0RadiusXY[1]     = new TH2F("fHistV0RadiusXYSec","y vs x decay radius",250,-125.0,125.0,250,-125.0,125.0);
 534     fHistV0RadiusXYVSY[1]  = new TH2F("fHistV0RadiusXYVSYSec","2D decay radius vs rap",100,-1,1,100,0.0,100.0);
 535     fHistArmenteros[1]     = new TH2F("fHistArmenterosSec"," pi+pi- armenteros",nbMass,-1.,1.,500,0.,0.5);
 536
 537     //***************************************** K0s **********************************//
 538
 539     //***************************************** Lambda **********************************//
 540     fHistPiPMass[1]               = new TH1F("fHistPiPMassSec"," p+pi- InvMass distribution",2*nbMass,0.,2.);
 541     fHistPiPMassVSPt[1]           = new TH2F("fHistPiPMassVSPtSec","p+pi- InvMass distribution",nbMass,1.05,1.25,300,0.0,30.0);
 542     fHistPiPMassVSPtMCTruth[1]    = new TH2F("fHistPiPMassVSPtMCTruthSec","p+pi- InvMass distribution vs pt MCTruth",nbMass,1.05,1.25,300,0.0,30.0);
 543     fHistPiPPtVSY[1]              = new TH2F("fHistPiPPtVSYSec","p{t} vs y",100,-1,1,100,0.0,20);
 544     fHistPiPDecayLengthVsPt[1]    = new TH2F("fHistPiPDecayLengthVsPtSec","#Lambda decay length vs pt",200,0.0,20.0,220,0.0,11.0);
 545     fHistPiPDecayLengthVsCtau[1]  = new TH2F("fHistPiPDecayLengthVsCtauSec","AL ctau vs mass",nbMass,1.05,1.25,250,0.0,50.0);
 546     fHistPiPDecayLengthVsMass[1]  = new TH2F("fHistPiPDecayLengthVsMassSec","#Lambda decay length vs mass",nbMass,1.05,1.25,220,0.0,110.0);
 547     fHistPiPMonitorCuts[1]        = new TH1F("fHistPiPMonitorCutsSec","#Lambda cut monitor",25,0.5,25.5);
 548     fHistPiPMonitorMCCuts[1]      = new TH1F("fHistPiPMonitorMCCutsSec","#Lambda cut monitor mc",25,0.5,25.5);
 549     //***************************************** AntiLambda **********************************//
 550     fHistPiAPMass[1]              = new TH1F("fHistPiAPMassSec"," ap-pi+ InvMass distribution",2*nbMass,0.,2.);
 551     fHistPiAPMassVSPt[1]          = new TH2F("fHistPiAPMassVSPtSec","p-pi+ InvMass distribution",nbMass,1.05,1.25,300,0.0,30.0);
 552     fHistPiAPMassVSPtMCTruth[1]   = new TH2F("fHistPiAPMassVSPtMCTruthSec","p-pi+ InvMass distribution vs pt MCTruth",nbMass,1.05,1.25,300,0.0,30.0);
 553     fHistPiAPPtVSY[1]             = new TH2F("fHistPiAPPtVSYSec","p{t} vs y",100,-1,1,100,0.0,20);
 554     fHistPiAPDecayLengthVsPt[1]   = new TH2F("fHistPiAPDecayLengthVsPtSec","#bar{#Lambda} decay length vs pt",200,0.0,20.0,220,0.0,11.0);
 555     fHistPiAPDecayLengthVsCtau[1] = new TH2F("fHistPiAPDecayLengthVsCtauSec","AL ctau vs mass",nbMass,1.05,1.25,250,0.0,50.0);
 556     fHistPiAPDecayLengthVsMass[1] = new TH2F("fHistPiAPDecayLengthVsMassSec","#bar{#Lambda} decay length vs mass",nbMass,1.05,1.25,220,0.0,110.0);
 557     fHistPiAPMonitorCuts[1]       = new TH1F("fHistPiAPMonitorCutsSec","#bar{#Lambda} cut monitor",25,0.5,25.5);
 558     fHistPiAPMonitorMCCuts[1]     = new TH1F("fHistPiAPMonitorMCCutsSec","#bar{#Lambda} cut monitor mc",25,0.5,25.5);
 559   }
 560
 561   //------- add K0s histos -----//
 562   fOutputContainer->Add(fHistPiPiMass[0]);
 563   fOutputContainer->Add(fHistPiPiMassVSPt[0]);
 564   fOutputContainer->Add(fHistPiPiMassVSPtMCTruth[0]);
 565   fOutputContainer->Add(fHistPiPiPtVSY[0]);
 566   fOutputContainer->Add(fHistPiPiDecayLengthVsPt[0]);
 567   fOutputContainer->Add(fHistPiPiDecayLengthVsCtau[0]);
 568   fOutputContainer->Add(fHistPiPiDecayLengthVsMass[0]);
 569   fOutputContainer->Add(fHistPiPiMonitorCuts[0]);
 570   fOutputContainer->Add(fHistPiPiMonitorMCCuts[0]);
 571
 572   //------- add (A)Lambda histos -----//
 573   for(Int_t j=0;j<mchist;j++){
 574
 575     fOutputContainer->Add(fHistArmenteros[j]);
 576     fOutputContainer->Add(fHistV0RadiusZ[j]);
 577     fOutputContainer->Add(fHistV0RadiusZVSPt[j]);
 578     fOutputContainer->Add(fHistV0RadiusXY[j]);
 579     fOutputContainer->Add(fHistV0RadiusXYVSY[j]);
 580
 581     fOutputContainer->Add(fHistPiPMass[j]);
 582     fOutputContainer->Add(fHistPiAPMass[j]);
 583
 584     fOutputContainer->Add(fHistPiPMassVSPt[j]);
 585     fOutputContainer->Add(fHistPiAPMassVSPt[j]);
 586
 587     fOutputContainer->Add(fHistPiPMassVSPtMCTruth[j]);
 588     fOutputContainer->Add(fHistPiAPMassVSPtMCTruth[j]);
 589
 590     fOutputContainer->Add(fHistPiPPtVSY[j]);
 591     fOutputContainer->Add(fHistPiAPPtVSY[j]);
 592
 593     fOutputContainer->Add(fHistPiPDecayLengthVsPt[j]);
 594     fOutputContainer->Add(fHistPiAPDecayLengthVsPt[j]);
 595
 596     fOutputContainer->Add(fHistPiPDecayLengthVsCtau[j]);
 597     fOutputContainer->Add(fHistPiAPDecayLengthVsCtau[j]);
 598
 599     fOutputContainer->Add(fHistPiPDecayLengthVsMass[j]);
 600     fOutputContainer->Add(fHistPiAPDecayLengthVsMass[j]);
 601
 602     fOutputContainer->Add(fHistPiPMonitorCuts[j]);
 603     fOutputContainer->Add(fHistPiAPMonitorCuts[j]);
 604
 605     fOutputContainer->Add(fHistPiPMonitorMCCuts[j]);
 606     fOutputContainer->Add(fHistPiAPMonitorMCCuts[j]);
 607   }
 608
 609   //----------------- for reco or data or mc data like MC reco only -----------------//
 610   if((fMCMode) || (!fMCTruthMode && !fMCMode)){//reco or data
 611
 612     fHistPiPiEtaDReco[0] = new TH2F("fHistPiPiEtaDRecoRaw","K0s daughters eta raw",300,-6,6,100,0,20);
 613     fOutputContainer->Add(fHistPiPiEtaDReco[0]);
 614     fHistPiPiEtaDReco[1] = new TH2F("fHistPiPiEtaDReco","K0s daughters eta after rap V0 cut pos",300,-3,3,300,-3.00,3.0);
 615     fOutputContainer->Add(fHistPiPiEtaDReco[1]);
 616     fHistPiPEtaDReco[0]  = new TH2F("fHistPiPEtaDRecoRaw","#Lambda daughters eta raw",300,-6,6,100,0,20);
 617     fOutputContainer->Add(fHistPiPEtaDReco[0]);
 618     fHistPiPEtaDReco[1]  = new TH2F("fHistPiPEtaDReco","#Lambda daughters eta after rap V0 cut neg",300,-3,3,300,-3.00,3.0);
 619     fOutputContainer->Add(fHistPiPEtaDReco[1]);
 620
 621     //********************************************** K0 *************************************************//
 622     // fHistPiPiMassVSAlpha[0] = new TH2F("fHistPiPiMassVSAlpha"," alpha armenteros vs pi+pi- InvMass distribution",nbMass,0.25,0.75,500,-1.,1.);
 623     fHistPiPiDCADaughters[0]                  = new TH2F("fHistPiPiDCADaughters","dca of K0 daughters",nbMass,0.25,0.75,250,0.0,2);
 624     fHistPiPiDCADaughterPosToPrimVtxVSMass[0] = new TH2F("fHistPiPiDCADaughterPosToPrimVtxVSMass","pi+ DCA daughter to prim vtx vsinvmass",nbMass,0.25,0.75,250,0.0,25.0);
 625     fHistPiPiDCAVSMass[0]                     = new TH2F("fHistPiPiDCAVSMass","pi+pi- dca  vs pt",nbMass,0.25,0.75,250,0.0,5.0);
 626     fHistPiPiCosPointAng[0]                   = new TH2F("fHistPiPiCosPointAng","K0 cosine of pointing angle vs mass ",nbMass,0.25,0.75,200,0.9,1.00);
 627     fHistPiPiRadiusXY[0]                      = new TH2F("fHistPiPiRadiusXY","pi+pi- phi dist vs mass",nbMass,0.25,0.75,200,0.0,4.0);
 628     // fHistPiPiPtDaughters[0] = new TH2F("fHistPiPiPtDaughters","p_{t} pos vs p_{t} neg of daughters",400,0.0,20.0,400,0,20.0);
 629
 630     //***************************************** Lambda *********************************************//
 631     fHistPiPDCADaughters[0] = new TH2F("fHistPiPDCADaughters","dca of #Lambda daughters",nbMass,1.05,1.25,250,0.0,2.0);
 632     fHistPiPDCADaughterPosToPrimVtxVSMass[0]  = new TH2F("fHistPiPDCADaughterPosToPrimVtxVSMass","p DCA daughter to prim vtx vs invmass",nbMass,1.05,1.25,250,0.0,25.0);
 633     fHistPiPDCAVSMass[0]                      = new TH2F("fHistPiPDCAVSMass","ppi- dca  vs pt",nbMass,1.05,1.25,250,0.0,5.0);
 634     fHistPiPCosPointAng[0]                    = new TH2F("fHistPiPCosPointAng","#Lambda cosine of pointing angle vs mass ",nbMass,1.05,1.25,200,0.9,1.00);
 635     fHistPiPRadiusXY[0]                       = new TH2F("fHistPiPRadiusXY","pi-p+ phi dist vs mass",nbMass,1.05,1.25,200,0.0,4.0);
 636     // fHistPiPPtDaughters[0] = new TH2F("fHistPiPPtDaughters","p_{t} pos vs p_{t} neg of daughters",400,0.0,20.0,400,0,20.0);
 637
 638     //************************************** Antilambda **********************************************//
 639     fHistPiAPDCADaughters[0]                  = new TH2F("fHistPiAPDCADaughters","dca of #bar{#Lambda} daughters",nbMass,1.05,1.25,250,0.0,2.0);
 640     fHistPiAPDCADaughterPosToPrimVtxVSMass[0] = new TH2F("fHistPiAPDCADaughterPosToPrimVtxVSMass","pi+ DCA daughter to prim vtx vs invmass",nbMass,1.05,1.25,250,0.0,25.0);
 641     fHistPiAPDCAVSMass[0]                     = new TH2F("fHistPiAPDCAVSMass","pi+p- dca  vs pt",nbMass,1.05,1.25,250,0.0,5.0);
 642     fHistPiAPCosPointAng [0]                  = new TH2F("fHistPiAPCosPointAng","#bar{#Lambda} cosine of pointing angle vs mass",nbMass,1.05,1.25,200,0.9,1.00);
 643     fHistPiAPRadiusXY[0]                      = new TH2F("fHistPiAPRadiusXY","pi+p- phi dist vs mass",nbMass,1.05,1.25,200,0.0,4.0);
 644     // fHistPiAPPtDaughters[0] = new TH2F("fHistPiAPPtDaughters","p_{t} pos vs p_{t} neg of daughters",400,0.0,20.0,400,0,20.0);
 645
 646     //**********************************************TPC*****************************************************//
 647
 648     fHistDedxSecProt[0]         = new TH2F("fHistDedxSecProt","proton", nbPt, 0, 20, 100, 0, 400);
 649     fHistDedxSecPiPlus[0]       = new TH2F("fHistDedxSecPiPlus","pi plus", nbPt, 0, 20, 100, 0, 400);
 650     fHistDedxSecAProt[0]        = new TH2F("fHistDedxSecAProt","antiproton", nbPt, 0, 20, 100, 0, 400);
 651     fHistDedxSecPiMinus[0]      = new TH2F("fHistDedxSecPiMinus","pi minus", nbPt, 0, 20, 100, 0, 400);
 652     fHistDedxProt[0]            = new TH2F("fHistDedxProt","proton", nbPt, 0, 20, 100, 0, 400);
 653     fHistDedxPiPlus[0]          = new TH2F("fHistDedxPiPlus","pi plus", nbPt, 0, 20, 100, 0, 400);
 654     fHistDedxAProt[0]           = new TH2F("fHistDedxAProt","antiproton", nbPt, 0, 20, 100, 0, 400);
 655     fHistDedxPiMinus[0]         = new TH2F("fHistDedxPiMinus","pi minus", nbPt, 0, 20, 100, 0, 400);
 656     //K0s
 657     fHistNclsITSPosK0[0]        = new TH1F("fHistNclsITSPosK0","fHistNclsITSPos K0",10,-0.5,9.5);
 658     fHistNclsITSNegK0[0]        = new TH1F("fHistNclsITSNegK0","fHistNclsITSNeg K0",10,-0.5,9.5);
 659
 660     fHistNclsTPCPosK0[0]        = new TH1F("fHistNclsTPCPosK0","fHistNclsTPCPos K0",200,-0.5,199.5);
 661     fHistNclsTPCNegK0[0]        = new TH1F("fHistNclsTPCNegK0","fHistNclsTPCNeg K0",200,-0.5,199.5);
 662
 663     fHistChi2PerNclsITSPosK0[0] = new TH1F("fHistChi2PerNclsITSPosK0","chi2 per cluster ITS pi+pi-",250,0.0,50.0);
 664     fHistChi2PerNclsITSNegK0[0] = new TH1F("fHistChi2PerNclsITSNegK0","chi2 per cluster ITS pi+pi- neg",250,0.0,50.0);
 665
 666     //Lambda
 667     fHistNclsITSPosL[0]         = new TH1F("fHistNclsITSPosL","fHistNclsITSPos #Lambda",10,-0.5,9.5);
 668     fHistNclsITSNegL[0]         = new TH1F("fHistNclsITSNegL","fHistNclsITSNeg #Lambda",10,-0.5,9.5);
 669
 670     fHistNclsTPCPosL[0]         = new TH1F("fHistNclsTPCPosL","fHistNclsTPCPos #Lambda",200,-0.5,199.5);
 671     fHistNclsTPCNegL[0]         = new TH1F("fHistNclsTPCNegL","fHistNclsTPCNeg #Lambda",200,-0.5,199.5);
 672
 673     fHistChi2PerNclsITSPosL[0]  = new TH1F("fHistChi2PerNclsITSPosL","chi2 per cluster ITS pi-p+ pos",250,0.0,50.0);
 674     fHistChi2PerNclsITSNegL[0]  = new TH1F("fHistChi2PerNclsITSNegL","chi2 per cluster ITS pi-p+ neg",250,0.0,50.0);
 675
 676     //------------for 2d pt dep studies ------------------------//
 677     fHistNclsITSPos[0]        = new TH2F("fHistNclsITSPos","fHistNclsITSPos  vs pt pos",nbMass,1.05,1.25,7,-0.5,6.5);
 678     fHistNclsITSNeg[0]        = new TH2F("fHistNclsITSNeg","fHistNclsITSNeg vs pt neg",nbMass,0.25,0.75,7,-0.5,6.5);
 679     fHistNclsITSPos[1]        = new TH2F("fHistNclsITSPosSec","fHistNclsITSPos  vs pt pos",nbMass,1.05,1.25,7,-0.5,6.5);
 680     fHistNclsITSNeg[1]        = new TH2F("fHistNclsITSNegSec","fHistNclsITSNeg vs pt neg",nbMass,0.25,0.75,7,-0.5,6.5);
 681
 682     fHistNclsTPCPos[0]        = new TH2F("fHistNclsTPCPos","K0s mass vs phi neg",nbMass,1.05,1.25,200,0.0,200.0);
 683     fHistNclsTPCNeg[0]        = new TH2F("fHistNclsTPCNeg","L mass vs phi neg",nbMass,0.25,0.75,200,0.0,200.0);
 684     fHistNclsTPCPos[1]        = new TH2F("fHistNclsTPCPosSec","K0s mass vs phi neg",nbMass,1.05,1.25,200,0.0,200.0);
 685     fHistNclsTPCNeg[1]        = new TH2F("fHistNclsTPCNegSec","L mass vs phi neg",nbMass,0.25,0.75,200,0.0,200.0);
 686
 687     fHistChi2PerNclsITSPos[0] = new TH2F("fHistChi2PerNclsITSPos","chi2 per cluster ITS K0s neg",nbMass,1.05,1.25,250,0.0,25.0);
 688     fHistChi2PerNclsITSNeg[0] = new TH2F("fHistChi2PerNclsITSNeg","chi2 per cluster ITS L neg",nbMass,0.25,0.75,250,0.0,25.0);
 689     fHistChi2PerNclsITSPos[1] = new TH2F("fHistChi2PerNclsITSPosSec","chi2 per cluster ITS K0s neg",nbMass,1.05,1.25,250,0.0,25.0);
 690     fHistChi2PerNclsITSNeg[1] = new TH2F("fHistChi2PerNclsITSNegSec","chi2 per cluster ITS L neg",nbMass,0.25,0.75,250,0.0,25.0);
 691
 692     fHistNCRowsTPCPos[0]      = new TH2F("fHistNCRowsTPCPos","n crossed rows vs K0s neg",nbMass,1.05,1.25,200,0.0,200.0);
 693     fHistNCRowsTPCNeg[0]      = new TH2F("fHistNCRowsTPCNeg","n crossed rows vs L neg",nbMass,0.25,0.75,200,0.0,200.0);
 694     fHistNCRowsTPCPos[1]      = new TH2F("fHistNCRowsTPCPosSec","n crossed rows vs K0s neg",nbMass,1.05,1.25,200,0.0,200.0);
 695     fHistNCRowsTPCNeg[1]      = new TH2F("fHistNCRowsTPCNegSec","n crossed rows vs L neg",nbMass,0.25,0.75,200,0.0,200.0);
 696
 697     fHistNclsITS[0]           = new TH2F("fHistNclsITS","fHistNclsITS pos vs neg L",10,-0.5,9.5,10,-0.5,9.5);
 698     fHistNclsITS[1]           = new TH2F("fHistNclsITSSec","fHistNclsITS pos vs neg K0",10,-0.5,9.5,10,-0.5,9.5);
 699
 700     fHistNclsTPC[0]           = new TH2F("fHistNclsTPC","ncls TPC neg vs crossed rows neg L",200,-0.5,199.5,200,-0.5,199.5);
 701     fHistNclsTPC[1]           = new TH2F("fHistNclsTPCSec","ncls TPC neg vs crossed rows neg K0",200,-0.5,199.5,200,-0.5,199.5);
 702
 703     fHistRatioFoundOverFinableTPCK0[0] = new TH2F("fHistRatioFoundOverFinableTPCK0","ncls found over findable K0s",nbMass,0.25,0.75,200,0.0,2.0);
 704     fHistRatioFoundOverFinableTPCK0[1] = new TH2F("fHistRatioFoundOverFinableTPCK0Sec","ncls found over findable K0s",nbMass,0.25,0.75,200,0.0,2.0);
 705
 706     fHistRatioFoundOverFinableTPCL[0]  = new TH2F("fHistRatioFoundOverFinableTPCL","ncls found over findable L neg",nbMass,1.05,1.25,200,0.0,2.0);
 707     fHistRatioFoundOverFinableTPCL[1]  = new TH2F("fHistRatioFoundOverFinableTPCLSec","ncls found over findable L neg",nbMass,1.05,1.25,200,0.0,2.0);
 708
 709     //--- add histos ----//
 710     fOutputContainer->Add(fHistNclsTPCNeg[0]);
 711     fOutputContainer->Add(fHistChi2PerNclsITSNeg[0]);
 712     fOutputContainer->Add(fHistNCRowsTPCNeg[0]);
 713     fOutputContainer->Add(fHistRatioFoundOverFinableTPCK0[0]);
 714     fOutputContainer->Add(fHistNclsITSPos[0]);
 715     fOutputContainer->Add(fHistNclsITSNeg[0]);
 716     fOutputContainer->Add(fHistNclsTPCPos[0]);
 717     fOutputContainer->Add(fHistChi2PerNclsITSPos[0]);
 718     fOutputContainer->Add(fHistNCRowsTPCPos[0]);
 719     fOutputContainer->Add(fHistRatioFoundOverFinableTPCL[0]);
 720     fOutputContainer->Add(fHistNclsITS[0]) ;
 721     fOutputContainer->Add(fHistNclsTPC[0]);
 722
 723     fOutputContainer->Add(fHistNclsTPCNeg[1]);
 724     fOutputContainer->Add(fHistChi2PerNclsITSNeg[1]);
 725     fOutputContainer->Add(fHistNCRowsTPCNeg[1]);
 726     fOutputContainer->Add(fHistRatioFoundOverFinableTPCK0[1]);
 727     fOutputContainer->Add(fHistNclsTPCPos[1]);
 728     fOutputContainer->Add(fHistChi2PerNclsITSPos[1]);
 729     fOutputContainer->Add(fHistNCRowsTPCPos[1]);
 730     fOutputContainer->Add(fHistRatioFoundOverFinableTPCL[1]);
 731     fOutputContainer->Add(fHistNclsITSPos[1]);
 732     fOutputContainer->Add(fHistNclsITSNeg[1]);
 733     fOutputContainer->Add(fHistNclsITS[1]);
 734     fOutputContainer->Add(fHistNclsTPC[1]);
 735
 736
 737     if(mchist==2){// for MC reco
 738
 739       //********************************************** K0 *************************************************//
 740
 741       //********************************************** Lambda ***********************************************//
 742       fHistPiPDCADaughters[1] = new TH2F("fHistPiPDCADaughtersSec","dca of #Lambda daughters",nbMass,1.05,1.25,250,0.0,2.0);
 743       fHistPiPDCADaughterPosToPrimVtxVSMass[1] = new TH2F("fHistPiPDCADaughterPosToPrimVtxVSMassSec","p DCA daughter to prim vtx vs invmass",nbMass,1.05,1.25,250,0.0,25.0);
 744       fHistPiPDCAVSMass[1] = new TH2F("fHistPiPDCAVSMassSec","ppi- dca  vs pt",nbMass,1.05,1.25,250,0.0,5.0);
 745       fHistPiPCosPointAng[1]  = new TH2F("fHistPiPCosPointAngSec","#Lambda cosine of pointing angle vs mass",nbMass,1.05,1.25,200,0.9,1.00);
 746
 747       fHistPiPRadiusXY[1] = new TH2F("fHistPiPRadiusXYSec","pi-p+ phi dist vs mass",nbMass,1.05,1.25,200,0.0,4.0);
 748
 749       //******************************************* Antilambda **********************************************//
 750       fHistPiAPDCADaughters[1] = new TH2F("fHistPiAPDCADaughtersSec","dca of #bar{#Lambda} daughters",nbMass,1.05,1.25,250,0.0,2.0);
 751       fHistPiAPDCADaughterPosToPrimVtxVSMass[1] = new TH2F("fHistPiAPDCADaughterPosToPrimVtxVSMassSec","pi+ DCA daughter to prim vtx vs invmass",nbMass,1.05,1.25,250,0.0,25.0);
 752       fHistPiAPDCAVSMass[1] = new TH2F("fHistPiAPDCAVSMassSec","pi+p- dca  vs pt",nbMass,1.05,1.25,250,0.0,5.0);
 753       fHistPiAPCosPointAng[1] = new TH2F("fHistPiAPCosPointAngSec","#bar{#Lambda} cosine of pointing angle vs mass",nbMass,1.05,1.25,200,0.9,1.00);
 754
 755       fHistPiAPRadiusXY[1] = new TH2F("fHistPiAPRadiusXYSec","pi+p- phi dist vs mass",nbMass,1.05,1.25,200,0.0,4.0);
 756
 757
 758       //******************************************* TPC ****************************************************//
 759
 760       fHistDedxSecProt[1] = new TH2F("fHistDedxSecProtSec","proton", nbPt, 0, 20, 100, 0, 400);
 761       fHistDedxSecPiPlus[1] = new TH2F("fHistDedxSecPiPlusSec","pi plus", nbPt, 0, 20, 100, 0, 400);
 762       fHistDedxSecAProt[1] = new TH2F("fHistDedxSecAProtSec","antiproton", nbPt, 0, 20, 100, 0, 400);
 763       fHistDedxSecPiMinus[1] = new TH2F("fHistDedxSecPiMinusSec","pi minus", nbPt, 0, 20, 100, 0, 400);
 764       fHistDedxProt[1] = new TH2F("fHistDedxProtSec","proton", nbPt, 0, 20, 100, 0, 400);
 765       fHistDedxPiPlus[1] = new TH2F("fHistDedxPiPlusSec","pi plus", nbPt, 0, 20, 100, 0, 400);
 766       fHistDedxAProt[1] = new TH2F("fHistDedxAProtSec","antiproton", nbPt, 0, 20, 100, 0, 400);
 767       fHistDedxPiMinus[1] = new TH2F("fHistDedxPiMinusSec","pi minus", nbPt, 0, 20, 100, 0, 400);
 768       fHistNclsITSPosK0[1] = new TH1F("fHistNclsITSPosK0Sec","fHistNclsITSPos K0 ",10,-0.5,9.5);
 769       fHistNclsITSNegK0[1] = new TH1F("fHistNclsITSNegK0Sec","fHistNclsITSNeg K0",10,-0.5,9.5);
 770       fHistNclsTPCPosK0[1] = new TH1F("fHistNclsTPCPosK0Sec","fHistNclsTPCPos K0",200,-0.5,199.5);
 771       fHistNclsTPCNegK0[1] = new TH1F("fHistNclsTPCNegK0Sec","fHistNclsTPCNeg K0",200,-0.5,199.5);
 772       fHistChi2PerNclsITSPosK0[1] = new TH1F("fHistChi2PerNclsITSPosK0Sec","chi2 per cluster ITS pi+pi- pos",250,0.0,50.0);
 773       fHistChi2PerNclsITSNegK0[1] = new TH1F("fHistChi2PerNclsITSNegK0Sec","chi2 per cluster ITS pi+pi- neg",250,0.0,50.0);
 774       fHistNclsITSPosL[1] = new TH1F("fHistNclsITSPosLSec","fHistNclsITSPos #Lambda",10,-0.5,9.5);
 775       fHistNclsITSNegL[1] = new TH1F("fHistNclsITSNegLSec","fHistNclsITSNeg #Lambda",10,-0.5,9.5);
 776       fHistNclsTPCPosL[1] = new TH1F("fHistNclsTPCPosLSec","fHistNclsTPCPos #Lambda",200,-0.5,199.5);
 777       fHistNclsTPCNegL[1] = new TH1F("fHistNclsTPCNegLSec","fHistNclsTPCNeg #Lambda",200,-0.5,199.5);
 778       fHistChi2PerNclsITSPosL[1] = new TH1F("fHistChi2PerNclsITSPosLSec","chi2 per cluster ITS pi-p+ pos",250,0.0,50.0);
 779       fHistChi2PerNclsITSNegL[1] = new TH1F("fHistChi2PerNclsITSNegLSec","chi2 per cluster ITS pi-p+ neg",250,0.0,50.0);
 780     }
 781
 782
 783     //--------------- add K0s histos --------------//
 784     fOutputContainer->Add(fHistPiPiDCADaughters[0]);
 785     fOutputContainer->Add(fHistPiPiDCADaughterPosToPrimVtxVSMass[0]);
 786     /*
 787       fOutputContainer->Add( fHistPiPiPtDaughters[0]);
 788
 789     */
 790     fOutputContainer->Add(fHistPiPiDCAVSMass[0]);
 791     fOutputContainer->Add(fHistPiPiCosPointAng[0]);
 792     /*
 793       fOutputContainer->Add(fHistPiPiRadiusXY[0]);
 794     */
 795     fOutputContainer->Add(fHistNclsITSPosK0[0]);
 796     fOutputContainer->Add(fHistNclsITSNegK0[0]);
 797     fOutputContainer->Add(fHistNclsTPCPosK0[0]);
 798     fOutputContainer->Add(fHistNclsTPCNegK0[0]);
 799     fOutputContainer->Add(fHistChi2PerNclsITSPosK0[0]);
 800     fOutputContainer->Add(fHistChi2PerNclsITSNegK0[0]);
 801
 802     //--------------- add (A)Lambda histos --------------//
 803     for(Int_t j=0;j<mchist;j++){
 804       fOutputContainer->Add(fHistPiPDCADaughters[j]);
 805       fOutputContainer->Add(fHistPiAPDCADaughters[j]);
 806
 807       fOutputContainer->Add( fHistPiPDCADaughterPosToPrimVtxVSMass[j]);
 808       fOutputContainer->Add( fHistPiAPDCADaughterPosToPrimVtxVSMass[j]);
 809
 810       /*
 811         fOutputContainer->Add( fHistPiPPtDaughters[j]);
 812         fOutputContainer->Add( fHistPiAPPtDaughters[j]);
 813       */
 814       fOutputContainer->Add(fHistPiPDCAVSMass[j]);
 815       fOutputContainer->Add(fHistPiAPDCAVSMass[j]);
 816
 817       fOutputContainer->Add(fHistPiPCosPointAng[j]);
 818       fOutputContainer->Add(fHistPiAPCosPointAng[j]);
 819       /*
 820         fOutputContainer->Add(fHistPiPRadiusXY[j]);
 821         fOutputContainer->Add(fHistPiAPRadiusXY[j]);
 822       */
 823
 824       //others
 825       fOutputContainer->Add(fHistDedxSecProt[j]);
 826       fOutputContainer->Add(fHistDedxSecAProt[j]);
 827       fOutputContainer->Add(fHistDedxSecPiPlus[j]);
 828       fOutputContainer->Add(fHistDedxSecPiMinus[j]);
 829
 830       fOutputContainer->Add(fHistDedxProt[j]);
 831       fOutputContainer->Add(fHistDedxAProt[j]);
 832       fOutputContainer->Add(fHistDedxPiPlus[j]);
 833       fOutputContainer->Add(fHistDedxPiMinus[j]);
 834
 835
 836       fOutputContainer->Add(fHistNclsITSPosL[j]);
 837       fOutputContainer->Add(fHistNclsITSNegL[j]);
 838
 839
 840       fOutputContainer->Add(fHistNclsTPCPosL[j]);
 841       fOutputContainer->Add(fHistNclsTPCNegL[j]);
 842
 843
 844       fOutputContainer->Add(fHistChi2PerNclsITSPosL[j]);
 845       fOutputContainer->Add(fHistChi2PerNclsITSNegL[j]);
 846     }
 847
 848   }
 849
 850   //----------------------------- MC reco or MC truth only --------------------------//
 851   if((fMCMode && fMCTruthMode) || fMCTruthMode){
 852     if(fAnapp){
 853       fHistPrimVtxZESDVSNContributorsMC = new TH2F("fHistPrimVtxZESDVSNContributorsMC","prim vtx pos z ESD vs no. of contributers MC",250,-50,50,500,0.0,500.0);
 854       fOutputContainer->Add(fHistPrimVtxZESDVSNContributorsMC);
 855       fHistPrimVtxZESDTPCVSNContributorsMC = new TH2F("fHistPrimVtxZESDTPCVSNContributorsMC","prim vtx pos z TPC vs no. of contributers MC",250,-50,50,500,0.0,500.0);
 856       fOutputContainer->Add(fHistPrimVtxZESDTPCVSNContributorsMC);
 857       fHistPrimVtxZESDSPDVSNContributorsMC = new TH2F("fHistPrimVtxZESDSPDVSNContributorsMC","prim vtx pos z SPD vs no. of contributers MC",250,-50,50,500,0.0,500.0);
 858       fOutputContainer->Add(fHistPrimVtxZESDSPDVSNContributorsMC);
 859     }
 860
 861     fHistMCVertexZ= new TH1F("fHistMCVertexZ"," z vertex distr in cm MC",500,-50,50);
 862     fOutputContainer->Add(fHistMCVertexZ);
 863     fHistPiPiPDGCode = new TH1F("fHistPiPiPDGCode","PDG code of K0s mothers",3500,0,3500);
 864     fOutputContainer->Add(fHistPiPiPDGCode);
 865     fHistPiPPDGCode = new TH1F("fHistPiPPDGCode","PDG code of #Lambda  mothers",3500,0,3500);
 866     fOutputContainer->Add(fHistPiPPDGCode);
 867     fHistPiAPPDGCode = new TH1F("fHistPiAPPDGCode","PDG code of #bar{#Lambda} mothers",3500,0,3500);
 868     fOutputContainer->Add(fHistPiAPPDGCode);
 869     fHistPiPCosPointAngXiVsPt= new TH2F("fHistPiPCosPointAngXiVsPt","pi-p cos of pointing angle vs pt from xi",200,0.0,20.0,250,0.99,1.00);
 870     fOutputContainer->Add(fHistPiPCosPointAngXiVsPt);
 871     fHistPiAPCosPointAngXiVsPt= new TH2F("fHistPiAPCosPointAngXiVsPt","pi+p- cos of pointing angle vs pt from xi",200,0.0,20.0,250,0.98,1.00);
 872     fOutputContainer->Add(fHistPiAPCosPointAngXiVsPt);
 873     fHistPiPiEtaDMC[0] = new TH2F("fHistPiPiEtaDMCRaw","K0s daughters etaMC raw",300,-6,6,100,0,20);//
 874     fOutputContainer->Add(fHistPiPiEtaDMC[0]);
 875     fHistPiPiEtaDMC[1] = new TH2F("fHistPiPiEtaDMC","K0s daughters etaMC after rap V0 cut",300,-6,6,100,0,20);
 876     fOutputContainer->Add(fHistPiPiEtaDMC[1]);
 877     fHistPiPEtaDMC[0] = new TH2F("fHistPiPEtaDMCRaw","#Lambda daughters etaMC raw",300,-6,6,100,0,20);
 878     fOutputContainer->Add(fHistPiPEtaDMC[0]);
 879     fHistPiPEtaDMC[1] = new TH2F("fHistPiPEtaDMC","#Lambda daughters etaMC after rap V0 cut",300,-6,6,100,0,20);
 880     fOutputContainer->Add(fHistPiPEtaDMC[1]);
 881
 882     //********************************************** K0 *************************************************//
 883
 884     fHistPiPiDecayLengthResolution[0] = new TH2F("fHistPiPiDecayLengthResolution","K0s decay length resolution MC",220,0.0,110.0,220,0.0,110);
 885
 886     //********************************************** Lambda **********************************************//
 887     fHistPiPDecayLengthResolution[0] = new TH2F("fHistPiPDecayLengthResolution","Lambda decay length resolution MC",220,0.0,110.0,220,0.0,110);
 888     fHistPiPDecayLengthResolution[1] = new TH2F("fHistPiPDecayLengthResolutionSec","Lambda sec decay length resolution MC",220,0.0,110.0,220,0.0,110);
 889
 890     fHistPiPMassVSPtSecSigma[0] = new TH2F("fHistPiPMassVSPtSecSigmaMC"," pi-p+ InvMass distribution secondaries from sigma MC",nbMass,1.05,1.25,200,0.,20);
 891     fHistPiPMassVSPtSecSigma[1] = new TH2F("fHistPiPMassVSPtSecSigma"," pi-p+ InvMass distribution secondaries from Sigma reco",nbMass,1.05,1.25,200,0.,20);
 892
 893     fHistPiPMassVSPtSecXi[0] = new TH2F("fHistPiPMassVSPtSecXiMC"," pi-p+ InvMass distribution secondaries from  xi MC",nbMass,1.05,1.25,200,0.,20);
 894     fHistPiPMassVSPtSecXi[1] = new TH2F("fHistPiPMassVSPtSecXi"," pi-p+ InvMass distribution secondaries from  xi  reco",nbMass,1.05,1.25,200,0.,20);
 895
 896     fHistPiPMassVSPtSecXiMCTruth = new TH2F("fHistPiPMassVSPtSecXiMCTruth","Lambda mass reco vs pt sec Lambda from xi MC truth pt",nbMass,1.05,1.25,200,0.0,20.0);
 897
 898     fHistPiPMassVSYSecXi[0] = new TH2F("fHistPiPMassVSYSecXiMC"," pi-p+ InvMass distribution secondaries from xi MC",nbMass,1.05,1.25,100,-2.,2);
 899     fHistPiPMassVSYSecXi[1] = new TH2F("fHistPiPMassVSYSecXi"," pi-p+ InvMass distribution secondaries from xi reco",nbMass,1.05,1.25,100,-2.,2);
 900
 901     fHistPiPXi0PtVSLambdaPt[0]= new TH2F("fHistPiPXi0PtVSLambdaPtMC"," pt xi 0 vs pt lambda MC truth",200,0.0,20.0,200,0.0,20.0);
 902     fHistPiPXi0PtVSLambdaPt[1]= new TH2F("fHistPiPXi0PtVSLambdaPt"," pt xi 0 truth vs pt lambda reco",200,0.0,20.0,200,0.0,20.0);
 903
 904     fHistPiPXiMinusPtVSLambdaPt[0]= new TH2F("fHistPiPXiMinusPtVSLambdaPtMC","pt xi- vs pt lambda MC truth",200,0.0,20.0,200,0.0,20.0);
 905     fHistPiPXiMinusPtVSLambdaPt[1]= new TH2F("fHistPiPXiMinusPtVSLambdaPt","pt xi- truth vs pt lambda reco",200,0.0,20.0,200,0.0,20.0);
 906
 907     fHistPiPOmegaPtVSLambdaPt[0] = new TH2F("fHistPiPOmegaPtVSLambdaPtMC","pt omega vs pt lambda MC truth",200,0.0,20.0,200,0.0,20.0);
 908     fHistPiPOmegaPtVSLambdaPt[1] = new TH2F("fHistPiPOmegaPtVSLambdaPt","pt omega vs pt lambda MC reco",200,0.0,20.0,200,0.0,20.0);
 909
 910     fHistPiPMassVSPtSecOmega[0] = new TH2F("fHistPiPMassVSPtSecOmegaMC","Lambda mass vs pt omega MCtruth",nbMass,1.05,1.25,200,0.0,20.0);
 911     fHistPiPMassVSPtSecOmega[1] = new TH2F("fHistPiPMassVSPtSecOmega","Lambda mass vs pt omega MCreco",nbMass,1.05,1.25,200,0.0,20.0);
 912     fHistPiPMassVSPtSecOmegaMCTruth= new TH2F("fHistPiPMassVSPtSecOmegaMCTruth","Lambda mass vs pt sec Lambda from Omega MC truth pt",nbMass,1.05,1.25,200,0.0,20.0);
 913
 914     //******************************************* Antilambda **********************************************//
 915     fHistPiAPDecayLengthResolution[0] = new TH2F("fHistPiAPDecayLengthResolution","ALambda decay length resolution MC",220,0.0,110.0,220,0.0,110);
 916     fHistPiAPDecayLengthResolution[1] = new TH2F("fHistPiAPDecayLengthResolutionSec","ALambda sec decay length resolution MC",220,0.0,110.0,220,0.0,110);
 917
 918     fHistPiAPMassVSPtSecSigma[0] = new TH2F("fHistPiAPMassVSPtSecSigmaMC"," pi+p- InvMass distribution secondaries from Sigma MC",nbMass,1.05,1.25,200,0.,20);
 919     fHistPiAPMassVSPtSecSigma[1] = new TH2F("fHistPiAPMassVSPtSecSigma"," pi+p- InvMass distribution secondaries from  Sigma  reco",nbMass,1.05,1.25,200,0.,20);
 920
 921     fHistPiAPMassVSPtSecXi[0] = new TH2F("fHistPiAPMassVSPtSecXiMC"," pi+p- InvMass distribution secondaries from xi MC",nbMass,1.05,1.25,200,0.,20);
 922     fHistPiAPMassVSPtSecXi[1] = new TH2F("fHistPiAPMassVSPtSecXi"," pi+p- InvMass distribution secondaries from  Xi reco",nbMass,1.05,1.25,200,0.,20);
 923
 924     fHistPiAPMassVSPtSecXiMCTruth = new TH2F("fHistPiAPMassVSPtSecXiMCTruth","ALambda mass reco vs pt sec Lambda from xi MC truth pt",nbMass,1.05,1.25,200,0.0,20.0);
 925
 926     fHistPiAPMassVSYSecXi[0] = new TH2F("fHistPiAPMassVSYSecXiMC"," pi+p- InvMass distribution secondaries from  xi MC",nbMass,1.05,1.25,100,-2,2);
 927     fHistPiAPMassVSYSecXi[1] = new TH2F("fHistPiAPMassVSYSecXi"," pi+p- InvMass distribution secondaries from xi reco",nbMass,1.05,1.25,100,-2.,2);
 928
 929
 930     fHistPiAPXi0PtVSLambdaPt[0]= new TH2F("fHistPiAPXi0PtVSLambdaPtMC"," pt xi 0 vs pt Alambda MC truth",200,0.0,20.0,200,0.0,20.0);
 931     fHistPiAPXi0PtVSLambdaPt[1]= new TH2F("fHistPiAPXi0PtVSLambdaPt"," pt xi 0 truth vs pt Alambda reco",200,0.0,20.0,200,0.0,20.0);
 932
 933     fHistPiAPXiMinusPtVSLambdaPt[0]= new TH2F("fHistPiAPXiMinusPtVSLambdaPtMC","pt xi- vs pt Alambda MC truth",200,0.0,20.0,200,0.0,20.0);
 934     fHistPiAPXiMinusPtVSLambdaPt[1]= new TH2F("fHistPiAPXiMinusPtVSLambdaPt","pt xi- truth vs pt Alambda reco",200,0.0,20.0,200,0.0,20.0);
 935
 936     fHistPiAPOmegaPtVSLambdaPt[0] = new TH2F("fHistPiAPOmegaPtVSLambdaPtMC","pt omega vs pt alambda MC truth",200,0.0,20.0,200,0.0,20.0);
 937     fHistPiAPOmegaPtVSLambdaPt[1] = new TH2F("fHistPiAPOmegaPtVSLambdaPt","pt omega vs pt alambda MC reco",200,0.0,20.0,200,0.0,20.0);
 938
 939     fHistPiAPMassVSPtSecOmega[0] = new TH2F("fHistPiAPMassVSPtSecOmegaMC","ALambda mass vs pt omega MCtruth",nbMass,1.05,1.25,200,0.0,20.0);
 940     fHistPiAPMassVSPtSecOmega[1] = new TH2F("fHistPiAPMassVSPtSecOmega","ALambda mass vs pt omega MCreco",nbMass,1.05,1.25,200,0.0,20.0);
 941     fHistPiAPMassVSPtSecOmegaMCTruth= new TH2F("fHistPiAPMassVSPtSecOmegaMCTruth","ALambda mass vs pt sec Lambda from Omega MC truth pt",nbMass,1.05,1.25,200,0.0,20.0);
 942
 943     fOutputContainer->Add(fHistPiPMassVSPtSecXiMCTruth);
 944     fOutputContainer->Add(fHistPiPMassVSPtSecOmegaMCTruth);
 945
 946     fOutputContainer->Add(fHistPiAPMassVSPtSecXiMCTruth);
 947     fOutputContainer->Add(fHistPiAPMassVSPtSecOmegaMCTruth);
 948
 949     //--- add K0s histo ------//
 950     fOutputContainer->Add(fHistPiPiDecayLengthResolution[0]);
 951
 952     //--- add (A)Lambda histo ------//
 953     for(Int_t j=0;j<2;j++){
 954
 955       fOutputContainer->Add(fHistPiPDecayLengthResolution[j]);
 956       fOutputContainer->Add(fHistPiAPDecayLengthResolution[j]);
 957
 958       fOutputContainer->Add(fHistPiPMassVSPtSecXi[j]);
 959       fOutputContainer->Add(fHistPiAPMassVSPtSecXi[j]);
 960
 961       fOutputContainer->Add(fHistPiPMassVSYSecXi[j]);
 962       fOutputContainer->Add(fHistPiAPMassVSYSecXi[j]);
 963
 964       fOutputContainer->Add(fHistPiPXi0PtVSLambdaPt[j]);
 965       fOutputContainer->Add(fHistPiAPXi0PtVSLambdaPt[j]);
 966
 967       fOutputContainer->Add(fHistPiPXiMinusPtVSLambdaPt[j]);
 968       fOutputContainer->Add(fHistPiAPXiMinusPtVSLambdaPt[j]);
 969
 970       fOutputContainer->Add(fHistPiPMassVSPtSecSigma[j]);
 971       fOutputContainer->Add(fHistPiAPMassVSPtSecSigma[j]);
 972
 973       fOutputContainer->Add(fHistPiPOmegaPtVSLambdaPt[j]);
 974       fOutputContainer->Add(fHistPiAPOmegaPtVSLambdaPt[j]);
 975       fOutputContainer->Add(fHistPiPMassVSPtSecOmega[j]);
 976       fOutputContainer->Add(fHistPiAPMassVSPtSecOmega[j]);
 977     }
 978   }
 979
 980   /*
 981   //shift q/pt
 982   fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5);
 983   */
 984
 985 }
 986
 987 //________________________________________________________________________
 988 void AliAnalysisTaskV0ForRAA::UserExec(Option_t *) {
 989   //user exec
 990
 991   //-- esd handler --//
 992   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
 993     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 994   if (!esdH) {
 995     printf("ERROR: Could not get ESDInputHandler");
 996     return;
 997   }
 998   fESD = esdH->GetEvent();
 999   if(!fESD) {
1000     printf("ERROR: fESD not available \n");
1001     return ;
1002   }
1003
1004   //-- mc handler --//
1005   if(fMCMode || fMCTruthMode){
1006     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>
1007       (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1008     if(!mcH) {
1009       printf("ERROR: Could not get MCInputHandler");
1010       return;
1011     }
1012     fMCev = mcH->MCEvent();
1013     if (!fMCev) {
1014       printf("ERROR: fMCev not available \n");
1015       return ;
1016     }
1017   }
1018
1019   //--  AliPIDResponse --//
1020   fESDpid = esdH->GetPIDResponse();
1021
1022   //-- Count events before cuts --//
1023   fHistNEvents->Fill(0);
1024
1025   //-- Check object existence --//
1026   const AliESDVertex *    vtxESD    = fESD->GetPrimaryVertexTracks();
1027   const AliESDVertex *    vtxESDTPC = fESD->GetPrimaryVertexTPC();
1028   const AliESDVertex *    vtxESDSPD = fESD->GetPrimaryVertexSPD();
1029   const AliMultiplicity * multESD   = fESD->GetMultiplicity();
1030
1031   if ( !vtxESD ){
1032     AliError("No Tracks Vertex");
1033     return;
1034   }
1035
1036   if ( !vtxESDTPC ){
1037     AliError("No TPC Vertex");
1038     return ;
1039   }
1040
1041   if ( !vtxESDSPD ){
1042     AliError("No SPD Vertex");
1043     return ;
1044   }
1045
1046   if ( !multESD ){
1047     AliError("No Multiplicity");
1048     return ;
1049   }
1050
1051
1052   // ----------- MC vertex -----------------------------------//
1053
1054   Int_t nContr =0;
1055
1056   if(fMCTruthMode){
1057     Double_t vVertexPrim[3];
1058     fMCev->GetPrimaryVertex()->GetXYZ(vVertexPrim);
1059     fHistMCVertexZ->Fill(vVertexPrim[2]);
1060
1061     if(fMCMode && fAnapp){
1062       if (vtxESD->GetStatus()){
1063         nContr=vtxESD->GetNContributors();
1064         fHistPrimVtxZESDVSNContributorsMC->Fill(vVertexPrim[2],nContr);
1065         fHistPrimVtxZESDTPCVSNContributorsMC->Fill(vVertexPrim[2],nContr);
1066       }
1067       else {
1068         if(vtxESDSPD->GetStatus()){
1069           nContr=vtxESDSPD->GetNContributors();
1070           fHistPrimVtxZESDTPCVSNContributorsMC->Fill(vVertexPrim[2],nContr);
1071           fHistPrimVtxZESDSPDVSNContributorsMC->Fill(vVertexPrim[2],nContr);
1072         }
1073         else{
1074           fHistPrimVtxZESDVSNContributorsMC->Fill(vVertexPrim[2],nContr);//add for correction ESD and ESDPSD!!!!
1075           fHistPrimVtxZESDTPCVSNContributorsMC->Fill(vVertexPrim[2],nContr);
1076         }
1077       }
1078     }
1079   }
1080
1081
1082
1083   // -- Check fo centrality
1084   Bool_t process = kTRUE;
1085   Int_t centBin = -1;
1086   if(fUseCentrality) {
1087     centBin = CalculateCentralityBin();
1088     if(!fUseCentralityRange){
1089       if(centBin!= fUseCentralityBin) process=kFALSE;
1090     }
1091     else if(centBin < fUseCentralityBin || centBin > fUseCentralityBin+fUseCentralityRange)
1092       process = kFALSE;
1093   }
1094
1095   AliESDVZERO* esdV0 = fESD->GetVZEROData();
1096   Float_t multV0 = esdV0->GetMTotV0A() + esdV0->GetMTotV0C();
1097
1098   if(fAnapp){// pp Analysis
1099     // SDD test for 2.76TeV pp
1100     // select events with SDD
1101     //   TString trCl = fESD->GetFiredTriggerClasses();
1102     //if(!(trCl.Contains("ALLNOTRD")) && fSelSDD) return;
1103     UInt_t maskSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1104     if(maskSel& AliVEvent::kFastOnly && fSelSDD) return;
1105     if(!(maskSel& AliVEvent::kFastOnly) && fSelNoSDD) return;
1106
1107
1108     //-- Monitor event cuts --//
1109     fHistNEvents->Fill(1);
1110
1111     Int_t ntracks = fESD->GetNumberOfTracks();
1112     for(Int_t i=0;i<ntracks;i++){//check sdd event selection
1113       AliESDtrack *tr=   fESD->GetTrack(i);
1114
1115       Bool_t sdd0 = tr->HasPointOnITSLayer(0);
1116       Bool_t sdd1 = tr->HasPointOnITSLayer(1);
1117       Bool_t sdd2 = tr->HasPointOnITSLayer(2);
1118       Bool_t sdd3 = tr->HasPointOnITSLayer(3);
1119       Bool_t sdd4 = tr->HasPointOnITSLayer(4);
1120       Bool_t sdd5 = tr->HasPointOnITSLayer(5);
1121
1122       fHistITSLayerHits->Fill(Int_t(sdd0)*(-1),ntracks);
1123       fHistITSLayerHits->Fill(Int_t(sdd1)*1,ntracks);
1124       fHistITSLayerHits->Fill(Int_t(sdd2)*2,ntracks);
1125       fHistITSLayerHits->Fill(Int_t(sdd3)*3,ntracks);
1126       fHistITSLayerHits->Fill(Int_t(sdd4)*4,ntracks);
1127       fHistITSLayerHits->Fill(Int_t(sdd5)*5,ntracks);
1128     }
1129
1130     //vertex selection
1131     if (vtxESD->GetStatus()){
1132       fHistNEvents->Fill(2);
1133       fHistESDVertexZ->Fill(vtxESD->GetZv());
1134       if(fabs(vtxESD->GetZv()) < fVertexZCut){
1135         fHistMuliplicityRaw->Fill(multV0);
1136         fHistNEvents->Fill(3);
1137         fHistNPrim->Fill(nContr);
1138
1139         Process();
1140
1141         fHistMuliplicity->Fill(multV0);
1142
1143         nContr = vtxESD->GetNContributors();
1144
1145         fHistPrimVtxZESDVSNContributors->Fill(vtxESD->GetZv(),nContr);
1146         fHistPrimVtxZESDTPCVSNContributors->Fill(vtxESDTPC->GetZv(),nContr);
1147
1148         fHistPrimVtxZESD->Fill(vtxESD->GetZv());
1149         fHistPrimVtxZESDTPC->Fill(vtxESDTPC->GetZv());
1150
1151         // -- count events after processing
1152         fHistNEvents->Fill(4);
1153       }
1154     }
1155     else{
1156       if(vtxESDSPD->GetStatus()){
1157         fHistNEvents->Fill(2);
1158
1159         fHistESDVertexZ->Fill(vtxESDSPD->GetZv());
1160         if(fabs(vtxESDSPD->GetZv()) < fVertexZCut){
1161
1162           fHistMuliplicityRaw->Fill(multV0);
1163           fHistNEvents->Fill(3);
1164           fHistNPrim->Fill(nContr);
1165
1166           Process();
1167
1168           fHistMuliplicity->Fill(multV0);
1169
1170           nContr = vtxESDSPD->GetNContributors();
1171
1172           fHistPrimVtxZESDTPCVSNContributors->Fill(vtxESDTPC->GetZv(),nContr);
1173           fHistPrimVtxZESDSPDVSNContributors->Fill(vtxESDSPD->GetZv(),nContr);
1174
1175           fHistPrimVtxZESDTPC->Fill(vtxESDTPC->GetZv());
1176           fHistPrimVtxZESDSPD->Fill(vtxESDSPD->GetZv());
1177           // -- count events after processing
1178           fHistNEvents->Fill(4);
1179         }
1180       }
1181       //else return;
1182     }
1183   }
1184   else{// PbPb analysis
1185     //-- Monitor event cuts --//
1186     fHistNEvents->Fill(1);
1187
1188     if(vtxESD->GetStatus()){
1189       Double_t vtxZ = vtxESD->GetZv();
1190       fHistESDVertexZ->Fill(vtxZ);
1191       if(process){
1192         fHistNEvents->Fill(2);
1193         if(fabs(vtxZ) < fVertexZCut){
1194           nContr = vtxESD->GetNContributors();
1195           fHistMuliplicityRaw->Fill(multV0);
1196           fHistNEvents->Fill(3);
1197           fHistNPrim->Fill(nContr);
1198           Process();
1199           fHistMuliplicity->Fill(multV0);
1200           fHistPrimVtxZESD->Fill(vtxZ);
1201           fHistPrimVtxZESDVSNContributors->Fill(vtxZ,nContr);
1202           // -- count events after processing --//
1203           fHistCentBin->Fill(centBin);
1204           fHistNEvents->Fill(4);
1205         }
1206       }
1207       if(fabs(vtxZ) < fVertexZCut) fHistCentBinRaw->Fill(centBin);
1208     }
1209   }
1210   PostData(1,fOutputContainer);
1211
1212 }
1213
1214 //________________________________________________________________________
1215 void AliAnalysisTaskV0ForRAA::Terminate(Option_t *) {
1216   //terminate
1217 }
1218 //________________________________________________________________________
1219 void AliAnalysisTaskV0ForRAA::Process(){
1220   //run the analysis
1221
1222   Int_t ntracks = fESD->GetNumberOfTracks();
1223   Int_t count = 0;
1224
1225   //-- count number of tracks --//
1226
1227   if(!(!fMCMode && fMCTruthMode)){
1228     for(Int_t i=0;i<ntracks;i++){
1229       AliESDtrack *track = (AliESDtrack*)fESD->GetTrack(i);
1230       if(!fESDTrackCuts->AcceptTrack(track)) continue;
1231       if( track->Eta() > fEtaCutMCDaughtersVal) continue;
1232       count++;
1233     }
1234     fHistMultiplicityPrimary->Fill(count);
1235   }
1236
1237   //-- check number of V0s in case of data or mc data like analysis--//
1238   Int_t nV0 = fESD->GetNumberOfV0s();
1239   if(!fMCTruthMode) if(nV0 < 1) return;
1240
1241   //-- run analysis --//
1242   if(fMCTruthMode)  V0MCTruthLoop();
1243   else  V0RecoLoop(0,0,0,0,0.0,0,0.0,0.0);
1244
1245 }
1246 //________________________________________________________________________
1247 void AliAnalysisTaskV0ForRAA::V0MCTruthLoop(){
1248   //loop over MC truth particles
1249
1250   //-- get MC stack --//
1251   AliStack *stack = fMCev->Stack();
1252
1253   /*
1254   //histo fo user defined shift in charge/pt
1255   if(fShift){
1256   fHistUserPtShift->Fill(fDeltaInvP);
1257   }
1258   */
1259   /*
1260     AliKFVertex primVtxStart(*(fESD->GetPrimaryVertex()));
1261     Int_t nTracksPrim=primVtxStart.GetNContributors();
1262     fHistNPrim->Fill(nTracksPrim);
1263   */
1264   /*
1265   // MC
1266
1267   Int_t mcPrimaries = stack->GetNprimary();
1268   Int_t mcParticles    = stack->GetNtrack();
1269
1270   fHistMultiplicityPrimary->Fill(mcPrimaries);
1271   fHistMCMultiplicityTracks->Fill(mcParticles);
1272
1273   // number of V0
1274   fHistNV0->Fill(nV0);
1275   if(nTracksPrim>0) {
1276   fHistNV0WithVertex->Fill(nV0);
1277   }
1278   */
1279
1280   //-- MC truht loop for V0s --//
1281   for (Int_t iMc = 0; iMc < (stack->GetNtrack()); iMc++){//MC truth loop
1282     Int_t fillMCtruth= int(fMCTruthMode);
1283     if(fMCTruthMode){
1284       fHistPiPiMonitorMCCuts[0]->Fill(1*fillMCtruth);
1285       fHistPiPMonitorMCCuts[0]->Fill(1*fillMCtruth);
1286       fHistPiAPMonitorMCCuts[0]->Fill(1*fillMCtruth);
1287     }
1288     TParticle *p0 = stack->Particle(iMc);
1289     if(!p0) continue;
1290
1291     if(fMCTruthMode){
1292       fHistPiPiMonitorMCCuts[0]->Fill(2*fillMCtruth);
1293       fHistPiPMonitorMCCuts[0]->Fill(2*fillMCtruth);
1294       fHistPiAPMonitorMCCuts[0]->Fill(2*fillMCtruth);
1295     }
1296
1297
1298
1299     Int_t pdgCode = p0->GetPdgCode();
1300
1301     //-------------- only K0s and Lambda ----------//
1302     if( (pdgCode != 310 ) && ( fabs(pdgCode) != 3122 ) ) continue; //accept only K0s and (A)Lambda
1303     Int_t fillFlagK0 = (3122- fabs(pdgCode))/(3122-310)*fillMCtruth;//K0s flag
1304     Int_t fillFlagL  = (fabs(pdgCode) - 310)/(3122-310)*(pdgCode+3122)/(2*3122)*fillMCtruth;// Lambda flag
1305     Int_t fillFlagAL = (fabs(pdgCode) - 310)/(3122-310)*(pdgCode-3122)/(-2*3122)*fillMCtruth; // AntiLambda flag
1306
1307     fHistPiPiMonitorMCCuts[0]->Fill(3*fillFlagK0);
1308     fHistPiPMonitorMCCuts[0]->Fill(3*fillFlagL);
1309     fHistPiAPMonitorMCCuts[0]->Fill(3*fillFlagAL);
1310
1311     if(p0->GetNDaughters() !=2) continue; // 2 body decay
1312     fHistPiPiMonitorMCCuts[0]->Fill(4*fillFlagK0);
1313     fHistPiPMonitorMCCuts[0]->Fill(4*fillFlagL);
1314     fHistPiAPMonitorMCCuts[0]->Fill(4*fillFlagAL);
1315
1316     //-------------- unique ID check-------------- //
1317     Int_t uniqueID =  p0->GetUniqueID();
1318     if(uniqueID==13) continue; //reject part. from material
1319
1320     fHistPiPiMonitorMCCuts[0]->Fill(5*fillFlagK0);
1321     fHistPiPMonitorMCCuts[0]->Fill(5*fillFlagL);
1322     fHistPiAPMonitorMCCuts[0]->Fill(5*fillFlagAL);
1323
1324     //-------------- daughters --------------------//
1325     Int_t id0  = p0->GetDaughter(0);
1326     Int_t id1  = p0->GetDaughter(1);
1327     if(id0<0 || id1 <0) continue;//check if daughters exists in stack
1328
1329     fHistPiPiMonitorMCCuts[0]->Fill(6*fillFlagK0);
1330     fHistPiPMonitorMCCuts[0]->Fill(6*fillFlagL);
1331     fHistPiAPMonitorMCCuts[0]->Fill(6*fillFlagAL);
1332
1333     Int_t pdgCodeD0 = stack->Particle(id0)->GetPdgCode();
1334     Int_t pdgCodeD1 = stack->Particle(id1)->GetPdgCode();
1335
1336     if(pdgCodeD0 == pdgCodeD1) continue;//pdg code shall be different
1337     if(pdgCodeD0*pdgCodeD1>0) continue;//and not of the same sign
1338
1339     fHistPiPiMonitorMCCuts[0]->Fill(7*fillFlagK0);
1340     fHistPiPMonitorMCCuts[0]->Fill(7*fillFlagL);
1341     fHistPiAPMonitorMCCuts[0]->Fill(7*fillFlagAL);
1342
1343     if((fabs(pdgCodeD0) != 211 ) && ( fabs(pdgCodeD0) != 2212 )) continue;//accept only pions and protons as daughters
1344     if((fabs(pdgCodeD1) != 211 ) && ( fabs(pdgCodeD1) != 2212 )) continue;//accept only pions and protons as daughters
1345
1346     fHistPiPiMonitorMCCuts[0]->Fill(8*fillFlagK0);
1347     fHistPiPMonitorMCCuts[0]->Fill(8*fillFlagL);
1348     fHistPiAPMonitorMCCuts[0]->Fill(8*fillFlagAL);
1349
1350     TParticle *p00 =stack->Particle(id0);
1351     TParticle *p01 =stack->Particle(id1);
1352     Double_t etaMC00   = p00->Eta();
1353     Double_t etaMC01   = p01->Eta();
1354
1355     //----------- unique ID check daughters-------- //
1356     Int_t uniqueIDdaughter0 = p00->GetUniqueID();
1357     Int_t uniqueIDdaughter1 = p01->GetUniqueID();
1358     if (uniqueIDdaughter0 !=4 || uniqueIDdaughter1 !=4 ) continue;// daughters shall come from decay
1359
1360     fHistPiPiMonitorMCCuts[0]->Fill(9*fillFlagK0);
1361     fHistPiPMonitorMCCuts[0]->Fill(9*fillFlagL);
1362     fHistPiAPMonitorMCCuts[0]->Fill(9*fillFlagAL);
1363
1364     fHistPiPiMonitorMCCuts[1]->Fill(9*fillFlagK0);
1365     fHistPiPMonitorMCCuts[1]->Fill(9*fillFlagL);
1366     fHistPiAPMonitorMCCuts[1]->Fill(9*fillFlagAL);
1367
1368     //------------ check label reco -------------------//
1369     if(fCheckNegLabelReco || fOnlyFoundRecoV0){//optional, set flag in AddTask
1370       Bool_t found =kFALSE;
1371       Int_t label0=0,label1=0;
1372       AliESDv0 * v0MIsMC=NULL;
1373       AliESDtrack *tr0 = NULL;
1374       AliESDtrack *tr1 = NULL;
1375       for(Int_t recL=0;recL < fESD->GetNumberOfV0s();recL++){
1376         v0MIsMC = fESD->GetV0(recL);
1377         if(!v0MIsMC) continue;
1378         tr0 = fESD->GetTrack(v0MIsMC->GetPindex());
1379         tr1 = fESD->GetTrack(v0MIsMC->GetNindex());
1380         if(tr0 && tr1){
1381           label0 = tr0->GetLabel();
1382           label1 = tr1->GetLabel();
1383           if((fabs(label0) == id0 && fabs(label1) == id1) ||
1384              (fabs(label0) == id1 && fabs(label1) == id0)){
1385             found =kTRUE;
1386             break;
1387           }
1388         }
1389       }
1390       if(fCheckNegLabelReco && !fOnlyFoundRecoV0) {
1391         if(found && (label0 <0 || label1 < 0)) continue;
1392       }
1393       else{
1394         if(!found) continue;
1395         if(fCheckNegLabelReco && found && (label0 <0 || label1 < 0)) continue;
1396       }
1397
1398     }
1399     //-----------get geometric properties --------------//
1400     // DCA of mother to prim vertex = production vertex
1401
1402     //-- primary and secondary vetex --//
1403     Double_t vVertexPrimMC[3];
1404     fMCev->GetPrimaryVertex()->GetXYZ(vVertexPrimMC);
1405     // Double_t x0=p0->Vx(),y0=p0->Vy(),z0=p0->Vz();//mother production vertex
1406
1407     Double_t x=p00->Vx(),y=p00->Vy(),z=p00->Vz();//daughter vertex =V0 decay vertex
1408     Double_t rx = x - vVertexPrimMC[0];
1409     Double_t ry = y - vVertexPrimMC[1];
1410     Double_t rz = z - vVertexPrimMC[2];
1411     Double_t sdeclength = rx*rx+ry*ry;//+rz*rz;//=p00->Rho();
1412     Double_t declength =0.0;
1413     if(sdeclength>0) declength = sqrt(sdeclength);
1414     Double_t declength3d = sqrt( rx*rx+ry*ry+rz*rz);
1415
1416     //-- decay radii --//
1417     Double_t rMC2D  = sqrt(x*x+y*y);
1418     const  Double_t xyzMC[3] = {x,y,z};
1419     // Double_t rMC = p00->R();
1420
1421     //-- phi --//
1422     Double_t pi = TMath::Pi();
1423     Double_t phi = p0->Phi();
1424     if(phi>pi) phi -=2*pi;
1425
1426     //-------------------- V0 variables ----------------//
1427     Double_t rapidity = p0->Y();
1428     Double_t massV0MC = p0->GetMass();
1429     Double_t ptV0MC =  p0->Pt();
1430     Double_t pV0MC =  p0->P();
1431
1432
1433     //----------------- mother variables-----------------//
1434     Int_t indexMother1  = p0->GetMother(0);
1435     Int_t isSecd=0; //primary:0, secondary:1
1436     Int_t pdgMother =0;
1437     Int_t uniqueIDmother=0;
1438     Double_t ptXiMother=0.0;
1439     Double_t rapXiMother = 0.0;
1440
1441
1442     //------check mother and fill mother histos---------//
1443     Bool_t isPrim= stack->IsPhysicalPrimary(iMc);
1444
1445     if(!isPrim){//prim or sec
1446       isSecd=1;// is secondary V0s
1447
1448       if(indexMother1 >-1){// && !isPrim){//secondary V0s
1449         //     isSecd=1;// is secondary V0s
1450
1451         //-- check for mother --//
1452         TParticle *mother = stack->Particle(indexMother1);
1453         if(!mother) {
1454           Printf("no mother pointer!");continue;
1455         }
1456         pdgMother = mother->GetPdgCode();
1457         fHistPiPiMonitorMCCuts[1]->Fill(11*fillFlagK0);
1458         fHistPiPMonitorMCCuts[1]->Fill(11*fillFlagL);
1459         fHistPiAPMonitorMCCuts[1]->Fill(11*fillFlagAL);
1460
1461         //-- check for injejcted --//
1462         Bool_t notinjectedMother = kTRUE;
1463         notinjectedMother = fMCev->IsFromBGEvent(indexMother1);
1464
1465         if(fSelectInjected && !notinjectedMother ) continue;//reject injected, optional, set flag in AddTask
1466         fHistPiPiMonitorMCCuts[1]->Fill(10*fillMCtruth);
1467         fHistPiPMonitorMCCuts[1]->Fill(10*fillMCtruth);
1468         fHistPiAPMonitorMCCuts[1]->Fill(10*fillMCtruth);
1469
1470         Bool_t isPrimMother= stack->IsPhysicalPrimary(indexMother1);
1471         if(!isPrimMother) continue;
1472         fHistPiPiMonitorMCCuts[1]->Fill(12*fillFlagK0);
1473         fHistPiPMonitorMCCuts[1]->Fill(12*fillFlagL);
1474         fHistPiAPMonitorMCCuts[1]->Fill(12*fillFlagAL);
1475
1476         uniqueIDmother =  mother->GetUniqueID();
1477
1478         if(uniqueIDmother==13){
1479           continue;
1480         }
1481         fHistPiPiMonitorMCCuts[1]->Fill(13*fillFlagK0);
1482         fHistPiPMonitorMCCuts[1]->Fill(13*fillFlagL);
1483         fHistPiAPMonitorMCCuts[1]->Fill(13*fillFlagAL);
1484
1485
1486         //-- fill secondary V0s histos and pdg histos --//
1487         ptXiMother=mother->Pt();
1488         rapXiMother=mother->Y();
1489
1490
1491         //-- K0s --//
1492         if(pdgCode==310){
1493           if(fabs(pdgMother)==311 || fabs(pdgMother)==313 || fabs(pdgMother)==323 ) isSecd=0; // from K0L,  K0 and K* as primary
1494           else fHistPiPiPDGCode->Fill(fabs(pdgMother));
1495         }
1496
1497         //-- Lambda --//
1498         if(pdgCode==3122){
1499           fHistPiPPDGCode->Fill(fabs(pdgMother));
1500           if (//sigma family
1501               ( TMath::Abs(pdgMother) == 3112) || //sigma minus
1502               ( TMath::Abs(pdgMother) == 3222) || //sigma plus
1503               ( TMath::Abs(pdgMother) == 3224) || //sigma *plus
1504               ( TMath::Abs(pdgMother) == 3114) || //sigma *minus
1505               ( TMath::Abs(pdgMother) == 3214) || //sigma *0 counts as primary????
1506               ( TMath::Abs(pdgMother) == 3212)    //sigma 0 counts as primary
1507               )
1508             {
1509               isSecd=0;
1510             }
1511
1512           if( pdgMother == 3322) //xi0
1513             {
1514               if(!fRapCutV0 || fabs(rapidity)<fRap){
1515                 if(!fEtaCutMCDaughters  ||  (fabs(etaMC00)<fEtaCutMCDaughtersVal|| fabs(etaMC01)<fEtaCutMCDaughtersVal)){
1516                   fHistPiPMassVSPtSecXi[0]->Fill(massV0MC,ptV0MC);
1517                   fHistPiPMassVSYSecXi[0]->Fill(massV0MC,rapidity);
1518                   if(!fRapCutV0 || fabs(rapXiMother)<fRap) fHistPiPXi0PtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1519                 }
1520               }
1521             }
1522           if(pdgMother == 3312) //xi minus
1523             {
1524               if(!fRapCutV0 || fabs(rapidity)<fRap){
1525                 if(!fEtaCutMCDaughters  ||  (fabs(etaMC00)<fEtaCutMCDaughtersVal|| fabs(etaMC01)<fEtaCutMCDaughtersVal)){
1526                   fHistPiPMassVSPtSecXi[0]->Fill(massV0MC,ptV0MC);
1527                   fHistPiPMassVSYSecXi[0]->Fill(massV0MC,rapidity);
1528                   if(!fRapCutV0 || fabs(rapXiMother)<fRap)  fHistPiPXiMinusPtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1529                 }
1530               }
1531             }
1532           if(pdgMother == 3334)//omega-
1533             {
1534               if(!fRapCutV0 || fabs(rapXiMother)<fRap) fHistPiPOmegaPtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1535               fHistPiPMassVSPtSecOmega[0]->Fill(massV0MC,ptV0MC);
1536             }
1537         }
1538
1539         //-- AntiLambda --//
1540         if(pdgCode==-3122 ){
1541           fHistPiAPPDGCode->Fill(fabs(pdgMother));
1542           if (//sigma family
1543               ( TMath::Abs(pdgMother) == 3112) ||//sigma minus
1544               ( TMath::Abs(pdgMother) == 3222) ||//sigma plus
1545               ( TMath::Abs(pdgMother) == 3224) ||//sigma *plus
1546               ( TMath::Abs(pdgMother) == 3114) ||//sigma *minus
1547               ( TMath::Abs(pdgMother) == 3214) || //sigma *0
1548               ( TMath::Abs(pdgMother) == 3212)    //sigma 0 counts as primary
1549               )
1550             {
1551               isSecd=0;
1552             }
1553
1554           if( pdgMother == -3322) //xi0
1555             {
1556               if(!fRapCutV0 || fabs(rapidity)<fRap){
1557                 if(!fEtaCutMCDaughters  ||  (fabs(etaMC00)<fEtaCutMCDaughtersVal|| fabs(etaMC01)<fEtaCutMCDaughtersVal)){
1558                   fHistPiAPMassVSPtSecXi[0]->Fill(massV0MC,ptV0MC);
1559                   fHistPiAPMassVSYSecXi[0]->Fill(massV0MC,rapidity);
1560                   if(!fRapCutV0 || fabs(rapXiMother)<fRap) fHistPiAPXi0PtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1561                 }
1562               }
1563             }
1564           if(pdgMother == -3312) //xi plus
1565             {
1566               if(!fRapCutV0 || fabs(rapidity)<fRap){
1567                 if(!fEtaCutMCDaughters  ||  (fabs(etaMC00)<fEtaCutMCDaughtersVal|| fabs(etaMC01)<fEtaCutMCDaughtersVal)){
1568                   fHistPiAPMassVSPtSecXi[0]->Fill(massV0MC,ptV0MC);
1569                   fHistPiAPMassVSYSecXi[0]->Fill(massV0MC,rapidity);
1570                   if(!fRapCutV0 || fabs(rapXiMother)<fRap) fHistPiAPXiMinusPtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1571                 }
1572               }
1573             }
1574           if(pdgMother == -3334)//omega+
1575             {
1576               fHistPiAPOmegaPtVSLambdaPt[0]->Fill(ptV0MC,ptXiMother);
1577               if(!fRapCutV0 || fabs(rapXiMother)<fRap) fHistPiAPMassVSPtSecOmega[0]->Fill(massV0MC,ptV0MC);
1578             }
1579         }
1580       }
1581     }//end secondaries
1582     else{
1583       //-- check for injejcted --//
1584       Bool_t notinjected = kTRUE;
1585       notinjected = fMCev->IsFromBGEvent(iMc);
1586
1587       if(fSelectInjected && !notinjected ) continue;//removes injected, optional, set flag in AddTask
1588       fHistPiPiMonitorMCCuts[0]->Fill(10*fillMCtruth);
1589       fHistPiPMonitorMCCuts[0]->Fill(10*fillMCtruth);
1590       fHistPiAPMonitorMCCuts[0]->Fill(10*fillMCtruth);
1591     }
1592
1593
1594     //-------------- MC truth mode  --------------//
1595     //---------- fill MC truth histos --------------//
1596     if(fMCTruthMode && !fMCMode){//MC true ana
1597       fHistPiPiMonitorMCCuts[isSecd]->Fill(14*fillFlagK0);
1598       fHistPiPMonitorMCCuts[isSecd]->Fill(14*fillFlagL);
1599       fHistPiAPMonitorMCCuts[isSecd]->Fill(14*fillFlagAL);
1600
1601       //-- DCA daughters --//
1602       // values of one daugher, should be the same
1603       /*
1604       //to primary vertex
1605       trackPos->GetImpactParameters(tdcaPosToVertex[0],tdcaPosToVertex[1]);
1606       trackNeg->GetImpactParameters(tdcaNegToVertex[0],tdcaNegToVertex[1]);
1607
1608       Double_t dcaPosToVertex = TMath::Sqrt(tdcaPosToVertex[0]*tdcaPosToVertex[0]+tdcaPosToVertex[1]*tdcaPosToVertex[1]);
1609       Double_t dcaNegToVertex = TMath::Sqrt(tdcaNegToVertex[0]*tdcaNegToVertex[0]+tdcaNegToVertex[1]*tdcaNegToVertex[1]);
1610       fHistDCADaughtersToPrimVtx[isSecd]->Fill(dcaPosToVertex,dcaNegToVertex);
1611       */
1612
1613
1614       //-- armenteros values --//
1615       TVector3 vecPip;
1616       TVector3 vecPin;
1617
1618       Double_t ptPlus=0, ptMinus=0;
1619       Double_t pt00 = p00->Pt();
1620       Double_t pt01 = p01->Pt();
1621       Double_t  phiPosMC=0.0;
1622
1623       if(p00->GetPdgCode()<0)
1624         {
1625           vecPip.SetXYZ(p01->Px(),p01->Py(),p01->Pz());
1626           vecPin.SetXYZ(p00->Px(),p00->Py(),p00->Pz());
1627           ptMinus = pt00;
1628           ptPlus = pt01;
1629           phiPosMC = p01->Phi();
1630         }
1631       else{
1632         vecPin.SetXYZ(p01->Px(),p01->Py(),p01->Pz());
1633         vecPip.SetXYZ(p00->Px(),p00->Py(),p00->Pz());
1634         ptMinus = pt01;
1635         ptPlus = pt00;
1636         phiPosMC = p00->Phi();
1637       }
1638
1639       TVector3 momTot(p0->Px(),p0->Py(),p0->Pz());
1640       Double_t lQlNeg = fabs(vecPin.Dot(momTot)/momTot.Mag());
1641       Double_t lQlPos = fabs(vecPip.Dot(momTot)/momTot.Mag());
1642       Double_t alfa =0.0;
1643       Double_t den = lQlPos + lQlNeg;
1644       if(den>0) alfa = (lQlPos - lQlNeg)/den;
1645       TVector3 qtvec= vecPin.Cross(momTot);//vecPip.Mag()*sqrt(1-pow(thetapip,2));
1646       Float_t qt = qtvec.Mag()/momTot.Mag();
1647
1648
1649
1650       if(pdgCode == 310) {
1651         fHistPiPiEtaDMC[isSecd]->Fill(etaMC00,ptV0MC);
1652         fHistPiPiEtaDMC[isSecd]->Fill(etaMC01,ptV0MC);
1653       }
1654       if(fabs(pdgCode) == 3122) {
1655         fHistPiPEtaDMC[isSecd]->Fill(etaMC00,ptV0MC);
1656         fHistPiPEtaDMC[isSecd]->Fill(etaMC01,ptV0MC);
1657       }
1658
1659       //-- rapidity and eta cut --//
1660       if(fRapCutV0 && fabs(rapidity)>fRap) continue;//optional, set flag in AddTask
1661       fHistPiPiMonitorMCCuts[isSecd]->Fill(15*fillFlagK0);
1662       fHistPiPMonitorMCCuts[isSecd]->Fill(15*fillFlagL);
1663       fHistPiAPMonitorMCCuts[isSecd]->Fill(15*fillFlagAL);
1664
1665       if(fEtaCutMCDaughters) { if(fabs(etaMC00)>fEtaCutMCDaughtersVal || fabs(etaMC01)>fEtaCutMCDaughtersVal ) continue; }//optional, set flag in AddTask
1666       fHistPiPiMonitorMCCuts[isSecd]->Fill(16*fillFlagK0);
1667       fHistPiPMonitorMCCuts[isSecd]->Fill(16*fillFlagL);
1668       fHistPiAPMonitorMCCuts[isSecd]->Fill(16*fillFlagAL);
1669
1670
1671
1672       //-- Fill Particle histos --//
1673       if(pdgCode==310){//K0s
1674         fHistPiPiMonitorMCCuts[isSecd]->Fill(17);
1675
1676         fHistPiPiEtaDMC[1]->Fill(etaMC00,ptV0MC);
1677         fHistPiPiEtaDMC[1]->Fill(etaMC01,ptV0MC);
1678
1679         fHistPiPiMass[isSecd]->Fill(massV0MC);
1680         fHistPiPiMassVSPt[isSecd]->Fill(massV0MC,ptV0MC);
1681         // fHistPiPiPtDaughters[isSecd]->Fill(ptMinus,ptPlus);
1682         fHistPiPiPtVSY[isSecd]->Fill(rapidity,ptV0MC);
1683
1684         Double_t ctTK0s=0.0,ctK0s=0.0;
1685         if(pV0MC>0.0) ctK0s=declength3d*0.497614/pV0MC;
1686         if(ptV0MC>0.0) ctTK0s=declength*0.497614/ptV0MC;
1687         fHistPiPiDecayLengthResolution[0]->Fill(declength3d,declength);
1688         fHistPiPiDecayLengthVsPt[isSecd]->Fill(ptV0MC,declength);
1689         fHistPiPiDecayLengthVsCtau[isSecd]->Fill(massV0MC,ctTK0s);
1690         fHistPiPiDecayLengthVsMass[isSecd]->Fill(massV0MC,declength);
1691         //all V0s histo
1692         fHistArmenteros[isSecd]->Fill(alfa,qt);
1693         fHistV0RadiusZ[isSecd]->Fill(rMC2D,xyzMC[2]);
1694         fHistV0RadiusXY[isSecd]->Fill(xyzMC[0],xyzMC[1]);
1695         fHistV0RadiusZVSPt[isSecd]->Fill(ptV0MC,xyzMC[2]);
1696       }
1697       if(pdgCode==3122){ //Lambda
1698         fHistPiPMonitorMCCuts[isSecd]->Fill(17);
1699
1700         fHistPiPEtaDMC[1]->Fill(etaMC00,ptV0MC);
1701         fHistPiPEtaDMC[1]->Fill(etaMC01,ptV0MC);
1702
1703         fHistPiPMassVSPt[isSecd]->Fill(massV0MC,ptV0MC);
1704         fHistPiPMass[isSecd]->Fill(massV0MC);
1705         //  fHistPiPPtDaughters[isSecd]->Fill(ptMinus,ptPlus);
1706         fHistPiPPtVSY[isSecd]->Fill(rapidity,ptV0MC);
1707
1708         Double_t ctTL=0.0, ctL=0.0;
1709         if(pV0MC>0.0) ctL=declength3d*1.115683/pV0MC;
1710         if(ptV0MC>0.0) ctTL=declength*1.115683/ptV0MC;
1711         fHistPiPDecayLengthResolution[0]->Fill(declength3d,declength);
1712         fHistPiPDecayLengthVsPt[isSecd]->Fill(ptV0MC,declength);
1713         fHistPiPDecayLengthVsCtau[isSecd]->Fill(massV0MC,ctTL);
1714         fHistPiPDecayLengthVsMass[isSecd]->Fill(massV0MC,declength);
1715         //all V0s hito
1716         fHistArmenteros[isSecd]->Fill(alfa,qt);
1717         fHistV0RadiusZ[isSecd]->Fill(rMC2D,xyzMC[2]);
1718         fHistV0RadiusXY[isSecd]->Fill(xyzMC[0],xyzMC[1]);
1719         fHistV0RadiusZVSPt[isSecd]->Fill(ptV0MC,xyzMC[2]);
1720       }
1721
1722       if(pdgCode==-3122){ //AntiLambda
1723         fHistPiAPMonitorMCCuts[isSecd]->Fill(17);
1724
1725         fHistPiPEtaDMC[1]->Fill(etaMC00,ptV0MC);
1726         fHistPiPEtaDMC[1]->Fill(etaMC01,ptV0MC);
1727
1728         fHistPiAPMassVSPt[isSecd]->Fill(massV0MC,ptV0MC);
1729         fHistPiAPMass[isSecd]->Fill(massV0MC);
1730         //  fHistPiAPPtDaughters[isSecd]->Fill(ptMinus,ptPlus);
1731         fHistPiAPPtVSY[isSecd]->Fill(rapidity,ptV0MC);
1732
1733         Double_t ctTL=0.0, ctL=0.0;
1734         if(pV0MC>0.0) ctL=declength3d*1.115683/pV0MC;
1735         if(ptV0MC>0.0) ctTL=declength*1.115683/ptV0MC;
1736         fHistPiAPDecayLengthResolution[0]->Fill(declength3d,declength);
1737         fHistPiAPDecayLengthVsPt[isSecd]->Fill(ptV0MC,declength);
1738         fHistPiAPDecayLengthVsCtau[isSecd]->Fill(massV0MC,ctL);
1739         fHistPiAPDecayLengthVsMass[isSecd]->Fill(massV0MC,ctTL);
1740         //all V0s histo
1741         fHistArmenteros[isSecd]->Fill(alfa,qt);
1742         fHistV0RadiusZ[isSecd]->Fill(rMC2D,xyzMC[2]);
1743         fHistV0RadiusXY[isSecd]->Fill(xyzMC[0],xyzMC[1]);
1744         fHistV0RadiusZVSPt[isSecd]->Fill(ptV0MC,xyzMC[2]);
1745         fHistV0RadiusXYVSY[isSecd]->Fill(rapidity,rMC2D);
1746       }
1747     }//MC true ana
1748     else{//MC reco ana
1749       V0RecoLoop(id0,id1,isSecd,pdgCode,ptV0MC,pdgMother,ptXiMother,declength);
1750     }
1751
1752   }//end MC stack loop
1753 }
1754 //________________________________________________________________________
1755 void AliAnalysisTaskV0ForRAA::V0RecoLoop(Int_t id0,Int_t id1,Int_t isSecd,Int_t what,Double_t ptV0MC, Int_t pdgMother,Double_t ptXiMother,Double_t declengthV0MC){
1756   //loop over reconstructed particles MC or Data
1757
1758
1759   //--------------------- define variables -----------------------//
1760   Double_t pp[3];
1761   Double_t pm[3];
1762   Double_t xr[3];
1763
1764   Double_t massPi=0.13957018;
1765   Double_t massP=0.93827203;
1766
1767   TLorentzVector positivesMIP;
1768   TLorentzVector negativesMIAP;
1769   TLorentzVector positivesMIPi;
1770   TLorentzVector negativesMIPi;
1771
1772   /*
1773     AliKFParticle::SetField(fESD->GetMagneticField());
1774     AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
1775     AliKFVertex primVtxImproved = primVtx;
1776
1777     AliKFParticle* negPiKF=NULL;
1778     AliKFParticle* posPiKF=NULL;
1779     AliKFParticle* posPKF=NULL;
1780     AliKFParticle* negAPKF=NULL;
1781   */
1782
1783   AliESDtrack* trackPos=NULL;
1784   AliESDtrack* trackNeg=NULL;
1785   AliESDtrack* trackPosTest = NULL;
1786   AliESDtrack* trackNegTest =NULL;
1787
1788   Double_t primaryVtxPosition[3];
1789   primaryVtxPosition[0] = fESD->GetPrimaryVertex()->GetXv();
1790   primaryVtxPosition[1] = fESD->GetPrimaryVertex()->GetYv();
1791   primaryVtxPosition[2] = fESD->GetPrimaryVertex()->GetZv();
1792
1793   Int_t nV0 = fESD->GetNumberOfV0s();
1794   AliESDv0 * v0MIs=NULL;
1795
1796
1797   Bool_t stopLoop = kFALSE;
1798   //------------------------ V0 reco loop --------------------//
1799   for(Int_t iV0MI = 0; iV0MI < nV0; iV0MI++) {//V0 loop
1800
1801     if(stopLoop) break; //in case of MC if MC truth partner is found
1802
1803     //-- get V0 info --//
1804     v0MIs = fESD->GetV0(iV0MI);
1805     if(!v0MIs ) continue;
1806
1807     fHistPiPiMonitorCuts[isSecd]->Fill(1);
1808     fHistPiPMonitorCuts[isSecd]->Fill(1);
1809     fHistPiAPMonitorCuts[isSecd]->Fill(1);
1810
1811     //------------ get references of daughters --------------//
1812     //-- esd tracks --//
1813     trackPosTest = fESD->GetTrack(v0MIs->GetPindex());
1814     trackNegTest = fESD->GetTrack(v0MIs->GetNindex());
1815
1816     if( trackPosTest->GetSign() == trackNegTest->GetSign()) continue;//reject same signed daughters
1817
1818     fHistPiPiMonitorCuts[isSecd]->Fill(2);
1819     fHistPiPMonitorCuts[isSecd]->Fill(2);
1820     fHistPiAPMonitorCuts[isSecd]->Fill(2);
1821
1822     Bool_t onthefly = v0MIs->GetOnFlyStatus();
1823     if(fOntheFly!=onthefly) continue;// select on-the-fly or offline V0 in AddTask
1824
1825     fHistPiPiMonitorCuts[isSecd]->Fill(3);
1826     fHistPiPMonitorCuts[isSecd]->Fill(3);
1827     fHistPiAPMonitorCuts[isSecd]->Fill(3);
1828
1829
1830     //-- for MC mode --//
1831     if(fMCMode){
1832       //check MC labels (and find partners for MC truth V0 daughters for fMCTruthMode=kTRUE)
1833       if(!GetMCTruthPartner(trackPosTest,trackNegTest,id0,id1)) continue;
1834       else stopLoop = kTRUE;
1835     }
1836
1837     fHistPiPiMonitorCuts[isSecd]->Fill(4);
1838     fHistPiPMonitorCuts[isSecd]->Fill(4);
1839     fHistPiAPMonitorCuts[isSecd]->Fill(4);
1840
1841
1842     //--  get eta from V0 daughters --//
1843     Double_t posDaughterEta=0.0;
1844     Double_t negDaughterEta=0.0;
1845     Double_t posDaughterPhi=0.0;
1846     Double_t negDaughterPhi=0.0;
1847
1848     Double_t eta00 = trackPosTest->Eta();
1849     Double_t eta01 = trackNegTest->Eta();
1850
1851     Int_t indexPos = 0,indexNeg=0;
1852
1853
1854     //---------- check sign assignment for daughters --------//
1855     Bool_t switchSign = kFALSE;
1856
1857     if( trackPosTest->GetSign() >0){//pos
1858       indexPos = v0MIs->GetPindex();
1859       indexNeg = v0MIs->GetNindex();
1860
1861       v0MIs->GetPPxPyPz(pp[0],pp[1],pp[2]);
1862       v0MIs->GetNPxPyPz(pm[0],pm[1],pm[2]);
1863
1864       posDaughterEta = v0MIs->GetParamP()->Eta();
1865       negDaughterEta = v0MIs->GetParamN()->Eta();
1866       posDaughterPhi = v0MIs->GetParamP()->Phi();
1867       negDaughterPhi = v0MIs->GetParamN()->Phi();
1868       /*
1869               if (negPiKF) delete negPiKF; negPiKF=NULL;
1870               if (posPiKF) delete posPiKF; posPiKF=NULL;
1871               if (posPKF) delete posPKF; posPKF=NULL;
1872               if (negAPKF) delete negAPKF; negAPKF=NULL;
1873
1874               negPiKF = new AliKFParticle( *(v0MIs->GetParamN()) ,-211);
1875               posPiKF = new AliKFParticle( *(v0MIs->GetParamP()) ,211);
1876               posPKF = new AliKFParticle( *(v0MIs->GetParamP()) ,2212);
1877               negAPKF = new AliKFParticle( *(v0MIs->GetParamN()) ,-2212);
1878       */
1879
1880     }
1881     if( trackPosTest->GetSign() <0){//neg
1882
1883       indexPos = v0MIs->GetNindex();
1884       indexNeg = v0MIs->GetPindex();
1885
1886       v0MIs->GetNPxPyPz(pp[0],pp[1],pp[2]);
1887       v0MIs->GetPPxPyPz(pm[0],pm[1],pm[2]);
1888
1889       posDaughterEta = v0MIs->GetParamN()->Eta();
1890       negDaughterEta = v0MIs->GetParamP()->Eta();
1891       posDaughterPhi = v0MIs->GetParamN()->Phi();
1892       negDaughterPhi = v0MIs->GetParamP()->Phi();
1893       /*
1894         if (negPiKF) delete negPiKF; negPiKF=NULL;
1895         if (posPiKF) delete posPiKF; posPiKF=NULL;
1896         if (posPKF) delete posPKF; posPKF=NULL;
1897         if (negAPKF) delete negAPKF; negAPKF=NULL;
1898
1899
1900         negPiKF = new AliKFParticle( *(v0MIs->GetParamP()) ,-211);
1901         posPiKF = new AliKFParticle( *(v0MIs->GetParamN()) ,211);
1902         posPKF = new AliKFParticle( *(v0MIs->GetParamN()) ,2212);
1903         negAPKF = new AliKFParticle( *(v0MIs->GetParamP()) ,-2212);
1904       */
1905
1906       switchSign = kTRUE;
1907       eta01 = trackPosTest->Eta();
1908       eta00 = trackNegTest->Eta();
1909
1910     }
1911
1912     trackPos =fESD->GetTrack(indexPos);
1913     trackNeg =fESD->GetTrack(indexNeg);
1914
1915     // ------------- calc masses and 4 vectors -------------- //
1916     //K0
1917     positivesMIPi.SetXYZM(pp[0],pp[1],pp[2],massPi);
1918     negativesMIPi.SetXYZM(pm[0],pm[1],pm[2],massPi);
1919     TLorentzVector v0K0=positivesMIPi+negativesMIPi;
1920
1921     //Lambda
1922     positivesMIP.SetXYZM(pp[0],pp[1],pp[2],massP);
1923     TLorentzVector v0Lambda=positivesMIP+negativesMIPi;
1924
1925     //Anitlambda
1926     negativesMIAP.SetXYZM(pm[0],pm[1],pm[2],massP);
1927     TLorentzVector v0ALambda=positivesMIPi+negativesMIAP;
1928
1929
1930     TVector3 ppTrack(pp);
1931     TVector3 pmTrack(pm);
1932
1933     //    Double_t posDaughterPt = ppTrack.Pt();
1934     //    Double_t negDaughterPt = pmTrack.Pt();
1935
1936     Double_t posDaughterP = ppTrack.Mag();
1937     Double_t negDaughterP = pmTrack.Mag();
1938
1939     //---------------------AliKFParticle ---------------------//
1940     /*
1941         Double_t chi2K0C=0.0;
1942         Double_t chi2LambdaC=0.0;
1943         Double_t chi2ALambdaC=0.0;
1944
1945
1946         AliKFParticle v0K0KF;
1947         v0K0KF +=(*negPiKF);
1948         v0K0KF +=(*posPiKF);
1949         //v0K0C.SetVtxGuess(xr[0],xr[1],xr[2]);
1950         v0K0KF.SetProductionVertex(primVtxImproved);
1951
1952         AliKFParticle v0LambdaKF;
1953         v0LambdaKF +=(*negPiKF);
1954         v0LambdaKF +=(*posPKF);
1955         //v0LambdaC.SetVtxGuess(xr[0],xr[1],xr[2]);
1956         v0LambdaKF.SetProductionVertex(primVtxImproved);
1957
1958         AliKFParticle v0ALambdaKF;
1959         v0ALambdaKF +=(*negAPKF);
1960         v0ALambdaKF +=(*posPiKF);
1961         //v0ALambdaC.SetVtxGuess(xr[0],xr[1],xr[2]);
1962         v0ALambdaKF.SetProductionVertex(primVtxImproved);
1963
1964         if( v0K0KF.GetNDF() != 0) {
1965         chi2K0C = v0K0KF.GetChi2()/v0K0KF.GetNDF();
1966         }
1967
1968         Double_t chi2LambdaC=100000.;
1969         if( v0LambdaKF.GetNDF() != 0) {
1970         chi2LambdaC = v0LambdaKF.GetChi2()/v0LambdaKF.GetNDF();
1971         }
1972
1973         Double_t chi2ALambdaC=100000.;
1974         if( v0ALambdaKF.GetNDF() != 0) {
1975         chi2ALambdaC = v0ALambdaKF.GetChi2()/v0ALambdaKF.GetNDF();
1976         }
1977     */
1978
1979
1980     //------------------- DCA and decay radius V0 -------------//
1981     v0MIs->GetXYZ(xr[0],xr[1],xr[2]);
1982
1983     Double_t dim2V0Radius= sqrt( pow(xr[0] - primaryVtxPosition[0],2.0)
1984                                  +pow(xr[1] - primaryVtxPosition[1],2.0));
1985
1986     Double_t decayLength = sqrt( pow(xr[0] - primaryVtxPosition[0],2.0)
1987                                  +pow(xr[1] - primaryVtxPosition[1],2.0)
1988                                  +pow(xr[2] - primaryVtxPosition[2],2.0));
1989
1990     Double_t  dcaV0ToPrimVertex= v0MIs->GetD(primaryVtxPosition[0],primaryVtxPosition[1]);
1991     ////v0K0KF.GetDistanceFromVertexXY(tPrimaryVtxPosition);
1992
1993
1994     //-------------------- general V0 cuts -------------------//
1995
1996     //-- eta cut --//
1997     if( fabs(posDaughterEta) > fEtaCutMCDaughtersVal || fabs(negDaughterEta) > fEtaCutMCDaughtersVal) continue;
1998     fHistPiPiMonitorCuts[isSecd]->Fill(6);
1999     fHistPiPMonitorCuts[isSecd]->Fill(6);
2000     fHistPiAPMonitorCuts[isSecd]->Fill(6);
2001
2002     //-- pt cut --//
2003     //if( fabs(posDaughterPt)<fMinPt || fabs(negDaughterPt) < fMinPt ) continue;
2004     // fHistPiPiMonitorCuts[isSecd]->Fill(7);
2005     // fHistPiPMonitorCuts[isSecd]->Fill(7);
2006     // fHistPiAPMonitorCuts[isSecd]->Fill(7);
2007
2008
2009     //-- radius xy min cut --//
2010     if(dim2V0Radius < fDecayRadXYMin) continue;
2011     fHistPiPiMonitorCuts[isSecd]->Fill(7);
2012     fHistPiPMonitorCuts[isSecd]->Fill(7);
2013     fHistPiAPMonitorCuts[isSecd]->Fill(7);
2014
2015     //-- radius xy max cut --//
2016     if(dim2V0Radius > fDecayRadXYMax) continue;
2017     fHistPiPiMonitorCuts[isSecd]->Fill(8);
2018     fHistPiPMonitorCuts[isSecd]->Fill(8);
2019     fHistPiAPMonitorCuts[isSecd]->Fill(8);
2020
2021     //-- decay length min ->ctau --//
2022     if(decayLength > fDecayLengthMax) continue;
2023     fHistPiPiMonitorCuts[isSecd]->Fill(9);
2024     fHistPiPMonitorCuts[isSecd]->Fill(9);
2025     fHistPiAPMonitorCuts[isSecd]->Fill(9);
2026
2027     //-- decay length min cut --//
2028     if(decayLength < fDecayLengthMin) continue;
2029     fHistPiPiMonitorCuts[isSecd]->Fill(10);
2030     fHistPiPMonitorCuts[isSecd]->Fill(10);
2031     fHistPiAPMonitorCuts[isSecd]->Fill(10);
2032
2033
2034     //----------------------- V0 variables --------------------//
2035     //-- armenteros --//
2036     TVector3 momTot = ppTrack + pmTrack;
2037     Double_t lQlNeg = fabs(pmTrack.Dot(momTot)/momTot.Mag());
2038     Double_t lQlPos = fabs(ppTrack.Dot(momTot)/momTot.Mag());
2039
2040     Double_t alfa =0.0;
2041     Double_t den = lQlPos + lQlNeg;
2042     if(den>0) alfa = (lQlPos - lQlNeg)/den;
2043     TVector3 qtvec= pmTrack.Cross(momTot);
2044     Double_t qt = qtvec.Mag()/momTot.Mag();
2045
2046     //-- momenta --//
2047     //pT
2048     Double_t ptK0s = v0K0.Pt();
2049     Double_t ptLambda = v0Lambda.Pt();
2050     Double_t ptALambda = v0ALambda.Pt();
2051     //p
2052     Double_t pK0s = v0K0.P();
2053     Double_t pLambda = v0Lambda.P();
2054     Double_t pALambda = v0ALambda.P();
2055
2056     //-- masses --//
2057     Double_t massK0s = v0K0.M();
2058     Double_t massLambda = v0Lambda.M();
2059     Double_t massALambda = v0ALambda.M();
2060
2061     //-- rapidity --//
2062     Double_t rapK0s = v0MIs->Y(310);
2063     Double_t rapL   = v0MIs->Y(3122);
2064     Double_t rapAL  = v0MIs->Y(3122);
2065
2066     //-- other variables --//
2067     Double_t opAng =   fabs(ppTrack.Angle(pmTrack));
2068     /*
2069       Double_t phiK0 =  v0K0.Phi();
2070       Double_t phiL =  v0Lambda.Phi();
2071       Double_t phiAL =  v0ALambda.Phi();
2072     */
2073     Double_t cosOPAng = v0MIs->GetV0CosineOfPointingAngle();
2074
2075     /*
2076     //introduce more histo
2077     Double_t errOnMassK0s = v0MIs->ChangeMassHypothesis(310);
2078     Double_t errOnMassLambda = 0.0;
2079     Double_t errOnMassALambda = 0.0;
2080     if(!switchSign){
2081     errOnMassLambda  = v0MIs->ChangeMassHypothesis(3122);
2082     errOnMassALambda = v0MIs->ChangeMassHypothesis(-3122);
2083     }
2084     else{
2085     errOnMassLambda  = v0MIs->ChangeMassHypothesis(-3122);
2086     errOnMassALambda = v0MIs->ChangeMassHypothesis(3122);
2087     }
2088     */
2089
2090     //-- DCAbetween the daughters --//
2091     Double_t dcaDaughters = v0MIs->GetDcaV0Daughters();
2092
2093     //-- DCA daughters to primary vertex --//
2094     Double_t dcaPosToVertex=0.0,dcaNegToVertex=0.0;
2095     Double_t dzPos=(primaryVtxPosition[0]-xr[0])*ppTrack.Y() - (primaryVtxPosition[1]-xr[1])*ppTrack.X();
2096     dcaPosToVertex=TMath::Sqrt(dzPos*dzPos/(pow(ppTrack.X(),2)+pow(ppTrack.Y(),2)));
2097     Double_t dzNeg=(primaryVtxPosition[0]-xr[0])*pmTrack.Y() - (primaryVtxPosition[1]-xr[1])*pmTrack.X();
2098     dcaNegToVertex=TMath::Sqrt(dzNeg*dzNeg/(pow(pmTrack.X(),2)+pow(pmTrack.Y(),2)));
2099
2100     // Double_t dcaPosToVertex[3];dcaNegToVertex[3];
2101     // trackPos->GetImpactParameters(dcaPosToVertex[0],dcaPosToVertex[1]);
2102     // trackNeg->GetImpactParameters(dcaNegToVertex[0],dcaNegToVertex[1]);
2103
2104     // dcaPosToVertex = TMath::Sqrt(dcaPosToVertex[0]*dcaPosToVertex[0]+dcaPosToVertex[1]*dcaPosToVertex[1]);
2105     // dcaNegToVertex = TMath::Sqrt(dcaNegToVertex[0]*dcaNegToVertex[0]+dcaNegToVertex[1]*dcaNegToVertex[1]);
2106
2107     //dcaPosToVertex  =   posPKF->GetDistanceFromVertexXY(primaryVtxPosition);
2108     //dcaNegToVertex  =   negPiKF->GetDistanceFromVertexXY(primaryVtxPosition);
2109
2110
2111     //------------------------ detector values -------------------------------//
2112     //-- TPC ITS values pos --//
2113     Int_t nclsTPCPos =  trackPos->GetNcls(1);
2114     Int_t nclsTPCFindablePos =  trackPos->GetTPCNclsF();
2115     Int_t nclsITSPos =  trackPos->GetNcls(0);
2116     Double_t chi2PerClusterITSPos = -1.0;
2117     if(nclsITSPos>0) chi2PerClusterITSPos = trackPos->GetITSchi2()/Double_t(nclsITSPos);
2118     Double_t crossedRowsTPCPos = trackPos->GetTPCCrossedRows();
2119
2120     //-- TPC ITS values neg --//
2121     Int_t nclsTPCNeg =  trackNeg->GetNcls(1);
2122     Int_t nclsTPCFindableNeg =  trackNeg->GetTPCNclsF();
2123     Int_t nclsITSNeg =  trackNeg->GetNcls(0);
2124     Double_t chi2PerClusterITSNeg = -1.0;
2125     if(nclsITSNeg>0) chi2PerClusterITSNeg =trackNeg->GetITSchi2()/Double_t(nclsITSNeg);
2126     Double_t crossedRowsTPCNeg = trackNeg->GetTPCCrossedRows();
2127
2128     //-- ratio TPC crossed rows to clusters --//
2129     Double_t ratio = 10.0;
2130     if(nclsTPCFindableNeg >0.0) ratio        = double(crossedRowsTPCNeg)/ double(nclsTPCFindableNeg);
2131
2132     Double_t ratioPos = 10.0;
2133     if(nclsTPCFindablePos >0.0) ratioPos     = double(crossedRowsTPCPos)/ double(nclsTPCFindablePos);
2134
2135     Double_t ratioFoFi = 10.0;
2136     if(nclsTPCFindableNeg >0.0) ratioFoFi    = double(nclsTPCNeg)/ double(nclsTPCFindableNeg);
2137
2138     Double_t ratioFoFiPos = 10.0;
2139     if(nclsTPCFindablePos >0.0) ratioFoFiPos = double(nclsTPCPos)/ double(nclsTPCFindablePos);
2140
2141
2142     //-------------- cuts on det. values ---------------------------------------//
2143     if(fMoreNclsThanRows && (crossedRowsTPCPos < nclsTPCPos || crossedRowsTPCNeg < nclsTPCNeg  )) continue;
2144     fHistPiPiMonitorCuts[isSecd]->Fill(11);
2145     fHistPiPMonitorCuts[isSecd]->Fill(11);
2146     fHistPiAPMonitorCuts[isSecd]->Fill(11);
2147
2148     if(fMoreNclsThanFindable && (nclsTPCFindablePos < nclsTPCPos || nclsTPCFindableNeg < nclsTPCNeg  )) continue;
2149     fHistPiPiMonitorCuts[isSecd]->Fill(12);
2150     fHistPiPMonitorCuts[isSecd]->Fill(12);
2151     fHistPiAPMonitorCuts[isSecd]->Fill(12);
2152
2153     if(fMoreNclsThanFindableMax && ( nclsTPCPos < (nclsTPCFindablePos -60.0) || nclsTPCNeg < (nclsTPCFindableNeg - 60.0) )) continue;
2154     // if(chi2PerClusterITSNeg > fChi2PerClusterITS || chi2PerClusterITSPos > fChi2PerClusterITS ) continue;
2155     fHistPiPiMonitorCuts[isSecd]->Fill(13);
2156     fHistPiPMonitorCuts[isSecd]->Fill(13);
2157     fHistPiAPMonitorCuts[isSecd]->Fill(13);
2158
2159     if(ratio > fRatioMaxCRowsOverFindable || ratioPos > fRatioMaxCRowsOverFindable) continue;
2160
2161     if(ratioFoFi < fRatioFoundOverFindable || ratioFoFiPos < fRatioFoundOverFindable) continue;
2162     fHistPiPiMonitorCuts[isSecd]->Fill(14);
2163     fHistPiPMonitorCuts[isSecd]->Fill(14);
2164     fHistPiAPMonitorCuts[isSecd]->Fill(14);
2165
2166
2167
2168     //------------------- cuts on V0 -------------------------------------------------------------------//
2169     //-- cut flags for V0 specific cuts --//
2170     Bool_t cutOKK0s = kTRUE;
2171     Bool_t cutOKLambda = kTRUE;
2172     Bool_t cutOKALambda = kTRUE;
2173
2174     // ----------------- for MC mode ------------------------ //
2175     Bool_t fillK0sMC = kTRUE;
2176     Bool_t fillLambdaMC = kTRUE;
2177     Bool_t fillALambdaMC = kTRUE;
2178
2179     if(fMCMode && fMCTruthMode) {
2180       if(what == 310) {
2181         fillLambdaMC = kFALSE;
2182         fillALambdaMC = kFALSE;
2183       }
2184       else if(what == 3122){
2185         fillALambdaMC = kFALSE;
2186         fillK0sMC = kFALSE;
2187       }
2188       else if(what == -3122){
2189         fillLambdaMC = kFALSE;
2190         fillK0sMC = kFALSE;
2191       }
2192     }
2193
2194     //-------------------------- K0 cuts -----------------------------//
2195
2196     if(dcaV0ToPrimVertex > fDCAToVertexK0)  cutOKK0s = kFALSE;
2197     if(fabs(xr[2])>fDCAZ) continue;//like decay radius z component
2198     fHistPiPiMonitorCuts[isSecd]->Fill(15);
2199
2200     Double_t ctK0 = 0.0,ctTK0 = 0.0;
2201     if(fabs(pK0s)>0.0)  ctK0 = decayLength*0.497614/pK0s;//ctau in 3D
2202     if(fabs(ptK0s)>0.0)  ctTK0 = dim2V0Radius*0.497614/ptK0s;// "transverse" ctau
2203     if(ctK0 > fCtauK0s &&  fabs(ptK0s) <fCtauPtCutK0) cutOKK0s = kFALSE;
2204     else  fHistPiPiMonitorCuts[isSecd]->Fill(16);
2205
2206     if((cosOPAng < fCosPointAngK && fabs(ptK0s) < fCPAPtCutK0)|| cosOPAng<0.99)
2207       cutOKK0s = kFALSE;//cosine of pointing angle
2208     else  fHistPiPiMonitorCuts[isSecd]->Fill(17);
2209
2210     if(dcaDaughters > fDCADaughtersK0 )cutOKK0s = kFALSE;// dca daughters
2211     else  fHistPiPiMonitorCuts[isSecd]->Fill(18);
2212
2213     if(dcaNegToVertex < fDCADaughtersToVtxSmall || dcaPosToVertex < fDCADaughtersToVtxSmall)
2214       cutOKK0s = kFALSE; //dca neg and pos daughter to vertex
2215     else  fHistPiPiMonitorCuts[isSecd]->Fill(19);
2216
2217     if(fRapCutV0 && fabs(rapK0s) > fRap) cutOKK0s = kFALSE;//rapidity
2218     else  fHistPiPiMonitorCuts[isSecd]->Fill(20);
2219     /*
2220         if(chi2K0C > fChiCutKf) cutOKK0s = kFALSE;
2221         if(opAng < fOpengAngleDaughters && fabs(ptK0s) < fOpAngPtCut )  cutOKK0s = kFALSE;//opening angle cut
2222         else
2223     */
2224
2225     if(fabs(1.115 - massLambda) < fExcludeLambdaFromK0s) cutOKK0s = kFALSE;//exclude Lambdas from K0s
2226   if(fabs(1.115 - massALambda) < fExcludeLambdaFromK0s) cutOKK0s = kFALSE;//exclude AntiLambdas from K0s
2227     fHistPiPiMonitorCuts[isSecd]->Fill(21);
2228
2229     Bool_t ptbinokK0s=kFALSE;
2230     if( ptK0s < fQtCut &&  ptK0s > fQtCutPtLow ) ptbinokK0s=kTRUE;// cut below pt threshhold
2231     Double_t qtval  = fArmQtSlope*fabs(alfa);
2232     if(fArmCutK0 && ptbinokK0s && qt < qtval) cutOKK0s = kFALSE; //cut in armenteros-podolanski diagram
2233     else  fHistPiPiMonitorCuts[isSecd]->Fill(22);
2234
2235     if( ptK0s > fPtTPCCut){
2236       if(fESDTrackCuts){
2237         if(!fESDTrackCuts->AcceptTrack(trackPosTest)) cutOKK0s = kFALSE; //esd track cuts for pion
2238         if(!fESDTrackCuts->AcceptTrack(trackNegTest)) cutOKK0s = kFALSE; //esd track cuts for pion
2239       }
2240       else  fHistPiPiMonitorCuts[isSecd]->Fill(5);
2241     }
2242     else{
2243       if(fESDTrackCutsLowPt){
2244         if(!fESDTrackCutsLowPt->AcceptTrack(trackPosTest))  cutOKK0s = kFALSE; //esd track cuts at low pt
2245         if(!fESDTrackCutsLowPt->AcceptTrack(trackNegTest))  cutOKK0s = kFALSE; //esd track cuts at low pt
2246       }
2247     }
2248
2249
2250
2251     //-------------------------- Lambda cuts -------------------------//
2252     if(dcaV0ToPrimVertex > fDCAToVertexL) cutOKLambda = kFALSE;
2253     if(fabs(xr[2])>fDCAZ) continue;//like decay radius z component
2254     fHistPiPMonitorCuts[isSecd]->Fill(15);
2255
2256     Double_t ctL = 0.0,ctTL=0.0;
2257     if(fabs(pLambda)>0.0)  ctL  = decayLength*1.115683/fabs(pLambda);//ctau in 3D
2258     if(fabs(ptLambda)>0.0) ctTL = dim2V0Radius*1.115683/fabs(ptLambda);// "transverse" ctau
2259
2260     if(ctL > fCtauL && fabs(ptLambda) <fCtauPtCutL)  cutOKLambda = kFALSE;
2261     else  fHistPiPMonitorCuts[isSecd]->Fill(16);
2262
2263     if((cosOPAng<fCosPointAngL && fabs(ptLambda) < fCPAPtCutL)|| cosOPAng<0.99)
2264       cutOKLambda = kFALSE;//cosine of pointing angle
2265     else fHistPiPMonitorCuts[isSecd]->Fill(17);
2266
2267     if(dcaDaughters > fDCADaughtersL )cutOKLambda = kFALSE;// dca daughters
2268     else  fHistPiPMonitorCuts[isSecd]->Fill(18);
2269
2270     if( dcaNegToVertex < fDCADaughtersToVtxSmall || dcaPosToVertex < fDCADaughtersToVtxLarge)
2271       cutOKLambda = kFALSE;//dca negand pos  daughter to vertex
2272     else  fHistPiPMonitorCuts[isSecd]->Fill(19);
2273
2274     if(fRapCutV0 && fabs(rapL) > fRap) cutOKLambda = kFALSE;//rapidity
2275     else  fHistPiPMonitorCuts[isSecd]->Fill(20);
2276
2277     if(fabs(0.497 - massK0s) < fExcludeK0sFromLambda) cutOKLambda = kFALSE;//exclude K0s from Lambda
2278
2279     /*
2280          if(chi2LambdaC > fChiCutKf) cutOKLambda = kFALSE;
2281          else  fHistPiPMonitorCuts[isSecd]->Fill(20);
2282     */
2283
2284     if(opAng < fOpengAngleDaughters && fabs(ptLambda) < fOpAngPtCut )  cutOKLambda = kFALSE;//opening angle cut
2285     else  fHistPiPMonitorCuts[isSecd]->Fill(21);
2286
2287
2288     if(alfa<fAlfaCut  || (fArmCutL && qt>qtval)) cutOKLambda = kFALSE; //cut in armenteros-podolanski diagram
2289     else  fHistPiPMonitorCuts[isSecd]->Fill(22);
2290
2291     if(ptLambda > fPtTPCCut){
2292       if(fESDTrackCuts && fESDTrackCutsCharged){
2293         if(!fESDTrackCutsCharged->AcceptTrack(trackPosTest)) cutOKLambda = kFALSE;//esd track cuts for proton
2294         if(!fESDTrackCuts->AcceptTrack(trackNegTest)) cutOKLambda = kFALSE;//esd track cuts for pion
2295       }
2296       else  fHistPiPMonitorCuts[isSecd]->Fill(5);
2297     }
2298     else{
2299       if(fESDTrackCutsLowPt){
2300         if(!fESDTrackCutsLowPt->AcceptTrack(trackPosTest))  cutOKLambda = kFALSE; //esd track cuts at low pt
2301         if(!fESDTrackCutsLowPt->AcceptTrack(trackNegTest))  cutOKLambda = kFALSE; //esd track cuts at low pt
2302       }
2303     }
2304
2305
2306     //--------------------------- ALambda cuts --------------------------//
2307
2308     if(dcaV0ToPrimVertex > fDCAToVertexL) cutOKALambda = kFALSE;
2309     if(fabs(xr[2])>fDCAZ) continue;//like decay radius z component
2310     fHistPiAPMonitorCuts[isSecd]->Fill(15);
2311
2312     Double_t ctAL = 0.0,ctTAL=0.0;
2313     if(fabs(pALambda)>0.0)  ctAL  = decayLength*1.115683/fabs(pALambda);//ctau in 3D
2314     if(fabs(ptALambda)>0.0) ctTAL = dim2V0Radius*1.115683/fabs(ptALambda);// "transverse" ctau
2315     if(ctAL > fCtauL &&  fabs(ptALambda) <fCtauPtCutL)  cutOKALambda = kFALSE;
2316     else  fHistPiAPMonitorCuts[isSecd]->Fill(16);
2317
2318     if((cosOPAng<fCosPointAngL && fabs(ptALambda) < fCPAPtCutL)|| cosOPAng<0.99)
2319       cutOKALambda = kFALSE;//cosine of pointing angle
2320     else fHistPiAPMonitorCuts[isSecd]->Fill(17);
2321
2322     if(dcaDaughters > fDCADaughtersAL )cutOKALambda = kFALSE;// dca daughters
2323     else  fHistPiAPMonitorCuts[isSecd]->Fill(18);
2324
2325     if( dcaPosToVertex < fDCADaughtersToVtxSmall || dcaNegToVertex < fDCADaughtersToVtxLarge)
2326       cutOKALambda = kFALSE;//dca pos and neg daughter to vertex
2327     else fHistPiAPMonitorCuts[isSecd]->Fill(19);
2328
2329     if(fRapCutV0 && fabs(rapAL) > fRap) cutOKALambda = kFALSE;//rapidity
2330     else fHistPiAPMonitorCuts[isSecd]->Fill(20);
2331
2332     if(fabs(0.497 - massK0s) < fExcludeK0sFromLambda) cutOKALambda = kFALSE;//exclude K0s from (A)Lambda
2333     /*
2334       if(chi2ALambdaC > fChiCutKf) cutOKALambda = kFALSE;
2335       else  fHistPiAPMonitorCuts[isSecd]->Fill(20);
2336     */
2337
2338     if(opAng < fOpengAngleDaughters && fabs(ptALambda) < fOpAngPtCut )  cutOKALambda = kFALSE;//opening angle cut
2339     else  fHistPiAPMonitorCuts[isSecd]->Fill(21);
2340
2341
2342     if((fArmCutL && qt>qtval) || alfa > -1.0*fAlfaCut) cutOKALambda = kFALSE; //cut in armenteros-podolanski diagram
2343     else  fHistPiAPMonitorCuts[isSecd]->Fill(22);
2344
2345     if(ptALambda > fPtTPCCut){
2346       if(fESDTrackCuts && fESDTrackCutsCharged){
2347         if(!fESDTrackCuts->AcceptTrack(trackPosTest)) cutOKALambda = kFALSE;//esd track cuts for pion
2348         if(!fESDTrackCutsCharged->AcceptTrack(trackNegTest)) cutOKALambda = kFALSE;//esd track cuts for antiproton
2349       }
2350       else  fHistPiAPMonitorCuts[isSecd]->Fill(5);
2351     }
2352     else{
2353       if(fESDTrackCutsLowPt){
2354         if(!fESDTrackCutsLowPt->AcceptTrack(trackPosTest))  cutOKALambda = kFALSE; //esd track cuts at low pt
2355         if(!fESDTrackCutsLowPt->AcceptTrack(trackNegTest))  cutOKALambda = kFALSE; //esd track cuts at low pt
2356       }
2357     }
2358
2359
2360     //--------------------- cut on PID: TPC dEdx ----------------------------//
2361     //-- dEdx --//
2362     Float_t nSigmaTPCtrackPosToPion = 0.0;
2363     Float_t nSigmaTPCtrackNegToPion = 0.0;
2364     Float_t nSigmaTPCtrackPosToProton = 0.0;
2365     Float_t nSigmaTPCtrackNegToProton = 0.0;
2366
2367
2368     if(fESDpid){
2369       nSigmaTPCtrackPosToPion   = fabs(fESDpid->NumberOfSigmasTPC(trackPos,AliPID::kPion));
2370       nSigmaTPCtrackNegToPion   = fabs(fESDpid->NumberOfSigmasTPC(trackNeg,AliPID::kPion));
2371       nSigmaTPCtrackPosToProton = fabs(fESDpid->NumberOfSigmasTPC(trackPos,AliPID::kProton));
2372       nSigmaTPCtrackNegToProton = fabs(fESDpid->NumberOfSigmasTPC(trackNeg,AliPID::kProton));
2373     }
2374
2375     Bool_t pipidEdx = kTRUE;
2376     Bool_t pipdEdx  = kTRUE;
2377     Bool_t piapdEdx = kTRUE;
2378
2379     Double_t tpcsigPos = trackPos->GetTPCsignal();
2380     Double_t tpcsigNeg = trackNeg->GetTPCsignal();
2381
2382     /*
2383       Double_t tpcsigNPos= trackPos->GetTPCsignalN();
2384       Double_t tpcsigNNeg= trackNeg->GetTPCsignalN();
2385     */
2386
2387     Double_t tpcMomPos = trackPos->GetInnerParam()->GetP();
2388     Double_t tpcMomNeg = trackNeg->GetInnerParam()->GetP();
2389
2390     //-- dedx cut --//
2391     if(fUsePID){
2392
2393       if(fabs(posDaughterP)<fPPIDcut && tpcsigPos < 5.0){//no zero dedx values!
2394         pipidEdx =kFALSE;//k0s
2395         piapdEdx =kFALSE;//antilambda
2396       }
2397       if(fabs(posDaughterP)<fPPIDcut && (nSigmaTPCtrackPosToProton > fNSigma || tpcsigPos < 5.0)) pipdEdx =kFALSE;//lambda
2398
2399       if(fabs(negDaughterP)<fPPIDcut &&  tpcsigNeg < 5.0){//no zero dedx values!
2400         pipidEdx =kFALSE;//k0s
2401         pipdEdx =kFALSE;//lambda
2402       }
2403       if(fabs(negDaughterP)<fPPIDcut && (nSigmaTPCtrackNegToProton > fNSigma || tpcsigNeg< 5.0)) piapdEdx =kFALSE;//antilambda
2404
2405     }
2406
2407     //-------------------- V0 ana -------------------------//
2408     //-- cut flags for furhter histos--//
2409     Bool_t k0sOK=kFALSE;
2410     Bool_t lambdaOK=kFALSE;
2411     Bool_t alambdaOK=kFALSE;
2412
2413     //--  Check for K0 --//
2414     if(ptK0s > fMinPt){//min pt cut
2415       if( cutOKK0s  && fillK0sMC && isSecd == 0 ){//cuts ok
2416         fHistDedxPiPlus[0]->Fill(tpcMomPos,tpcsigPos);
2417         fHistDedxPiMinus[0]->Fill(tpcMomNeg,tpcsigNeg);
2418         fHistPiPiMonitorCuts[0]->Fill(23);
2419         if(pipidEdx){//PID
2420           fHistPiPiMonitorCuts[0]->Fill(24);
2421           k0sOK=kTRUE;
2422           if(massK0s>0.25 && massK0s<0.75 ){//mass window
2423             if(!fMCMode){
2424               ptV0MC = ptK0s;
2425               declengthV0MC = dim2V0Radius;
2426             }
2427             fHistPiPiMonitorCuts[0]->Fill(25);
2428             fHistPiPiMass[0]->Fill(massK0s);
2429             fHistPiPiMassVSPt[0]->Fill(massK0s,ptK0s);
2430             fHistPiPiMassVSPtMCTruth[0]->Fill(massK0s,ptV0MC);
2431             fHistPiPiRadiusXY[0]->Fill(massK0s,posDaughterPhi);
2432             fHistPiPiCosPointAng[0]->Fill(massK0s,cosOPAng);
2433             if(massK0s>0.46 && massK0s<0.53)  fHistPiPiDecayLengthVsPt[0]->Fill(ptV0MC,dim2V0Radius);
2434             if(fMCMode)  fHistPiPiDecayLengthResolution[0]->Fill(declengthV0MC,dim2V0Radius);
2435             fHistPiPiDecayLengthVsCtau[0]->Fill(massK0s,ctTK0);
2436             fHistPiPiDecayLengthVsMass[0]->Fill(massK0s,dim2V0Radius);
2437
2438             fHistPiPiDCADaughters[0]->Fill(massK0s,dcaDaughters);
2439             fHistPiPiDCADaughterPosToPrimVtxVSMass[0]->Fill(massK0s,dcaPosToVertex);
2440             fHistPiPiDCAVSMass[0]->Fill(massK0s,dcaV0ToPrimVertex);
2441             // fHistPiPiPtDaughters[0]->Fill(posDaughterPt,negDaughterPt);
2442             fHistPiPiPtVSY[0]->Fill(rapK0s,ptK0s);
2443           }
2444           fHistArmenteros[0]->Fill(alfa,qt);
2445           fHistDedxSecPiPlus[0]->Fill(tpcMomPos,tpcsigPos);
2446           fHistDedxSecPiMinus[0]->Fill(tpcMomNeg,tpcsigNeg);
2447
2448           fHistV0RadiusZ[0]->Fill(dim2V0Radius,xr[2]);
2449           fHistV0RadiusXY[0]->Fill(xr[0],xr[1]);
2450           fHistV0RadiusZVSPt[0]->Fill(ptK0s,dim2V0Radius);
2451           fHistV0RadiusXYVSY[0]->Fill(rapK0s,dim2V0Radius);
2452
2453           //-- detector values --/
2454           fHistNclsITSPosK0[0]->Fill(nclsITSPos);
2455           fHistNclsITSNegK0[0]->Fill(nclsITSNeg);
2456           fHistNclsTPCPosK0[0]->Fill(nclsTPCPos);
2457           fHistNclsTPCNegK0[0]->Fill(nclsTPCNeg);
2458           fHistChi2PerNclsITSPosK0[0]->Fill(chi2PerClusterITSPos);
2459           fHistChi2PerNclsITSNegK0[0]->Fill(chi2PerClusterITSNeg);
2460
2461
2462           fHistNclsITSNeg[0]->Fill(massK0s,nclsITSPos);
2463           fHistNclsTPCNeg[0]->Fill(massK0s,nclsTPCPos);
2464           fHistChi2PerNclsITSNeg[0]->Fill(massK0s,chi2PerClusterITSPos);
2465           fHistNCRowsTPCNeg[0]->Fill(massK0s,crossedRowsTPCPos);
2466
2467
2468           fHistNclsITSNeg[1]->Fill(massK0s,nclsITSNeg);
2469           fHistNclsTPCNeg[1]->Fill(massK0s,nclsTPCNeg);
2470           fHistChi2PerNclsITSNeg[1]->Fill(massK0s,chi2PerClusterITSNeg);
2471           fHistNCRowsTPCNeg[1]->Fill(massK0s,crossedRowsTPCNeg);
2472
2473           fHistRatioFoundOverFinableTPCK0[0]->Fill(massK0s,ratioFoFi);
2474           fHistRatioFoundOverFinableTPCK0[1]->Fill(massK0s,ratioFoFiPos);
2475
2476           fHistNclsITS[1]->Fill(nclsITSPos,nclsITSNeg);
2477           fHistNclsTPC[1]->Fill(crossedRowsTPCNeg,nclsTPCNeg);
2478         }
2479       }
2480     }
2481
2482     //-----  Check for Lambda -------//
2483     if(ptLambda >fMinPt){//min pt cut
2484       if(cutOKLambda && fillLambdaMC){// cuts ok
2485         fHistDedxProt[isSecd]->Fill(tpcMomPos,tpcsigPos);
2486         fHistDedxPiMinus[isSecd]->Fill(tpcMomNeg,tpcsigNeg);
2487         fHistPiPMonitorCuts[isSecd]->Fill(23);
2488         if(pipdEdx){//PID
2489           fHistPiPMonitorCuts[isSecd]->Fill(24);
2490           lambdaOK=kTRUE;
2491           if( massLambda>1.05 && massLambda<1.25 ){//mass window
2492             if(!fMCMode) {
2493               ptV0MC = ptLambda;
2494               declengthV0MC = dim2V0Radius;
2495             }
2496             fHistPiPMonitorCuts[isSecd]->Fill(25);
2497             fHistPiPMass[isSecd]->Fill(massLambda);
2498             fHistPiPMassVSPtMCTruth[isSecd]->Fill(massLambda,ptV0MC);
2499             fHistPiPMassVSPt[isSecd]->Fill(massLambda,ptLambda);
2500             fHistPiPRadiusXY[isSecd]->Fill(massLambda,posDaughterPhi);
2501             fHistPiPCosPointAng[isSecd]->Fill(massLambda,cosOPAng);
2502             fHistPiPPtVSY[isSecd]->Fill(rapL,ptLambda);
2503             //   fHistPiPPtDaughters[isSecd]->Fill(posDaughterPt,negDaughterPt);
2504             fHistPiPDCADaughters[isSecd]->Fill(massLambda,dcaDaughters);
2505             fHistPiPDCAVSMass[isSecd]->Fill(massLambda,dcaV0ToPrimVertex);
2506             fHistPiPDCADaughterPosToPrimVtxVSMass[isSecd]->Fill(massLambda,dcaPosToVertex);
2507             if( massLambda>1.108 && massLambda<1.123 ) fHistPiPDecayLengthVsPt[isSecd]->Fill(ptV0MC,dim2V0Radius);
2508             if(fMCMode)  fHistPiPDecayLengthResolution[isSecd]->Fill(declengthV0MC,dim2V0Radius);
2509             fHistPiPDecayLengthVsMass[isSecd]->Fill(massLambda,dim2V0Radius);
2510             fHistPiPDecayLengthVsCtau[isSecd]->Fill(massLambda,ctTL);
2511
2512             //-- secondaries --//
2513             if(isSecd == 1){
2514               if(fabs(pdgMother) == 3112 || fabs(pdgMother) == 3114 || fabs(pdgMother) == 3222 || fabs(pdgMother) == 3224 || fabs(pdgMother) == 3214 ){
2515                 fHistPiPMassVSPtSecSigma[1]->Fill(massLambda,ptLambda);
2516               }
2517             }
2518             if(pdgMother == 3322 || pdgMother == 3312){//Xi0 and xi minus
2519               fHistPiPCosPointAngXiVsPt->Fill(ptLambda,cosOPAng);
2520               fHistPiPMassVSPtSecXi[1]->Fill(massLambda,ptLambda);
2521               fHistPiPMassVSPtSecXiMCTruth->Fill(massLambda,ptV0MC);
2522               fHistPiPMassVSYSecXi[1]->Fill(massLambda,rapL);
2523               fHistPiPXi0PtVSLambdaPt[1]->Fill(ptLambda,ptXiMother);
2524             }
2525             if(pdgMother == 3334){//Omega
2526               fHistPiPMassVSPtSecOmega[1]->Fill(massLambda,ptLambda);
2527               fHistPiPMassVSPtSecOmegaMCTruth->Fill(massLambda,ptV0MC);
2528               fHistPiPOmegaPtVSLambdaPt[1]->Fill(ptLambda,ptXiMother);
2529             }
2530           }
2531         }
2532
2533         if(ptLambda > 0.4) fHistArmenteros[isSecd]->Fill(alfa,qt);
2534         fHistV0RadiusZ[isSecd]->Fill(dim2V0Radius,xr[2]);
2535         fHistV0RadiusXY[isSecd]->Fill(xr[0],xr[1]);
2536         fHistV0RadiusZVSPt[isSecd]->Fill(ptLambda,dim2V0Radius);
2537         fHistV0RadiusXYVSY[isSecd]->Fill(rapL,dim2V0Radius);
2538         fHistDedxSecProt[isSecd]->Fill(tpcMomPos,tpcsigPos);
2539         fHistDedxSecPiMinus[isSecd]->Fill(tpcMomNeg,tpcsigNeg);
2540
2541         //-- detector values --//
2542         fHistNclsITSPosL[isSecd]->Fill(nclsITSPos);
2543         fHistNclsITSNegL[isSecd]->Fill(nclsITSNeg);
2544         fHistNclsTPCPosL[isSecd]->Fill(nclsTPCPos);
2545         fHistNclsTPCNegL[isSecd]->Fill(nclsTPCNeg);
2546         fHistChi2PerNclsITSPosL[isSecd]->Fill(chi2PerClusterITSPos);
2547         fHistChi2PerNclsITSNegL[isSecd]->Fill(chi2PerClusterITSNeg);
2548
2549         fHistNclsITSPos[0]->Fill(massLambda,nclsITSPos);
2550         fHistNclsTPCPos[0]->Fill(massLambda,nclsTPCPos);
2551         fHistChi2PerNclsITSPos[0]->Fill(massLambda,chi2PerClusterITSPos);
2552         fHistNCRowsTPCPos[0]->Fill(massLambda,crossedRowsTPCPos);
2553
2554         fHistNclsITSPos[1]->Fill(massLambda,nclsITSNeg);
2555         fHistNclsTPCPos[1]->Fill(massLambda,nclsTPCNeg);
2556         fHistChi2PerNclsITSPos[1]->Fill(massLambda,chi2PerClusterITSNeg);
2557         fHistNCRowsTPCPos[1]->Fill(massLambda,crossedRowsTPCNeg);
2558
2559         fHistRatioFoundOverFinableTPCL[0]->Fill(massLambda,ratioFoFi);
2560         fHistRatioFoundOverFinableTPCL[1]->Fill(massLambda,ratioFoFiPos);
2561
2562         fHistNclsITS[0]->Fill(nclsITSPos,nclsITSNeg);
2563         fHistNclsTPC[0]->Fill(crossedRowsTPCNeg,nclsTPCNeg);
2564       }
2565     }
2566
2567     if(ptALambda >fMinPt){//min pt cut
2568       //-- Check for AntiLambda --//
2569       if(cutOKALambda && fillALambdaMC){// cuts ok
2570         fHistDedxAProt[isSecd]->Fill(tpcMomNeg,tpcsigNeg);
2571         fHistDedxPiPlus[isSecd]->Fill(tpcMomPos,tpcsigPos);
2572         fHistPiAPMonitorCuts[isSecd]->Fill(23);
2573         if(piapdEdx){//PID
2574           fHistPiAPMonitorCuts[isSecd]->Fill(24);
2575           alambdaOK=kTRUE;
2576           if( massALambda>1.05 && massALambda<1.25){// mass window
2577             if(!fMCMode) {
2578               ptV0MC = ptALambda;
2579               declengthV0MC = dim2V0Radius;
2580             }
2581             fHistPiAPMonitorCuts[isSecd]->Fill(25);
2582             fHistPiAPMass[isSecd]->Fill(massALambda);
2583             fHistPiAPMassVSPtMCTruth[isSecd]->Fill(massALambda,ptV0MC);
2584             fHistPiAPMassVSPt[isSecd]->Fill(massALambda,ptALambda);
2585             fHistPiAPRadiusXY[isSecd]->Fill(massALambda,posDaughterPhi);
2586             fHistPiAPCosPointAng[isSecd]->Fill(massALambda,cosOPAng);
2587             fHistPiAPPtVSY[isSecd]->Fill(rapAL,ptALambda);
2588             //  fHistPiAPPtDaughters[isSecd]->Fill(posDaughterPt,negDaughterPt);
2589             fHistPiAPDCADaughters[isSecd]->Fill(massALambda,dcaDaughters);
2590             fHistPiAPDCAVSMass[isSecd]->Fill(massALambda,dcaV0ToPrimVertex);
2591             fHistPiAPDCADaughterPosToPrimVtxVSMass[isSecd]->Fill(massALambda,dcaPosToVertex);
2592             if( massALambda>1.108 && massALambda<1.123 )  fHistPiAPDecayLengthVsPt[isSecd]->Fill(ptV0MC,dim2V0Radius);
2593             if(fMCMode)  fHistPiAPDecayLengthResolution[isSecd]->Fill(declengthV0MC,dim2V0Radius);
2594             fHistPiAPDecayLengthVsMass[isSecd]->Fill(massALambda,dim2V0Radius);
2595             fHistPiAPDecayLengthVsCtau[isSecd]->Fill(massALambda,ctTAL);
2596
2597             //-- secondaries --//
2598             if(isSecd == 1){
2599               if(fabs(pdgMother) == 3112 || fabs(pdgMother) == 3114 || fabs(pdgMother) == 3222 || fabs(pdgMother) == 3224 || fabs(pdgMother) == 3214 ){
2600                 fHistPiAPMassVSPtSecSigma[1]->Fill(massLambda,ptLambda);
2601
2602               }
2603             }
2604             if(pdgMother == -3322 || pdgMother == -3312){//Xi0 and xiplus
2605               fHistPiAPCosPointAngXiVsPt->Fill(ptALambda,cosOPAng);
2606               fHistPiAPMassVSPtSecXi[1]->Fill(massALambda,ptALambda);
2607               fHistPiAPMassVSPtSecXiMCTruth->Fill(massALambda,ptV0MC);
2608               fHistPiAPMassVSYSecXi[1]->Fill(massALambda,rapAL);
2609               fHistPiAPXi0PtVSLambdaPt[1]->Fill(ptALambda,ptXiMother);
2610             }
2611             if(pdgMother == -3334){//Omega
2612               fHistPiAPMassVSPtSecOmega[1]->Fill(massALambda,ptALambda);
2613               fHistPiAPMassVSPtSecOmegaMCTruth->Fill(massALambda,ptV0MC);
2614               fHistPiAPOmegaPtVSLambdaPt[1]->Fill(ptALambda,ptXiMother);
2615             }
2616           }
2617
2618         }
2619         if(ptALambda > 0.4) fHistArmenteros[isSecd]->Fill(alfa,qt);
2620         fHistDedxSecAProt[isSecd]->Fill(tpcMomNeg,tpcsigNeg);
2621         fHistDedxSecPiPlus[isSecd]->Fill(tpcMomPos,tpcsigPos);
2622         fHistV0RadiusZ[isSecd]->Fill(dim2V0Radius,xr[2]);
2623         fHistV0RadiusXY[isSecd]->Fill(xr[0],xr[1]);
2624         fHistV0RadiusZVSPt[isSecd]->Fill(ptALambda,dim2V0Radius);
2625         fHistV0RadiusXYVSY[isSecd]->Fill(rapAL,dim2V0Radius);
2626       }
2627     }
2628
2629
2630     //-- fill detector histos general --//
2631     if(lambdaOK || alambdaOK || k0sOK){
2632       fHistPiPiEtaDReco[1]->Fill(posDaughterEta,eta00);
2633       fHistPiPEtaDReco[1]->Fill(negDaughterEta,eta01);
2634     }
2635
2636     /*
2637     //-- AliKFParticle --//
2638     if (negPiKF) delete negPiKF; negPiKF=NULL;
2639     if (posPiKF) delete posPiKF; posPiKF=NULL;
2640     if (posPKF) delete posPKF; posPKF=NULL;
2641     if (negAPKF) delete negAPKF; negAPKF=NULL;
2642     */
2643
2644
2645   }//end V0 reco loop
2646 }
2647
2648 //________________________________________________________________________
2649
2650 Int_t AliAnalysisTaskV0ForRAA::CalculateCentralityBin(){
2651   //find centrality bin for centrality selection
2652
2653   if (fUseCentrality == 0) return -1;
2654
2655   AliCentrality *esdCentrality = fESD->GetCentrality();
2656
2657   Float_t centralityVZERO  = esdCentrality->GetCentralityPercentile("V0M");
2658   Float_t centralitySPD    = esdCentrality->GetCentralityPercentile("CL1");
2659
2660   Int_t centralityVZEROBin = -1;
2661   Int_t centralitySPDBin   = -1;
2662
2663   //-- SPD centrality --//
2664   if ( fUseCentrality == 2 ){
2665     if      ( centralitySPD >=  0. && centralitySPD <   5.) centralitySPDBin =  0;
2666     else if ( centralitySPD >=  5. && centralitySPD <  10.) centralitySPDBin =  5;
2667     else if ( centralitySPD >= 10. && centralitySPD <  20.) centralitySPDBin = 10;
2668     else if ( centralitySPD >= 20. && centralitySPD <  30.) centralitySPDBin = 20;
2669     else if ( centralitySPD >= 30. && centralitySPD <  40.) centralitySPDBin = 30;
2670     else if ( centralitySPD >= 40. && centralitySPD <  50.) centralitySPDBin = 40;
2671     else if ( centralitySPD >= 50. && centralitySPD <  60.) centralitySPDBin = 50;
2672     else if ( centralitySPD >= 60. && centralitySPD <  70.) centralitySPDBin = 60;
2673     else if ( centralitySPD >= 70. && centralitySPD <  80.) centralitySPDBin = 70;
2674     else if ( centralitySPD >= 80. && centralitySPD <  90.) centralitySPDBin = 80;
2675     else if ( centralitySPD >= 90. && centralitySPD <  99.) centralitySPDBin = 90;
2676     else if ( centralitySPD >= 99. ) centralitySPDBin = 100;
2677     else if ( fabs(centralitySPD)< 0.0001 ) centralitySPDBin = 100;
2678     return centralitySPDBin;
2679   }
2680   //-- V0 centrality --//
2681   if ( fUseCentrality == 1 ){
2682     if      ( centralityVZERO >  0. && centralityVZERO <   5.) centralityVZEROBin =  0;
2683     else if ( centralityVZERO >=  5. && centralityVZERO <  10.) centralityVZEROBin =  5;
2684     else if ( centralityVZERO >= 10. && centralityVZERO <  20.) centralityVZEROBin = 10;
2685     else if ( centralityVZERO >= 20. && centralityVZERO <  30.) centralityVZEROBin = 20;
2686     else if ( centralityVZERO >= 30. && centralityVZERO <  40.) centralityVZEROBin = 30;
2687     else if ( centralityVZERO >= 40. && centralityVZERO <  50.) centralityVZEROBin = 40;
2688     else if ( centralityVZERO >= 50. && centralityVZERO <  60.) centralityVZEROBin = 50;
2689     else if ( centralityVZERO >= 60. && centralityVZERO <  70.) centralityVZEROBin = 60;
2690     else if ( centralityVZERO >= 70. && centralityVZERO <  80.) centralityVZEROBin = 70;
2691     else if ( centralityVZERO >= 80. && centralityVZERO <  90.) centralityVZEROBin = 80;
2692     else if ( centralityVZERO >= 90. && centralityVZERO <  99.) centralityVZEROBin = 90;
2693     else if ( centralityVZERO >= 99. ) centralityVZEROBin = 100;
2694     else if ( fabs(centralityVZERO)< 0.0001 ) centralityVZEROBin = 100;
2695     return centralityVZEROBin;
2696   }
2697   return -1;
2698
2699 }
2700
2701 //__________________________________________________________________________________________________________
2702 Bool_t  AliAnalysisTaskV0ForRAA::GetMCTruthPartner(AliESDtrack *pos,AliESDtrack *neg,Int_t id0,Int_t id1){
2703   //-- get daughter label and check it --//
2704   Int_t labelP = fabs(pos->GetLabel());
2705   Int_t labelN = fabs(neg->GetLabel());
2706
2707   if (labelN==labelP)  return kFALSE;
2708
2709   if(fMCTruthMode){
2710     if ((labelP!=id0) && (labelP!=id1))  return kFALSE;
2711     if ((labelN!=id0) && (labelN!=id1))  return kFALSE;
2712   }
2713
2714   return kTRUE;
2715 }
2716