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   fHistCentralityClass10(0),
  99   fHistCentralityPercentile(0),
 100   fHistTriggerStat(0),
 101   fHistTriggerStatAfterEventSelection(0),
 102   fHistLambdaNeutronPtGen(0),
 103   fHistLambdaNeutronPtGenMinBias(0),
 104   fHistLambdaNeutronPtGenCentral(0),
 105   fHistLambdaNeutronPtGenSemiCentral(0),
 106   fHistAntiLambdaNeutronPtGen(0),
 107   fHistAntiLambdaNeutronPtGenMinBias(0),
 108   fHistAntiLambdaNeutronPtGenCentral(0),
 109   fHistAntiLambdaNeutronPtGenSemiCentral(0),
 110   fHistLambdaNeutronInvaMassGen(0),
 111   fHistAntiLambdaNeutronInvaMassGen(0),
 112   fHistLambdaNeutronDecayLengthGen(0),
 113   fHistAntiLambdaNeutronDecayLengthGen(0),
 114   fHistLambdaNeutronPtAso(0),
 115   fHistLambdaNeutronPtAsoCuts(0),
 116   fHistAntiLambdaNeutronPtAso(0),
 117   fHistAntiLambdaNeutronPtAsoCuts(0),
 118   fHistLambdaNeutronInvaMassAso(0),
 119   fHistAntiLambdaNeutronInvaMassAso(0),
 120   fHistLambdaNeutronDecayLengthAso(0),
 121   fHistAntiLambdaNeutronDecayLengthAso(0),
 122   fof(0),
 123   fHistArmenterosPodolanskiDeuteronPion(0),
 124   fHistArmenterosPodolanskiAntiDeuteronPion(0),
 125   fHistDeDxQA(0),
 126 //fHistdEdxMC(0),
 127   fNTriggers(5),
 128   fMCtrue(0),
 129   fOnlyQA(0),
 130   fTriggerFired(),
 131 //fV0object(NULL),
 132   fItrk(0),
 133   fOutputContainer(NULL)
 134 {
 135   // default Constructor
 136
 137   // Define input and output slots here
 138 }
 139
 140 //________________________________________________________________________
 141 AliAnalysisTaskLambdaNAOD::AliAnalysisTaskLambdaNAOD(const char *name)
 142   : AliAnalysisTaskSE(name),
 143     fAnalysisType("ESD"),
 144     fEventHandler(0),
 145     fESDtrackCutsV0(0),
 146     fESDpid(0),
 147     fPIDResponse(0),
 148     fTreeV0(0),
 149     fHistCentralityClass10(0),
 150     fHistCentralityPercentile(0),
 151     fHistTriggerStat(0),
 152     fHistTriggerStatAfterEventSelection(0),
 153     fHistLambdaNeutronPtGen(0),
 154     fHistLambdaNeutronPtGenMinBias(0),
 155     fHistLambdaNeutronPtGenCentral(0),
 156     fHistLambdaNeutronPtGenSemiCentral(0),
 157     fHistAntiLambdaNeutronPtGen(0),
 158     fHistAntiLambdaNeutronPtGenMinBias(0),
 159     fHistAntiLambdaNeutronPtGenCentral(0),
 160     fHistAntiLambdaNeutronPtGenSemiCentral(0),
 161     fHistLambdaNeutronInvaMassGen(0),
 162     fHistAntiLambdaNeutronInvaMassGen(0),
 163     fHistLambdaNeutronDecayLengthGen(0),
 164     fHistAntiLambdaNeutronDecayLengthGen(0),
 165     fHistLambdaNeutronPtAso(0),
 166     fHistLambdaNeutronPtAsoCuts(0),
 167     fHistAntiLambdaNeutronPtAso(0),
 168     fHistAntiLambdaNeutronPtAsoCuts(0),
 169     fHistLambdaNeutronInvaMassAso(0),
 170     fHistAntiLambdaNeutronInvaMassAso(0),
 171     fHistLambdaNeutronDecayLengthAso(0),
 172     fHistAntiLambdaNeutronDecayLengthAso(0),
 173     fof(0),
 174     fHistArmenterosPodolanskiDeuteronPion(0),
 175     fHistArmenterosPodolanskiAntiDeuteronPion(0),
 176     fHistDeDxQA(0),
 177     //fHistdEdxMC(0),
 178     fNTriggers(5),
 179     fMCtrue(0),
 180     fOnlyQA(0),
 181     fTriggerFired(),
 182     //fV0object(NULL),
 183     fItrk(0),
 184     fOutputContainer(NULL)
 185 {
 186   // Constructor
 187
 188   // Define input and output slots here
 189   // Input slot #0 works with a TChain
 190   DefineInput(0, TChain::Class());
 191   // Output slot #0 writes into a TH1 container
 192   DefineOutput(1, TObjArray::Class());
 193   DefineOutput(2, TTree::Class());
 194
 195   //ESD Track Cuts for v0's
 196   fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
 197   fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
 198   fESDtrackCutsV0->SetMinNClustersTPC(60);
 199   fESDtrackCutsV0->SetMaxChi2PerClusterTPC(6); //set to 6 for sytematics
 200   //fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
 201   fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
 202   //  fESDtrackCutsV0->SetMinNClustersITS(1); // to be tested if this cut is not too strong
 203   fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
 204
 205
 206   fMCtrue = kTRUE;
 207   fOnlyQA = kTRUE;
 208
 209 }
 210
 211 //____________________________________________________________
 212 AliAnalysisTaskLambdaNAOD::~AliAnalysisTaskLambdaNAOD(){
 213   //
 214   // Destructor
 215   //
 216   //if (fESD) delete fESD;
 217   if (fESDtrackCutsV0) delete fESDtrackCutsV0;
 218
 219 }
 220
 221 //____________________________________________________________
 222 const Int_t AliAnalysisTaskLambdaNAOD::fgkPdgCode[] = {
 223   211,                //PionPlus
 224   -211,               //PionMinus
 225   2212,               //Proton
 226   -2212,              //Anti-Proton
 227   1000010020,         //Deuteron
 228   -1000010020,        //Anti-Deuteron
 229   1000020030,         //Helium3
 230   -1000020030,        //Anti-Helium3
 231   3122,               //Lambda
 232   -3122,              //Anti-Lambda
 233   1010000020,         //LambdaNeutron
 234   -1010000020,        //Anti-Lambda-Neutron
 235   1010010030,         //Hypertriton
 236   -1010010030         //Anti-Hypertriton
 237 };
 238
 239 //____________________________________________________________
 240 const Double_t AliAnalysisTaskLambdaNAOD::fgkMass[] = {
 241   0.13957,           //Pion
 242   0.93827,           //Proton
 243   1.875612,          //Deuteron
 244   2.808921,          //Triton
 245   2.80941            //Helium3 Quelle: Wolfram Alpha
 246 };
 247
 248 //________________________________________________________________________
 249 void AliAnalysisTaskLambdaNAOD::UserCreateOutputObjects(){
 250
 251   // New PID object
 252   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 253   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 254   fPIDResponse = inputHandler->GetPIDResponse();
 255
 256   // Create histograms for MC
 257   //generated histogramms
 258   fHistLambdaNeutronPtGen = new TH1F("fHistLambdaNeutronPtGen", "Generated  Lambdan", 201, 0., 10.1);
 259   fHistLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
 260   fHistLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 261
 262   fHistLambdaNeutronPtGenMinBias = new TH1F("fHistLambdaNeutronPtGenMinBias", "Generated  Lambdan", 201, 0., 10.1);
 263   fHistLambdaNeutronPtGenMinBias->GetYaxis()->SetTitle("Counts");
 264   fHistLambdaNeutronPtGenMinBias->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 265
 266   fHistLambdaNeutronPtGenCentral = new TH1F("fHistLambdaNeutronPtGenCentral", "Generated  Lambdan", 201, 0., 10.1);
 267   fHistLambdaNeutronPtGenCentral->GetYaxis()->SetTitle("Counts");
 268   fHistLambdaNeutronPtGenCentral->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 269
 270   fHistLambdaNeutronPtGenSemiCentral = new TH1F("fHistLambdaNeutronPtGenSemiCentral", "Generated  Lambdan", 201, 0., 10.1);
 271   fHistLambdaNeutronPtGenSemiCentral->GetYaxis()->SetTitle("Counts");
 272   fHistLambdaNeutronPtGenSemiCentral->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 273
 274   fHistAntiLambdaNeutronPtGen = new TH1F("fHistAntiLambdaNeutronPtGen", "Generated  #bar{#Lambdan}", 201, 0., 10.1);
 275   fHistAntiLambdaNeutronPtGen->GetYaxis()->SetTitle("Counts");
 276   fHistAntiLambdaNeutronPtGen->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 277
 278   fHistAntiLambdaNeutronPtGenMinBias = new TH1F("fHistAntiLambdaNeutronPtGenMinBias", "Generated  #bar{#Lambdan}", 201, 0., 10.1);
 279   fHistAntiLambdaNeutronPtGenMinBias->GetYaxis()->SetTitle("Counts");
 280   fHistAntiLambdaNeutronPtGenMinBias->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 281
 282   fHistAntiLambdaNeutronPtGenCentral = new TH1F("fHistAntiLambdaNeutronPtGenCentral", "Generated  #bar{#Lambdan}", 201, 0., 10.1);
 283   fHistAntiLambdaNeutronPtGenCentral->GetYaxis()->SetTitle("Counts");
 284   fHistAntiLambdaNeutronPtGenCentral->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 285
 286   fHistAntiLambdaNeutronPtGenSemiCentral = new TH1F("fHistAntiLambdaNeutronPtGenSemiCentral", "Generated  #bar{#Lambdan}", 201, 0., 10.1);
 287   fHistAntiLambdaNeutronPtGenSemiCentral->GetYaxis()->SetTitle("Counts");
 288   fHistAntiLambdaNeutronPtGenSemiCentral->GetXaxis()->SetTitle("#it{p}_{T}   (GeV/#it{c})");
 289
 290   fHistLambdaNeutronInvaMassGen = new TH1F("fHistLambdaNeutronInvaMassGen", "Generated #Lambdan ", 100, 2.0, 2.1);
 291   fHistLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
 292   fHistLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
 293
 294   fHistAntiLambdaNeutronInvaMassGen = new TH1F("fHistAntiLambdaNeutronInvaMassGen", "Generated  #bar{#Lambdan}", 100, 2.0, 2.1);
 295   fHistAntiLambdaNeutronInvaMassGen->GetYaxis()->SetTitle("Counts");
 296   fHistAntiLambdaNeutronInvaMassGen->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+}) (GeV/#it{c}^{2})");
 297
 298   fHistLambdaNeutronDecayLengthGen = new TH1F("fHistLambdaNeutronDecayLengthGen", "Generated  #Lambdan", 401, 0., 400.1);
 299   fHistLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
 300   fHistLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 301
 302   fHistAntiLambdaNeutronDecayLengthGen = new TH1F("fHistAntiLambdaNeutronDecayLengthGen", "Generated #bar{#Lambdan}", 401, 0., 400.1);
 303   fHistAntiLambdaNeutronDecayLengthGen->GetYaxis()->SetTitle("Counts");
 304   fHistAntiLambdaNeutronDecayLengthGen->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 305
 306   //assoziated (reconstracted) histogramms
 307   fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated   #Lambdan", 201, 0., 10.1);
 308   fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 309   fHistLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 310
 311   fHistLambdaNeutronPtAsoCuts = new TH1F("fHistLambdaNeutronPtAsoCuts", "Associated   #Lambdan Cuts", 201, 0., 10.1);
 312   fHistLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
 313   fHistLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 314
 315   fHistAntiLambdaNeutronPtAsoCuts = new TH1F("fHistAntiLambdaNeutronPtAsoCuts", " associated  #bar{#Lambdan} Cuts", 201, 0., 10.1);
 316   fHistAntiLambdaNeutronPtAsoCuts->GetYaxis()->SetTitle("Counts");
 317   fHistAntiLambdaNeutronPtAsoCuts->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 318
 319   fHistAntiLambdaNeutronPtAso = new TH1F("fHistAntiLambdaNeutronPtAso", " associated  #bar{#Lambdan}", 201, 0., 10.1);
 320   fHistAntiLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 321   fHistAntiLambdaNeutronPtAso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 322
 323   fHistLambdaNeutronInvaMassAso = new TH1F("fHistLambdaNeutronInvaMassAso", "Associated  #Lambdan", 100, 2.0, 2.1);
 324   fHistLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
 325   fHistLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (d #pi^{-}) (GeV/#it{c}^{2})");
 326
 327   fHistAntiLambdaNeutronInvaMassAso = new TH1F("fHistAntiLambdaNeutronInvaMassAso", " Associated  #bar{#Lambdan}", 100, 2.0, 2.1);
 328   fHistAntiLambdaNeutronInvaMassAso->GetYaxis()->SetTitle("Counts");
 329   fHistAntiLambdaNeutronInvaMassAso->GetXaxis()->SetTitle("Invariant mass (#bar{d} #pi^{+})  (GeV/#it{c}^{2})");
 330
 331   fHistLambdaNeutronDecayLengthAso = new TH1F("fHistLambdaNeutronDecayLengthAso", "Associated  #Lambdan", 401, 0., 400.1);
 332   fHistLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
 333   fHistLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 334
 335   fHistAntiLambdaNeutronDecayLengthAso = new TH1F("fHistAntiLambdaNeutronDecayLengthAso", "Associated #bar{#Lambdan}", 401, 0., 400.1);
 336   fHistAntiLambdaNeutronDecayLengthAso->GetYaxis()->SetTitle("Counts");
 337   fHistAntiLambdaNeutronDecayLengthAso->GetXaxis()->SetTitle("#it{decay length}  (cm)");
 338
 339   //fHistLambdaNeutronPtAso = new TH1F("fHistLambdaNeutronPtAso", "Associated   #Lambdan", 201, 0., 10.1);
 340   //fHistLambdaNeutronPtAso->GetYaxis()->SetTitle("Counts");
 341   //Fhistlambdaneutronptaso->GetXaxis()->SetTitle("#it{p}_{T}  (GeV/#it{c})");
 342
 343   //Tree
 344   //------------ Tree and branch definitions ----------------//
 345   fTreeV0 = new TTree("tree", "fTreeV0");
 346   fTreeV0->Branch("fItrk", &fItrk, "fItrk/I");
 347
 348   //fTreeV0->Branch("fV0object",&fV0object,"fV0object[fItrk]");
 349   fTreeV0->Branch("fV0finder",fV0finder,"fV0finder[fItrk]/I");
 350   fTreeV0->Branch("fkMB",fkMB,"fkMB[fItrk]/I");
 351   fTreeV0->Branch("fkCentral",fkCentral,"fkCentral[fItrk]/I");
 352   fTreeV0->Branch("fkSemiCentral",fkSemiCentral,"fkSemiCentral[fItrk]/I");
 353   fTreeV0->Branch("fkEMCEJE",fkEMCEJE,"fkEMCEJE[fItrk]/I");
 354   fTreeV0->Branch("fkEMCEGA",fkEMCEGA,"fkEMCEGA[fItrk]/I");
 355
 356   fTreeV0->Branch("fCentralityClass10",fCentralityClass10,"fCentralityClass10[fItrk]/I");
 357   fTreeV0->Branch("fCentralityPercentile",fCentralityPercentile,"fCentralityPercentile[fItrk]/I");
 358   fTreeV0->Branch("fMultiplicity",fMultiplicity,"fMultipicity[fItrk]/I");
 359   fTreeV0->Branch("fRefMultiplicity",fRefMultiplicity,"fRefMultipicity[fItrk]/I");
 360
 361   fTreeV0->Branch("fPtotN",fPtotN,"fPtotN[fItrk]/D");
 362   fTreeV0->Branch("fPtotP",fPtotP,"fPtotP[fItrk]/D");
 363   fTreeV0->Branch("fMotherPt",fMotherPt,"fMotherPt[fItrk]/D");
 364   fTreeV0->Branch("fdEdxN",fdEdxN,"fdEdxN[fItrk]/D");
 365   fTreeV0->Branch("fdEdxP",fdEdxP,"fdEdxP[fItrk]/D");
 366   fTreeV0->Branch("fSignN",fSignN,"fSignN[fItrk]/D");
 367   fTreeV0->Branch("fSignP",fSignP,"fSignP[fItrk]/D");
 368
 369   fTreeV0->Branch("fSigmadEdxPionPos",fSigmadEdxPionPos,"fSigmadEdxPionPos[fItrk]/D");
 370   fTreeV0->Branch("fSigmadEdxPionNeg",fSigmadEdxPionNeg,"fSigmadEdxPionNeg[fItrk]/D");
 371
 372   fTreeV0->Branch("fDCAv0",fDCAv0,"fDCAv0[fItrk]/F"); //Dca v0 Daughters
 373   fTreeV0->Branch("fCosinePAv0",fCosinePAv0,"fCosinePAv0[fItrk]/F"); //Cosine of Pionting Angle
 374   fTreeV0->Branch("fDecayRadiusTree",fDecayRadiusTree,"fDecayRadiusTree[fItrk]/F"); //decay radius
 375
 376   fTreeV0->Branch("fChargeComboDeuteronPionTree",fChargeComboDeuteronPionTree,"fChargeComboDeuteronPionTree[fItrk]/I");
 377   fTreeV0->Branch("fInvaMassDeuteronPionTree",fInvaMassDeuteronPionTree,"fInvaMassDeuteronPionTree[fItrk]/D"); //invariant mass
 378
 379   fTreeV0->Branch("fIsCorrectlyAssociated",fIsCorrectlyAssociated,"fIsCorrectlyAssociated[fItrk]/O"); //associated hypertriton
 380
 381   fTreeV0->Branch("fAmenterosAlphaTree",fAmenterosAlphaTree,"fAmenterosAlphaTree[fItrk]/D");
 382   fTreeV0->Branch("fAmenterosQtTree",fAmenterosQtTree,"fAmenterosQtTree[fItrk]/D");
 383   fTreeV0->Branch("fRotationTree",fRotationTree,"fRotationTree[fItrk]/I");
 384
 385   fTreeV0->Branch("fImpactParameterDeuteronPos",fImpactParameterDeuteronPos,"fImpactParameterDeuteronPos[fItrk]/D");
 386   fTreeV0->Branch("fImpactParameterDeuteronNeg",fImpactParameterDeuteronNeg,"fImpactParameterDeuteronNeg[fItrk]/D");
 387   fTreeV0->Branch("fImpactParameterPionPos",fImpactParameterPionPos,"fImpactParameterPionPos[fItrk]/D");
 388   fTreeV0->Branch("fImpactParameterPionNeg",fImpactParameterPionNeg,"fImpactParameterPionNeg[fItrk]/D");
 389
 390   fTreeV0->Branch("fImpactParameterDeuteronPosAliKF",fImpactParameterDeuteronPosAliKF,"fImpactParameterDeuteronPosAliKF[fItrk]/D");
 391   fTreeV0->Branch("fImpactParameterDeuteronNegAliKF",fImpactParameterDeuteronNegAliKF,"fImpactParameterDeuteronNegAliKF[fItrk]/D");
 392   fTreeV0->Branch("fImpactParameterPionPosAliKF",fImpactParameterPionPosAliKF,"fImpactParameterPionPosAliKF[fItrk]/D");
 393   fTreeV0->Branch("fImpactParameterPionNegAliKF",fImpactParameterPionNegAliKF,"fImpactParameterPionNegAliKF[fItrk]/D");
 394
 395   fTreeV0->Branch("fMinNClustersTPCPos",fMinNClustersTPCPos,"fMinNClustersTPCPos[fItrk]/s");
 396   fTreeV0->Branch("fMinNClustersTPCNeg",fMinNClustersTPCNeg,"fMinNClustersTPCNeg[fItrk]/s");
 397
 398   fTreeV0->Branch("fMaxChi2PerClusterTPCPos",fMaxChi2PerClusterTPCPos,"fMaxChi2PerClusterTPCPos[fItrk]/D");
 399   fTreeV0->Branch("fMaxChi2PerClusterTPCNeg",fMaxChi2PerClusterTPCNeg,"fMaxChi2PerClusterTPCNeg[fItrk]/D");
 400
 401   //Armenteros-Podolanski
 402   fHistArmenterosPodolanskiDeuteronPion= new TH2F("fHistArmenterosPodolanskiDeuteronPion", "Armenteros-Podolanski d #pi^{-}", 200,-1.0,1.0, 500,0,1);
 403   fHistArmenterosPodolanskiDeuteronPion->GetXaxis()->SetTitle("#alpha");
 404   fHistArmenterosPodolanskiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
 405   fHistArmenterosPodolanskiDeuteronPion->SetMarkerStyle(kFullCircle);
 406
 407   fHistArmenterosPodolanskiAntiDeuteronPion= new TH2F("fHistArmenterosPodolanskiAntiDeuteronPion", "Armenteros-Podolanski #bar{d} #pi^{+}", 200,-1.0,1.0, 500,0,1);
 408   fHistArmenterosPodolanskiAntiDeuteronPion->GetXaxis()->SetTitle("#alpha");
 409   fHistArmenterosPodolanskiAntiDeuteronPion->GetYaxis()->SetTitle("q_{t}");
 410   fHistArmenterosPodolanskiAntiDeuteronPion->SetMarkerStyle(kFullCircle);
 411
 412   // control histograms
 413   fof = new TH2F("fof", "OnTheFlyStatus ",5,0.5,5.5,2,-0.5,1.5);
 414   fof->GetYaxis()->SetBinLabel(1,"offline");
 415   fof->GetYaxis()->SetBinLabel(2,"onTheFly");
 416   fof->GetXaxis()->SetBinLabel(1,"total");
 417   fof->GetXaxis()->SetBinLabel(2,"dcaCut");
 418   fof->GetXaxis()->SetBinLabel(3,"cosCut");
 419   fof->GetXaxis()->SetBinLabel(4,"nucleonPID");
 420   fof->GetXaxis()->SetBinLabel(5,"pionPID");
 421   //fof->GetXaxis()->SetBinLabel(6,"decayRadiusCut");
 422   //fof->SetMarkerStyle(kFullCircle);
 423
 424   //histogram to count number of events
 425   fHistCentralityClass10  = new TH1F("fHistCentralityClass10", "centrality with class10", 11, -0.5, 10.5);
 426   fHistCentralityClass10->GetXaxis()->SetTitle("Centrality");
 427   fHistCentralityClass10->GetYaxis()->SetTitle("Entries");
 428
 429   fHistCentralityPercentile  = new TH1F("fHistCentralityPercentile", "centrality with percentile", 101, -0.1, 100.1);
 430   fHistCentralityPercentile->GetXaxis()->SetTitle("Centrality");
 431   fHistCentralityPercentile->GetYaxis()->SetTitle("Entries");
 432
 433   //trigger statitics histogram
 434   fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
 435   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 436   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 437     fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 438
 439   fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", fNTriggers,-0.5,fNTriggers-0.5);
 440   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 441     fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 442
 443   //QA dE/dx
 444   fHistDeDxQA = new TH3F("fHistDeDxQA", "QA dE/dx", 400, -6.0, 6.0, 500, 0.0, 2000, 15, -0.5, 14.5);
 445   fHistDeDxQA->GetYaxis()->SetTitle("TPC Signal");
 446   fHistDeDxQA->GetXaxis()->SetTitle("p (GeV/c)");
 447   fHistDeDxQA->GetZaxis()->SetBinLabel(0,"all pos v0 tracks");
 448   fHistDeDxQA->GetZaxis()->SetBinLabel(1,"all neg v0 tracks");
 449   fHistDeDxQA->GetZaxis()->SetBinLabel(2,"all neg deuteron");
 450   fHistDeDxQA->GetZaxis()->SetBinLabel(3,"all pos deuteron");
 451   fHistDeDxQA->GetZaxis()->SetBinLabel(6,"all selected tracks");
 452   fHistDeDxQA->GetZaxis()->SetBinLabel(7,"neg deuteron for deuteron+pion");
 453   fHistDeDxQA->GetZaxis()->SetBinLabel(8,"pos pion for deuteron+pion");
 454   fHistDeDxQA->GetZaxis()->SetBinLabel(9,"pos deuteron for deuteron+pion");
 455   fHistDeDxQA->GetZaxis()->SetBinLabel(10,"neg pion for deuteron+pion");
 456
 457   //Define and fill the OutputContainer
 458   fOutputContainer = new TObjArray(1);
 459   fOutputContainer->SetOwner(kTRUE);
 460   fOutputContainer->SetName(GetName()) ;
 461   fOutputContainer->Add(fHistArmenterosPodolanskiDeuteronPion);
 462   fOutputContainer->Add(fHistArmenterosPodolanskiAntiDeuteronPion);
 463   fOutputContainer->Add(fof);
 464   fOutputContainer->Add(fHistDeDxQA);
 465   fOutputContainer->Add(fHistCentralityClass10);
 466   fOutputContainer->Add(fHistCentralityPercentile);
 467   fOutputContainer->Add(fHistTriggerStat);
 468   fOutputContainer->Add(fHistTriggerStatAfterEventSelection);
 469   fOutputContainer->Add(fHistLambdaNeutronPtGen);
 470   fOutputContainer->Add(fHistLambdaNeutronPtGenMinBias);
 471   fOutputContainer->Add(fHistLambdaNeutronPtGenCentral);
 472   fOutputContainer->Add(fHistLambdaNeutronPtGenSemiCentral);
 473   fOutputContainer->Add(fHistAntiLambdaNeutronPtGen);
 474   fOutputContainer->Add(fHistAntiLambdaNeutronPtGenMinBias);
 475   fOutputContainer->Add(fHistAntiLambdaNeutronPtGenCentral);
 476   fOutputContainer->Add(fHistAntiLambdaNeutronPtGenSemiCentral);
 477   fOutputContainer->Add(fHistLambdaNeutronInvaMassGen);
 478   fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassGen);
 479   fOutputContainer->Add(fHistLambdaNeutronDecayLengthGen);
 480   fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthGen);
 481   fOutputContainer->Add(fHistLambdaNeutronPtAso);
 482   fOutputContainer->Add(fHistLambdaNeutronPtAsoCuts);
 483   fOutputContainer->Add(fHistAntiLambdaNeutronPtAso);
 484   fOutputContainer->Add(fHistAntiLambdaNeutronPtAsoCuts);
 485   fOutputContainer->Add(fHistLambdaNeutronInvaMassAso);
 486   fOutputContainer->Add(fHistAntiLambdaNeutronInvaMassAso);
 487   fOutputContainer->Add(fHistLambdaNeutronDecayLengthAso);
 488   fOutputContainer->Add(fHistAntiLambdaNeutronDecayLengthAso);
 489 }
 490
 491 //________________________________________________________________________
 492 void AliAnalysisTaskLambdaNAOD::UserExec(Option_t *){
 493
 494   // Main loop
 495   // Called for each event
 496
 497   //cout << "katze1" << endl;
 498
 499   //Data
 500   //get Event-Handler for the trigger information
 501   fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 502   if (!fEventHandler) {
 503     AliError("Could not get InputHandler");
 504     //return -1;
 505     return;
 506   }
 507
 508   // Monte Carlo
 509   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 510   if (!eventHandler){
 511     //Printf("ERROR: Could not retrieve MC event handler");
 512     fMCtrue = kFALSE;
 513   }
 514
 515   AliMCEvent* mcEvent = 0x0;
 516   AliStack* stack = 0x0;
 517   if (eventHandler) mcEvent = eventHandler->MCEvent();
 518   if (!mcEvent){
 519     //Printf("ERROR: Could not retrieve MC event");
 520     if (fMCtrue) return;
 521   }
 522
 523   if (fMCtrue){
 524     stack = mcEvent->Stack();
 525     if (!stack) return;
 526   }
 527
 528   //look for the generated particles
 529   //MCGenerated(stack);
 530
 531   //cout << "katze2" << endl;
 532   if (SetupEvent() < 0) {
 533     PostData(1, fOutputContainer);
 534     return;
 535     }
 536
 537   //DATA
 538   AliESDEvent *fESDevent = 0x0;
 539   AliAODEvent *fAODevent = 0x0;
 540
 541
 542   if(!fMCtrue){
 543     if(!fPIDResponse) {
 544       AliError("Cannot get pid response");
 545       return;
 546     }
 547   }
 548
 549   Int_t centralityClass10 = -1;
 550   Int_t centralityPercentile = -1;
 551   Double_t vertex[3]          = {-100.0, -100.0, -100.0};
 552
 553   //Initialisation of the event and basic event cuts:
 554   //1.) ESDs:
 555   if (fAnalysisType == "ESD") {
 556
 557     fESDevent=dynamic_cast<AliESDEvent*>(InputEvent());
 558     if (!fESDevent) {
 559       AliWarning("ERROR: fESDevent not available \n");
 560       return;
 561     }
 562
 563     //Basic event cuts:
 564     //1.) vertex existence
 565     /*const AliESDVertex *vertexESD = fESDevent->GetPrimaryVertexTracks();
 566     if (vertexESD->GetNContributors()<1)
 567       {
 568         // SPD vertex
 569         vertexESD = fESDevent->GetPrimaryVertexSPD();
 570         if(vertexESD->GetNContributors()<1) {
 571           PostData(1, fOutputContainer);
 572           return;
 573         }
 574       }
 575
 576     vertexESD->GetXYZ(vertex);
 577
 578     //2. vertex position within 10 cm
 579     if (TMath::Abs(vertexESD->GetZv()) > 10) return;*/
 580
 581     const AliESDVertex *vertexTracks = fESDevent->GetPrimaryVertexTracks();
 582     if (vertexTracks->GetNContributors()<1) vertexTracks = 0x0;
 583
 584     const AliESDVertex *vertexSPD    = fESDevent->GetPrimaryVertexSPD();
 585     if (vertexSPD->GetNContributors()<1) vertexSPD = 0x0;
 586
 587     //cout << "before" << endl;
 588     if(!fMCtrue){
 589       if(!vertexTracks || !vertexSPD) return;
 590        //cout << "after" <<endl;
 591
 592       //if (vertexTracks && vertexSPD){
 593       //cout << "Vertex: " << TMath::Abs(vertexTracks->GetZv()) << endl;
 594       if (TMath::Abs(vertexTracks->GetZv()) > 10 || TMath::Abs(vertexSPD->GetZv()) > 10) return;
 595       //}
 596     }
 597
 598     //centrality selection in PbPb
 599     if (fESDevent->GetEventSpecie() == 4) //species == 4 == PbPb
 600       { // PbPb
 601         AliCentrality *esdCentrality = fESDevent->GetCentrality();
 602         centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 603         centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
 604         //cout << "************************************EventSpecie " << fESDevent->GetEventSpecie() << " centrality "<< centrality << endl;
 605         if(!fMCtrue){
 606           if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
 607         }
 608       }
 609
 610     //cout << "EventSpecie " << fESDevent->GetEventSpecie() << " centrality "<< centrality << endl;
 611     //count centrality classes
 612     fHistCentralityClass10->Fill(centralityClass10);
 613     fHistCentralityPercentile->Fill(centralityPercentile);
 614
 615     if(fTriggerFired[0] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(0);
 616     if(fTriggerFired[1] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(1);
 617     if(fTriggerFired[2] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(2);
 618     if(fTriggerFired[3] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(3);
 619     if(fTriggerFired[4] == kTRUE)fHistTriggerStatAfterEventSelection->Fill(4);
 620
 621  }
 622
 623 //2.) AODs:
 624  else if (fAnalysisType == "AOD") {
 625
 626    fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 627    if (!fAODevent) {
 628      AliWarning("ERROR: lAODevent not available \n");
 629      return;
 630    }
 631
 632    const AliAODVertex *vertexAOD = fAODevent->GetPrimaryVertex();
 633     if (vertexAOD->GetNContributors()<1)
 634       {
 635         // SPD vertex
 636         vertexAOD = fAODevent->GetPrimaryVertexSPD();
 637         if(vertexAOD->GetNContributors()<1) {
 638           PostData(1, fOutputContainer);
 639           return;
 640         }
 641       }
 642     vertexAOD->GetXYZ(vertex);
 643
 644     //2. vertex position within 10 cm
 645     if (TMath::Abs(vertex[2]) > 10) return;
 646
 647     //centrality selection
 648    AliCentrality *aodCentrality = fAODevent->GetCentrality();
 649    centralityClass10 = aodCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 650    centralityPercentile = aodCentrality->GetCentralityPercentile("V0M");
 651    if (centralityPercentile < 0. || centralityPercentile > 80.) return; //select only events with centralities between 0 and 80 %
 652
 653    // count number of events
 654    fHistCentralityClass10->Fill(centralityClass10);
 655
 656  } else {
 657
 658    Printf("Analysis type (ESD or AOD) not specified \n");
 659    return;
 660
 661  }
 662
 663
 664  //v0 loop
 665
 666
 667   fItrk = 0;
 668
 669  Int_t runNumber = 0;
 670  runNumber = (InputEvent())->GetRunNumber();
 671
 672  Int_t    nTrackMultiplicity             = -1;
 673  nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
 674  //const AliMulitiplicity *spdMultiplicity = 0x0;
 675  //if (fAnalysisType == "ESD")spdMultiplicity= fESDevent->GetMultiplicity();
 676  Int_t refMultTpc = -1;
 677  //const AliMulitiplicity *spdMultiplicity = 0x0;
 678  Int_t multi1 =-1;
 679  Int_t multi2 =-1;
 680  Int_t multi3 =-1;
 681
 682  // EstimateMultiplicity(&multi1,&multi2,&multi3, 0.9,kTRUE, kTRUE);
 683
 684  if (fAnalysisType == "ESD"){
 685    fESDevent->EstimateMultiplicity(multi1,multi2,multi3,0.9,kTRUE,kTRUE);
 686    refMultTpc = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent, kTRUE);
 687    //   const AliMultiplicity *spdMultiplicity = fESDevent->GetMultiplicity();
 688    //cout << "nTrackMultiplicity: " << nTrackMultiplicity << " multi1: " << multi1 << " multi2: " << multi2 << " multi3: " << multi3 << " refMultTpc: " << refMultTpc << endl;
 689  }
 690
 691  //look for the generated particles
 692  MCGenerated(stack,refMultTpc);
 693
 694  //AliKFParticle settings
 695  AliKFParticle::SetField(fInputEvent->GetMagneticField());
 696  AliKFVertex primVtx(*(fInputEvent->GetPrimaryVertex()));
 697
 698  for (Int_t ivertex=0; ivertex<(InputEvent()->GetNumberOfV0s()); ivertex++) { //beginn v0 loop
 699
 700    AliESDv0 * v0ESD = 0x0;
 701    AliAODv0 * v0AOD = 0x0;
 702
 703    Bool_t v0ChargesCorrect = kTRUE;
 704    Bool_t testTrackCuts = kFALSE;
 705    //Bool_t testFilterBit = kFALSE;
 706
 707    Int_t of = 7;
 708    Int_t onFlyStatus = 5;
 709
 710    Float_t dcaV0 = -1;
 711    Float_t cosPointing= 2;
 712    Float_t decayRadius= -1;
 713
 714    AliVTrack * trackN = 0x0;
 715    AliVTrack * trackP = 0x0;
 716
 717    Double_t ptotN = -1000;
 718    Double_t ptotP = -1000;
 719
 720    Int_t chargeComboDeuteronPion = -1;
 721
 722    Double_t momentumPion[3]={0,0,0};
 723    Double_t momentumPionRot[3]={0,0,0};
 724    Double_t momentumDeuteron[3]={0,0,0};
 725    Double_t momentumDeuteronRot[3]={0,0,0};
 726    Double_t momentumMother[3] = {0,0,0};
 727
 728    Double_t transversMomentumPion = 0;
 729    Double_t transversMomentumDeuteron = 0;
 730
 731    Double_t transversMomentumMother = 0;
 732    Double_t longitudinalMomentumMother = 0;
 733
 734    Double_t totalEnergyMother = 0;
 735    Double_t energyPion = 0;
 736    Double_t energyDeuteron = 0;
 737
 738    Double_t rapidity = 2;
 739
 740    TVector3 vecPion(0,0,0);
 741    TVector3 vecDeuteron(0,0,0);
 742    TVector3 vecMother(0,0,0);
 743
 744    Double_t alpha = 2;
 745    Double_t qt = -1;
 746
 747    Double_t thetaPion = 0;
 748    Double_t thetaDeuteron = 0;
 749
 750    Double_t phi=0;
 751    Double_t invaMassDeuteronPion = 0;
 752
 753    Double_t e12 = 0;
 754    Double_t r12   = 0;
 755    Double_t d22 = 0;
 756    Double_t dr22   = 0;
 757
 758    //Tree variables
 759    fV0finder[fItrk]         = -1;
 760    fkMB[fItrk]              = -1;
 761    fkCentral[fItrk]         = -1;
 762    fkSemiCentral[fItrk]     = -1;
 763    fkEMCEJE[fItrk]          = -1;
 764    fkEMCEGA[fItrk]          = -1;
 765
 766    fCentralityClass10[fItrk]  = -1;
 767    fCentralityPercentile[fItrk]  = -1;
 768    fMultiplicity[fItrk]     = -1;
 769    fRefMultiplicity[fItrk]     = -1;
 770
 771    fPtotN[fItrk]            = -1000;
 772    fPtotP[fItrk]            = -1000;
 773    fMotherPt[fItrk]         = -1000;
 774
 775    fdEdxN[fItrk]            = -1;
 776    fdEdxP[fItrk]            = -1;
 777    fSignN[fItrk]            = 0;
 778    fSignP[fItrk]            = 0;
 779
 780    fSigmadEdxPionPos[fItrk] = -1;
 781    fSigmadEdxPionNeg[fItrk] = -1;
 782
 783    fDCAv0[fItrk]            = -1;
 784    fCosinePAv0[fItrk]       = -2;
 785    fDecayRadiusTree[fItrk]  = -1;
 786
 787    fInvaMassDeuteronPionTree[fItrk] = 0;
 788    fChargeComboDeuteronPionTree[fItrk] = -1;
 789
 790    fAmenterosAlphaTree[fItrk] = 2;
 791    fAmenterosQtTree[fItrk] = -1;
 792
 793    fImpactParameterDeuteronPos[fItrk]            = -1;
 794    fImpactParameterDeuteronNeg[fItrk]            = -1;
 795
 796    fImpactParameterPionPos[fItrk]            = -1;
 797    fImpactParameterPionNeg[fItrk]            = -1;
 798
 799    fImpactParameterDeuteronPosAliKF[fItrk]            = -1;
 800    fImpactParameterDeuteronNegAliKF[fItrk]            = -1;
 801
 802    fImpactParameterPionPosAliKF[fItrk]            = -1;
 803    fImpactParameterPionNegAliKF[fItrk]            = -1;
 804
 805    fMinNClustersTPCPos[fItrk]            = 0;
 806    fMinNClustersTPCNeg[fItrk]            = 0;
 807
 808    fMaxChi2PerClusterTPCPos[fItrk]            = -1;
 809    fMaxChi2PerClusterTPCNeg[fItrk]            = -1;
 810
 811    //Get v0 object
 812    if(fAnalysisType == "ESD")v0ESD = fESDevent->GetV0(ivertex);
 813    if(fAnalysisType == "AOD")v0AOD = fAODevent->GetV0(ivertex);
 814
 815
 816    //distinguish between the two V0 finders: 0 offline V0 finder / 1 online V0 finder
 817
 818    if(fAnalysisType == "ESD") of = v0ESD->GetOnFlyStatus();
 819    if(fAnalysisType == "AOD") of = v0AOD->GetOnFlyStatus();
 820
 821    if(of)  onFlyStatus= 1;
 822    if(!of) onFlyStatus= 0;
 823
 824    if(onFlyStatus==0)fof->Fill(1,0);
 825    if(onFlyStatus==1)fof->Fill(1,1);
 826
 827    //Get dca, cos of pointing angle and decay radius
 828
 829
 830    if(fAnalysisType == "ESD")
 831      {
 832        dcaV0 = v0ESD->GetDcaV0Daughters();
 833        cosPointing = v0ESD->GetV0CosineOfPointingAngle();
 834        decayRadius = v0ESD->GetRr();
 835      }
 836
 837    if(fAnalysisType == "AOD")
 838      {
 839        dcaV0 = v0AOD->DcaV0Daughters();
 840        cosPointing = v0AOD->CosPointingAngle(vertex);
 841        decayRadius = v0AOD->DecayLengthV0(vertex);
 842      }
 843
 844
 845       // select coresponding tracks
 846       if(fAnalysisType == "ESD")
 847         {
 848           trackN = fESDevent->GetTrack(v0ESD->GetIndex(0));
 849           trackP = fESDevent->GetTrack(v0ESD->GetIndex(1));
 850
 851           if (trackN->Charge() > 0) // switch because of bug in V0 interface
 852             {
 853               trackN = fESDevent->GetTrack(v0ESD->GetIndex(1));
 854               trackP = fESDevent->GetTrack(v0ESD->GetIndex(0));
 855               v0ChargesCorrect = kFALSE;
 856             }
 857           //Track-Cuts
 858           testTrackCuts = TrackCuts(trackN,testTrackCuts);
 859           if(testTrackCuts == kFALSE) continue;
 860           testTrackCuts = TrackCuts(trackP,testTrackCuts);
 861           if(testTrackCuts == kFALSE) continue;
 862         }
 863
 864       if(fAnalysisType == "AOD")
 865         {
 866           trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
 867           trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
 868
 869           if (trackN->Charge() > 0) // switch because of bug in V0 interface
 870             {
 871               trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
 872               trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
 873               v0ChargesCorrect = kFALSE;
 874             }
 875           //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
 876           //testFilterBit = FilterBit(trackN,testFilterBit);
 877           //if(testFilterBit == kFALSE) continue;
 878           //testFilterBit = FilterBit(trackP,testFilterBit);
 879           //if(testFilterBit == kFALSE) continue;
 880           //Track-Cuts
 881           /*testTrackCuts = TrackCuts(trackN,testTrackCuts);
 882           if(testTrackCuts == kFALSE) continue;
 883           testTrackCuts = TrackCuts(trackP,testTrackCuts);
 884           if(testTrackCuts == kFALSE) continue;*/
 885         }
 886
 887       //Get the total momentum for each track ---> at the inner readout of the TPC???? // momenta a always positive
 888             if(fAnalysisType == "AOD") {
 889         ptotN = trackN->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
 890         ptotP = trackP->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
 891       }
 892
 893       if(fAnalysisType == "ESD") {
 894         ptotN =  MomentumInnerParam(trackN,ptotN);
 895         ptotP =  MomentumInnerParam(trackP,ptotP);
 896       }
 897
 898       //fill QA dEdx with all V0 candidates
 899       if(trackP) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),0);
 900       if(trackN) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),1);
 901
 902       if (dcaV0 > 3 || dcaV0 < 1e-20) continue;
 903
 904       //Check how much the dca cut reduces the statistic (background) for the different VO-finder
 905       if(onFlyStatus==0)fof->Fill(2,0);
 906       if(onFlyStatus==1)fof->Fill(2,1);
 907
 908       if (cosPointing < 0.9) continue;
 909
 910       //Check how much the cos-of-the-pointing-angle-cut reduces the statistic (background) for the different VO-finder
 911       if(onFlyStatus==0)fof->Fill(3,0);
 912       if(onFlyStatus==1)fof->Fill(3,1);
 913
 914       //deuteron PID via specific energy loss in the TPC
 915       Bool_t isDeuteron[3] = {kFALSE,kFALSE,kFALSE}; //0 = posDeuteron, 1 = negDeuteron, 2 = trackN is deuteron
 916
 917       DeuteronPID(trackP,trackN,ptotP,ptotN,runNumber,isDeuteron);
 918
 919       //if(!isDeuteron[0] && !isDeuteron[1]) continue;
 920       if((isDeuteron[0]==kFALSE) && (isDeuteron[1]==kFALSE)) continue;
 921
 922
 923       //Check how much the nuclei PID cuts reduce the statistics (background) for the two V0-finders
 924       if(onFlyStatus==0)fof->Fill(4,0);
 925       if(onFlyStatus==1)fof->Fill(4,1);
 926
 927       //Fill the QA dEdx with deuterons and helium3 after the nuclei PID cut
 928       if(isDeuteron[1]) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),2);
 929       if(isDeuteron[0]) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),3);
 930
 931       //deuteron PID via specific energy loss in the TPC
 932       Bool_t isPion[2] = {kFALSE,kFALSE}; //0 = posPion, 1 = negPion
 933
 934       PionPID(trackP,trackN,ptotP,ptotN,runNumber,isPion);
 935
 936       //if(isDeuteron[0] && !isPion[1]) continue; //pos deuteron and neg Pion
 937       //if(isDeuteron[1] && !isPion[0]) continue; //neg deuteron and pos Pion
 938
 939       //Check how much the pion PID cuts reduce the statistics (background) for the two V0-finders
 940       if(onFlyStatus==0)fof->Fill(5,0);
 941       if(onFlyStatus==1)fof->Fill(5,1);
 942
 943
 944       //Save the different charge combinations to differentiat between particles and anit-particles:
 945       // -/+ = Anti-deuteron + positive Pion
 946       // -/- = Anti-deuteron + negative Pion
 947       // +/- = deuteron + negative Pion
 948       // +/+ = deuteron + positive Pion
 949
 950       //Like-sign
 951       if (trackN->Charge()<0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 1; // -/-
 952       if (trackN->Charge()>0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 3; // +/+
 953
 954       //unlinke-sign
 955       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+
 956       //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+ //should not exist because of charge correction
 957       //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/- //should not exist because of charge correction
 958       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/-
 959
 960       //if(chargeComboDeuteronPion==0) cout << "chargeN: " << trackN->Charge() << " chargeP: " << trackP->Charge() << endl;
 961
 962
 963       //Fill the QA dEdx with all selected tracks after the PID cuts
 964       fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),6);
 965       fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),6);
 966
 967
 968
 969       //get the momenta of the daughters
 970
 971       if(fAnalysisType == "ESD")
 972         {
 973           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
 974
 975             v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 976             v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 977             if (!v0ChargesCorrect) {
 978
 979               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 980               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 981             }
 982           }
 983
 984           if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
 985
 986             v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 987             v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 988             if (!v0ChargesCorrect){
 989
 990               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 991               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
 992             }
 993           }
 994
 995           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
 996             if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
 997
 998               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
 999               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1000               if (!v0ChargesCorrect) {
1001
1002                 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1003                 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1004               }
1005             }
1006             if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1007
1008               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1009               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1010               if (!v0ChargesCorrect) {
1011
1012                 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
1013                 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
1014               }
1015             }
1016           }
1017
1018
1019           //get the momenta of the mother
1020           v0ESD->GetPxPyPz(momentumMother[0],momentumMother[1],momentumMother[2]);
1021         }
1022
1023
1024       if(fAnalysisType == "AOD")
1025         {
1026           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
1027
1028             momentumPion[0] = v0AOD->MomPosX();
1029             momentumPion[1] = v0AOD->MomPosY();
1030             momentumPion[2] = v0AOD->MomPosZ();
1031
1032             momentumDeuteron[0] = v0AOD->MomNegX();
1033             momentumDeuteron[1] = v0AOD->MomNegY();
1034             momentumDeuteron[2] = v0AOD->MomNegZ();
1035
1036             if (!v0ChargesCorrect){
1037
1038               momentumPion[0] = v0AOD->MomNegX();
1039               momentumPion[1] = v0AOD->MomNegY();
1040               momentumPion[2] = v0AOD->MomNegZ();
1041
1042               momentumDeuteron[0] = v0AOD->MomPosX();
1043               momentumDeuteron[1] = v0AOD->MomPosY();
1044               momentumDeuteron[2] = v0AOD->MomPosZ();
1045             }
1046           }
1047
1048           if (chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
1049
1050             momentumPion[0] = v0AOD->MomNegX();
1051             momentumPion[1] = v0AOD->MomNegY();
1052             momentumPion[2] = v0AOD->MomNegZ();
1053
1054             momentumDeuteron[0] = v0AOD->MomPosX();
1055             momentumDeuteron[1] = v0AOD->MomPosY();
1056             momentumDeuteron[2] = v0AOD->MomPosZ();
1057
1058             if (!v0ChargesCorrect){
1059
1060               momentumPion[0] = v0AOD->MomPosX();
1061               momentumPion[1] = v0AOD->MomPosY();
1062               momentumPion[2] = v0AOD->MomPosZ();
1063
1064               momentumDeuteron[0] = v0AOD->MomNegX();
1065               momentumDeuteron[1] = v0AOD->MomNegY();
1066               momentumDeuteron[2] = v0AOD->MomNegZ();
1067             }
1068           }
1069
1070           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN correponds to deuteron or pion, trackP corresponds to deuteron or pion
1071             if(isDeuteron[2]==kTRUE){ //trackN correponds to deuteron, trackP corresponds to pion
1072
1073               momentumPion[0] = v0AOD->MomPosX();
1074               momentumPion[1] = v0AOD->MomPosY();
1075               momentumPion[2] = v0AOD->MomPosZ();
1076
1077               momentumDeuteron[0] = v0AOD->MomNegX();
1078               momentumDeuteron[1] = v0AOD->MomNegY();
1079               momentumDeuteron[2] = v0AOD->MomNegZ();
1080
1081               if (!v0ChargesCorrect){
1082
1083                 momentumPion[0] = v0AOD->MomNegX();
1084                 momentumPion[1] = v0AOD->MomNegY();
1085                 momentumPion[2] = v0AOD->MomNegZ();
1086
1087                 momentumDeuteron[0] = v0AOD->MomPosX();
1088                 momentumDeuteron[1] = v0AOD->MomPosY();
1089                 momentumDeuteron[2] = v0AOD->MomPosZ();
1090               }
1091             }
1092
1093             if(isDeuteron[2]==kFALSE){ //trackP corresponds to deuteron, trackN corresponds to pion
1094
1095               momentumPion[0] = v0AOD->MomNegX();
1096               momentumPion[1] = v0AOD->MomNegY();
1097               momentumPion[2] = v0AOD->MomNegZ();
1098
1099               momentumDeuteron[0] = v0AOD->MomPosX();
1100               momentumDeuteron[1] = v0AOD->MomPosY();
1101               momentumDeuteron[2] = v0AOD->MomPosZ();
1102
1103               if (!v0ChargesCorrect){
1104
1105                 momentumPion[0] = v0AOD->MomPosX();
1106                 momentumPion[1] = v0AOD->MomPosY();
1107                 momentumPion[2] = v0AOD->MomPosZ();
1108
1109                 momentumDeuteron[0] = v0AOD->MomNegX();
1110                 momentumDeuteron[1] = v0AOD->MomNegY();
1111                 momentumDeuteron[2] = v0AOD->MomNegZ();
1112               }
1113             }
1114
1115             //get the momenta of the mother
1116             momentumMother[0] = v0AOD->MomV0X();
1117             momentumMother[1] = v0AOD->MomV0Y();
1118             momentumMother[2] = v0AOD->MomV0Z();
1119           }
1120         }
1121
1122       //Rapidity - cut
1123       transversMomentumPion      = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]);
1124       transversMomentumDeuteron  = TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]);
1125
1126       transversMomentumMother=TMath::Sqrt((momentumDeuteron[0]+momentumPion[0])*(momentumDeuteron[0]+momentumPion[0])+(momentumDeuteron[1]+momentumPion[1])*(momentumDeuteron[1]+momentumPion[1]));
1127
1128       longitudinalMomentumMother = TMath::Sqrt((momentumDeuteron[2]+momentumPion[2])*(momentumDeuteron[2]+momentumPion[2]));
1129
1130       energyDeuteron =  TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]+momentumDeuteron[2]*momentumDeuteron[2]+fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]);
1131
1132       energyPion = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]+momentumPion[2]*momentumPion[2]+fgkMass[kMassPion]*fgkMass[kMassPion]);
1133
1134       totalEnergyMother = energyPion + energyDeuteron;
1135
1136       //rapidity cut +- 1
1137       rapidity = 0.5 * TMath::Log((totalEnergyMother+longitudinalMomentumMother)/( totalEnergyMother-longitudinalMomentumMother));
1138
1139       if(rapidity > 1.0 || rapidity < -1) continue;
1140
1141
1142       //Armanteros-Podolanski
1143       vecPion.SetXYZ(momentumPion[0],momentumPion[1],momentumPion[2]);
1144       vecDeuteron.SetXYZ(momentumDeuteron[0],momentumDeuteron[1],momentumDeuteron[2]);
1145       vecMother.SetXYZ(momentumMother[0],momentumMother[1],momentumMother[2]);
1146
1147       thetaPion = acos((vecPion * vecMother)/(vecPion.Mag() * vecMother.Mag()));
1148       thetaDeuteron = acos((vecDeuteron * vecMother)/(vecDeuteron.Mag() * vecMother.Mag()));
1149
1150       alpha = ((vecPion.Mag())*cos(thetaPion)-(vecDeuteron.Mag())*cos(thetaDeuteron))/
1151         ((vecDeuteron.Mag())*cos(thetaDeuteron)+(vecPion.Mag())*cos(thetaPion)) ;
1152       qt = vecDeuteron.Mag()*sin(thetaDeuteron);
1153
1154
1155       //Rotation for background calculation
1156       //Int_t rotation=1; // =1 signal, =2 Rotation of the pion , =3 Rotation of the deuteron
1157
1158       //Double_t fStartAnglePhi=TMath::Pi();
1159       //Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1160       //phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1161
1162       for(Int_t rotation=1;rotation<4;rotation++){ //loop for rotation
1163
1164         Double_t fStartAnglePhi=TMath::Pi();
1165         Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1166         phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1167
1168         //calculate new rotated momenta
1169         momentumPionRot[0]=TMath::Cos(phi)*momentumPion[0]-TMath::Sin(phi)*momentumPion[1];
1170         momentumPionRot[1]=TMath::Sin(phi)*momentumPion[0]+TMath::Cos(phi)*momentumPion[1];
1171
1172         momentumDeuteronRot[0]=TMath::Cos(phi)*momentumDeuteron[0]-TMath::Sin(phi)*momentumDeuteron[1];
1173         momentumDeuteronRot[1]=TMath::Sin(phi)*momentumDeuteron[0]+TMath::Cos(phi)*momentumDeuteron[1];
1174
1175         //invariant mass calculations
1176         fIsCorrectlyAssociated[fItrk] = kFALSE;
1177
1178         //factor for the invariant mass calculation, which only include the pion
1179         e12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1180           +momentumPion[0]*momentumPion[0]
1181           +momentumPion[1]*momentumPion[1]
1182           +momentumPion[2]*momentumPion[2];
1183
1184         r12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1185           +momentumPionRot[0]*momentumPionRot[0]
1186           +momentumPionRot[1]*momentumPionRot[1]
1187           +momentumPion[2]*momentumPion[2];
1188
1189         //factor for the invariant mass calculation, which only include the deuterons
1190         d22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1191           +momentumDeuteron[0]*momentumDeuteron[0]
1192           +momentumDeuteron[1]*momentumDeuteron[1]
1193           +momentumDeuteron[2]*momentumDeuteron[2];
1194         dr22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1195           +momentumDeuteronRot[0]*momentumDeuteronRot[0]
1196           +momentumDeuteronRot[1]*momentumDeuteronRot[1]
1197           +momentumDeuteron[2]*momentumDeuteron[2];
1198
1199         if(rotation == 1){ //signal
1200           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1201                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1202                                                         + 2.*(TMath::Sqrt(e12*d22)
1203                                                               -momentumPion[0]*momentumDeuteron[0]
1204                                                               -momentumPion[1]*momentumDeuteron[1]
1205                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1206
1207           if (fMCtrue){
1208
1209             Int_t labelN = 0;
1210             Int_t labelP = 0;
1211             labelN = trackN->GetLabel();
1212             labelP = trackP->GetLabel();
1213
1214             TParticle *tparticleDaughterN = stack->Particle(TMath::Abs(labelN));
1215             TParticle *tparticleDaughterP = stack->Particle(TMath::Abs(labelP));
1216
1217             Int_t labelMotherN = tparticleDaughterN->GetFirstMother();
1218             Int_t labelMotherP = tparticleDaughterP->GetFirstMother();
1219
1220             TParticle *tparticleMotherN = stack->Particle(TMath::Abs(labelMotherN));
1221             TParticle *tparticleMotherP = stack->Particle(TMath::Abs(labelMotherP));
1222
1223             //LambdaNeuteron
1224             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
1225
1226               //cout << "number of daughters: " << tparticleMotherN->GetNDaughters() << endl;;
1227               if(tparticleMotherN->GetNDaughters() > 2.) continue;
1228
1229               Int_t labelSecondDaughter = tparticleMotherP->GetDaughter(1);
1230               Int_t labelFirstDaughter = labelSecondDaughter-1;
1231
1232               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1233               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1234
1235               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGDeuteron]){//check first daughter PDG
1236
1237                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionMinus]){//check second daughter PDG
1238
1239                   if(invaMassDeuteronPion < 2.02) continue;
1240
1241                   fHistLambdaNeutronPtAso->Fill(transversMomentumMother);
1242                   fHistLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1243                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1244                   fHistLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1245                   //CUTS
1246
1247                   if((dcaV0 < 0.5) &&
1248                      (cosPointing > 0.999) &&
1249                      (decayRadius < 50.0) &&
1250                      (ptotP > 0.2)){//trackP->GetTPCsignal()>180 &&&&  decayRadius > 1.5 && decayRadius< 50
1251
1252                     fHistLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1253                   }
1254                 }//end check second daughter PDG
1255               }//end check first daughter PDG
1256             }//end LambdaNeutron
1257
1258             //Anti-LambdaNeutron
1259             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
1260
1261               Int_t labelSecondDaughter = tparticleMotherN->GetDaughter(1);
1262               Int_t labelFirstDaughter = labelSecondDaughter-1;
1263
1264               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1265               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1266
1267               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGAntiDeuteron]) {//check first daughter PDG
1268
1269                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionPlus]) {//check second daughter PDG
1270
1271                   if(invaMassDeuteronPion < 2.02) continue;
1272
1273                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1274                   fHistAntiLambdaNeutronPtAso->Fill(transversMomentumMother);
1275                   fHistAntiLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1276                   fHistAntiLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1277
1278                   //CUTS
1279                   if(dcaV0 < 1. && cosPointing > 0.99 && ptotN > 0.2)//&& trackN->GetTPCsignal()>180 && decayRadius > 1.5 && decayRadius< 50
1280                     {
1281                       fHistAntiLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1282                     }
1283                 }//end check second daughter PDG
1284               }//end check first daughter PDG
1285             }//end Anti-LambdaNeutron
1286           }//end MC
1287         }//end rotation == 1, signal
1288
1289         if(rotation == 2){ // rotation of the pion
1290           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1291                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1292                                                         + 2.*(TMath::Sqrt(r12*d22)
1293                                                               -momentumPionRot[0]*momentumDeuteron[0]
1294                                                               -momentumPionRot[1]*momentumDeuteron[1]
1295                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1296         }//end rotation == 2, rotation of the pion
1297
1298         if(rotation == 3){// Rotation of the deuteron
1299           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1300                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1301                                                         + 2.*(TMath::Sqrt(e12*dr22)
1302                                                               -momentumPion[0]*momentumDeuteronRot[0]
1303                                                               -momentumPion[1]*momentumDeuteronRot[1]
1304                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1305         }//end rotation == 3, rotation of the deuteron
1306
1307         //fill the THnSparse and the tree variables
1308
1309         //tree variables which are independent of the particle-species
1310         /*fV0finder[fItrk]         = onFlyStatus;
1311         fkMB[fItrk]              = fTriggerFired[0];
1312         fkCentral[fItrk]         = fTriggerFired[1];
1313         fkSemiCentral[fItrk]     = fTriggerFired[2];
1314         fkEMCEJE[fItrk]          = fTriggerFired[3];
1315         fkEMCEGA[fItrk]          = fTriggerFired[4];
1316
1317         fPtotN[fItrk]            = ptotN;
1318         fPtotP[fItrk]            = ptotP;
1319         fMotherPt[fItrk]         = transversMomentumMother;
1320
1321         fdEdxN[fItrk]            = trackN->GetTPCsignal();
1322         fdEdxP[fItrk]            = trackP->GetTPCsignal();
1323         fSignN[fItrk]            = trackN->Charge();
1324         fSignP[fItrk]            = trackP->Charge();
1325
1326         fDCAv0[fItrk]            = dcaV0;
1327         fCosinePAv0[fItrk]       = cosPointing;
1328         fDecayRadiusTree[fItrk]  = decayRadius;
1329
1330         fAmenterosAlphaTree[fItrk] = alpha;
1331         fAmenterosQtTree[fItrk] = qt;
1332         fRotationTree[fItrk] = rotation;*/
1333
1334
1335         if (isDeuteron[0] == kTRUE)  //pos deuteron
1336           {
1337             /*  fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1338             fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1339
1340             fItrk++;*/
1341
1342             if(invaMassDeuteronPion < 2.1)
1343               {
1344                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotP*(trackP->Charge()), trackP->GetTPCsignal(),9);
1345                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotN*(trackN->Charge()), trackN->GetTPCsignal(),10);
1346               }
1347             fHistArmenterosPodolanskiDeuteronPion->Fill(alpha*(-1),qt);
1348           }
1349
1350         if (isDeuteron[1] == kTRUE)
1351           {
1352             /*fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1353             fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1354
1355             fItrk++;*/
1356
1357             if(invaMassDeuteronPion < 2.1)
1358               {
1359                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),7);
1360                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),8);
1361               }
1362             fHistArmenterosPodolanskiAntiDeuteronPion->Fill(alpha,qt);
1363           }
1364
1365         UShort_t numberOfTPCclustersPos = 0;
1366         UShort_t numberOfTPCclustersNeg = 0;
1367
1368         numberOfTPCclustersPos = TPCclusters(trackP,numberOfTPCclustersPos);
1369         numberOfTPCclustersNeg = TPCclusters(trackN,numberOfTPCclustersNeg);
1370
1371         Double_t numberOfChi2clustersTPCPos = 10.;
1372         Double_t numberOfChi2clustersTPCNeg = 10.;
1373
1374         numberOfChi2clustersTPCPos = TPCchi2(trackP,numberOfChi2clustersTPCPos,numberOfTPCclustersPos);
1375         numberOfChi2clustersTPCNeg = TPCchi2(trackN,numberOfChi2clustersTPCNeg,numberOfTPCclustersNeg);
1376
1377         //ImpactParameters
1378         Double_t dcaPosToVertex = 0;
1379         Double_t dcaNegToVertex =  0;
1380
1381         dcaPosToVertex = ImpactParameter(trackP,dcaPosToVertex);
1382         dcaNegToVertex = ImpactParameter(trackN,dcaNegToVertex);
1383
1384         //ImpactParameter with AliKF
1385         //Get the primary vertex from the ITS for the following calculations, because it is the most precise determination of the primary vertex
1386         Double_t itsVertex[3]={0.,0.,0.};
1387
1388         if(fAnalysisType == "ESD"){
1389           itsVertex[0]=fESDevent->GetPrimaryVertexSPD()->GetX();
1390           itsVertex[1]=fESDevent->GetPrimaryVertexSPD()->GetY();
1391           itsVertex[2]=fESDevent->GetPrimaryVertexSPD()->GetZ();
1392         }
1393
1394         Double_t impactParameterPionPos = -1;
1395         Double_t impactParameterPionNeg = -1;
1396
1397         Double_t impactParameterDeuteronPos = -1;
1398         Double_t impactParameterDeuteronNeg = -1;
1399
1400         Double_t impactParameterPionPosAliKF = -1;
1401         Double_t impactParameterPionNegAliKF = -1;
1402
1403         Double_t impactParameterDeuteronPosAliKF = -1;
1404         Double_t impactParameterDeuteronNegAliKF = -1;
1405
1406         AliKFParticle* negPionKF;
1407         AliKFParticle* posPionKF;
1408         AliKFParticle* negDeuteronKF;
1409         AliKFParticle* posDeuteronKF;
1410
1411         if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN corresponds to anti-deuteron, tackP corresponds to pion
1412
1413           posPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1414           impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1415           impactParameterPionPos = dcaPosToVertex;
1416           negDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1417           impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1418           impactParameterDeuteronNeg = dcaNegToVertex;
1419         }
1420
1421         if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP corresponds to deuteron, tackN corresponds to pion
1422
1423           negPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1424           impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1425           impactParameterPionNeg = dcaNegToVertex;
1426           posDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1427           impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1428           impactParameterDeuteronPos = dcaPosToVertex;
1429         }
1430
1431         if(chargeComboDeuteronPion==1){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
1432           if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
1433
1434             negPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1435             impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1436             impactParameterPionNeg = dcaPosToVertex;
1437             negDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1438             impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1439             impactParameterDeuteronNeg = dcaNegToVertex;
1440           }
1441           if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1442
1443             negPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionPlus]);
1444             impactParameterPionNegAliKF = negPionKF->GetDistanceFromVertex(itsVertex);
1445             impactParameterPionNeg = dcaNegToVertex;
1446             negDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1447             impactParameterDeuteronNegAliKF = negDeuteronKF->GetDistanceFromVertex(itsVertex);
1448             impactParameterDeuteronNeg = dcaPosToVertex;
1449           }
1450         }
1451
1452         if(chargeComboDeuteronPion==3){ //trackN corresponds to deuteron or pion, trackP corresponds to deuteron or pion
1453           if(isDeuteron[2]==kTRUE){ //trackN corresponds to deuteron, trackP corresponds to pion
1454
1455             posPionKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1456             impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1457             impactParameterPionPos = dcaPosToVertex;
1458             posDeuteronKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionMinus]);
1459             impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1460             impactParameterDeuteronPos  = dcaNegToVertex;
1461           }
1462           if(isDeuteron[2]==kFALSE){ //trackP corresponds tp deuteron, trackN corresponds to pion
1463
1464             posPionKF = new AliKFParticle (*(trackN),fgkPdgCode[kPDGPionPlus]);
1465             impactParameterPionPosAliKF = posPionKF->GetDistanceFromVertex(itsVertex);
1466             impactParameterPionPos  = dcaNegToVertex;
1467             posDeuteronKF = new AliKFParticle (*(trackP),fgkPdgCode[kPDGPionPlus]);
1468             impactParameterDeuteronPosAliKF = posDeuteronKF->GetDistanceFromVertex(itsVertex);
1469             impactParameterDeuteronPos = dcaPosToVertex;
1470           }
1471         }
1472
1473         if((isDeuteron[0] == kTRUE) || (isDeuteron[1] == kTRUE)){
1474           //tree variables which are independent of the particle-species
1475           fV0finder[fItrk]         = onFlyStatus;
1476           fkMB[fItrk]              = fTriggerFired[0];
1477           fkCentral[fItrk]         = fTriggerFired[1];
1478           fkSemiCentral[fItrk]     = fTriggerFired[2];
1479           fkEMCEJE[fItrk]          = fTriggerFired[3];
1480           fkEMCEGA[fItrk]          = fTriggerFired[4];
1481
1482           fCentralityClass10[fItrk]= centralityClass10;
1483           fCentralityPercentile[fItrk]= centralityPercentile;
1484           fMultiplicity[fItrk]     = multi1;
1485           fRefMultiplicity[fItrk]  = refMultTpc;
1486
1487           fPtotN[fItrk]            = ptotN;
1488           fPtotP[fItrk]            = ptotP;
1489           fMotherPt[fItrk]         = transversMomentumMother;
1490
1491           fdEdxN[fItrk]            = trackN->GetTPCsignal();
1492           fdEdxP[fItrk]            = trackP->GetTPCsignal();
1493           fSignN[fItrk]            = trackN->Charge();
1494           fSignP[fItrk]            = trackP->Charge();
1495
1496           fSigmadEdxPionPos[fItrk] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion));
1497           fSigmadEdxPionNeg[fItrk] = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion));
1498
1499           fDCAv0[fItrk]            = dcaV0;
1500           fCosinePAv0[fItrk]       = cosPointing;
1501           fDecayRadiusTree[fItrk]  = decayRadius;
1502
1503           fAmenterosAlphaTree[fItrk] = alpha;
1504           fAmenterosQtTree[fItrk] = qt;
1505           fRotationTree[fItrk] = rotation;
1506
1507           fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1508           fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1509
1510           fImpactParameterDeuteronPos[fItrk] = impactParameterDeuteronPos;
1511           fImpactParameterDeuteronNeg[fItrk] = impactParameterDeuteronNeg;
1512
1513           fImpactParameterPionPos[fItrk] = impactParameterPionPos;
1514           fImpactParameterPionNeg[fItrk] = impactParameterPionNeg;
1515
1516           fImpactParameterDeuteronPosAliKF[fItrk] = impactParameterDeuteronPosAliKF;
1517           fImpactParameterDeuteronNegAliKF[fItrk] = impactParameterDeuteronNegAliKF;
1518
1519           fImpactParameterPionPosAliKF[fItrk] = impactParameterPionPosAliKF;
1520           fImpactParameterPionNegAliKF[fItrk] = impactParameterPionNegAliKF;
1521
1522           fMinNClustersTPCPos[fItrk]            = numberOfTPCclustersPos;
1523           fMinNClustersTPCNeg[fItrk]            = numberOfTPCclustersNeg;
1524
1525           //  cout << "numberOfTPCclustersPos: " << numberOfTPCclustersPos << " fMinNClustersTPCPos[fItrk]: "  << fMinNClustersTPCPos[fItrk] << endl;
1526
1527           fMaxChi2PerClusterTPCPos[fItrk]            = numberOfChi2clustersTPCPos;
1528           fMaxChi2PerClusterTPCNeg[fItrk]            = numberOfChi2clustersTPCNeg;
1529
1530           fItrk++;
1531         }
1532         //cout << "fItrk: " << fItrk-1 << " fChargeComboDeuteronPionTree[fItrk]: " << fChargeComboDeuteronPionTree[fItrk-1] << " fV0finder[fItrk]: " << fV0finder[fItrk-1] << " fInvaMassDeuteronPionTree[fItrk]: " << fInvaMassDeuteronPionTree[fItrk-1] << " fMaxChi2PerClusterTPCPos[fItrk]: " << fMaxChi2PerClusterTPCPos[fItrk-1] << " numberOfTPCclustersPos: " << numberOfTPCclustersPos << " fMinNClustersTPCPos[fItrk]: "  << fMinNClustersTPCPos[fItrk-1] << " fCentralityPercentile[fItrk]: " << fCentralityPercentile[fItrk-1]  << " fMotherPt[fItrk]: " << fMotherPt[fItrk-1] << endl;
1533
1534       }//end rotation loop
1535
1536  } //end of v0 loop
1537
1538   //fill the tree
1539   fTreeV0->Fill();
1540
1541   // Post output data.
1542   PostData(1, fOutputContainer);
1543   PostData(2, fTreeV0);
1544 }
1545
1546 //________________________________________________________________________
1547 void AliAnalysisTaskLambdaNAOD::Terminate(const Option_t *)
1548 {
1549   // Draw result to the screen
1550   // Called once at the end of the query
1551
1552   //get output data and darw 'fHistPt'
1553   if (!GetOutputData(0)) return;
1554   //TH1F *hist=(TH1F*)(((TObjArray*)GetOutputData(0))->FindObject("fHistPt"));
1555   //if (hist) hist->Draw();
1556 }
1557
1558 //_____________________________________________________________________________
1559 Int_t AliAnalysisTaskLambdaNAOD::Initialize() {
1560
1561
1562   // -- Reset Event
1563   // ----------------
1564   ResetEvent();
1565
1566   return 0;
1567   }
1568 //________________________________________________________________________
1569 Int_t AliAnalysisTaskLambdaNAOD::SetupEvent() {
1570   // Setup Reading of event
1571
1572   ResetEvent();
1573
1574   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1575   // -- Get Trigger
1576   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1577
1578   //Bool_t isTriggered = IsTriggered();
1579   IsTriggered();
1580   return 0;
1581 }
1582 //_____________________________________________________________________________
1583 void AliAnalysisTaskLambdaNAOD::ResetEvent() {
1584
1585   // Reset event
1586
1587   // -- Reset QA
1588   for (Int_t ii = 0; ii < fNTriggers; ++ii)
1589     fTriggerFired[ii] = kFALSE;
1590
1591   return;
1592 }
1593 //_______________________________________________________________________
1594 void AliAnalysisTaskLambdaNAOD::BinLogAxis(const THnSparse *h, Int_t axisNumber) {
1595   //
1596   // Method for the correct logarithmic binning of histograms
1597   //
1598   TAxis *axis = h->GetAxis(axisNumber);
1599   int bins = axis->GetNbins();
1600
1601   Double_t from = axis->GetXmin();
1602   Double_t to = axis->GetXmax();
1603   Double_t *newBins = new Double_t[bins + 1];
1604
1605   newBins[0] = from;
1606   Double_t factor = pow(to/from, 1./bins);
1607
1608   for (int i = 1; i <= bins; i++) {
1609    newBins[i] = factor * newBins[i-1];
1610   }
1611   axis->Set(bins, newBins);
1612   delete [] newBins;
1613
1614   }
1615 //________________________________________________________________________
1616 Bool_t AliAnalysisTaskLambdaNAOD::IsTriggered() {
1617   // Check if Event is triggered and fill Trigger Histogram
1618
1619   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
1620   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
1621   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
1622   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
1623   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
1624
1625   Bool_t isTriggered = kFALSE;
1626
1627   for (Int_t ii=0; ii<fNTriggers; ++ii) {
1628     if(fTriggerFired[ii]) {
1629       isTriggered = kTRUE;
1630       fHistTriggerStat->Fill(ii);
1631     }
1632   }
1633
1634   return isTriggered;
1635   }
1636 //______________________________________________________________________________
1637 Bool_t AliAnalysisTaskLambdaNAOD::DeuteronPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isDeuteron[3]) {
1638
1639      //define the arrays for the Bethe-Bloch-Parameters
1640       Double_t parDeuteron[5] = {0,0,0,0,0};
1641
1642       if(runNumber < 166500) //LHC10h
1643         {
1644           parDeuteron[0] = 1.45802; // ALEPH parameters for deuterons (pass2)
1645           parDeuteron[1] = 27.4992;
1646           parDeuteron[2] = 4.00313e-15;
1647           parDeuteron[3] = 2.48485;
1648           parDeuteron[4] = 8.31768;
1649         }
1650
1651       if(runNumber > 166500) //LHC11h
1652         {
1653           parDeuteron[0] = 1.11243; // ALEPH parameters for deuterons (pass2)
1654           parDeuteron[1] = 26.1144;
1655           parDeuteron[2] = 4.00313e-15;
1656           parDeuteron[3] = 2.72969 ;
1657           parDeuteron[4] = 9.15038;
1658         }
1659
1660       if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1661         parDeuteron[0] = 4.449; // ALEPH parameters for deuterons (pass2)
1662         parDeuteron[1] = 6.91865;
1663         parDeuteron[2] = 0.0183501;
1664         parDeuteron[3] = 2.49437;
1665         parDeuteron[4] = 2.62616;
1666       }
1667
1668
1669       //define expected signals for the various species
1670       Double_t expSignalDeuteronN = 0;
1671       Double_t expSignalDeuteronP = 0;
1672
1673       isDeuteron[0] = kFALSE;
1674       isDeuteron[1] = kFALSE;
1675       isDeuteron[2] = kFALSE;
1676
1677       //for data
1678       if(!fMCtrue){
1679
1680         expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1681         expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1682
1683         if(trackP->GetTPCsignal() >= 100 && //needed to reduce the size of the tree -- below the deuterons merge with the other particle spezies; standard value 110, variation for systematics
1684            trackP->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1685            (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.32 && //for systemactics -- normal value 0.24
1686            ptotP > 0.2 ){
1687
1688           if(trackP->Charge() >0)               isDeuteron[0] = kTRUE; //pos deuteron
1689           if(trackP->Charge() <0)               isDeuteron[1] = kTRUE; //neg deuteron
1690         }
1691
1692         if(trackN->GetTPCsignal() >= 100 && //needed to reduce the size of the tree -- below the deuterons merge with the other particle spezies
1693            trackN->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1694            (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.32 && //for systemactics -- normal value 0.24
1695            ptotN > 0.2){
1696
1697           isDeuteron[2] = kTRUE;
1698
1699           if(trackN->Charge() >0)        isDeuteron[0] = kTRUE; //pos deuteron
1700           if(trackN->Charge() <0)        isDeuteron[1] = kTRUE; //neg deuteron
1701         }
1702       }
1703
1704       //for MC
1705       if(fMCtrue)
1706         {
1707           if(runNumber < 166500) //2010
1708             {
1709               expSignalDeuteronN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1710               expSignalDeuteronP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1711             }
1712           if(runNumber > 166500) //2011
1713             {
1714               expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1715               expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1716             }
1717
1718           if(trackP->GetTPCsignal() < 1200 &&
1719              (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1720              ptotP > 0.2 )
1721             {
1722               if(trackP->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1723               if(trackP->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1724             }
1725
1726           if(trackN->GetTPCsignal() < 1200 &&
1727              (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1728              ptotN > 0.2 )
1729             {
1730               isDeuteron[2] = kTRUE;
1731
1732               if(trackN->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1733               if(trackN->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1734             }
1735
1736         }
1737
1738
1739
1740      return isDeuteron;
1741
1742 }
1743 //______________________________________________________________________________
1744 Bool_t AliAnalysisTaskLambdaNAOD::PionPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isPion[2]) {
1745
1746   //Pion PID via specific energy loss in the TPC
1747   //define the array for the Bethe-Bloch-Parameters (only necessary for local MC)
1748   Double_t parPion[5]  = {0,0,0,0,0};
1749
1750     if(runNumber < 166500) { //LHC10h
1751         parPion[0]   = 1.45802; // ALEPH parameters for pions (pass2)
1752         parPion[1]   = 27.4992;
1753         parPion[2]   = 4.00313e-15;
1754         parPion[3]   = 2.48485;
1755         parPion[4]   = 8.31768;
1756     }
1757
1758     if(runNumber > 166500) {  //LHC11h
1759         parPion[0]   = 1.11243; // ALEPH parameters for pions (pass2)
1760         parPion[1]   = 26.1144;
1761         parPion[2]   = 4.00313e-15;
1762         parPion[3]   = 2.72969;
1763         parPion[4]   = 9.15038;
1764     }
1765
1766     if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1767       parPion[0] = 1.064; // ALEPH parameters for deuterons (pass2)
1768       parPion[1] = 26.3;
1769       parPion[2] = 4.00313e-15;
1770       parPion[3] = 2.703 ;
1771       parPion[4] = 9.967;
1772     }
1773
1774     Double_t expSignalPionP = 0;
1775     Double_t expSignalPionN = 0;
1776
1777     //for MC
1778     if(fMCtrue){
1779       if(runNumber < 166500){ //2010
1780         expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1781         expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1782       }
1783       if(runNumber > 166500){ //2011
1784         expSignalPionP = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1785         expSignalPionN = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1786       }
1787
1788     }
1789
1790     isPion[0] = kFALSE;
1791     isPion[1] = kFALSE;
1792     //data
1793     if(!fMCtrue){
1794       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion))<4){ //for systematics -- normal value 3
1795
1796         if(trackP->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1797         if(trackP->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1798      }
1799       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion))<4){ //for systematics -- normal value 3
1800
1801         if(trackN->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1802         if(trackN->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1803       }
1804     }
1805
1806     //MC
1807     if(fMCtrue){
1808       if(TMath::Abs(trackP->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.2
1809          && ptotP>0.00001){
1810         if(trackP->Charge()>0) isPion[0] = kTRUE; //pos pion
1811         if(trackP->Charge()<0) isPion[1] = kTRUE; //neg pion
1812       }
1813      if(TMath::Abs(trackN->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.2
1814         && ptotN>0.00001){
1815         if(trackN->Charge()>0) isPion[0] = kTRUE; //pos pion
1816         if(trackN->Charge()<0) isPion[1] = kTRUE; //neg pion
1817       }
1818     }
1819
1820     return isPion;
1821
1822 }
1823 //______________________________________________________________________________
1824 Bool_t AliAnalysisTaskLambdaNAOD::TrackCuts(AliVTrack *track, Bool_t testTrackCuts) {
1825
1826   testTrackCuts = kFALSE;
1827
1828   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1829   //if(!esdtrack) testTrackCuts = kFALSE;
1830   if (fESDtrackCutsV0->AcceptTrack(esdtrack)) testTrackCuts = kTRUE;
1831 //if( testTrackCuts == kTRUE) cout <<   "testTrackCuts im test: " << testTrackCuts << endl;
1832
1833
1834   return testTrackCuts;
1835 }
1836 //______________________________________________________________________________
1837 /*Bool_t AliAnalysisTaskLambdaNAOD::FilterBit(AliVTrack *track, Bool_t testFilterBit) {
1838
1839   testFilterBit = kFALSE;
1840
1841   AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1842   //if(!aodtrack) testFilterBit = kFALSE;
1843   if (aodtrack->TestFilterBit(7)) testFilterBit = kTRUE;
1844
1845   //if(testFilterBit == kTRUE) cout <<   "testFilterBit im test: " << testFilterBit<< endl;
1846
1847   return testFilterBit;
1848   }*/
1849 //______________________________________________________________________
1850 void AliAnalysisTaskLambdaNAOD::MCGenerated(AliStack* stack, Int_t multiplicity)
1851 {
1852
1853   // Monte Carlo for genenerated particles
1854   if (fMCtrue) //MC loop
1855     {
1856
1857       Int_t stackN = 0;
1858
1859       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1860         {
1861
1862           const TParticle *tparticleMother = stack->Particle(stackN);
1863           Long_t pdgCodeMother = tparticleMother->GetPdgCode();
1864
1865           //check which particle the mother is
1866
1867           //LambdaNeutron
1868           if(pdgCodeMother == fgkPdgCode[kPDGLambdaNeutron]) //check mother PDG
1869             {
1870               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGDeuteron],fgkPdgCode[kPDGPionMinus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1871             }
1872
1873           //Anti-LambdaNeutron
1874           if(pdgCodeMother == fgkPdgCode[kPDGAntiLambdaNeutron]) //check mother PDG
1875             {
1876               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGAntiDeuteron],fgkPdgCode[kPDGPionPlus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1877             }
1878
1879
1880         }//end loop over stack
1881
1882
1883     }//end MC
1884 }
1885 //_____________________________________________
1886 void AliAnalysisTaskLambdaNAOD::MCTwoBodyDecay(AliStack* stack, const TParticle *tparticleMother, Long_t PDGMother, Long_t PDGFirstDaughter, Long_t PDGSecondDaughter, Double_t massFirstDaughter, Double_t massSecondDaughter, Int_t multiplicity){ //function that calculates the invariant mass and the transverse momentum for MC
1887
1888   Double_t momentumFirstDaughterGen[3]={0,0,0};
1889   Double_t momentumSecondDaughterGen[3]={0,0,0};
1890
1891   Double_t energyFirstDaughterGen = 0;
1892   Double_t energySecondDaughterGen = 0;
1893
1894   Double_t transversMomentumMotherGen = 0;
1895   Double_t longitudinalMomentumMotherGen = 0;
1896   Double_t totalEnergyMotherGen = 0;
1897
1898   Double_t rapidityGen = 2;
1899
1900   //Int_t labelFirstDaughter = tparticleMother->GetDaughter(0);
1901   Int_t labelSecondDaughter = tparticleMother->GetDaughter(1);
1902   Int_t labelFirstDaughter = labelSecondDaughter -1;
1903
1904   TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1905   TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1906
1907   if(tparticleFirstDaughter->GetPdgCode() == PDGFirstDaughter) //check first daughter PDG
1908     {
1909       if(tparticleSecondDaughter->GetPdgCode() == PDGSecondDaughter) //check second daughter PDG
1910         {
1911
1912           momentumFirstDaughterGen[0] = tparticleFirstDaughter->Px();
1913           momentumFirstDaughterGen[1] = tparticleFirstDaughter->Py();
1914           momentumFirstDaughterGen[2] = tparticleFirstDaughter->Pz();
1915
1916           momentumSecondDaughterGen[0] = tparticleSecondDaughter->Px();
1917           momentumSecondDaughterGen[1] = tparticleSecondDaughter->Py();
1918           momentumSecondDaughterGen[2] = tparticleSecondDaughter->Pz();
1919
1920           energyFirstDaughterGen  = tparticleFirstDaughter->Energy();
1921           energySecondDaughterGen = tparticleSecondDaughter->Energy();
1922
1923           transversMomentumMotherGen=TMath::Sqrt((momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])*(momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])+(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1])*(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1]));
1924
1925           //rapidity cut +- 1
1926           //longitudinal momentum of mother
1927           longitudinalMomentumMotherGen = TMath::Sqrt((momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2])*(momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2]));
1928
1929           //total energy of mother
1930           totalEnergyMotherGen = energyFirstDaughterGen + energySecondDaughterGen;
1931
1932           //rapidity
1933           rapidityGen = 0.5 * TMath::Log( (totalEnergyMotherGen+longitudinalMomentumMotherGen)/(totalEnergyMotherGen-longitudinalMomentumMotherGen));
1934
1935           if(rapidityGen > 1.0 || rapidityGen < -1 ) return;
1936
1937           //cout << "decay lenght first daughter :" << tparticleFirstDaughter->Rho()  << " decay lenght second daughter :" << tparticleSecondDaughter->Rho()<< endl;
1938           //fHistLambdaNeutronDecayLengthGen->Fill(tparticleFirstDaughter->Rho());
1939
1940           //calculate the invariant mass
1941           Double_t firstDaughterPart = 0;
1942           Double_t secondDaughterPart = 0;
1943           Double_t invaMass = 0;
1944
1945           firstDaughterPart = massFirstDaughter*massFirstDaughter
1946             +momentumFirstDaughterGen[0]*momentumFirstDaughterGen[0]
1947             +momentumFirstDaughterGen[1]*momentumFirstDaughterGen[1]
1948             +momentumFirstDaughterGen[2]*momentumFirstDaughterGen[2];
1949           secondDaughterPart = massSecondDaughter*massSecondDaughter
1950             +momentumSecondDaughterGen[0]*momentumSecondDaughterGen[0]
1951             +momentumSecondDaughterGen[1]*momentumSecondDaughterGen[1]
1952             +momentumSecondDaughterGen[2]*momentumSecondDaughterGen[2];
1953
1954           invaMass = TMath::Sqrt(TMath::Max(massFirstDaughter*massFirstDaughter
1955                                                         +massSecondDaughter*massSecondDaughter
1956                                                         + 2.*(TMath::Sqrt(firstDaughterPart*secondDaughterPart)
1957                                                               -momentumFirstDaughterGen[0]*momentumSecondDaughterGen[0]
1958                                                               -momentumFirstDaughterGen[1]*momentumSecondDaughterGen[1]
1959                                                               -momentumFirstDaughterGen[2]*momentumSecondDaughterGen[2]), 0.));
1960
1961 #if 0
1962
1963                                                               switch(PDGMother) {
1964                                                               case fgkPdgCode[kPDGLambdaNeutron] :
1965                                                                 fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1966                                                                 if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1967                                                                 if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1968                                                                 if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1969                                                                 fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1970                                                                 break;
1971                                                               case fgkPdgCode[kPDGAntiLambdaNeutron] :
1972                                                                 fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1973                                                                 if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1974                                                                 if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1975                                                                 if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1976                                                                 fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1977                                                                 break;
1978                                                               default :
1979                                                                 printf("should not happen!!!! \n");
1980                                                               }
1981
1982
1983 #else
1984           //LambdaNeutron
1985           if(PDGMother == fgkPdgCode[kPDGLambdaNeutron])
1986             {
1987               fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1988               if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1989               if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1990               if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1991               fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1992               fHistLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1993             }
1994           //Anti-LambdaNeutron
1995           if(PDGMother == fgkPdgCode[kPDGAntiLambdaNeutron])
1996             {
1997               fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1998               if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1999               if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
2000               if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
2001               fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
2002               fHistAntiLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
2003             }
2004
2005 #endif
2006         }//end of check second daughter PDG
2007     }//end of check first daughter PDG
2008
2009 }
2010 //_____________________________________________
2011 Double_t AliAnalysisTaskLambdaNAOD::MomentumInnerParam(AliVTrack *track, Double_t ptot){ //function to get the momentum at the inner wall of the TPC
2012
2013   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2014   ptot = esdtrack->GetInnerParam()->GetP();
2015
2016   return ptot;
2017 }
2018 //_____________________________________________
2019 UShort_t AliAnalysisTaskLambdaNAOD::TPCclusters(AliVTrack *track, UShort_t numberOfTPCclusters){ //function to get the number of clusters used for each track
2020
2021   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2022   numberOfTPCclusters = esdtrack->GetTPCNcls();
2023
2024   return numberOfTPCclusters;
2025 }
2026 //_____________________________________________
2027 Double_t AliAnalysisTaskLambdaNAOD::TPCchi2(AliVTrack *track, Double_t numberOfChi2clustersTPC, UShort_t numberOfTPCclusters){ //function to get the chi2 per clusters used for each track
2028
2029   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2030   numberOfChi2clustersTPC = esdtrack->GetTPCchi2()/numberOfTPCclusters;
2031
2032   return numberOfChi2clustersTPC;
2033 }
2034 //_____________________________________________
2035 Double_t AliAnalysisTaskLambdaNAOD::ImpactParameter(AliVTrack *track, Double_t dcaToVertex){ //function to get the number of clusters used for each track
2036
2037   Float_t tdcaToVertex[2] = {-1,-1};
2038
2039   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2040   esdtrack->GetImpactParameters(tdcaToVertex[0],tdcaToVertex[1]);
2041   dcaToVertex = TMath::Sqrt(tdcaToVertex[0]*tdcaToVertex[0]+tdcaToVertex[1]*tdcaToVertex[1]);
2042
2043   return dcaToVertex;
2044 }
2045 //________________________________________________________________________
2046 /*void AliAnalysisTaskLambdaNAOD::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
2047
2048      // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
2049      if (!kfParticle) return;
2050      Double_t dx = 0.;
2051      Double_t dy = 0.;
2052      Double_t dz = 0.;
2053
2054      if (ev){
2055        dx = ev->GetPrimaryVertex()->GetX()-0.;
2056        dy = ev->GetPrimaryVertex()->GetY()-0.;
2057        dz = ev->GetPrimaryVertex()->GetZ()-0.;
2058      }
2059
2060      kfParticle->X() = kfParticle->GetX() - dx;
2061      kfParticle->Y() = kfParticle->GetY() - dy;
2062      kfParticle->Z() = kfParticle->GetZ() - dz;
2063
2064
2065      // Rotate the kf particle
2066      Double_t c = cos(angle);
2067      Double_t s = sin(angle);
2068
2069      Double_t mA[8][ 8];
2070      for( Int_t i=0; i<8; i++ ){
2071        for( Int_t j=0; j<8; j++){
2072          mA[i][j] = 0;
2073        }
2074      }
2075      for( int i=0; i<8; i++ ){
2076        mA[i][i] = 1;
2077      }
2078      mA[0][0] =  c;  mA[0][1] = s;
2079      mA[1][0] = -s;  mA[1][1] = c;
2080      mA[3][3] =  c;  mA[3][4] = s;
2081      mA[4][3] = -s;  mA[4][4] = c;
2082
2083      Double_t mAC[8][8];
2084      Double_t mAp[8];
2085
2086      for( Int_t i=0; i<8; i++ ){
2087        mAp[i] = 0;
2088        for( Int_t k=0; k<8; k++){
2089          mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
2090        }
2091      }
2092
2093      for( Int_t i=0; i<8; i++){
2094        kfParticle->Parameter(i) = mAp[i];
2095      }
2096
2097      for( Int_t i=0; i<8; i++ ){
2098        for( Int_t j=0; j<8; j++ ){
2099          mAC[i][j] = 0;
2100          for( Int_t k=0; k<8; k++ ){
2101            mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
2102          }
2103        }
2104      }
2105
2106      for( Int_t i=0; i<8; i++ ){
2107        for( Int_t j=0; j<=i; j++ ){
2108          Double_t xx = 0;
2109          for( Int_t k=0; k<8; k++){
2110            xx+= mAC[i][k]*mA[j][k];
2111          }
2112          kfParticle->Covariance(i,j) = xx;
2113        }
2114      }
2115
2116      kfParticle->X() = kfParticle->GetX() + dx;
2117      kfParticle->Y() = kfParticle->GetY() + dy;
2118      kfParticle->Z() = kfParticle->GetZ() + dz;
2119
2120      }*/