PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskLambdaNAOD.cxx

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