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