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