PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHdibaryonLPpi.cxx

   1 /**************************************************************************
   2  * Author : Benjamin Dönigus (benjamin.doenigus@cern.ch)                  *
   3  *                                                                        *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //-----------------------------------------------------------------
  16 //                 AliAnalysisTaskHdibaryonLPpi class
  17 //          used to search for the H-Dibaryon in weak
  18 //          (Lambda Proton Pion) and strong (Lambda Lambda) decays
  19 //-----------------------------------------------------------------
  20
  21 #include "Riostream.h"
  22 #include "TROOT.h"
  23 #include "TChain.h"
  24 #include "TStyle.h"
  25 #include "TSystem.h"
  26 #include "TTree.h"
  27 #include "TH1F.h"
  28 #include "TH1.h"
  29 #include "TH2D.h"
  30 #include "TH3F.h"
  31 #include "TCanvas.h"
  32 #include "TParticle.h"
  33 #include "TNtuple.h"
  34 #include "TObjString.h"
  35 #include "TLorentzVector.h"
  36 #include "TDatabasePDG.h"
  37
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliMCEventHandler.h"
  41 #include "AliMCEvent.h"
  42
  43 #include "AliESD.h"
  44 #include "AliESDpid.h"
  45 #include "AliESDEvent.h"
  46 #include "AliGenEventHeader.h"
  47 #include "AliESDInputHandler.h"
  48 #include "AliLog.h"
  49 #include "AliStack.h"
  50 #include "AliHeader.h"
  51
  52 #include "AliKFParticle.h"
  53 #include "AliKFVertex.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliESDv0Cuts.h"
  56 #include "AliAnalysisTaskHdibaryonLPpi.h"
  57
  58 #include "TFile.h"
  59 #include "TH2F.h"
  60 #include "TF1.h"
  61 #include "TCanvas.h"
  62 #include "TList.h"
  63 #include "TParallelCoord.h"
  64
  65 #include "AliMCParticle.h"
  66 #include "AliGenPythiaEventHeader.h"
  67
  68 #include "AliPID.h"
  69 #include "AliESDtrack.h"
  70 #include "AliCentrality.h"
  71
  72 #include "THnSparse.h"
  73
  74 #include "AliVertexerTracks.h"
  75
  76 using namespace std;
  77 ClassImp(AliAnalysisTaskHdibaryonLPpi)
  78 //________________________________________________________________________
  79 AliAnalysisTaskHdibaryonLPpi::AliAnalysisTaskHdibaryonLPpi() : AliAnalysisTaskSE()/*AliAnalysisTask(name, ""), fMCEvent(0)*/, fESD(0),   fESDtrackCutsV0(0),
  80   fESDCutsV0(0),
  81   fEsdTrackCuts(0),
  82   fBin(0),
  83   fEvent(0x0),
  84   fHistList(0),
  85   fHistMassDPi(0),
  86   fHistMassLPi(0),
  87   fHistMassLambdaPi(0),
  88   fHistMassLambda(0),
  89   fHistMassLambdaPPi(0),
  90   fHistMassLambdaP(0),
  91   fHistMassLambdaK(0),
  92   fHistMassK0onl(0),
  93   fHistMassK0offl(0),
  94   fHistMassK0onlC(0),
  95   fHistMassK0offlC(0),
  96   fHistMassPQonl(0),
  97   fHistMassPQoffl(0),
  98   fHistDC(0),
  99   fHistArmenterosPodolanski(0),
 100   fHistArmenterosPodolanskiCut(0),
 101   fHistHDibaryonInvaMassGen(0),
 102   fHistHDibaryonInvaMassGenRes(0),
 103   fHistAntiHDibaryonInvaMassGen(0),
 104   fHistHDibaryonInvaMassAso(0),
 105   fHistHDibaryonInvaMassAsoReso(0),
 106   fHistAntiHDibaryonInvaMassAso(0),
 107   fHistCheck(0),
 108   fHistHPointingAngle(0),
 109   fHistMassH(0),
 110   fHistMassLambdaFromH(0),
 111   fHistMassLambdaFromHtLorentz(0),
 112   fHistMassPpi(0),
 113   fHistMassPpiReso(0),
 114   fHistMassLpi(0),
 115   fHistMassLP(0),
 116   fHistProtonPIDBb(0),
 117   fHistPionPIDBb(0),
 118   fHistProtonPIDLambda(0),
 119   fHistPionPIDLambda(0),
 120   fHistMCdcaPvtxDvtx(0),
 121   fHistMCdcaPvtxLvtx(0),
 122   fHistMCdcaDvtxLvtx(0),
 123   fHistMCangleLH(0),
 124   fHistMCdecayAngle(0),
 125   fHistMCpointingAngle(0),
 126   fHistMCap(0),
 127   fHistMCdcaPvtxDvtxReso(0),
 128   fHistMCdcaPvtxLvtxReso(0),
 129   fHistMCdcaDvtxLvtxReso(0),
 130   fHistMCangleLHReso(0),
 131   fHistMCdecayAngleReso(0),
 132   fHistMCpointingAngleReso(0),
 133   fHistMCapReso(0),
 134   fHistCentrality(0),
 135   fHistCentralityAC(0),
 136   fHistMultiplicity(0),
 137   fHistHilf1(0),
 138   fHistHilf2(0),
 139   fHistHilf3(0),
 140   fHistHilf4(0),
 141   fHistHilf5(0),
 142   fHistHilf6(0),
 143   fHistPtvsEtaGen(0),
 144   fHistPtvsEtaAso(0),
 145   fHistPtvsYGen(0),
 146   fHistPtvsYAso(0),
 147   fHistRap(0),
 148   fHistCount(0),
 149   fPIDtpcESD(0),
 150   fHistTriggerStat(0),
 151   fHistTriggerStatAfterEventSelection(0),
 152   fHistMassHcentMult(0),
 153   fHistNdim(0)
 154 {
 155   // DefaultConstructor
 156
 157 }
 158
 159 //________________________________________________________________________
 160 AliAnalysisTaskHdibaryonLPpi::AliAnalysisTaskHdibaryonLPpi(const char *name) : AliAnalysisTaskSE(name)/*AliAnalysisTask(name, ""), fMCEvent(0)*/, fESD(0),   fESDtrackCutsV0(0),
 161   fESDCutsV0(0),
 162   fEsdTrackCuts(0),
 163   fBin(0),
 164   fEvent(0x0),
 165   fHistList(0),
 166   fHistMassDPi(0),
 167   fHistMassLPi(0),
 168   fHistMassLambdaPi(0),
 169   fHistMassLambda(0),
 170   fHistMassLambdaPPi(0),
 171   fHistMassLambdaP(0),
 172   fHistMassLambdaK(0),
 173   fHistMassK0onl(0),
 174   fHistMassK0offl(0),
 175   fHistMassK0onlC(0),
 176   fHistMassK0offlC(0),
 177   fHistMassPQonl(0),
 178   fHistMassPQoffl(0),
 179   fHistDC(0),
 180   fHistArmenterosPodolanski(0),
 181   fHistArmenterosPodolanskiCut(0),
 182   fHistHDibaryonInvaMassGen(0),
 183   fHistHDibaryonInvaMassGenRes(0),
 184   fHistAntiHDibaryonInvaMassGen(0),
 185   fHistHDibaryonInvaMassAso(0),
 186   fHistHDibaryonInvaMassAsoReso(0),
 187   fHistAntiHDibaryonInvaMassAso(0),
 188   fHistCheck(0),
 189   fHistHPointingAngle(0),
 190   fHistMassH(0),
 191   fHistMassLambdaFromH(0),
 192   fHistMassLambdaFromHtLorentz(0),
 193   fHistMassPpi(0),
 194   fHistMassPpiReso(0),
 195   fHistMassLpi(0),
 196   fHistMassLP(0),
 197   fHistProtonPIDBb(0),
 198   fHistPionPIDBb(0),
 199   fHistProtonPIDLambda(0),
 200   fHistPionPIDLambda(0),
 201   fHistMCdcaPvtxDvtx(0),
 202   fHistMCdcaPvtxLvtx(0),
 203   fHistMCdcaDvtxLvtx(0),
 204   fHistMCangleLH(0),
 205   fHistMCdecayAngle(0),
 206   fHistMCpointingAngle(0),
 207   fHistMCap(0),
 208   fHistMCdcaPvtxDvtxReso(0),
 209   fHistMCdcaPvtxLvtxReso(0),
 210   fHistMCdcaDvtxLvtxReso(0),
 211   fHistMCangleLHReso(0),
 212   fHistMCdecayAngleReso(0),
 213   fHistMCpointingAngleReso(0),
 214   fHistMCapReso(0),
 215   fHistCentrality(0),
 216   fHistCentralityAC(0),
 217   fHistMultiplicity(0),
 218   fHistHilf1(0),
 219   fHistHilf2(0),
 220   fHistHilf3(0),
 221   fHistHilf4(0),
 222   fHistHilf5(0),
 223   fHistHilf6(0),
 224   fHistPtvsEtaGen(0),
 225   fHistPtvsEtaAso(0),
 226   fHistPtvsYGen(0),
 227   fHistPtvsYAso(0),
 228   fHistRap(0),
 229   fHistCount(0),
 230   fPIDtpcESD(0),
 231   fHistTriggerStat(0),
 232   fHistTriggerStatAfterEventSelection(0),
 233   fHistMassHcentMult(0),
 234   fHistNdim(0)
 235
 236 {
 237   // Constructor
 238
 239   // Define input and output slots here
 240   // Input from a TChain
 241   DefineInput(0, TChain::Class());
 242   // Output to TList container
 243   DefineOutput(1, TList::Class()); //full
 244
 245   //MC info contol
 246   if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())SetHasMC();
 247
 248   //V0 cuts
 249
 250   fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
 251   fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
 252   fESDtrackCutsV0->SetMinNClustersTPC(80);
 253   fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
 254   fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
 255   fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
 256   fESDtrackCutsV0->SetPtRange(0.2,1.5);
 257   fESDtrackCutsV0->SetMinDCAToVertexXY(3);
 258   fESDtrackCutsV0->SetMinDCAToVertexZ(3);
 259
 260   fESDCutsV0 = new AliESDv0Cuts("AliESDCutsV0","AliESDCutsV0");
 261   fESDCutsV0->SetMaxDcaV0Daughters(1.0);
 262   fESDCutsV0->SetMinDcaNegToVertex(2.0);
 263   fESDCutsV0->SetMinDcaPosToVertex(2.0);
 264
 265   //ESD Track cuts
 266   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 267   fEsdTrackCuts->SetMinNClustersTPC(80);
 268   fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 269   fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
 270   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 271   fEsdTrackCuts->SetEtaRange(-0.9,0.9);
 272 }
 273
 274 //____________________________________________________________
 275 AliAnalysisTaskHdibaryonLPpi::~AliAnalysisTaskHdibaryonLPpi(){
 276   //
 277   // Destructor
 278   //
 279   if(fHistList){
 280     fHistList->Clear();
 281     delete fHistList;
 282   }
 283   if(fEsdTrackCuts) delete fEsdTrackCuts;
 284   if(fESDtrackCutsV0) delete fESDtrackCutsV0;
 285   if(fESDCutsV0) delete fESDCutsV0;
 286
 287 }
 288
 289 //________________________________________________________________________
 290 void AliAnalysisTaskHdibaryonLPpi::UserCreateOutputObjects()
 291 {
 292   // Create histograms
 293   // Called once
 294
 295  fHistList = new TList();
 296
 297  fHistMassDPi = new TH1F("fHistMassDPi", "Invariant mass distribution p+#pi^{-} ", 500, 1.0, 1.25);
 298  fHistMassDPi->GetXaxis()->SetTitle("Invariant mass p+#pi^{-} (GeV/c^{2})");
 299  fHistMassDPi->GetYaxis()->SetTitle("Entries");
 300  fHistMassDPi->SetMarkerStyle(kFullCircle);
 301
 302  fHistMassLPi = new TH1F("fHistMassLPi", "Offline Invariant mass distribution p+#pi^{-} ", 500, 1.0, 1.25);
 303  fHistMassLPi->GetXaxis()->SetTitle("Invariant mass p+#pi^{-} (GeV/c^{2})");
 304  fHistMassLPi->GetYaxis()->SetTitle("Entries");
 305  fHistMassLPi->SetMarkerStyle(kFullCircle);
 306
 307  fHistMassLambdaPi = new TH1F("fHistMassLambdaPi", "Invariant mass distribution #Lambda+#pi^{-} ", 500, 1.2, 1.5);
 308  fHistMassLambdaPi->GetXaxis()->SetTitle("Invariant mass #Lambda+#pi^{-} (GeV/c^{2})");
 309  fHistMassLambdaPi->GetYaxis()->SetTitle("Entries");
 310  fHistMassLambdaPi->SetMarkerStyle(kFullCircle);
 311
 312  fHistMassLambda = new TH1F("fHistMassLambda", "Invariant mass distribution of #Lambda for further analyis", 500, 1.0, 1.2);
 313  fHistMassLambda->GetXaxis()->SetTitle("Invariant mass p+#pi^{+} (GeV/c^{2})");
 314  fHistMassLambda->GetYaxis()->SetTitle("Entries");
 315  fHistMassLambda->SetMarkerStyle(kFullCircle);
 316
 317  fHistMassLambdaPPi = new TH1F("fHistMassLambdaPPi", "Invariant mass distribution #Lambdap#pi^{-} ", 300, 2.2, 2.5);
 318  fHistMassLambdaPPi->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})");
 319  fHistMassLambdaPPi->GetYaxis()->SetTitle("Entries");
 320  fHistMassLambdaPPi->SetMarkerStyle(kFullCircle);
 321
 322  fHistMassLambdaP = new TH1F("fHistMassLambdaP", "Invariant mass distribution #Lambdap ", 300, 2.2, 2.5);
 323  fHistMassLambdaP->GetXaxis()->SetTitle("Invariant mass #Lambdap (GeV/c^{2})");
 324  fHistMassLambdaP->GetYaxis()->SetTitle("Entries");
 325  fHistMassLambdaP->SetMarkerStyle(kFullCircle);
 326
 327  fHistMassLambdaK = new TH1F("fHistMassLambdaK", "Invariant mass distribution #LambdaK^{-} ", 300, 1.6, 1.9);
 328  fHistMassLambdaK->GetXaxis()->SetTitle("Invariant mass #LambdaK^{-} (GeV/c^{2})");
 329  fHistMassLambdaK->GetYaxis()->SetTitle("Entries");
 330  fHistMassLambdaK->SetMarkerStyle(kFullCircle);
 331
 332  fHistMassK0onl = new TH1F("fHistMassK0onl", "Invariant mass distribution K_{s}^{0} online V0 finder ", 400, 0.2, 1.0);
 333  fHistMassK0onl->GetXaxis()->SetTitle("Invariant mass #pi^{+}#pi^{-} (GeV/c^{2})");
 334  fHistMassK0onl->GetYaxis()->SetTitle("Entries");
 335  fHistMassK0onl->SetMarkerStyle(kFullCircle);
 336
 337  fHistMassK0offl = new TH1F("fHistMassK0offl", "Invariant mass distribution  K_{s}^{0} offline V0 finder ", 400, 0.2, 1.0);
 338  fHistMassK0offl->GetXaxis()->SetTitle("Invariant mass  #pi^{+}#pi^{-} (GeV/c^{2})");
 339  fHistMassK0offl->GetYaxis()->SetTitle("Entries");
 340  fHistMassK0offl->SetMarkerStyle(kFullCircle);
 341
 342  fHistMassK0onlC = new TH1F("fHistMassK0onlC", "Invariant mass distribution K_{s}^{0} online V0 finder ", 400, 0.2, 1.0);
 343  fHistMassK0onlC->GetXaxis()->SetTitle("Invariant mass #pi^{+}#pi^{-} (GeV/c^{2})");
 344  fHistMassK0onlC->GetYaxis()->SetTitle("Entries");
 345  fHistMassK0onlC->SetMarkerStyle(kFullCircle);
 346
 347  fHistMassK0offlC = new TH1F("fHistMassK0offlC", "Invariant mass distribution  K_{s}^{0} offline V0 finder ", 400, 0.2, 1.0);
 348  fHistMassK0offlC->GetXaxis()->SetTitle("Invariant mass  #pi^{+}#pi^{-} (GeV/c^{2})");
 349  fHistMassK0offlC->GetYaxis()->SetTitle("Entries");
 350  fHistMassK0offlC->SetMarkerStyle(kFullCircle);
 351
 352  fHistMassPQonl = new TH1F("fHistMassPQonl", "Invariant mass distribution K_{s}^{0}p using online V0 finder ", 500, 1.3, 2.3);
 353  fHistMassPQonl->GetXaxis()->SetTitle("Invariant mass K_{s}^{0}p (GeV/c^{2})");
 354  fHistMassPQonl->GetYaxis()->SetTitle("Entries");
 355  fHistMassPQonl->SetMarkerStyle(kFullCircle);
 356
 357  fHistMassPQoffl = new TH1F("fHistMassPQoffl", "Invariant mass distribution  K_{s}^{0}p using offline V0 finder ", 500, 1.3, 2.3);
 358  fHistMassPQoffl->GetXaxis()->SetTitle("Invariant mass K_{s}^{0}p (GeV/c^{2})");
 359  fHistMassPQoffl->GetYaxis()->SetTitle("Entries");
 360  fHistMassPQoffl->SetMarkerStyle(kFullCircle);
 361
 362  fHistDC = new TH1F("fHistDC", "Proper decay length", 500, 0.0, 25);
 363  fHistDC->GetXaxis()->SetTitle("c#tau (cm)");
 364  fHistDC->GetYaxis()->SetTitle("Entries");
 365  fHistDC->SetMarkerStyle(kFullCircle);
 366
 367  fHistArmenterosPodolanski = new TH2F("fHistArmenterosPodolanski", "Armenteros-Podolanski plot", 200,-1.0,1.0, 500,0,1);
 368  fHistArmenterosPodolanski->GetXaxis()->SetTitle("#alpha");
 369  fHistArmenterosPodolanski->GetYaxis()->SetTitle("q_{t}");
 370  fHistArmenterosPodolanski->SetMarkerStyle(kFullCircle);
 371
 372  fHistArmenterosPodolanskiCut = new TH2F("fHistArmenterosPodolanskiCut", "Armenteros-Podolanski plot after cut", 200,-1.0,1.0, 500,0,1);
 373  fHistArmenterosPodolanskiCut->GetXaxis()->SetTitle("#alpha");
 374  fHistArmenterosPodolanskiCut->GetYaxis()->SetTitle("q_{t}");
 375  fHistArmenterosPodolanskiCut->SetMarkerStyle(kFullCircle);
 376
 377  fHistHDibaryonInvaMassGen = new TH1F("fHistHDibaryonInvaMassGen", "Generated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
 378  fHistHDibaryonInvaMassGen->GetYaxis()->SetTitle("Counts");
 379  fHistHDibaryonInvaMassGen->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 380
 381  fHistHDibaryonInvaMassGenRes = new TH1F("fHistHDibaryonInvaMassGenRes", "Generated  #Lambda p #pi^{-} with particles in rapidity!!!", 200, 2.1, 2.3);
 382  fHistHDibaryonInvaMassGenRes->GetYaxis()->SetTitle("Counts");
 383  fHistHDibaryonInvaMassGenRes->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 384
 385  fHistAntiHDibaryonInvaMassGen = new TH1F("fHistAntiHDibaryonInvaMassGen", "Generated  #bar{#Lambda} #bar{p} #pi^{+}", 200, 2.1, 2.3);
 386  fHistAntiHDibaryonInvaMassGen->GetYaxis()->SetTitle("Counts");
 387  fHistAntiHDibaryonInvaMassGen->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 388
 389  fHistHDibaryonInvaMassAso = new TH1F("fHistHDibaryonInvaMassAso", "Associated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
 390  fHistHDibaryonInvaMassAso->GetYaxis()->SetTitle("Counts");
 391  fHistHDibaryonInvaMassAso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 392
 393  fHistHDibaryonInvaMassAsoReso = new TH1F("fHistHDibaryonInvaMassAsoReso", "Associated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
 394  fHistHDibaryonInvaMassAsoReso->GetYaxis()->SetTitle("Counts");
 395  fHistHDibaryonInvaMassAsoReso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 396
 397  fHistAntiHDibaryonInvaMassAso = new TH1F("fHistAntiHDibaryonInvaMassAso", "Associated  #bar{#Lambda} #bar{p} #pi^{+}", 200, 2.1, 2.3);
 398  fHistAntiHDibaryonInvaMassAso->GetYaxis()->SetTitle("Counts");
 399  fHistAntiHDibaryonInvaMassAso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
 400
 401  fHistCheck = new TH2F("fHistCheck", "Check online/offline", 200, -0.5, 1.5, 200, -0.5, 1.5);
 402  fHistCheck->GetXaxis()->SetTitle("offline");
 403  fHistCheck->GetYaxis()->SetTitle("online");
 404  fHistCheck->SetMarkerStyle(kFullCircle);
 405
 406  fHistHPointingAngle= new TH1F("fHistHPointingAngle", "Pointing angle distribution for #Lambdap#pi^{-}", 200, 0., 2*TMath::Pi());
 407  fHistHPointingAngle->GetXaxis()->SetTitle("Pointing angle distribution for #Lambdap#pi^{-}");
 408  fHistHPointingAngle->GetYaxis()->SetTitle("Entries");
 409  fHistHPointingAngle->SetMarkerStyle(kFullCircle);
 410
 411  fHistMassH= new TH1F("fHistMassH", "Invariant mass distribution #Lambdap#pi^{-} (GeV/c^{2}) after pointing angle cut", 3000, 2.2, 2.5);
 412  fHistMassH->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})");
 413  fHistMassH->GetYaxis()->SetTitle("Entries");
 414  fHistMassH->SetMarkerStyle(kFullCircle);
 415
 416  fHistMassLambdaFromH= new TH1F("fHistMassLambdaFromH", "Invariant mass distribution #Lambda (GeV/c^{2}) asking for Mother to be H", 300, 1.0, 1.3);
 417  fHistMassLambdaFromH->GetXaxis()->SetTitle("Invariant mass #Lambda (GeV/c^{2})");
 418  fHistMassLambdaFromH->GetYaxis()->SetTitle("Entries");
 419  fHistMassLambdaFromH->SetMarkerStyle(kFullCircle);
 420
 421  fHistMassLambdaFromHtLorentz= new TH1F("fHistMassLambdaFromHtLorentz", "Invariant mass distribution #Lambda (GeV/c^{2}) asking for Mother to be H", 300, 1.0, 1.3);
 422  fHistMassLambdaFromHtLorentz->GetXaxis()->SetTitle("Invariant mass #Lambda (GeV/c^{2})");
 423  fHistMassLambdaFromHtLorentz->GetYaxis()->SetTitle("Entries");
 424  fHistMassLambdaFromHtLorentz->SetMarkerStyle(kFullCircle);
 425
 426  fHistMassPpi= new TH1F("fHistMassPpi", "Invariant mass distribution of the p#pi^{-} used for combing with Lambda", 300, 1.0, 1.6);
 427  fHistMassPpi->GetXaxis()->SetTitle("Invariant mass p#pi^{-} (GeV/c^{2})");
 428  fHistMassPpi->GetYaxis()->SetTitle("Entries");
 429  fHistMassPpi->SetMarkerStyle(kFullCircle);
 430
 431  fHistMassPpiReso= new TH1F("fHistMassPpiReso", "Invariant mass distribution of the p#pi^{-} used for combing with Lambda", 300, 1.0, 1.6);
 432  fHistMassPpiReso->GetXaxis()->SetTitle("Invariant mass p#pi^{-} (GeV/c^{2})");
 433  fHistMassPpiReso->GetYaxis()->SetTitle("Entries");
 434  fHistMassPpiReso->SetMarkerStyle(kFullCircle);
 435
 436  fHistMassLpi= new TH1F("fHistMassLpi", "Invariant mass distribution of the #Lambda#pi^{-} used for combing with p", 300, 1.1, 1.7);
 437  fHistMassLpi->GetXaxis()->SetTitle("Invariant mass #Lambda#pi^{-} (GeV/c^{2})");
 438  fHistMassLpi->GetYaxis()->SetTitle("Entries");
 439  fHistMassLpi->SetMarkerStyle(kFullCircle);
 440
 441  fHistMassLP= new TH1F("fHistMassLP", "Invariant mass distribution of the #Lambda p used for combing with #pi^{-}", 300, 2.0, 2.3);
 442  fHistMassLP->GetXaxis()->SetTitle("Invariant mass #Lambda p (GeV/c^{2})");
 443  fHistMassLP->GetYaxis()->SetTitle("Entries");
 444  fHistMassLP->SetMarkerStyle(kFullCircle);
 445
 446  fHistProtonPIDBb = new TH2F("fHistProtonPIDBb", "dE/dx after p PID", 100, 0., 10, 100, 0, 100);
 447  fHistProtonPIDBb->GetYaxis()->SetTitle("TPC Signal");
 448  fHistProtonPIDBb->GetXaxis()->SetTitle("P (GeV/c)");
 449  fHistProtonPIDBb->SetOption("scat");
 450  fHistProtonPIDBb->SetMarkerStyle(kFullCircle);
 451
 452  fHistPionPIDBb = new TH2F("fHistPionPIDBb", "dE/dx after K PID", 100, 0., 10, 100, 0, 100);
 453  fHistPionPIDBb->GetYaxis()->SetTitle("TPC Signal");
 454  fHistPionPIDBb->GetXaxis()->SetTitle("P (GeV/c)");
 455  fHistPionPIDBb->SetOption("scat");
 456  fHistPionPIDBb->SetMarkerStyle(kFullCircle);
 457
 458  fHistProtonPIDLambda = new TH2F("fHistProtonPIDLambda", "dE/dx after p PID from V0", 100, 0., 10, 100, 0, 100);
 459  fHistProtonPIDLambda->GetYaxis()->SetTitle("TPC Signal");
 460  fHistProtonPIDLambda->GetXaxis()->SetTitle("P (GeV/c)");
 461  fHistProtonPIDLambda->SetOption("scat");
 462  fHistProtonPIDLambda->SetMarkerStyle(kFullCircle);
 463
 464  fHistPionPIDLambda = new TH2F("fHistPionPIDLambda", "dE/dx after #pi PID from V0", 100, 0, 10, 100, 0, 100);
 465  fHistPionPIDLambda->GetYaxis()->SetTitle("TPC Signal");
 466  fHistPionPIDLambda->GetXaxis()->SetTitle("P (GeV/c)");
 467  fHistPionPIDLambda->SetOption("scat");
 468  fHistPionPIDLambda->SetMarkerStyle(kFullCircle);
 469
 470  fHistMCdcaPvtxDvtx= new TH1F("fHistMCdcaPvtxDvtx", "MC True DCA Primary Vertex - H Decay Vertex", 300, -0.1, 11.9);
 471  fHistMCdcaPvtxDvtx->GetXaxis()->SetTitle("dca prim. vtx- decay vtx (cm)");
 472  fHistMCdcaPvtxDvtx->GetYaxis()->SetTitle("Entries");
 473  fHistMCdcaPvtxDvtx->SetMarkerStyle(kFullCircle);
 474
 475  fHistMCdcaPvtxLvtx= new TH1F("fHistMCdcaPvtxLvtx", "MC True DCA Primary Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
 476  fHistMCdcaPvtxLvtx->GetXaxis()->SetTitle("dca prim. vtx-#Lambda decay vtx (cm)");
 477  fHistMCdcaPvtxLvtx->GetYaxis()->SetTitle("Entries");
 478  fHistMCdcaPvtxLvtx->SetMarkerStyle(kFullCircle);
 479
 480  fHistMCdcaDvtxLvtx= new TH1F("fHistMCdcaDvtxLvtx", "MC True DCA H Decay Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
 481  fHistMCdcaDvtxLvtx->GetXaxis()->SetTitle("dca H deacy vtx-#Lambda decay vtx (cm)");
 482  fHistMCdcaDvtxLvtx->GetYaxis()->SetTitle("Entries");
 483  fHistMCdcaDvtxLvtx->SetMarkerStyle(kFullCircle);
 484
 485  fHistMCangleLH= new TH1F("fHistMCangleLH", "MC True Angle between #Lambda and H", 300, 0., 2*TMath::Pi());
 486  fHistMCangleLH->GetXaxis()->SetTitle("Angle (#Lambda-H)");
 487  fHistMCangleLH->GetYaxis()->SetTitle("Entries");
 488  fHistMCangleLH->SetMarkerStyle(kFullCircle);
 489
 490  fHistMCdecayAngle= new TH1F("fHistMCdecayAngle", "MC True Angle between decay products", 300, 0., 2*TMath::Pi());
 491  fHistMCdecayAngle->GetXaxis()->SetTitle("Angle (#Lambda-p#pi)");
 492  fHistMCdecayAngle->GetYaxis()->SetTitle("Entries");
 493  fHistMCdecayAngle->SetMarkerStyle(kFullCircle);
 494
 495  fHistMCpointingAngle= new TH1F("fHistMCpointingAngle", "MC True Pointing Angle", 3000, 0., 2*TMath::Pi());
 496  fHistMCpointingAngle->GetXaxis()->SetTitle("Pointing Angle");
 497  fHistMCpointingAngle->GetYaxis()->SetTitle("Entries");
 498  fHistMCpointingAngle->SetMarkerStyle(kFullCircle);
 499
 500  fHistMCap = new TH2F("fHistMCap", "True MC Armenteros-Podolanski", 200,-1.0,1.0, 500,0,1);
 501  fHistMCap->GetYaxis()->SetTitle("#alpha");
 502  fHistMCap->GetXaxis()->SetTitle("q_{t}");
 503  fHistMCap->SetOption("scat");
 504  fHistMCap->SetMarkerStyle(kFullCircle);
 505
 506  fHistMCdcaPvtxDvtxReso= new TH1F("fHistMCdcaPvtxDvtxReso", "MC True DCA Primary Vertex - H Decay Vertex", 300, -0.1, 11.9);
 507  fHistMCdcaPvtxDvtxReso->GetXaxis()->SetTitle("dca prim. vtx- decay vtx (cm)");
 508  fHistMCdcaPvtxDvtxReso->GetYaxis()->SetTitle("Entries");
 509  fHistMCdcaPvtxDvtxReso->SetMarkerStyle(kFullCircle);
 510
 511  fHistMCdcaPvtxLvtxReso= new TH1F("fHistMCdcaPvtxLvtxReso", "MC True DCA Primary Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
 512  fHistMCdcaPvtxLvtxReso->GetXaxis()->SetTitle("dca prim. vtx-#Lambda decay vtx (cm)");
 513  fHistMCdcaPvtxLvtxReso->GetYaxis()->SetTitle("Entries");
 514  fHistMCdcaPvtxLvtxReso->SetMarkerStyle(kFullCircle);
 515
 516  fHistMCdcaDvtxLvtxReso= new TH1F("fHistMCdcaDvtxLvtxReso", "MC True DCA H Decay Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
 517  fHistMCdcaDvtxLvtxReso->GetXaxis()->SetTitle("dca H deacy vtx-#Lambda decay vtx (cm)");
 518  fHistMCdcaDvtxLvtxReso->GetYaxis()->SetTitle("Entries");
 519  fHistMCdcaDvtxLvtxReso->SetMarkerStyle(kFullCircle);
 520
 521  fHistMCangleLHReso= new TH1F("fHistMCangleLHReso", "MC True Angle between #Lambda and H", 300, 0., 2*TMath::Pi());
 522  fHistMCangleLHReso->GetXaxis()->SetTitle("Angle (#Lambda-H)");
 523  fHistMCangleLHReso->GetYaxis()->SetTitle("Entries");
 524  fHistMCangleLHReso->SetMarkerStyle(kFullCircle);
 525
 526  fHistMCdecayAngleReso= new TH1F("fHistMCdecayAngleReso", "MC True Angle between decay products", 300, 0., 2*TMath::Pi());
 527  fHistMCdecayAngleReso->GetXaxis()->SetTitle("Angle (#Lambda-p#pi)");
 528  fHistMCdecayAngleReso->GetYaxis()->SetTitle("Entries");
 529  fHistMCdecayAngleReso->SetMarkerStyle(kFullCircle);
 530
 531  fHistMCpointingAngleReso= new TH1F("fHistMCpointingAngleReso", "MC True Pointing Angle", 300, 0., 2*TMath::Pi());
 532  fHistMCpointingAngleReso->GetXaxis()->SetTitle("Pointing Angle");
 533  fHistMCpointingAngleReso->GetYaxis()->SetTitle("Entries");
 534  fHistMCpointingAngleReso->SetMarkerStyle(kFullCircle);
 535
 536  fHistMCapReso = new TH2F("fHistMCapReso", "True MC Armenteros-Podolanski", 200,-1.0,1.0, 500,0,1);
 537  fHistMCapReso->GetYaxis()->SetTitle("#alpha");
 538  fHistMCapReso->GetXaxis()->SetTitle("q_{t}");
 539  fHistMCapReso->SetOption("scat");
 540  fHistMCapReso->SetMarkerStyle(kFullCircle);
 541
 542  fHistCentrality = new TH1F("Centrality ", "Centrality", 11, -0.5, 10.5);
 543  fHistCentrality ->GetXaxis()->SetTitle("Centrality");
 544  fHistCentrality ->GetYaxis()->SetTitle("Entries");
 545
 546  fHistCentralityAC = new TH1F("CentralityAC ", "CentralityAC", 11, -0.5, 10.5);
 547  fHistCentralityAC ->GetXaxis()->SetTitle("Centrality");
 548  fHistCentralityAC ->GetYaxis()->SetTitle("Entries");
 549
 550  fHistMultiplicity = new TH1F("Multiplicity ", "Multiplicity", 100, 0, 20000);
 551  fHistMultiplicity ->GetXaxis()->SetTitle("Centrality");
 552  fHistMultiplicity ->GetYaxis()->SetTitle("Entries");
 553
 554  fHistList->Add(fHistMassDPi);
 555  fHistList->Add(fHistMassLPi);
 556  fHistList->Add(fHistMassLambdaPi);
 557  fHistList->Add(fHistMassLambda);
 558  fHistList->Add(fHistMassLambdaPPi);
 559  fHistList->Add(fHistMassLambdaP);
 560  fHistList->Add(fHistMassLambdaK);
 561  fHistList->Add(fHistMassK0onl);
 562  fHistList->Add(fHistMassK0offl);
 563  fHistList->Add(fHistMassK0onlC);
 564  fHistList->Add(fHistMassK0offlC);
 565  fHistList->Add(fHistMassPQonl);
 566  fHistList->Add(fHistMassPQoffl);
 567  fHistList->Add(fHistDC);
 568  fHistList->Add(fHistArmenterosPodolanski);
 569  fHistList->Add(fHistArmenterosPodolanskiCut);
 570  fHistList->Add(fHistHDibaryonInvaMassGen);
 571  fHistList->Add(fHistHDibaryonInvaMassGenRes);
 572  fHistList->Add(fHistAntiHDibaryonInvaMassGen);
 573  fHistList->Add(fHistHDibaryonInvaMassAso);
 574  fHistList->Add(fHistHDibaryonInvaMassAsoReso);
 575  fHistList->Add(fHistAntiHDibaryonInvaMassAso);
 576  fHistList->Add(fHistCheck);
 577  fHistList->Add(fHistHPointingAngle);
 578  fHistList->Add(fHistMassH);
 579  fHistList->Add(fHistMassPpi);
 580  fHistList->Add(fHistMassPpiReso);
 581  fHistList->Add(fHistMassLpi);
 582  fHistList->Add(fHistMassLP);
 583  fHistList->Add(fHistMassLambdaFromH);
 584  fHistList->Add(fHistMassLambdaFromHtLorentz);
 585  fHistList->Add(fHistProtonPIDBb);
 586  fHistList->Add(fHistPionPIDBb);
 587  fHistList->Add(fHistProtonPIDLambda);
 588  fHistList->Add(fHistPionPIDLambda);
 589  fHistList->Add(fHistMCdcaPvtxDvtx);
 590  fHistList->Add(fHistMCdcaPvtxLvtx);
 591  fHistList->Add(fHistMCdcaDvtxLvtx);
 592  fHistList->Add(fHistMCangleLH);
 593  fHistList->Add(fHistMCdecayAngle);
 594  fHistList->Add(fHistMCpointingAngle);
 595  fHistList->Add(fHistMCap);
 596  fHistList->Add(fHistMCdcaPvtxDvtxReso);
 597  fHistList->Add(fHistMCdcaPvtxLvtxReso);
 598  fHistList->Add(fHistMCdcaDvtxLvtxReso);
 599  fHistList->Add(fHistMCangleLHReso);
 600  fHistList->Add(fHistMCdecayAngleReso);
 601  fHistList->Add(fHistMCpointingAngleReso);
 602  fHistList->Add(fHistMCapReso);
 603  fHistList->Add(fHistCentrality);
 604  fHistList->Add(fHistCentralityAC);
 605  fHistList->Add(fHistMultiplicity);
 606
 607  fHistHilf1= new TH1F("fHistHilf1", "HD", 300, 0., 10);
 608  fHistHilf1->GetXaxis()->SetTitle("");
 609  fHistHilf1->GetYaxis()->SetTitle("Entries");
 610
 611  fHistHilf2= new TH1F("fHistHilf2", "HD", 300, 0., 10);
 612  fHistHilf2->GetXaxis()->SetTitle("");
 613  fHistHilf2->GetYaxis()->SetTitle("Entries");
 614
 615  fHistHilf3= new TH1F("fHistHilf3", "HD", 300, 0., 10);
 616  fHistHilf3->GetXaxis()->SetTitle("");
 617  fHistHilf3->GetYaxis()->SetTitle("Entries");
 618
 619  fHistHilf4= new TH1F("fHistHilf4", "HD", 300, 0., 10);
 620  fHistHilf4->GetXaxis()->SetTitle("");
 621  fHistHilf4->GetYaxis()->SetTitle("Entries");
 622
 623  fHistHilf5= new TH1F("fHistHilf5", "HD", 300, 0., 10);
 624  fHistHilf5->GetXaxis()->SetTitle("");
 625  fHistHilf5->GetYaxis()->SetTitle("Entries");
 626
 627  fHistHilf6= new TH1F("fHistHilf6", "HD", 300, 0., 10);
 628  fHistHilf6->GetXaxis()->SetTitle("");
 629  fHistHilf6->GetYaxis()->SetTitle("Entries");
 630
 631  fHistPtvsEtaGen = new TH2F("fHistPtvsEtaGen", "p_{t} vs #eta from generated H", 200,0.0,10.0, 200,-1,1);
 632  fHistPtvsEtaGen->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 633  fHistPtvsEtaGen->GetYaxis()->SetTitle("#eta");
 634  fHistPtvsEtaGen->SetOption("scat");
 635  fHistPtvsEtaGen->SetMarkerStyle(kFullCircle);
 636
 637  fHistPtvsEtaAso = new TH2F("fHistPtvsEtaAso", "p_{t} vs #eta from associated H", 200,0.0,10.0, 200,-1,1);
 638  fHistPtvsEtaAso->GetYaxis()->SetTitle("p_{t} (GeV/c)");
 639  fHistPtvsEtaAso->GetXaxis()->SetTitle("#eta");
 640  fHistPtvsEtaAso->SetOption("scat");
 641  fHistPtvsEtaAso->SetMarkerStyle(kFullCircle);
 642
 643  fHistPtvsYGen = new TH2F("fHistPtvsYGen", "p_{t} vs rapidity from generated H", 200,0.0,10.0, 200,-1,1);
 644  fHistPtvsYGen->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 645  fHistPtvsYGen->GetYaxis()->SetTitle("y");
 646  fHistPtvsYGen->SetOption("scat");
 647  fHistPtvsYGen->SetMarkerStyle(kFullCircle);
 648
 649  fHistPtvsYAso = new TH2F("fHistPtvsYAso", "p_{t} vs rapidity from associated H", 200,0.0,10.0, 200,-1,1);
 650  fHistPtvsYAso->GetXaxis()->SetTitle("p_{t} (GeV/c)");
 651  fHistPtvsYAso->GetYaxis()->SetTitle("y");
 652  fHistPtvsYAso->SetOption("scat");
 653  fHistPtvsYAso->SetMarkerStyle(kFullCircle);
 654
 655  fHistRap= new TH1F("fHistRap", "Rapidity", 400, -2., 2);
 656  fHistRap->GetXaxis()->SetTitle("Y");
 657  fHistRap->GetYaxis()->SetTitle("Entries");
 658  fHistPtvsEtaAso->SetMarkerStyle(kFullCircle);
 659
 660  fHistList->Add(fHistHilf1);
 661  fHistList->Add(fHistHilf2);
 662  fHistList->Add(fHistHilf3);
 663  fHistList->Add(fHistHilf4);
 664  fHistList->Add(fHistHilf5);
 665  fHistList->Add(fHistHilf6);
 666  fHistList->Add(fHistPtvsEtaGen);
 667  fHistList->Add(fHistPtvsEtaAso);
 668  fHistList->Add(fHistPtvsYGen);
 669  fHistList->Add(fHistPtvsYAso);
 670  fHistList->Add(fHistRap);
 671
 672  fHistCount = new TH1F("fHistCount","test",17,0,17);
 673  fHistCount->GetXaxis()->SetBinLabel(1,"Events");
 674  fHistCount->GetXaxis()->SetBinLabel(2,"MC All");
 675  fHistCount->GetXaxis()->SetBinLabel(3,"MC from Primary Vtx");
 676  fHistCount->GetXaxis()->SetBinLabel(4,"Horst");
 677  fHistCount->GetXaxis()->SetBinLabel(5,"Lambda Candidates");
 678  fHistCount->GetXaxis()->SetBinLabel(6,"Sigma Candidates");
 679  fHistCount->GetXaxis()->SetBinLabel(7,"Horst");
 680  fHistCount->GetXaxis()->SetBinLabel(8,"Horst");
 681  fHistCount->GetXaxis()->SetBinLabel(9,"Horst");
 682  fHistCount->GetXaxis()->SetBinLabel(10,"MC All #bar{Lambda}(1520)s");
 683  fHistCount->GetXaxis()->SetBinLabel(11,"");
 684  fHistCount->GetXaxis()->SetBinLabel(12,"H-Dibaryon");
 685  fHistCount->GetXaxis()->SetBinLabel(13,"Hypertriton 2-Body");
 686  fHistCount->GetXaxis()->SetBinLabel(14,"Hypertriton 3-Body");
 687  fHistCount->GetXaxis()->SetBinLabel(15,"");
 688  fHistCount->GetXaxis()->SetBinLabel(16,"");
 689  fHistCount->GetXaxis()->SetBinLabel(17,"Lambdas!!!");
 690  fHistCount->SetStats(0);
 691  fHistCount->SetFillColor(38);
 692  fHistList->Add(fHistCount);
 693
 694  //trigger statistics histogram
 695   fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", 4,-0.5, 3.5);
 696   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral" };
 697   for ( Int_t ii=0; ii < 3; ii++ )
 698     fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 699
 700   fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", 4,-0.5, 3.5);
 701   for ( Int_t ii=0; ii < 3; ii++ )
 702     fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 703   fHistList->Add(fHistTriggerStat);
 704   fHistList->Add(fHistTriggerStatAfterEventSelection);
 705
 706   fHistMassHcentMult =  new TH3F("fHistMassHcentMult", "Inv. Mass vs. centrality vs. multiplicity", 100, 2.2, 2.3, 5, 0, 4, 300, 0, 6000);
 707   fHistMassHcentMult->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})"); // inv. mass
 708   fHistMassHcentMult->GetYaxis()->SetTitle("Centrality"); // triggertype
 709   fHistMassHcentMult->GetZaxis()->SetTitle("Multiplicity"); // refTPC
 710   fHistList->Add(fHistMassHcentMult);
 711
 712
 713   const Double_t kz = 2*TMath::Pi();
 714   Int_t binsD01[16]={  300, 200, 100, 100, 100, 100, 100, 100, 200, 200, 200, 200, 400, 200, 200, 3};
 715   Double_t xminD01[16]={2.0, 1.0, 0., -1, 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.,  0., 0};
 716   Double_t xmaxD01[16]={2.3, 1.2, kz, 1, 1, 10, 10, 5, 5, 5, 5, 100, 1, 100,  4000, 1};
 717
 718
 719   fHistNdim = new THnSparseF("fHistNdim","THnS;InvMass, InvMassLambda, pointingAngle, armPoAlpha, armPoQt, pTL, pTH, d0p, d0n, dcaHd, dca, decayL, cosPA, centr, multi, mcinf;InvMassH", 16,binsD01,xminD01,xmaxD01);
 720   fHistList->Add(fHistNdim);
 721 }
 722
 723  //________________________________________________________________________
 724 void AliAnalysisTaskHdibaryonLPpi::UserExec(Option_t *)
 725 {
 726   // Main loop
 727   // Called for each event
 728
 729   //define improtant masses
 730   Double_t pionMass     = 0.13957;
 731   Double_t protonMass   = 0.93827;
 732
 733   //define PDGCodes
 734   Long_t pdgPionPlus          =         211;
 735   Long_t pdgPionMinus         =        -211;
 736   Long_t pdgProton            =        2212;
 737   Long_t pdgAntiProton        =       -2212;
 738   Long_t pdgLambda            =        3122;
 739   Long_t pdgAntiLambda        =       -3122;
 740   Long_t pdgHDibaryon         =  1020000020;
 741   Long_t pdgAntiHDibaryon     = -1020000020;
 742
 743   AliStack* stack(NULL);
 744   if(HasMC()){
 745
 746     if(!fMCEvent)return;
 747
 748     AliHeader *head = fMCEvent->Header();
 749     if(!head)return;
 750     AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)head->GenEventHeader();
 751     if(!pyheader)return;
 752
 753     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 754       if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()){
 755         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->TreeK()))return;
 756         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->TreeTR()))return;
 757         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->Init("local")))return;
 758         stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 759         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack()->TreeK()))return;
 760       }
 761     }
 762
 763     if(!stack)return;
 764   }
 765
 766   // -------------------------------------------------------
 767   // Loop for Inv. Mass via ESD tracks
 768   // -------------------------------------------------------
 769
 770   fHistCount->Fill(0);
 771
 772   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 773
 774   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 775   if (vertex->GetNContributors()<1)
 776     {
 777       // SPD vertex
 778       vertex = fESD->GetPrimaryVertexSPD();
 779        if(vertex->GetNContributors()<1) {
 780          return;
 781        }
 782     }
 783   if (TMath::Abs(vertex->GetZv()) > 10) return;
 784
 785   Int_t centrality = -5;
 786   Double_t centrPerc = -5;
 787
 788   if (fESD->GetEventSpecie() == 4)
 789     { // PbPb
 790       AliCentrality *esdCentrality = fESD->GetCentrality();
 791       centrality = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 792       centrPerc = esdCentrality->GetCentralityPercentile("V0M");
 793       if (centrality < 0. || centrality > 8.) return; //0 bis 80 %
 794       //  cout<<"Centrality: "<< centrality << endl;
 795     }
 796
 797   fHistCentrality->Fill(centrality);
 798
 799   //*****************//
 800   //*   Centrality  *//
 801   //*****************//
 802
 803   //  Float_t percentile=centrality->GetCentralityPercentile("V0M");
 804
 805   Bool_t isSelectedCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
 806   Bool_t isSelectedSemiCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
 807   Bool_t isSelectedMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 808
 809   Int_t triggertype = 17;
 810
 811   if(isSelectedCentral){
 812     fHistTriggerStat->Fill(1);
 813     triggertype=1;
 814   }
 815
 816   if(isSelectedSemiCentral){
 817     fHistTriggerStat->Fill(2);
 818     triggertype=2;
 819   }
 820
 821   if(isSelectedMB){
 822     fHistTriggerStat->Fill(0);
 823     triggertype=3;
 824   }
 825
 826   //  if(isSelectedCentral || isSelectedSemiCentral || isSelectedMB){
 827
 828   //*******************************
 829
 830   Int_t runNumber = 0;
 831   //  itrk = 0;
 832   runNumber = fESD->GetRunNumber();
 833
 834   if (!fPIDtpcESD) fPIDtpcESD = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 835   if (!fPIDtpcESD) {
 836     fPIDtpcESD = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
 837     fPIDtpcESD->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
 838   }
 839
 840   Double_t pionK=1;
 841   Double_t pK=1;
 842
 843   TObjArray* listCrossV0   = fESDCutsV0->GetAcceptedV0s(fESD);
 844   Int_t nGoodV0s      = listCrossV0->GetEntries();
 845
 846   const AliESDVertex *esdVer = fESD->GetPrimaryVertex();
 847   AliESDVertex *esdVer1 = new AliESDVertex(*esdVer);
 848
 849   AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
 850   TObjArray *trkArray = new TObjArray(2);
 851   AliVertexerTracks *vertexer1 = new AliVertexerTracks(fESD->GetMagneticField());
 852   TObjArray *trkArray1 = new TObjArray(2);
 853
 854   AliKFParticle::SetField(fESD->GetMagneticField());
 855
 856   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
 857
 858   Int_t refMultTpc = AliESDtrackCuts::GetReferenceMultiplicity(fESD, kTRUE);
 859   //cout<<"Multiplicity: "<< refMultTpc << endl;
 860   fHistMultiplicity->Fill(refMultTpc);
 861
 862   Double_t mn[3] = {0,0,0};
 863   Double_t mp[3] = {0,0,0};
 864   Double_t mm[3] = {0,0,0};
 865   Double_t dd[3] = {0,0,0};
 866   Double_t dd1[3] = {0,0,0};
 867   const Double_t cProtonMass=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
 868   const Double_t cPionMass=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 869   const Double_t cElectronMass=TDatabasePDG::Instance()->GetParticle(11)->Mass();
 870   const Double_t cLambdaMass=TDatabasePDG::Instance()->GetParticle(3122)->Mass();
 871   Double_t decayLength=0;
 872   Double_t decayLengthH=0;
 873
 874   //V0 Loop for Lambda and Anti-Lambda
 875   for(Int_t iV0MI = 0; iV0MI < nGoodV0s ; iV0MI++) {
 876     AliESDv0 * fV0MIs = fESD->GetV0(iV0MI);
 877     Int_t    lOnFlyStatus = 0;
 878
 879     lOnFlyStatus = fV0MIs->GetOnFlyStatus();
 880     Double_t lInvMassLambda=0;
 881     Double_t lInvMassLambdaPi=0;
 882     Double_t lPtLambda=0;
 883     Double_t lPzLambda=0;
 884     Double_t lPLambda=0;
 885     Int_t onl=0;
 886     Int_t offl=0;
 887
 888     TLorentzVector posE;
 889     TLorentzVector negE;
 890     TLorentzVector photon;
 891
 892     if (lOnFlyStatus){
 893       onl=1;
 894       //      return;
 895     }
 896     if (!lOnFlyStatus){
 897       offl=1;
 898       //return;
 899     }
 900
 901     //    fHistMultiplicity->Fill(refMultTpc);
 902     fHistCentralityAC->Fill(centrality);
 903     fHistCheck->Fill(offl,onl);
 904
 905     AliESDtrack* trackPosTest = fESD->GetTrack(fV0MIs->GetPindex());
 906     AliESDtrack* trackNegTest = fESD->GetTrack(fV0MIs->GetNindex());
 907
 908     //    if (!
 909     if (!fEsdTrackCuts->AcceptTrack(trackPosTest)) continue;
 910     if (!fESDtrackCutsV0->AcceptTrack(trackNegTest)) continue;
 911
 912
 913       //PID via specific energy loss in the TPC
 914       //define the arrays for the Bethe-Bloch-Parameters
 915       Double_t parProton[5]   = {0,0,0,0,0};
 916
 917       if(runNumber < 166500) //LHC10h
 918         {
 919           parProton[0]   = 1.45802; // ALEPH parameters for protons (pass2)
 920           parProton[1]   = 27.4992;
 921           parProton[2]   = 4.00313e-15;
 922           parProton[3]   = 2.48485;
 923           parProton[4]   = 8.31768;
 924         }
 925
 926       if(runNumber > 166500) //LHC11h
 927         {
 928           parProton[0]   = 1.11243; // ALEPH parameters for protons (pass2)
 929           parProton[1]   = 26.1144;
 930           parProton[2]   = 4.00313e-15;
 931           parProton[3]   = 2.72969;
 932           parProton[4]   = 9.15038;
 933         }
 934
 935       //Get the total momentum for each track at the inner readout of the TPC
 936       Double_t ptotN = trackNegTest->GetInnerParam()->GetP();
 937       Double_t ptotP = trackPosTest->GetInnerParam()->GetP();
 938
 939       //define expected signals for the various species
 940       Double_t expSignalPionP = 0;
 941       Double_t expSignalPionN = 0;
 942       Double_t expSignalProtonN = 0;
 943       Double_t expSignalProtonP = 0;
 944
 945       //for data
 946       if(!HasMC())
 947         {
 948           expSignalProtonN = AliExternalTrackParam::BetheBlochAleph(ptotN/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 949           expSignalProtonP = AliExternalTrackParam::BetheBlochAleph(ptotP/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 950         }
 951
 952       //for MC
 953       if(HasMC())
 954         {
 955           expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(pionMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 956           expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(pionMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 957
 958           expSignalProtonN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 959           expSignalProtonP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
 960         }
 961
 962       // PID cut on the nuclei (proton, deuteron, triton, helium3)
 963       Bool_t corrParticle = kFALSE;
 964
 965       Bool_t posProton = kFALSE;
 966       //data
 967       if(!HasMC())
 968         {
 969           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackPosTest, AliPID::kProton)) < 3)
 970             {
 971               posProton = kTRUE;
 972               corrParticle = kTRUE;
 973             }
 974         }
 975       //MC
 976       if(HasMC())
 977         {
 978           if(//trackPosTest->GetTPCsignal() < 1200 &&
 979              TMath::Abs(trackPosTest->GetTPCsignal() - expSignalProtonP)/expSignalProtonP < 0.4)
 980             {
 981               posProton = kTRUE;
 982               corrParticle = kTRUE;
 983             }
 984         }
 985
 986       Bool_t negProton = kFALSE;
 987       //data
 988       if(!HasMC())
 989         {
 990           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackNegTest, AliPID::kProton)) < 3)
 991             {
 992               negProton = kTRUE;
 993               corrParticle = kTRUE;
 994             }
 995         }
 996      //MC
 997       if(HasMC())
 998         {
 999           if(//trackNegTest->GetTPCsignal() < 1200 &&
1000              TMath::Abs(trackNegTest->GetTPCsignal() - expSignalProtonN)/expSignalProtonN < 0.4)
1001             {
1002               negProton = kTRUE;
1003               corrParticle = kTRUE;
1004             }
1005         }
1006
1007      //PID cut for pions
1008       //data: 3sigma cut on the pions
1009
1010       Bool_t negPion = kFALSE;
1011       Bool_t posPion = kFALSE;
1012
1013       if (!HasMC())
1014         {
1015           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackPosTest, AliPID::kPion)) < 4) posPion=kTRUE; //pos daughter has to be a pion
1016           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackNegTest, AliPID::kPion)) < 4) negPion=kTRUE; // negative daughter has to be a pion
1017         }
1018
1019       //MC: like the nuclei via the specific energyloss in the TPC
1020       if (HasMC())
1021         {
1022           if (TMath::Abs(trackPosTest->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.4) posPion=kTRUE;
1023           if (TMath::Abs(trackNegTest->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.4) negPion=kTRUE;
1024         }
1025
1026       if (!(posProton==kTRUE)) continue;
1027       if (!(negPion==kTRUE)) continue;
1028
1029     //To avoid ghosts
1030
1031     if( !(trackPosTest->GetStatus() & AliESDtrack::kTPCrefit)){
1032       continue;
1033     }
1034
1035     if( !(trackNegTest->GetStatus() & AliESDtrack::kTPCrefit)){
1036       continue;
1037     }
1038
1039     if( trackPosTest->GetSign() >0 && trackNegTest->GetSign() <0){
1040       fV0MIs->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1041       fV0MIs->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
1042     }
1043
1044     if( trackPosTest->GetSign() <0 && trackNegTest->GetSign() >0){
1045       fV0MIs->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1046       fV0MIs->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
1047     }
1048
1049     fV0MIs->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
1050
1051     TVector3 vecN(mn[0],mn[1],mn[2]);
1052     TVector3 vecP(mp[0],mp[1],mp[2]);
1053     TVector3 vecM(mm[0],mm[1],mm[2]);
1054
1055     Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
1056     Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
1057
1058     Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
1059       ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
1060     Double_t qt = vecP.Mag()*sin(thetaP);
1061
1062     fHistArmenterosPodolanski->Fill(alfa,qt);    //Armenteros-Podolanski calculation
1063
1064     TLorentzVector k0;
1065     TLorentzVector k0daugh1;
1066     TLorentzVector k0daugh2;
1067     TLorentzVector proton;
1068     TLorentzVector pq;
1069
1070     k0daugh1.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1071     k0daugh2.SetXYZM(mp[0],mp[1],mp[2],cPionMass);
1072     k0=k0daugh1+k0daugh2;
1073
1074     fV0MIs->ChangeMassHypothesis(3122);
1075     lInvMassLambda = fV0MIs->GetEffMass();
1076     lPtLambda = fV0MIs->Pt();
1077     lPzLambda = fV0MIs->Pz();
1078     lPLambda = fV0MIs->P();
1079
1080     trkArray->AddAt(trackPosTest,0);
1081     trkArray->AddAt(trackNegTest,1);
1082
1083     vertexer->SetVtxStart(esdVer1);
1084     AliESDVertex *decayVertex = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1085
1086     dd[0]=fESD->GetPrimaryVertexSPD()->GetX()-decayVertex->GetX();
1087     dd[1]=fESD->GetPrimaryVertexSPD()->GetY()-decayVertex->GetY();
1088     dd[2]=fESD->GetPrimaryVertexSPD()->GetZ()-decayVertex->GetZ();
1089
1090     decayLength=sqrt(dd[0]*dd[0]+dd[1]*dd[1]+dd[2]*dd[2]);
1091
1092     if (decayVertex) delete decayVertex;
1093
1094     TLorentzVector negPio1;
1095     TLorentzVector posProt1;
1096     TLorentzVector posP;
1097     TLorentzVector posProt;
1098     TLorentzVector negK;
1099     TLorentzVector negPio;
1100     TLorentzVector negPi;
1101     TLorentzVector omega;
1102     TLorentzVector threeSum;
1103     TLorentzVector fourSum;
1104     TLorentzVector ppK;
1105     TLorentzVector posPiK;
1106     TLorentzVector negPiK;
1107     TLorentzVector kaon;
1108     TLorentzVector lambda;
1109     TLorentzVector lambdaH;
1110     TLorentzVector hDibaryon;
1111     TVector3 h;
1112     TVector3 h1;
1113
1114     Int_t mcStatus=0;
1115
1116     h.SetXYZ(-dd[0],-dd[1],-dd[2]);
1117
1118     if (onl==1)fHistMassDPi->Fill(lInvMassLambda);
1119
1120     if (offl==1||onl==1){
1121       fHistMassLPi->Fill(lInvMassLambda);
1122
1123       negE.SetXYZM(mn[0],mn[1],mn[2],cElectronMass);
1124       posE.SetXYZM(mp[0],mp[1],mp[2],cElectronMass);
1125       photon=posE+negE;
1126
1127       negPiK.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1128       posPiK.SetXYZM(mp[0],mp[1],mp[2],cPionMass);
1129       kaon=posPiK+negPiK;
1130
1131       negPi.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1132       posP.SetXYZM(mp[0],mp[1],mp[2],cProtonMass);
1133       lambda=negPi+posP;
1134       lambdaH.SetXYZM(mm[0],mm[1],mm[2],cLambdaMass);
1135
1136       if (lInvMassLambda>1.1113&&lInvMassLambda<1.1202){
1137
1138         if (!HasMC()){
1139           if (qt<-2.21*alfa*alfa+2.945*alfa-0.887) continue;
1140           if (qt>-2.21*alfa*alfa+2.945*alfa-0.873) continue;
1141           if (photon.M()<0.005) continue;
1142           if (kaon.M()>0.495 && kaon.M()<0.500 ) continue;
1143         }
1144
1145
1146         fHistMassLambda->Fill(lInvMassLambda);
1147         //
1148         Bool_t isCorrectlyAssociatedLambda = kFALSE;
1149         Bool_t isPartialCorrectlyAssociatedLambda = kFALSE;
1150         Int_t labelAssociatedH=-1;
1151         Int_t labelLambda=-1;
1152         //
1153         if (HasMC()) {
1154           Int_t labelPosTest = trackPosTest->GetLabel();
1155           TParticle *tparticleDaughter = stack->Particle(TMath::Abs(labelPosTest));
1156           Int_t labelMother = tparticleDaughter->GetFirstMother();
1157           TParticle *tparticleMother = stack->Particle(TMath::Abs(labelMother));
1158
1159           Int_t labelOma = tparticleMother->GetFirstMother();
1160           TParticle *tparticleOma = stack->Particle(TMath::Abs(labelOma));
1161
1162           if ((tparticleOma->GetPdgCode() < 0) && TMath::Abs(tparticleMother->GetPdgCode())==pdgLambda){// check mother to be Lambda
1163             Int_t labelFirstDaughter  = tparticleMother->GetDaughter(1);// Proton
1164             Int_t labelThirdDaughter  = tparticleMother->GetDaughter(0);// Pion
1165
1166             TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1167             TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1168
1169             if((tparticleFirstDaughter->GetPdgCode() == pdgProton && tparticleThirdDaughter->GetPdgCode()== pdgPionMinus) ||
1170                (tparticleFirstDaughter->GetPdgCode() == pdgPionMinus && tparticleThirdDaughter->GetPdgCode()== pdgProton)){ //daughter PDGs
1171               isPartialCorrectlyAssociatedLambda = kTRUE;
1172               labelLambda=labelMother;
1173             }
1174           }
1175
1176           //H-Dibaryon
1177           if(tparticleOma->GetPdgCode() == pdgHDibaryon){ //check grandmother to be H PDG
1178             if (TMath::Abs(tparticleMother->GetPdgCode())==pdgLambda){// check mother to be Lambda
1179               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(1);// Proton
1180               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(0);// Pion
1181
1182               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1183               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1184
1185               if((tparticleFirstDaughter->GetPdgCode() == pdgProton && tparticleThirdDaughter->GetPdgCode()== pdgPionMinus) ||
1186                  (tparticleFirstDaughter->GetPdgCode() == pdgPionMinus && tparticleThirdDaughter->GetPdgCode()== pdgProton)){ //daughter PDGs
1187                 isCorrectlyAssociatedLambda = kTRUE;
1188                 labelAssociatedH=labelOma;
1189                 fHistMassLambdaFromH->Fill(lInvMassLambda);
1190                 fHistMassLambdaFromHtLorentz->Fill(lambda.M());
1191               }
1192             }
1193           }
1194         }
1195
1196         fHistProtonPIDLambda->Fill(trackPosTest->GetInnerParam()->GetP(), trackPosTest->GetTPCsignal());
1197         fHistPionPIDLambda->Fill(trackNegTest->GetInnerParam()->GetP(), trackNegTest->GetTPCsignal());
1198
1199         //---------------------------------------------------------
1200         // Proton track loop
1201         //---------------------------------------------------------
1202         fHistArmenterosPodolanskiCut->Fill(alfa,qt);
1203
1204         for (Int_t iTracksP = 0; iTracksP < fESD->GetNumberOfTracks(); iTracksP++) {
1205           AliESDtrack* trackP = dynamic_cast<AliESDtrack*> (fESD->GetTrack(iTracksP));
1206           if (trackP->GetSign()<0) continue;
1207
1208           if (iTracksP==fV0MIs->GetPindex())continue;
1209           if (iTracksP==fV0MIs->GetNindex())continue;
1210
1211           if (!fEsdTrackCuts->AcceptTrack(trackP)) continue;
1212
1213           AliKFParticle protonKF( *(trackP), 2212);
1214
1215           if (!trackP->GetInnerParam()) continue;
1216
1217           if (HasMC()) {
1218             pK=0.65;
1219           }
1220
1221           if (!HasMC()){
1222             if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackP, AliPID::kProton)) > 3) continue;
1223           }
1224
1225           fHistProtonPIDBb->Fill(trackP->GetInnerParam()->GetP(), trackP->GetTPCsignal());
1226
1227           posProt.SetXYZM(trackP->Px(),trackP->Py(),trackP->Pz(),cProtonMass);
1228
1229           //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1230           //Pion Track loop!!!!!!!!!!!!!!!!!!!!!
1231           //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1232           for (Int_t iTracksN = iTracksP+1; iTracksN < fESD->GetNumberOfTracks(); iTracksN++) {
1233             AliESDtrack* trackN = dynamic_cast<AliESDtrack*> (fESD->GetTrack(iTracksN));
1234
1235             if (iTracksN==fV0MIs->GetPindex())continue;
1236             if (iTracksN==fV0MIs->GetNindex())continue;
1237             if (trackN->GetSign()>0) continue;
1238
1239             if (!fEsdTrackCuts->AcceptTrack(trackN)) continue;
1240             if (!fESDtrackCutsV0->AcceptTrack(trackN)) continue;
1241
1242             negPi.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1243             posP.SetXYZM(mp[0],mp[1],mp[2],cProtonMass);
1244             negPio.SetXYZM(trackN->Px(),trackN->Py(),trackN->Pz(),cPionMass);
1245
1246             threeSum=negPi+posP+negPio;
1247             lInvMassLambdaPi=threeSum.M();
1248
1249             fHistDC->Fill(decayLength*lInvMassLambda/lPLambda);
1250
1251             AliKFParticle posPionKF( *(trackN) ,-211);
1252
1253             if (!trackN->GetInnerParam()) continue;
1254             if (HasMC()) {
1255               pionK=0.7;
1256             }
1257
1258             if (!HasMC()){
1259               if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackN, AliPID::kPion)) > 3) continue;
1260             }
1261             fHistPionPIDBb->Fill(trackN->GetInnerParam()->GetP(), trackN->GetTPCsignal());
1262
1263             trkArray1->AddAt(trackP,0);
1264             trkArray1->AddAt(trackN,1);
1265
1266             vertexer1->SetVtxStart(esdVer1);
1267             AliESDVertex *decayVertex1 = (AliESDVertex*)vertexer1->VertexForSelectedESDTracks(trkArray1);
1268
1269             dd1[0]=fESD->GetPrimaryVertexSPD()->GetX()-decayVertex1->GetX();
1270             dd1[1]=fESD->GetPrimaryVertexSPD()->GetY()-decayVertex1->GetY();
1271             dd1[2]=fESD->GetPrimaryVertexSPD()->GetZ()-decayVertex1->GetZ();
1272
1273             decayLengthH=sqrt(dd1[0]*dd1[0]+dd1[1]*dd1[1]+dd1[2]*dd1[2]);
1274
1275             Double_t bz = fESD->GetMagneticField();
1276
1277             trackP->PropagateToDCA(decayVertex1, bz, 10);
1278             trackN->PropagateToDCA(decayVertex1, bz, 10);
1279
1280             Double_t xthiss(0.0);
1281             Double_t xpp(0.0);
1282             Double_t dca = trackN->GetDCA(trackP,bz,xthiss,xpp);
1283
1284             if (decayVertex1) delete decayVertex1;
1285             h1.SetXYZ(-dd1[0],-dd1[1],-dd1[2]);
1286
1287             //      if (dca>1) continue;
1288             if (dca>0.1) continue;
1289
1290             fourSum=threeSum+posProt;
1291
1292             posProt1.SetXYZM(trackP->Px(),trackP->Py(),trackP->Pz(),cProtonMass);
1293             negPio1.SetXYZM(trackN->Px(),trackN->Py(),trackN->Pz(),cPionMass);
1294             hDibaryon=lambdaH+posProt1+negPio1;
1295             Double_t hPointingAngle = hDibaryon.Angle(h);
1296             Double_t pointingAngleH = hDibaryon.Angle(h1);
1297             Double_t decayAngleH = h.Angle(h1);
1298             TVector3 vecDist(dd[0]-dd1[0],dd[1]-dd1[1],dd[2]-dd1[2]);
1299             fHistMassLambdaPPi->Fill(hDibaryon.M());
1300             fHistHPointingAngle->Fill(pointingAngleH);
1301
1302             fHistMassHcentMult->Fill(hDibaryon.M(),triggertype,refMultTpc);
1303
1304             Double_t rapidity = hDibaryon.Rapidity();
1305             if(rapidity > 1.0 || rapidity < -1.0) continue;
1306
1307             Double_t vec[16]={hDibaryon.M(), lInvMassLambda, pointingAngleH, alfa, qt, lPtLambda, hDibaryon.Pt(), posPionKF.GetDistanceFromVertex(primVtx), protonKF.GetDistanceFromVertex(primVtx), dca, protonKF.GetDistanceFromVertex(posPionKF), TMath::Cos(pointingAngleH), centrPerc, refMultTpc, mcStatus};
1308             fHistNdim->Fill(vec);
1309
1310             fHistRap->Fill(rapidity);
1311             //if (pointingAngleH > 0.1) continue;
1312             if (pointingAngleH > 0.05) continue;
1313
1314             ///////////////////////////
1315             //MC part for Associated H
1316             ///////////////////////////
1317
1318             if (HasMC() && isCorrectlyAssociatedLambda) {
1319               Int_t labelP = trackP->GetLabel();
1320               TParticle *tparticleDaughter = stack->Particle(TMath::Abs(labelP));
1321               Int_t labelMother = tparticleDaughter->GetFirstMother();
1322               TParticle *tparticleMother = stack->Particle(TMath::Abs(labelMother));
1323
1324               Int_t labelProton = trackP->GetLabel();
1325               Int_t labelPion = trackN->GetLabel();
1326
1327               //H-Dibaryon
1328               if(tparticleMother->GetPdgCode() == pdgHDibaryon && labelAssociatedH==labelMother){ //check mother PDG
1329                 Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1330                 Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1331                 Int_t labelSecondDaughter = labelFirstDaughter +1;
1332
1333                 TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1334                 TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1335
1336                 TLorentzVector ppi;
1337                 TLorentzVector lpi;
1338                 TLorentzVector lP;
1339
1340                 ppi=posProt+negPio;
1341                 lpi=lambdaH+negPio;
1342                 lP=lambdaH+posProt;
1343
1344                 if((tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter)&&(tparticleSecondDaughter->GetPdgCode() == pdgProton||tparticleSecondDaughter->GetPdgCode() == pdgPionMinus) && labelProton==labelSecondDaughter) fHistMassPpi->Fill(ppi.M());
1345
1346                 if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter) fHistMassLpi->Fill(lpi.M());
1347
1348                 if(tparticleSecondDaughter->GetPdgCode() == pdgProton && labelProton==labelSecondDaughter) fHistMassLP->Fill(lP.M());
1349
1350                 if(tparticleSecondDaughter->GetPdgCode() == pdgProton && labelProton==labelSecondDaughter){//check second daughter PDG
1351                   if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter){//check second daughter PDG
1352
1353                     fHistHDibaryonInvaMassAso->Fill(hDibaryon.M());
1354
1355                     Double_t distance01=vecDist.Mag();
1356                     fHistMCdcaPvtxDvtx->Fill(decayLengthH);
1357                     fHistMCdcaPvtxLvtx->Fill(decayLength);
1358                     fHistMCdcaDvtxLvtx->Fill(distance01);
1359                     fHistMCangleLH->Fill(hPointingAngle);
1360                     fHistMCdecayAngle->Fill(decayAngleH);
1361                     fHistMCpointingAngle->Fill(pointingAngleH);
1362                     fHistMCap->Fill(alfa,qt);
1363
1364                     fHistHilf1->Fill(posPionKF.GetDistanceFromVertex(primVtx));
1365                     fHistHilf2->Fill(protonKF.GetDistanceFromVertex(primVtx));
1366                     fHistHilf3->Fill(protonKF.GetDistanceFromVertex(posPionKF));
1367                     fHistHilf6->Fill(dca);
1368                     fHistPtvsYAso->Fill(hDibaryon.Pt(),hDibaryon.Rapidity());
1369                     fHistPtvsEtaAso->Fill(hDibaryon.Pt(),hDibaryon.Eta());
1370                     mcStatus=1;
1371                   }//end check for third daughter PDG
1372                 }//end check second daughter PDG
1373               }//end H-Dibaryon
1374             }//end MC
1375
1376             //      cout<<"Trigger: "<<triggertype<<endl;
1377             fHistMassH->Fill(hDibaryon.M());
1378             fHistMassHcentMult->Fill(hDibaryon.M(),triggertype,refMultTpc);
1379             ppK=lambdaH+posProt;
1380             fHistMassLambdaP->Fill(ppK.M());
1381
1382             //  fHistNdim = new THnSparseF("fHistNdim","THnS;InvMass, InvMassLambda, pointingAngle, armPoAlpha, armPoQt, pTL, pTH, d0p, d0n, dcaHd, dca, decayL, cosPA, centr, multi, mcinf;InvMassH", 16,binsD01,xminD01,xmaxD01);
1383
1384             //      Double_t vec[16]={hDibaryon.M(), lInvMassLambda, pointingAngleH, alfa, qt, lPtLambda, hDibaryon.Pt(), posPionKF.GetDistanceFromVertex(primVtx), protonKF.GetDistanceFromVertex(primVtx), dca, protonKF.GetDistanceFromVertex(posPionKF), TMath::Cos(pointingAngleH), centrPerc, refMultTpc, mcStatus};
1385             //      fHistNdim->Fill(vec);
1386
1387           }
1388         }
1389       }
1390     }
1391   }
1392
1393   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394   //Pure MC Part!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1395   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1396
1397   // Monte Carlo for genenerated particles
1398   if (HasMC()) //MC loop
1399     {
1400
1401       Int_t stackN = 0;
1402
1403       Double_t momentumPionGen[3]={0,0,0};
1404       Double_t momentumNucleonGen[3]={0,0,0};
1405       Double_t momentumLambdaGen[3]={0,0,0};
1406
1407       Double_t energyPionGen = 0;
1408       Double_t energyNucleonGen = 0;
1409       Double_t energyLambdaGen = 0;
1410
1411       Double_t transversMomentumMotherGen = 0;
1412       Double_t longitudinalMomentumMotherGen = 0;
1413       Double_t totalEnergyMotherGen = 0;
1414
1415       Double_t rapidityGen = 2;
1416
1417       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1418         {
1419
1420           TParticle *tparticleMother = stack->Particle(stackN);
1421
1422           if(tparticleMother->GetPdgCode() == pdgLambda) fHistCount->Fill(16);
1423
1424           //H-Dibaryon
1425           if(tparticleMother->GetPdgCode() == pdgHDibaryon) //check mother PDG
1426             {
1427               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1428               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1429               Int_t labelSecondDaughter = labelFirstDaughter +1;
1430
1431               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1432               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1433               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1434
1435               if(tparticleFirstDaughter->GetPdgCode() == pdgLambda) //check first daughter PDG
1436                 {
1437                   if(tparticleSecondDaughter->GetPdgCode() == pdgProton)//check second daughter PDG
1438                     {
1439                       if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus)//check second daughter PDG
1440                         {
1441                           momentumLambdaGen[0] = tparticleFirstDaughter->Px();
1442                           momentumLambdaGen[1] = tparticleFirstDaughter->Py();
1443                           momentumLambdaGen[2] = tparticleFirstDaughter->Pz();
1444
1445                           momentumNucleonGen[0] = tparticleSecondDaughter->Px();
1446                           momentumNucleonGen[1] = tparticleSecondDaughter->Py();
1447                           momentumNucleonGen[2] = tparticleSecondDaughter->Pz();
1448
1449                           momentumPionGen[0] = tparticleThirdDaughter->Px();
1450                           momentumPionGen[1] = tparticleThirdDaughter->Py();
1451                           momentumPionGen[2] = tparticleThirdDaughter->Pz();
1452
1453                           TLorentzVector lorentzVectorLambda;
1454                           TLorentzVector lorentzVectorProton;
1455                           TLorentzVector lorentzVectorPion;
1456                           TLorentzVector lorentzVectorHDibaryon;
1457
1458                           lorentzVectorLambda.SetXYZM(momentumLambdaGen[0],momentumLambdaGen[1],momentumLambdaGen[2],1.115);
1459                           lorentzVectorProton.SetXYZM(momentumNucleonGen[0],momentumNucleonGen[1],momentumNucleonGen[2],protonMass);
1460                           lorentzVectorPion.SetXYZM(momentumPionGen[0],momentumPionGen[1],momentumPionGen[2],pionMass);
1461
1462                           lorentzVectorHDibaryon = lorentzVectorLambda + lorentzVectorProton + lorentzVectorPion;
1463                           rapidityGen=lorentzVectorHDibaryon.Rapidity();
1464                           transversMomentumMotherGen = lorentzVectorHDibaryon.Pt();
1465                           longitudinalMomentumMotherGen = lorentzVectorHDibaryon.Pz();
1466                           totalEnergyMotherGen = lorentzVectorHDibaryon.Energy();
1467
1468                           if(rapidityGen > 1.0 || rapidityGen < -1 ) continue;
1469                           //lorentzVectorLambda
1470                           fHistHDibaryonInvaMassGen->Fill(lorentzVectorHDibaryon.M());
1471                           if (lorentzVectorLambda.Rapidity()  > 1.0 || lorentzVectorLambda.Rapidity() < -1) continue;
1472                           if (lorentzVectorProton.Rapidity()  > 1.0 ||      lorentzVectorProton.Rapidity() < -1) continue;
1473
1474                           if (lorentzVectorPion.Rapidity()  > 1.0 ||  lorentzVectorPion.Rapidity() < -1) continue;
1475                           fHistHDibaryonInvaMassGenRes->Fill(lorentzVectorHDibaryon.M());
1476                           fHistPtvsEtaGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Eta());
1477                           fHistPtvsYGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Rapidity());
1478                           fHistPtvsEtaGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Eta());
1479                           fHistCount->Fill(11);
1480                         }//end of check third daughter PDG
1481                     }//end of check second daughter PDG
1482                 }//end of check first daughter PDG
1483             }//end of H-Dibaryon
1484
1485           //Anti-H-Dibaryon
1486           if(tparticleMother->GetPdgCode() == pdgAntiHDibaryon) //check mother PDG
1487             {
1488               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1489               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1490               Int_t labelSecondDaughter = labelFirstDaughter +1;
1491
1492               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1493               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1494               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1495
1496               if(tparticleFirstDaughter->GetPdgCode() == pdgAntiLambda) //check first daughter PDG
1497                 {
1498                   if(tparticleSecondDaughter->GetPdgCode() == pdgAntiProton)//check second daughter PDG
1499                     {
1500                       if(tparticleThirdDaughter->GetPdgCode() == pdgPionPlus)//check second daughter PDG
1501                         {
1502                           momentumLambdaGen[0] = tparticleFirstDaughter->Px();
1503                           momentumLambdaGen[1] = tparticleFirstDaughter->Py();
1504                           momentumLambdaGen[2] = tparticleFirstDaughter->Pz();
1505
1506                           momentumNucleonGen[0] = tparticleSecondDaughter->Px();
1507                           momentumNucleonGen[1] = tparticleSecondDaughter->Py();
1508                           momentumNucleonGen[2] = tparticleSecondDaughter->Pz();
1509
1510                           momentumPionGen[0] = tparticleThirdDaughter->Px();
1511                           momentumPionGen[1] = tparticleThirdDaughter->Py();
1512                           momentumPionGen[2] = tparticleThirdDaughter->Pz();
1513
1514                           energyLambdaGen  = tparticleFirstDaughter->Energy();
1515                           energyNucleonGen = tparticleSecondDaughter->Energy();
1516                           energyPionGen    = tparticleThirdDaughter->Energy();
1517
1518                           TLorentzVector lorentzVectorLambda;
1519                           TLorentzVector lorentzVectorProton;
1520                           TLorentzVector lorentzVectorPion;
1521                           TLorentzVector lorentzVectorHDibaryon;
1522
1523                           lorentzVectorLambda.SetXYZM(momentumLambdaGen[0],momentumLambdaGen[1],momentumLambdaGen[2],1.115);
1524                           lorentzVectorProton.SetXYZM(momentumNucleonGen[0],momentumNucleonGen[1],momentumNucleonGen[2],protonMass);
1525                           lorentzVectorPion.SetXYZM(momentumPionGen[0],momentumPionGen[1],momentumPionGen[2],pionMass);
1526
1527                           lorentzVectorHDibaryon = lorentzVectorLambda + lorentzVectorProton + lorentzVectorPion;
1528
1529                           rapidityGen=lorentzVectorHDibaryon.Rapidity();
1530                           if(rapidityGen > 1.0 || rapidityGen < -1 ) continue;
1531                           fHistAntiHDibaryonInvaMassGen->Fill(lorentzVectorHDibaryon.M());
1532                         }//end of check third daughter PDG
1533                     }//end of check second daughter PDG
1534                 }//end of check first daughter PDG
1535             }//end of Anti-H-Dibaryon
1536         }
1537     }//end MC
1538
1539   // Post output data.
1540   PostData(1,fHistList);
1541   //PostData(0,fHistList);
1542
1543   if (listCrossV0) delete listCrossV0;
1544   if (esdVer1) delete esdVer1;
1545   if (vertexer) delete vertexer;
1546   if (vertexer1) delete vertexer1;
1547   if (trkArray) delete trkArray;
1548   if (trkArray1) delete trkArray1;
1549 }
1550
1551 //________________________________________________________________________
1552 void AliAnalysisTaskHdibaryonLPpi::Terminate(Option_t *)
1553 {
1554   // Draw result to the screen
1555   // Called once at the end of the query
1556
1557 }
1558