reduced size of the tree
[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   
450   if (SetupEvent() < 0) {
451     PostData(1, fOutputContainer);
452     return;
453     }
454   
455   //DATA
456   AliESDEvent *fESDevent = 0x0;
457   AliAODEvent *fAODevent = 0x0;
458   
459   
460   if(!fMCtrue){
461     if(!fPIDResponse) {
462       AliError("Cannot get pid response");
463       return;
464     }
465   }
466  
467
468   Int_t centrality = -1;
469   Double_t vertex[3]          = {-100.0, -100.0, -100.0};
470
471   //Initialisation of the event and basic event cuts:
472   //1.) ESDs: 
473   if (fAnalysisType == "ESD") {
474         
475     fESDevent=dynamic_cast<AliESDEvent*>(InputEvent());
476     if (!fESDevent) {
477       AliWarning("ERROR: fESDevent not available \n");
478       return;
479     }    
480     
481     //Basic event cuts: 
482     //1.) vertex existence
483     const AliESDVertex *vertexESD = fESDevent->GetPrimaryVertexTracks();
484     if (vertexESD->GetNContributors()<1) 
485       {
486         // SPD vertex
487         vertexESD = fESDevent->GetPrimaryVertexSPD();
488         if(vertexESD->GetNContributors()<1) {
489           PostData(1, fOutputContainer);
490           return;
491         }  
492       }
493
494     vertexESD->GetXYZ(vertex);
495
496     //2. vertex position within 10 cm
497     if (TMath::Abs(vertexESD->GetZv()) > 10) return;
498     
499     //centrality selection in PbPb
500     if (fESDevent->GetEventSpecie() == 4) //species == 4 == PbPb
501       { // PbPb
502         AliCentrality *esdCentrality = fESDevent->GetCentrality();
503         centrality = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
504         if(!fMCtrue){
505           if (centrality < 0. || centrality > 8. ) return; //select only events with centralities between 0 and 80 %
506         }
507       }
508
509     //cout << "centrality "<< centrality << endl;
510     // count number of events
511     fHistNumberOfEvents->Fill(centrality);
512
513  }
514
515 //2.) AODs: 
516  else if (fAnalysisType == "AOD") {
517    
518    fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
519    if (!fAODevent) {
520      AliWarning("ERROR: lAODevent not available \n");
521      return;
522    }
523       
524    const AliAODVertex *vertexAOD = fAODevent->GetPrimaryVertex();
525     if (vertexAOD->GetNContributors()<1) 
526       {
527         // SPD vertex
528         vertexAOD = fAODevent->GetPrimaryVertexSPD();
529         if(vertexAOD->GetNContributors()<1) {
530           PostData(1, fOutputContainer);
531           return;
532         }  
533       }
534     vertexAOD->GetXYZ(vertex);
535
536     //2. vertex position within 10 cm
537     if (TMath::Abs(vertex[2]) > 10) return;
538    
539     //centrality selection
540    AliCentrality *aodCentrality = fAODevent->GetCentrality();
541    centrality = aodCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
542    if (centrality < 0 || centrality > 8) return; //select only events with centralities between 0 and 80 %
543     
544    // count number of events
545    fHistNumberOfEvents->Fill(centrality);
546
547  } else {
548    
549    Printf("Analysis type (ESD or AOD) not specified \n");
550    return;
551    
552  }
553  
554   
555 //v0 loop
556
557
558  fItrk = 0;
559
560  Int_t runNumber = 0;
561  runNumber = (InputEvent())->GetRunNumber();
562  
563  Int_t    nTrackMultiplicity             = -1;
564  nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
565   
566
567  for (Int_t ivertex=0; ivertex<(InputEvent()->GetNumberOfV0s()); ivertex++) { //beginn v0 loop
568          
569    AliESDv0 * v0ESD = 0x0;
570    AliAODv0 * v0AOD = 0x0;
571
572    Bool_t v0ChargesCorrect = kTRUE;
573    Bool_t testTrackCuts = kFALSE;
574    //Bool_t testFilterBit = kFALSE;
575
576    Int_t of = 7;
577    Int_t onFlyStatus = 5;
578    
579    Float_t dcaV0 = -1;
580    Float_t cosPointing= 2;
581    Float_t decayRadius= -1;
582    
583    AliVTrack * trackN = 0x0;
584    AliVTrack * trackP = 0x0;
585    
586    Double_t ptotN = -1000;
587    Double_t ptotP = -1000;
588       
589    Int_t chargeComboDeuteronPion = -1; 
590    
591    Double_t momentumPion[3]={0,0,0};
592    Double_t momentumPionRot[3]={0,0,0};
593    Double_t momentumDeuteron[3]={0,0,0};
594    Double_t momentumDeuteronRot[3]={0,0,0};
595    Double_t momentumMother[3] = {0,0,0};
596    
597    Double_t transversMomentumPion = 0;
598    Double_t transversMomentumDeuteron = 0;
599    
600    Double_t transversMomentumMother = 0;
601    Double_t longitudinalMomentumMother = 0;
602    
603    Double_t totalEnergyMother = 0;
604    Double_t energyPion = 0;
605    Double_t energyDeuteron = 0;
606    
607    Double_t rapidity = 2;
608    
609    TVector3 vecPion(0,0,0);
610    TVector3 vecDeuteron(0,0,0);
611    TVector3 vecMother(0,0,0);
612    
613    Double_t alpha = 2;
614    Double_t qt = -1;
615    
616    Double_t thetaPion = 0;
617    Double_t thetaDeuteron = 0;
618    
619    Double_t phi=0;
620    Double_t invaMassDeuteronPion = 0;
621    
622    Double_t e12 = 0;
623    Double_t r12   = 0;
624    Double_t d22 = 0;
625    Double_t dr22   = 0;
626
627    //Tree variables
628    fV0finder[fItrk]         = -1;
629    fkMB[fItrk]              = -1;
630    fkCentral[fItrk]         = -1;
631    fkSemiCentral[fItrk]     = -1;
632    //fkEMCEJE[fItrk]          = -1;
633    //fkEMCEGA[fItrk]          = -1;
634    
635    fPtotN[fItrk]            = -1000;
636    fPtotP[fItrk]            = -1000;
637    fMotherPt[fItrk]         = -1000;
638    
639    fdEdxN[fItrk]            = -1;
640    fdEdxP[fItrk]            = -1;
641    fSignN[fItrk]            = 0;
642    fSignP[fItrk]            = 0;
643    
644    fDCAv0[fItrk]            = -1;
645    fCosinePAv0[fItrk]       = -2;
646    fDecayRadiusTree[fItrk]  = -1;
647    
648    fInvaMassDeuteronPionTree[fItrk] = 0;
649    fChargeComboDeuteronPionTree[fItrk] = -1;
650    
651    //fAmenterosAlphaTree[fItrk] = 2;
652    //fAmenterosQtTree[fItrk] = -1;
653    
654    //Get v0 object
655    if(fAnalysisType == "ESD")v0ESD = fESDevent->GetV0(ivertex);
656    if(fAnalysisType == "AOD")v0AOD = fAODevent->GetV0(ivertex);
657    
658    
659    //distinguish between the two V0 finders: 0 offline V0 finder / 1 online V0 finder
660    
661    if(fAnalysisType == "ESD") of = v0ESD->GetOnFlyStatus();
662    if(fAnalysisType == "AOD") of = v0AOD->GetOnFlyStatus();
663    
664    if(of)  onFlyStatus= 1;
665    if(!of) onFlyStatus= 0;
666    
667    if(onFlyStatus==0)fof->Fill(1,0);
668    if(onFlyStatus==1)fof->Fill(1,1);
669    
670    //Get dca, cos of pointing angle and decay radius
671    
672    
673    if(fAnalysisType == "ESD")
674      { 
675        dcaV0 = v0ESD->GetDcaV0Daughters();
676        cosPointing = v0ESD->GetV0CosineOfPointingAngle();
677        decayRadius = v0ESD->GetRr();
678      }
679    
680    if(fAnalysisType == "AOD")
681      { 
682        dcaV0 = v0AOD->DcaV0Daughters();
683        cosPointing = v0AOD->CosPointingAngle(vertex);
684        decayRadius = v0AOD->DecayLengthV0(vertex);
685      }
686
687
688       // select coresponding tracks
689       if(fAnalysisType == "ESD")
690         { 
691           trackN = fESDevent->GetTrack(v0ESD->GetIndex(0));       
692           trackP = fESDevent->GetTrack(v0ESD->GetIndex(1));
693
694           if (trackN->Charge() > 0) // switch because of bug in V0 interface
695             { 
696               trackN = fESDevent->GetTrack(v0ESD->GetIndex(1));
697               trackP = fESDevent->GetTrack(v0ESD->GetIndex(0));
698               v0ChargesCorrect = kFALSE;
699             }
700           //Track-Cuts
701           testTrackCuts = TrackCuts(trackN,testTrackCuts);
702           if(testTrackCuts == kFALSE) continue;
703           testTrackCuts = TrackCuts(trackP,testTrackCuts);
704           if(testTrackCuts == kFALSE) continue;
705         }
706       
707       if(fAnalysisType == "AOD")
708         { 
709           trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
710           trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
711
712           if (trackN->Charge() > 0) // switch because of bug in V0 interface
713             { 
714               trackN = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(1));
715               trackP = dynamic_cast<AliVTrack*>(v0AOD->GetDaughter(0));
716               v0ChargesCorrect = kFALSE;
717             }
718           //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  
719           //testFilterBit = FilterBit(trackN,testFilterBit);
720           //if(testFilterBit == kFALSE) continue;
721           //testFilterBit = FilterBit(trackP,testFilterBit);
722           //if(testFilterBit == kFALSE) continue;       
723           //Track-Cuts
724           /*testTrackCuts = TrackCuts(trackN,testTrackCuts);
725           if(testTrackCuts == kFALSE) continue;
726           testTrackCuts = TrackCuts(trackP,testTrackCuts);
727           if(testTrackCuts == kFALSE) continue;*/
728         }
729       
730       //Get the total momentum for each track ---> at the inner readout of the TPC???? // momenta a always positive
731             if(fAnalysisType == "AOD") {
732         ptotN = trackN->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
733         ptotP = trackP->P(); //GetInnerParam()->GetP(); ******************FIXME: InnerParam
734       }
735
736       if(fAnalysisType == "ESD") {
737         ptotN =  MomentumInnerParam(trackN,ptotN); 
738         ptotP =  MomentumInnerParam(trackP,ptotP); 
739       }
740       
741       //fill QA dEdx with all V0 candidates      
742       if(trackP) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),0);
743       if(trackN) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),1);
744
745       if (dcaV0 > 3 || dcaV0 < 1e-20) continue;
746       
747       //Check how much the dca cut reduces the statistic (background) for the different VO-finder
748       if(onFlyStatus==0)fof->Fill(2,0);
749       if(onFlyStatus==1)fof->Fill(2,1);
750     
751       if (cosPointing < 0.9) continue;
752         
753       //Check how much the cos-of-the-pointing-angle-cut reduces the statistic (background) for the different VO-finder
754       if(onFlyStatus==0)fof->Fill(3,0);
755       if(onFlyStatus==1)fof->Fill(3,1);
756
757       //deuteron PID via specific energy loss in the TPC
758       Bool_t isDeuteron[3] = {kFALSE,kFALSE,kFALSE}; //0 = posDeuteron, 1 = negDeuteron, 2 = trackN is deuteron
759
760       DeuteronPID(trackP,trackN,ptotP,ptotN,runNumber,isDeuteron);
761       
762       //if(!isDeuteron[0] && !isDeuteron[1]) continue;
763       if((isDeuteron[0]==kFALSE) && (isDeuteron[1]==kFALSE)) continue;
764
765
766       //Check how much the nuclei PID cuts reduce the statistics (background) for the two V0-finders
767       if(onFlyStatus==0)fof->Fill(4,0);
768       if(onFlyStatus==1)fof->Fill(4,1);
769       
770       //Fill the QA dEdx with deuterons and helium3 after the nuclei PID cut
771       if(isDeuteron[1]) fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),2);
772       if(isDeuteron[0]) fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),3);
773
774       //deuteron PID via specific energy loss in the TPC
775       Bool_t isPion[2] = {kFALSE,kFALSE}; //0 = posPion, 1 = negPion
776       
777       PionPID(trackP,trackN,ptotP,ptotN,runNumber,isPion);
778
779       //if(isDeuteron[0] && !isPion[1]) continue; //pos deuteron and neg Pion
780       //if(isDeuteron[1] && !isPion[0]) continue; //neg deuteron and pos Pion
781
782       //Check how much the pion PID cuts reduce the statistics (background) for the two V0-finders
783       if(onFlyStatus==0)fof->Fill(5,0);
784       if(onFlyStatus==1)fof->Fill(5,1);
785
786
787       //Save the different charge combinations to differentiat between particles and anit-particles:
788       // -/+ = Anti-deuteron + positive Pion
789       // -/- = Anti-deuteron + negative Pion
790       // +/- = deuteron + negative Pion
791       // +/+ = deuteron + positive Pion
792
793       //Like-sign
794       if (trackN->Charge()<0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 1; // -/-
795       if (trackN->Charge()>0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 3; // +/+
796
797       //unlinke-sign
798       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+
799       //if (trackN->Charge()>0 && trackP->Charge()<0 && isDeuteron[1]==kTRUE && isPion[0]==kTRUE ) chargeComboDeuteronPion = 0; // -/+ //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; // +/- //dürfte wegen charge correctur am anfang nicht existiere
801       if (trackN->Charge()<0 && trackP->Charge()>0 && isDeuteron[0]==kTRUE && isPion[1]==kTRUE) chargeComboDeuteronPion = 2; // +/-
802
803       //if(chargeComboDeuteronPion==0) cout << "chargeN: " << trackN->Charge() << " chargeP: " << trackP->Charge() << endl;
804
805
806       //Fill the QA dEdx with all selected tracks after the PID cuts
807       fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),6);
808       fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),6);
809
810  
811             
812       //get the momenta of the daughters
813
814       if(fAnalysisType == "ESD")
815         {
816           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN gehört zu anti-deuteron, tackP gehört zu pion
817             
818             v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
819             v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
820             if (!v0ChargesCorrect) {
821
822               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
823               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
824             }
825           }
826
827           if(chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP gehört zu deuteron, tackN gehört zu pion
828             
829             v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
830             v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
831             if (!v0ChargesCorrect){ 
832
833               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
834               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
835             }
836           }
837
838           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN gehört zu deuteron oder pion, trackP gehört zu deuteron oder pion
839             if(isDeuteron[2]==kTRUE){ //trackN gehört zu deuteron, trackP gehört zum pion
840
841               v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
842               v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
843               if (!v0ChargesCorrect) {
844                 
845                 v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
846                 v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
847               }
848             }
849             if(isDeuteron[2]==kFALSE){ //trackP gehört zum deuteron, trackN gehört zum pion
850
851               v0ESD->GetNPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
852               v0ESD->GetPPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
853               if (!v0ChargesCorrect) {
854                 
855                 v0ESD->GetPPxPyPz(momentumPion[0], momentumPion[1], momentumPion[2]);
856                 v0ESD->GetNPxPyPz(momentumDeuteron[0], momentumDeuteron[1], momentumDeuteron[2]);
857               }
858             } 
859           }      
860
861  
862           //get the momenta of the mother
863           v0ESD->GetPxPyPz(momentumMother[0],momentumMother[1],momentumMother[2]);
864         }
865     
866
867       if(fAnalysisType == "AOD")
868         {
869           if(chargeComboDeuteronPion==0){ //anti-deuteron (isDeuteron[1]==kTRUE), positives pion (isPion[0]==kTRUE), trackN gehört zu anti-deuteron, tackP gehört zu pion 
870             
871             momentumPion[0] = v0AOD->MomPosX();
872             momentumPion[1] = v0AOD->MomPosY();
873             momentumPion[2] = v0AOD->MomPosZ();
874             
875             momentumDeuteron[0] = v0AOD->MomNegX();
876             momentumDeuteron[1] = v0AOD->MomNegY();
877             momentumDeuteron[2] = v0AOD->MomNegZ();
878             
879             if (!v0ChargesCorrect){ 
880
881               momentumPion[0] = v0AOD->MomNegX();
882               momentumPion[1] = v0AOD->MomNegY();
883               momentumPion[2] = v0AOD->MomNegZ();
884               
885               momentumDeuteron[0] = v0AOD->MomPosX();
886               momentumDeuteron[1] = v0AOD->MomPosY();
887               momentumDeuteron[2] = v0AOD->MomPosZ();
888             }
889           }
890           
891           if (chargeComboDeuteronPion==2){ //deuteron (isDeuteron[0]==kTRUE), negative pion (isPion[1]==kTRUE), trackP gehört zu deuteron, tackN gehört zu pion
892
893             momentumPion[0] = v0AOD->MomNegX();
894             momentumPion[1] = v0AOD->MomNegY();
895             momentumPion[2] = v0AOD->MomNegZ();
896             
897             momentumDeuteron[0] = v0AOD->MomPosX();
898             momentumDeuteron[1] = v0AOD->MomPosY();
899             momentumDeuteron[2] = v0AOD->MomPosZ();
900             
901             if (!v0ChargesCorrect){
902
903               momentumPion[0] = v0AOD->MomPosX();
904               momentumPion[1] = v0AOD->MomPosY();
905               momentumPion[2] = v0AOD->MomPosZ();
906               
907               momentumDeuteron[0] = v0AOD->MomNegX();
908               momentumDeuteron[1] = v0AOD->MomNegY();
909               momentumDeuteron[2] = v0AOD->MomNegZ();
910             }
911           }
912
913           if(chargeComboDeuteronPion==1 || chargeComboDeuteronPion==3){ //trackN gehört zu deuteron oder pion, trackP gehört zu deuteron oder pion
914             if(isDeuteron[2]==kTRUE){ //trackN gehört zu deuteron, trackP gehört zum pion
915           
916               momentumPion[0] = v0AOD->MomPosX();
917               momentumPion[1] = v0AOD->MomPosY();
918               momentumPion[2] = v0AOD->MomPosZ();
919               
920               momentumDeuteron[0] = v0AOD->MomNegX();
921               momentumDeuteron[1] = v0AOD->MomNegY();
922               momentumDeuteron[2] = v0AOD->MomNegZ();
923               
924               if (!v0ChargesCorrect){ 
925                 
926                 momentumPion[0] = v0AOD->MomNegX();
927                 momentumPion[1] = v0AOD->MomNegY();
928                 momentumPion[2] = v0AOD->MomNegZ();
929                 
930                 momentumDeuteron[0] = v0AOD->MomPosX();
931                 momentumDeuteron[1] = v0AOD->MomPosY();
932                 momentumDeuteron[2] = v0AOD->MomPosZ();
933               }
934             }
935            
936             if(isDeuteron[2]==kFALSE){ //trackP gehört zum deuteron, trackN gehört zum pion
937               
938               momentumPion[0] = v0AOD->MomNegX();
939               momentumPion[1] = v0AOD->MomNegY();
940               momentumPion[2] = v0AOD->MomNegZ();
941               
942               momentumDeuteron[0] = v0AOD->MomPosX();
943               momentumDeuteron[1] = v0AOD->MomPosY();
944               momentumDeuteron[2] = v0AOD->MomPosZ();
945               
946               if (!v0ChargesCorrect){
947                 
948                 momentumPion[0] = v0AOD->MomPosX();
949                 momentumPion[1] = v0AOD->MomPosY();
950                 momentumPion[2] = v0AOD->MomPosZ();
951                 
952                 momentumDeuteron[0] = v0AOD->MomNegX();
953                 momentumDeuteron[1] = v0AOD->MomNegY();
954                 momentumDeuteron[2] = v0AOD->MomNegZ();
955               }
956             }
957
958             //get the momenta of the mother
959             momentumMother[0] = v0AOD->MomV0X();
960             momentumMother[1] = v0AOD->MomV0Y();
961             momentumMother[2] = v0AOD->MomV0Z();
962           }
963         }
964       
965       //Rapidity - cut    
966       transversMomentumPion      = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]);
967       transversMomentumDeuteron  = TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]);
968   
969       transversMomentumMother=TMath::Sqrt((momentumDeuteron[0]+momentumPion[0])*(momentumDeuteron[0]+momentumPion[0])+(momentumDeuteron[1]+momentumPion[1])*(momentumDeuteron[1]+momentumPion[1]));
970       
971       longitudinalMomentumMother = TMath::Sqrt((momentumDeuteron[2]+momentumPion[2])*(momentumDeuteron[2]+momentumPion[2]));
972
973       energyDeuteron =  TMath::Sqrt(momentumDeuteron[0]*momentumDeuteron[0]+momentumDeuteron[1]*momentumDeuteron[1]+momentumDeuteron[2]*momentumDeuteron[2]+fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]);
974     
975       energyPion = TMath::Sqrt(momentumPion[0]*momentumPion[0]+momentumPion[1]*momentumPion[1]+momentumPion[2]*momentumPion[2]+fgkMass[kMassPion]*fgkMass[kMassPion]);
976     
977       totalEnergyMother = energyPion + energyDeuteron;
978                                     
979       //rapidity cut +- 1
980       rapidity = 0.5 * TMath::Log((totalEnergyMother+longitudinalMomentumMother)/( totalEnergyMother-longitudinalMomentumMother));
981     
982       if(rapidity > 1.0 || rapidity < -1) continue;
983
984     
985       //Armanteros-Podolanski
986       vecPion.SetXYZ(momentumPion[0],momentumPion[1],momentumPion[2]);
987       vecDeuteron.SetXYZ(momentumDeuteron[0],momentumDeuteron[1],momentumDeuteron[2]);
988       vecMother.SetXYZ(momentumMother[0],momentumMother[1],momentumMother[2]);
989   
990       thetaPion = acos((vecPion * vecMother)/(vecPion.Mag() * vecMother.Mag()));
991       thetaDeuteron = acos((vecDeuteron * vecMother)/(vecDeuteron.Mag() * vecMother.Mag()));
992
993       alpha = ((vecPion.Mag())*cos(thetaPion)-(vecDeuteron.Mag())*cos(thetaDeuteron))/
994         ((vecDeuteron.Mag())*cos(thetaDeuteron)+(vecPion.Mag())*cos(thetaPion)) ;
995       qt = vecDeuteron.Mag()*sin(thetaDeuteron);
996
997
998       //Rotation for background calculation
999       //Int_t rotation=1; // =1 signal, =2 Rotation of the pion , =3 Rotation of the deuteron
1000
1001       Double_t fStartAnglePhi=TMath::Pi();
1002       Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1003       phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1004      
1005       for(Int_t rotation=1;rotation<4;rotation++){ //loop for rotation
1006         
1007         //Double_t fStartAnglePhi=TMath::Pi();
1008         //Double_t fConeAnglePhi=TMath::Pi(); //-0.174;
1009         //phi  = fStartAnglePhi+(2*gRandom->Rndm()-1)*fConeAnglePhi;
1010           
1011         //calculate new rotated momenta
1012         momentumPionRot[0]=TMath::Cos(phi)*momentumPion[0]-TMath::Sin(phi)*momentumPion[1];
1013         momentumPionRot[1]=TMath::Sin(phi)*momentumPion[0]+TMath::Cos(phi)*momentumPion[1];
1014         
1015         momentumDeuteronRot[0]=TMath::Cos(phi)*momentumDeuteron[0]-TMath::Sin(phi)*momentumDeuteron[1];
1016         momentumDeuteronRot[1]=TMath::Sin(phi)*momentumDeuteron[0]+TMath::Cos(phi)*momentumDeuteron[1];
1017           
1018         //invariant mass calculations  
1019         fIsCorrectlyAssociated[fItrk] = kFALSE;
1020         
1021         //factor for the invariant mass calculation, which only include the pion
1022         e12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1023           +momentumPion[0]*momentumPion[0]
1024           +momentumPion[1]*momentumPion[1]
1025           +momentumPion[2]*momentumPion[2];
1026         
1027         r12   = fgkMass[kMassPion]*fgkMass[kMassPion]
1028           +momentumPionRot[0]*momentumPionRot[0]
1029           +momentumPionRot[1]*momentumPionRot[1]
1030           +momentumPion[2]*momentumPion[2];
1031         
1032         //factor for the invariant mass calculation, which only include the deuterons
1033         d22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1034           +momentumDeuteron[0]*momentumDeuteron[0]
1035           +momentumDeuteron[1]*momentumDeuteron[1]
1036           +momentumDeuteron[2]*momentumDeuteron[2];
1037         dr22   = fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1038           +momentumDeuteronRot[0]*momentumDeuteronRot[0]
1039           +momentumDeuteronRot[1]*momentumDeuteronRot[1]
1040           +momentumDeuteron[2]*momentumDeuteron[2];
1041
1042         if(rotation == 1){ //signal
1043           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1044                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron] 
1045                                                         + 2.*(TMath::Sqrt(e12*d22)
1046                                                               -momentumPion[0]*momentumDeuteron[0]
1047                                                               -momentumPion[1]*momentumDeuteron[1]
1048                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1049
1050           if (fMCtrue){
1051             
1052             Int_t labelN = 0;
1053             Int_t labelP = 0;
1054             labelN = trackN->GetLabel();
1055             labelP = trackP->GetLabel();
1056             
1057             TParticle *tparticleDaughterN = stack->Particle(TMath::Abs(labelN));
1058             TParticle *tparticleDaughterP = stack->Particle(TMath::Abs(labelP));
1059             
1060             Int_t labelMotherN = tparticleDaughterN->GetFirstMother();
1061             Int_t labelMotherP = tparticleDaughterP->GetFirstMother();
1062             
1063             TParticle *tparticleMotherN = stack->Particle(TMath::Abs(labelMotherN));
1064             TParticle *tparticleMotherP = stack->Particle(TMath::Abs(labelMotherP));
1065             
1066             //LambdaNeuteron
1067             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
1068               
1069               //cout << "number of daughters: " << tparticleMotherN->GetNDaughters() << endl;;
1070               if(tparticleMotherN->GetNDaughters() > 2.) continue;
1071               
1072               Int_t labelSecondDaughter = tparticleMotherP->GetDaughter(1);
1073               Int_t labelFirstDaughter = labelSecondDaughter-1;
1074           
1075               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1076               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1077               
1078               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGDeuteron]){//check first daughter PDG
1079                 
1080                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionMinus]){//check second daughter PDG
1081                   
1082                   if(invaMassDeuteronPion < 2.02) continue;
1083                   
1084                   fHistLambdaNeutronPtAso->Fill(transversMomentumMother);
1085                   fHistLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1086                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1087                   fHistLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1088                   //CUTS
1089                   
1090                   if((dcaV0 < 0.5) && 
1091                      (cosPointing > 0.999) && 
1092                      (decayRadius < 50.0) && 
1093                      (ptotP > 0.2)){//trackP->GetTPCsignal()>180 &&&&  decayRadius > 1.5 && decayRadius< 50 
1094                     
1095                     fHistLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1096                   }
1097                 }//end check second daughter PDG
1098               }//end check first daughter PDG
1099             }//end LambdaNeutron
1100             
1101             //Anti-LambdaNeutron
1102             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
1103               
1104               Int_t labelSecondDaughter = tparticleMotherN->GetDaughter(1);
1105               Int_t labelFirstDaughter = labelSecondDaughter-1;
1106               
1107               TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1108               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1109               
1110               if(tparticleFirstDaughter->GetPdgCode() == fgkPdgCode[kPDGAntiDeuteron]) {//check first daughter PDG
1111                 
1112                 if(tparticleSecondDaughter->GetPdgCode() == fgkPdgCode[kPDGPionPlus]) {//check second daughter PDG
1113                   
1114                   if(invaMassDeuteronPion < 2.02) continue;
1115                   
1116                   fIsCorrectlyAssociated[fItrk] = kTRUE;
1117                   fHistAntiLambdaNeutronPtAso->Fill(transversMomentumMother);
1118                   fHistAntiLambdaNeutronInvaMassAso->Fill(invaMassDeuteronPion);
1119                   fHistAntiLambdaNeutronDecayLengthAso->Fill((decayRadius*2.054)/(tparticleMotherP->P()));
1120                   
1121                   //CUTS
1122                   if(dcaV0 < 1. && cosPointing > 0.99 && ptotN > 0.2)//&& trackN->GetTPCsignal()>180 && decayRadius > 1.5 && decayRadius< 50
1123                     {
1124                       fHistAntiLambdaNeutronPtAsoCuts->Fill(transversMomentumMother);
1125                     }
1126                 }//end check second daughter PDG
1127               }//end check first daughter PDG
1128             }//end Anti-LambdaNeutron
1129           }//end MC
1130         }//end rotation == 1, signal
1131         
1132         if(rotation == 2){ // rotation of the pion
1133           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1134                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron]
1135                                                         + 2.*(TMath::Sqrt(r12*d22)
1136                                                               -momentumPionRot[0]*momentumDeuteron[0]
1137                                                               -momentumPionRot[1]*momentumDeuteron[1]
1138                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1139         }//end rotation == 2, rotation of the pion
1140         
1141         if(rotation == 3){// Rotation of the deuteron
1142           invaMassDeuteronPion = TMath::Sqrt(TMath::Max(fgkMass[kMassPion]*fgkMass[kMassPion]
1143                                                         +fgkMass[kMassDeuteron]*fgkMass[kMassDeuteron] 
1144                                                         + 2.*(TMath::Sqrt(e12*dr22)
1145                                                               -momentumPion[0]*momentumDeuteronRot[0]
1146                                                               -momentumPion[1]*momentumDeuteronRot[1]
1147                                                               -momentumPion[2]*momentumDeuteron[2]), 0.));
1148         }//end rotation == 3, rotation of the deuteron
1149         
1150         //fill the THnSparse and the tree variables
1151         
1152         //tree variables which are independent of the particle-species
1153         fV0finder[fItrk]         = onFlyStatus;
1154         fkMB[fItrk]              = fTriggerFired[0];
1155         fkCentral[fItrk]         = fTriggerFired[1];
1156         fkSemiCentral[fItrk]     = fTriggerFired[2];
1157         //fkEMCEJE[fItrk]          = fTriggerFired[3];
1158         //fkEMCEGA[fItrk]          = fTriggerFired[4];
1159         
1160         fPtotN[fItrk]            = trackN->P(); //InnerParam???
1161         fPtotP[fItrk]            = trackP->P(); //InnerParam???
1162         fMotherPt[fItrk]         = transversMomentumMother;
1163         
1164         fdEdxN[fItrk]            = trackN->GetTPCsignal();
1165         fdEdxP[fItrk]            = trackP->GetTPCsignal();
1166         fSignN[fItrk]            = trackN->Charge();
1167         fSignP[fItrk]            = trackP->Charge();
1168         
1169         fDCAv0[fItrk]            = dcaV0;
1170         fCosinePAv0[fItrk]       = cosPointing;
1171         fDecayRadiusTree[fItrk]  = decayRadius;
1172         
1173         //fAmenterosAlphaTree[fItrk] = alpha;
1174         //fAmenterosQtTree[fItrk] = qt;
1175         fRotationTree[fItrk] = rotation;
1176         
1177             
1178         if (isDeuteron[0] == kTRUE)  //pos deuteron
1179           {
1180             //fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1181             //fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1182             
1183             //fItrk++;
1184             
1185             if(invaMassDeuteronPion < 2.1) 
1186               {
1187                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotP*(trackP->Charge()), trackP->GetTPCsignal(),9);
1188                 if(chargeComboDeuteronPion == 2)fHistDeDxQA->Fill(ptotN*(trackN->Charge()), trackN->GetTPCsignal(),10);
1189               }
1190             //fHistArmenterosPodolanskiDeuteronPion->Fill(alpha*(-1),qt);
1191           }
1192
1193         if (isDeuteron[1] == kTRUE) 
1194           {
1195             fInvaMassDeuteronPionTree[fItrk] = invaMassDeuteronPion;
1196             fChargeComboDeuteronPionTree[fItrk] = chargeComboDeuteronPion;
1197             
1198             fItrk++;
1199             
1200             if(invaMassDeuteronPion < 2.1) 
1201               {
1202                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotN*trackN->Charge(), trackN->GetTPCsignal(),7);
1203                 if(chargeComboDeuteronPion == 0)fHistDeDxQA->Fill(ptotP*trackP->Charge(), trackP->GetTPCsignal(),8);
1204               }
1205             //fHistArmenterosPodolanskiAntiDeuteronPion->Fill(alpha,qt);
1206           }
1207         
1208       }//end rotation loop
1209       
1210  } //end of v0 loop
1211
1212   //fill the tree
1213   fTreeV0->Fill();
1214  
1215   // Post output data.
1216   PostData(1, fOutputContainer);
1217   PostData(2, fTreeV0);
1218 }      
1219
1220 //________________________________________________________________________
1221 void AliAnalysisTaskLambdaNAOD::Terminate(const Option_t *)
1222 {
1223   // Draw result to the screen
1224   // Called once at the end of the query
1225
1226   //get output data and darw 'fHistPt'
1227   if (!GetOutputData(0)) return;
1228   //TH1F *hist=(TH1F*)(((TObjArray*)GetOutputData(0))->FindObject("fHistPt"));
1229   //if (hist) hist->Draw();
1230 }
1231
1232 //_____________________________________________________________________________
1233 Int_t AliAnalysisTaskLambdaNAOD::Initialize() {
1234
1235  
1236   // -- Reset Event
1237   // ----------------
1238   ResetEvent();
1239
1240   return 0;
1241   }
1242 //________________________________________________________________________
1243 Int_t AliAnalysisTaskLambdaNAOD::SetupEvent() {
1244   // Setup Reading of event
1245
1246   ResetEvent();
1247
1248   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --   
1249   // -- Get Trigger 
1250   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --   
1251
1252   //Bool_t isTriggered = IsTriggered();
1253   IsTriggered();
1254   return 0;
1255 }
1256 //_____________________________________________________________________________
1257 void AliAnalysisTaskLambdaNAOD::ResetEvent() {
1258
1259   // Reset event
1260
1261   // -- Reset QA
1262   for (Int_t ii = 0; ii < fNTriggers; ++ii)
1263     fTriggerFired[ii] = kFALSE;
1264   
1265   return;
1266 }
1267 //_______________________________________________________________________
1268 void AliAnalysisTaskLambdaNAOD::BinLogAxis(const THnSparse *h, Int_t axisNumber) {
1269   //
1270   // Method for the correct logarithmic binning of histograms
1271   //
1272   TAxis *axis = h->GetAxis(axisNumber);
1273   int bins = axis->GetNbins();
1274
1275   Double_t from = axis->GetXmin();
1276   Double_t to = axis->GetXmax();
1277   Double_t *newBins = new Double_t[bins + 1];
1278    
1279   newBins[0] = from;
1280   Double_t factor = pow(to/from, 1./bins);
1281   
1282   for (int i = 1; i <= bins; i++) {
1283    newBins[i] = factor * newBins[i-1];
1284   }
1285   axis->Set(bins, newBins);
1286   delete [] newBins;
1287   
1288   }
1289 //________________________________________________________________________
1290 Bool_t AliAnalysisTaskLambdaNAOD::IsTriggered() {
1291   // Check if Event is triggered and fill Trigger Histogram
1292
1293   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
1294   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
1295   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
1296   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
1297   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
1298
1299   Bool_t isTriggered = kFALSE;
1300
1301   for (Int_t ii=0; ii<fNTriggers; ++ii) {
1302     if(fTriggerFired[ii]) {
1303       isTriggered = kTRUE;
1304       fHistTriggerStat->Fill(ii);
1305     }
1306   }
1307
1308   return isTriggered;
1309   }
1310 //______________________________________________________________________________
1311 Bool_t AliAnalysisTaskLambdaNAOD::DeuteronPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isDeuteron[3]) {
1312  
1313      //define the arrays for the Bethe-Bloch-Parameters
1314       Double_t parDeuteron[5] = {0,0,0,0,0};
1315
1316       if(runNumber < 166500) //LHC10h
1317         {
1318           parDeuteron[0] = 1.45802; // ALEPH parameters for deuterons (pass2)
1319           parDeuteron[1] = 27.4992;
1320           parDeuteron[2] = 4.00313e-15;
1321           parDeuteron[3] = 2.48485;
1322           parDeuteron[4] = 8.31768; 
1323         }
1324       
1325       if(runNumber > 166500) //LHC11h
1326         { 
1327           parDeuteron[0] = 1.11243; // ALEPH parameters for deuterons (pass2)
1328           parDeuteron[1] = 26.1144;
1329           parDeuteron[2] = 4.00313e-15;
1330           parDeuteron[3] = 2.72969 ;
1331           parDeuteron[4] = 9.15038; 
1332         }
1333  
1334
1335       //define expected signals for the various species
1336       Double_t expSignalDeuteronN = 0;
1337       Double_t expSignalDeuteronP = 0;
1338   
1339       isDeuteron[0] = kFALSE;
1340       isDeuteron[1] = kFALSE;
1341       isDeuteron[2] = kFALSE;
1342       
1343       //for data
1344       if(!fMCtrue){
1345        
1346         expSignalDeuteronN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1347         expSignalDeuteronP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1348         
1349         if(trackP->GetTPCsignal() >= 110 && ///??????????????????????????????????
1350            trackP->GetTPCsignal() < 1200 && 
1351            (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1352            ptotP > 0.2 ){
1353           
1354           if(trackP->Charge() >0)               isDeuteron[0] = kTRUE; //pos deuteron
1355           if(trackP->Charge() <0)               isDeuteron[1] = kTRUE; //neg deuteron
1356         }
1357         
1358         if(trackN->GetTPCsignal() >= 110 && ///??????????????????????????????????
1359            trackN->GetTPCsignal() < 1200 && 
1360            (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1361            ptotN > 0.2){ 
1362           
1363           isDeuteron[2] = kTRUE;
1364           
1365           if(trackN->Charge() >0)        isDeuteron[0] = kTRUE; //pos deuteron
1366           if(trackN->Charge() <0)        isDeuteron[1] = kTRUE; //neg deuteron
1367         }
1368       }
1369       
1370       //for MC
1371       if(fMCtrue)
1372         {
1373           if(runNumber < 166500) //2010
1374             { 
1375               expSignalDeuteronN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1376               expSignalDeuteronP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]); 
1377             }
1378           if(runNumber > 166500) //2011
1379             { 
1380               expSignalDeuteronN = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]);
1381               expSignalDeuteronP = 0.8*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassDeuteron]),parDeuteron[0],parDeuteron[1],parDeuteron[2],parDeuteron[3],parDeuteron[4]); 
1382             }
1383          
1384           if(trackP->GetTPCsignal() < 1200 && 
1385              (TMath::Abs(trackP->GetTPCsignal() - expSignalDeuteronP)/expSignalDeuteronP) < 0.2 &&
1386              ptotP > 0.2 )
1387             {
1388               if(trackP->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1389               if(trackP->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron         
1390             }
1391
1392           if(trackN->GetTPCsignal() < 1200 && 
1393              (TMath::Abs(trackN->GetTPCsignal() - expSignalDeuteronN)/expSignalDeuteronN) < 0.2 &&
1394              ptotN > 0.2 )
1395             {
1396               isDeuteron[2] = kTRUE;
1397
1398               if(trackN->Charge() >0)           isDeuteron[0] = kTRUE; //pos deuteron
1399               if(trackN->Charge() <0)           isDeuteron[1] = kTRUE; //neg deuteron         
1400             }
1401
1402         }
1403
1404
1405
1406      return isDeuteron;
1407   
1408 }
1409 //______________________________________________________________________________
1410 Bool_t AliAnalysisTaskLambdaNAOD::PionPID(AliVTrack *trackP, AliVTrack *trackN, Double_t ptotP, Double_t ptotN, Int_t runNumber, Bool_t isPion[2]) {
1411     
1412   //Pion PID via specific energy loss in the TPC
1413   //define the array for the Bethe-Bloch-Parameters (only necessary for local MC)
1414   Double_t parPion[5]  = {0,0,0,0,0};
1415   
1416     if(runNumber < 166500) { //LHC10h
1417         parPion[0]   = 1.45802; // ALEPH parameters for pions (pass2)
1418         parPion[1]   = 27.4992;
1419         parPion[2]   = 4.00313e-15;
1420         parPion[3]   = 2.48485;
1421         parPion[4]   = 8.31768; 
1422     }
1423       
1424     if(runNumber > 166500) {  //LHC11h
1425         parPion[0]   = 1.11243; // ALEPH parameters for pions (pass2)
1426         parPion[1]   = 26.1144;
1427         parPion[2]   = 4.00313e-15;
1428         parPion[3]   = 2.72969;
1429         parPion[4]   = 9.15038; 
1430     }
1431
1432
1433     Double_t expSignalPionP = 0;
1434     Double_t expSignalPionN = 0;
1435
1436     //for MC
1437     if(fMCtrue){
1438       if(runNumber < 166500){ //2010
1439         expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1440         expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1441       }
1442       if(runNumber > 166500){ //2011
1443         expSignalPionP = AliExternalTrackParam::BetheBlochAleph(ptotP/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1444         expSignalPionN = AliExternalTrackParam::BetheBlochAleph(ptotN/(fgkMass[kMassPion]),parPion[0],parPion[1],parPion[2],parPion[3],parPion[4]);
1445       }
1446       
1447     }
1448      
1449     isPion[0] = kFALSE;
1450     isPion[1] = kFALSE;
1451     //data
1452     if(!fMCtrue){
1453       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackP, AliPID::kPion))<3){
1454
1455         if(trackP->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1456         if(trackP->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1457      }
1458       if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackN, AliPID::kPion))<3){
1459
1460         if(trackN->Charge()>0 ) isPion[0] = kTRUE; //pos pion
1461         if(trackN->Charge()<0 ) isPion[1] = kTRUE; //neg pion
1462       }
1463     }
1464
1465     //MC
1466     if(fMCtrue){
1467       if(TMath::Abs(trackP->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.2
1468          && ptotP>0.00001){
1469         if(trackP->Charge()>0) isPion[0] = kTRUE; //pos pion
1470         if(trackP->Charge()<0) isPion[1] = kTRUE; //neg pion
1471       }
1472      if(TMath::Abs(trackN->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.2
1473         && ptotN>0.00001){
1474         if(trackN->Charge()>0) isPion[0] = kTRUE; //pos pion
1475         if(trackN->Charge()<0) isPion[1] = kTRUE; //neg pion
1476       }
1477     }
1478
1479     return isPion;
1480
1481 }
1482 //______________________________________________________________________________
1483 Bool_t AliAnalysisTaskLambdaNAOD::TrackCuts(AliVTrack *track, Bool_t testTrackCuts) {
1484
1485   testTrackCuts = kFALSE;
1486
1487   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1488   //if(!esdtrack) testTrackCuts = kFALSE;
1489   if (fESDtrackCutsV0->AcceptTrack(esdtrack)) testTrackCuts = kTRUE;
1490 //if( testTrackCuts == kTRUE) cout <<   "testTrackCuts im test: " << testTrackCuts << endl;
1491
1492
1493   return testTrackCuts;
1494 }
1495 //______________________________________________________________________________
1496 /*Bool_t AliAnalysisTaskLambdaNAOD::FilterBit(AliVTrack *track, Bool_t testFilterBit) {
1497
1498   testFilterBit = kFALSE;
1499
1500   AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1501   //if(!aodtrack) testFilterBit = kFALSE;
1502   if (aodtrack->TestFilterBit(7)) testFilterBit = kTRUE;
1503
1504   //if(testFilterBit == kTRUE) cout <<   "testFilterBit im test: " << testFilterBit<< endl;
1505
1506   return testFilterBit;
1507   }*/
1508 //______________________________________________________________________
1509 void AliAnalysisTaskLambdaNAOD::MCGenerated(AliStack* stack) 
1510
1511
1512   // Monte Carlo for genenerated particles
1513   if (fMCtrue) //MC loop  
1514     {
1515  
1516       Int_t stackN = 0;
1517
1518       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1519         {
1520
1521           const TParticle *tparticleMother = stack->Particle(stackN);
1522           Long_t pdgCodeMother = tparticleMother->GetPdgCode();
1523
1524           //check which particle the mother is 
1525
1526           //LambdaNeutron
1527           if(pdgCodeMother == fgkPdgCode[kPDGLambdaNeutron]) //check mother PDG 
1528             {
1529               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGDeuteron],fgkPdgCode[kPDGPionMinus],fgkMass[kMassDeuteron],fgkMass[kMassPion]);
1530             }
1531
1532           //Anti-LambdaNeutron
1533           if(pdgCodeMother == fgkPdgCode[kPDGAntiLambdaNeutron]) //check mother PDG 
1534             {
1535               MCTwoBodyDecay(stack,tparticleMother,pdgCodeMother,fgkPdgCode[kPDGAntiDeuteron],fgkPdgCode[kPDGPionPlus],fgkMass[kMassDeuteron],fgkMass[kMassPion]);
1536             }
1537
1538               
1539         }//end loop over stack
1540       
1541
1542     }//end MC
1543 }
1544 //_____________________________________________
1545 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
1546
1547   Double_t momentumFirstDaughterGen[3]={0,0,0};
1548   Double_t momentumSecondDaughterGen[3]={0,0,0};
1549   
1550   Double_t energyFirstDaughterGen = 0;
1551   Double_t energySecondDaughterGen = 0;
1552   
1553   Double_t transversMomentumMotherGen = 0;
1554   Double_t longitudinalMomentumMotherGen = 0;
1555   Double_t totalEnergyMotherGen = 0;
1556   
1557   Double_t rapidityGen = 2;
1558   
1559   //Int_t labelFirstDaughter = tparticleMother->GetDaughter(0);
1560   Int_t labelSecondDaughter = tparticleMother->GetDaughter(1);
1561   Int_t labelFirstDaughter = labelSecondDaughter -1;
1562
1563   TParticle *tparticleFirstDaughter = stack->Particle(TMath::Abs(labelFirstDaughter));
1564   TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1565
1566   if(tparticleFirstDaughter->GetPdgCode() == PDGFirstDaughter) //check first daughter PDG
1567     {
1568       if(tparticleSecondDaughter->GetPdgCode() == PDGSecondDaughter) //check second daughter PDG
1569         {
1570           
1571           momentumFirstDaughterGen[0] = tparticleFirstDaughter->Px();
1572           momentumFirstDaughterGen[1] = tparticleFirstDaughter->Py();
1573           momentumFirstDaughterGen[2] = tparticleFirstDaughter->Pz();
1574           
1575           momentumSecondDaughterGen[0] = tparticleSecondDaughter->Px();
1576           momentumSecondDaughterGen[1] = tparticleSecondDaughter->Py();
1577           momentumSecondDaughterGen[2] = tparticleSecondDaughter->Pz();
1578           
1579           energyFirstDaughterGen  = tparticleFirstDaughter->Energy();
1580           energySecondDaughterGen = tparticleSecondDaughter->Energy();
1581           
1582           transversMomentumMotherGen=TMath::Sqrt((momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])*(momentumFirstDaughterGen[0]+momentumSecondDaughterGen[0])+(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1])*(momentumFirstDaughterGen[1]+momentumSecondDaughterGen[1]));
1583           
1584           //rapidity cut +- 1
1585           //longitudinal momentum of mother
1586           longitudinalMomentumMotherGen = TMath::Sqrt((momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2])*(momentumFirstDaughterGen[2]+momentumSecondDaughterGen[2]));
1587           
1588           //total energy of mother
1589           totalEnergyMotherGen = energyFirstDaughterGen + energySecondDaughterGen;
1590           
1591           //rapidity
1592           rapidityGen = 0.5 * TMath::Log( (totalEnergyMotherGen+longitudinalMomentumMotherGen)/(totalEnergyMotherGen-longitudinalMomentumMotherGen));
1593           
1594           if(rapidityGen > 1.0 || rapidityGen < -1 ) return;
1595
1596           //cout << "decay lenght first daughter :" << tparticleFirstDaughter->Rho()  << " decay lenght second daughter :" << tparticleSecondDaughter->Rho()<< endl;
1597           //fHistLambdaNeutronDecayLengthGen->Fill(tparticleFirstDaughter->Rho());
1598
1599           //calculate the invariant mass
1600           Double_t firstDaughterPart = 0;
1601           Double_t secondDaughterPart = 0;
1602           Double_t invaMass = 0;
1603                       
1604           firstDaughterPart = massFirstDaughter*massFirstDaughter
1605             +momentumFirstDaughterGen[0]*momentumFirstDaughterGen[0]
1606             +momentumFirstDaughterGen[1]*momentumFirstDaughterGen[1]
1607             +momentumFirstDaughterGen[2]*momentumFirstDaughterGen[2];
1608           secondDaughterPart = massSecondDaughter*massSecondDaughter
1609             +momentumSecondDaughterGen[0]*momentumSecondDaughterGen[0]
1610             +momentumSecondDaughterGen[1]*momentumSecondDaughterGen[1]
1611             +momentumSecondDaughterGen[2]*momentumSecondDaughterGen[2];
1612           
1613           invaMass = TMath::Sqrt(TMath::Max(massFirstDaughter*massFirstDaughter
1614                                                         +massSecondDaughter*massSecondDaughter 
1615                                                         + 2.*(TMath::Sqrt(firstDaughterPart*secondDaughterPart)
1616                                                               -momentumFirstDaughterGen[0]*momentumSecondDaughterGen[0]
1617                                                               -momentumFirstDaughterGen[1]*momentumSecondDaughterGen[1]
1618                                                               -momentumFirstDaughterGen[2]*momentumSecondDaughterGen[2]), 0.));
1619
1620 #if 0
1621
1622                                                               switch(PDGMother) {
1623                                                               case fgkPdgCode[kPDGLambdaNeutron] : 
1624                                                                fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1625                                                                fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1626                                                                break;
1627                                                               case fgkPdgCode[kPDGAntiLambdaNeutron] :
1628                                                                fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1629                                                                fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1630                                                                break;
1631                                                               default :
1632                                                               printf("should not happen!!!! \n");
1633                                                               }
1634
1635
1636 #else
1637           //LambdaNeutron
1638           if(PDGMother == fgkPdgCode[kPDGLambdaNeutron])
1639             {  
1640               fHistLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1641               fHistLambdaNeutronInvaMassGen->Fill(invaMass);
1642               fHistLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1643             }
1644           //Anti-LambdaNeutron
1645           if(PDGMother == fgkPdgCode[kPDGAntiLambdaNeutron])
1646             {
1647               fHistAntiLambdaNeutronPtGen->Fill(transversMomentumMotherGen);
1648               fHistAntiLambdaNeutronInvaMassGen->Fill(invaMass);
1649               fHistAntiLambdaNeutronDecayLengthGen->Fill(((tparticleFirstDaughter->Rho())*2.054)/(tparticleMother->P()));
1650             }         
1651
1652 #endif    
1653         }//end of check second daughter PDG
1654     }//end of check first daughter PDG
1655
1656 }
1657 //_____________________________________________
1658 Double_t AliAnalysisTaskLambdaNAOD::MomentumInnerParam(AliVTrack *track, Double_t ptot){ //function to get the momentum at the inner wall of the TPC
1659  
1660   AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1661   ptot = esdtrack->GetInnerParam()->GetP();
1662
1663   return ptot;
1664