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