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]/I");
 370   fTreeV0->Branch("fSigmadEdxPionNeg",fSigmadEdxPionNeg,"fSigmadEdxPionNeg[fItrk]/I");
 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/D");
 386   fTreeV0->Branch("fImpactParameterDeuteronNeg",fImpactParameterDeuteronNeg,"fImpactParameterDeuteronNeg/D");
 387   fTreeV0->Branch("fImpactParameterPionPos",fImpactParameterPionPos,"fImpactParameterPionPos/D");
 388   fTreeV0->Branch("fImpactParameterPionNeg",fImpactParameterPionNeg,"fImpactParameterPionNeg/D");
 389
 390   fTreeV0->Branch("fImpactParameterDeuteronPosAliKF",fImpactParameterDeuteronPosAliKF,"fImpactParameterDeuteronPosAliKF/D");
 391   fTreeV0->Branch("fImpactParameterDeuteronNegAliKF",fImpactParameterDeuteronNegAliKF,"fImpactParameterDeuteronNegAliKF/D");
 392   fTreeV0->Branch("fImpactParameterPionPosAliKF",fImpactParameterPionPosAliKF,"fImpactParameterPionPosAliKF/D");
 393   fTreeV0->Branch("fImpactParameterPionNegAliKF",fImpactParameterPionNegAliKF,"fImpactParameterPionNegAliKF/D");
 394
 395   fTreeV0->Branch("fMinNClustersTPCPos",fMinNClustersTPCPos,"fMinNClustersTPCPos/I");
 396   fTreeV0->Branch("fMinNClustersTPCNeg",fMinNClustersTPCNeg,"fMinNClustersTPCNeg/I");
 397
 398   fTreeV0->Branch("fMaxChi2PerClusterTPCPos",fMaxChi2PerClusterTPCPos,"fMaxChi2PerClusterTPCPos/F");
 399   fTreeV0->Branch("fMaxChi2PerClusterTPCNeg",fMaxChi2PerClusterTPCNeg,"fMaxChi2PerClusterTPCNeg/F");
 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]            = 0;
 794    fImpactParameterDeuteronNeg[fItrk]            = 0;
 795
 796    fImpactParameterPionPos[fItrk]            = 0;
 797    fImpactParameterPionNeg[fItrk]            = 0;
 798
 799    fImpactParameterDeuteronPosAliKF[fItrk]            = 0;
 800    fImpactParameterDeuteronNegAliKF[fItrk]            = 0;
 801
 802    fImpactParameterPionPosAliKF[fItrk]            = 0;
 803    fImpactParameterPionNegAliKF[fItrk]            = 0;
 804
 805    fMinNClustersTPCPos[fItrk]            = 0;
 806    fMinNClustersTPCNeg[fItrk]            = 0;
 807
 808    fMaxChi2PerClusterTPCPos[fItrk]            = 0;
 809    fMaxChi2PerClusterTPCNeg[fItrk]            = 0;
 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         Int_t numberOfTPCclustersPos = 0.;
1366         Int_t numberOfTPCclustersNeg = 0.;
1367
1368         numberOfTPCclustersPos = TPCclusters(trackP,numberOfTPCclustersPos);
1369         numberOfTPCclustersNeg = TPCclusters(trackN,numberOfTPCclustersNeg);
1370
1371         Float_t numberOfChi2clustersTPCPos = 10;
1372         Float_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           fMaxChi2PerClusterTPCPos[fItrk]            = numberOfChi2clustersTPCPos;
1526           fMaxChi2PerClusterTPCNeg[fItrk]            = numberOfChi2clustersTPCNeg;
1527
1528           fItrk++;
1529         }
1530
1531       }//end rotation loop
1532
1533  } //end of v0 loop
1534
1535   //fill the tree
1536   fTreeV0->Fill();
1537
1538   // Post output data.
1539   PostData(1, fOutputContainer);
1540   PostData(2, fTreeV0);
1541 }
1542
1543 //________________________________________________________________________
1544 void AliAnalysisTaskLambdaNAOD::Terminate(const Option_t *)
1545 {
1546   // Draw result to the screen
1547   // Called once at the end of the query
1548
1549   //get output data and darw 'fHistPt'
1550   if (!GetOutputData(0)) return;
1551   //TH1F *hist=(TH1F*)(((TObjArray*)GetOutputData(0))->FindObject("fHistPt"));
1552   //if (hist) hist->Draw();
1553 }
1554
1555 //_____________________________________________________________________________
1556 Int_t AliAnalysisTaskLambdaNAOD::Initialize() {
1557
1558
1559   // -- Reset Event
1560   // ----------------
1561   ResetEvent();
1562
1563   return 0;
1564   }
1565 //________________________________________________________________________
1566 Int_t AliAnalysisTaskLambdaNAOD::SetupEvent() {
1567   // Setup Reading of event
1568
1569   ResetEvent();
1570
1571   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1572   // -- Get Trigger
1573   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
1574
1575   //Bool_t isTriggered = IsTriggered();
1576   IsTriggered();
1577   return 0;
1578 }
1579 //_____________________________________________________________________________
1580 void AliAnalysisTaskLambdaNAOD::ResetEvent() {
1581
1582   // Reset event
1583
1584   // -- Reset QA
1585   for (Int_t ii = 0; ii < fNTriggers; ++ii)
1586     fTriggerFired[ii] = kFALSE;
1587
1588   return;
1589 }
1590 //_______________________________________________________________________
1591 void AliAnalysisTaskLambdaNAOD::BinLogAxis(const THnSparse *h, Int_t axisNumber) {
1592   //
1593   // Method for the correct logarithmic binning of histograms
1594   //
1595   TAxis *axis = h->GetAxis(axisNumber);
1596   int bins = axis->GetNbins();
1597
1598   Double_t from = axis->GetXmin();
1599   Double_t to = axis->GetXmax();
1600   Double_t *newBins = new Double_t[bins + 1];
1601
1602   newBins[0] = from;
1603   Double_t factor = pow(to/from, 1./bins);
1604
1605   for (int i = 1; i <= bins; i++) {
1606    newBins[i] = factor * newBins[i-1];
1607   }
1608   axis->Set(bins, newBins);
1609   delete [] newBins;
1610
1611   }
1612 //________________________________________________________________________
1613 Bool_t AliAnalysisTaskLambdaNAOD::IsTriggered() {
1614   // Check if Event is triggered and fill Trigger Histogram
1615
1616   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
1617   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
1618   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
1619   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
1620   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
1621
1622   Bool_t isTriggered = kFALSE;
1623
1624   for (Int_t ii=0; ii<fNTriggers; ++ii) {
1625     if(fTriggerFired[ii]) {
1626       isTriggered = kTRUE;
1627       fHistTriggerStat->Fill(ii);
1628     }
1629   }
1630
1631   return isTriggered;
1632   }
1633 //______________________________________________________________________________
1634 Bool_t AliAnalysisTaskLambdaNAOD::DeuteronPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isDeuteron[3]) {
1635
1636      //define the arrays for the Bethe-Bloch-Parameters
1637       Double_t parDeuteron[5] = {0,0,0,0,0};
1638
1639       if(runNumber < 166500) //LHC10h
1640         {
1641           parDeuteron[0] = 1.45802; // ALEPH parameters for deuterons (pass2)
1642           parDeuteron[1] = 27.4992;
1643           parDeuteron[2] = 4.00313e-15;
1644           parDeuteron[3] = 2.48485;
1645           parDeuteron[4] = 8.31768;
1646         }
1647
1648       if(runNumber > 166500) //LHC11h
1649         {
1650           parDeuteron[0] = 1.11243; // ALEPH parameters for deuterons (pass2)
1651           parDeuteron[1] = 26.1144;
1652           parDeuteron[2] = 4.00313e-15;
1653           parDeuteron[3] = 2.72969 ;
1654           parDeuteron[4] = 9.15038;
1655         }
1656
1657       if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1658         parDeuteron[0] = 1.064; // ALEPH parameters for deuterons (pass2)
1659         parDeuteron[1] = 26.3;
1660         parDeuteron[2] = 4.00313e-15;
1661         parDeuteron[3] = 2.703 ;
1662         parDeuteron[4] = 9.967;
1663       }
1664
1665
1666       //define expected signals for the various species
1667       Double_t expSignalDeuteronN = 0;
1668       Double_t expSignalDeuteronP = 0;
1669
1670       isDeuteron[0] = kFALSE;
1671       isDeuteron[1] = kFALSE;
1672       isDeuteron[2] = kFALSE;
1673
1674       //for data
1675       if(!fMCtrue){
1676
1677         expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1678         expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1679
1680         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
1681            trackP->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1682            (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.3 &&
1683            ptotP > 0.2 ){
1684
1685           if(trackP->Charge() >0)               isDeuteron[0] = kTRUE; //pos deuteron
1686           if(trackP->Charge() <0)               isDeuteron[1] = kTRUE; //neg deuteron
1687         }
1688
1689         if(trackN->GetTPCsignal() >= 100 && //needed to reduce the size of the tree -- below the deuterons merge with the other particle spezies
1690            trackN->GetTPCsignal() < 1200 && //needed to reduce the size of the tree -- anyway above the TPC does not work properly (?)
1691            (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.3 &&
1692            ptotN > 0.2){
1693
1694           isDeuteron[2] = kTRUE;
1695
1696           if(trackN->Charge() >0)        isDeuteron[0] = kTRUE; //pos deuteron
1697           if(trackN->Charge() <0)        isDeuteron[1] = kTRUE; //neg deuteron
1698         }
1699       }
1700
1701       //for MC
1702       if(fMCtrue)
1703         {
1704           if(runNumber < 166500) //2010
1705             {
1706               expSignalDeuteronN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1707               expSignalDeuteronP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1708             }
1709           if(runNumber > 166500) //2011
1710             {
1711               expSignalDeuteronN = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1712               expSignalDeuteronP = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1713             }
1714
1715           if(trackP->GetTPCsignal() < 1200 &&
1716              (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.5 &&
1717              ptotP > 0.2 )
1718             {
1719               if(trackP->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1720               if(trackP->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1721             }
1722
1723           if(trackN->GetTPCsignal() < 1200 &&
1724              (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.5 &&
1725              ptotN > 0.2 )
1726             {
1727               isDeuteron[2] = kTRUE;
1728
1729               if(trackN->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1730               if(trackN->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron
1731             }
1732
1733         }
1734
1735
1736
1737      return isDeuteron;
1738
1739 }
1740 //______________________________________________________________________________
1741 Bool_t AliAnalysisTaskLambdaNAOD::PionPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isPion[2]) {
1742
1743   //Pion PID via specific energy loss in the TPC
1744   //define the array for the Bethe-Bloch-Parameters (only necessary for local MC)
1745   Double_t parPion[5]  = {0,0,0,0,0};
1746
1747     if(runNumber < 166500) { //LHC10h
1748         parPion[0]   = 1.45802; // ALEPH parameters for pions (pass2)
1749         parPion[1]   = 27.4992;
1750         parPion[2]   = 4.00313e-15;
1751         parPion[3]   = 2.48485;
1752         parPion[4]   = 8.31768;
1753     }
1754
1755     if(runNumber > 166500) {  //LHC11h
1756         parPion[0]   = 1.11243; // ALEPH parameters for pions (pass2)
1757         parPion[1]   = 26.1144;
1758         parPion[2]   = 4.00313e-15;
1759         parPion[3]   = 2.72969;
1760         parPion[4]   = 9.15038;
1761     }
1762
1763     if(fMCtrue && runNumber > 166500){ //LHC11h MC (local production)
1764       parPion[0] = 1.064; // ALEPH parameters for deuterons (pass2)
1765       parPion[1] = 26.3;
1766       parPion[2] = 4.00313e-15;
1767       parPion[3] = 2.703 ;
1768       parPion[4] = 9.967;
1769     }
1770
1771     Double_t expSignalPionP = 0;
1772     Double_t expSignalPionN = 0;
1773
1774     //for MC
1775     if(fMCtrue){
1776       if(runNumber < 166500){ //2010
1777         expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1778         expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1779       }
1780       if(runNumber > 166500){ //2011
1781         expSignalPionP = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1782         expSignalPionN = 1.1*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1783       }
1784
1785     }
1786
1787     isPion[0] = kFALSE;
1788     isPion[1] = kFALSE;
1789     //data
1790     if(!fMCtrue){
1791       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion))<4){
1792
1793         if(trackP->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1794         if(trackP->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1795      }
1796       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion))<4){
1797
1798         if(trackN->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1799         if(trackN->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1800       }
1801     }
1802
1803     //MC
1804     if(fMCtrue){
1805       if(TMath::Abs(trackP->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.5
1806          && ptotP>0.00001){
1807         if(trackP->Charge()>0) isPion[0] = kTRUE; //pos pion
1808         if(trackP->Charge()<0) isPion[1] = kTRUE; //neg pion
1809       }
1810      if(TMath::Abs(trackN->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.5
1811         && ptotN>0.00001){
1812         if(trackN->Charge()>0) isPion[0] = kTRUE; //pos pion
1813         if(trackN->Charge()<0) isPion[1] = kTRUE; //neg pion
1814       }
1815     }
1816
1817     return isPion;
1818
1819 }
1820 //______________________________________________________________________________
1821 Bool_t AliAnalysisTaskLambdaNAOD::TrackCuts(AliVTrack *track, Bool_t testTrackCuts) {
1822
1823   testTrackCuts = kFALSE;
1824
1825   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1826   //if(!esdtrack) testTrackCuts = kFALSE;
1827   if (fESDtrackCutsV0->AcceptTrack(esdtrack)) testTrackCuts = kTRUE;
1828 //if( testTrackCuts == kTRUE) cout <<   "testTrackCuts im test: " << testTrackCuts << endl;
1829
1830
1831   return testTrackCuts;
1832 }
1833 //______________________________________________________________________________
1834 /*Bool_t AliAnalysisTaskLambdaNAOD::FilterBit(AliVTrack *track, Bool_t testFilterBit) {
1835
1836   testFilterBit = kFALSE;
1837
1838   AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1839   //if(!aodtrack) testFilterBit = kFALSE;
1840   if (aodtrack->TestFilterBit(7)) testFilterBit = kTRUE;
1841
1842   //if(testFilterBit == kTRUE) cout <<   "testFilterBit im test: " << testFilterBit<< endl;
1843
1844   return testFilterBit;
1845   }*/
1846 //______________________________________________________________________
1847 void AliAnalysisTaskLambdaNAOD::MCGenerated(AliStack* stack, Int_t multiplicity)
1848 {
1849
1850   // Monte Carlo for genenerated particles
1851   if (fMCtrue) //MC loop
1852     {
1853
1854       Int_t stackN = 0;
1855
1856       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1857         {
1858
1859           const TParticle *tparticleMother = stack->Particle(stackN);
1860           Long_t pdgCodeMother = tparticleMother->GetPdgCode();
1861
1862           //check which particle the mother is
1863
1864           //LambdaNeutron
1865           if(pdgCodeMother == fgkPdgCode[kPDGLambdaNeutron]) //check mother PDG
1866             {
1867               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGDeuteron],fgkPdgCode[kPDGPionMinus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1868             }
1869
1870           //Anti-LambdaNeutron
1871           if(pdgCodeMother == fgkPdgCode[kPDGAntiLambdaNeutron]) //check mother PDG
1872             {
1873               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGAntiDeuteron],fgkPdgCode[kPDGPionPlus],fgkMass[kMassDeuteron],fgkMass[kMassPion],multiplicity);
1874             }
1875
1876
1877         }//end loop over stack
1878
1879
1880     }//end MC
1881 }
1882 //_____________________________________________
1883 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
1884
1885   Double_t momentumFirstDaughterGen[3]={0,0,0};
1886   Double_t momentumSecondDaughterGen[3]={0,0,0};
1887
1888   Double_t energyFirstDaughterGen = 0;
1889   Double_t energySecondDaughterGen = 0;
1890
1891   Double_t transversMomentumMotherGen = 0;
1892   Double_t longitudinalMomentumMotherGen = 0;
1893   Double_t totalEnergyMotherGen = 0;
1894
1895   Double_t rapidityGen = 2;
1896
1897   //Int_t labelFirstDaughter = tparticleMother->GetDaughter(0);
1898   Int_t labelSecondDaughter = tparticleMother->GetDaughter(1);
1899   Int_t labelFirstDaughter = labelSecondDaughter -1;
1900
1901   TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1902   TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1903
1904   if(tparticleFirstDaughter->GetPdgCode() == PDGFirstDaughter) //check first daughter PDG
1905     {
1906       if(tparticleSecondDaughter->GetPdgCode() == PDGSecondDaughter) //check second daughter PDG
1907         {
1908
1909           momentumFirstDaughterGen[0] = tparticleFirstDaughter->Px();
1910           momentumFirstDaughterGen[1] = tparticleFirstDaughter->Py();
1911           momentumFirstDaughterGen[2] = tparticleFirstDaughter->Pz();
1912
1913           momentumSecondDaughterGen[0] = tparticleSecondDaughter->Px();
1914           momentumSecondDaughterGen[1] = tparticleSecondDaughter->Py();
1915           momentumSecondDaughterGen[2] = tparticleSecondDaughter->Pz();
1916
1917           energyFirstDaughterGen  = tparticleFirstDaughter->Energy();
1918           energySecondDaughterGen = tparticleSecondDaughter->Energy();
1919
1920           transversMomentumMotherGen=TMath::Sqrt((momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])*(momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])+(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1])*(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1]));
1921
1922           //rapidity cut +- 1
1923           //longitudinal momentum of mother
1924           longitudinalMomentumMotherGen = TMath::Sqrt((momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2])*(momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2]));
1925
1926           //total energy of mother
1927           totalEnergyMotherGen = energyFirstDaughterGen + energySecondDaughterGen;
1928
1929           //rapidity
1930           rapidityGen = 0.5 * TMath::Log( (totalEnergyMotherGen+longitudinalMomentumMotherGen)/(totalEnergyMotherGen-longitudinalMomentumMotherGen));
1931
1932           if(rapidityGen > 1.0 || rapidityGen < -1 ) return;
1933
1934           //cout << "decay lenght first daughter :" << tparticleFirstDaughter->Rho()  << " decay lenght second daughter :" << tparticleSecondDaughter->Rho()<< endl;
1935           //fHistLambdaNeutronDecayLengthGen->Fill(tparticleFirstDaughter->Rho());
1936
1937           //calculate the invariant mass
1938           Double_t firstDaughterPart = 0;
1939           Double_t secondDaughterPart = 0;
1940           Double_t invaMass = 0;
1941
1942           firstDaughterPart = massFirstDaughter*massFirstDaughter
1943             +momentumFirstDaughterGen[0]*momentumFirstDaughterGen[0]
1944             +momentumFirstDaughterGen[1]*momentumFirstDaughterGen[1]
1945             +momentumFirstDaughterGen[2]*momentumFirstDaughterGen[2];
1946           secondDaughterPart = massSecondDaughter*massSecondDaughter
1947             +momentumSecondDaughterGen[0]*momentumSecondDaughterGen[0]
1948             +momentumSecondDaughterGen[1]*momentumSecondDaughterGen[1]
1949             +momentumSecondDaughterGen[2]*momentumSecondDaughterGen[2];
1950
1951           invaMass = TMath::Sqrt(TMath::Max(massFirstDaughter*massFirstDaughter
1952                                                         +massSecondDaughter*massSecondDaughter
1953                                                         + 2.*(TMath::Sqrt(firstDaughterPart*secondDaughterPart)
1954                                                               -momentumFirstDaughterGen[0]*momentumSecondDaughterGen[0]
1955                                                               -momentumFirstDaughterGen[1]*momentumSecondDaughterGen[1]
1956                                                               -momentumFirstDaughterGen[2]*momentumSecondDaughterGen[2]), 0.));
1957
1958 #if 0
1959
1960                                                               switch(PDGMother) {
1961                                                               case fgkPdgCode[kPDGLambdaNeutron] :
1962                                                                 fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1963                                                                 if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1964                                                                 if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1965                                                                 if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1966                                                                 fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1967                                                                 break;
1968                                                               case fgkPdgCode[kPDGAntiLambdaNeutron] :
1969                                                                 fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1970                                                                 if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1971                                                                 if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1972                                                                 if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1973                                                                 fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1974                                                                 break;
1975                                                               default :
1976                                                                 printf("should not happen!!!! \n");
1977                                                               }
1978
1979
1980 #else
1981           //LambdaNeutron
1982           if(PDGMother == fgkPdgCode[kPDGLambdaNeutron])
1983             {
1984               fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1985               if(multiplicity < 2500)fHistLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1986               if(multiplicity > 1500 && multiplicity < 2750)fHistLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1987               if(multiplicity > 300 && multiplicity < 2000)fHistLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1988               fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1989               fHistLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1990             }
1991           //Anti-LambdaNeutron
1992           if(PDGMother == fgkPdgCode[kPDGAntiLambdaNeutron])
1993             {
1994               fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1995               if(multiplicity < 2500)fHistAntiLambdaNeutronPtGenMinBias->Fill(transversMomentumMotherGen);
1996               if(multiplicity > 1500 && multiplicity < 2750)fHistAntiLambdaNeutronPtGenCentral->Fill(transversMomentumMotherGen);
1997               if(multiplicity > 300 && multiplicity < 2000)fHistAntiLambdaNeutronPtGenSemiCentral->Fill(transversMomentumMotherGen);
1998               fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1999               fHistAntiLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
2000             }
2001
2002 #endif
2003         }//end of check second daughter PDG
2004     }//end of check first daughter PDG
2005
2006 }
2007 //_____________________________________________
2008 Double_t AliAnalysisTaskLambdaNAOD::MomentumInnerParam(AliVTrack *track, Double_t ptot){ //function to get the momentum at the inner wall of the TPC
2009
2010   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2011   ptot = esdtrack->GetInnerParam()->GetP();
2012
2013   return ptot;
2014 }
2015 //_____________________________________________
2016 Int_t AliAnalysisTaskLambdaNAOD::TPCclusters(AliVTrack *track, Int_t numberOfTPCclusters){ //function to get the number of clusters used for each track
2017
2018   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2019   numberOfTPCclusters = esdtrack->GetTPCNcls();
2020
2021   return numberOfTPCclusters;
2022 }
2023 //_____________________________________________
2024 Int_t AliAnalysisTaskLambdaNAOD::TPCchi2(AliVTrack *track, Float_t numberOfChi2clustersTPC, Int_t numberOfTPCclusters){ //function to get the chi2 per clusters used for each track
2025
2026   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2027   numberOfChi2clustersTPC = esdtrack->GetTPCchi2()/numberOfTPCclusters;
2028
2029   return numberOfChi2clustersTPC;
2030 }
2031 //_____________________________________________
2032 Double_t AliAnalysisTaskLambdaNAOD::ImpactParameter(AliVTrack *track, Double_t dcaToVertex){ //function to get the number of clusters used for each track
2033
2034   Float_t tdcaToVertex[2] = {-1,-1};
2035
2036   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
2037   esdtrack->GetImpactParameters(tdcaToVertex[0],tdcaToVertex[1]);
2038   dcaToVertex = TMath::Sqrt(tdcaToVertex[0]*tdcaToVertex[0]+tdcaToVertex[1]*tdcaToVertex[1]);
2039
2040   return dcaToVertex;
2041 }
2042 //________________________________________________________________________
2043 /*void AliAnalysisTaskLambdaNAOD::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
2044
2045      // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
2046      if (!kfParticle) return;
2047      Double_t dx = 0.;
2048      Double_t dy = 0.;
2049      Double_t dz = 0.;
2050
2051      if (ev){
2052        dx = ev->GetPrimaryVertex()->GetX()-0.;
2053        dy = ev->GetPrimaryVertex()->GetY()-0.;
2054        dz = ev->GetPrimaryVertex()->GetZ()-0.;
2055      }
2056
2057      kfParticle->X() = kfParticle->GetX() - dx;
2058      kfParticle->Y() = kfParticle->GetY() - dy;
2059      kfParticle->Z() = kfParticle->GetZ() - dz;
2060
2061
2062      // Rotate the kf particle
2063      Double_t c = cos(angle);
2064      Double_t s = sin(angle);
2065
2066      Double_t mA[8][ 8];
2067      for( Int_t i=0; i<8; i++ ){
2068        for( Int_t j=0; j<8; j++){
2069          mA[i][j] = 0;
2070        }
2071      }
2072      for( int i=0; i<8; i++ ){
2073        mA[i][i] = 1;
2074      }
2075      mA[0][0] =  c;  mA[0][1] = s;
2076      mA[1][0] = -s;  mA[1][1] = c;
2077      mA[3][3] =  c;  mA[3][4] = s;
2078      mA[4][3] = -s;  mA[4][4] = c;
2079
2080      Double_t mAC[8][8];
2081      Double_t mAp[8];
2082
2083      for( Int_t i=0; i<8; i++ ){
2084        mAp[i] = 0;
2085        for( Int_t k=0; k<8; k++){
2086          mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
2087        }
2088      }
2089
2090      for( Int_t i=0; i<8; i++){
2091        kfParticle->Parameter(i) = mAp[i];
2092      }
2093
2094      for( Int_t i=0; i<8; i++ ){
2095        for( Int_t j=0; j<8; j++ ){
2096          mAC[i][j] = 0;
2097          for( Int_t k=0; k<8; k++ ){
2098            mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
2099          }
2100        }
2101      }
2102
2103      for( Int_t i=0; i<8; i++ ){
2104        for( Int_t j=0; j<=i; j++ ){
2105          Double_t xx = 0;
2106          for( Int_t k=0; k<8; k++){
2107            xx+= mAC[i][k]*mA[j][k];
2108          }
2109          kfParticle->Covariance(i,j) = xx;
2110        }
2111      }
2112
2113      kfParticle->X() = kfParticle->GetX() + dx;
2114      kfParticle->Y() = kfParticle->GetY() + dy;
2115      kfParticle->Z() = kfParticle->GetZ() + dz;
2116
2117      }*/