b3fb58d9c5d3193c8d58aaf6fbed230608b3fee6
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAnalysisTaskHdibaryonLPpi.cxx
1 /**************************************************************************
2  * Author : Benjamin Dönigus (benjamin.doenigus@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 //                 AliAnalysisTaskHdibaryonLPpi class
17 //          used to search for the H-Dibaryon in weak 
18 //          (Lambda Proton Pion) and strong (Lambda Lambda) decays
19 //-----------------------------------------------------------------
20
21 #include "Riostream.h"
22 #include "TROOT.h"
23 #include "TChain.h"
24 #include "TStyle.h"
25 #include "TSystem.h"
26 #include "TTree.h"
27 #include "TH1F.h"
28 #include "TH1.h"
29 #include "TH2D.h"
30 #include "TH3F.h"
31 #include "TCanvas.h"
32 #include "TParticle.h"
33 #include "TNtuple.h"
34 #include "TObjString.h"
35 #include "TLorentzVector.h"
36 #include "TDatabasePDG.h"
37
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40 #include "AliMCEventHandler.h"
41 #include "AliMCEvent.h"
42
43 #include "AliESD.h"
44 #include "AliESDpid.h"
45 #include "AliESDEvent.h"
46 #include "AliGenEventHeader.h"
47 #include "AliESDInputHandler.h"
48 #include "AliLog.h"
49 #include "AliStack.h"
50 #include "AliHeader.h"
51
52 #include "AliKFParticle.h"
53 #include "AliKFVertex.h"
54 #include "AliESDtrackCuts.h"
55 #include "AliESDv0Cuts.h"
56 #include "AliAnalysisTaskHdibaryonLPpi.h"
57
58 #include "TFile.h"
59 #include "TH2F.h"
60 #include "TF1.h"
61 #include "TCanvas.h"
62 #include "TList.h"
63 #include "TParallelCoord.h"
64
65 #include "AliMCParticle.h"
66 #include "AliGenPythiaEventHeader.h"
67
68 #include "AliPID.h"
69 #include "AliESDtrack.h"
70 #include "AliCentrality.h"
71
72 #include "THnSparse.h"
73
74 #include "AliVertexerTracks.h"
75
76 using namespace std;
77 ClassImp(AliAnalysisTaskHdibaryonLPpi)
78 //________________________________________________________________________
79 AliAnalysisTaskHdibaryonLPpi::AliAnalysisTaskHdibaryonLPpi() : AliAnalysisTaskSE()/*AliAnalysisTask(name, ""), fMCEvent(0)*/, fESD(0),   fESDtrackCutsV0(0),
80   fESDCutsV0(0),
81   fEsdTrackCuts(0),
82   fBin(0),
83   fEvent(0x0),
84   fHistList(0),  
85   fHistMassDPi(0), 
86   fHistMassLPi(0),
87   fHistMassLambdaPi(0),
88   fHistMassLambda(0),
89   fHistMassLambdaPPi(0),
90   fHistMassLambdaP(0),
91   fHistMassLambdaK(0),
92   fHistMassK0onl(0),
93   fHistMassK0offl(0),
94   fHistMassK0onlC(0),
95   fHistMassK0offlC(0),
96   fHistMassPQonl(0), 
97   fHistMassPQoffl(0),
98   fHistDC(0),
99   fHistArmenterosPodolanski(0),
100   fHistArmenterosPodolanskiCut(0), 
101   fHistHDibaryonInvaMassGen(0),
102   fHistHDibaryonInvaMassGenRes(0),
103   fHistAntiHDibaryonInvaMassGen(0),
104   fHistHDibaryonInvaMassAso(0),
105   fHistHDibaryonInvaMassAsoReso(0),
106   fHistAntiHDibaryonInvaMassAso(0),
107   fHistCheck(0),
108   fHistHPointingAngle(0),
109   fHistMassH(0),
110   fHistMassLambdaFromH(0),
111   fHistMassLambdaFromHtLorentz(0),
112   fHistMassPpi(0),
113   fHistMassPpiReso(0),
114   fHistMassLpi(0),
115   fHistMassLP(0),
116   fHistProtonPIDBb(0),
117   fHistPionPIDBb(0),
118   fHistProtonPIDLambda(0),
119   fHistPionPIDLambda(0),
120   fHistMCdcaPvtxDvtx(0),
121   fHistMCdcaPvtxLvtx(0),
122   fHistMCdcaDvtxLvtx(0),
123   fHistMCangleLH(0),
124   fHistMCdecayAngle(0),
125   fHistMCpointingAngle(0),
126   fHistMCap(0),
127   fHistMCdcaPvtxDvtxReso(0),
128   fHistMCdcaPvtxLvtxReso(0),
129   fHistMCdcaDvtxLvtxReso(0),
130   fHistMCangleLHReso(0),
131   fHistMCdecayAngleReso(0),
132   fHistMCpointingAngleReso(0),
133   fHistMCapReso(0),
134   fHistCentrality(0),
135   fHistCentralityAC(0), 
136   fHistMultiplicity(0),
137   fHistHilf1(0),
138   fHistHilf2(0), 
139   fHistHilf3(0),
140   fHistHilf4(0),
141   fHistHilf5(0), 
142   fHistHilf6(0),
143   fHistPtvsEtaGen(0),
144   fHistPtvsEtaAso(0),
145   fHistPtvsYGen(0),
146   fHistPtvsYAso(0), 
147   fHistRap(0),
148   fHistCount(0),
149   fPIDtpcESD(0),
150   fHistTriggerStat(0),
151   fHistTriggerStatAfterEventSelection(0), 
152   fHistMassHcentMult(0),
153   fHistNdim(0)
154 {
155   // DefaultConstructor
156
157 }
158
159 //________________________________________________________________________
160 AliAnalysisTaskHdibaryonLPpi::AliAnalysisTaskHdibaryonLPpi(const char *name) : AliAnalysisTaskSE(name)/*AliAnalysisTask(name, ""), fMCEvent(0)*/, fESD(0),   fESDtrackCutsV0(0),
161   fESDCutsV0(0),
162   fEsdTrackCuts(0),
163   fBin(0),
164   fEvent(0x0),
165   fHistList(0),  
166   fHistMassDPi(0), 
167   fHistMassLPi(0),
168   fHistMassLambdaPi(0),
169   fHistMassLambda(0),
170   fHistMassLambdaPPi(0),
171   fHistMassLambdaP(0),
172   fHistMassLambdaK(0),
173   fHistMassK0onl(0),
174   fHistMassK0offl(0),
175   fHistMassK0onlC(0),
176   fHistMassK0offlC(0),
177   fHistMassPQonl(0), 
178   fHistMassPQoffl(0),
179   fHistDC(0),
180   fHistArmenterosPodolanski(0),
181   fHistArmenterosPodolanskiCut(0), 
182   fHistHDibaryonInvaMassGen(0),
183   fHistHDibaryonInvaMassGenRes(0),
184   fHistAntiHDibaryonInvaMassGen(0),
185   fHistHDibaryonInvaMassAso(0),
186   fHistHDibaryonInvaMassAsoReso(0),
187   fHistAntiHDibaryonInvaMassAso(0),
188   fHistCheck(0),
189   fHistHPointingAngle(0),
190   fHistMassH(0),
191   fHistMassLambdaFromH(0),
192   fHistMassLambdaFromHtLorentz(0),
193   fHistMassPpi(0),
194   fHistMassPpiReso(0),
195   fHistMassLpi(0),
196   fHistMassLP(0),
197   fHistProtonPIDBb(0),
198   fHistPionPIDBb(0),
199   fHistProtonPIDLambda(0),
200   fHistPionPIDLambda(0),
201   fHistMCdcaPvtxDvtx(0),
202   fHistMCdcaPvtxLvtx(0),
203   fHistMCdcaDvtxLvtx(0),
204   fHistMCangleLH(0),
205   fHistMCdecayAngle(0),
206   fHistMCpointingAngle(0),
207   fHistMCap(0),
208   fHistMCdcaPvtxDvtxReso(0),
209   fHistMCdcaPvtxLvtxReso(0),
210   fHistMCdcaDvtxLvtxReso(0),
211   fHistMCangleLHReso(0),
212   fHistMCdecayAngleReso(0),
213   fHistMCpointingAngleReso(0),
214   fHistMCapReso(0),
215   fHistCentrality(0),
216   fHistCentralityAC(0), 
217   fHistMultiplicity(0), 
218   fHistHilf1(0),
219   fHistHilf2(0), 
220   fHistHilf3(0),
221   fHistHilf4(0),
222   fHistHilf5(0), 
223   fHistHilf6(0),
224   fHistPtvsEtaGen(0),
225   fHistPtvsEtaAso(0),
226   fHistPtvsYGen(0),
227   fHistPtvsYAso(0), 
228   fHistRap(0),
229   fHistCount(0),
230   fPIDtpcESD(0),
231   fHistTriggerStat(0),
232   fHistTriggerStatAfterEventSelection(0),
233   fHistMassHcentMult(0),
234   fHistNdim(0)
235
236 {
237   // Constructor
238
239   // Define input and output slots here
240   // Input from a TChain
241   DefineInput(0, TChain::Class());
242   // Output to TList container
243   DefineOutput(1, TList::Class()); //full
244
245   //MC info contol
246   if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())SetHasMC();
247
248   //V0 cuts
249
250   fESDtrackCutsV0 = new AliESDtrackCuts("AliESDtrackCutsV0","AliESDtrackCutsV0");
251   fESDtrackCutsV0->SetAcceptKinkDaughters(kFALSE);
252   fESDtrackCutsV0->SetMinNClustersTPC(80);
253   fESDtrackCutsV0->SetMaxChi2PerClusterTPC(5);
254   fESDtrackCutsV0->SetRequireTPCRefit(kTRUE);
255   fESDtrackCutsV0->SetEtaRange(-0.9,0.9);
256   fESDtrackCutsV0->SetPtRange(0.2,1.5);
257   fESDtrackCutsV0->SetMinDCAToVertexXY(2); //war inzwischen 1 & 3
258   fESDtrackCutsV0->SetMinDCAToVertexZ(2); //war inzwischen 1 & 3
259
260   fESDCutsV0 = new AliESDv0Cuts("AliESDCutsV0","AliESDCutsV0");
261   fESDCutsV0->SetMaxDcaV0Daughters(1.0);
262   fESDCutsV0->SetMinDcaNegToVertex(2); //1.5
263     fESDCutsV0->SetMinDcaPosToVertex(2); //1.5
264
265   //ESD Track cuts
266   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");    
267   fEsdTrackCuts->SetMinNClustersTPC(80);
268   fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
269   fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
270   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
271   fEsdTrackCuts->SetEtaRange(-0.9,0.9);
272 }
273
274 //____________________________________________________________
275 AliAnalysisTaskHdibaryonLPpi::~AliAnalysisTaskHdibaryonLPpi(){
276   //
277   // Destructor
278   //
279   if(fHistList){ 
280     fHistList->Clear();
281     delete fHistList;
282   }
283   if(fEsdTrackCuts) delete fEsdTrackCuts;
284   if(fESDtrackCutsV0) delete fESDtrackCutsV0;
285   if(fESDCutsV0) delete fESDCutsV0;
286
287 }
288
289 //________________________________________________________________________
290 void AliAnalysisTaskHdibaryonLPpi::UserCreateOutputObjects()
291 {
292   // Create histograms
293   // Called once
294
295  fHistList = new TList();
296  fHistList->SetOwner();
297
298  fHistMassDPi = new TH1F("fHistMassDPi", "Invariant mass distribution p+#pi^{-} ", 500, 1.0, 1.25);
299  fHistMassDPi->GetXaxis()->SetTitle("Invariant mass p+#pi^{-} (GeV/c^{2})");
300  fHistMassDPi->GetYaxis()->SetTitle("Entries");
301  fHistMassDPi->SetMarkerStyle(kFullCircle); 
302
303  fHistMassLPi = new TH1F("fHistMassLPi", "Offline Invariant mass distribution p+#pi^{-} ", 500, 1.0, 1.25);
304  fHistMassLPi->GetXaxis()->SetTitle("Invariant mass p+#pi^{-} (GeV/c^{2})");
305  fHistMassLPi->GetYaxis()->SetTitle("Entries");
306  fHistMassLPi->SetMarkerStyle(kFullCircle);
307
308  fHistMassLambdaPi = new TH1F("fHistMassLambdaPi", "Invariant mass distribution #Lambda+#pi^{-} ", 500, 1.2, 1.5);
309  fHistMassLambdaPi->GetXaxis()->SetTitle("Invariant mass #Lambda+#pi^{-} (GeV/c^{2})");
310  fHistMassLambdaPi->GetYaxis()->SetTitle("Entries");
311  fHistMassLambdaPi->SetMarkerStyle(kFullCircle);
312  
313  fHistMassLambda = new TH1F("fHistMassLambda", "Invariant mass distribution of #Lambda for further analyis", 500, 1.0, 1.2);
314  fHistMassLambda->GetXaxis()->SetTitle("Invariant mass p+#pi^{+} (GeV/c^{2})");
315  fHistMassLambda->GetYaxis()->SetTitle("Entries");
316  fHistMassLambda->SetMarkerStyle(kFullCircle);
317  
318  fHistMassLambdaPPi = new TH1F("fHistMassLambdaPPi", "Invariant mass distribution #Lambdap#pi^{-} ", 300, 2.2, 2.5);
319  fHistMassLambdaPPi->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})");
320  fHistMassLambdaPPi->GetYaxis()->SetTitle("Entries");
321  fHistMassLambdaPPi->SetMarkerStyle(kFullCircle);
322
323  fHistMassLambdaP = new TH1F("fHistMassLambdaP", "Invariant mass distribution #Lambdap ", 300, 2.2, 2.5);
324  fHistMassLambdaP->GetXaxis()->SetTitle("Invariant mass #Lambdap (GeV/c^{2})");
325  fHistMassLambdaP->GetYaxis()->SetTitle("Entries");
326  fHistMassLambdaP->SetMarkerStyle(kFullCircle);
327
328  fHistMassLambdaK = new TH1F("fHistMassLambdaK", "Invariant mass distribution #LambdaK^{-} ", 300, 1.6, 1.9);
329  fHistMassLambdaK->GetXaxis()->SetTitle("Invariant mass #LambdaK^{-} (GeV/c^{2})");
330  fHistMassLambdaK->GetYaxis()->SetTitle("Entries");
331  fHistMassLambdaK->SetMarkerStyle(kFullCircle);
332
333  fHistMassK0onl = new TH1F("fHistMassK0onl", "Invariant mass distribution K_{s}^{0} online V0 finder ", 400, 0.2, 1.0);
334  fHistMassK0onl->GetXaxis()->SetTitle("Invariant mass #pi^{+}#pi^{-} (GeV/c^{2})");
335  fHistMassK0onl->GetYaxis()->SetTitle("Entries");
336  fHistMassK0onl->SetMarkerStyle(kFullCircle);
337
338  fHistMassK0offl = new TH1F("fHistMassK0offl", "Invariant mass distribution  K_{s}^{0} offline V0 finder ", 400, 0.2, 1.0);
339  fHistMassK0offl->GetXaxis()->SetTitle("Invariant mass  #pi^{+}#pi^{-} (GeV/c^{2})");
340  fHistMassK0offl->GetYaxis()->SetTitle("Entries");
341  fHistMassK0offl->SetMarkerStyle(kFullCircle);
342  
343  fHistMassK0onlC = new TH1F("fHistMassK0onlC", "Invariant mass distribution K_{s}^{0} online V0 finder ", 400, 0.2, 1.0);
344  fHistMassK0onlC->GetXaxis()->SetTitle("Invariant mass #pi^{+}#pi^{-} (GeV/c^{2})");
345  fHistMassK0onlC->GetYaxis()->SetTitle("Entries");
346  fHistMassK0onlC->SetMarkerStyle(kFullCircle);
347
348  fHistMassK0offlC = new TH1F("fHistMassK0offlC", "Invariant mass distribution  K_{s}^{0} offline V0 finder ", 400, 0.2, 1.0);
349  fHistMassK0offlC->GetXaxis()->SetTitle("Invariant mass  #pi^{+}#pi^{-} (GeV/c^{2})");
350  fHistMassK0offlC->GetYaxis()->SetTitle("Entries");
351  fHistMassK0offlC->SetMarkerStyle(kFullCircle);
352
353  fHistMassPQonl = new TH1F("fHistMassPQonl", "Invariant mass distribution K_{s}^{0}p using online V0 finder ", 500, 1.3, 2.3);
354  fHistMassPQonl->GetXaxis()->SetTitle("Invariant mass K_{s}^{0}p (GeV/c^{2})");
355  fHistMassPQonl->GetYaxis()->SetTitle("Entries");
356  fHistMassPQonl->SetMarkerStyle(kFullCircle);
357
358  fHistMassPQoffl = new TH1F("fHistMassPQoffl", "Invariant mass distribution  K_{s}^{0}p using offline V0 finder ", 500, 1.3, 2.3);
359  fHistMassPQoffl->GetXaxis()->SetTitle("Invariant mass K_{s}^{0}p (GeV/c^{2})");
360  fHistMassPQoffl->GetYaxis()->SetTitle("Entries");
361  fHistMassPQoffl->SetMarkerStyle(kFullCircle);
362
363  fHistDC = new TH1F("fHistDC", "Proper decay length", 500, 0.0, 25);
364  fHistDC->GetXaxis()->SetTitle("c#tau (cm)");
365  fHistDC->GetYaxis()->SetTitle("Entries");
366  fHistDC->SetMarkerStyle(kFullCircle);
367   
368  fHistArmenterosPodolanski = new TH2F("fHistArmenterosPodolanski", "Armenteros-Podolanski plot", 200,-1.0,1.0, 500,0,1);
369  fHistArmenterosPodolanski->GetXaxis()->SetTitle("#alpha");
370  fHistArmenterosPodolanski->GetYaxis()->SetTitle("q_{t}");
371  fHistArmenterosPodolanski->SetMarkerStyle(kFullCircle);
372
373  fHistArmenterosPodolanskiCut = new TH2F("fHistArmenterosPodolanskiCut", "Armenteros-Podolanski plot after cut", 200,-1.0,1.0, 500,0,1);
374  fHistArmenterosPodolanskiCut->GetXaxis()->SetTitle("#alpha");
375  fHistArmenterosPodolanskiCut->GetYaxis()->SetTitle("q_{t}");
376  fHistArmenterosPodolanskiCut->SetMarkerStyle(kFullCircle);
377
378  fHistHDibaryonInvaMassGen = new TH1F("fHistHDibaryonInvaMassGen", "Generated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
379  fHistHDibaryonInvaMassGen->GetYaxis()->SetTitle("Counts");
380  fHistHDibaryonInvaMassGen->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
381
382  fHistHDibaryonInvaMassGenRes = new TH1F("fHistHDibaryonInvaMassGenRes", "Generated  #Lambda p #pi^{-} with particles in rapidity!!!", 200, 2.1, 2.3);
383  fHistHDibaryonInvaMassGenRes->GetYaxis()->SetTitle("Counts");
384  fHistHDibaryonInvaMassGenRes->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
385
386  fHistAntiHDibaryonInvaMassGen = new TH1F("fHistAntiHDibaryonInvaMassGen", "Generated  #bar{#Lambda} #bar{p} #pi^{+}", 200, 2.1, 2.3);
387  fHistAntiHDibaryonInvaMassGen->GetYaxis()->SetTitle("Counts");
388  fHistAntiHDibaryonInvaMassGen->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
389
390  fHistHDibaryonInvaMassAso = new TH1F("fHistHDibaryonInvaMassAso", "Associated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
391  fHistHDibaryonInvaMassAso->GetYaxis()->SetTitle("Counts");
392  fHistHDibaryonInvaMassAso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
393
394  fHistHDibaryonInvaMassAsoReso = new TH1F("fHistHDibaryonInvaMassAsoReso", "Associated  #Lambda p #pi^{-}", 200, 2.1, 2.3);
395  fHistHDibaryonInvaMassAsoReso->GetYaxis()->SetTitle("Counts");
396  fHistHDibaryonInvaMassAsoReso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
397
398  fHistAntiHDibaryonInvaMassAso = new TH1F("fHistAntiHDibaryonInvaMassAso", "Associated  #bar{#Lambda} #bar{p} #pi^{+}", 200, 2.1, 2.3);
399  fHistAntiHDibaryonInvaMassAso->GetYaxis()->SetTitle("Counts");
400  fHistAntiHDibaryonInvaMassAso->GetXaxis()->SetTitle("Invariant mass  (GeV/c^{2})");
401
402  fHistCheck = new TH2F("fHistCheck", "Check online/offline", 200, -0.5, 1.5, 200, -0.5, 1.5);
403  fHistCheck->GetXaxis()->SetTitle("offline");
404  fHistCheck->GetYaxis()->SetTitle("online");
405  fHistCheck->SetMarkerStyle(kFullCircle);
406
407  fHistHPointingAngle= new TH1F("fHistHPointingAngle", "Pointing angle distribution for #Lambdap#pi^{-}", 200, 0., 2*TMath::Pi());
408  fHistHPointingAngle->GetXaxis()->SetTitle("Pointing angle distribution for #Lambdap#pi^{-}");
409  fHistHPointingAngle->GetYaxis()->SetTitle("Entries");
410  fHistHPointingAngle->SetMarkerStyle(kFullCircle);
411
412  fHistMassH= new TH1F("fHistMassH", "Invariant mass distribution #Lambdap#pi^{-} (GeV/c^{2}) after pointing angle cut", 3000, 2.2, 2.5);
413  fHistMassH->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})");
414  fHistMassH->GetYaxis()->SetTitle("Entries");
415  fHistMassH->SetMarkerStyle(kFullCircle);
416  
417  fHistMassLambdaFromH= new TH1F("fHistMassLambdaFromH", "Invariant mass distribution #Lambda (GeV/c^{2}) asking for Mother to be H", 300, 1.0, 1.3);
418  fHistMassLambdaFromH->GetXaxis()->SetTitle("Invariant mass #Lambda (GeV/c^{2})");
419  fHistMassLambdaFromH->GetYaxis()->SetTitle("Entries");
420  fHistMassLambdaFromH->SetMarkerStyle(kFullCircle);
421
422  fHistMassLambdaFromHtLorentz= new TH1F("fHistMassLambdaFromHtLorentz", "Invariant mass distribution #Lambda (GeV/c^{2}) asking for Mother to be H", 300, 1.0, 1.3);
423  fHistMassLambdaFromHtLorentz->GetXaxis()->SetTitle("Invariant mass #Lambda (GeV/c^{2})");
424  fHistMassLambdaFromHtLorentz->GetYaxis()->SetTitle("Entries");
425  fHistMassLambdaFromHtLorentz->SetMarkerStyle(kFullCircle);
426  
427  fHistMassPpi= new TH1F("fHistMassPpi", "Invariant mass distribution of the p#pi^{-} used for combing with Lambda", 300, 1.0, 1.6);
428  fHistMassPpi->GetXaxis()->SetTitle("Invariant mass p#pi^{-} (GeV/c^{2})");
429  fHistMassPpi->GetYaxis()->SetTitle("Entries");
430  fHistMassPpi->SetMarkerStyle(kFullCircle);
431
432  fHistMassPpiReso= new TH1F("fHistMassPpiReso", "Invariant mass distribution of the p#pi^{-} used for combing with Lambda", 300, 1.0, 1.6);
433  fHistMassPpiReso->GetXaxis()->SetTitle("Invariant mass p#pi^{-} (GeV/c^{2})");
434  fHistMassPpiReso->GetYaxis()->SetTitle("Entries");
435  fHistMassPpiReso->SetMarkerStyle(kFullCircle);
436
437  fHistMassLpi= new TH1F("fHistMassLpi", "Invariant mass distribution of the #Lambda#pi^{-} used for combing with p", 300, 1.1, 1.7);
438  fHistMassLpi->GetXaxis()->SetTitle("Invariant mass #Lambda#pi^{-} (GeV/c^{2})");
439  fHistMassLpi->GetYaxis()->SetTitle("Entries");
440  fHistMassLpi->SetMarkerStyle(kFullCircle);
441
442  fHistMassLP= new TH1F("fHistMassLP", "Invariant mass distribution of the #Lambda p used for combing with #pi^{-}", 300, 2.0, 2.3);
443  fHistMassLP->GetXaxis()->SetTitle("Invariant mass #Lambda p (GeV/c^{2})");
444  fHistMassLP->GetYaxis()->SetTitle("Entries");
445  fHistMassLP->SetMarkerStyle(kFullCircle);
446
447  fHistProtonPIDBb = new TH2F("fHistProtonPIDBb", "dE/dx after p PID", 100, 0., 10, 100, 0, 100);
448  fHistProtonPIDBb->GetYaxis()->SetTitle("TPC Signal");
449  fHistProtonPIDBb->GetXaxis()->SetTitle("P (GeV/c)");
450  fHistProtonPIDBb->SetOption("scat");
451  fHistProtonPIDBb->SetMarkerStyle(kFullCircle);
452
453  fHistPionPIDBb = new TH2F("fHistPionPIDBb", "dE/dx after K PID", 100, 0., 10, 100, 0, 100);
454  fHistPionPIDBb->GetYaxis()->SetTitle("TPC Signal");
455  fHistPionPIDBb->GetXaxis()->SetTitle("P (GeV/c)");
456  fHistPionPIDBb->SetOption("scat");
457  fHistPionPIDBb->SetMarkerStyle(kFullCircle);
458
459  fHistProtonPIDLambda = new TH2F("fHistProtonPIDLambda", "dE/dx after p PID from V0", 100, 0., 10, 100, 0, 100);
460  fHistProtonPIDLambda->GetYaxis()->SetTitle("TPC Signal");
461  fHistProtonPIDLambda->GetXaxis()->SetTitle("P (GeV/c)");
462  fHistProtonPIDLambda->SetOption("scat");
463  fHistProtonPIDLambda->SetMarkerStyle(kFullCircle);
464
465  fHistPionPIDLambda = new TH2F("fHistPionPIDLambda", "dE/dx after #pi PID from V0", 100, 0, 10, 100, 0, 100);
466  fHistPionPIDLambda->GetYaxis()->SetTitle("TPC Signal");
467  fHistPionPIDLambda->GetXaxis()->SetTitle("P (GeV/c)");
468  fHistPionPIDLambda->SetOption("scat");
469  fHistPionPIDLambda->SetMarkerStyle(kFullCircle);
470
471  fHistMCdcaPvtxDvtx= new TH1F("fHistMCdcaPvtxDvtx", "MC True DCA Primary Vertex - H Decay Vertex", 300, -0.1, 11.9);
472  fHistMCdcaPvtxDvtx->GetXaxis()->SetTitle("dca prim. vtx- decay vtx (cm)");
473  fHistMCdcaPvtxDvtx->GetYaxis()->SetTitle("Entries");
474  fHistMCdcaPvtxDvtx->SetMarkerStyle(kFullCircle);
475
476  fHistMCdcaPvtxLvtx= new TH1F("fHistMCdcaPvtxLvtx", "MC True DCA Primary Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
477  fHistMCdcaPvtxLvtx->GetXaxis()->SetTitle("dca prim. vtx-#Lambda decay vtx (cm)");
478  fHistMCdcaPvtxLvtx->GetYaxis()->SetTitle("Entries");
479  fHistMCdcaPvtxLvtx->SetMarkerStyle(kFullCircle);
480
481  fHistMCdcaDvtxLvtx= new TH1F("fHistMCdcaDvtxLvtx", "MC True DCA H Decay Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
482  fHistMCdcaDvtxLvtx->GetXaxis()->SetTitle("dca H deacy vtx-#Lambda decay vtx (cm)");
483  fHistMCdcaDvtxLvtx->GetYaxis()->SetTitle("Entries");
484  fHistMCdcaDvtxLvtx->SetMarkerStyle(kFullCircle);
485
486  fHistMCangleLH= new TH1F("fHistMCangleLH", "MC True Angle between #Lambda and H", 300, 0., 2*TMath::Pi());
487  fHistMCangleLH->GetXaxis()->SetTitle("Angle (#Lambda-H)");
488  fHistMCangleLH->GetYaxis()->SetTitle("Entries");
489  fHistMCangleLH->SetMarkerStyle(kFullCircle);
490
491  fHistMCdecayAngle= new TH1F("fHistMCdecayAngle", "MC True Angle between decay products", 300, 0., 2*TMath::Pi());
492  fHistMCdecayAngle->GetXaxis()->SetTitle("Angle (#Lambda-p#pi)");
493  fHistMCdecayAngle->GetYaxis()->SetTitle("Entries");
494  fHistMCdecayAngle->SetMarkerStyle(kFullCircle);
495
496  fHistMCpointingAngle= new TH1F("fHistMCpointingAngle", "MC True Pointing Angle", 3000, 0., 2*TMath::Pi());
497  fHistMCpointingAngle->GetXaxis()->SetTitle("Pointing Angle");
498  fHistMCpointingAngle->GetYaxis()->SetTitle("Entries");
499  fHistMCpointingAngle->SetMarkerStyle(kFullCircle);
500
501  fHistMCap = new TH2F("fHistMCap", "True MC Armenteros-Podolanski", 200,-1.0,1.0, 500,0,1);
502  fHistMCap->GetYaxis()->SetTitle("#alpha");
503  fHistMCap->GetXaxis()->SetTitle("q_{t}");
504  fHistMCap->SetOption("scat");
505  fHistMCap->SetMarkerStyle(kFullCircle);
506
507  fHistMCdcaPvtxDvtxReso= new TH1F("fHistMCdcaPvtxDvtxReso", "MC True DCA Primary Vertex - H Decay Vertex", 300, -0.1, 11.9);
508  fHistMCdcaPvtxDvtxReso->GetXaxis()->SetTitle("dca prim. vtx- decay vtx (cm)");
509  fHistMCdcaPvtxDvtxReso->GetYaxis()->SetTitle("Entries");
510  fHistMCdcaPvtxDvtxReso->SetMarkerStyle(kFullCircle);
511
512  fHistMCdcaPvtxLvtxReso= new TH1F("fHistMCdcaPvtxLvtxReso", "MC True DCA Primary Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
513  fHistMCdcaPvtxLvtxReso->GetXaxis()->SetTitle("dca prim. vtx-#Lambda decay vtx (cm)");
514  fHistMCdcaPvtxLvtxReso->GetYaxis()->SetTitle("Entries");
515  fHistMCdcaPvtxLvtxReso->SetMarkerStyle(kFullCircle);
516
517  fHistMCdcaDvtxLvtxReso= new TH1F("fHistMCdcaDvtxLvtxReso", "MC True DCA H Decay Vertex - Lambda Decay Vertex", 300, -0.1, 11.9);
518  fHistMCdcaDvtxLvtxReso->GetXaxis()->SetTitle("dca H deacy vtx-#Lambda decay vtx (cm)");
519  fHistMCdcaDvtxLvtxReso->GetYaxis()->SetTitle("Entries");
520  fHistMCdcaDvtxLvtxReso->SetMarkerStyle(kFullCircle);
521
522  fHistMCangleLHReso= new TH1F("fHistMCangleLHReso", "MC True Angle between #Lambda and H", 300, 0., 2*TMath::Pi());
523  fHistMCangleLHReso->GetXaxis()->SetTitle("Angle (#Lambda-H)");
524  fHistMCangleLHReso->GetYaxis()->SetTitle("Entries");
525  fHistMCangleLHReso->SetMarkerStyle(kFullCircle);
526
527  fHistMCdecayAngleReso= new TH1F("fHistMCdecayAngleReso", "MC True Angle between decay products", 300, 0., 2*TMath::Pi());
528  fHistMCdecayAngleReso->GetXaxis()->SetTitle("Angle (#Lambda-p#pi)");
529  fHistMCdecayAngleReso->GetYaxis()->SetTitle("Entries");
530  fHistMCdecayAngleReso->SetMarkerStyle(kFullCircle);
531
532  fHistMCpointingAngleReso= new TH1F("fHistMCpointingAngleReso", "MC True Pointing Angle", 300, 0., 2*TMath::Pi());
533  fHistMCpointingAngleReso->GetXaxis()->SetTitle("Pointing Angle");
534  fHistMCpointingAngleReso->GetYaxis()->SetTitle("Entries");
535  fHistMCpointingAngleReso->SetMarkerStyle(kFullCircle);
536
537  fHistMCapReso = new TH2F("fHistMCapReso", "True MC Armenteros-Podolanski", 200,-1.0,1.0, 500,0,1);
538  fHistMCapReso->GetYaxis()->SetTitle("#alpha");
539  fHistMCapReso->GetXaxis()->SetTitle("q_{t}");
540  fHistMCapReso->SetOption("scat");
541  fHistMCapReso->SetMarkerStyle(kFullCircle);
542
543  fHistCentrality = new TH1F("Centrality ", "Centrality", 11, -0.5, 10.5);
544  fHistCentrality ->GetXaxis()->SetTitle("Centrality");
545  fHistCentrality ->GetYaxis()->SetTitle("Entries");
546
547  fHistCentralityAC = new TH1F("CentralityAC ", "CentralityAC", 11, -0.5, 10.5);
548  fHistCentralityAC ->GetXaxis()->SetTitle("Centrality");
549  fHistCentralityAC ->GetYaxis()->SetTitle("Entries");
550
551  fHistMultiplicity = new TH1F("Multiplicity ", "Multiplicity", 100, 0, 20000);
552  fHistMultiplicity ->GetXaxis()->SetTitle("Centrality");
553  fHistMultiplicity ->GetYaxis()->SetTitle("Entries");
554
555  fHistList->Add(fHistMassDPi);
556  fHistList->Add(fHistMassLPi);
557  fHistList->Add(fHistMassLambdaPi);
558  fHistList->Add(fHistMassLambda);
559  fHistList->Add(fHistMassLambdaPPi);
560  fHistList->Add(fHistMassLambdaP);
561  fHistList->Add(fHistMassLambdaK);
562  fHistList->Add(fHistMassK0onl);
563  fHistList->Add(fHistMassK0offl);
564  fHistList->Add(fHistMassK0onlC);
565  fHistList->Add(fHistMassK0offlC);
566  fHistList->Add(fHistMassPQonl); 
567  fHistList->Add(fHistMassPQoffl);
568  fHistList->Add(fHistDC); 
569  fHistList->Add(fHistArmenterosPodolanski);
570  fHistList->Add(fHistArmenterosPodolanskiCut);
571  fHistList->Add(fHistHDibaryonInvaMassGen);
572  fHistList->Add(fHistHDibaryonInvaMassGenRes);
573  fHistList->Add(fHistAntiHDibaryonInvaMassGen);
574  fHistList->Add(fHistHDibaryonInvaMassAso);
575  fHistList->Add(fHistHDibaryonInvaMassAsoReso);
576  fHistList->Add(fHistAntiHDibaryonInvaMassAso);
577  fHistList->Add(fHistCheck);
578  fHistList->Add(fHistHPointingAngle);
579  fHistList->Add(fHistMassH);
580  fHistList->Add(fHistMassPpi);
581  fHistList->Add(fHistMassPpiReso);
582  fHistList->Add(fHistMassLpi);
583  fHistList->Add(fHistMassLP);
584  fHistList->Add(fHistMassLambdaFromH);
585  fHistList->Add(fHistMassLambdaFromHtLorentz);
586  fHistList->Add(fHistProtonPIDBb);
587  fHistList->Add(fHistPionPIDBb);
588  fHistList->Add(fHistProtonPIDLambda);
589  fHistList->Add(fHistPionPIDLambda);
590  fHistList->Add(fHistMCdcaPvtxDvtx);
591  fHistList->Add(fHistMCdcaPvtxLvtx);
592  fHistList->Add(fHistMCdcaDvtxLvtx);
593  fHistList->Add(fHistMCangleLH);
594  fHistList->Add(fHistMCdecayAngle);
595  fHistList->Add(fHistMCpointingAngle);
596  fHistList->Add(fHistMCap);
597  fHistList->Add(fHistMCdcaPvtxDvtxReso);
598  fHistList->Add(fHistMCdcaPvtxLvtxReso);
599  fHistList->Add(fHistMCdcaDvtxLvtxReso);
600  fHistList->Add(fHistMCangleLHReso);
601  fHistList->Add(fHistMCdecayAngleReso);
602  fHistList->Add(fHistMCpointingAngleReso);
603  fHistList->Add(fHistMCapReso);
604  fHistList->Add(fHistCentrality);
605  fHistList->Add(fHistCentralityAC);
606  fHistList->Add(fHistMultiplicity);
607  
608  fHistHilf1= new TH1F("fHistHilf1", "HD", 300, 0., 10);
609  fHistHilf1->GetXaxis()->SetTitle("");
610  fHistHilf1->GetYaxis()->SetTitle("Entries");
611
612  fHistHilf2= new TH1F("fHistHilf2", "HD", 300, 0., 10);
613  fHistHilf2->GetXaxis()->SetTitle("");
614  fHistHilf2->GetYaxis()->SetTitle("Entries");
615
616  fHistHilf3= new TH1F("fHistHilf3", "HD", 300, 0., 10);
617  fHistHilf3->GetXaxis()->SetTitle("");
618  fHistHilf3->GetYaxis()->SetTitle("Entries");
619
620  fHistHilf4= new TH1F("fHistHilf4", "HD", 300, 0., 10);
621  fHistHilf4->GetXaxis()->SetTitle("");
622  fHistHilf4->GetYaxis()->SetTitle("Entries");
623
624  fHistHilf5= new TH1F("fHistHilf5", "HD", 300, 0., 10);
625  fHistHilf5->GetXaxis()->SetTitle("");
626  fHistHilf5->GetYaxis()->SetTitle("Entries");
627
628  fHistHilf6= new TH1F("fHistHilf6", "HD", 300, 0., 10);
629  fHistHilf6->GetXaxis()->SetTitle("");
630  fHistHilf6->GetYaxis()->SetTitle("Entries");
631  
632  fHistPtvsEtaGen = new TH2F("fHistPtvsEtaGen", "p_{t} vs #eta from generated H", 200,0.0,10.0, 200,-1,1);
633  fHistPtvsEtaGen->GetXaxis()->SetTitle("p_{t} (GeV/c)");
634  fHistPtvsEtaGen->GetYaxis()->SetTitle("#eta");
635  fHistPtvsEtaGen->SetOption("scat");
636  fHistPtvsEtaGen->SetMarkerStyle(kFullCircle);
637
638  fHistPtvsEtaAso = new TH2F("fHistPtvsEtaAso", "p_{t} vs #eta from associated H", 200,0.0,10.0, 200,-1,1);
639  fHistPtvsEtaAso->GetYaxis()->SetTitle("p_{t} (GeV/c)");
640  fHistPtvsEtaAso->GetXaxis()->SetTitle("#eta");
641  fHistPtvsEtaAso->SetOption("scat");
642  fHistPtvsEtaAso->SetMarkerStyle(kFullCircle);
643
644  fHistPtvsYGen = new TH2F("fHistPtvsYGen", "p_{t} vs rapidity from generated H", 200,0.0,10.0, 200,-1,1);
645  fHistPtvsYGen->GetXaxis()->SetTitle("p_{t} (GeV/c)");
646  fHistPtvsYGen->GetYaxis()->SetTitle("y");
647  fHistPtvsYGen->SetOption("scat");
648  fHistPtvsYGen->SetMarkerStyle(kFullCircle);
649
650  fHistPtvsYAso = new TH2F("fHistPtvsYAso", "p_{t} vs rapidity from associated H", 200,0.0,10.0, 200,-1,1);
651  fHistPtvsYAso->GetXaxis()->SetTitle("p_{t} (GeV/c)");
652  fHistPtvsYAso->GetYaxis()->SetTitle("y");
653  fHistPtvsYAso->SetOption("scat");
654  fHistPtvsYAso->SetMarkerStyle(kFullCircle);
655
656  fHistRap= new TH1F("fHistRap", "Rapidity", 400, -2., 2);
657  fHistRap->GetXaxis()->SetTitle("Y");
658  fHistRap->GetYaxis()->SetTitle("Entries");
659  fHistPtvsEtaAso->SetMarkerStyle(kFullCircle);
660
661  fHistList->Add(fHistHilf1);
662  fHistList->Add(fHistHilf2); 
663  fHistList->Add(fHistHilf3);
664  fHistList->Add(fHistHilf4);
665  fHistList->Add(fHistHilf5); 
666  fHistList->Add(fHistHilf6);
667  fHistList->Add(fHistPtvsEtaGen);
668  fHistList->Add(fHistPtvsEtaAso);
669  fHistList->Add(fHistPtvsYGen);
670  fHistList->Add(fHistPtvsYAso);
671  fHistList->Add(fHistRap);
672
673  fHistCount = new TH1F("fHistCount","test",17,0,17);
674  fHistCount->GetXaxis()->SetBinLabel(1,"Events");
675  fHistCount->GetXaxis()->SetBinLabel(2,"MC All");
676  fHistCount->GetXaxis()->SetBinLabel(3,"MC from Primary Vtx");
677  fHistCount->GetXaxis()->SetBinLabel(4,"Horst");
678  fHistCount->GetXaxis()->SetBinLabel(5,"Lambda Candidates");
679  fHistCount->GetXaxis()->SetBinLabel(6,"Sigma Candidates");
680  fHistCount->GetXaxis()->SetBinLabel(7,"Horst");
681  fHistCount->GetXaxis()->SetBinLabel(8,"Horst");
682  fHistCount->GetXaxis()->SetBinLabel(9,"Horst");
683  fHistCount->GetXaxis()->SetBinLabel(10,"MC All #bar{Lambda}(1520)s");
684  fHistCount->GetXaxis()->SetBinLabel(11,"");
685  fHistCount->GetXaxis()->SetBinLabel(12,"H-Dibaryon");
686  fHistCount->GetXaxis()->SetBinLabel(13,"Hypertriton 2-Body");
687  fHistCount->GetXaxis()->SetBinLabel(14,"Hypertriton 3-Body");
688  fHistCount->GetXaxis()->SetBinLabel(15,"");
689  fHistCount->GetXaxis()->SetBinLabel(16,"");
690  fHistCount->GetXaxis()->SetBinLabel(17,"Lambdas!!!");
691  fHistCount->SetStats(0);
692  fHistCount->SetFillColor(38);
693  fHistList->Add(fHistCount);
694
695  //trigger statistics histogram
696   fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", 4,-0.5, 3.5);
697   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral" };
698   for ( Int_t ii=0; ii < 3; ii++ )
699     fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
700
701   fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", 4,-0.5, 3.5);
702   for ( Int_t ii=0; ii < 3; ii++ )
703     fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
704   fHistList->Add(fHistTriggerStat);
705   fHistList->Add(fHistTriggerStatAfterEventSelection);
706
707   fHistMassHcentMult =  new TH3F("fHistMassHcentMult", "Inv. Mass vs. centrality vs. multiplicity", 100, 2.2, 2.3, 5, 0, 4, 300, 0, 6000);
708   fHistMassHcentMult->GetXaxis()->SetTitle("Invariant mass #Lambdap#pi^{-} (GeV/c^{2})"); // inv. mass
709   fHistMassHcentMult->GetYaxis()->SetTitle("Centrality"); // triggertype
710   fHistMassHcentMult->GetZaxis()->SetTitle("Multiplicity"); // refTPC
711   fHistList->Add(fHistMassHcentMult);
712
713   
714   const Double_t kz = 2*TMath::Pi();
715   Int_t binsD01[16]={  300, 200, 100, 100, 100, 100, 100, 100, 200, 200, 200, 200, 400, 200, 200, 3};
716   Double_t xminD01[16]={2.0, 1.0, 0., -1, 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.,  0., 0};
717   Double_t xmaxD01[16]={2.3, 1.2, kz, 1, 1, 10, 10, 5, 5, 5, 5, 100, 1, 100,  4000, 1};
718   
719
720   fHistNdim = new THnSparseF("fHistNdim","THnS;InvMass, InvMassLambda, pointingAngle, armPoAlpha, armPoQt, pTL, pTH, d0p, d0n, dcaHd, dca, decayL, cosPA, centr, multi, mcinf;InvMassH", 16,binsD01,xminD01,xmaxD01);
721   fHistList->Add(fHistNdim);
722
723  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
724  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
725  fPIDtpcESD = inputHandler->GetPIDResponse();
726
727 // Post output data (if histograms are not used later, PostData is at least called here)
728   PostData(1, fHistList);
729
730 }
731
732  //________________________________________________________________________
733 void AliAnalysisTaskHdibaryonLPpi::UserExec(Option_t *) 
734
735   // Main loop
736   // Called for each event
737
738   //define improtant masses
739   Double_t pionMass     = 0.13957;
740   Double_t protonMass   = 0.93827;
741
742   //define PDGCodes
743   Long_t pdgPionPlus          =         211;
744   Long_t pdgPionMinus         =        -211;
745   Long_t pdgProton            =        2212;
746   Long_t pdgAntiProton        =       -2212;
747   Long_t pdgLambda            =        3122;
748   Long_t pdgAntiLambda        =       -3122;
749   Long_t pdgHDibaryon         =  1020000020;
750   Long_t pdgAntiHDibaryon     = -1020000020;
751
752   AliStack* stack(NULL);
753   if(HasMC()){
754
755     if(!fMCEvent)return;
756
757     AliHeader *head = fMCEvent->Header();
758     if(!head)return;
759     AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)head->GenEventHeader();
760     if(!pyheader)return;
761     
762     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
763       if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()){
764         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->TreeK()))return;
765         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->TreeTR()))return;
766         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->Init("local")))return;
767         stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
768         if(!(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack()->TreeK()))return;
769       }
770     }
771         
772     if(!stack)return;
773   }
774
775   // -------------------------------------------------------
776   // Loop for Inv. Mass via ESD tracks
777   // -------------------------------------------------------
778
779   fHistCount->Fill(0);
780   
781   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
782
783     if (!fESD) {
784       //Printf("ERROR: fESD not available");
785       return;
786     }
787
788
789   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
790   if (vertex->GetNContributors()<1) 
791     {
792       // SPD vertex
793       vertex = fESD->GetPrimaryVertexSPD();
794        if(vertex->GetNContributors()<1) {
795          return;
796        }  
797     }
798   if (TMath::Abs(vertex->GetZv()) > 10) return;
799   
800   Int_t centrality = -5;
801   Double_t centrPerc = -5;
802
803   if (fESD->GetEventSpecie() == 4) 
804     { // PbPb
805       AliCentrality *esdCentrality = fESD->GetCentrality();
806       centrality = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
807       centrPerc = esdCentrality->GetCentralityPercentile("V0M");
808 //      if (centrality < 0. || centrality > 8.) return; //0 bis 80 %
809       if (centrality > 8) return; //0 bis 80 %
810       //  cout<<"Centrality: "<< centrality << endl;
811     }
812
813
814   fHistCentrality->Fill(centrality);   
815
816   //*****************//  
817   //*   Centrality  *//
818   //*****************//
819  
820   //  Float_t percentile=centrality->GetCentralityPercentile("V0M");
821
822   Bool_t isSelectedCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
823   Bool_t isSelectedSemiCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
824   Bool_t isSelectedMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
825  
826   Int_t triggertype = 17;
827
828   if(isSelectedCentral){
829     fHistTriggerStat->Fill(1);
830     triggertype=1;
831   }
832
833   if(isSelectedSemiCentral){
834     fHistTriggerStat->Fill(2);
835     triggertype=2;
836   }
837
838   if(isSelectedMB){
839     fHistTriggerStat->Fill(0);
840     triggertype=3;
841   }
842
843   //  if(isSelectedCentral || isSelectedSemiCentral || isSelectedMB){
844
845   //*******************************
846
847   Int_t runNumber = 0;
848   //  itrk = 0;
849   runNumber = fESD->GetRunNumber();
850 /*  
851   if (!fPIDtpcESD) fPIDtpcESD = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
852   if (!fPIDtpcESD) {
853     fPIDtpcESD = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
854     fPIDtpcESD->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
855   }
856 */
857   Double_t pionK=1;
858   Double_t pK=1;
859  
860   TObjArray* listCrossV0   = fESDCutsV0->GetAcceptedV0s(fESD);  
861   Int_t nGoodV0s      = listCrossV0->GetEntries();
862
863   const AliESDVertex *esdVer = fESD->GetPrimaryVertex();
864   AliESDVertex *esdVer1 = new AliESDVertex(*esdVer);
865    
866   AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
867   TObjArray *trkArray = new TObjArray(2);
868   AliVertexerTracks *vertexer1 = new AliVertexerTracks(fESD->GetMagneticField());
869   TObjArray *trkArray1 = new TObjArray(2);
870
871   AliKFParticle::SetField(fESD->GetMagneticField());
872         
873   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
874
875   Int_t refMultTpc = AliESDtrackCuts::GetReferenceMultiplicity(fESD, kTRUE);
876   //cout<<"Multiplicity: "<< refMultTpc << endl;
877   fHistMultiplicity->Fill(refMultTpc);
878   
879   Double_t mn[3] = {0,0,0};
880   Double_t mp[3] = {0,0,0}; 
881   Double_t mm[3] = {0,0,0};
882   Double_t dd[3] = {0,0,0};
883   Double_t dd1[3] = {0,0,0};
884   const Double_t cProtonMass=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
885   const Double_t cPionMass=TDatabasePDG::Instance()->GetParticle(211)->Mass();
886   const Double_t cElectronMass=TDatabasePDG::Instance()->GetParticle(11)->Mass();
887   const Double_t cLambdaMass=TDatabasePDG::Instance()->GetParticle(3122)->Mass();
888   Double_t decayLength=0;
889   Double_t decayLengthH=0;
890
891   //V0 Loop for Lambda and Anti-Lambda
892   for(Int_t iV0MI = 0; iV0MI < nGoodV0s ; iV0MI++) {
893     AliESDv0 * fV0MIs = fESD->GetV0(iV0MI);
894     Int_t    lOnFlyStatus = 0;
895
896     lOnFlyStatus = fV0MIs->GetOnFlyStatus();
897     Double_t lInvMassLambda=0;
898     Double_t lInvMassLambdaPi=0;
899     Double_t lPtLambda=0;
900     Double_t lPzLambda=0;
901     Double_t lPLambda=0;
902     Int_t onl=0;
903     Int_t offl=0;
904
905     TLorentzVector posE;
906     TLorentzVector negE;
907     TLorentzVector photon;
908
909     if (lOnFlyStatus){
910       onl=1; 
911       //      return;
912     }
913     if (!lOnFlyStatus){
914       offl=1;
915       //return;
916     }
917
918     //    fHistMultiplicity->Fill(refMultTpc); 
919     fHistCentralityAC->Fill(centrality);
920     fHistCheck->Fill(offl,onl);
921
922     AliESDtrack* trackPosTest = fESD->GetTrack(fV0MIs->GetPindex());
923     AliESDtrack* trackNegTest = fESD->GetTrack(fV0MIs->GetNindex());
924
925     //    if (!
926     if (!fEsdTrackCuts->AcceptTrack(trackPosTest)) continue;
927     if (!fESDtrackCutsV0->AcceptTrack(trackNegTest)) continue;
928
929
930       //PID via specific energy loss in the TPC
931       //define the arrays for the Bethe-Bloch-Parameters
932       Double_t parProton[5]   = {0,0,0,0,0};
933
934       if(runNumber < 166500) //LHC10h
935         {
936           parProton[0]   = 1.45802; // ALEPH parameters for protons (pass2)
937           parProton[1]   = 27.4992;
938           parProton[2]   = 4.00313e-15;
939           parProton[3]   = 2.48485;
940           parProton[4]   = 8.31768; 
941         }
942       
943       if(runNumber > 166500) //LHC11h
944         { 
945           parProton[0]   = 1.11243; // ALEPH parameters for protons (pass2)
946           parProton[1]   = 26.1144;
947           parProton[2]   = 4.00313e-15;
948           parProton[3]   = 2.72969;
949           parProton[4]   = 9.15038; 
950         }
951  
952       //Get the total momentum for each track at the inner readout of the TPC
953       Double_t ptotN = trackNegTest->GetInnerParam()->GetP();
954       Double_t ptotP = trackPosTest->GetInnerParam()->GetP();
955
956       //define expected signals for the various species
957       Double_t expSignalPionP = 0;
958       Double_t expSignalPionN = 0;
959       Double_t expSignalProtonN = 0;
960       Double_t expSignalProtonP = 0;
961
962       //for data
963       if(!HasMC())
964         {
965           expSignalProtonN = AliExternalTrackParam::BetheBlochAleph(ptotN/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
966           expSignalProtonP = AliExternalTrackParam::BetheBlochAleph(ptotP/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
967         }
968     
969       //for MC
970       if(HasMC())
971         {
972           expSignalPionP = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotP/(pionMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]); 
973           expSignalPionN = 0.7*AliExternalTrackParam::BetheBlochAleph(ptotN/(pionMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]); 
974
975           expSignalProtonN = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotN/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);
976           expSignalProtonP = 0.65*AliExternalTrackParam::BetheBlochAleph(ptotP/(protonMass),parProton[0],parProton[1],parProton[2],parProton[3],parProton[4]);   
977         }
978
979       // PID cut on the nuclei (proton, deuteron, triton, helium3)
980       Bool_t corrParticle = kFALSE;
981
982       Bool_t posProton = kFALSE;
983       //data
984       if(!HasMC())
985         {
986           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackPosTest, AliPID::kProton)) < 3)
987             {
988               posProton = kTRUE;
989               corrParticle = kTRUE;
990             }
991         }
992       //MC
993       if(HasMC())
994         {
995           if(//trackPosTest->GetTPCsignal() < 1200 && 
996              TMath::Abs(trackPosTest->GetTPCsignal() - expSignalProtonP)/expSignalProtonP < 0.4)
997             {
998               posProton = kTRUE;
999               corrParticle = kTRUE;
1000             }
1001         }
1002
1003       Bool_t negProton = kFALSE;
1004       //data
1005       if(!HasMC())
1006         {
1007           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackNegTest, AliPID::kProton)) < 3)
1008             {
1009               negProton = kTRUE;
1010               corrParticle = kTRUE;
1011             }
1012         }
1013      //MC
1014       if(HasMC())
1015         {
1016           if(//trackNegTest->GetTPCsignal() < 1200 && 
1017              TMath::Abs(trackNegTest->GetTPCsignal() - expSignalProtonN)/expSignalProtonN < 0.4)
1018             {
1019               negProton = kTRUE;
1020               corrParticle = kTRUE;
1021             }
1022         }
1023
1024      //PID cut for pions
1025       //data: 3sigma cut on the pions
1026
1027       Bool_t negPion = kFALSE;     
1028       Bool_t posPion = kFALSE;
1029
1030       if (!HasMC())
1031         {
1032           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackPosTest, AliPID::kPion)) < 4) posPion=kTRUE; //pos daughter has to be a pion
1033           if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackNegTest, AliPID::kPion)) < 4) negPion=kTRUE; // negative daughter has to be a pion
1034         }
1035       
1036       //MC: like the nuclei via the specific energyloss in the TPC
1037       if (HasMC())
1038         {
1039           if (TMath::Abs(trackPosTest->GetTPCsignal() - expSignalPionP)/expSignalPionP < 0.4) posPion=kTRUE;
1040           if (TMath::Abs(trackNegTest->GetTPCsignal() - expSignalPionN)/expSignalPionN < 0.4) negPion=kTRUE;
1041         }
1042
1043       if (!(posProton==kTRUE)) continue;
1044       if (!(negPion==kTRUE)) continue;
1045
1046     //To avoid ghosts
1047
1048     if( !(trackPosTest->GetStatus() & AliESDtrack::kTPCrefit)){
1049       continue;
1050     }
1051
1052     if( !(trackNegTest->GetStatus() & AliESDtrack::kTPCrefit)){
1053       continue;
1054     }
1055
1056     if( trackPosTest->GetSign() >0 && trackNegTest->GetSign() <0){
1057       fV0MIs->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1058       fV0MIs->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter 
1059     }
1060
1061     if( trackPosTest->GetSign() <0 && trackNegTest->GetSign() >0){
1062       fV0MIs->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
1063       fV0MIs->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
1064     }
1065
1066     fV0MIs->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
1067
1068     TVector3 vecN(mn[0],mn[1],mn[2]);
1069     TVector3 vecP(mp[0],mp[1],mp[2]);
1070     TVector3 vecM(mm[0],mm[1],mm[2]);
1071     
1072     Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
1073     Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
1074     
1075     Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
1076       ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
1077     Double_t qt = vecP.Mag()*sin(thetaP);
1078
1079     fHistArmenterosPodolanski->Fill(alfa,qt);    //Armenteros-Podolanski calculation
1080
1081     TLorentzVector k0;
1082     TLorentzVector k0daugh1;
1083     TLorentzVector k0daugh2;
1084     TLorentzVector proton;
1085     TLorentzVector pq;
1086
1087     k0daugh1.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1088     k0daugh2.SetXYZM(mp[0],mp[1],mp[2],cPionMass);
1089     k0=k0daugh1+k0daugh2;
1090
1091     fV0MIs->ChangeMassHypothesis(3122);
1092     lInvMassLambda = fV0MIs->GetEffMass();
1093     lPtLambda = fV0MIs->Pt();
1094     lPzLambda = fV0MIs->Pz();
1095     lPLambda = fV0MIs->P();
1096
1097     trkArray->AddAt(trackPosTest,0);
1098     trkArray->AddAt(trackNegTest,1);
1099
1100     vertexer->SetVtxStart(esdVer1);
1101     AliESDVertex *decayVertex = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
1102
1103     dd[0]=fESD->GetPrimaryVertexSPD()->GetX()-decayVertex->GetX();
1104     dd[1]=fESD->GetPrimaryVertexSPD()->GetY()-decayVertex->GetY();
1105     dd[2]=fESD->GetPrimaryVertexSPD()->GetZ()-decayVertex->GetZ();
1106     
1107     decayLength=sqrt(dd[0]*dd[0]+dd[1]*dd[1]+dd[2]*dd[2]);
1108
1109       if (decayVertex == NULL) cout << "Lambda decay vtx pointer NULL" << endl;
1110     if (decayVertex) delete decayVertex;
1111
1112     TLorentzVector negPio1;
1113     TLorentzVector posProt1;        
1114     TLorentzVector posP;
1115     TLorentzVector posProt;
1116     TLorentzVector negK;
1117     TLorentzVector negPio;
1118     TLorentzVector negPi;
1119     TLorentzVector omega;
1120     TLorentzVector threeSum;
1121     TLorentzVector fourSum;
1122     TLorentzVector ppK;
1123     TLorentzVector posPiK;
1124     TLorentzVector negPiK;
1125     TLorentzVector kaon;
1126     TLorentzVector lambda;
1127     TLorentzVector lambdaH;
1128     TLorentzVector hDibaryon;
1129     TVector3 h;
1130     TVector3 h1;
1131
1132     Int_t mcStatus=0;
1133
1134     h.SetXYZ(-dd[0],-dd[1],-dd[2]);
1135
1136
1137     if (offl==1)fHistMassDPi->Fill(lInvMassLambda);
1138
1139     if (onl==1){
1140       fHistMassLPi->Fill(lInvMassLambda);
1141       
1142       negE.SetXYZM(mn[0],mn[1],mn[2],cElectronMass);
1143       posE.SetXYZM(mp[0],mp[1],mp[2],cElectronMass);
1144       photon=posE+negE;
1145       
1146       negPiK.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1147       posPiK.SetXYZM(mp[0],mp[1],mp[2],cPionMass);
1148       kaon=posPiK+negPiK;
1149       
1150       negPi.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1151       posP.SetXYZM(mp[0],mp[1],mp[2],cProtonMass);
1152       lambda=negPi+posP;
1153       lambdaH.SetXYZM(mm[0],mm[1],mm[2],cLambdaMass);
1154
1155       if (lInvMassLambda>1.1113&&lInvMassLambda<1.1202){
1156         
1157         if (!HasMC()){
1158           if (qt<-2.21*alfa*alfa+2.945*alfa-0.887) continue;
1159           if (qt>-2.21*alfa*alfa+2.945*alfa-0.873) continue;
1160           if (photon.M()<0.005) continue;
1161           if (kaon.M()>0.495 && kaon.M()<0.500 ) continue;
1162         }
1163         
1164
1165         fHistMassLambda->Fill(lInvMassLambda);
1166         //
1167         Bool_t isCorrectlyAssociatedLambda = kFALSE;
1168         Bool_t isPartialCorrectlyAssociatedLambda = kFALSE;
1169         Int_t labelAssociatedH=-1; 
1170         Int_t labelLambda=-1;
1171         //
1172         if (HasMC()) {
1173           Int_t labelPosTest = trackPosTest->GetLabel();
1174           TParticle *tparticleDaughter = stack->Particle(TMath::Abs(labelPosTest));
1175           Int_t labelMother = tparticleDaughter->GetFirstMother();
1176           TParticle *tparticleMother = stack->Particle(TMath::Abs(labelMother));
1177           
1178           Int_t labelOma = tparticleMother->GetFirstMother();
1179           TParticle *tparticleOma = stack->Particle(TMath::Abs(labelOma));
1180
1181           if ((tparticleOma->GetPdgCode() < 0) && TMath::Abs(tparticleMother->GetPdgCode())==pdgLambda){// check mother to be Lambda 
1182             Int_t labelFirstDaughter  = tparticleMother->GetDaughter(1);// Proton
1183             Int_t labelThirdDaughter  = tparticleMother->GetDaughter(0);// Pion
1184             
1185             TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1186             TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1187             
1188             if((tparticleFirstDaughter->GetPdgCode() == pdgProton && tparticleThirdDaughter->GetPdgCode()== pdgPionMinus) || 
1189                (tparticleFirstDaughter->GetPdgCode() == pdgPionMinus && tparticleThirdDaughter->GetPdgCode()== pdgProton)){ //daughter PDGs
1190               isPartialCorrectlyAssociatedLambda = kTRUE;
1191               labelLambda=labelMother;
1192             }
1193           }
1194           
1195           //H-Dibaryon
1196           if(tparticleOma->GetPdgCode() == pdgHDibaryon){ //check grandmother to be H PDG
1197             if (TMath::Abs(tparticleMother->GetPdgCode())==pdgLambda){// check mother to be Lambda 
1198               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(1);// Proton
1199               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(0);// Pion
1200                         
1201               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1202               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1203               
1204               if((tparticleFirstDaughter->GetPdgCode() == pdgProton && tparticleThirdDaughter->GetPdgCode()== pdgPionMinus) || 
1205                  (tparticleFirstDaughter->GetPdgCode() == pdgPionMinus && tparticleThirdDaughter->GetPdgCode()== pdgProton)){ //daughter PDGs
1206                 isCorrectlyAssociatedLambda = kTRUE;
1207                 labelAssociatedH=labelOma;
1208                 fHistMassLambdaFromH->Fill(lInvMassLambda);
1209                 fHistMassLambdaFromHtLorentz->Fill(lambda.M());
1210               }
1211             }
1212           }
1213         }
1214
1215         fHistProtonPIDLambda->Fill(trackPosTest->GetInnerParam()->GetP(), trackPosTest->GetTPCsignal());
1216         fHistPionPIDLambda->Fill(trackNegTest->GetInnerParam()->GetP(), trackNegTest->GetTPCsignal());
1217
1218         //---------------------------------------------------------
1219         // Proton track loop
1220         //---------------------------------------------------------
1221         fHistArmenterosPodolanskiCut->Fill(alfa,qt);
1222
1223         for (Int_t iTracksP = 0; iTracksP < fESD->GetNumberOfTracks(); iTracksP++) {
1224           AliESDtrack* trackP = dynamic_cast<AliESDtrack*> (fESD->GetTrack(iTracksP));
1225           if (trackP->GetSign()<0) continue;
1226           
1227           if (iTracksP==fV0MIs->GetPindex())continue;
1228           if (iTracksP==fV0MIs->GetNindex())continue;
1229                 
1230           if (!fEsdTrackCuts->AcceptTrack(trackP)) continue;
1231                 
1232           AliKFParticle protonKF( *(trackP), 2212);
1233
1234           if (!trackP->GetInnerParam()) continue;
1235           
1236           if (HasMC()) {
1237             pK=0.65;
1238           }
1239
1240           if (!HasMC()){
1241             if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackP, AliPID::kProton)) > 3) continue;
1242           }             
1243           
1244           fHistProtonPIDBb->Fill(trackP->GetInnerParam()->GetP(), trackP->GetTPCsignal());
1245           
1246           posProt.SetXYZM(trackP->Px(),trackP->Py(),trackP->Pz(),cProtonMass);
1247
1248           //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1249           //Pion Track loop!!!!!!!!!!!!!!!!!!!!!
1250           //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1251           for (Int_t iTracksN = iTracksP+1; iTracksN < fESD->GetNumberOfTracks(); iTracksN++) {
1252             AliESDtrack* trackN = dynamic_cast<AliESDtrack*> (fESD->GetTrack(iTracksN));
1253
1254             if (iTracksN==fV0MIs->GetPindex())continue;
1255             if (iTracksN==fV0MIs->GetNindex())continue;
1256             if (trackN->GetSign()>0) continue;
1257           
1258             if (!fEsdTrackCuts->AcceptTrack(trackN)) continue;
1259             if (!fESDtrackCutsV0->AcceptTrack(trackN)) continue;
1260
1261             Double_t bz = fESD->GetMagneticField();
1262            
1263             Double_t xthiss(0.0);
1264             Double_t xpp(0.0);
1265             Double_t dca = trackN->GetDCA(trackP,bz,xthiss,xpp);
1266             if (dca>0.5) continue;
1267
1268             negPi.SetXYZM(mn[0],mn[1],mn[2],cPionMass);
1269             posP.SetXYZM(mp[0],mp[1],mp[2],cProtonMass);
1270             negPio.SetXYZM(trackN->Px(),trackN->Py(),trackN->Pz(),cPionMass);
1271               
1272             threeSum=negPi+posP+negPio;
1273             lInvMassLambdaPi=threeSum.M();         
1274             
1275             fHistDC->Fill(decayLength*lInvMassLambda/lPLambda);
1276             
1277             AliKFParticle posPionKF( *(trackN) ,-211);
1278             
1279             if (!trackN->GetInnerParam()) continue;
1280             if (HasMC()) {
1281               pionK=0.7;
1282             }
1283             
1284           //  if (!HasMC()){
1285               if (TMath::Abs(fPIDtpcESD->NumberOfSigmasTPC(trackN, AliPID::kPion)) > 3) continue;
1286           //  }
1287             fHistPionPIDBb->Fill(trackN->GetInnerParam()->GetP(), trackN->GetTPCsignal());
1288
1289             trkArray1->AddAt(trackP,0);
1290             trkArray1->AddAt(trackN,1);
1291             
1292             vertexer1->SetVtxStart(esdVer1);
1293             AliESDVertex *decayVertex1 = (AliESDVertex*)vertexer1->VertexForSelectedESDTracks(trkArray1);
1294             
1295             dd1[0]=fESD->GetPrimaryVertexSPD()->GetX()-decayVertex1->GetX();
1296             dd1[1]=fESD->GetPrimaryVertexSPD()->GetY()-decayVertex1->GetY();
1297             dd1[2]=fESD->GetPrimaryVertexSPD()->GetZ()-decayVertex1->GetZ();
1298             
1299             decayLengthH=sqrt(dd1[0]*dd1[0]+dd1[1]*dd1[1]+dd1[2]*dd1[2]);
1300
1301               //            Double_t bz = fESD->GetMagneticField();
1302             
1303             trackP->PropagateToDCA(decayVertex1, bz, 10);
1304             trackN->PropagateToDCA(decayVertex1, bz, 10);
1305
1306               //            Double_t xthiss(0.0);
1307               //            Double_t xpp(0.0);
1308               //            Double_t dca = trackN->GetDCA(trackP,bz,xthiss,xpp);
1309
1310               if (decayVertex1 == NULL) cout << "Secondary decay vtx pointer NULL" << endl;
1311             if (decayVertex1) delete decayVertex1;
1312             h1.SetXYZ(-dd1[0],-dd1[1],-dd1[2]);
1313
1314             //      if (dca>1) continue;
1315               //            if (dca>0.1) continue;
1316
1317             fourSum=threeSum+posProt;
1318
1319             posProt1.SetXYZM(trackP->Px(),trackP->Py(),trackP->Pz(),cProtonMass);
1320             negPio1.SetXYZM(trackN->Px(),trackN->Py(),trackN->Pz(),cPionMass);
1321             hDibaryon=lambdaH+posProt1+negPio1;
1322             Double_t hPointingAngle = hDibaryon.Angle(h);
1323             Double_t pointingAngleH = hDibaryon.Angle(h1);
1324             Double_t decayAngleH = h.Angle(h1);
1325             TVector3 vecDist(dd[0]-dd1[0],dd[1]-dd1[1],dd[2]-dd1[2]);
1326             fHistMassLambdaPPi->Fill(hDibaryon.M());
1327             fHistHPointingAngle->Fill(pointingAngleH);
1328
1329             fHistMassHcentMult->Fill(hDibaryon.M(),triggertype,refMultTpc);
1330
1331             Double_t rapidity = hDibaryon.Rapidity();
1332             if(rapidity > 1.0 || rapidity < -1.0) continue;
1333
1334               //Double_t vec[16]={hDibaryon.M(), lInvMassLambda, pointingAngleH, alfa, qt, lPtLambda, hDibaryon.Pt(), posPionKF.GetDistanceFromVertex(primVtx), protonKF.GetDistanceFromVertex(primVtx), dca, protonKF.GetDistanceFromVertex(posPionKF), TMath::Cos(pointingAngleH), centrPerc, refMultTpc, mcStatus};
1335               //fHistNdim->Fill(vec);
1336
1337             fHistRap->Fill(rapidity);
1338             //if (pointingAngleH > 0.1) continue;
1339             if (pointingAngleH > 0.05) continue;
1340
1341             ///////////////////////////
1342             //MC part for Associated H
1343             ///////////////////////////
1344
1345             if (HasMC() && isCorrectlyAssociatedLambda) {
1346               Int_t labelP = trackP->GetLabel();
1347               TParticle *tparticleDaughter = stack->Particle(TMath::Abs(labelP));
1348               Int_t labelMother = tparticleDaughter->GetFirstMother();
1349               TParticle *tparticleMother = stack->Particle(TMath::Abs(labelMother));
1350               
1351               Int_t labelProton = trackP->GetLabel();
1352               Int_t labelPion = trackN->GetLabel();
1353               
1354               //H-Dibaryon
1355               if(tparticleMother->GetPdgCode() == pdgHDibaryon && labelAssociatedH==labelMother){ //check mother PDG
1356                 Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1357                 Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1358                 Int_t labelSecondDaughter = labelFirstDaughter +1;
1359
1360                 TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1361                 TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1362                 
1363                 TLorentzVector ppi;
1364                 TLorentzVector lpi;
1365                 TLorentzVector lP;
1366                 
1367                 ppi=posProt+negPio;
1368                 lpi=lambdaH+negPio;
1369                 lP=lambdaH+posProt;
1370
1371                 if((tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter)&&(tparticleSecondDaughter->GetPdgCode() == pdgProton||tparticleSecondDaughter->GetPdgCode() == pdgPionMinus) && labelProton==labelSecondDaughter) fHistMassPpi->Fill(ppi.M());
1372                 
1373                 if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter) fHistMassLpi->Fill(lpi.M());
1374                 
1375                 if(tparticleSecondDaughter->GetPdgCode() == pdgProton && labelProton==labelSecondDaughter) fHistMassLP->Fill(lP.M());
1376                 
1377                 if(tparticleSecondDaughter->GetPdgCode() == pdgProton && labelProton==labelSecondDaughter){//check second daughter PDG
1378                   if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus && labelPion==labelThirdDaughter){//check second daughter PDG
1379
1380                     fHistHDibaryonInvaMassAso->Fill(hDibaryon.M()); 
1381                     
1382                     Double_t distance01=vecDist.Mag();
1383                     fHistMCdcaPvtxDvtx->Fill(decayLengthH);
1384                     fHistMCdcaPvtxLvtx->Fill(decayLength);
1385                     fHistMCdcaDvtxLvtx->Fill(distance01);
1386                     fHistMCangleLH->Fill(hPointingAngle);
1387                     fHistMCdecayAngle->Fill(decayAngleH);
1388                     fHistMCpointingAngle->Fill(pointingAngleH);
1389                     fHistMCap->Fill(alfa,qt);
1390
1391                     fHistHilf1->Fill(posPionKF.GetDistanceFromVertex(primVtx));
1392                     fHistHilf2->Fill(protonKF.GetDistanceFromVertex(primVtx));
1393                     fHistHilf3->Fill(protonKF.GetDistanceFromVertex(posPionKF));
1394                     fHistHilf6->Fill(dca);
1395                     fHistPtvsYAso->Fill(hDibaryon.Pt(),hDibaryon.Rapidity());
1396                     fHistPtvsEtaAso->Fill(hDibaryon.Pt(),hDibaryon.Eta());
1397                     mcStatus=1;
1398                   }//end check for third daughter PDG
1399                 }//end check second daughter PDG
1400               }//end H-Dibaryon
1401             }//end MC
1402             
1403             //      cout<<"Trigger: "<<triggertype<<endl;
1404             fHistMassH->Fill(hDibaryon.M());
1405             fHistMassHcentMult->Fill(hDibaryon.M(),triggertype,refMultTpc);
1406             ppK=lambdaH+posProt;
1407             fHistMassLambdaP->Fill(ppK.M());
1408
1409               //fHistNdim = new THnSparseF("fHistNdim","THnS;InvMass, InvMassLambda, pointingAngle, armPoAlpha, armPoQt, pTL, pTH, d0p, d0n, dcaHd, dca, decayL, cosPA, centr, multi, mcinf;InvMassH", 16,binsD01,xminD01,xmaxD01);
1410
1411             Double_t vec[16]={hDibaryon.M(), lInvMassLambda, pointingAngleH, alfa, qt, lPtLambda, hDibaryon.Pt(), posPionKF.GetDistanceFromVertex(primVtx), protonKF.GetDistanceFromVertex(primVtx), dca, protonKF.GetDistanceFromVertex(posPionKF), TMath::Cos(pointingAngleH), centrPerc, static_cast<Double_t>(refMultTpc), static_cast<Double_t>(mcStatus)};
1412                       fHistNdim->Fill(vec);
1413
1414           }
1415         }
1416       }
1417     }    
1418   }
1419   
1420   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1421   //Pure MC Part!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1422   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1423
1424   // Monte Carlo for genenerated particles
1425   if (HasMC()) //MC loop  
1426     {
1427
1428       Int_t stackN = 0;
1429
1430       Double_t momentumPionGen[3]={0,0,0};
1431       Double_t momentumNucleonGen[3]={0,0,0};
1432       Double_t momentumLambdaGen[3]={0,0,0};
1433
1434       Double_t energyPionGen = 0;
1435       Double_t energyNucleonGen = 0;
1436       Double_t energyLambdaGen = 0;
1437
1438       Double_t transversMomentumMotherGen = 0;
1439       Double_t longitudinalMomentumMotherGen = 0;
1440       Double_t totalEnergyMotherGen = 0;
1441       
1442       Double_t rapidityGen = 2;
1443
1444       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1445         {
1446
1447           TParticle *tparticleMother = stack->Particle(stackN);
1448
1449           if(tparticleMother->GetPdgCode() == pdgLambda) fHistCount->Fill(16); 
1450
1451           //H-Dibaryon
1452           if(tparticleMother->GetPdgCode() == pdgHDibaryon) //check mother PDG
1453             {
1454               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1455               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1456               Int_t labelSecondDaughter = labelFirstDaughter +1;
1457
1458               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1459               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1460               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1461              
1462               if(tparticleFirstDaughter->GetPdgCode() == pdgLambda) //check first daughter PDG
1463                 {
1464                   if(tparticleSecondDaughter->GetPdgCode() == pdgProton)//check second daughter PDG
1465                     {
1466                       if(tparticleThirdDaughter->GetPdgCode() == pdgPionMinus)//check second daughter PDG
1467                         {                                     
1468                           momentumLambdaGen[0] = tparticleFirstDaughter->Px();
1469                           momentumLambdaGen[1] = tparticleFirstDaughter->Py();
1470                           momentumLambdaGen[2] = tparticleFirstDaughter->Pz();
1471                 
1472                           momentumNucleonGen[0] = tparticleSecondDaughter->Px();
1473                           momentumNucleonGen[1] = tparticleSecondDaughter->Py();
1474                           momentumNucleonGen[2] = tparticleSecondDaughter->Pz();
1475
1476                           momentumPionGen[0] = tparticleThirdDaughter->Px();
1477                           momentumPionGen[1] = tparticleThirdDaughter->Py();
1478                           momentumPionGen[2] = tparticleThirdDaughter->Pz();
1479
1480                           TLorentzVector lorentzVectorLambda;
1481                           TLorentzVector lorentzVectorProton;
1482                           TLorentzVector lorentzVectorPion;
1483                           TLorentzVector lorentzVectorHDibaryon;
1484                           
1485                           lorentzVectorLambda.SetXYZM(momentumLambdaGen[0],momentumLambdaGen[1],momentumLambdaGen[2],1.115);
1486                           lorentzVectorProton.SetXYZM(momentumNucleonGen[0],momentumNucleonGen[1],momentumNucleonGen[2],protonMass);
1487                           lorentzVectorPion.SetXYZM(momentumPionGen[0],momentumPionGen[1],momentumPionGen[2],pionMass);
1488                           
1489                           lorentzVectorHDibaryon = lorentzVectorLambda + lorentzVectorProton + lorentzVectorPion;
1490                           rapidityGen=lorentzVectorHDibaryon.Rapidity();
1491                           transversMomentumMotherGen = lorentzVectorHDibaryon.Pt();
1492                           longitudinalMomentumMotherGen = lorentzVectorHDibaryon.Pz();
1493                           totalEnergyMotherGen = lorentzVectorHDibaryon.Energy();
1494
1495                           if(rapidityGen > 1.0 || rapidityGen < -1 ) continue;
1496                           //lorentzVectorLambda
1497                           fHistHDibaryonInvaMassGen->Fill(lorentzVectorHDibaryon.M()); 
1498                           if (lorentzVectorLambda.Rapidity()  > 1.0 || lorentzVectorLambda.Rapidity() < -1) continue;
1499                           if (lorentzVectorProton.Rapidity()  > 1.0 ||      lorentzVectorProton.Rapidity() < -1) continue;
1500
1501                           if (lorentzVectorPion.Rapidity()  > 1.0 ||  lorentzVectorPion.Rapidity() < -1) continue;
1502                           fHistHDibaryonInvaMassGenRes->Fill(lorentzVectorHDibaryon.M());
1503                           fHistPtvsEtaGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Eta());
1504                           fHistPtvsYGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Rapidity());
1505                           fHistPtvsEtaGen->Fill(lorentzVectorHDibaryon.Pt(),lorentzVectorHDibaryon.Eta());
1506                           fHistCount->Fill(11);
1507                         }//end of check third daughter PDG
1508                     }//end of check second daughter PDG
1509                 }//end of check first daughter PDG
1510             }//end of H-Dibaryon
1511
1512           //Anti-H-Dibaryon
1513           if(tparticleMother->GetPdgCode() == pdgAntiHDibaryon) //check mother PDG
1514             {
1515               Int_t labelFirstDaughter  = tparticleMother->GetDaughter(0);
1516               Int_t labelThirdDaughter  = tparticleMother->GetDaughter(1);
1517               Int_t labelSecondDaughter = labelFirstDaughter +1;
1518
1519               TParticle *tparticleFirstDaughter  = stack->Particle(TMath::Abs(labelFirstDaughter));
1520               TParticle *tparticleSecondDaughter = stack->Particle(TMath::Abs(labelSecondDaughter));
1521               TParticle *tparticleThirdDaughter  = stack->Particle(TMath::Abs(labelThirdDaughter));
1522
1523               if(tparticleFirstDaughter->GetPdgCode() == pdgAntiLambda) //check first daughter PDG
1524                 {
1525                   if(tparticleSecondDaughter->GetPdgCode() == pdgAntiProton)//check second daughter PDG
1526                     {
1527                       if(tparticleThirdDaughter->GetPdgCode() == pdgPionPlus)//check second daughter PDG
1528                         {                
1529                           momentumLambdaGen[0] = tparticleFirstDaughter->Px();
1530                           momentumLambdaGen[1] = tparticleFirstDaughter->Py();
1531                           momentumLambdaGen[2] = tparticleFirstDaughter->Pz();
1532                 
1533                           momentumNucleonGen[0] = tparticleSecondDaughter->Px();
1534                           momentumNucleonGen[1] = tparticleSecondDaughter->Py();
1535                           momentumNucleonGen[2] = tparticleSecondDaughter->Pz();
1536
1537                           momentumPionGen[0] = tparticleThirdDaughter->Px();
1538                           momentumPionGen[1] = tparticleThirdDaughter->Py();
1539                           momentumPionGen[2] = tparticleThirdDaughter->Pz();
1540                           
1541                           energyLambdaGen  = tparticleFirstDaughter->Energy();
1542                           energyNucleonGen = tparticleSecondDaughter->Energy();
1543                           energyPionGen    = tparticleThirdDaughter->Energy();
1544                            
1545                           TLorentzVector lorentzVectorLambda;
1546                           TLorentzVector lorentzVectorProton;
1547                           TLorentzVector lorentzVectorPion;
1548                           TLorentzVector lorentzVectorHDibaryon;
1549                           
1550                           lorentzVectorLambda.SetXYZM(momentumLambdaGen[0],momentumLambdaGen[1],momentumLambdaGen[2],1.115);
1551                           lorentzVectorProton.SetXYZM(momentumNucleonGen[0],momentumNucleonGen[1],momentumNucleonGen[2],protonMass);
1552                           lorentzVectorPion.SetXYZM(momentumPionGen[0],momentumPionGen[1],momentumPionGen[2],pionMass);
1553                           
1554                           lorentzVectorHDibaryon = lorentzVectorLambda + lorentzVectorProton + lorentzVectorPion;
1555
1556                           rapidityGen=lorentzVectorHDibaryon.Rapidity();
1557                           if(rapidityGen > 1.0 || rapidityGen < -1 ) continue;
1558                           fHistAntiHDibaryonInvaMassGen->Fill(lorentzVectorHDibaryon.M()); 
1559                         }//end of check third daughter PDG
1560                     }//end of check second daughter PDG
1561                 }//end of check first daughter PDG
1562             }//end of Anti-H-Dibaryon
1563         }      
1564     }//end MC
1565
1566   // Post output data.
1567   PostData(1,fHistList);
1568   //PostData(0,fHistList);
1569
1570     if (listCrossV0 == NULL) return;
1571
1572   if (listCrossV0) delete listCrossV0;
1573   if (esdVer1) delete esdVer1;
1574   if (vertexer) delete vertexer;
1575   if (vertexer1) delete vertexer1;
1576   if (trkArray) delete trkArray;
1577   if (trkArray1) delete trkArray1;
1578 }
1579
1580 //________________________________________________________________________
1581 void AliAnalysisTaskHdibaryonLPpi::Terminate(Option_t *) 
1582 {
1583   // Draw result to the screen
1584   // Called once at the end of the query
1585
1586 }
1587