PWGLF/STRANGENESS/Correlations/AliAnalysisTaskLambdaOverK0sJets.cxx

   1 /*************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: X. Sanchez Castro                                              *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include <TCanvas.h>
  17 #include <TTree.h>
  18 #include <TFile.h>
  19 #include <TH1F.h>
  20 #include <TH2F.h>
  21 #include <TH3F.h>
  22 #include <TPDGCode.h>
  23 #include <TDatabasePDG.h>
  24 #include <TClonesArray.h>
  25 #include <TROOT.h>
  26
  27 #include "AliOADBContainer.h"
  28
  29 #include "AliAODMCHeader.h"
  30 #include "AliAODMCParticle.h"
  31 #include "AliGenHijingEventHeader.h"
  32
  33 #include "AliAODEvent.h"
  34 #include "AliAODv0.h"
  35 #include "AliAODcascade.h"
  36
  37 #include "AliCFContainer.h"
  38 #include "AliCentrality.h"
  39
  40 #include "AliPID.h"
  41 #include "AliPIDResponse.h"
  42 #include "AliAODPid.h"
  43
  44 #include "AliInputEventHandler.h"
  45 #include "AliAnalysisManager.h"
  46
  47 #include "AliAnalysisTaskLambdaOverK0sJets.h"
  48
  49 extern TROOT *gROOT;
  50
  51
  52 ClassImp(AliAnalysisTaskLambdaOverK0sJets)
  53
  54 // Global variables:
  55 static Int_t    nbins = 100;                 // Number of bins for l, pt, mass for V0
  56 static Int_t    nbinsPhi = 120;              // Number of bins for Phi
  57 static Int_t    nbinsdPhi = 42;              // Number of bins for dPhi
  58 static Int_t    nbinsdEta = 40;              // Number of bins for dEta
  59 static Int_t    nbinPtLP = 200;
  60 static Int_t    nbinsVtx = 20;
  61
  62 static Double_t pMin = 0.0;                  // Lower cut for transverse momentum
  63 static Double_t pMax = 10.;                  // Max cut for transverse momentum for V0
  64 static Double_t ptMaxLP = 50.;               // Max cut for transverse momentum LP
  65
  66 static Double_t lMin = 0.0;                  // Limits in the histo for fidutial volume
  67 static Double_t lMax = 100.;                 // Limits in the fidutial volume
  68
  69 //
  70 //
  71 //
  72
  73 AliAnalysisTaskLambdaOverK0sJets::AliAnalysisTaskLambdaOverK0sJets(const char *name) :
  74   AliAnalysisTaskSE(name),
  75
  76   fAOD(0),
  77   fIsMC(kFALSE),
  78   fUsePID(kFALSE),
  79   fCentMin(0.),
  80   fCentMax(90.),
  81   fTrigPtMin(8.),
  82   fTrigPtMax(20.),
  83   fTrigEtaMax(0.8),
  84   fSeparateInjPart(kTRUE),
  85   fEndOfHijingEvent(-1),
  86   fPIDResponse(0),
  87
  88   fMinPtDaughter(0.160),
  89   fMaxEtaDaughter(0.8),
  90   fMaxDCADaughter(1.0),
  91   fYMax(0.5),
  92   fDCAToPrimVtx(0.1),
  93   fMinCPA(0.998),
  94   fNSigma(3.0),
  95   fMinCtau(0.),
  96   fMaxCtau(3.),
  97
  98   fOutput(0),
  99
 100   fEvents(0),
 101   fCentrality(0),
 102   fPrimaryVertexX(0),
 103   fPrimaryVertexY(0),
 104   fPrimaryVertexZ(0),
 105   fNumberPileUp(0),
 106   fCentMult(0),
 107   fdEdx(0),
 108   fdEdxPid(0),
 109
 110   fTriggerMCPtCent(0),
 111   fTriggerPtCent(0),
 112   fTriggerEtaPhi(0),
 113   fCheckTriggerFromV0Daug(0),
 114   fTriggerComingFromDaug(0),
 115   fTriggerIsV0(0),
 116
 117   fInjectedParticles(0),
 118
 119   fK0sMCPt(0),
 120   fK0sMCPtRap(0),
 121   fK0sMCPtEta(0),
 122   fK0sMCPtLt(0),
 123   fK0sAssocPt(0),
 124   fK0sAssocPtLt(0),
 125   fK0sAssocPtLtArm(0),
 126   fK0sAssocPtRap(0),
 127   fK0sAssocPtEta(0),
 128
 129   fLambdaMCPt(0),
 130   fLambdaMCPtRap(0),
 131   fLambdaMCPtEta(0),
 132   fLambdaMCPtLt(0),
 133   fLambdaAssocPt(0),
 134   fLambdaAssocPtLt(0),
 135   fLambdaAssocPtLtArm(0),
 136   fLambdaAssocPtRap(0),
 137   fLambdaAssocPtEta(0),
 138
 139   fHistArmenterosPodolanski(0),
 140   fHistArmPodBckg(0),
 141
 142   fK0sMass(0),
 143   fK0sPtLtSB(0),
 144   fK0sPtvsEta (0),
 145   fK0sPtvsRap(0),
 146   fK0sEtaPhi(0),
 147   fK0sMassPtPhi(0),
 148
 149   fK0sMassPtvsPtL(0),
 150   fK0sSiPtL(0),
 151   fK0sDaughtersPt(0),
 152   fK0sdPhiPtAssocPtL(0),
 153   fK0sDCADaugToPrimVtx(0),
 154
 155   fK0sBckgDecLength(0),
 156   fK0sBckgDCADaugToPrimVtx(0),
 157   fK0sdEdxPosDaug(0),
 158   fK0sdEdxNegDaug(0),
 159   fK0sBckgEtaPhi(0),
 160   fK0sBckgPhiRadio(0),
 161   fK0sBckgDCANegDaugToPrimVtx(0),
 162   fK0sBckgDCAPosDaugToPrimVtx(0),
 163   fK0sMassCascade(0),
 164
 165   fLambdaMass(0),
 166   fLambdaPtLtSB(0),
 167   fLambdaPtvsEta(0),
 168   fLambdaPtvsRap(0),
 169   fLambdaEtaPhi(0),
 170   fLambdaMassPtPhi(0),
 171
 172   fLambdadEdx(0),
 173   fCPA(0),
 174   fDCA(0),
 175
 176   fLambdaMassPtvsPtL(0),
 177   fLambdaSiPtL(0),
 178   fLambdaDaughtersPt(0),
 179   fLambdadPhiPtAssocPtL(0),
 180   fLambdaDCADaugToPrimVtx(0),
 181
 182   fLambdaBckgDecLength(0),
 183   fLambdaBckgDCADaugToPrimVtx(0),
 184   fLambdadEdxPosDaug(0),
 185   fLambdadEdxNegDaug(0),
 186   fLambdaBckgEtaPhi(0),
 187   fLambdaBckgPhiRadio(0),
 188   fLambdaBckgDCANegDaugToPrimVtx(0),
 189   fLambdaBckgDCAPosDaugToPrimVtx(0),
 190   fLambdaMassCascade(0)
 191
 192
 193 {
 194   // Dummy Constructor
 195
 196   for (Int_t i=0; i<kN1; i++){
 197
 198     fK0sdPhidEtaMC[i] = 0;
 199     fK0sdPhidEtaMCCent[i] = 0;
 200
 201     fK0sdPhidEtaPtL[i] = 0;
 202     fK0sdPhidEtaPtLCent[i] = 0;
 203     fK0sdPhidEtaPtLBckg[i] = 0;
 204     fK0sdPhidEtaPtLCentBckg[i] = 0;
 205
 206     fK0sdPhidEtaPtL2[i] = 0;
 207     fK0sdPhidEtaPtLCent2[i] = 0;
 208     fK0sdPhidEtaPtLBckg2[i] = 0;
 209     fK0sdPhidEtaPtLCentBckg2[i] = 0;
 210
 211     fLambdadPhidEtaMC[i] = 0;
 212     fLambdadPhidEtaMCCent[i] = 0;
 213
 214     fLambdadPhidEtaPtL[i] = 0;
 215     fLambdadPhidEtaPtLCent[i] = 0;
 216     fLambdadPhidEtaPtLBckg[i] = 0;
 217     fLambdadPhidEtaPtLCentBckg[i] = 0;
 218
 219     fLambdadPhidEtaPtL2[i] = 0;
 220     fLambdadPhidEtaPtLCent2[i] = 0;
 221     fLambdadPhidEtaPtLBckg2[i] = 0;
 222     fLambdadPhidEtaPtLCentBckg2[i] = 0;
 223   }
 224
 225   // Constructor. Initialization of pointers
 226   DefineOutput(1, TList::Class());
 227
 228 }
 229
 230 //___________________________________________________________________________________________
 231
 232 void AliAnalysisTaskLambdaOverK0sJets::UserCreateOutputObjects()
 233 {
 234   // Creating the histograms that are needed for the output
 235
 236   fOutput = new TList();
 237   fOutput->SetOwner();
 238
 239   // ====== General characteristics of the event and tracks ====== //
 240
 241   // Counter for the number of events in each step:
 242   fEvents=new TH1F("fEvents","Number of events",14,-0.5,13.5);
 243   fEvents->GetXaxis()->SetBinLabel(1,"calls to UserExec()");
 244   fEvents->GetXaxis()->SetBinLabel(2,"AOD available");
 245   fEvents->GetXaxis()->SetBinLabel(3,"CINT1B");
 246   fEvents->GetXaxis()->SetBinLabel(4,"V0M Cent");
 247   fEvents->GetXaxis()->SetBinLabel(5,"Vtx > 3 part");
 248   fEvents->GetXaxis()->SetBinLabel(6,"|VtxZ| < 10 cm");
 249   fEvents->GetXaxis()->SetBinLabel(7,"Mult && Cent");
 250   fEvents->GetXaxis()->SetBinLabel(8,"Bad ID Trigger");
 251   fEvents->GetXaxis()->SetBinLabel(9,"V0 is LP");
 252   fEvents->GetXaxis()->SetBinLabel(10,"Trigger is V0 daug");
 253   fEvents->GetXaxis()->SetBinLabel(11,"Trigger is V0 daug && 2nd check");
 254   fEvents->GetXaxis()->SetBinLabel(12,"Triggered");
 255   fEvents->GetXaxis()->SetBinLabel(13,"NOT Triggered");
 256   fEvents->GetYaxis()->SetTitle("Counts");
 257   fOutput->Add(fEvents);
 258
 259   // Centrality:
 260   fCentrality = new TH1F("fCentrality","Centrality",100,0.,100.);
 261   fCentrality->GetXaxis()->SetTitle("Centrality (%)");
 262   fCentrality->GetYaxis()->SetTitle("Events");
 263   fOutput->Add(fCentrality);
 264
 265   // Primary Vertex:
 266   fPrimaryVertexX = new TH1F("fPrimaryVertexX", "Primary Vertex Position X",100,-0.5,0.5);
 267   fPrimaryVertexX->GetXaxis()->SetTitle("Primary Vertex Position X (cm)");
 268   fPrimaryVertexX->GetYaxis()->SetTitle("Events");
 269   fOutput->Add(fPrimaryVertexX);
 270
 271   fPrimaryVertexY = new TH1F("fPrimaryVertexY", "Primary Vertex Position Y",100,-0.5,0.5);
 272   fPrimaryVertexY->GetXaxis()->SetTitle("Primary Vertex Position Y (cm)");
 273   fPrimaryVertexY->GetYaxis()->SetTitle("Events");
 274   fOutput->Add(fPrimaryVertexY);
 275
 276   fPrimaryVertexZ = new TH1F("fPrimaryVertexZ", "Primary Vertex Position Z",200,-20,20);
 277   fPrimaryVertexZ->GetXaxis()->SetTitle("Primary Vertex Position Z (cm)");
 278   fPrimaryVertexZ->GetYaxis()->SetTitle("Events");
 279   fOutput->Add(fPrimaryVertexZ);
 280
 281   // Pile-up
 282   fNumberPileUp = new TH2F("fNumberPileUp","Pile up: SPD vs TRacks",5,-0.5,4.5,5,-0.5,4.5);
 283   fNumberPileUp->GetXaxis()->SetTitle("Pile-up wi the SPD");
 284   fNumberPileUp->GetYaxis()->SetTitle("Pile-up wi tracks");
 285   fOutput->Add(fNumberPileUp);
 286
 287   // Multiplicity:
 288   fCentMult = new TH2F("fCentMult","Centrality vs Multiplicity",100,0.,100.,1100,0.,3300);
 289   fCentMult->GetXaxis()->SetTitle("Centrality");
 290   fCentMult->GetYaxis()->SetTitle("Multiplicity");
 291   fOutput->Add(fCentMult);
 292
 293   // Energy loss:
 294   fdEdx = new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
 295   fdEdx->GetXaxis()->SetTitle("TPC Momentum (GeV/c)");
 296   fdEdx->GetYaxis()->SetTitle("a.u.");
 297   fOutput->Add(fdEdx);
 298
 299   fdEdxPid = new TH2F("fdEdxPid","dE/dx with PID (protons)",50,0.2,3,50,0.,6.);
 300   fdEdxPid->GetXaxis()->SetTitle("TPC Momentum (GeV/c)");
 301   fdEdxPid->GetYaxis()->SetTitle("a.u.");
 302   fOutput->Add(fdEdxPid);
 303
 304
 305   // ====== Trigger Particle characteristics ====== //
 306
 307   // MC: Pt Trigger particle vs centrality:
 308   fTriggerMCPtCent = new TH3F("fTriggerMCPtCent","Trigger particle MC",2*nbinPtLP,pMin,2*ptMaxLP,100,-1.,1.,100,0.,100.);
 309   fTriggerMCPtCent->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 310   fTriggerMCPtCent->GetYaxis()->SetTitle("#eta");
 311   fTriggerMCPtCent->GetZaxis()->SetTitle("centrality (%)");
 312   fOutput->Add(fTriggerMCPtCent);
 313
 314   // Pt Trigger particle vs centrality:
 315   fTriggerPtCent = new TH2F("fTriggerPtCent","Trigger particle",nbinPtLP,pMin,ptMaxLP,100,0.,100.);
 316   fTriggerPtCent->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 317   fTriggerPtCent->GetYaxis()->SetTitle("centrality (%)");
 318   fOutput->Add(fTriggerPtCent);
 319
 320   // Phi vs pseudorapidity:
 321   fTriggerEtaPhi = new TH2F("fTriggerEtaPhi","Trigger particle",nbinsPhi,0.,2.*TMath::Pi(),100,-1.,1.);
 322   fTriggerEtaPhi->GetXaxis()->SetTitle("#phi (rad)");
 323   fTriggerEtaPhi->GetYaxis()->SetTitle("#eta");
 324   fOutput->Add(fTriggerEtaPhi);
 325
 326   // Check if Trigger particle comes from a V0 daughter:
 327   fCheckTriggerFromV0Daug =
 328     new TH1F("fCheckTriggerFromV0Daug","Trigger particle from a V0 daughter",7,-2.5,4.5);
 329   fCheckTriggerFromV0Daug->GetXaxis()->SetTitle("Flag");
 330   fCheckTriggerFromV0Daug->GetXaxis()->SetBinLabel(2,"NOT V0 daug & 2nd check");
 331   fCheckTriggerFromV0Daug->GetXaxis()->SetBinLabel(3,"NOT V0 daug");
 332   fCheckTriggerFromV0Daug->GetXaxis()->SetBinLabel(4,"V0 daug");
 333   fCheckTriggerFromV0Daug->GetXaxis()->SetBinLabel(5,"V0 daug & V0 LP");
 334   fCheckTriggerFromV0Daug->GetXaxis()->SetBinLabel(6,"V0 daug & V0 LP & 2nd Check");
 335   fCheckTriggerFromV0Daug->GetYaxis()->SetTitle("Counts");
 336   fOutput->Add(fCheckTriggerFromV0Daug);
 337
 338   fTriggerComingFromDaug = new TH1F("fTriggerComingFromDaug","Trigger particle from a V0 daughter",240, 0, 12);
 339   fTriggerComingFromDaug->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 340   fTriggerComingFromDaug->GetYaxis()->SetTitle("Counts");
 341   fOutput->Add(fTriggerComingFromDaug);
 342
 343   fTriggerIsV0 = new TH1F("fTriggerIsV0","V0 candidate is a LP",nbinPtLP,pMin,ptMaxLP);
 344   fTriggerIsV0->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 345   fTriggerIsV0->GetYaxis()->SetTitle("Counts");
 346   fOutput->Add(fTriggerIsV0);
 347
 348
 349   // ====== MC-true and  MC-Association information ====== //
 350   if(fIsMC){
 351
 352   fInjectedParticles = new TH1F("fInjectedParticles","Injected particles",2,0.,2.);
 353   fInjectedParticles->GetXaxis()->SetBinLabel(1,"Injected");
 354   fInjectedParticles->GetXaxis()->SetBinLabel(2,"Natural");
 355   fInjectedParticles->GetYaxis()->SetTitle("Counts");
 356   fOutput->Add(fInjectedParticles);
 357
 358   // K0s MC-true:
 359   fK0sMCPt       = new TH1F("fK0sMCPt", "K^{0}_{S} MC",nbins,pMin,pMax);
 360   fK0sMCPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 361   fOutput->Add(fK0sMCPt);
 362
 363   fK0sMCPtRap    = new TH2F("fK0sMCPtRap", "K^{0}_{S} MC",nbins,pMin,pMax,30,-1.5,1.5);
 364   fK0sMCPtRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 365   fK0sMCPtRap->GetYaxis()->SetTitle("y");
 366   fOutput->Add(fK0sMCPtRap);
 367
 368   fK0sMCPtEta    = new TH2F("fK0sMCPtEta", "K^{0}_{S} MC",nbins,pMin,pMax,30,-1.5,1.5);
 369   fK0sMCPtEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 370   fK0sMCPtEta->GetYaxis()->SetTitle("#eta");
 371   fOutput->Add(fK0sMCPtEta);
 372
 373   fK0sMCPtLt =
 374     new TH3F("fK0sMCPtLt", "K^{0}_{S} MC",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 375   fK0sMCPtLt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 376   fK0sMCPtLt->GetYaxis()->SetTitle("L_{T} (cm)");
 377   fK0sMCPtLt->GetZaxis()->SetTitle("centrality");
 378   fOutput->Add(fK0sMCPtLt);
 379
 380   // K0s MC-Association:
 381   fK0sAssocPt =
 382     new TH1F("fK0sAssocPt","K^{0}_{S} Assoc: L_{T} vs p_{T}",nbins,pMin,pMax);
 383   fK0sAssocPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 384   fOutput->Add(fK0sAssocPt);
 385
 386   /*
 387   fK0sAssocPtLt =
 388     new TH3F("fK0sAssocPtLt","K^{0}_{S} Assoc: L_{T} vs p_{T}",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 389   fK0sAssocPtLt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 390   fK0sAssocPtLt->GetYaxis()->SetTitle("L_{T} (cm)");
 391   fK0sAssocPtLt->GetZaxis()->SetTitle("centrality");
 392   fOutput->Add(fK0sAssocPtLt);
 393   */
 394
 395   fK0sAssocPtLt =
 396     new TH3F("fK0sAssocPtLt","K^{0}_{S} Assoc:  p_{T} vs y vs centrlaity",nbins,pMin,pMax,30,-1.5,1.5,100,0.,100.);
 397   fK0sAssocPtLt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 398   fK0sAssocPtLt->GetYaxis()->SetTitle("y");
 399   fK0sAssocPtLt->GetZaxis()->SetTitle("centrality");
 400   fOutput->Add(fK0sAssocPtLt);
 401
 402   /*
 403   fK0sAssocPtLtArm =
 404     new TH3F("fK0sAssocPtLtArm","K^{0}_{S} Assoc: L_{T} vs p_{T}",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 405   fK0sAssocPtLtArm->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 406   fK0sAssocPtLtArm->GetYaxis()->SetTitle("L_{T} (cm)");
 407   fK0sAssocPtLtArm->GetZaxis()->SetTitle("centrality");
 408   fOutput->Add(fK0sAssocPtLtArm);
 409   */
 410
 411   fK0sAssocPtLtArm =
 412     new TH3F("fK0sAssocPtLtArm","K^{0}_{S} Assoc: p_{T} vs y vs centrality",nbins,pMin,pMax,30,-1.5,1.5,100,0.,100.);
 413   fK0sAssocPtLtArm->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 414   fK0sAssocPtLtArm->GetYaxis()->SetTitle("y");
 415   fK0sAssocPtLtArm->GetZaxis()->SetTitle("centrality");
 416   fOutput->Add(fK0sAssocPtLtArm);
 417
 418   fK0sAssocPtRap    = new TH2F("fK0sAssocPtRap","K^{0}_{S} Assoc",nbins,pMin,pMax,30,-1.5,1.5);
 419   fK0sAssocPtRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 420   fK0sAssocPtRap->GetYaxis()->SetTitle("y");
 421   fOutput->Add(fK0sAssocPtRap);
 422
 423   fK0sAssocPtEta    = new TH2F("fK0sAssocPtEta","K^{0}_{S} Assoc",nbins,pMin,pMax,30,-1.5,1.5);
 424   fK0sAssocPtEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 425   fK0sAssocPtEta->GetYaxis()->SetTitle("#eta");
 426   fOutput->Add(fK0sAssocPtEta);
 427
 428   // Lambda MC-true:
 429   fLambdaMCPt = new TH1F("fLambdaMCPt","#Lambda MC",nbins,pMin,pMax);
 430   fLambdaMCPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 431   fOutput->Add(fLambdaMCPt);
 432
 433   fLambdaMCPtRap = new TH2F("fLambdaMCPtRap","#Lambda MC",nbins,pMin,pMax,30,-1.5,1.5);
 434   fLambdaMCPtRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 435   fLambdaMCPtRap->GetYaxis()->SetTitle("y");
 436   fOutput->Add(fLambdaMCPtRap);
 437
 438   fLambdaMCPtEta = new TH2F("fLambdaMCPtEta","#Lambda MC",nbins,pMin,pMax,30,-1.5,1.5);
 439   fLambdaMCPtEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 440   fLambdaMCPtEta->GetYaxis()->SetTitle("#eta");
 441   fOutput->Add(fLambdaMCPtEta);
 442
 443   fLambdaMCPtLt =
 444     new TH3F("fLambdaMCPtLt", "#Lambda MC",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 445   fLambdaMCPtLt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 446   fLambdaMCPtLt->GetYaxis()->SetTitle("L_{T} (cm)");
 447   fLambdaMCPtLt->GetZaxis()->SetTitle("centrality");
 448   fOutput->Add(fLambdaMCPtLt);
 449
 450   // Lambda MC-Association:
 451   fLambdaAssocPt =
 452     new TH1F("fLambdaAssocPt","#Lambda Assoc: L_{T} vs p_{T}",nbins,pMin,pMax);
 453   fLambdaAssocPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 454   fOutput->Add(fLambdaAssocPt);
 455
 456   fLambdaAssocPtLt =
 457     new TH3F("fLambdaAssocPtLt","#Lambda Assoc: L_{T} vs p_{T}",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 458   fLambdaAssocPtLt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 459   fLambdaAssocPtLt->GetYaxis()->SetTitle("L_{T} (cm)");
 460   fLambdaAssocPtLt->GetZaxis()->SetTitle("centrality");
 461   fOutput->Add(fLambdaAssocPtLt);
 462
 463   fLambdaAssocPtLtArm =
 464     new TH3F("fLambdaAssocPtLtArm","#Lambda Assoc: L_{T} vs p_{T}",nbins,pMin,pMax,2*nbins,lMin,lMax,100,0.,100.);
 465   fLambdaAssocPtLtArm->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 466   fLambdaAssocPtLtArm->GetYaxis()->SetTitle("L_{T} (cm)");
 467   fLambdaAssocPtLtArm->GetZaxis()->SetTitle("centrality");
 468   fOutput->Add(fLambdaAssocPtLtArm);
 469
 470   fLambdaAssocPtRap = new TH2F("fLambdaAssocPtRap", "#Lambda Assoc",nbins,pMin,pMax,30,-1.5,1.5);
 471   fLambdaAssocPtRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 472   fLambdaAssocPtRap->GetYaxis()->SetTitle("y");
 473   fOutput->Add(fLambdaAssocPtRap);
 474
 475   fLambdaAssocPtEta = new TH2F("fLambdaAssocPtEta", "#Lambda Assoc",nbins,pMin,pMax,30,-1.5,1.5);
 476   fLambdaAssocPtEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 477   fLambdaAssocPtEta->GetYaxis()->SetTitle("#eta");
 478   fOutput->Add(fLambdaAssocPtEta);
 479
 480   }
 481
 482   // ====== Reconstruction information in AOD ====== //
 483   fHistArmenterosPodolanski  =
 484     new TH3F("fHistArmenterosPodolanski","Armenteros-Podolanski phase space;#alpha;p_{t} arm",
 485              100,-1.0,1.0,50,0,0.5,4,-0.5,3.5);
 486   fHistArmenterosPodolanski->GetZaxis()->SetBinLabel(1,"All events");
 487   fHistArmenterosPodolanski->GetZaxis()->SetBinLabel(3,"Triggered events");
 488   fOutput->Add(fHistArmenterosPodolanski);
 489
 490   fHistArmPodBckg =
 491     new TH3F("fHistArmPodBckg","Background: Armenteros-Podolanski phase space;#alpha;p_{t} arm",
 492              100,-1.0,1.0,50,0,0.5,4,-0.5,3.5);
 493   fHistArmPodBckg->GetZaxis()->SetBinLabel(1,"K0s: All events");
 494   fHistArmPodBckg->GetZaxis()->SetBinLabel(2,"K0s: Trig events");
 495   fHistArmPodBckg->GetZaxis()->SetBinLabel(3,"Lambda: All events");
 496   fHistArmPodBckg->GetZaxis()->SetBinLabel(4,"Lambda: Trig events");
 497   fOutput->Add(fHistArmPodBckg);
 498
 499   // ****** K0s ******
 500   fK0sMass =
 501     new TH3F("fK0sMass", "K^{0}_{s}: mass vs p_{T}",nbins/2,0.448,0.548,nbins,pMin,pMax,6,-0.5,5.5);
 502   fK0sMass->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 503   fK0sMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 504   fK0sMass->GetZaxis()->SetBinLabel(1,"All events");
 505   fK0sMass->GetZaxis()->SetBinLabel(2,"Central events");
 506   fK0sMass->GetZaxis()->SetBinLabel(3,"Arm-Pod Cut");
 507   fK0sMass->GetZaxis()->SetBinLabel(4,"Triggered events");
 508   fK0sMass->GetZaxis()->SetBinLabel(5,"Triggered central events");
 509   fOutput->Add(fK0sMass);
 510
 511   fK0sPtLtSB =
 512     new TH2F("fK0sPtLtSB","K^{0}_{s}: L_{T} vs p_{T}, side-band subtracted",
 513              nbins,pMin,pMax,2*nbins,lMin,lMax);
 514   fK0sPtLtSB->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 515   fK0sPtLtSB->GetYaxis()->SetTitle("L_{T} (cm)");
 516   fOutput->Add(fK0sPtLtSB);
 517
 518   fK0sPtvsEta =
 519     new TH3F("fK0sPtvsEta","K^{0}_{s}: p_{T} vs #eta",nbins,pMin,pMax,30,-1.5,1.5,4,-0.5,3.5);
 520   fK0sPtvsEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 521   fK0sPtvsEta->GetYaxis()->SetTitle("#eta");
 522   fK0sPtvsEta->GetZaxis()->SetBinLabel(1,"All events");
 523   fK0sPtvsEta->GetZaxis()->SetBinLabel(2,"Triggered events");
 524   fK0sPtvsEta->GetZaxis()->SetBinLabel(3,"All ev wi Inv Mass cut");
 525   fK0sPtvsEta->GetZaxis()->SetBinLabel(4,"Trig ev wi Inv Mass cut");
 526   fOutput->Add(fK0sPtvsEta);
 527
 528   fK0sPtvsRap =
 529     new TH3F("fK0sPtvsRap","K^{0}_{s}: p_{T} vs y",nbins,pMin,pMax,30,-1.5,1.5,4,-0.5,3.5);
 530   fK0sPtvsRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 531   fK0sPtvsRap->GetYaxis()->SetTitle("y");
 532   fK0sPtvsRap->GetZaxis()->SetBinLabel(1,"All events");
 533   fK0sPtvsRap->GetZaxis()->SetBinLabel(2,"Triggered events");
 534   fK0sPtvsRap->GetZaxis()->SetBinLabel(3,"All ev wi Inv Mass cut");
 535   fK0sPtvsRap->GetZaxis()->SetBinLabel(4,"Trig ev wi Inv Mass cut");
 536   fOutput->Add(fK0sPtvsRap);
 537
 538   fK0sEtaPhi =
 539     new TH2F("fK0sEtaPhi","K^{0}_{s}: #phi vs #eta",
 540              nbinsPhi,0.,2.*TMath::Pi(),100,-1.,1.);
 541   fK0sEtaPhi->GetXaxis()->SetTitle("#phi (rad)");
 542   fK0sEtaPhi->GetYaxis()->SetTitle("#eta");
 543   fOutput->Add(fK0sEtaPhi);
 544
 545
 546   fK0sMassPtPhi  =
 547     new TH3F("fK0sMassPtPhi","K^{0}_{s}: mass vs pt vs #phi",
 548              nbins/2,0.448,0.548,nbins,pMin,pMax,nbinsPhi,0.,2.*TMath::Pi());
 549   fK0sMassPtPhi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 550   fK0sMassPtPhi->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 551   fK0sMassPtPhi->GetZaxis()->SetTitle("#phi (rad)");
 552   fOutput->Add(fK0sMassPtPhi);
 553
 554
 555   // Correlations:
 556   fK0sMassPtvsPtL =
 557     new TH3F("fK0sMassPtvsPtL", "K^{0}_{s}: mass vs p_{T} vs p_{T,l}",
 558              nbins/2,0.448,0.548,
 559              nbins,pMin,pMax,
 560              nbinPtLP,pMin,ptMaxLP);
 561   fK0sMassPtvsPtL->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 562   fK0sMassPtvsPtL->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 563   fK0sMassPtvsPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 564   fOutput->Add(fK0sMassPtvsPtL);
 565
 566   /*
 567   fK0sSiPtL =
 568     new TH3F("fK0sSiPtL","K^{0}_{s}: L_{T} vs p_{T} vs p_{t,L}, side-band subtracted",
 569              nbins,pMin,pMax,2*nbins,lMin,lMax,
 570              nbinPtLP,pMin,ptMaxLP);
 571   fK0sSiPtL->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 572   fK0sSiPtL->GetYaxis()->SetTitle("L_{T} (cm)");
 573   fK0sSiPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 574   fOutput->Add(fK0sSiPtL);
 575
 576   fK0sDaughtersPt =
 577     new TH2F("fK0sDaughtersPt","K^{0}_{s}: p_{T} of daughters",
 578              nbins/2,pMin,pMax/2.,nbins/2,pMin,pMax/2.);
 579   fK0sDaughtersPt->GetXaxis()->SetTitle("Neg. Daught p_{T} (GeV/c)");
 580   fK0sDaughtersPt->GetYaxis()->SetTitle("Pos. Daught p_{T} (GeV/c)");
 581   fOutput->Add(fK0sDaughtersPt);
 582   */
 583
 584   fK0sdPhiPtAssocPtL =
 585     new TH3F("fK0sdPhiPtAssocPtL","K^{0}_{s}: #Delta#phi vs  p_{T,a} vs p_{T,l}",
 586              nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),nbins,pMin,pMax,nbinPtLP,pMin,ptMaxLP);
 587   fK0sdPhiPtAssocPtL->GetXaxis()->SetTitle("#Delta#phi (rad)");
 588   fK0sdPhiPtAssocPtL->GetYaxis()->SetTitle("p_{T,a} (GeV/c)");
 589   fK0sdPhiPtAssocPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 590   fOutput->Add(fK0sdPhiPtAssocPtL);
 591
 592
 593   fK0sDCADaugToPrimVtx
 594     = new TH3F("fK0sDCADaugToPrimVtx","K^{0}_{S} Bckg: dca daughter vs. p_{T,l}",
 595                90,0.,3.3,90,0.,3.3,nbinPtLP,pMin,ptMaxLP);
 596   fK0sDCADaugToPrimVtx->GetXaxis()->SetTitle("DCA Pos daug (cm)");
 597   fK0sDCADaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg daug (cm)");
 598   fK0sDCADaugToPrimVtx->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 599   fOutput->Add(fK0sDCADaugToPrimVtx);
 600
 601   //Double_t kVtxBins[] = {-10.,-7.,-4.,-2.,0.,2.,4.,7.,10.};
 602   char hNameHist[256];
 603   for(Int_t k=0;k<kN1;k++){
 604
 605     // Monte-Carlo level:
 606     sprintf(hNameHist, "fK0sdPhidEtaMC_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 607     fK0sdPhidEtaMC[k] = new TH3F(hNameHist,"K^{0}_{S} MC: #Delta#phi vs #Delta#eta vs p_{T,l}",
 608                                   nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 609                                   nbinsdEta,-2.,2.,
 610                                   nbinsVtx,-10.,10.);
 611     fK0sdPhidEtaMC[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 612     fK0sdPhidEtaMC[k]->GetYaxis()->SetTitle("#Delta#eta");
 613     fK0sdPhidEtaMC[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 614     fOutput->Add(fK0sdPhidEtaMC[k]);
 615
 616     sprintf(hNameHist, "fK0sdPhidEtaMC_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 617     fK0sdPhidEtaMCCent[k] = new TH3F(hNameHist,"K^{0}_{S} MC: #Delta#phi vs #Delta#eta vs p_{T,l}",
 618                                       nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 619                                       nbinsdEta,-2.,2.,
 620                                       nbinsVtx,-10.,10.);
 621     fK0sdPhidEtaMCCent[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 622     fK0sdPhidEtaMCCent[k]->GetYaxis()->SetTitle("#Delta#eta");
 623     fK0sdPhidEtaMCCent[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 624     fOutput->Add(fK0sdPhidEtaMCCent[k]);
 625
 626     // Reconstruction level:
 627     sprintf(hNameHist, "fK0sdPhidEtaPtL_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 628     fK0sdPhidEtaPtL[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 629                                   nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 630                                   nbinsdEta,-2.,2.,
 631                                   nbinsVtx,-10.,10.);
 632     fK0sdPhidEtaPtL[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 633     fK0sdPhidEtaPtL[k]->GetYaxis()->SetTitle("#Delta#eta");
 634     fK0sdPhidEtaPtL[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 635     fOutput->Add(fK0sdPhidEtaPtL[k]);
 636
 637     sprintf(hNameHist, "fK0sdPhidEtaPtL_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 638     fK0sdPhidEtaPtLCent[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 639                                       nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 640                                       nbinsdEta,-2.,2.,
 641                                       nbinsVtx,-10.,10.);
 642     fK0sdPhidEtaPtLCent[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 643     fK0sdPhidEtaPtLCent[k]->GetYaxis()->SetTitle("#Delta#eta");
 644     fK0sdPhidEtaPtLCent[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 645     fOutput->Add(fK0sdPhidEtaPtLCent[k]);
 646
 647     sprintf(hNameHist, "fK0sdPhidEtaPtL_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 648     fK0sdPhidEtaPtLBckg[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 649                                       nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 650                                       nbinsdEta,-2.,2.,
 651                                       nbinsVtx,-10.,10.);
 652     fK0sdPhidEtaPtLBckg[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 653     fK0sdPhidEtaPtLBckg[k]->GetYaxis()->SetTitle("#Delta#eta");
 654     fK0sdPhidEtaPtLBckg[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 655     fOutput->Add(fK0sdPhidEtaPtLBckg[k]);
 656
 657     sprintf(hNameHist, "fK0sdPhidEtaPtL_Cent_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 658     fK0sdPhidEtaPtLCentBckg[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 659                                           nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 660                                           nbinsdEta,-2.,2.,
 661                                           nbinsVtx,-10.,10.);
 662     fK0sdPhidEtaPtLCentBckg[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 663     fK0sdPhidEtaPtLCentBckg[k]->GetYaxis()->SetTitle("#Delta#eta");
 664     fK0sdPhidEtaPtLCentBckg[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 665     fOutput->Add(fK0sdPhidEtaPtLCentBckg[k]);
 666
 667     // To be done with Armenteros Podolanski cut
 668     sprintf(hNameHist, "fK0sdPhidEtaPtL2_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 669     fK0sdPhidEtaPtL2[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 670                                    nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 671                                    nbinsdEta,-2.,2.,
 672                                    nbinsVtx,-10.,10.);
 673     fK0sdPhidEtaPtL2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 674     fK0sdPhidEtaPtL2[k]->GetYaxis()->SetTitle("#Delta#eta");
 675     fK0sdPhidEtaPtL2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 676     fOutput->Add(fK0sdPhidEtaPtL2[k]);
 677
 678     sprintf(hNameHist, "fK0sdPhidEtaPtL2_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 679     fK0sdPhidEtaPtLCent2[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 680                                        nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 681                                        nbinsdEta,-2.,2.,
 682                                        nbinsVtx,-10.,10.);
 683     fK0sdPhidEtaPtLCent2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 684     fK0sdPhidEtaPtLCent2[k]->GetYaxis()->SetTitle("#Delta#eta");
 685     fK0sdPhidEtaPtLCent2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 686     fOutput->Add(fK0sdPhidEtaPtLCent2[k]);
 687
 688     sprintf(hNameHist, "fK0sdPhidEtaPtL2_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 689     fK0sdPhidEtaPtLBckg2[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 690                                        nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 691                                        nbinsdEta,-2.,2.,
 692                                        nbinsVtx,-10.,10.);
 693     fK0sdPhidEtaPtLBckg2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 694     fK0sdPhidEtaPtLBckg2[k]->GetYaxis()->SetTitle("#Delta#eta");
 695     fK0sdPhidEtaPtLBckg2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 696     fOutput->Add(fK0sdPhidEtaPtLBckg2[k]);
 697
 698     sprintf(hNameHist, "fK0sdPhidEtaPtL2_Cent_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 699     fK0sdPhidEtaPtLCentBckg2[k] = new TH3F(hNameHist,"K^{0}_{S}: #Delta#phi vs #Delta#eta vs p_{T,l}",
 700                                            nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 701                                            nbinsdEta,-2.,2.,
 702                                            nbinsVtx,-10.,10.);
 703     fK0sdPhidEtaPtLCentBckg2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 704     fK0sdPhidEtaPtLCentBckg2[k]->GetYaxis()->SetTitle("#Delta#eta");
 705     fK0sdPhidEtaPtLCentBckg2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 706     fOutput->Add(fK0sdPhidEtaPtLCentBckg2[k]);
 707   }
 708
 709   // Correlations (side-band):
 710   fK0sBckgDecLength
 711     = new TH2F("fK0sBckgDecLength","K^{0}_{S} Bckg: c#tau vs. p_{T,l}",
 712                100,0.,15.,nbinPtLP,pMin,ptMaxLP);
 713   fK0sBckgDecLength->GetXaxis()->SetTitle("c#tau (cm)");
 714   fK0sBckgDecLength->GetYaxis()->SetTitle("p_{T,l} (GeV/c)");
 715   fOutput->Add(fK0sBckgDecLength);
 716
 717   fK0sBckgDCADaugToPrimVtx
 718     = new TH3F("fK0sBckgDCADaugToPrimVtx","K^{0}_{S} Bckg: dca daughter vs. p_{T,l}",
 719                90,0.,3.3,90,0.,3.3,nbinPtLP,pMin,ptMaxLP);
 720   fK0sBckgDCADaugToPrimVtx->GetXaxis()->SetTitle("DCA Pos daug (cm)");
 721   fK0sBckgDCADaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg daug (cm)");
 722   fK0sBckgDCADaugToPrimVtx->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 723   fOutput->Add(fK0sBckgDCADaugToPrimVtx);
 724
 725   fK0sBckgEtaPhi =
 726     new TH2F("fK0sBckgEtaPhi","K^{0}_{s} Bckg: #phi vs #eta",
 727              nbinsPhi,0.,2.*TMath::Pi(),100,-1.,1.);
 728   fK0sBckgEtaPhi->GetXaxis()->SetTitle("#phi (rad)");
 729   fK0sBckgEtaPhi->GetYaxis()->SetTitle("#eta");
 730   fOutput->Add(fK0sBckgEtaPhi);
 731
 732   fK0sBckgPhiRadio
 733     = new TH2F("fK0sBckgPhiRadio","K^{0}_{S} Bckg: #phi vs l_{T}",
 734                nbinsPhi,0.,2.*TMath::Pi(),2*nbins,lMin,lMax);
 735   fK0sBckgPhiRadio->GetXaxis()->SetTitle("#phi (rad)");
 736   fK0sBckgPhiRadio->GetYaxis()->SetTitle("l_{T} (cm)");
 737   fOutput->Add(fK0sBckgPhiRadio);
 738
 739   fK0sBckgDCANegDaugToPrimVtx
 740     = new TH2F("fK0sBckgDCANegDaugToPrimVtx","K^{0}_{S} Bckg: dca NegDaughter",
 741                6,-0.5,5.5,90,0.,3.3);
 742   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetTitle("MC Production");
 743   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(1,"Rec");
 744   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(2,"Primary");
 745   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(3,"Weak Decay");
 746   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(4,"Gamma conv.");
 747   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(5,"Unidentified mother");
 748   fK0sBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(6,"Other");
 749   fK0sBckgDCANegDaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg Daug (cm)");
 750   fOutput->Add(fK0sBckgDCANegDaugToPrimVtx);
 751
 752   fK0sBckgDCAPosDaugToPrimVtx
 753     = new TH2F("fK0sBckgDCAPosDaugToPrimVtx","K^{0}_{S} Bckg: dca PosDaughter",
 754                6,-0.5,5.5,90,0.,3.3);
 755   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetTitle("MC Production");
 756   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(1,"Rec");
 757   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(2,"Primary");
 758   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(3,"Weak Decay");
 759   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(4,"Gamma conv.");
 760   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(5,"Unidentified mother");
 761   fK0sBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(6,"Other");
 762   fK0sBckgDCAPosDaugToPrimVtx->GetYaxis()->SetTitle("DCA Pos Daug (cm)");
 763   fOutput->Add(fK0sBckgDCAPosDaugToPrimVtx);
 764
 765   fK0sMassCascade
 766     = new TH2F("fK0sMassCascade","Cascade Reconstruction in K^{0}_{S} Bckg",650, 1.2, 2.5,4,0.5,4.5);
 767   fK0sMassCascade->GetXaxis()->SetTitle("Invariant Mass (GeV/c^{2})");
 768   fK0sMassCascade->GetYaxis()->SetTitle("Cascade type");
 769   fOutput->Add(fK0sMassCascade);
 770
 771
 772   // ****** Lambda ******
 773   fLambdaMass =
 774     new TH3F("fLambdaMass","Mass vs p_{T} for \\Lambda",nbins,1.065,1.165,nbins,pMin,pMax,6,-0.5,5.5);
 775   fLambdaMass->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 776   fLambdaMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 777   fLambdaMass->GetZaxis()->SetBinLabel(1,"All events");
 778   fLambdaMass->GetZaxis()->SetBinLabel(2,"Central events");
 779   fLambdaMass->GetZaxis()->SetBinLabel(3,"Arm-Pod cut");
 780   fLambdaMass->GetZaxis()->SetBinLabel(4,"Triggered events");
 781   fLambdaMass->GetZaxis()->SetBinLabel(5,"Triggered central events");
 782   fOutput->Add(fLambdaMass);
 783
 784   fLambdaPtLtSB =
 785     new TH2F("fLambdaPtLtSB","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
 786              nbins,pMin,pMax,nbins,2*lMin,lMax);
 787   fLambdaPtLtSB->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 788   fLambdaPtLtSB->GetYaxis()->SetTitle("L_{T} (cm)");
 789   fOutput->Add(fLambdaPtLtSB);
 790
 791   fLambdaPtvsEta =
 792     new TH3F("fLambdaPtvsEta","\\Lambda: p_{T} vs #eta",nbins,pMin,pMax,30,-1.5,1.5,4,-0.5,3.5);
 793   fLambdaPtvsEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 794   fLambdaPtvsEta->GetYaxis()->SetTitle("#eta");
 795   fK0sPtvsEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 796   fLambdaPtvsEta->GetYaxis()->SetTitle("#eta");
 797   fLambdaPtvsEta->GetZaxis()->SetBinLabel(1,"All events");
 798   fLambdaPtvsEta->GetZaxis()->SetBinLabel(2,"Triggered events");
 799   fLambdaPtvsEta->GetZaxis()->SetBinLabel(3,"All ev wi Inv Mass cut");
 800   fLambdaPtvsEta->GetZaxis()->SetBinLabel(4,"Trig ev wi Inv Mass cut");
 801   fOutput->Add(fLambdaPtvsEta);
 802
 803   fLambdaPtvsRap =
 804     new TH3F("fLambdaPtvsRap","\\Lambda: p_{T} vs y",nbins,pMin,pMax,30,-1.5,1.5,4,-0.5,3.5);
 805   fLambdaPtvsRap->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 806   fLambdaPtvsRap->GetYaxis()->SetTitle("y");
 807   fLambdaPtvsRap->GetYaxis()->SetTitle("#eta");
 808   fLambdaPtvsRap->GetZaxis()->SetBinLabel(1,"All events");
 809   fLambdaPtvsRap->GetZaxis()->SetBinLabel(2,"Triggered events");
 810   fLambdaPtvsRap->GetZaxis()->SetBinLabel(3,"All ev wi Inv Mass cut");
 811   fLambdaPtvsRap->GetZaxis()->SetBinLabel(4,"Trig ev wi Inv Mass cut");
 812   fOutput->Add(fLambdaPtvsRap);
 813
 814   fLambdaEtaPhi =
 815     new TH2F("fLambdaEtaPhi","#phi vs #eta for #Lambda",
 816              nbinsPhi,0.,2.*TMath::Pi(),100,-1.,1.);
 817   fLambdaEtaPhi->GetXaxis()->SetTitle("#phi (rad)");
 818   fLambdaEtaPhi->GetYaxis()->SetTitle("#eta");
 819   fOutput->Add(fLambdaEtaPhi);
 820
 821
 822  fLambdaMassPtPhi  =
 823     new TH3F("fLambdaMassPtPhi","#Lambda: mass vs pt vs #phi",
 824              nbins/2,0.448,0.548,nbins,pMin,pMax,nbinsPhi,0.,2.*TMath::Pi());
 825   fLambdaMassPtPhi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 826   fLambdaMassPtPhi->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 827   fLambdaMassPtPhi->GetZaxis()->SetTitle("#phi (rad)");
 828   fOutput->Add(fLambdaMassPtPhi);
 829
 830   fLambdadEdx = new TH2F("fLambdadEdx","dE/dx for protons of #Lambda",50,0.2,3,50,0.,6.);
 831   fOutput->Add(fLambdadEdx);
 832
 833   //----------------------
 834   fCPA = new TH1F("fCPA","Cosine of the pointing angle",30,0.9978,1.);
 835   fOutput->Add(fCPA);
 836
 837   fDCA = new TH1F("fDCA","DCA between the daughters",50,0.,1.1);
 838   fOutput->Add(fDCA);
 839   //----------------------
 840
 841   // Correlations:
 842   fLambdaMassPtvsPtL =
 843     new TH3F("fLambdaMassPtvsPtL", "#Lambda: mass vs p_{T} vs p_{T,l}",
 844              nbins,1.065,1.165,
 845              nbins,pMin,pMax,
 846              nbinPtLP,pMin,ptMaxLP);
 847   fLambdaMassPtvsPtL->GetXaxis()->SetTitle("Mass (GeV/c^2)");
 848   fLambdaMassPtvsPtL->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 849   fLambdaMassPtvsPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 850   fOutput->Add(fLambdaMassPtvsPtL);
 851
 852   /*
 853   fLambdaSiPtL =
 854     new TH3F("fLambdaSiPtL","#Lambda: L_{T} vs p_{T} vs p_{t,L}, side-band subtracted",
 855              nbins,pMin,pMax,2*nbins,lMin,lMax,
 856              nbinPtLP,pMin,ptMaxLP);
 857   fLambdaSiPtL->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 858   fLambdaSiPtL->GetYaxis()->SetTitle("L_{T} (cm)");
 859   fLambdaSiPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 860   fOutput->Add(fLambdaSiPtL);
 861
 862   fLambdaDaughtersPt =
 863     new TH2F("fLambdaDaughtersPt","#Lambda: p_{T} of daughters",
 864              nbins/2,pMin,pMax/2.,nbins/2,pMin,pMax/2.);
 865   fLambdaDaughtersPt->GetXaxis()->SetTitle("Neg. Daught p_{T} (GeV/c)");
 866   fLambdaDaughtersPt->GetYaxis()->SetTitle("Pos. Daught p_{T} (GeV/c)");
 867   fOutput->Add(fLambdaDaughtersPt);
 868   */
 869
 870   fLambdadPhiPtAssocPtL =
 871     new TH3F("fLambdadPhiPtAssocPtL","#Lambda: #Delta#phi vs  p_{T,a} vs p_{T,l}",
 872              nbinsPhi,-TMath::PiOver2(),3*TMath::PiOver2(),nbins,pMin,pMax,nbinPtLP,pMin,ptMaxLP);
 873   fLambdadPhiPtAssocPtL->GetXaxis()->SetTitle("#Delta#phi (rad)");
 874   fLambdadPhiPtAssocPtL->GetYaxis()->SetTitle("p_{T,a} (GeV/c)");
 875   fLambdadPhiPtAssocPtL->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 876   fOutput->Add(fLambdadPhiPtAssocPtL);
 877
 878   fLambdaDCADaugToPrimVtx
 879     = new TH3F("fLambdaDCADaugToPrimVtx","#Lambda Bckg: dca daughter vs. p_{T,l}",
 880                90,0.,3.3,90,0.,3.3,nbinPtLP,pMin,ptMaxLP);
 881   fLambdaDCADaugToPrimVtx->GetXaxis()->SetTitle("DCA Pos daug (cm)");
 882   fLambdaDCADaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg daug (cm)");
 883   fLambdaDCADaugToPrimVtx->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
 884   fOutput->Add(fLambdaDCADaugToPrimVtx);
 885
 886
 887   for(Int_t k=0;k<kN1;k++){
 888
 889     // Monte-Carlo level:
 890     sprintf(hNameHist, "fLambdadPhidEtaMC_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 891     fLambdadPhidEtaMC[k] = new TH3F(hNameHist,"#Lambda MC: #Delta#phi vs #Delta#eta vs p_{T,l}",
 892                                   nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 893                                   nbinsdEta,-2.,2.,
 894                                   nbinsVtx,-10.,10.);
 895     fLambdadPhidEtaMC[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 896     fLambdadPhidEtaMC[k]->GetYaxis()->SetTitle("#Delta#eta");
 897     fLambdadPhidEtaMC[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 898     fOutput->Add(fLambdadPhidEtaMC[k]);
 899
 900     sprintf(hNameHist, "fLambdadPhidEtaMC_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 901     fLambdadPhidEtaMCCent[k] = new TH3F(hNameHist,"#Lambda MC: #Delta#phi vs #Delta#eta vs p_{T,l}",
 902                                       nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 903                                       nbinsdEta,-2.,2.,
 904                                       nbinsVtx,-10.,10.);
 905     fLambdadPhidEtaMCCent[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 906     fLambdadPhidEtaMCCent[k]->GetYaxis()->SetTitle("#Delta#eta");
 907     fLambdadPhidEtaMCCent[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 908     fOutput->Add(fLambdadPhidEtaMCCent[k]);
 909
 910     // Reconstruction level:
 911     sprintf(hNameHist, "fLambdadPhidEtaPtL_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 912     fLambdadPhidEtaPtL[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi vs #Delta#eta vs p_{T,l}",
 913                                      nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 914                                      nbinsdEta,-2.,2.,
 915                                      nbinsVtx,-10.,10.);
 916     fLambdadPhidEtaPtL[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 917     fLambdadPhidEtaPtL[k]->GetYaxis()->SetTitle("#Delta#eta");
 918     fLambdadPhidEtaPtL[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 919     fOutput->Add(fLambdadPhidEtaPtL[k]);
 920
 921     sprintf(hNameHist, "fLambdadPhidEtaPtL_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 922     fLambdadPhidEtaPtLCent[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 923                                          nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 924                                          nbinsdEta,-2.,2.,
 925                                          nbinsVtx,-10.,10.);
 926     fLambdadPhidEtaPtLCent[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 927     fLambdadPhidEtaPtLCent[k]->GetYaxis()->SetTitle("#Delta#eta");
 928     fLambdadPhidEtaPtLCent[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 929     fOutput->Add(fLambdadPhidEtaPtLCent[k]);
 930
 931     sprintf(hNameHist, "fLambdadPhidEtaPtL_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 932     fLambdadPhidEtaPtLBckg[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 933                                          nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 934                                          nbinsdEta,-2.,2.,
 935                                          nbinsVtx,-10.,10.);
 936     fLambdadPhidEtaPtLBckg[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 937     fLambdadPhidEtaPtLBckg[k]->GetYaxis()->SetTitle("#Delta#eta");
 938     fLambdadPhidEtaPtLBckg[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 939     fOutput->Add(fLambdadPhidEtaPtLBckg[k]);
 940
 941     sprintf(hNameHist, "fLambdadPhidEtaPtL_Cent_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 942     fLambdadPhidEtaPtLCentBckg[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 943                                              nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 944                                              nbinsdEta,-2.,2.,
 945                                              nbinsVtx,-10.,10.);
 946     fLambdadPhidEtaPtLCentBckg[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 947     fLambdadPhidEtaPtLCentBckg[k]->GetYaxis()->SetTitle("#Delta#eta");
 948     fLambdadPhidEtaPtLCentBckg[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 949     fOutput->Add(fLambdadPhidEtaPtLCentBckg[k]);
 950
 951
 952     //  To be used when the Armenteros Podolanski plot is applied
 953     sprintf(hNameHist, "fLambdadPhidEtaPtL2_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 954     fLambdadPhidEtaPtL2[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi vs #Delta#eta vs p_{T,l}",
 955                                       nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 956                                       nbinsdEta,-2.,2.,
 957                                       nbinsVtx,-10.,10.);
 958     fLambdadPhidEtaPtL2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 959     fLambdadPhidEtaPtL2[k]->GetYaxis()->SetTitle("#Delta#eta");
 960     fLambdadPhidEtaPtL2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 961     fOutput->Add(fLambdadPhidEtaPtL2[k]);
 962
 963     sprintf(hNameHist, "fLambdadPhidEtaPtL2_Cent_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 964     fLambdadPhidEtaPtLCent2[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 965                                           nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 966                                           nbinsdEta,-2.,2.,
 967                                           nbinsVtx,-10.,10.);
 968     fLambdadPhidEtaPtLCent2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 969     fLambdadPhidEtaPtLCent2[k]->GetYaxis()->SetTitle("#Delta#eta");
 970     fLambdadPhidEtaPtLCent2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 971     fOutput->Add(fLambdadPhidEtaPtLCent2[k]);
 972
 973     sprintf(hNameHist, "fLambdadPhidEtaPtL2_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 974     fLambdadPhidEtaPtLBckg2[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 975                                           nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 976                                           nbinsdEta,-2.,2.,
 977                                           nbinsVtx,-10.,10.);
 978     fLambdadPhidEtaPtLBckg2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 979     fLambdadPhidEtaPtLBckg2[k]->GetYaxis()->SetTitle("#Delta#eta");
 980     fLambdadPhidEtaPtLBckg2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 981     fOutput->Add(fLambdadPhidEtaPtLBckg2[k]);
 982
 983     sprintf(hNameHist, "fLambdadPhidEtaPtL2_Cent_Bckg_%.1f_%.1f",kPtBinV0[k],kPtBinV0[k+1]);
 984     fLambdadPhidEtaPtLCentBckg2[k] = new TH3F(hNameHist,"#Lambda: #Delta#phi  vs #Delta#eta vs p_{T,l}",
 985                                               nbinsdPhi,-TMath::PiOver2(),3*TMath::PiOver2(),
 986                                               nbinsdEta,-2.,2.,
 987                                               nbinsVtx,-10.,10.);
 988     fLambdadPhidEtaPtLCentBckg2[k]->GetXaxis()->SetTitle("#Delta#phi (rad)");
 989     fLambdadPhidEtaPtLCentBckg2[k]->GetYaxis()->SetTitle("#Delta#eta");
 990     fLambdadPhidEtaPtLCentBckg2[k]->GetZaxis()->SetTitle("Vertex Z (cm)");
 991     fOutput->Add(fLambdadPhidEtaPtLCentBckg2[k]);
 992
 993   }
 994
 995
 996   // Correlations (side-band):
 997   fLambdaBckgDecLength
 998     = new TH2F("fLambdaBckgDecLength","#Lambda Bckg: c#tau vs. p_{T,l}",
 999                100,0.,25.,nbinPtLP,pMin,ptMaxLP);
1000   fLambdaBckgDecLength->GetXaxis()->SetTitle("c#tau (cm)");
1001   fLambdaBckgDecLength->GetYaxis()->SetTitle("p_{T,l} (GeV/c)");
1002   fOutput->Add(fLambdaBckgDecLength);
1003
1004   fLambdaBckgDCADaugToPrimVtx
1005     = new TH3F("fLambdaBckgDCADaugToPrimVtx","#Lambda Bckg: dca daughter vs. p_{T,l}",
1006                90,0.,3.3,90,0.,3.3,nbinPtLP,pMin,ptMaxLP);
1007   fLambdaBckgDCADaugToPrimVtx->GetXaxis()->SetTitle("DCA Pos daug (cm)");
1008   fLambdaBckgDCADaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg daug (cm)");
1009   fLambdaBckgDCADaugToPrimVtx->GetZaxis()->SetTitle("p_{T,l} (GeV/c)");
1010   fOutput->Add(fLambdaBckgDCADaugToPrimVtx);
1011
1012   fLambdaBckgEtaPhi =
1013     new TH2F("fLambdaBckgEtaPhi","#Lambda Bckg: #phi vs #eta",
1014              nbinsPhi,0.,2.*TMath::Pi(),100,-1.,1.);
1015   fLambdaBckgEtaPhi->GetXaxis()->SetTitle("#phi (rad)");
1016   fLambdaBckgEtaPhi->GetYaxis()->SetTitle("#eta");
1017   fOutput->Add(fLambdaBckgEtaPhi);
1018
1019   fLambdaBckgPhiRadio
1020     = new TH2F("fLambdaBckgPhiRadio","#Lambda Bckg: #phi vs l_{T}",
1021                nbinsPhi,0.,2.*TMath::Pi(),2*nbins,lMin,lMax);
1022   fLambdaBckgPhiRadio->GetXaxis()->SetTitle("#phi (rad)");
1023   fLambdaBckgPhiRadio->GetYaxis()->SetTitle("l_{T} (cm)");
1024   fOutput->Add(fLambdaBckgPhiRadio);
1025
1026
1027   fLambdaBckgDCANegDaugToPrimVtx
1028     = new TH2F("fLambdaBckgDCANegDaugToPrimVtx","#Lambda Bckg: dca NegDaughter",
1029                6,-0.5,5.5,90,0.,3.3);
1030   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetTitle("MC Production");
1031   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(1,"Rec");
1032   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(2,"Primary");
1033   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(3,"Weak Decay");
1034   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(4,"Gamma conv.");
1035   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(5,"Unidentified mother");
1036   fLambdaBckgDCANegDaugToPrimVtx->GetXaxis()->SetBinLabel(6,"Other");
1037   fLambdaBckgDCANegDaugToPrimVtx->GetYaxis()->SetTitle("DCA Neg Daug (cm)");
1038   fOutput->Add(fLambdaBckgDCANegDaugToPrimVtx);
1039
1040
1041   fLambdaBckgDCAPosDaugToPrimVtx
1042     = new TH2F("fLambdaBckgDCAPosDaugToPrimVtx","#Lambda Bckg: dca PosDaughter",
1043                6,-0.5,5.5,90,0.,3.3);
1044   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetTitle("MC Production");
1045   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(1,"Rec");
1046   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(2,"Primary");
1047   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(3,"Weak Decay");
1048   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(4,"Gamma conv.");
1049   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(5,"Unidentified mother");
1050   fLambdaBckgDCAPosDaugToPrimVtx->GetXaxis()->SetBinLabel(6,"Other");
1051   fLambdaBckgDCAPosDaugToPrimVtx->GetYaxis()->SetTitle("DCA Pos Daug (cm)");
1052   fOutput->Add(fLambdaBckgDCAPosDaugToPrimVtx);
1053
1054
1055   fLambdaMassCascade
1056     = new TH2F("fLambdaMassCascade","Cascade Reconstruction in #Lambda Bckg",650, 1.2, 2.5,4,0.5,4.5);
1057   fLambdaMassCascade->GetXaxis()->SetTitle("Invariant Mass (GeV/c^{2})");
1058   fLambdaMassCascade->GetYaxis()->SetTitle("Cascade type");
1059   fOutput->Add(fLambdaMassCascade);
1060
1061
1062   //----------------------
1063
1064   PostData(1, fOutput);
1065
1066 }
1067
1068 //___________________________________________________________________________________________
1069
1070 static Bool_t AcceptTrack(const AliAODTrack *t)
1071 {
1072   // Track criteria for primaries particles
1073   if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
1074
1075   Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
1076   if (nCrossedRowsTPC < 70) return kFALSE;
1077   Int_t findable=t->GetTPCNclsF();
1078   if(findable <= 0) return kFALSE;
1079   if(nCrossedRowsTPC/findable < 0.8) return kFALSE;
1080
1081   return kTRUE;
1082 }
1083
1084
1085 //___________________________________________________________________________________________
1086
1087 static Bool_t AcceptTrackV0(const AliAODTrack *t)
1088 {
1089   // Track criteria for daughter particles of V0 candidate
1090
1091   if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
1092   //if (t->GetKinkIndex(0)>0) return kFALSE;
1093
1094   Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
1095   if (nCrossedRowsTPC < 70) return kFALSE;
1096   Int_t findable=t->GetTPCNclsF();
1097   if(findable <= 0) return kFALSE;
1098   if(nCrossedRowsTPC/findable < 0.8) return kFALSE;
1099
1100   return kTRUE;
1101 }
1102
1103
1104 //___________________________________________________________________________________________
1105
1106 Bool_t AliAnalysisTaskLambdaOverK0sJets::AcceptV0(AliAODVertex *vtx, const AliAODv0 *v1)
1107 {
1108   // Selection for accepting V0 candidates
1109
1110   if (v1->GetOnFlyStatus()) return kFALSE;
1111
1112   //if (v1->Pt() < pMin) return kFALSE; ***
1113
1114   const AliAODTrack *ntrack1=(AliAODTrack *)v1->GetDaughter(1);
1115   const AliAODTrack *ptrack1=(AliAODTrack *)v1->GetDaughter(0);
1116
1117   if( !ntrack1 || !ptrack1 ) return kFALSE;
1118   if( !AcceptTrackV0(ntrack1) ) return kFALSE;
1119   if( !AcceptTrackV0(ptrack1) ) return kFALSE;
1120
1121   if( ntrack1->Charge() == ptrack1->Charge())
1122     return kFALSE;
1123
1124   // Daughters: pseudo-rapidity cut
1125   if ( TMath::Abs(ntrack1->Eta()) > fMaxEtaDaughter  ||
1126        TMath::Abs(ptrack1->Eta()) > fMaxEtaDaughter  )
1127     return kFALSE;
1128
1129   // Daughters: transverse momentum cut
1130   if ( ( ntrack1->Pt() < fMinPtDaughter ) ||
1131        ( ptrack1->Pt() < fMinPtDaughter )  )
1132     return kFALSE;
1133
1134   // Daughters: Impact parameter of daughter to prim vtx
1135   Float_t xy = v1->DcaNegToPrimVertex();
1136   if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE;
1137   xy = v1->DcaPosToPrimVertex();
1138   if (TMath::Abs(xy)<fDCAToPrimVtx) return kFALSE;
1139
1140   // Daughters: DCA
1141   Double_t dca = v1->DcaV0Daughters();
1142   if (dca>fMaxDCADaughter) return kFALSE;
1143
1144   // V0: Cosine of the pointing angle
1145   Double_t cpa=v1->CosPointingAngle(vtx);
1146   if (cpa<fMinCPA) return kFALSE;
1147
1148   // V0: Fiducial volume
1149   Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
1150   Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
1151   if (r2<0.9*0.9) return kFALSE;
1152   if (r2>lMax*lMax) return kFALSE;
1153
1154   return kTRUE;
1155 }
1156
1157 //___________________________________________________________________________________________
1158
1159 static Double_t dPHI(Double_t phi1, Double_t phi2)
1160 {
1161   // Calculate the phi difference between two tracks
1162   Double_t deltaPhi = phi1 - phi2;
1163
1164   if (deltaPhi<-TMath::PiOver2())    deltaPhi = deltaPhi + 2*(TMath::Pi());
1165   if (deltaPhi>(3*TMath::PiOver2()))  deltaPhi = deltaPhi - 2*(TMath::Pi());
1166   return deltaPhi;
1167 }
1168
1169 //___________________________________________________________________________________________
1170
1171 static Double_t MyRapidity(Double_t rE, Double_t rPz)
1172 {
1173   // Local method for rapidity
1174   return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1175 }
1176
1177 //___________________________________________________________________________________________
1178
1179 static Int_t EqualPt(AliAODTrack *trk, const AliAODTrack *nTrk, const AliAODTrack *pTrk)
1180 {
1181   // Local method to compaire the momentum between two tracks
1182
1183   double const kEpsilon = 0.000001;
1184   Int_t    isSamePt = 0;
1185
1186   Double_t p[3];     trk->GetPxPyPz(p);
1187   Double_t pNegTrk[3]; nTrk->GetPxPyPz(pNegTrk);
1188   Double_t pPosTrk[3]; pTrk->GetPxPyPz(pPosTrk);
1189
1190   if( (  fabs(p[0]-pNegTrk[0])<kEpsilon &&
1191          fabs(p[1]-pNegTrk[1])<kEpsilon &&
1192          fabs(p[2]-pNegTrk[2])<kEpsilon )
1193       ||
1194       (  fabs(p[0]-pPosTrk[0])<kEpsilon &&
1195          fabs(p[1]-pPosTrk[1])<kEpsilon &&
1196          fabs(p[2]-pPosTrk[2])<kEpsilon ) )
1197     isSamePt = 1;
1198
1199   return isSamePt;
1200
1201 }
1202
1203 //___________________________________________________________________________________________
1204
1205 void AliAnalysisTaskLambdaOverK0sJets::RecCascade(AliAODTrack *trk1,const AliAODTrack *trk2,const AliAODTrack *trkBch,TString histo)
1206 {
1207   // Local method to reconstruct cascades candidates from the combinations of three tracks
1208   // The input tracks correspond to the trigger particle and the daughter tracks of the V0 candidate (correlation step)
1209   // The trigger particle track will be always consider as a possible daughter of the V0 which coming from the Cascade decay.
1210   // The daughters of the V0 candidates are switched to be the bachelor track for the Cascade reconstruction.
1211
1212   Double_t lMassBach=0., lPtot2Bach=0., lEBach=0.;
1213   Double_t lMassLambda=0., lPtot2Lambda=0., lELambda = 0.;
1214   Double_t pLambda[3] = {0.,0.,0.};
1215   Double_t pCascade[3] = {0.,0.,0.};
1216   Double_t lMassCascade = 0., lPtot2Cascade=0.;
1217
1218   // Two loops are done to consider the posibility to reconstruct a Xi or an Omega
1219   for(Int_t i=0;i<2;i++){
1220
1221     // 0. Check the charge for both tracks: trk1 & trk2.
1222     //    Usefull in the Lambda step.
1223     if( trk1->Charge() == trk2->Charge() )
1224       continue;
1225
1226     // 1. Bachelor: Allocation for the track
1227     if(i==0) // Xi
1228       lMassBach = TDatabasePDG::Instance()->GetParticle(kPiMinus)->Mass();
1229     else if(i==1) //Omega
1230       lMassBach = TDatabasePDG::Instance()->GetParticle(kKMinus)->Mass();
1231
1232     lPtot2Bach = TMath::Power(trkBch->P(),2);
1233
1234     lEBach = TMath::Sqrt(lPtot2Bach + lMassBach*lMassBach);
1235
1236     // 2. Lambda: Kinematical properties
1237     lMassLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1238
1239     pLambda[0] = trk1->Px() + trk2->Px();
1240     pLambda[1] = trk1->Py() + trk2->Py();
1241     pLambda[2] = trk1->Pz() + trk2->Pz();
1242
1243     lPtot2Lambda = pLambda[0]*pLambda[0] +  pLambda[1]*pLambda[1] +  pLambda[2]*pLambda[2];
1244
1245     lELambda = TMath::Sqrt(lPtot2Lambda + lMassLambda*lMassLambda);
1246
1247     // 3. Cascade: Reconstruction
1248     pCascade[0] = pLambda[0] + trkBch->Px();
1249     pCascade[1] = pLambda[1] + trkBch->Py();
1250     pCascade[2] = pLambda[2] + trkBch->Pz();
1251
1252     lPtot2Cascade = pCascade[0]*pCascade[0] + pCascade[1]*pCascade[1] + pCascade[2]*pCascade[2];
1253
1254     lMassCascade = TMath::Sqrt( TMath::Power(lEBach+lELambda,2) - lPtot2Cascade );
1255
1256     // 4. Filling histograms
1257     if( histo.Contains("K0s") ) {
1258       if(i==0)
1259         fK0sMassCascade->Fill(lMassCascade,1);
1260       else if(i==1)
1261         fK0sMassCascade->Fill(lMassCascade,3);
1262     }
1263     else if( histo.Contains("Lambda") ) {
1264       if(i==0)
1265         fLambdaMassCascade->Fill(lMassCascade,1);
1266       else if(i==1)
1267         fLambdaMassCascade->Fill(lMassCascade,3);
1268     }
1269
1270
1271   }
1272
1273 }
1274
1275 //___________________________________________________________________________________________
1276
1277 TArrayD* AliAnalysisTaskLambdaOverK0sJets::V0Loop(AliAODTrack *trkTrig, V0LoopStep_t step, Bool_t isTriggered)
1278 {
1279   // Three options for the 'step' variable:
1280   // 1) TriggerCheck
1281   // 2) Reconstruction and Correlation
1282   // 3) Mixed event
1283
1284   Double_t ptTrig  = -100.;
1285   Double_t phiTrig = -100.;
1286   Double_t etaTrig = -100.;
1287
1288   if( step==kTriggerCheck || isTriggered ){
1289     ptTrig  = trkTrig->Pt();
1290     phiTrig = trkTrig->Phi();
1291     etaTrig = trkTrig->Eta();
1292   }
1293
1294   // ---- 1) TriggerCheck: Variables used to crosscheck if trigger particle is a V0 daughter ---- //
1295   Int_t    isTrigFromV0daug = 0, isV0LP = 0;
1296   Double_t ptV0LP = -1000.;
1297
1298
1299   // *************************************************
1300   // Centrality selection
1301   AliCentrality *cent = fAOD->GetCentrality();
1302   Double_t centrality = cent->GetCentralityPercentile("V0M");
1303
1304   // *************************************************
1305   // MC Event
1306   TClonesArray *stackMC = 0x0;
1307   Double_t mcXv=0., mcYv=0., mcZv=0.;
1308
1309   if(fIsMC){
1310     TList *lst = fAOD->GetList();
1311     stackMC = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
1312     if (!stackMC) {
1313       Printf("ERROR: stack not available");
1314     }
1315
1316     AliAODMCHeader *mcHdr =
1317       (AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
1318
1319     mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();
1320   }
1321
1322   // *************************************************
1323   // V0 loop - AOD
1324   const AliAODVertex *vtx=fAOD->GetPrimaryVertex();
1325   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
1326   Int_t nV0sTot = fAOD->GetNumberOfV0s();
1327   for (Int_t iV0 = 0; iV0 < nV0sTot; iV0++) {
1328
1329     AliAODv0 *v0=fAOD->GetV0(iV0);
1330     if (!v0) continue;
1331     if (!AcceptV0(fAOD->GetPrimaryVertex(),v0)) continue;
1332
1333     const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
1334     const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
1335
1336     // Decay vertex
1337     Double_t xyz[3]; v0->GetSecondaryVtx(xyz);
1338     Double_t dx=xyz[0]-xv, dy=xyz[1]-yv;
1339
1340     // Decay lenght and pt
1341     Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
1342     Double_t pt=TMath::Sqrt(v0->Pt2V0());
1343
1344     Double_t dlK = 0.4977*lt/pt;
1345     Double_t dlL = 1.1157*lt/pt;
1346
1347     // ctau
1348     Bool_t ctK=kTRUE; if (dlK > fMaxCtau*2.68 || dlK < fMinCtau*2.68) ctK=kFALSE;
1349     Bool_t ctL=kTRUE; if (dlL > fMaxCtau*7.89 || dlL < fMinCtau*7.89) ctL=kFALSE;
1350
1351     // Armenteros variables:
1352     Double_t lAlphaV0      =  v0->AlphaV0();
1353     Double_t lPtArmV0      =  v0->PtArmV0();
1354
1355     // MC Association:
1356     Bool_t lComeFromSigma     = kFALSE;
1357     Bool_t lCheckMcK0Short    = kFALSE;
1358     Bool_t lCheckMcLambda     = kFALSE;
1359     Bool_t lCheckMcAntiLambda = kFALSE;
1360
1361     Int_t lMCAssocNegDaug = 0;
1362     Int_t lMCAssocPosDaug = 0;
1363
1364     // ********* MC - Association *********
1365     // In case of injected-MC, the correlations might be done with only natural particles
1366     Bool_t isNaturalPart = kTRUE;
1367     if(step==kCorrelation){ //Correlation
1368
1369       if(fIsMC){
1370         if(!stackMC) goto noas;
1371
1372         Int_t ntrkMC=stackMC->GetEntriesFast();
1373
1374         Int_t nlab = TMath::Abs(ntrack->GetLabel());//** UInt_t
1375         Int_t plab = TMath::Abs(ptrack->GetLabel());
1376
1377         // To avoid futher problems
1378         if ( (nlab<0 || plab<0) ||
1379              (nlab>=ntrkMC || plab>=ntrkMC) )
1380           goto noas;
1381
1382         AliAODMCParticle *nPart=(AliAODMCParticle*)stackMC->UncheckedAt(nlab);
1383         AliAODMCParticle *pPart=(AliAODMCParticle*)stackMC->UncheckedAt(plab);
1384         // MC origin of daughters: Primaries?
1385         if( nPart->IsPhysicalPrimary() ) lMCAssocNegDaug = 1;
1386         if( pPart->IsPhysicalPrimary() ) lMCAssocPosDaug = 1;
1387
1388         if(!nPart || !pPart)   goto noas;
1389
1390         if ( TMath::Abs(nPart->Eta()) > fMaxEtaDaughter ||
1391              TMath::Abs(pPart->Eta()) > fMaxEtaDaughter )
1392           goto noas;
1393
1394         // Daughter momentum cut
1395         if ( ( nPart->Pt() < fMinPtDaughter ) ||
1396              ( pPart->Pt() < fMinPtDaughter )  )
1397           goto noas;
1398
1399         // ----------------------------------------
1400
1401         Int_t lPDGCodeNegDaughter = nPart->GetPdgCode();
1402         Int_t lPDGCodePosDaughter = pPart->GetPdgCode();
1403
1404         Int_t ipMother = pPart->GetMother();
1405         Int_t inMother = nPart->GetMother();
1406
1407         if(inMother<0 || inMother>=ntrkMC) lMCAssocNegDaug = 5;
1408         if(ipMother<0 || ipMother>=ntrkMC) lMCAssocPosDaug = 5;
1409
1410         if(inMother<0 || inMother>=ntrkMC) {  goto noas;}
1411         if(inMother != ipMother) { // did the negative daughter decay ?
1412           AliAODMCParticle *negMotherOfMotherPart = (AliAODMCParticle*)stackMC->UncheckedAt(inMother);
1413           if (negMotherOfMotherPart->GetMother() != ipMother)
1414             goto noas;
1415         }
1416
1417         if (ipMother<0 || ipMother>=ntrkMC)
1418           goto noas;
1419
1420         AliAODMCParticle *p0=(AliAODMCParticle*)stackMC->UncheckedAt(ipMother);
1421         if(!p0)
1422           goto noas;
1423
1424         // ----------------------------------------
1425
1426         if ( (ipMother>=fEndOfHijingEvent) &&
1427              (fEndOfHijingEvent!=-1)     &&
1428              (p0->GetMother()<0) )
1429           isNaturalPart = kFALSE;
1430         else  isNaturalPart = kTRUE;
1431
1432         // ----------------------------------------
1433
1434         if(fSeparateInjPart && !isNaturalPart) goto noas;
1435
1436         Int_t lPDGCodeV0 = p0->GetPdgCode();
1437
1438         // MC origin of daughters:
1439         //Decay from Weak Decay?
1440         if( (TMath::Abs(lPDGCodeV0) == kK0Short) || (TMath::Abs(lPDGCodeV0) == kLambda0) ||
1441             (TMath::Abs(lPDGCodeV0) == kSigmaMinus) || (TMath::Abs(lPDGCodeV0) == kSigmaPlus) ||
1442             (TMath::Abs(lPDGCodeV0) == kSigma0) || (TMath::Abs(lPDGCodeV0) == kXiMinus) ||
1443             (TMath::Abs(lPDGCodeV0) == kOmegaMinus) )
1444           { lMCAssocNegDaug = 2;          lMCAssocPosDaug = 2; }
1445         // Gamma conversion
1446         else if( TMath::Abs(lPDGCodeV0) == kGamma )
1447           { lMCAssocNegDaug = 3;          lMCAssocPosDaug = 3; }
1448         // Unidentied mother:
1449         else
1450           { lMCAssocNegDaug = 4;          lMCAssocPosDaug = 4; }
1451
1452
1453         Int_t lIndexMotherOfMother   = p0->GetMother();
1454         Int_t lPdgcodeMotherOfMother = 0;
1455         if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
1456         else {
1457           AliAODMCParticle *lMCAODMotherOfMother=(AliAODMCParticle*)stackMC->UncheckedAt(lIndexMotherOfMother);
1458           if (!lMCAODMotherOfMother) {lPdgcodeMotherOfMother=0;}
1459           lPdgcodeMotherOfMother = lMCAODMotherOfMother->GetPdgCode();
1460         }
1461
1462         // Daughter momentum cut: ! FIX it in case of AOD ! //MC or REc
1463         if ( (nPart->Pt()  < fMinPtDaughter ) ||
1464              (pPart->Pt()  < fMinPtDaughter ) )
1465           goto noas;
1466
1467         if( (lPDGCodeV0 != kK0Short) &&
1468             (lPDGCodeV0 != kLambda0) &&
1469             (lPDGCodeV0 != kLambda0Bar) )
1470           goto noas;
1471
1472
1473         lComeFromSigma     = kFALSE;
1474         lCheckMcK0Short    = kFALSE;
1475         lCheckMcLambda     = kFALSE;
1476         lCheckMcAntiLambda = kFALSE;
1477
1478         // ----------------------------------------
1479
1480         // K0s
1481         if( (lPDGCodePosDaughter==+211) && (lPDGCodeNegDaughter==-211) &&
1482             (inMother==ipMother) && (lPDGCodeV0==310) ) {
1483
1484           if ( ((AliAODMCParticle*)stackMC->UncheckedAt(ipMother))->IsPrimary()  )
1485             lCheckMcK0Short  = kTRUE;
1486
1487         }
1488         // Lambda
1489         else if( (lPDGCodePosDaughter==+2212) && (lPDGCodeNegDaughter==-211)  &&
1490                  (inMother==ipMother) && (lPDGCodeV0==3122)  ){
1491
1492           if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
1493                ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
1494                ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
1495                ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
1496                ) lComeFromSigma = kTRUE;
1497           else lComeFromSigma = kFALSE;
1498
1499           if ( ((AliAODMCParticle*)stackMC->UncheckedAt(ipMother))->IsPrimary() ||
1500                ( !(((AliAODMCParticle*)stackMC->UncheckedAt(ipMother))->IsPrimary() )
1501                  && (lComeFromSigma==kTRUE) )
1502                ) lCheckMcLambda  = kTRUE;
1503
1504         }
1505         // AntiLambda
1506         else if( (lPDGCodePosDaughter==211) && (lPDGCodeNegDaughter==-2212) &&
1507                  (inMother==ipMother) && (lPDGCodeV0==-3122) ) {
1508
1509
1510           if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
1511                ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
1512                ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
1513                ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
1514                ) lComeFromSigma = kTRUE;
1515           else lComeFromSigma = kFALSE;
1516
1517           if ( ((AliAODMCParticle*)stackMC->UncheckedAt(ipMother))->IsPrimary() ||
1518                ( (!((AliAODMCParticle*)stackMC->UncheckedAt(ipMother))->IsPrimary())
1519                  && (lComeFromSigma==kTRUE) )
1520                ) lCheckMcAntiLambda  = kTRUE;
1521
1522         }
1523
1524         //  ----------------------------------------
1525
1526         if ((p0->Pt())<pMin) goto noas;
1527         if (TMath::Abs(p0->Y())>fYMax ) goto noas;
1528
1529         Double_t dxAs = mcXv - p0->Xv(),  dyAs = mcYv - p0->Yv(),  dzAs = mcZv - p0->Zv();
1530         Double_t l = TMath::Sqrt(dxAs*dxAs + dyAs*dyAs + dzAs*dzAs);
1531
1532         dxAs = mcXv - pPart->Xv(); dyAs = mcYv - pPart->Yv();
1533         //Double_t ltAs = TMath::Sqrt(dxAs*dxAs + dyAs*dyAs);
1534         Double_t ptAs = p0->Pt();
1535         Double_t rapAs = p0->Y();
1536         Double_t etaAs = p0->Eta();
1537
1538         if (l < 0.01) { // Primary V0
1539
1540           if(ctK && lCheckMcK0Short){
1541             fK0sAssocPt->Fill(ptAs);
1542             fK0sAssocPtLt->Fill(ptAs,rapAs,centrality);
1543             fK0sAssocPtRap->Fill(ptAs,rapAs);
1544             fK0sAssocPtEta->Fill(ptAs,etaAs);
1545
1546             if(lPtArmV0 > TMath::Abs(0.2*lAlphaV0) )
1547               fK0sAssocPtLtArm->Fill(ptAs,rapAs,centrality);
1548
1549           }
1550           else if(ctL && lCheckMcLambda) {
1551             fLambdaAssocPt->Fill(ptAs);
1552             fLambdaAssocPtLt->Fill(ptAs,rapAs,centrality);
1553             fLambdaAssocPtRap->Fill(ptAs,rapAs);
1554             fLambdaAssocPtEta->Fill(ptAs,etaAs);
1555
1556             if(lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1557               fLambdaAssocPtLtArm->Fill(ptAs,rapAs,centrality);
1558           }
1559           /*
1560           else if (ctL && lCheckMcAntiLambda){
1561             fAntiLambdaAssocPt->Fill(ptAs);
1562             fAntiLambdaAssocPtLt->Fill(ptAs,ltAs,centrality);
1563             fAntiLambdaAssocPtRap->Fill(ptAs,rapAs);
1564             fAntiLambdaAssocPtEta->Fill(ptAs,etaAs);
1565           }
1566           */
1567
1568         }
1569
1570         // After the kinematical selection of K0s and Lambdas
1571         // it might be that the daugthers are not identified through MC Association
1572         if(lMCAssocNegDaug==0)
1573           lMCAssocNegDaug = 5;
1574         if(lMCAssocPosDaug==0)
1575           lMCAssocPosDaug = 5;
1576
1577
1578       } // End MC-Association
1579
1580       // ************************************
1581
1582     }// End Correlation Step
1583
1584   noas:
1585
1586     Double_t pPos = -100.;
1587     Double_t dedxPos = -1000.;
1588
1589     Double_t nsigPosPion   = 0.;
1590     Double_t nsigPosProton = 0.;
1591     Double_t nsigNegPion   = 0.;
1592     Double_t nsigNegProton = 0.;
1593
1594     if(fUsePID) {
1595       const AliAODPid *pidNeg = ntrack->GetDetPid();
1596       const AliAODPid *pidPos = ptrack->GetDetPid();
1597
1598       if (pidNeg && pidPos) {
1599         pPos = pidPos->GetTPCmomentum();
1600         dedxPos = pidPos->GetTPCsignal()/47.;
1601       }
1602
1603       nsigPosPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(ptrack,AliPID::kPion));
1604       nsigPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton));
1605       nsigNegPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(ntrack,AliPID::kPion));
1606       nsigNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(ntrack,AliPID::kProton));
1607     }
1608
1609     Double_t dcaNeg = -100.;
1610     Double_t dcaPos = -100.;
1611     Double_t lPtNeg = -100.;
1612     Double_t lPtPos = -100.;
1613
1614     Double_t dca   = -100;
1615     Double_t cpa   = -100;
1616     Double_t lEta  = -100.;
1617     Double_t lPhi  = -100.;
1618
1619     Double_t radio = -100.;
1620     Double_t dPhi  = -100.;
1621     Double_t dEta  = -100.;
1622
1623     /*
1624     // Good regions
1625     lPhi  = v0->Phi();
1626     if(lPhi>0. && lPhi<1.8) continue;
1627     if(lPhi>2.2 && lPhi<3.8) continue;
1628     if(lPhi>4.2 && lPhi<5.8) continue;
1629     */
1630
1631     /*
1632     // Bad regions
1633     lPhi  = v0->Phi();
1634     if(lPhi>1.8 && lPhi<2.2) continue;
1635     if(lPhi>3.8 && lPhi<4.2) continue;
1636     if(lPhi>5.8 && lPhi<6.2) continue;
1637     */
1638
1639     if(step==kCorrelation){//Correlation
1640       dcaNeg = v0->DcaNegToPrimVertex();
1641       dcaPos = v0->DcaPosToPrimVertex();
1642       lPtNeg = ntrack->Pt();
1643       lPtPos = ptrack->Pt();
1644
1645       dca   = v0->DcaV0Daughters();
1646       cpa   = v0->CosPointingAngle(fAOD->GetPrimaryVertex());
1647       lEta  = v0->PseudoRapV0();
1648       lPhi  = v0->Phi();
1649
1650       lPhi  = ( (lPhi < 0) ? lPhi + 2*TMath::Pi() : lPhi );
1651       dPhi  = dPHI(phiTrig,lPhi);
1652       dEta  = etaTrig - v0->PseudoRapV0();
1653       radio = TMath::Sqrt(dPhi*dPhi + dEta*dEta);
1654     }
1655
1656     // Comparing the pt of the trigger particle wrt the v0-candidate's daughter:
1657     // It is used as well for the side-band subtraction
1658     Int_t isSameTrk = -1;
1659     isSameTrk = EqualPt(trkTrig,ntrack,ptrack);
1660
1661     // *******************
1662     // Disentangle the V0 candidate
1663
1664     Double_t massK0s = 0., mK0s = 0., sK0s = 0.;
1665     Double_t massLambda = 0., mLambda = 0., sLambda = 0.;
1666     Double_t massAntiLambda = 0.;
1667
1668     Bool_t isCandidate2K0s = kFALSE;
1669     massK0s = v0->MassK0Short();
1670     mK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1671     sK0s = 0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
1672     if (ctK && (TMath::Abs(v0->RapK0Short())<fYMax) &&
1673         (nsigNegPion < fNSigma) && (nsigPosPion < fNSigma) &&
1674         (TMath::Abs(mK0s-massK0s) < 3*sK0s) )
1675       isCandidate2K0s = kTRUE;
1676
1677     Bool_t isCandidate2Lambda = kFALSE;
1678     massLambda = v0->MassLambda();
1679     mLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1680     //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
1681     //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
1682     sLambda=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
1683     if (ctL && (TMath::Abs(v0->RapLambda())<fYMax)  &&
1684         (nsigNegPion < fNSigma)   &&
1685         (nsigPosProton < fNSigma) &&
1686         (TMath::Abs(mLambda-massLambda) < 3*sLambda) )
1687       isCandidate2Lambda = kTRUE;
1688
1689     Bool_t isCandidate2LambdaBar = kFALSE;
1690     massAntiLambda = v0->MassAntiLambda();
1691     if (ctL && (TMath::Abs(v0->RapLambda())<fYMax)  &&
1692         (nsigPosPion < fNSigma)   &&
1693         (nsigNegProton < fNSigma) &&
1694         (TMath::Abs(mLambda-massAntiLambda) < 3*sLambda) )
1695       isCandidate2LambdaBar = kTRUE;
1696
1697
1698     // *******************
1699     //   K0s selection
1700     // *******************
1701     if (ctK && (TMath::Abs(v0->RapK0Short())<fYMax) &&
1702         (nsigNegPion < fNSigma) && (nsigPosPion < fNSigma)
1703         //&& (lPtArmV0 > TMath::Abs(0.2*lAlphaV0) )
1704         && !isCandidate2Lambda && !isCandidate2LambdaBar) {
1705
1706       switch(step) {
1707       case kTriggerCheck:
1708
1709         if (isCandidate2K0s){
1710           if(pt>ptTrig){
1711             isV0LP = 1;
1712             ptV0LP = pt;
1713           }
1714
1715           isTrigFromV0daug = isSameTrk;
1716           if(isTrigFromV0daug)
1717             Printf("  The LP has the same momentum in X and Y as one of the K0s daughters *** iV0 %d",iV0);
1718         }
1719
1720         break; // End K0s selection for TriggerCheck
1721       case kCorrelation:
1722
1723         fK0sMass->Fill(massK0s,pt,0);
1724         fK0sPtvsEta->Fill(pt,lEta,0);
1725         fK0sPtvsRap->Fill(pt,v0->RapK0Short(),0);
1726         if(centrality<10) fK0sMass->Fill(massK0s,pt,1);
1727
1728         if( lPtArmV0 > TMath::Abs(0.2*lAlphaV0) )
1729           fK0sMass->Fill(massK0s,pt,2);
1730
1731         if(pt>2. && pt<5.)
1732           fK0sMassPtPhi->Fill(massK0s,pt,lPhi);
1733
1734         if(isTriggered && isNaturalPart){
1735
1736           fK0sMass->Fill(massK0s,pt,3);
1737           fK0sPtvsEta->Fill(pt,lEta,1);
1738           fK0sPtvsRap->Fill(pt,v0->RapK0Short(),1);
1739
1740           fK0sMassPtvsPtL->Fill(massK0s,pt,ptTrig);
1741           if(centrality<10) fK0sMass->Fill(massK0s,pt,4);
1742
1743         }
1744         // Invariant Mass cut
1745         if (TMath::Abs(mK0s-massK0s) < 3*sK0s) {
1746
1747           fK0sPtLtSB->Fill(pt,lt);
1748           fK0sEtaPhi->Fill(lPhi,lEta);
1749
1750           fK0sPtvsEta->Fill(pt,lEta,2);
1751           fK0sPtvsRap->Fill(pt,v0->RapK0Short(),2);
1752           fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0,0);
1753
1754           if(isTriggered && isNaturalPart){
1755
1756             fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0,2);
1757             fK0sPtvsEta->Fill(pt,lEta,3);
1758             fK0sPtvsRap->Fill(pt,v0->RapK0Short(),3);
1759             //fK0sSiPtL->Fill(pt,lt,ptTrig);
1760             //fK0sDaughtersPt->Fill(lPtNeg,lPtPos);
1761
1762             fK0sdPhiPtAssocPtL->Fill(dPhi,pt,ptTrig);
1763
1764             // *** pt bin and centrality
1765             for(Int_t k=0;k<kN1;k++)
1766               if( (pt>kPtBinV0[k]) && (pt<kPtBinV0[k+1]) ){
1767
1768                 fK0sdPhidEtaPtL[k]->Fill(dPhi,dEta,zv);
1769                 if(lPtArmV0 > TMath::Abs(0.2*lAlphaV0))
1770                   fK0sdPhidEtaPtL2[k]->Fill(dPhi,dEta,zv);
1771
1772                 if(centrality<10){
1773                   fK0sdPhidEtaPtLCent[k]->Fill(dPhi,dEta,zv);
1774                   if(lPtArmV0 > TMath::Abs(0.2*lAlphaV0))
1775                     fK0sdPhidEtaPtLCent2[k]->Fill(dPhi,dEta,zv);
1776                 }
1777
1778                 if(radio < 0.4)
1779                   fK0sDCADaugToPrimVtx->Fill(dcaPos,dcaNeg,ptTrig);
1780
1781               } // End selection pt bin
1782
1783           } // End triggered selection
1784
1785         } // End selection in mass
1786
1787         if( TMath::Abs(mK0s-massK0s + 6.5*sK0s) < 1.5*sK0s ||
1788             TMath::Abs(mK0s-massK0s - 6.5*sK0s) < 1.5*sK0s  ) {
1789
1790           fK0sEtaPhi->Fill(lPhi,lEta,-1);
1791           fK0sPtLtSB->Fill(pt,lt,-1);
1792
1793           //  if(radio < 0.4) // Under the peak
1794           fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,0);
1795
1796           if( !isSameTrk && isTriggered && isNaturalPart){
1797
1798             //fK0sSiPtL->Fill(pt,lt,ptTrig,-1);
1799             fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,1);
1800
1801             // Pt bin
1802             for(Int_t k=0;k<kN1;k++)
1803               if( (pt>kPtBinV0[k]) && (pt<kPtBinV0[k+1]) ){
1804
1805                 fK0sdPhidEtaPtLBckg[k]->Fill(dPhi,dEta,zv);
1806                 if(lPtArmV0 > TMath::Abs(0.2*lAlphaV0))
1807                   fK0sdPhidEtaPtLBckg2[k]->Fill(dPhi,dEta,zv);
1808
1809                 if(centrality<10){
1810                   fK0sdPhidEtaPtLCentBckg[k]->Fill(dPhi,dEta,zv);
1811                   if(lPtArmV0 > TMath::Abs(0.2*lAlphaV0))
1812                     fK0sdPhidEtaPtLCentBckg2[k]->Fill(dPhi,dEta,zv);
1813                 }
1814
1815
1816                 if(radio < 0.4){ // Under the correlation peak
1817                   //fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,1);
1818                   fK0sBckgDecLength->Fill(dlK,ptTrig);
1819                   fK0sBckgDCADaugToPrimVtx->Fill(dcaPos,dcaNeg,ptTrig);
1820                   fK0sBckgEtaPhi->Fill(lPhi,lEta);
1821                   fK0sBckgPhiRadio->Fill(lPhi,lt);
1822
1823                   // MC Association of daughter particles
1824                   for(Int_t ii=0;ii<5;ii++){
1825                     if(lMCAssocNegDaug == ii){
1826                       fK0sBckgDCANegDaugToPrimVtx->Fill(ii,dcaNeg);
1827                       //if eta selection for the trigger particle
1828                     }
1829                     if(lMCAssocPosDaug == ii){
1830                       fK0sBckgDCAPosDaugToPrimVtx->Fill(ii,dcaPos);
1831                       //if eta selection for the trigger particle
1832                     }
1833
1834                   }
1835
1836                   RecCascade(trkTrig,ntrack,ptrack,"K0s");
1837                   RecCascade(trkTrig,ptrack,ntrack,"K0s");
1838
1839                 }// End salection in the correlation peak
1840
1841               } // End selection pt bin
1842
1843           } // End triggered selection
1844
1845         }// End selection in outside the mass cut
1846
1847         break; // End K0s selection for Corrleation
1848       default:
1849         Printf( " Selection of 'step' is not set properly");
1850         break;
1851
1852       }// End switch
1853
1854     } // End K0s selection
1855
1856     // *******************
1857     // Lambda selection
1858     // *******************
1859     if ( ctL && (TMath::Abs(v0->RapLambda())<fYMax)  &&
1860          (nsigNegPion < fNSigma) && (nsigPosProton < fNSigma)
1861          //&& (lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1862          && !isCandidate2K0s && !isCandidate2LambdaBar) {
1863
1864
1865       switch(step) {
1866       case kTriggerCheck:
1867
1868         if (isCandidate2Lambda){
1869           if(pt>ptTrig) {
1870             isV0LP = 1;
1871             ptV0LP = pt;
1872           }
1873
1874           isTrigFromV0daug = isSameTrk;
1875           if(isTrigFromV0daug)
1876             Printf("  The LP has the same momentum in X and Y as one of the Lambda daughters *** iV0 %d",iV0);
1877         }
1878
1879         break; // End Lambda selection for TriggerCheck
1880       case kCorrelation:
1881
1882         if(fUsePID)
1883           fLambdadEdx->Fill(pPos,dedxPos);
1884
1885         fLambdaMass->Fill(massLambda,pt,0);
1886         fLambdaPtvsEta->Fill(pt,lEta,0);
1887         fLambdaPtvsRap->Fill(pt,v0->RapLambda(),0);
1888         if(centrality<10) fLambdaMass->Fill(massLambda,pt,1);
1889
1890         if( lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1891           fLambdaMass->Fill(massLambda,pt,2);
1892
1893
1894         if(pt>2. && pt<5.)
1895           fLambdaMassPtPhi->Fill(massLambda,pt,lPhi);
1896
1897         if( isTriggered && isNaturalPart ){
1898
1899           fLambdaMass->Fill(massLambda,pt,3);
1900           fLambdaPtvsEta->Fill(pt,lEta,1);
1901           fLambdaPtvsRap->Fill(pt,v0->RapLambda(),1);
1902
1903           fLambdaMassPtvsPtL->Fill(massLambda,pt,ptTrig);
1904           if(centrality<10)  fLambdaMass->Fill(massLambda,pt,4);
1905
1906         }
1907         // Invariant Mass cut
1908         if (TMath::Abs(mLambda-massLambda) < 3*sLambda) {
1909
1910           fDCA->Fill(dca,1);
1911           fCPA->Fill(cpa,1);
1912           fLambdaPtLtSB->Fill(pt,lt);
1913           fLambdaEtaPhi->Fill(lPhi,lEta);
1914
1915           fLambdaPtvsEta->Fill(pt,lEta,2);
1916           fLambdaPtvsRap->Fill(pt,v0->RapLambda(),2);
1917           fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0,0);
1918
1919           if( isTriggered && isNaturalPart ){
1920             fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0,2);
1921             fLambdaPtvsEta->Fill(pt,lEta,3);
1922             fLambdaPtvsRap->Fill(pt,v0->RapLambda(),3);
1923             //fLambdaSiPtL->Fill(pt,lt,ptTrig); //***
1924             //fLambdaDaughtersPt->Fill(lPtNeg,lPtPos);
1925
1926             fLambdadPhiPtAssocPtL->Fill(dPhi,pt,ptTrig);
1927
1928             // *** pt bin and centrality
1929             for(Int_t k=0;k<kN1;k++)
1930               if( (pt>kPtBinV0[k]) && (pt<kPtBinV0[k+1]) ){
1931                 fLambdadPhidEtaPtL[k]->Fill(dPhi,dEta,zv);
1932
1933                 if( lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1934                   fLambdadPhidEtaPtL2[k]->Fill(dPhi,dEta,zv);
1935
1936                   if(centrality<10){
1937                     fLambdadPhidEtaPtLCent[k]->Fill(dPhi,dEta,zv);
1938                     if( lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1939                       fLambdadPhidEtaPtLCent2[k]->Fill(dPhi,dEta,zv);
1940                   }
1941
1942                 if(radio < 0.4)
1943                   fLambdaDCADaugToPrimVtx->Fill(dcaPos,dcaNeg,ptTrig);
1944
1945               } // End selection pt bin
1946
1947           } // End triggered selection
1948
1949         } // End selection in mass
1950
1951         if( (TMath::Abs(mLambda-massLambda + 6.5*sLambda) < 1.5*sLambda) ||
1952             (TMath::Abs(mLambda-massLambda - 6.5*sLambda) < 1.5*sLambda) ){
1953
1954           fCPA->Fill(cpa,-1);
1955           fDCA->Fill(dca,-1);
1956           fLambdaPtLtSB->Fill(pt,lt,-1);
1957           fLambdaEtaPhi->Fill(lPhi,lEta,-1);
1958
1959           //   if(radio < 0.4)
1960           fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,2);
1961
1962           if( !isSameTrk && isTriggered && isNaturalPart ){
1963
1964             //fLambdaSiPtL->Fill(pt,lt,ptTrig,-1); //***
1965             fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,3);
1966
1967             // Pt bin
1968             for(Int_t k=0;k<kN1;k++)
1969               if( (pt>kPtBinV0[k]) && (pt<kPtBinV0[k+1]) ){
1970
1971                 fLambdadPhidEtaPtLBckg[k]->Fill(dPhi,dEta,zv);
1972                 if( lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1973                   fLambdadPhidEtaPtLBckg2[k]->Fill(dPhi,dEta,zv);
1974
1975                 if(centrality<10){
1976                   fLambdadPhidEtaPtLCentBckg[k]->Fill(dPhi,dEta,zv);
1977                   if( lPtArmV0 < TMath::Abs(0.2*lAlphaV0) )
1978                     fLambdadPhidEtaPtLCentBckg2[k]->Fill(dPhi,dEta,zv);
1979                 }
1980
1981                 if(radio < 0.4){ // Under the peak
1982                   //fHistArmPodBckg->Fill(lAlphaV0,lPtArmV0,3);
1983                   fLambdaBckgDecLength->Fill(dlL,ptTrig);
1984                   fLambdaBckgDCADaugToPrimVtx->Fill(dcaPos,dcaNeg,ptTrig);
1985                   fLambdaBckgEtaPhi->Fill(lPhi,lEta);
1986                   fLambdaBckgPhiRadio->Fill(lPhi,lt);
1987
1988                   // MC Association of daughter particles
1989                   for(Int_t ii=0;ii<5;ii++){
1990                     if(lMCAssocNegDaug == ii)
1991                       fLambdaBckgDCANegDaugToPrimVtx->Fill(ii,dcaNeg);
1992
1993                     if(lMCAssocPosDaug == ii)
1994                       fLambdaBckgDCAPosDaugToPrimVtx->Fill(ii,dcaPos);
1995                   }
1996
1997                   RecCascade(trkTrig,ntrack,ptrack,"Lambda");
1998                   RecCascade(trkTrig,ptrack,ntrack,"Lambda");
1999
2000                 }// End salection in the correlation peak
2001
2002               }// End selection pt bi
2003
2004           }// End triggered selection
2005
2006         }// End selection in outside the mass cut
2007
2008         break; // End Lambda selection for Correlation
2009       default:
2010         Printf(" Selection of 'step' is not set properly");
2011         break;
2012
2013       }// End switch
2014
2015     } // End Lambda selection
2016
2017     // *******************
2018     // AntiLambda selection
2019     // *******************
2020     if ( ctL && (TMath::Abs(v0->RapLambda())<fYMax)  &&
2021          (nsigNegPion < fNSigma) && (nsigPosProton < fNSigma)
2022          && !isCandidate2K0s && !isCandidate2LambdaBar) {
2023
2024       switch(step) {
2025       case kTriggerCheck:
2026
2027         if (isCandidate2LambdaBar){
2028           if(pt>ptTrig) {
2029             isV0LP = 1;
2030             ptV0LP = pt;
2031           }
2032
2033           isTrigFromV0daug = isSameTrk;
2034           if(isTrigFromV0daug)
2035             Printf("  The LP has the same momentum in X and Y as one of the AntiLambda daughters *** iV0 %d",iV0);
2036
2037         }
2038         break; // End Lambda selection for Correlation
2039       case kCorrelation:
2040         break;
2041       default:
2042         Printf( " Selection of 'step' is not set properly");
2043         break;
2044       }// End switch
2045
2046     } // End AntiLambda selection
2047
2048   } // End V0 loop
2049
2050   // ----------------------------
2051
2052   TArrayD* aV0Loop = new TArrayD(3);
2053   aV0Loop->AddAt(isV0LP,0);
2054   aV0Loop->AddAt(isTrigFromV0daug,1);
2055   aV0Loop->AddAt(ptV0LP,2);
2056
2057   return aV0Loop;
2058   delete aV0Loop;
2059 }
2060
2061 //___________________________________________________________________________________________
2062
2063 TArrayD* AliAnalysisTaskLambdaOverK0sJets::TriggerParticle()
2064 {
2065   // Obtain the trigger particle of the event to perform the correlations in phi and eta
2066
2067   Double_t ptTrigger  = -1000.;
2068   Int_t    iTrigger   = -1000;
2069   Double_t phiTrigger = -1000.;
2070   Double_t etaTrigger = -1000.;
2071   Int_t    isTriggerFromV0Daug = -1000;
2072   Int_t    isV0LP  = -1000;
2073   Double_t ptV0LP  = -100.;
2074   Int_t    isSndCheck = 0;
2075
2076   Int_t    idSndTrigger = 0;
2077
2078   // ----------------------------
2079   // 1. First trigger particle
2080
2081   Int_t nTrk= fAOD->GetNumberOfTracks();
2082   for (Int_t i=0; i<nTrk; i++) {
2083
2084     AliAODTrack *t = fAOD->GetTrack(i);
2085     if (t->IsMuonTrack()) continue;
2086     if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;
2087     if (!AcceptTrack(t)) continue;
2088
2089     Double_t xyz[3];
2090     // position of first point on track or dca
2091     if (t->GetPosition(xyz)) continue;
2092     if (TMath::Abs(xyz[0])>3.) continue;
2093     if (TMath::Abs(xyz[1])>3.) continue;
2094
2095     Double_t pt=t->Pt(),pz=t->Pz();
2096     if (TMath::Abs(pz/pt)>0.8) continue;
2097     if (TMath::Abs(t->Eta())>0.8 ) continue;
2098     /// Do I need another cuts for the primary tracks?
2099
2100     if(pt>ptTrigger) {
2101       ptTrigger = pt;
2102       iTrigger = i;
2103       //rapTrigger   = 0.5*TMath::Log((t->E()+t->Pz())/(t->E()-t->Pz()+1.e-13));
2104       phiTrigger = t->Phi();
2105       etaTrigger = t->Eta();
2106     }
2107
2108   }
2109
2110   TArrayD* v0Audit;
2111   AliAODTrack *trackLP;
2112
2113   if(iTrigger<0)
2114     goto endLP;
2115
2116   // ----------------------------
2117
2118   // 2. Checking if the trigger particle
2119   // might be a daughter from the V0-candidate
2120
2121   trackLP = (AliAODTrack*)fAOD->GetTrack(iTrigger);
2122
2123   v0Audit = V0Loop(trackLP,kTriggerCheck,kFALSE);
2124   isV0LP               = TMath::FloorNint( v0Audit->At(0) );
2125   isTriggerFromV0Daug  = TMath::FloorNint( v0Audit->At(1) );
2126   ptV0LP               = v0Audit->At(2);
2127
2128   fCheckTriggerFromV0Daug->Fill(isTriggerFromV0Daug);
2129
2130   if(isV0LP) fTriggerIsV0->Fill(ptV0LP);
2131
2132   if(isTriggerFromV0Daug && isV0LP)
2133     fCheckTriggerFromV0Daug->Fill(2);
2134
2135   if(isTriggerFromV0Daug)
2136     fTriggerComingFromDaug->Fill(ptTrigger);
2137   else
2138     goto endLP;
2139
2140   // ----------------------------
2141   // 3. Second trigger particle: when isV0LP is true
2142
2143   isSndCheck = 1;
2144   ptTrigger = -1000.;
2145
2146   for (Int_t i=0; i<nTrk; i++) {
2147     if(i == iTrigger) continue;
2148
2149     AliAODTrack *t=fAOD->GetTrack(i);
2150     if (t->IsMuonTrack()) continue;
2151     if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;
2152     if (!AcceptTrack(t)) continue;
2153
2154     Double_t xyz[3];
2155     // position of first point on track or dca
2156     if (t->GetPosition(xyz)) continue;
2157     if (TMath::Abs(xyz[0])>3.) continue;
2158     if (TMath::Abs(xyz[1])>3.) continue;
2159
2160     Double_t pt=t->Pt(),pz=t->Pz();
2161     if (TMath::Abs(pz/pt)>0.8) continue;
2162     if (TMath::Abs(t->Eta())>0.8 ) continue;
2163     /// Do I need another cuts for the primary tracks?
2164
2165     if(pt>ptTrigger) {
2166       ptTrigger = pt;
2167       idSndTrigger = i;
2168       //rapTrigger   = 0.5*TMath::Log((t->E()+t->Pz())/(t->E()-t->Pz()+1.e-13));
2169       phiTrigger = t->Phi();
2170       etaTrigger = t->Eta();
2171     }
2172
2173   }
2174
2175   iTrigger = idSndTrigger;
2176   trackLP =  (AliAODTrack*)fAOD->GetTrack(iTrigger);
2177
2178   v0Audit = V0Loop(trackLP,kTriggerCheck,kFALSE);
2179   //isV0LP           = TMath::FloorNint( v0Audit->At(0) );
2180   isTriggerFromV0Daug = TMath::FloorNint( v0Audit->At(1) );
2181   //ptV0LP           = v0Audit->At(2);
2182
2183
2184   // ----------------------------
2185   // 4. Fill the values for the array
2186
2187   endLP:
2188
2189   Double_t idTrigger = 1.*iTrigger;
2190   Double_t flagTriggerFromV0daug = 1.*isTriggerFromV0Daug;
2191   Double_t flagV0LP = 1.*isV0LP;
2192   Double_t flagSndLoop = 1.*isSndCheck;
2193
2194   TArrayD* myTrigger = new TArrayD(4);
2195   myTrigger->AddAt(idTrigger,0);
2196   myTrigger->AddAt(flagTriggerFromV0daug,1);
2197   myTrigger->AddAt(flagV0LP,2);
2198   myTrigger->AddAt(flagSndLoop,3);
2199
2200   return myTrigger;
2201   delete myTrigger;
2202   delete v0Audit;
2203   trackLP->Delete();
2204 }
2205
2206 //___________________________________________________________________________________________
2207
2208 void AliAnalysisTaskLambdaOverK0sJets::UserExec(Option_t *)
2209 {
2210   //Main loop for the Analysis
2211
2212   fAOD = (AliAODEvent *)InputEvent();
2213   fEvents->Fill(0); //event counter
2214
2215   if (!fAOD) {
2216     Printf("ERROR: aod not available");
2217     return;
2218   }
2219   fEvents->Fill(1);
2220
2221   // Physics selection
2222   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2223   AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
2224   UInt_t maskIsSelected = hdr->IsEventSelected();
2225   Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
2226   if (!isSelected) return;
2227   fEvents->Fill(2);
2228
2229   // Centrality selection
2230   AliCentrality *cent = fAOD->GetCentrality();
2231   Double_t centrality = cent->GetCentralityPercentile("V0M");
2232   fCentrality->Fill(centrality);
2233
2234   if (!cent->IsEventInCentralityClass(fCentMin,fCentMax,"V0M")) return;
2235   fEvents->Fill(3);
2236
2237   // Primary vertex
2238   const AliAODVertex *vtx = fAOD->GetPrimaryVertex();
2239   if (vtx->GetNContributors()<3) return;
2240   fEvents->Fill(4);
2241
2242   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
2243
2244   fPrimaryVertexX->Fill(xv);
2245   fPrimaryVertexY->Fill(yv);
2246   fPrimaryVertexZ->Fill(zv);
2247
2248   if (TMath::Abs(zv) > 10.) return;
2249   fEvents->Fill(5);
2250
2251   Int_t nVtxSPD = fAOD->GetNumberOfPileupVerticesSPD();
2252   Int_t nVtxTracks = fAOD->GetNumberOfPileupVerticesTracks();
2253
2254   fNumberPileUp->Fill(nVtxSPD,nVtxTracks);
2255
2256   // ***********************************************
2257   // Multiplicity with charged-tracks
2258
2259   fPIDResponse = hdr->GetPIDResponse();
2260
2261   Int_t nTrk = fAOD->GetNumberOfTracks();
2262   Int_t mult=0;
2263   Double_t nsig=0;
2264
2265   for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
2266
2267     AliAODTrack *track = fAOD->GetTrack(iTrk);
2268     if (track->IsMuonTrack()) continue;
2269     if (!track->IsOn(AliAODTrack::kTPCrefit)) continue;
2270
2271     Double_t xyz[3];
2272     if (track->GetPosition(xyz)) continue;
2273     if (TMath::Abs(xyz[0])>3.) continue;
2274     if (TMath::Abs(xyz[1])>3.) continue;
2275
2276     Double_t pt = track->Pt(), pz = track->Pz();
2277     if (TMath::Abs(pz/pt)>0.8) continue;
2278
2279     if( TMath::Abs(track->Eta())>0.8 ) continue;
2280
2281     mult++;
2282
2283     if(fUsePID){
2284       const AliAODPid *pid = track->GetDetPid();
2285       if (!pid) continue;
2286
2287       Double_t p = pid->GetTPCmomentum();
2288       Double_t dedx = pid->GetTPCsignal()/47.; //47?
2289       fdEdx->Fill(p,dedx,1);
2290
2291       nsig = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
2292       if ( TMath::Abs(nsig) < 3. )
2293         fdEdxPid->Fill(p,dedx,1);
2294     }
2295
2296   }
2297
2298   if( (mult<2300) && (centrality<2) )
2299     fEvents->Fill(6);
2300
2301   fCentMult->Fill(centrality,mult);
2302
2303   // **********************************************
2304   // Triggered Particle -  Trigger Particle
2305
2306   TArrayD* triggerArray = TriggerParticle();
2307   Int_t    iTrigger     = TMath::FloorNint( triggerArray->At(0) );
2308   Int_t    isTriggerFromV0daug = TMath::FloorNint( triggerArray->At(1) );
2309   Int_t    isV0LP    = TMath::FloorNint( triggerArray->At(2) );
2310   Int_t    is2ndLoop = TMath::FloorNint( triggerArray->At(3) );
2311
2312   Bool_t isTriggered = kFALSE;
2313   Double_t ptTrigger  = -100.;
2314   Double_t phiTrigger = -100.;
2315   Double_t etaTrigger = -100.;
2316   AliAODTrack *trkTrigger = 0x0;
2317
2318   if(iTrigger<0)
2319     fEvents->Fill(7);
2320   else {
2321     trkTrigger = (AliAODTrack*)fAOD->GetTrack(iTrigger);
2322     ptTrigger  = trkTrigger->Pt();
2323     phiTrigger = trkTrigger->Phi();
2324     etaTrigger = trkTrigger->Eta();
2325   }
2326
2327   // -------------------
2328
2329   if(!isTriggerFromV0daug && is2ndLoop)
2330     fCheckTriggerFromV0Daug->Fill(-1);
2331   else if(isTriggerFromV0daug && isV0LP && is2ndLoop)
2332     fCheckTriggerFromV0Daug->Fill(3);
2333
2334   // The highest-pt particle is a V0-candidate
2335   if(isV0LP)
2336     fEvents->Fill(8);
2337
2338   // The highest-pt particle is a V0-candidate's daughter
2339   if(isTriggerFromV0daug)
2340     fEvents->Fill(9);
2341
2342   // The sencond highest-pt particle is a V0-candidate's daughter
2343   if( isTriggerFromV0daug && is2ndLoop )
2344     fEvents->Fill(10);
2345
2346   if(!isTriggerFromV0daug){
2347     fTriggerEtaPhi->Fill(phiTrigger,etaTrigger);
2348     fTriggerPtCent->Fill(ptTrigger,centrality);
2349   }
2350
2351   if( (ptTrigger>=fTrigPtMin)  && (ptTrigger<=fTrigPtMax) &&
2352       (TMath::Abs(etaTrigger)<fTrigEtaMax) && !isTriggerFromV0daug ) {
2353     fEvents->Fill(11);
2354     isTriggered=kTRUE;
2355   }
2356   else
2357     fEvents->Fill(12);
2358
2359
2360   // ******************************************
2361   // Start loop over MC particles
2362
2363   TClonesArray *stack = 0x0;
2364   Double_t mcXv=0., mcYv=0., mcZv=0.;
2365
2366   if(fIsMC) {
2367
2368     TList *lst = fAOD->GetList();
2369     stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
2370     if (!stack) {
2371       Printf("ERROR: stack not available");
2372       return;
2373     }
2374
2375     AliAODMCHeader *mcHdr =
2376       (AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
2377
2378     mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();
2379
2380     Int_t nTrkMC = stack->GetEntriesFast();
2381     // -----------------------------------------
2382     // --------- Trigger particle --------------
2383     // -----------------------------------------
2384
2385     Double_t triggerMCPt = -1000.;
2386     Double_t triggerMCPhi = -1000.;
2387     Double_t triggerMCEta = -1000.;
2388     Bool_t   isTriggeredMC = kFALSE;
2389
2390     for (Int_t iTrkMC = 0; iTrkMC < nTrkMC; iTrkMC++){
2391
2392       AliAODMCParticle *p0 = (AliAODMCParticle*)stack->At(iTrkMC);
2393       if(!p0) continue;
2394       if(TMath::Abs(p0->Eta())>0.8) continue;
2395       if(p0->Pt()<0.15) continue;
2396
2397       // ----------------------------------------
2398
2399       // For injected MC: determine where HIJING event ends
2400       if (fEndOfHijingEvent==-1) {
2401         if ( ( p0->GetStatus() == 21 ) ||
2402              ( (p0->GetPdgCode() == 443) &&
2403                (p0->GetMother() == -1)   &&
2404                (p0->GetDaughter(0) ==  (iTrkMC+1))) ) {
2405           fEndOfHijingEvent = iTrkMC;
2406         }
2407       }
2408
2409       // ----------------------------------------
2410
2411       Int_t isNaturalPart = 1;
2412       if ( (iTrkMC>=fEndOfHijingEvent) &&
2413            (fEndOfHijingEvent!=-1)     &&
2414            (p0->GetMother()<0) )
2415         isNaturalPart = 0;
2416
2417       // ----------------------------------------
2418
2419       Int_t lPdgcodeCurrentPart = p0->GetPdgCode();
2420       if ( (lPdgcodeCurrentPart == kK0Short) ||
2421            (lPdgcodeCurrentPart == kLambda0) ||
2422            (lPdgcodeCurrentPart == kLambda0Bar) ) continue;
2423
2424       // ----------------------------------------
2425
2426       if(isNaturalPart == 0) continue;
2427       if( !p0->IsPhysicalPrimary() ) continue;
2428
2429       Double_t ptPrim = p0->Pt();
2430       if(ptPrim>triggerMCPt){
2431
2432         triggerMCPt  = p0->Pt();
2433         triggerMCPhi = p0->Phi();
2434         triggerMCEta = p0->Eta();
2435
2436       } // End trigger selection
2437
2438     } // End loop over charged particles
2439
2440
2441     // ----------------------------------------------------------------
2442     // Check if Primary particle has larger pt than strange particles
2443     // ----------------------------------------------------------------
2444
2445     Bool_t isV0LPMC = kFALSE;
2446     if(triggerMCPt>1.)
2447       for (Int_t iTrkMC = 0; iTrkMC < nTrkMC; iTrkMC++){
2448
2449         AliAODMCParticle *p0 = (AliAODMCParticle*)stack->At(iTrkMC);
2450         if(!p0) continue;
2451
2452         // ----------------------------------------
2453
2454         Int_t lPdgcodeCurrentPart = p0->GetPdgCode();
2455         if ( (lPdgcodeCurrentPart != kK0Short) &&
2456              (lPdgcodeCurrentPart != kLambda0) &&
2457              (lPdgcodeCurrentPart != kLambda0Bar) ) continue;
2458
2459         // ----------------------------------------
2460
2461         Int_t isNaturalPart = 1;
2462         if ( (iTrkMC>=fEndOfHijingEvent) &&
2463              (fEndOfHijingEvent!=-1)     &&
2464              (p0->GetMother()<0) )
2465           isNaturalPart = 0;
2466
2467         if(fSeparateInjPart)
2468           if(isNaturalPart == 0) continue;
2469
2470         // ----------------------------------------
2471
2472         Double_t lRapCurrentPart = MyRapidity(p0->E(),p0->Pz());
2473         Double_t lPtCurrentPart  = p0->Pt();
2474
2475         Int_t iCurrentMother = p0->GetMother();
2476         AliAODMCParticle *pCurrentMother = (AliAODMCParticle *)stack->At(iCurrentMother);
2477         Int_t lPdgCurrentMother = 0;
2478         if (iCurrentMother == -1) { lPdgCurrentMother = 0;}
2479         else { lPdgCurrentMother = pCurrentMother->GetPdgCode(); }
2480
2481         Int_t id0  = p0->GetDaughter(0);
2482         Int_t id1  = p0->GetDaughter(1);
2483
2484         //if ( id0 ==  id1 ) continue;
2485         if ( (id0 < 0 || id1 < 0) ||
2486              (id0 >=nTrkMC  || id1 >= nTrkMC) ) continue;
2487
2488         AliAODMCParticle *pDaughter0 = (AliAODMCParticle *)stack->UncheckedAt(id0);
2489         AliAODMCParticle *pDaughter1 = (AliAODMCParticle *)stack->UncheckedAt(id1);
2490         if (!pDaughter0 || !pDaughter1) continue;
2491
2492         if ( TMath::Abs(pDaughter0->Eta()) > fMaxEtaDaughter ||
2493              TMath::Abs(pDaughter1->Eta()) > fMaxEtaDaughter )
2494           continue;
2495
2496         // Daughter momentum cut: ! FIX it in case of AOD !
2497         if ( ( pDaughter0->Pt() < fMinPtDaughter ) ||
2498              ( pDaughter1->Pt() < fMinPtDaughter )  )
2499           continue;
2500
2501         if ((p0->Pt())<pMin) continue;
2502         if (TMath::Abs(lRapCurrentPart) > fYMax)  continue;
2503
2504         Double_t dx = mcXv-p0->Xv(),  dy = mcYv-p0->Yv(),  dz = mcZv-p0->Zv();
2505         Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2506
2507         //Cut in the 3D-distance of the secondary vertex to primary vertex
2508         if (l > 0.01) continue; // secondary V0
2509
2510         // is V0 the LP?
2511         if(lPtCurrentPart>triggerMCPt)
2512           isV0LPMC = kTRUE;
2513
2514       }
2515
2516     // If V0 is not the particle with highest pt
2517     // then the event is triggered
2518     if(!isV0LPMC)
2519       fTriggerMCPtCent->Fill(triggerMCPt,triggerMCEta,centrality);
2520
2521     if( (triggerMCPt>=fTrigPtMin)  && (triggerMCPt<=fTrigPtMax) &&
2522         (TMath::Abs(triggerMCEta)<fTrigEtaMax) &&
2523         (isV0LPMC==kFALSE) )
2524       isTriggeredMC=kTRUE;
2525
2526
2527     // -----------------------------------------
2528     // ---------- Strange particles ------------
2529     // -----------------------------------------
2530     Double_t dPhiMC = -100.;
2531     Double_t dEtaMC = -100.;
2532     for (Int_t iTrkMC = 0; iTrkMC < nTrkMC; iTrkMC++){
2533
2534       AliAODMCParticle *p0 = (AliAODMCParticle*)stack->At(iTrkMC);
2535       if(!p0) continue;
2536
2537       // ----------------------------------------
2538
2539       Int_t lPdgcodeCurrentPart = p0->GetPdgCode();
2540       if ( (lPdgcodeCurrentPart != kK0Short) &&
2541            (lPdgcodeCurrentPart != kLambda0) &&
2542            (lPdgcodeCurrentPart != kLambda0Bar) ) continue;
2543
2544       // ----------------------------------------
2545
2546       Int_t isNaturalPart = 1;
2547       if ( (iTrkMC>=fEndOfHijingEvent) &&
2548            (fEndOfHijingEvent!=-1)     &&
2549            (p0->GetMother()<0) )
2550         isNaturalPart = 0;
2551
2552       fInjectedParticles->Fill(isNaturalPart);
2553
2554       if(fSeparateInjPart)
2555         if(isNaturalPart == 0) continue;
2556
2557       // ----------------------------------------
2558
2559       Double_t lRapCurrentPart = MyRapidity(p0->E(),p0->Pz());
2560       Double_t lEtaCurrentPart = p0->Eta();
2561       Double_t lPhiCurrentPart = p0->Phi();
2562       Double_t lPtCurrentPart  = p0->Pt();
2563
2564       Int_t iCurrentMother = p0->GetMother();
2565       AliAODMCParticle *pCurrentMother = (AliAODMCParticle *)stack->At(iCurrentMother);
2566       Int_t lPdgCurrentMother = 0;
2567       if (iCurrentMother == -1) { lPdgCurrentMother = 0;}
2568       else { lPdgCurrentMother = pCurrentMother->GetPdgCode(); }
2569
2570       Int_t id0  = p0->GetDaughter(0);
2571       Int_t id1  = p0->GetDaughter(1);
2572
2573       //if ( id0 ==  id1 ) continue;
2574       if ( (id0 < 0 || id1 < 0) ||
2575            (id0 >=nTrkMC  || id1 >= nTrkMC) ) continue;
2576
2577       AliAODMCParticle *pDaughter0 = (AliAODMCParticle *)stack->UncheckedAt(id0);
2578       AliAODMCParticle *pDaughter1 = (AliAODMCParticle *)stack->UncheckedAt(id1);
2579       if (!pDaughter0 || !pDaughter1) continue;
2580
2581       if ( TMath::Abs(pDaughter0->Eta()) > fMaxEtaDaughter ||
2582            TMath::Abs(pDaughter1->Eta()) > fMaxEtaDaughter )
2583         continue;
2584
2585       // Daughter momentum cut: ! FIX it in case of AOD !
2586       if ( ( pDaughter0->Pt() < fMinPtDaughter ) ||
2587            ( pDaughter1->Pt() < fMinPtDaughter )  )
2588         continue;
2589
2590       if ((p0->Pt())<pMin) continue;
2591       if (TMath::Abs(lRapCurrentPart) > fYMax)  continue;
2592
2593       Double_t dx = mcXv-p0->Xv(),  dy = mcYv-p0->Yv(),  dz = mcZv-p0->Zv();
2594       Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2595
2596       //Cut in the 3D-distance of the secondary vertex to primary vertex
2597       if (l > 0.01) continue; // secondary V0
2598
2599       //Transverse distance to vertex
2600       dx = mcXv-pDaughter0->Xv(); dy = mcYv-pDaughter0->Yv();
2601       Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
2602
2603       // K0s
2604       if (lPdgcodeCurrentPart == kK0Short) {
2605         //if (multiplicity>=1){
2606         fK0sMCPt->Fill(lPtCurrentPart);
2607         fK0sMCPtRap->Fill(lPtCurrentPart,lRapCurrentPart);
2608         fK0sMCPtEta->Fill(lPtCurrentPart,lEtaCurrentPart);
2609         fK0sMCPtLt->Fill(lPtCurrentPart,lt,centrality);
2610         //}
2611
2612
2613         // Triggered Event
2614         if(isTriggeredMC){
2615
2616           dPhiMC = dPHI(triggerMCPhi,lPhiCurrentPart);
2617           dEtaMC = triggerMCEta - lEtaCurrentPart;
2618
2619           // Pt bin
2620           for(Int_t k=0;k<kN1;k++)
2621             if( (lPtCurrentPart>kPtBinV0[k]) && (lPtCurrentPart<kPtBinV0[k+1]) ){
2622               fK0sdPhidEtaMC[k]->Fill(dPhiMC,dEtaMC,zv);
2623
2624               if(centrality<10)
2625                 fK0sdPhidEtaMCCent[k]->Fill(dPhiMC,dEtaMC,zv);
2626
2627             } // End pt bin
2628
2629         } // End trigger event
2630
2631       } // End K0s selection
2632
2633       // Lambda
2634       if (lPdgcodeCurrentPart == kLambda0) {
2635         //if (multiplicity>=1){
2636         fLambdaMCPt->Fill(lPtCurrentPart);
2637         fLambdaMCPtRap->Fill(lPtCurrentPart,lRapCurrentPart);
2638         fLambdaMCPtEta->Fill(lPtCurrentPart,lEtaCurrentPart);
2639         fLambdaMCPtLt->Fill(lPtCurrentPart,lt,centrality);
2640         //}
2641
2642         // Triggered Event
2643         if(isTriggeredMC){
2644
2645           dPhiMC = dPHI(triggerMCPhi,lPhiCurrentPart);
2646           dEtaMC = triggerMCEta - lEtaCurrentPart;
2647
2648           // Pt bin
2649           for(Int_t k=0;k<kN1;k++)
2650             if( (lPtCurrentPart>kPtBinV0[k]) && (lPtCurrentPart<kPtBinV0[k+1]) ){
2651
2652               fLambdadPhidEtaMC[k]->Fill(dPhiMC,dEtaMC,zv);
2653
2654               if(centrality<10)
2655                 fLambdadPhidEtaMCCent[k]->Fill(dPhiMC,dEtaMC,zv);
2656
2657             } //End pt bin
2658
2659         } // End trigger event
2660
2661       } // End Lambda
2662       // AntiLambda
2663       /*
2664       if (lPdgcodeCurrentPart == kLambda0Bar) {
2665         //if (multiplicity>=1){
2666
2667         fAntiLambdaMCPt->Fill(lPtCurrentPart);
2668         fAntiLambdaMCPtRap->Fill(lPtCurrentPart,lRapCurrentPart);
2669         fAntiLambdaMCPtEta->Fill(lPtCurrentPart,lEtaCurrentPart);
2670         fAntiLambdaMCPtLt->Fill(lPtCurrentPart,lt,centrality);
2671
2672         //}
2673       }
2674       */
2675     } // End loop over MC
2676
2677   } // End MC condition
2678
2679   // *************************************************
2680   // V0 loop - AOD
2681
2682   TArrayD* v0Corr = V0Loop(trkTrigger,kCorrelation,isTriggered);
2683
2684   delete v0Corr;
2685   delete triggerArray;
2686
2687 }
2688
2689 //___________________________________________________________________________________________
2690
2691 void AliAnalysisTaskLambdaOverK0sJets::Terminate(Option_t *)
2692 {
2693   // The Terminate() function is the last function to be called during
2694   // a query. It always runs on the client, it can be used to present
2695   // the results graphically or save the results to file.
2696
2697   fOutput=(TList*)GetOutputData(1);
2698   if (fOutput) {
2699     Printf("\n\t *** DONE: fOutput available *** \n");
2700     return;
2701   }
2702   else{
2703     Printf("ERROR: fOutput not available");
2704     return;
2705   }
2706
2707 }