PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskLambdaNAOD.cxx

   1 /**************************************************************************
   2  * Author : Nicole Alice Martin (nicole.alice.martin@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 //               AliAnalysisTaskLambdaNAOD  class
  17 //  task for the investigation of (anti-)lambda-n bound state
  18 //          uses the V0 finders, based on AODs or ESDS
  19 //-----------------------------------------------------------------
  20
  21
  22 class TTree;
  23 class TParticle;
  24 class TVector3;
  25
  26 class AliESDVertex;
  27 class AliESDv0;
  28
  29 class AliAODVertex;
  30 class AliAODv0;
  31
  32 #include "TChain.h"
  33 #include <Riostream.h>
  34 #include <iostream>
  35 #include <fstream>
  36 #include <string>
  37 #include "TList.h"
  38 #include "TH1.h"
  39 #include "TH2.h"
  40 #include "TH3.h"
  41 #include "THnSparse.h"
  42 #include "TCanvas.h"
  43 #include "TMath.h"
  44 #include "TLegend.h"
  45 #include "TFile.h"
  46 #include <TRandom3.h>
  47 #include "TPDGCode.h"
  48 #include <TDatabasePDG.h>
  49
  50 #include "AliAnalysisManager.h"
  51 #include "AliAnalysisTask.h"
  52 #include "AliInputEventHandler.h"
  53 #include "AliLog.h"
  54 #include "AliCentrality.h"
  55 #include "AliV0vertexer.h"
  56 #include "AliPIDResponse.h"
  57 #include "AliMultiplicity.h"
  58 #include "AliVertexerTracks.h"
  59
  60 #include "AliVEvent.h"
  61 #include "AliVTrack.h"
  62
  63 #include "AliESDInputHandler.h"
  64 #include "AliESDEvent.h"
  65 #include "AliESDtrackCuts.h"
  66 #include "AliESDpid.h"
  67 #include "AliESDv0.h"
  68
  69 #include "AliAODInputHandler.h"
  70 #include "AliAODEvent.h"
  71 #include "AliAODTrack.h"
  72 #include "AliAODv0.h"
  73
  74 #include "AliCFContainer.h"
  75 #include "AliKFParticle.h"
  76 #include "AliKFVertex.h"
  77
  78 #include "AliMCEventHandler.h"
  79 #include "AliMCEvent.h"
  80 #include "AliStack.h"
  81
  82 #include "AliAnalysisTaskLambdaNAOD.h"
  83 using namespace std;
  84
  85
  86 ClassImp(AliAnalysisTaskLambdaNAOD)
  87
  88
  89 //________________________________________________________________________
  90 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD()
  91 : AliAnalysisTaskSE(),
  92   fAnalysisType("ESD"),
  93   fEventHandler(0),
  94   fESDtrackCutsV0(0),
  95   fESDpid(0),
  96   fPIDResponse(0),
  97   fTreeV0(0),
  98   fHistNumberOfEvents(0),
  99   fHistTriggerStat(0),
 100   fHistLambdaNeutronPtGen(0),
 101   fHistAntiLambdaNeutronPtGen(0),
 102   fHistLambdaNeutronInvaMassGen(0),
 103   fHistAntiLambdaNeutronInvaMassGen(0),
 104   fHistLambdaNeutronDecayLengthGen(0),
 105   fHistAntiLambdaNeutronDecayLengthGen(0),
 106   fHistLambdaNeutronPtAso(0),
 107   fHistLambdaNeutronPtAsoCuts(0),
 108   fHistAntiLambdaNeutronPtAso(0),
 109   fHistAntiLambdaNeutronPtAsoCuts(0),
 110   fHistLambdaNeutronInvaMassAso(0),
 111   fHistAntiLambdaNeutronInvaMassAso(0),
 112   fHistLambdaNeutronDecayLengthAso(0),
 113   fHistAntiLambdaNeutronDecayLengthAso(0),
 114   fof(0),
 115   fHistArmenterosPodolanskiDeuteronPion(0),
 116   fHistArmenterosPodolanskiAntiDeuteronPion(0),
 117   fHistDeDxQA(0),
 118 //fHistdEdxMC(0),
 119   fNTriggers(5),
 120   fMCtrue(0),
 121   fOnlyQA(0),
 122   fTriggerFired(),
 123 //fV0object(NULL),
 124   fItrk(0),
 125   fOutputContainer(NULL)
 126 {
 127   // default Constructor
 128
 129   // Define input and output slots here
 130 }
 131
 132 //________________________________________________________________________
 133 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD(const char *name)
 134   : AliAnalysisTaskSE(name),
 135     fAnalysisType("ESD"),
 136     fEventHandler(0),
 137     fESDtrackCutsV0(0),
 138     fESDpid(0),
 139     fPIDResponse(0),
 140     fTreeV0(0),
 141     fHistNumberOfEvents(0),
 142     fHistTriggerStat(0),
 143     fHistLambdaNeutronPtGen(0),
 144     fHistAntiLambdaNeutronPtGen(0),
 145     fHistLambdaNeutronInvaMassGen(0),
 146     fHistAntiLambdaNeutronInvaMassGen(0),
 147     fHistLambdaNeutronDecayLengthGen(0),
 148     fHistAntiLambdaNeutronDecayLengthGen(0),
 149     fHistLambdaNeutronPtAso(0),
 150     fHistLambdaNeutronPtAsoCuts(0),
 151     fHistAntiLambdaNeutronPtAso(0),
 152     fHistAntiLambdaNeutronPtAsoCuts(0),
 153     fHistLambdaNeutronInvaMassAso(0),
 154     fHistAntiLambdaNeutronInvaMassAso(0),
 155     fHistLambdaNeutronDecayLengthAso(0),
 156     fHistAntiLambdaNeutronDecayLengthAso(0),
 157     fof(0),
 158     fHistArmenterosPodolanskiDeuteronPion(0),
 159     fHistArmenterosPodolanskiAntiDeuteronPion(0),
 160     fHistDeDxQA(0),
 161     //fHistdEdxMC(0),
 162     fNTriggers(5),
 163     fMCtrue(0),
 164     fOnlyQA(0),
 165     fTriggerFired(),
 166     //fV0object(NULL),
 167     fItrk(0),
 168     fOutputContainer(NULL)
 169 {
 170   // Constructor
 171
 172   // Define input and output slots here
 173   // Input slot #0 works with a TChain
 174   DefineInput(0, TChain::Class());
 175   // Output slot #0 writes into a TH1 container
 176   DefineOutput(1, TObjArray::Class());
 177   DefineOutput(2, TTree::Class());
 178
 179   //ESD Track Cuts for v0's
 180   fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
 181   fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
 182   fESDtrackCutsV0->SetMinNClustersTPC(60);
 183   fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
 184   fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
 185   //  fESDtrackCutsV0->SetMinNClustersITS(1); // to be tested if this cut is not too strong
 186   fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
 187
 188
 189   fMCtrue = kTRUE;
 190   fOnlyQA = kTRUE;
 191
 192 }
 193
 194 //____________________________________________________________
 195 AliAnalysisTaskLambdaNAOD::~AliAnalysisTaskLambdaNAOD(){
 196   //
 197   // Destructor
 198   //
 199   //if (fESD) delete fESD;
 200   if (fESDtrackCutsV0) delete fESDtrackCutsV0;
 201
 202 }
 203
 204 //____________________________________________________________
 205 const Int_t AliAnalysisTaskLambdaNAOD::fgkPdgCode[] = {
 206   211,                //PionPlus
 207   -211,               //PionMinus
 208   2212,               //Proton
 209   -2212,              //Anti-Proton
 210   1000010020,         //Deuteron
 211   -1000010020,        //Anti-Deuteron
 212   1000020030,         //Helium3
 213   -1000020030,        //Anti-Helium3
 214   3122,               //Lambda
 215   -3122,              //Anti-Lambda
 216   1010000020,         //LambdaNeutron
 217   -1010000020,        //Anti-Lambda-Neutron
 218   1010010030,         //Hypertriton
 219   -1010010030         //Anti-Hypertriton
 220 };
 221
 222 //____________________________________________________________
 223 const Double_t AliAnalysisTaskLambdaNAOD::fgkMass[] = {
 224   0.13957,           //Pion
 225   0.93827,           //Proton
 226   1.875612,          //Deuteron
 227   2.808921,          //Triton
 228   2.80941            //Helium3 Quelle: Wolfram Alpha
 229 };
 230
 231 //________________________________________________________________________
 232 void AliAnalysisTaskLambdaNAOD::UserCreateOutputObjects(){
 233
 234   // New PID object
 235   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 236   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 237   fPIDResponse = inputHandler->GetPIDResponse();
 238
 239   // Create histograms for MC
 240   //generated histogramms
 241   fHistLambdaNeutronPtGen = new TH1F("fHistLambdaNeutronPtGen", "Generated  #Lambdan", 201, 0., 10.1);
 242   fHistLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
 243   fHistLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 244
 245   fHistAntiLambdaNeutronPtGen = new TH1F("fHistAntiLambdaNeutronPtGen", "Generated  #bar{#Lambdan}", 201, 0., 10.1);
 246   fHistAntiLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
 247   fHistAntiLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 248
 249   fHistLambdaNeutronInvaMassGen = new TH1F("fHistLambdaNeutronInvaMassGen", "Generated #Lambdan ", 100, 2.0, 2.1);
 250   fHistLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
 251   fHistLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
 252
 253   fHistAntiLambdaNeutronInvaMassGen = new TH1F("fHistAntiLambdaNeutronInvaMassGen", "Generated  #bar{#Lambdan}", 100, 2.0, 2.1);
 254   fHistAntiLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
 255   fHistAntiLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
 256
 257   fHistLambdaNeutronDecayLengthGen = new TH1F("fHistLambdaNeutronDecayLengthGen", "Generated  #Lambdan", 401, 0., 400.1);
 258   fHistLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
 259   fHistLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 260
 261   fHistAntiLambdaNeutronDecayLengthGen = new TH1F("fHistAntiLambdaNeutronDecayLengthGen", "Generated #bar{#Lambdan}", 401, 0., 400.1);
 262   fHistAntiLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
 263   fHistAntiLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 264
 265   //assoziated (reconstracted) histogramms
 266   fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated   #Lambdan", 201, 0., 10.1);
 267   fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 268   fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 269
 270   fHistLambdaNeutronPtAsoCuts = new TH1F("fHistLambdaNeutronPtAsoCuts", "Associated   #Lambdan Cuts", 201, 0., 10.1);
 271   fHistLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
 272   fHistLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 273
 274   fHistAntiLambdaNeutronPtAsoCuts = new TH1F("fHistAntiLambdaNeutronPtAsoCuts", " associated  #bar{#Lambdan} Cuts", 201, 0., 10.1);
 275   fHistAntiLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
 276   fHistAntiLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 277
 278   fHistAntiLambdaNeutronPtAso = new TH1F("fHistAntiLambdaNeutronPtAso", " associated  #bar{#Lambdan}", 201, 0., 10.1);
 279   fHistAntiLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 280   fHistAntiLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 281
 282   fHistLambdaNeutronInvaMassAso = new TH1F("fHistLambdaNeutronInvaMassAso", "Associated  #Lambdan", 100, 2.0, 2.1);
 283   fHistLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
 284   fHistLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
 285
 286   fHistAntiLambdaNeutronInvaMassAso = new TH1F("fHistAntiLambdaNeutronInvaMassAso", " Associated  #bar{#Lambdan}", 100, 2.0, 2.1);
 287   fHistAntiLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
 288   fHistAntiLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+})  (GeV/#it{c}^{2})");
 289
 290   fHistLambdaNeutronDecayLengthAso = new TH1F("fHistLambdaNeutronDecayLengthAso", "Associated  #Lambdan", 401, 0., 400.1);
 291   fHistLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
 292   fHistLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 293
 294   fHistAntiLambdaNeutronDecayLengthAso = new TH1F("fHistAntiLambdaNeutronDecayLengthAso", "Associated #bar{#Lambdan}", 401, 0., 400.1);
 295   fHistAntiLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
 296   fHistAntiLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 297
 298   fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated   #Lambdan", 201, 0., 10.1);
 299   fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 300   fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 301
 302   //Tree
 303   //------------ Tree and branch definitions ----------------//
 304   fTreeV0 = new TTree("tree", "fTreeV0");
 305   fTreeV0->Branch("fItrk", &fItrk, "fItrk/I");
 306
 307   //fTreeV0->Branch("fV0object",&fV0object,"fV0object[fItrk]");
 308   fTreeV0->Branch("fV0finder",fV0finder,"fV0finder[fItrk]/I");
 309   fTreeV0->Branch("fkMB",fkMB,"fkMB[fItrk]/I");
 310   fTreeV0->Branch("fkCentral",fkCentral,"fkCentral[fItrk]/I");
 311   fTreeV0->Branch("fkSemiCentral",fkSemiCentral,"fkSemiCentral[fItrk]/I");
 312   fTreeV0->Branch("fkEMCEJE",fkEMCEJE,"fkEMCEJE[fItrk]/I");
 313   fTreeV0->Branch("fkEMCEGA",fkEMCEGA,"fkEMCEGA[fItrk]/I");
 314
 315   fTreeV0->Branch("fPtotN",fPtotN,"fPtotN[fItrk]/D");
 316   fTreeV0->Branch("fPtotP",fPtotP,"fPtotP[fItrk]/D");
 317   fTreeV0->Branch("fMotherPt",fMotherPt,"fMotherPt[fItrk]/D");
 318   fTreeV0->Branch("fdEdxN",fdEdxN,"fdEdxN[fItrk]/D");
 319   fTreeV0->Branch("fdEdxP",fdEdxP,"fdEdxP[fItrk]/D");
 320   fTreeV0->Branch("fSignN",fSignN,"fSignN[fItrk]/D");
 321   fTreeV0->Branch("fSignP",fSignP,"fSignP[fItrk]/D");
 322
 323   fTreeV0->Branch("fDCAv0",fDCAv0,"fDCAv0[fItrk]/F"); //Dca v0 Daughters
 324   fTreeV0->Branch("fCosinePAv0",fCosinePAv0,"fCosinePAv0[fItrk]/F"); //Cosine of Pionting Angle
 325   fTreeV0->Branch("fDecayRadiusTree",fDecayRadiusTree,"fDecayRadiusTree[fItrk]/F"); //decay radius
 326
 327   fTreeV0->Branch("fChargeComboDeuteronPionTree",fChargeComboDeuteronPionTree,"fChargeComboDeuteronPionTree[fItrk]/I");
 328   fTreeV0->Branch("fInvaMassDeuteronPionTree",fInvaMassDeuteronPionTree,"fInvaMassDeuteronPionTree[fItrk]/D"); //invariant mass
 329
 330   fTreeV0->Branch("fIsCorrectlyAssociated",fIsCorrectlyAssociated,"fIsCorrectlyAssociated[fItrk]/O"); //associated hypertriton
 331
 332   fTreeV0->Branch("fAmenterosAlphaTree",fAmenterosAlphaTree,"fAmenterosAlphaTree[fItrk]/D");
 333   fTreeV0->Branch("fAmenterosQtTree",fAmenterosQtTree,"fAmenterosQtTree[fItrk]/D");
 334   fTreeV0->Branch("fRotationTree",fRotationTree,"fRotationTree[fItrk]/I");
 335
 336   //Armenteros-Podolanski
 337   fHistArmenterosPodolanskiDeuteronPion= new TH2F("fHistArmenterosPodolanskiDeuteronPion", "Armenteros-Podolanski d #pi^{-}", 200,-1.0,1.0, 500,0,1);
 338   fHistArmenterosPodolanskiDeuteronPion->GetXaxis()->SetTitle("#alpha");
 339   fHistArmenterosPodolanskiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
 340   fHistArmenterosPodolanskiDeuteronPion->SetMarkerStyle(kFullCircle);
 341
 342   fHistArmenterosPodolanskiAntiDeuteronPion= new TH2F("fHistArmenterosPodolanskiAntiDeuteronPion", "Armenteros-Podolanski #bar{d} #pi^{+}", 200,-1.0,1.0, 500,0,1);
 343   fHistArmenterosPodolanskiAntiDeuteronPion->GetXaxis()->SetTitle("#alpha");
 344   fHistArmenterosPodolanskiAntiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
 345   fHistArmenterosPodolanskiAntiDeuteronPion->SetMarkerStyle(kFullCircle);
 346
 347   // control histograms
 348   fof = new TH2F("fof", "OnTheFlyStatus ",5,0.5,5.5,2,-0.5,1.5);
 349   fof->GetYaxis()->SetBinLabel(1,"offline");
 350   fof->GetYaxis()->SetBinLabel(2,"onTheFly");
 351   fof->GetXaxis()->SetBinLabel(1,"total");
 352   fof->GetXaxis()->SetBinLabel(2,"dcaCut");
 353   fof->GetXaxis()->SetBinLabel(3,"cosCut");
 354   fof->GetXaxis()->SetBinLabel(4,"nucleonPID");
 355   fof->GetXaxis()->SetBinLabel(5,"pionPID");
 356   //fof->GetXaxis()->SetBinLabel(6,"decayRadiusCut");
 357   //fof->SetMarkerStyle(kFullCircle);
 358
 359   //histogram to count number of events
 360   fHistNumberOfEvents = new TH1F("fHistNumberOfEvents", "Number of events", 11, -0.5, 10.5);
 361   fHistNumberOfEvents ->GetXaxis()->SetTitle("Centrality");
 362   fHistNumberOfEvents ->GetYaxis()->SetTitle("Entries");
 363
 364   //trigger statitics histogram
 365   fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
 366
 367   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 368   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 369     fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 370
 371   //QA dE/dx
 372   fHistDeDxQA = new TH3F("fHistDeDxQA", "QA dE/dx", 400, -6.0, 6.0, 500, 0.0, 2000, 15, -0.5, 14.5);
 373   fHistDeDxQA->GetYaxis()->SetTitle("TPC Signal");
 374   fHistDeDxQA->GetXaxis()->SetTitle("p (GeV/c)");
 375   fHistDeDxQA->GetZaxis()->SetBinLabel(0,"all pos v0 tracks");
 376   fHistDeDxQA->GetZaxis()->SetBinLabel(1,"all neg v0 tracks");
 377   fHistDeDxQA->GetZaxis()->SetBinLabel(2,"all neg deuteron");
 378   fHistDeDxQA->GetZaxis()->SetBinLabel(3,"all pos deuteron");
 379   fHistDeDxQA->GetZaxis()->SetBinLabel(6,"all selected tracks");
 380   fHistDeDxQA->GetZaxis()->SetBinLabel(7,"neg deuteron for deuteron+pion");
 381   fHistDeDxQA->GetZaxis()->SetBinLabel(8,"pos pion for deuteron+pion");
 382   fHistDeDxQA->GetZaxis()->SetBinLabel(9,"pos deuteron for deuteron+pion");
 383   fHistDeDxQA->GetZaxis()->SetBinLabel(10,"neg pion for deuteron+pion");
 384
 385   //Define and fill the OutputContainer
 386   fOutputContainer = new TObjArray(1);
 387   fOutputContainer->SetOwner(kTRUE);
 388   fOutputContainer->SetName(GetName()) ;
 389   fOutputContainer->Add(fHistArmenterosPodolanskiDeuteronPion);
 390   fOutputContainer->Add(fHistArmenterosPodolanskiAntiDeuteronPion);
 391   fOutputContainer->Add(fof);
 392   fOutputContainer->Add(fHistDeDxQA);
 393   fOutputContainer->Add(fHistNumberOfEvents);
 394   fOutputContainer->Add(fHistTriggerStat);
 395   fOutputContainer->Add(fHistLambdaNeutronPtGen);
 396   fOutputContainer->Add(fHistAntiLambdaNeutronPtGen);
 397   fOutputContainer->Add(fHistLambdaNeutronInvaMassGen);
 398   fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassGen);
 399   fOutputContainer->Add(fHistLambdaNeutronDecayLengthGen);
 400   fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthGen);
 401   fOutputContainer->Add(fHistLambdaNeutronPtAso);
 402   fOutputContainer->Add(fHistLambdaNeutronPtAsoCuts);
 403   fOutputContainer->Add(fHistAntiLambdaNeutronPtAso);
 404   fOutputContainer->Add(fHistAntiLambdaNeutronPtAsoCuts);
 405   fOutputContainer->Add(fHistLambdaNeutronInvaMassAso);
 406   fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassAso);
 407   fOutputContainer->Add(fHistLambdaNeutronDecayLengthAso);
 408   fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthAso);
 409 }
 410
 411 //________________________________________________________________________
 412 void AliAnalysisTaskLambdaNAOD::UserExec(Option_t *){
 413
 414   // Main loop
 415   // Called for each event
 416
 417   //Data
 418   //get Event-Handler for the trigger information
 419   fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 420   if (!fEventHandler) {
 421     AliError("Could not get InputHandler");
 422     //return -1;
 423     return;
 424   }
 425
 426   // Monte Carlo
 427   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 428   if (!eventHandler){
 429     //Printf("ERROR: Could not retrieve MC event handler");
 430     fMCtrue = kFALSE;
 431   }
 432
 433   AliMCEvent* mcEvent = 0x0;
 434   AliStack* stack = 0x0;
 435   if (eventHandler) mcEvent = eventHandler->MCEvent();
 436   if (!mcEvent){
 437     //Printf("ERROR: Could not retrieve MC event");
 438     if (fMCtrue) return;
 439   }
 440
 441   if (fMCtrue){
 442     stack = mcEvent->Stack();
 443     if (!stack) return;
 444   }
 445
 446   //look for the generated particles
 447   MCGenerated(stack);
 448
 449
 450   if (SetupEvent() < 0) {
 451     PostData(1, fOutputContainer);
 452     return;
 453     }
 454
 455   //DATA
 456   AliESDEvent *fESDevent = 0x0;
 457   AliAODEvent *fAODevent = 0x0;
 458
 459
 460   if(!fMCtrue){
 461     if(!fPIDResponse) {
 462       AliError("Cannot get pid response");
 463       return;
 464     }
 465   }
 466
 467
 468   Int_t centrality = -1;
 469   Double_t vertex[3]          = {-100.0, -100.0, -100.0};
 470
 471   //Initialisation of the event and basic event cuts:
 472   //1.) ESDs:
 473   if (fAnalysisType == "ESD") {
 474
 475     fESDevent=dynamic_cast<AliESDEvent*>(InputEvent());
 476     if (!fESDevent) {
 477       AliWarning("ERROR: fESDevent not available \n");
 478       return;
 479     }
 480
 481     //Basic event cuts:
 482     //1.) vertex existence
 483     const AliESDVertex *vertexESD = fESDevent->GetPrimaryVertexTracks();
 484     if (vertexESD->GetNContributors()<1)
 485       {
 486         // SPD vertex
 487         vertexESD = fESDevent->GetPrimaryVertexSPD();
 488         if(vertexESD->GetNContributors()<1) {
 489           PostData(1, fOutputContainer);
 490           return;
 491         }
 492       }
 493
 494     vertexESD->GetXYZ(vertex);
 495
 496     //2. vertex position within 10 cm
 497     if (TMath::Abs(vertexESD->GetZv()) > 10) return;
 498
 499     //centrality selection in PbPb
 500     if (fESDevent->GetEventSpecie() == 4) //species == 4 == PbPb
 501       { // PbPb
 502         AliCentrality *esdCentrality = fESDevent->GetCentrality();
 503         centrality = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 504         if(!fMCtrue){
 505           if (centrality < 0. || centrality > 8. ) return; //select only events with centralities between 0 and 80 %
 506         }
 507       }
 508
 509     //cout << "centrality "<< centrality << endl;
 510     // count number of events
 511     fHistNumberOfEvents->Fill(centrality);
 512
 513  }
 514
 515 //2.) AODs:
 516  else if (fAnalysisType == "AOD") {
 517
 518    fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 519    if (!fAODevent) {
 520      AliWarning("ERROR: lAODevent not available \n");
 521      return;
 522    }
 523
 524    const AliAODVertex *vertexAOD = fAODevent->GetPrimaryVertex();
 525     if (vertexAOD->GetNContributors()<1)
 526       {
 527         // SPD vertex
 528         vertexAOD = fAODevent->GetPrimaryVertexSPD();
 529         if(vertexAOD->GetNContributors()<1) {
 530           PostData(1, fOutputContainer);
 531           return;
 532         }
 533       }
 534     vertexAOD->GetXYZ(vertex);
 535
 536     //2. vertex position within 10 cm
 537     if (TMath::Abs(vertex[2]) > 10) return;
 538
 539     //centrality selection
 540    AliCentrality *aodCentrality = fAODevent->GetCentrality();
 541    centrality = aodCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 542    if (centrality < 0 || centrality > 8) return; //select only events with centralities between 0 and 80 %
 543
 544    // count number of events
 545    fHistNumberOfEvents->Fill(centrality);
 546
 547  } else {
 548
 549    Printf("Analysis type (ESD or AOD) not specified \n");
 550    return;
 551
 552  }
 553
 554
 555 //v0 loop
 556
 557
 558  fItrk = 0;
 559
 560  Int_t runNumber = 0;
 561  runNumber = (InputEvent())->GetRunNumber();
 562
 563  Int_t    nTrackMultiplicity             = -1;
 564  nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
 565
 566
 567  for (Int_t ivertex=0; ivertex<(InputEvent()->GetNumberOfV0s()); ivertex++) { //beginn v0 loop
 568
 569    AliESDv0 * v0ESD = 0x0;
 570    AliAODv0 * v0AOD = 0x0;
 571
 572    Bool_t v0ChargesCorrect = kTRUE;
 573    Bool_t testTrackCuts = kFALSE;
 574    //Bool_t testFilterBit = kFALSE;
 575
 576    Int_t of = 7;
 577    Int_t onFlyStatus = 5;
 578
 579    Float_t dcaV0 = -1;
 580    Float_t cosPointing= 2;
 581    Float_t decayRadius= -1;
 582
 583    AliVTrack * trackN = 0x0;
 584    AliVTrack * trackP = 0x0;
 585
 586    Double_t ptotN = -1000;
 587    Double_t ptotP = -1000;
 588
 589    Int_t chargeComboDeuteronPion = -1;
 590
 591    Double_t momentumPion[3]={0,0,0};
 592    Double_t momentumPionRot[3]={0,0,0};
 593    Double_t momentumDeuteron[3]={0,0,0};
 594    Double_t momentumDeuteronRot[3]={0,0,0};
 595    Double_t momentumMother[3] = {0,0,0};
 596
 597    Double_t transversMomentumPion = 0;
 598    Double_t transversMomentumDeuteron = 0;
 599
 600    Double_t transversMomentumMother = 0;
 601    Double_t longitudinalMomentumMother = 0;
 602
 603    Double_t totalEnergyMother = 0;
 604    Double_t energyPion = 0;
 605    Double_t energyDeuteron = 0;
 606
 607    Double_t rapidity = 2;
 608
 609    TVector3 vecPion(0,0,0);
 610    TVector3 vecDeuteron(0,0,0);
 611    TVector3 vecMother(0,0,0);
 612
 613    Double_t alpha = 2;
 614    Double_t qt = -1;
 615
 616    Double_t thetaPion = 0;
 617    Double_t thetaDeuteron = 0;
 618
 619    Double_t phi=0;
 620    Double_t invaMassDeuteronPion = 0;
 621
 622    Double_t e12 = 0;
 623    Double_t r12   = 0;
 624    Double_t d22 = 0;
 625    Double_t dr22   = 0;
 626
 627    //Tree variables
 628    fV0finder[fItrk]         = -1;
 629    fkMB[fItrk]              = -1;
 630    fkCentral[fItrk]         = -1;
 631    fkSemiCentral[fItrk]     = -1;
 632    fkEMCEJE[fItrk]          = -1;
 633    fkEMCEGA[fItrk]          = -1;
 634
 635    fPtotN[fItrk]            = -1000;
 636    fPtotP[fItrk]            = -1000;
 637    fMotherPt[fItrk]         = -1000;
 638
 639    fdEdxN[fItrk]            = -1;
 640    fdEdxP[fItrk]            = -1;
 641    fSignN[fItrk]            = 0;
 642    fSignP[fItrk]            = 0;
 643
 644    fDCAv0[fItrk]            = -1;
 645    fCosinePAv0[fItrk]       = -2;
 646    fDecayRadiusTree[fItrk]  = -1;
 647
 648    fInvaMassDeuteronPionTree[fItrk] = 0;
 649    fChargeComboDeuteronPionTree[fItrk] = -1;
 650
 651    fAmenterosAlphaTree[fItrk] = 2;
 652    fAmenterosQtTree[fItrk] = -1;
 653
 654    //Get v0 object
 655    if(fAnalysisType == "ESD")v0ESD = fESDevent->GetV0(ivertex);
 656    if(fAnalysisType == "AOD")v0AOD = fAODevent->GetV0(ivertex);
 657
 658
 659    //distinguish between the two V0 finders: 0 offline V0 finder / 1 online V0 finder
 660
 661    if(fAnalysisType == "ESD") of = v0ESD->GetOnFlyStatus();
 662    if(fAnalysisType == "AOD") of = v0AOD->GetOnFlyStatus();
 663
 664    if(of)  onFlyStatus= 1;
 665    if(!of) onFlyStatus= 0;
 666
 667    if(onFlyStatus==0)fof->Fill(1,0);
 668    if(onFlyStatus==1)fof->Fill(1,1);
 669
 670    //Get dca, cos of pointing angle and decay radius
 671
 672
 673    if(fAnalysisType == "ESD")
 674      {
 675        dcaV0 = v0ESD->GetDcaV0Daughters();
 676        cosPointing = v0ESD->GetV0CosineOfPointingAngle();
 677        decayRadius = v0ESD->GetRr();
 678      }
 679
 680    if(fAnalysisType == "AOD")
 681      {
 682        dcaV0 = v0AOD->DcaV0Daughters();
 683        cosPointing = v0AOD->CosPointingAngle(vertex);
 684        decayRadius = v0AOD->DecayLengthV0(vertex);
 685      }
 686
 687
 688       // select coresponding tracks
 689       if(fAnalysisType == "ESD")
 690         {
 691           trackN = fESDevent->GetTrack(v0ESD->GetIndex(0));
 692           trackP = fESDevent->GetTrack(v0ESD->GetIndex(1));
 693
 694           if (trackN->Charge() > 0) // switch because of bug in V0 interface
 695             {
 696               trackN = fESDevent->GetTrack(v0ESD->GetIndex(1));
 697               trackP = fESDevent->GetTrack(v0ESD->GetIndex(0));
 698               v0ChargesCorrect = kFALSE;
 699             }
 700           //Track-Cuts
 701           testTrackCuts = TrackCuts(trackN,testTrackCuts);
 702           if(testTrackCuts == kFALSE) continue;
 703           testTrackCuts = TrackCuts(trackP,testTrackCuts);
 704           if(testTrackCuts == kFALSE) continue;
 705         }
 706
 707       if(fAnalysisType == "AOD")
 708         {
 709           trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
 710           trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
 711
 712           if (trackN->Charge() > 0) // switch because of bug in V0 interface
 713             {
 714               trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
 715               trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
 716               v0ChargesCorrect = kFALSE;
 717             }
 718           //Test-Filterbit - only use for track base analysis -- not for V0 candidates -- the AOD V0 candidates sould be a copy of the ESD candidates, BUT NOT for AOD0086 -> here there wqas a wider cos-cutto have more candidates in the AODs
 719           //testFilterBit = FilterBit(trackN,testFilterBit);
 720           //if(testFilterBit == kFALSE) continue;
 721           //testFilterBit = FilterBit(trackP,testFilterBit);
 722           //if(testFilterBit == kFALSE) continue;
 723           //Track-Cuts
 724           /*testTrackCuts = TrackCuts(trackN,testTrackCuts);
 725           if(testTrackCuts == kFALSE) continue;
 726           testTrackCuts = TrackCuts(trackP,testTrackCuts);
 727           if(testTrackCuts == kFALSE) continue;*/
 728         }
 729
 730       //Get the total momentum for each track ---> at the inner readout of the TPC???? // momenta a always positive
 731             if(fAnalysisType == "AOD") {
 732         ptotN = trackN->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
 733         ptotP = trackP->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
 734       }
 735
 736       if(fAnalysisType == "ESD") {
 737         ptotN =  MomentumInnerParam(trackN,ptotN);
 738         ptotP =  MomentumInnerParam(trackP,ptotP);
 739       }
 740
 741       //fill QA dEdx with all V0 candidates
 742       if(trackP) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),0);
 743       if(trackN) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),1);
 744
 745       if (dcaV0 > 3 || dcaV0 < 1e-20) continue;
 746
 747       //Check how much the dca cut reduces the statistic (background) for the different VO-finder
 748       if(onFlyStatus==0)fof->Fill(2,0);
 749       if(onFlyStatus==1)fof->Fill(2,1);
 750
 751       if (cosPointing < 0.9) continue;
 752
 753       //Check how much the cos-of-the-pointing-angle-cut reduces the statistic (background) for the different VO-finder
 754       if(onFlyStatus==0)fof->Fill(3,0);
 755       if(onFlyStatus==1)fof->Fill(3,1);
 756
 757       //deuteron PID via specific energy loss in the TPC
 758       Bool_t isDeuteron[3] = {kFALSE,kFALSE,kFALSE}; //0 = posDeuteron, 1 = negDeuteron, 2 = trackN is deuteron
 759
 760       DeuteronPID(trackP,trackN,ptotP,ptotN,runNumber,isDeuteron);
 761
 762       //if(!isDeuteron[0] && !isDeuteron[1]) continue;
 763       if((isDeuteron[0]==kFALSE) && (isDeuteron[1]==kFALSE)) continue;
 764
 765
 766       //Check how much the nuclei PID cuts reduce the statistics (background) for the two V0-finders
 767       if(onFlyStatus==0)fof->Fill(4,0);
 768       if(onFlyStatus==1)fof->Fill(4,1);
 769
 770       //Fill the QA dEdx with deuterons and helium3 after the nuclei PID cut
 771       if(isDeuteron[1]) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),2);
 772       if(isDeuteron[0]) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),3);
 773
 774       //deuteron PID via specific energy loss in the TPC
 775       Bool_t isPion[2] = {kFALSE,kFALSE}; //0 = posPion, 1 = negPion
 776
 777       PionPID(trackP,trackN,ptotP,ptotN,runNumber,isPion);
 778
 779       //if(isDeuteron[0] && !isPion[1]) continue; //pos deuteron and neg Pion
 780       //if(isDeuteron[1] && !isPion[0]) continue; //neg deuteron and pos Pion
 781
 782       //Check how much the pion PID cuts reduce the statistics (background) for the two V0-finders
 783       if(onFlyStatus==0)fof->Fill(5,0);
 784       if(onFlyStatus==1)fof->Fill(5,1);
 785
 786
 787       //Save the different charge combinations to differentiat between particles and anit-particles:
 788       // -/+ = Anti-deuteron + positive Pion
 789       // -/- = Anti-deuteron + negative Pion
 790       // +/- = deuteron + negative Pion
 791       // +/+ = deuteron + positive Pion
 792
 793       //Like-sign
 794       if (trackN->Charge()<0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 1; // -/-
 795       if (trackN->Charge()>0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 3; // +/+
 796
 797       //unlinke-sign
 798       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+
 799       //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+ //should not exist because of charge correction
 800       //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/- //should not exist because of charge correction
 801       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/-
 802
 803       //if(chargeComboDeuteronPion==0) cout << "chargeN: " << trackN->Charge() << " chargeP: " << trackP->Charge() << endl;
 804
 805
 806       //Fill the QA dEdx with all selected tracks after the PID cuts
 807       fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),6);
 808       fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),6);
 809
 810
 811
 812       //get the momenta of the daughters
 813
 814       if(fAnalysisType == "ESD")
 815         {
 816           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
 817
 818             v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 819             v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 820             if (!v0ChargesCorrect) {
 821
 822               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 823               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 824             }
 825           }
 826
 827           if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
 828
 829             v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 830             v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 831             if (!v0ChargesCorrect){
 832
 833               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 834               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 835             }
 836           }
 837
 838           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
 839             if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
 840
 841               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 842               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 843               if (!v0ChargesCorrect) {
 844
 845                 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 846                 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 847               }
 848             }
 849             if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
 850
 851               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 852               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 853               if (!v0ChargesCorrect) {
 854
 855                 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 856                 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 857               }
 858             }
 859           }
 860
 861
 862           //get the momenta of the mother
 863           v0ESD->GetPxPyPz(momentumMother[0],momentumMother[1],momentumMother[2]);
 864         }
 865
 866
 867       if(fAnalysisType == "AOD")
 868         {
 869           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
 870
 871             momentumPion[0] = v0AOD->MomPosX();
 872             momentumPion[1] = v0AOD->MomPosY();
 873             momentumPion[2] = v0AOD->MomPosZ();
 874
 875             momentumDeuteron[0] = v0AOD->MomNegX();
 876             momentumDeuteron[1] = v0AOD->MomNegY();
 877             momentumDeuteron[2] = v0AOD->MomNegZ();
 878
 879             if (!v0ChargesCorrect){
 880
 881               momentumPion[0] = v0AOD->MomNegX();
 882               momentumPion[1] = v0AOD->MomNegY();
 883               momentumPion[2] = v0AOD->MomNegZ();
 884
 885               momentumDeuteron[0] = v0AOD->MomPosX();
 886               momentumDeuteron[1] = v0AOD->MomPosY();
 887               momentumDeuteron[2] = v0AOD->MomPosZ();
 888             }
 889           }
 890
 891           if (chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
 892
 893             momentumPion[0] = v0AOD->MomNegX();
 894             momentumPion[1] = v0AOD->MomNegY();
 895             momentumPion[2] = v0AOD->MomNegZ();
 896
 897             momentumDeuteron[0] = v0AOD->MomPosX();
 898             momentumDeuteron[1] = v0AOD->MomPosY();
 899             momentumDeuteron[2] = v0AOD->MomPosZ();
 900
 901             if (!v0ChargesCorrect){
 902
 903               momentumPion[0] = v0AOD->MomPosX();
 904               momentumPion[1] = v0AOD->MomPosY();
 905               momentumPion[2] = v0AOD->MomPosZ();
 906
 907               momentumDeuteron[0] = v0AOD->MomNegX();
 908               momentumDeuteron[1] = v0AOD->MomNegY();
 909               momentumDeuteron[2] = v0AOD->MomNegZ();
 910             }
 911           }
 912
 913           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN correponds to deuteron or pion, trackP corresponds to deuteron or pion
 914             if(isDeuteron[2]==kTRUE){ //trackN correponds to deuteron, trackP corresponds to pion
 915
 916               momentumPion[0] = v0AOD->MomPosX();
 917               momentumPion[1] = v0AOD->MomPosY();
 918               momentumPion[2] = v0AOD->MomPosZ();
 919
 920               momentumDeuteron[0] = v0AOD->MomNegX();
 921               momentumDeuteron[1] = v0AOD->MomNegY();
 922               momentumDeuteron[2] = v0AOD->MomNegZ();
 923
 924               if (!v0ChargesCorrect){
 925
 926                 momentumPion[0] = v0AOD->MomNegX();
 927                 momentumPion[1] = v0AOD->MomNegY();
 928                 momentumPion[2] = v0AOD->MomNegZ();
 929
 930                 momentumDeuteron[0] = v0AOD->MomPosX();
 931                 momentumDeuteron[1] = v0AOD->MomPosY();
 932                 momentumDeuteron[2] = v0AOD->MomPosZ();
 933               }
 934             }
 935
 936             if(isDeuteron[2]==kFALSE){ //trackP corresponds to deuteron, trackN corresponds to pion
 937
 938               momentumPion[0] = v0AOD->MomNegX();
 939               momentumPion[1] = v0AOD->MomNegY();
 940               momentumPion[2] = v0AOD->MomNegZ();
 941
 942               momentumDeuteron[0] = v0AOD->MomPosX();
 943               momentumDeuteron[1] = v0AOD->MomPosY();
 944               momentumDeuteron[2] = v0AOD->MomPosZ();
 945
 946               if (!v0ChargesCorrect){
 947
 948                 momentumPion[0] = v0AOD->MomPosX();
 949                 momentumPion[1] = v0AOD->MomPosY();
 950                 momentumPion[2] = v0AOD->MomPosZ();
 951
 952                 momentumDeuteron[0] = v0AOD->MomNegX();
 953                 momentumDeuteron[1] = v0AOD->MomNegY();
 954                 momentumDeuteron[2] = v0AOD->MomNegZ();
 955               }
 956             }
 957
 958             //get the momenta of the mother
 959             momentumMother[0] = v0AOD->MomV0X();
 960             momentumMother[1] = v0AOD->MomV0Y();
 961             momentumMother[2] = v0AOD->MomV0Z();
 962           }
 963         }
 964
 965       //Rapidity - cut
 966       transversMomentumPion      = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]);
 967       transversMomentumDeuteron  = TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]);
 968
 969       transversMomentumMother=TMath::Sqrt((momentumDeuteron[0]+momentumPion[0])*(momentumDeuteron[0]+momentumPion[0])+(momentumDeuteron[1]+momentumPion[1])*(momentumDeuteron[1]+momentumPion[1]));
 970
 971       longitudinalMomentumMother = TMath::Sqrt((momentumDeuteron[2]+momentumPion[2])*(momentumDeuteron[2]+momentumPion[2]));
 972
 973       energyDeuteron =  TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]+momentumDeuteron[2]*momentumDeuteron[2]+fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]);
 974
 975       energyPion = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]+momentumPion[2]*momentumPion[2]+fgkMass[kMassPion]*fgkMass[kMassPion]);
 976
 977       totalEnergyMother = energyPion + energyDeuteron;
 978
 979       //rapidity cut +- 1
 980       rapidity = 0.5 * TMath::Log((totalEnergyMother+longitudinalMomentumMother)/( totalEnergyMother-longitudinalMomentumMother));
 981
 982       if(rapidity > 1.0 || rapidity < -1) continue;
 983
 984
 985       //Armanteros-Podolanski
 986       vecPion.SetXYZ(momentumPion[0],momentumPion[1],momentumPion[2]);
 987       vecDeuteron.SetXYZ(momentumDeuteron[0],momentumDeuteron[1],momentumDeuteron[2]);
 988       vecMother.SetXYZ(momentumMother[0],momentumMother[1],momentumMother[2]);
 989
 990       thetaPion = acos((vecPion * vecMother)/(vecPion.Mag() * vecMother.Mag()));
 991       thetaDeuteron = acos((vecDeuteron * vecMother)/(vecDeuteron.Mag() * vecMother.Mag()));
 992
 993       alpha = ((vecPion.Mag())*cos(thetaPion)-(vecDeuteron.Mag())*cos(thetaDeuteron))/
 994         ((vecDeuteron.Mag())*cos(thetaDeuteron)+(vecPion.Mag())*cos(thetaPion)) ;
 995       qt = vecDeuteron.Mag()*sin(thetaDeuteron);
 996
 997
 998       //Rotation for background calculation
 999       //Int_t rotation=1; // =1 signal, =2 Rotation of the pion , =3 Rotation of the deuteron
1000
1001       //Double_t fStartAnglePhi=TMath::Pi();
1002       //Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1003       //phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1004
1005       for(Int_t rotation=1;rotation<4;rotation++){ //loop for rotation
1006
1007         Double_t fStartAnglePhi=TMath::Pi();
1008         Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1009         phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1010
1011         //calculate new rotated momenta
1012         momentumPionRot[0]=TMath::Cos(phi)*momentumPion[0]-TMath::Sin(phi)*momentumPion[1];
1013         momentumPionRot[1]=TMath::Sin(phi)*momentumPion[0]+TMath::Cos(phi)*momentumPion[1];
1014
1015         momentumDeuteronRot[0]=TMath::Cos(phi)*momentumDeuteron[0]-TMath::Sin(phi)*momentumDeuteron[1];
1016         momentumDeuteronRot[1]=TMath::Sin(phi)*momentumDeuteron[0]+TMath::Cos(phi)*momentumDeuteron[1];
1017
1018         //invariant mass calculations
1019         fIsCorrectlyAssociated[fItrk] = kFALSE;
1020
1021         //factor for the invariant mass calculation, which only include the pion
1022         e12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1023           +momentumPion[0]*momentumPion[0]
1024           +momentumPion[1]*momentumPion[1]
1025           +momentumPion[2]*momentumPion[2];
1026
1027         r12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1028           +momentumPionRot[0]*momentumPionRot[0]
1029           +momentumPionRot[1]*momentumPionRot[1]
1030           +momentumPion[2]*momentumPion[2];
1031
1032         //factor for the invariant mass calculation, which only include the deuterons
1033         d22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1034           +momentumDeuteron[0]*momentumDeuteron[0]
1035           +momentumDeuteron[1]*momentumDeuteron[1]
1036           +momentumDeuteron[2]*momentumDeuteron[2];
1037         dr22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1038           +momentumDeuteronRot[0]*momentumDeuteronRot[0]
1039           +momentumDeuteronRot[1]*momentumDeuteronRot[1]
1040           +momentumDeuteron[2]*momentumDeuteron[2];
1041
1042         if(rotation == 1){ //signal
1043           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1044                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1045                                                         + 2.*(TMath::Sqrt(e12*d22)
1046                                                               -momentumPion[0]*momentumDeuteron[0]
1047                                                               -momentumPion[1]*momentumDeuteron[1]
1048                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1049
1050           if (fMCtrue){
1051
1052             Int_t labelN = 0;
1053             Int_t labelP = 0;
1054             labelN = trackN->GetLabel();
1055             labelP = trackP->GetLabel();
1056
1057             TParticle *tparticleDaughterN = stack->Particle(TMath::Abs(labelN));
1058             TParticle *tparticleDaughterP = stack->Particle(TMath::Abs(labelP));
1059
1060             Int_t labelMotherN = tparticleDaughterN->GetFirstMother();
1061             Int_t labelMotherP = tparticleDaughterP->GetFirstMother();
1062
1063             TParticle *tparticleMotherN = stack->Particle(TMath::Abs(labelMotherN));
1064             TParticle *tparticleMotherP = stack->Particle(TMath::Abs(labelMotherP));
1065
1066             //LambdaNeuteron
1067             if(tparticleMotherN->GetPdgCode() == fgkPdgCode[kPDGLambdaNeutron] && tparticleMotherP->GetPdgCode() == fgkPdgCode[kPDGLambdaNeutron] && onFlyStatus ==1 && labelMotherN == labelMotherP ){//check mother PDG  and fill the histogramms only for the only V0 finder
1068
1069               //cout << "number of daughters: " << tparticleMotherN->GetNDaughters() << endl;;
1070               if(tparticleMotherN->GetNDaughters() > 2.) continue;
1071
1072               Int_t labelSecondDaughter = tparticleMotherP->GetDaughter(1);
1073               Int_t labelFirstDaughter = labelSecondDaughter-1;
1074
1075               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1076               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1077
1078               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGDeuteron]){//check first daughter PDG
1079
1080                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionMinus]){//check second daughter PDG
1081
1082                   if(invaMassDeuteronPion < 2.02) continue;
1083
1084                   fHistLambdaNeutronPtAso->Fill(transversMomentumMother);
1085                   fHistLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1086                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1087                   fHistLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1088                   //CUTS
1089
1090                   if((dcaV0 < 0.5) &&
1091                      (cosPointing > 0.999) &&
1092                      (decayRadius < 50.0) &&
1093                      (ptotP > 0.2)){//trackP->GetTPCsignal()>180 &&&&  decayRadius > 1.5 && decayRadius< 50
1094
1095                     fHistLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1096                   }
1097                 }//end check second daughter PDG
1098               }//end check first daughter PDG
1099             }//end LambdaNeutron
1100
1101             //Anti-LambdaNeutron
1102             if(tparticleMotherN->GetPdgCode() == fgkPdgCode[kPDGAntiLambdaNeutron] && tparticleMotherP->GetPdgCode() == fgkPdgCode[kPDGAntiLambdaNeutron] && onFlyStatus ==1 && labelMotherN == labelMotherP){//check mother PDG  and fill the histogramms only for the only V0 finder
1103
1104               Int_t labelSecondDaughter = tparticleMotherN->GetDaughter(1);
1105               Int_t labelFirstDaughter = labelSecondDaughter-1;
1106
1107               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1108               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1109
1110               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGAntiDeuteron]) {//check first daughter PDG
1111
1112                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionPlus]) {//check second daughter PDG
1113
1114                   if(invaMassDeuteronPion < 2.02) continue;
1115
1116                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1117                   fHistAntiLambdaNeutronPtAso->Fill(transversMomentumMother);
1118                   fHistAntiLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1119                   fHistAntiLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1120
1121                   //CUTS
1122                   if(dcaV0 < 1. && cosPointing > 0.99 && ptotN > 0.2)//&& trackN->GetTPCsignal()>180 && decayRadius > 1.5 && decayRadius< 50
1123                     {
1124                       fHistAntiLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1125                     }
1126                 }//end check second daughter PDG
1127               }//end check first daughter PDG
1128             }//end Anti-LambdaNeutron
1129           }//end MC
1130         }//end rotation == 1, signal
1131
1132         if(rotation == 2){ // rotation of the pion
1133           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1134                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1135                                                         + 2.*(TMath::Sqrt(r12*d22)
1136                                                               -momentumPionRot[0]*momentumDeuteron[0]
1137                                                               -momentumPionRot[1]*momentumDeuteron[1]
1138                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1139         }//end rotation == 2, rotation of the pion
1140
1141         if(rotation == 3){// Rotation of the deuteron
1142           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1143                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1144                                                         + 2.*(TMath::Sqrt(e12*dr22)
1145                                                               -momentumPion[0]*momentumDeuteronRot[0]
1146                                                               -momentumPion[1]*momentumDeuteronRot[1]
1147                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1148         }//end rotation == 3, rotation of the deuteron
1149
1150         //fill the THnSparse and the tree variables
1151
1152         //tree variables which are independent of the particle-species
1153         fV0finder[fItrk]         = onFlyStatus;
1154         fkMB[fItrk]              = fTriggerFired[0];
1155         fkCentral[fItrk]         = fTriggerFired[1];
1156         fkSemiCentral[fItrk]     = fTriggerFired[2];
1157         fkEMCEJE[fItrk]          = fTriggerFired[3];
1158         fkEMCEGA[fItrk]          = fTriggerFired[4];
1159
1160         fPtotN[fItrk]            = ptotN;
1161         fPtotP[fItrk]            = ptotP;
1162         fMotherPt[fItrk]         = transversMomentumMother;
1163
1164         fdEdxN[fItrk]            = trackN->GetTPCsignal();
1165         fdEdxP[fItrk]            = trackP->GetTPCsignal();
1166         fSignN[fItrk]            = trackN->Charge();
1167         fSignP[fItrk]            = trackP->Charge();
1168
1169         fDCAv0[fItrk]            = dcaV0;
1170         fCosinePAv0[fItrk]       = cosPointing;
1171         fDecayRadiusTree[fItrk]  = decayRadius;
1172
1173         fAmenterosAlphaTree[fItrk] = alpha;
1174         fAmenterosQtTree[fItrk] = qt;
1175         fRotationTree[fItrk] = rotation;
1176
1177
1178         if (isDeuteron[0] == kTRUE)  //pos deuteron
1179           {
1180             fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1181             fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1182
1183             fItrk++;
1184
1185             if(invaMassDeuteronPion < 2.1)
1186               {
1187                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotP*(trackP->Charge()), trackP->GetTPCsignal(),9);
1188                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotN*(trackN->Charge()), trackN->GetTPCsignal(),10);
1189               }
1190             fHistArmenterosPodolanskiDeuteronPion->Fill(alpha*(-1),qt);
1191           }
1192
1193         if (isDeuteron[1] == kTRUE)
1194           {
1195             fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1196             fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1197
1198             fItrk++;
1199
1200             if(invaMassDeuteronPion < 2.1)
1201               {
1202                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),7);
1203                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),8);
1204               }
1205             fHistArmenterosPodolanskiAntiDeuteronPion->Fill(alpha,qt);
1206           }
1207
1208       }//end rotation loop
1209
1210  } //end of v0 loop
1211
1212   //fill the tree
1213   fTreeV0->Fill();
1214
1215   // Post output data.
1216   PostData(1, fOutputContainer);
1217   PostData(2, fTreeV0);
1218 }
1219
1220 //________________________________________________________________________
1221 void AliAnalysisTaskLambdaNAOD::Terminate(const Option_t *)
1222 {
1223   // Draw result to the screen
1224   // Called once at the end of the query
1225
1226   //get output data and darw 'fHistPt'
1227   if (!GetOutputData(0)) return;
1228   //TH1F *hist=(TH1F*)(((TObjArray*)GetOutputData(0))->FindObject("fHistPt"));
1229   //if (hist) hist->Draw();
1230 }
1231
1232 //_____________________________________________________________________________
1233 Int_t AliAnalysisTaskLambdaNAOD::Initialize() {
1234
1235
1236   // -- Reset Event
1237   // ----------------
1238   ResetEvent();
1239
1240   return 0;
1241   }
1242 //________________________________________________________________________
1243 Int_t AliAnalysisTaskLambdaNAOD::SetupEvent() {
1244   // Setup Reading of event
1245
1246   ResetEvent();
1247
1248   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1249   // -- Get Trigger
1250   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1251
1252   //Bool_t isTriggered = IsTriggered();
1253   IsTriggered();
1254   return 0;
1255 }
1256 //_____________________________________________________________________________
1257 void AliAnalysisTaskLambdaNAOD::ResetEvent() {
1258
1259   // Reset event
1260
1261   // -- Reset QA
1262   for (Int_t ii = 0; ii < fNTriggers; ++ii)
1263     fTriggerFired[ii] = kFALSE;
1264
1265   return;
1266 }
1267 //_______________________________________________________________________
1268 void AliAnalysisTaskLambdaNAOD::BinLogAxis(const THnSparse *h, Int_t axisNumber) {
1269   //
1270   // Method for the correct logarithmic binning of histograms
1271   //
1272   TAxis *axis = h->GetAxis(axisNumber);
1273   int bins = axis->GetNbins();
1274
1275   Double_t from = axis->GetXmin();
1276   Double_t to = axis->GetXmax();
1277   Double_t *newBins = new Double_t[bins + 1];
1278
1279   newBins[0] = from;
1280   Double_t factor = pow(to/from, 1./bins);
1281
1282   for (int i = 1; i <= bins; i++) {
1283    newBins[i] = factor * newBins[i-1];
1284   }
1285   axis->Set(bins, newBins);
1286   delete [] newBins;
1287
1288   }
1289 //________________________________________________________________________
1290 Bool_t AliAnalysisTaskLambdaNAOD::IsTriggered() {
1291   // Check if Event is triggered and fill Trigger Histogram
1292
1293   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
1294   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
1295   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
1296   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
1297   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
1298
1299   Bool_t isTriggered = kFALSE;
1300
1301   for (Int_t ii=0; ii<fNTriggers; ++ii) {
1302     if(fTriggerFired[ii]) {
1303       isTriggered = kTRUE;
1304       fHistTriggerStat->Fill(ii);
1305     }
1306   }
1307
1308   return isTriggered;
1309   }
1310 //______________________________________________________________________________
1311 Bool_t AliAnalysisTaskLambdaNAOD::DeuteronPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isDeuteron[3]) {
1312
1313      //define the arrays for the Bethe-Bloch-Parameters
1314       Double_t parDeuteron[5] = {0,0,0,0,0};
1315
1316       if(runNumber < 166500) //LHC10h
1317         {
1318           parDeuteron[0] = 1.45802; // ALEPH parameters for deuterons (pass2)
1319           parDeuteron[1] = 27.4992;
1320           parDeuteron[2] = 4.00313e-15;
1321           parDeuteron[3] = 2.48485;
1322           parDeuteron[4] = 8.31768;
1323         }
1324
1325       if(runNumber > 166500) //LHC11h
1326         {
1327           parDeuteron[0] = 1.11243; // ALEPH parameters for deuterons (pass2)
1328           parDeuteron[1] = 26.1144;
1329           parDeuteron[2] = 4.00313e-15;
1330           parDeuteron[3] = 2.72969 ;
1331           parDeuteron[4] = 9.15038;
1332         }
1333
1334
1335       //define expected signals for the various species
1336       Double_t expSignalDeuteronN = 0;
1337       Double_t expSignalDeuteronP = 0;
1338
1339       isDeuteron[0] = kFALSE;
1340       isDeuteron[1] = kFALSE;
1341       isDeuteron[2] = kFALSE;
1342
1343       //for data
1344       if(!fMCtrue){
1345
1346         expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1347         expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1348
1349         if(//trackP->GetTPCsignal() >= 110 && ///??????????????????????????????????
1350            //trackP->GetTPCsignal() < 1200 &&
1351            (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1352            ptotP > 0.2 ){
1353
1354           if(trackP->Charge() >0)               isDeuteron[0] = kTRUE; //pos deuteron
1355           if(trackP->Charge() <0)               isDeuteron[1] = kTRUE; //neg deuteron
1356         }
1357
1358         if(//trackN->GetTPCsignal() >= 110 && ///??????????????????????????????????
1359            //trackN->GetTPCsignal() < 1200 &&
1360            (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1361            ptotN > 0.2){
1362
1363           isDeuteron[2] = kTRUE;
1364
1365           if(trackN->Charge() >0)        isDeuteron[0] = kTRUE; //pos deuteron
1366           if(trackN->Charge() <0)        isDeuteron[1] = kTRUE; //neg deuteron
1367         }
1368       }
1369
1370       //for MC
1371       if(fMCtrue)
1372         {
1373           if(runNumber < 166500) //2010
1374             {
1375               expSignalDeuteronN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1376               expSignalDeuteronP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1377             }
1378           if(runNumber > 166500) //2011
1379             {
1380               expSignalDeuteronN = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1381               expSignalDeuteronP = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1382             }
1383
1384           if(trackP->GetTPCsignal() < 1200 &&
1385              (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1386              ptotP > 0.2 )
1387             {
1388               if(trackP->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1389               if(trackP->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1390             }
1391
1392           if(trackN->GetTPCsignal() < 1200 &&
1393              (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1394              ptotN > 0.2 )
1395             {
1396               isDeuteron[2] = kTRUE;
1397
1398               if(trackN->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1399               if(trackN->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1400             }
1401
1402         }
1403
1404
1405
1406      return isDeuteron;
1407
1408 }
1409 //______________________________________________________________________________
1410 Bool_t AliAnalysisTaskLambdaNAOD::PionPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isPion[2]) {
1411
1412   //Pion PID via specific energy loss in the TPC
1413   //define the array for the Bethe-Bloch-Parameters (only necessary for local MC)
1414   Double_t parPion[5]  = {0,0,0,0,0};
1415
1416     if(runNumber < 166500) { //LHC10h
1417         parPion[0]   = 1.45802; // ALEPH parameters for pions (pass2)
1418         parPion[1]   = 27.4992;
1419         parPion[2]   = 4.00313e-15;
1420         parPion[3]   = 2.48485;
1421         parPion[4]   = 8.31768;
1422     }
1423
1424     if(runNumber > 166500) {  //LHC11h
1425         parPion[0]   = 1.11243; // ALEPH parameters for pions (pass2)
1426         parPion[1]   = 26.1144;
1427         parPion[2]   = 4.00313e-15;
1428         parPion[3]   = 2.72969;
1429         parPion[4]   = 9.15038;
1430     }
1431
1432
1433     Double_t expSignalPionP = 0;
1434     Double_t expSignalPionN = 0;
1435
1436     //for MC
1437     if(fMCtrue){
1438       if(runNumber < 166500){ //2010
1439         expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1440         expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1441       }
1442       if(runNumber > 166500){ //2011
1443         expSignalPionP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1444         expSignalPionN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1445       }
1446
1447     }
1448
1449     isPion[0] = kFALSE;
1450     isPion[1] = kFALSE;
1451     //data
1452     if(!fMCtrue){
1453       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion))<3){
1454
1455         if(trackP->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1456         if(trackP->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1457      }
1458       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion))<3){
1459
1460         if(trackN->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1461         if(trackN->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1462       }
1463     }
1464
1465     //MC
1466     if(fMCtrue){
1467       if(TMath::Abs(trackP->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.2
1468          && ptotP>0.00001){
1469         if(trackP->Charge()>0) isPion[0] = kTRUE; //pos pion
1470         if(trackP->Charge()<0) isPion[1] = kTRUE; //neg pion
1471       }
1472      if(TMath::Abs(trackN->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.2
1473         && ptotN>0.00001){
1474         if(trackN->Charge()>0) isPion[0] = kTRUE; //pos pion
1475         if(trackN->Charge()<0) isPion[1] = kTRUE; //neg pion
1476       }
1477     }
1478
1479     return isPion;
1480
1481 }
1482 //______________________________________________________________________________
1483 Bool_t AliAnalysisTaskLambdaNAOD::TrackCuts(AliVTrack *track, Bool_t testTrackCuts) {
1484
1485   testTrackCuts = kFALSE;
1486
1487   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1488   //if(!esdtrack) testTrackCuts = kFALSE;
1489   if (fESDtrackCutsV0->AcceptTrack(esdtrack)) testTrackCuts = kTRUE;
1490 //if( testTrackCuts == kTRUE) cout <<   "testTrackCuts im test: " << testTrackCuts << endl;
1491
1492
1493   return testTrackCuts;
1494 }
1495 //______________________________________________________________________________
1496 /*Bool_t AliAnalysisTaskLambdaNAOD::FilterBit(AliVTrack *track, Bool_t testFilterBit) {
1497
1498   testFilterBit = kFALSE;
1499
1500   AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1501   //if(!aodtrack) testFilterBit = kFALSE;
1502   if (aodtrack->TestFilterBit(7)) testFilterBit = kTRUE;
1503
1504   //if(testFilterBit == kTRUE) cout <<   "testFilterBit im test: " << testFilterBit<< endl;
1505
1506   return testFilterBit;
1507   }*/
1508 //______________________________________________________________________
1509 void AliAnalysisTaskLambdaNAOD::MCGenerated(AliStack* stack)
1510 {
1511
1512   // Monte Carlo for genenerated particles
1513   if (fMCtrue) //MC loop
1514     {
1515
1516       Int_t stackN = 0;
1517
1518       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1519         {
1520
1521           const TParticle *tparticleMother = stack->Particle(stackN);
1522           Long_t pdgCodeMother = tparticleMother->GetPdgCode();
1523
1524           //check which particle the mother is
1525
1526           //LambdaNeutron
1527           if(pdgCodeMother == fgkPdgCode[kPDGLambdaNeutron]) //check mother PDG
1528             {
1529               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGDeuteron],fgkPdgCode[kPDGPionMinus],fgkMass[kMassDeuteron],fgkMass[kMassPion]);
1530             }
1531
1532           //Anti-LambdaNeutron
1533           if(pdgCodeMother == fgkPdgCode[kPDGAntiLambdaNeutron]) //check mother PDG
1534             {
1535               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGAntiDeuteron],fgkPdgCode[kPDGPionPlus],fgkMass[kMassDeuteron],fgkMass[kMassPion]);
1536             }
1537
1538
1539         }//end loop over stack
1540
1541
1542     }//end MC
1543 }
1544 //_____________________________________________
1545 void AliAnalysisTaskLambdaNAOD::MCTwoBodyDecay(AliStack* stack, const TParticle *tparticleMother, Long_t PDGMother, Long_t PDGFirstDaughter, Long_t PDGSecondDaughter, Double_t massFirstDaughter, Double_t massSecondDaughter){ //function that calculates the invariant mass and the transverse momentum for MC
1546
1547   Double_t momentumFirstDaughterGen[3]={0,0,0};
1548   Double_t momentumSecondDaughterGen[3]={0,0,0};
1549
1550   Double_t energyFirstDaughterGen = 0;
1551   Double_t energySecondDaughterGen = 0;
1552
1553   Double_t transversMomentumMotherGen = 0;
1554   Double_t longitudinalMomentumMotherGen = 0;
1555   Double_t totalEnergyMotherGen = 0;
1556
1557   Double_t rapidityGen = 2;
1558
1559   //Int_t labelFirstDaughter = tparticleMother->GetDaughter(0);
1560   Int_t labelSecondDaughter = tparticleMother->GetDaughter(1);
1561   Int_t labelFirstDaughter = labelSecondDaughter -1;
1562
1563   TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1564   TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1565
1566   if(tparticleFirstDaughter->GetPdgCode() == PDGFirstDaughter) //check first daughter PDG
1567     {
1568       if(tparticleSecondDaughter->GetPdgCode() == PDGSecondDaughter) //check second daughter PDG
1569         {
1570
1571           momentumFirstDaughterGen[0] = tparticleFirstDaughter->Px();
1572           momentumFirstDaughterGen[1] = tparticleFirstDaughter->Py();
1573           momentumFirstDaughterGen[2] = tparticleFirstDaughter->Pz();
1574
1575           momentumSecondDaughterGen[0] = tparticleSecondDaughter->Px();
1576           momentumSecondDaughterGen[1] = tparticleSecondDaughter->Py();
1577           momentumSecondDaughterGen[2] = tparticleSecondDaughter->Pz();
1578
1579           energyFirstDaughterGen  = tparticleFirstDaughter->Energy();
1580           energySecondDaughterGen = tparticleSecondDaughter->Energy();
1581
1582           transversMomentumMotherGen=TMath::Sqrt((momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])*(momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])+(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1])*(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1]));
1583
1584           //rapidity cut +- 1
1585           //longitudinal momentum of mother
1586           longitudinalMomentumMotherGen = TMath::Sqrt((momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2])*(momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2]));
1587
1588           //total energy of mother
1589           totalEnergyMotherGen = energyFirstDaughterGen + energySecondDaughterGen;
1590
1591           //rapidity
1592           rapidityGen = 0.5 * TMath::Log( (totalEnergyMotherGen+longitudinalMomentumMotherGen)/(totalEnergyMotherGen-longitudinalMomentumMotherGen));
1593
1594           if(rapidityGen > 1.0 || rapidityGen < -1 ) return;
1595
1596           //cout << "decay lenght first daughter :" << tparticleFirstDaughter->Rho()  << " decay lenght second daughter :" << tparticleSecondDaughter->Rho()<< endl;
1597           //fHistLambdaNeutronDecayLengthGen->Fill(tparticleFirstDaughter->Rho());
1598
1599           //calculate the invariant mass
1600           Double_t firstDaughterPart = 0;
1601           Double_t secondDaughterPart = 0;
1602           Double_t invaMass = 0;
1603
1604           firstDaughterPart = massFirstDaughter*massFirstDaughter
1605             +momentumFirstDaughterGen[0]*momentumFirstDaughterGen[0]
1606             +momentumFirstDaughterGen[1]*momentumFirstDaughterGen[1]
1607             +momentumFirstDaughterGen[2]*momentumFirstDaughterGen[2];
1608           secondDaughterPart = massSecondDaughter*massSecondDaughter
1609             +momentumSecondDaughterGen[0]*momentumSecondDaughterGen[0]
1610             +momentumSecondDaughterGen[1]*momentumSecondDaughterGen[1]
1611             +momentumSecondDaughterGen[2]*momentumSecondDaughterGen[2];
1612
1613           invaMass = TMath::Sqrt(TMath::Max(massFirstDaughter*massFirstDaughter
1614                                                         +massSecondDaughter*massSecondDaughter
1615                                                         + 2.*(TMath::Sqrt(firstDaughterPart*secondDaughterPart)
1616                                                               -momentumFirstDaughterGen[0]*momentumSecondDaughterGen[0]
1617                                                               -momentumFirstDaughterGen[1]*momentumSecondDaughterGen[1]
1618                                                               -momentumFirstDaughterGen[2]*momentumSecondDaughterGen[2]), 0.));
1619
1620 #if 0
1621
1622                                                               switch(PDGMother) {
1623                                                               case fgkPdgCode[kPDGLambdaNeutron] :
1624                                                                fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1625                                                                fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1626                                                                break;
1627                                                               case fgkPdgCode[kPDGAntiLambdaNeutron] :
1628                                                                fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1629                                                                fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1630                                                                break;
1631                                                               default :
1632                                                               printf("should not happen!!!! \n");
1633                                                               }
1634
1635
1636 #else
1637           //LambdaNeutron
1638           if(PDGMother == fgkPdgCode[kPDGLambdaNeutron])
1639             {
1640               fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1641               fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1642               fHistLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1643             }
1644           //Anti-LambdaNeutron
1645           if(PDGMother == fgkPdgCode[kPDGAntiLambdaNeutron])
1646             {
1647               fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1648               fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1649               fHistAntiLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1650             }
1651
1652 #endif
1653         }//end of check second daughter PDG
1654     }//end of check first daughter PDG
1655
1656 }
1657 //_____________________________________________
1658 Double_t AliAnalysisTaskLambdaNAOD::MomentumInnerParam(AliVTrack *track, Double_t ptot){ //function to get the momentum at the inner wall of the TPC
1659
1660   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1661   ptot = esdtrack->GetInnerParam()->GetP();
1662
1663   return ptot;
1664 }
1665 //________________________________________________________________________
1666 /*void AliAnalysisTaskLambdaNAOD::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
1667
1668      // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
1669      if (!kfParticle) return;
1670      Double_t dx = 0.;
1671      Double_t dy = 0.;
1672      Double_t dz = 0.;
1673
1674      if (ev){
1675        dx = ev->GetPrimaryVertex()->GetX()-0.;
1676        dy = ev->GetPrimaryVertex()->GetY()-0.;
1677        dz = ev->GetPrimaryVertex()->GetZ()-0.;
1678      }
1679
1680      kfParticle->X() = kfParticle->GetX() - dx;
1681      kfParticle->Y() = kfParticle->GetY() - dy;
1682      kfParticle->Z() = kfParticle->GetZ() - dz;
1683
1684
1685      // Rotate the kf particle
1686      Double_t c = cos(angle);
1687      Double_t s = sin(angle);
1688
1689      Double_t mA[8][ 8];
1690      for( Int_t i=0; i<8; i++ ){
1691        for( Int_t j=0; j<8; j++){
1692          mA[i][j] = 0;
1693        }
1694      }
1695      for( int i=0; i<8; i++ ){
1696        mA[i][i] = 1;
1697      }
1698      mA[0][0] =  c;  mA[0][1] = s;
1699      mA[1][0] = -s;  mA[1][1] = c;
1700      mA[3][3] =  c;  mA[3][4] = s;
1701      mA[4][3] = -s;  mA[4][4] = c;
1702
1703      Double_t mAC[8][8];
1704      Double_t mAp[8];
1705
1706      for( Int_t i=0; i<8; i++ ){
1707        mAp[i] = 0;
1708        for( Int_t k=0; k<8; k++){
1709          mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
1710        }
1711      }
1712
1713      for( Int_t i=0; i<8; i++){
1714        kfParticle->Parameter(i) = mAp[i];
1715      }
1716
1717      for( Int_t i=0; i<8; i++ ){
1718        for( Int_t j=0; j<8; j++ ){
1719          mAC[i][j] = 0;
1720          for( Int_t k=0; k<8; k++ ){
1721            mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
1722          }
1723        }
1724      }
1725
1726      for( Int_t i=0; i<8; i++ ){
1727        for( Int_t j=0; j<=i; j++ ){
1728          Double_t xx = 0;
1729          for( Int_t k=0; k<8; k++){
1730            xx+= mAC[i][k]*mA[j][k];
1731          }
1732          kfParticle->Covariance(i,j) = xx;
1733        }
1734      }
1735
1736      kfParticle->X() = kfParticle->GetX() + dx;
1737      kfParticle->Y() = kfParticle->GetY() + dy;
1738      kfParticle->Z() = kfParticle->GetZ() + dz;
1739
1740      }*/