New task for Lc->pKpi (Rossella)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSELambdac.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /////////////////////////////////////////////////////////////
17 //
18 // AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
19 // decay candidates with the MC truth.
20 // Authors: Renu Bala, bala@to.infn.it
21 // F. Prino, prino@to.infn.it
22 // G. Ortona, ortona@to.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TClonesArray.h>
26 #include <TNtuple.h>
27 #include <TCanvas.h>
28 #include <TList.h>
29 #include <TString.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TDatabasePDG.h>
33
34 #include "AliAnalysisManager.h"
35 #include "AliAODHandler.h"
36 #include "AliAODEvent.h"
37 #include "AliAODVertex.h"
38 #include "AliAODTrack.h"
39 #include "AliAODMCHeader.h"
40 #include "AliAODMCParticle.h"
41 #include "AliAODRecoDecayHF3Prong.h"
42 #include "AliAnalysisVertexingHF.h"
43 #include "AliAnalysisTaskSE.h"
44 #include "AliAnalysisTaskSELambdac.h"
45 #include "AliKFParticle.h"
46 #include "AliAODPidHF.h"
47 #include "AliRDHFCutsLctopKpi.h"
48 #include "AliRDHFCuts.h"
49
50 ClassImp(AliAnalysisTaskSELambdac)
51
52
53 //________________________________________________________________________
54 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac():
55 AliAnalysisTaskSE(),
56 fOutput(0), 
57 fHistNEvents(0),
58 fTPCSignal3Sigma(0),
59 fTPCSignal3SigmaReK(0),
60 fTPCSignal3SigmaRep(0),
61 fTPCSignal2Sigma(0),
62 fTPCSignal2SigmaReK(0),
63 fTPCSignal2SigmaRep(0),
64 fNtupleLambdac(0),
65 fUpmasslimit(2.486),
66 fLowmasslimit(2.086),
67 fNPtBins(0),
68 fRDCutsAnalysis(0),
69 fRDCutsProduction(0),
70 fListCuts(0),
71 fFillNtuple(kFALSE),
72 fReadMC(kFALSE),
73 fMCPid(kFALSE),
74 fRealPid(kFALSE),
75 fResPid(kTRUE),
76 fUseKF(kFALSE),
77 fVHF(0)
78 {
79    // Default constructor
80 }
81
82 //________________________________________________________________________
83 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac(const char *name,Bool_t fillNtuple,AliRDHFCutsLctopKpi *lccutsana,AliRDHFCutsLctopKpi *lccutsprod):
84 AliAnalysisTaskSE(name),
85 fOutput(0),
86 fHistNEvents(0),
87 fTPCSignal3Sigma(0),
88 fTPCSignal3SigmaReK(0),
89 fTPCSignal3SigmaRep(0),
90 fTPCSignal2Sigma(0),
91 fTPCSignal2SigmaReK(0),
92 fTPCSignal2SigmaRep(0),
93 fNtupleLambdac(0),
94 fUpmasslimit(2.486),
95 fLowmasslimit(2.086),
96 fNPtBins(0),
97 fRDCutsAnalysis(lccutsana),
98 fRDCutsProduction(lccutsprod),
99 fListCuts(0),
100 fFillNtuple(fillNtuple),
101 fReadMC(kFALSE),
102 fMCPid(kFALSE),
103 fRealPid(kTRUE),
104 fResPid(kFALSE),
105 fUseKF(kFALSE),
106 fVHF(0)
107 {
108   //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
109   //SetPtBinLimit(5, ptlim);
110    SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
111   // Default constructor
112    // Output slot #1 writes into a TList container
113   DefineOutput(1,TList::Class());  //My private output
114   DefineOutput(2,TList::Class());
115
116   if(fFillNtuple){
117     // Output slot #2 writes into a TNtuple container
118     DefineOutput(3,TNtuple::Class());  //My private output
119   }
120 }
121
122 //________________________________________________________________________
123 AliAnalysisTaskSELambdac::~AliAnalysisTaskSELambdac()
124 {
125   // Destructor
126   if (fOutput) {
127     delete fOutput;
128     fOutput = 0;
129   }
130   if (fVHF) {
131     delete fVHF;
132     fVHF = 0;
133   }
134   
135  if(fRDCutsAnalysis){
136     delete fRDCutsAnalysis;
137     fRDCutsAnalysis = 0;
138  }
139  if(fRDCutsProduction){
140     delete fRDCutsProduction;
141     fRDCutsProduction = 0;
142  }
143
144  if (fListCuts) {
145    delete fListCuts;
146    fListCuts = 0;
147  }
148 }  
149 //_________________________________________________________________
150 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t range){
151   fUpmasslimit = 2.286+range;
152   fLowmasslimit = 2.286-range;
153 }
154 //_________________________________________________________________
155 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t lowlimit, Float_t uplimit){
156   if(uplimit>lowlimit)
157     {
158       fUpmasslimit = lowlimit;
159       fLowmasslimit = uplimit;
160     }
161 }
162
163
164 //________________________________________________________________________
165 void AliAnalysisTaskSELambdac::SetPtBinLimit(Int_t n, Float_t* lim){
166   // define pt bins for analysis
167   if(n>kMaxPtBins){
168     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
169     fNPtBins=kMaxPtBins;
170     fArrayBinLimits[0]=0.;
171     fArrayBinLimits[1]=2.;
172     fArrayBinLimits[2]=3.;
173     fArrayBinLimits[3]=5.;
174     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
175   }else{
176     fNPtBins=n-1;
177     fArrayBinLimits[0]=lim[0];
178     for(Int_t i=1; i<fNPtBins+1; i++) 
179       if(lim[i]>fArrayBinLimits[i-1]){
180         fArrayBinLimits[i]=lim[i];
181       }
182       else {
183         fArrayBinLimits[i]=fArrayBinLimits[i-1];
184       }
185     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
186   }
187   if(fDebug > 1){
188     printf("Number of Pt bins = %d\n",fNPtBins);
189     for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);    
190   }
191 }
192 //_________________________________________________________________
193 Double_t  AliAnalysisTaskSELambdac::GetPtBinLimit(Int_t ibin){
194   if(ibin>fNPtBins)return -1;
195   return fArrayBinLimits[ibin];
196
197
198 //_________________________________________________________________
199 void AliAnalysisTaskSELambdac::Init()
200 {
201   // Initialization
202
203   if(fDebug > 1) printf("AnalysisTaskSELambdac::Init() \n");
204
205   fListCuts=new TList();
206   AliRDHFCutsLctopKpi *production = new AliRDHFCutsLctopKpi();
207     production=fRDCutsProduction;
208
209   AliRDHFCutsLctopKpi *analysis = new AliRDHFCutsLctopKpi();
210   analysis=fRDCutsAnalysis;
211
212    fListCuts->Add(analysis);
213    fListCuts->Add(production);
214    PostData(2,fListCuts);
215   return;
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskSELambdac::UserCreateOutputObjects()
220 {
221   // Create the output container
222   //
223   if(fDebug > 1) printf("AnalysisTaskSELambdac::UserCreateOutputObjects() \n");
224
225   // Several histograms are more conveniently managed in a TList
226   fOutput = new TList();
227   fOutput->SetOwner();
228   fOutput->SetName("OutputHistos");
229
230   TString hisname;
231   Int_t index=0;
232   Int_t indexLS=0;
233   for(Int_t i=0;i<fNPtBins;i++){
234
235     index=GetHistoIndex(i);
236     indexLS=GetLSHistoIndex(i);
237
238     hisname.Form("hMassPt%d",i);
239     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
240     fMassHist[index]->Sumw2();
241     hisname.Form("hCosPAllPt%d",i);
242     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
243     fCosPHist[index]->Sumw2();
244     hisname.Form("hDLenAllPt%d",i);
245     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
246     fDLenHist[index]->Sumw2();
247     hisname.Form("hSumd02AllPt%d",i);
248     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
249     fSumd02Hist[index]->Sumw2();
250     hisname.Form("hSigVertAllPt%d",i);
251     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
252     fSigVertHist[index]->Sumw2();
253     hisname.Form("hPtMaxAllPt%d",i);
254     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
255     fPtMaxHist[index]->Sumw2();
256
257     hisname.Form("hDCAAllPt%d",i);
258     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
259     fDCAHist[index]->Sumw2();
260
261
262
263     hisname.Form("hMassPt%dTC",i);
264     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
265     fMassHistTC[index]->Sumw2();
266
267
268
269
270     
271     hisname.Form("hCosPAllPt%dLS",i);
272     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
273     fCosPHistLS[index]->Sumw2();
274     hisname.Form("hDLenAllPt%dLS",i);
275     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
276     fDLenHistLS[index]->Sumw2();
277     hisname.Form("hSumd02AllPt%dLS",i);
278     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
279     fSumd02HistLS[index]->Sumw2();
280     hisname.Form("hSigVertAllPt%dLS",i);
281     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
282     fSigVertHistLS[index]->Sumw2();
283     hisname.Form("hPtMaxAllPt%dLS",i);
284     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
285     fPtMaxHistLS[index]->Sumw2();
286     
287     hisname.Form("hDCAAllPt%dLS",i);
288     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
289     fDCAHistLS[index]->Sumw2();
290     
291     hisname.Form("hLSPt%dLC",i);
292     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
293     fMassHistLS[indexLS]->Sumw2();
294     
295     hisname.Form("hLSPt%dTC",i);
296     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
297     fMassHistLSTC[indexLS]->Sumw2();
298
299
300     
301     index=GetSignalHistoIndex(i);    
302     indexLS++;
303     hisname.Form("hSigPt%d",i);
304     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
305     fMassHist[index]->Sumw2();
306     hisname.Form("hCosPSigPt%d",i);
307     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
308     fCosPHist[index]->Sumw2();
309     hisname.Form("hDLenSigPt%d",i);
310     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
311     fDLenHist[index]->Sumw2();
312     hisname.Form("hSumd02SigPt%d",i);
313     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
314     fSumd02Hist[index]->Sumw2();
315     hisname.Form("hSigVertSigPt%d",i);
316     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
317     fSigVertHist[index]->Sumw2();
318     hisname.Form("hPtMaxSigPt%d",i);
319     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
320     fPtMaxHist[index]->Sumw2();    
321
322     hisname.Form("hDCASigPt%d",i);
323     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
324     fDCAHist[index]->Sumw2();    
325
326
327     hisname.Form("hSigPt%dTC",i);
328     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
329     fMassHistTC[index]->Sumw2();
330
331     hisname.Form("hLSPt%dLCnw",i);
332     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
333     fMassHistLS[indexLS]->Sumw2();
334     hisname.Form("hLSPt%dTCnw",i);
335     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
336     fMassHistLSTC[indexLS]->Sumw2();
337
338
339     
340     hisname.Form("hCosPSigPt%dLS",i);
341     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
342     fCosPHistLS[index]->Sumw2();
343     hisname.Form("hDLenSigPt%dLS",i);
344     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
345     fDLenHistLS[index]->Sumw2();
346     hisname.Form("hSumd02SigPt%dLS",i);
347     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
348     fSumd02HistLS[index]->Sumw2();
349     hisname.Form("hSigVertSigPt%dLS",i);
350     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
351     fSigVertHistLS[index]->Sumw2();
352     hisname.Form("hPtMaxSigPt%dLS",i);
353     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
354     fPtMaxHistLS[index]->Sumw2();
355
356     hisname.Form("hDCASigPt%dLS",i);
357     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
358     fDCAHistLS[index]->Sumw2();
359     
360
361
362     index=GetBackgroundHistoIndex(i); 
363     indexLS++;
364     hisname.Form("hBkgPt%d",i);
365     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
366     fMassHist[index]->Sumw2();
367     hisname.Form("hCosPBkgPt%d",i);
368     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
369     fCosPHist[index]->Sumw2();
370     hisname.Form("hDLenBkgPt%d",i);
371     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
372     fDLenHist[index]->Sumw2();
373     hisname.Form("hSumd02BkgPt%d",i);
374     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
375     fSumd02Hist[index]->Sumw2();
376     hisname.Form("hSigVertBkgPt%d",i);
377     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
378     fSigVertHist[index]->Sumw2();
379     hisname.Form("hPtMaxBkgPt%d",i);
380     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
381     fPtMaxHist[index]->Sumw2();
382
383     hisname.Form("hDCABkgPt%d",i);
384     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
385     fDCAHist[index]->Sumw2();
386
387
388     hisname.Form("hBkgPt%dTC",i);
389     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
390     fMassHistTC[index]->Sumw2();
391
392     hisname.Form("hLSPt%dLCntrip",i);
393     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
394     fMassHistLS[indexLS]->Sumw2();
395     hisname.Form("hLSPt%dTCntrip",i);
396     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
397     fMassHistLSTC[indexLS]->Sumw2();
398
399     
400     hisname.Form("hCosPBkgPt%dLS",i);
401     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
402     fCosPHistLS[index]->Sumw2();
403     hisname.Form("hDLenBkgPt%dLS",i);
404     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
405     fDLenHistLS[index]->Sumw2();
406     hisname.Form("hSumd02BkgPt%dLS",i);
407     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
408     fSumd02HistLS[index]->Sumw2();
409     hisname.Form("hSigVertBkgPt%dLS",i);
410     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
411     fSigVertHistLS[index]->Sumw2();
412     hisname.Form("hPtMaxBkgPt%dLS",i);
413     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
414     fPtMaxHistLS[index]->Sumw2();
415     hisname.Form("hDCABkgPt%dLS",i);
416     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
417     fDCAHistLS[index]->Sumw2();
418     
419
420     indexLS++;
421     hisname.Form("hLSPt%dLCntripsinglecut",i);
422     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
423     fMassHistLS[indexLS]->Sumw2();
424     hisname.Form("hLSPt%dTCntripsinglecut",i);
425     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
426     fMassHistLSTC[indexLS]->Sumw2();
427
428     indexLS++;
429     hisname.Form("hLSPt%dLCspc",i);
430     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
431     fMassHistLS[indexLS]->Sumw2();
432     hisname.Form("hLSPt%dTCspc",i);
433     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
434     fMassHistLSTC[indexLS]->Sumw2();
435   }
436
437   fTPCSignal3Sigma=new TH2F("fTPCSignal3Sigma","fTPCSignal3Sigma",100,0, 10, 100, 0, 100);
438   fTPCSignal3SigmaReK=new TH2F("fTPCSignal3SigmaReK","fTPCSignal3SigmaRe",100,0, 10, 100, 0, 100);
439   fTPCSignal3SigmaRep=new TH2F("fTPCSignal3SigmaRep","fTPCSignal3SigmaRe",100,0, 10, 100, 0, 100);
440   fTPCSignal2Sigma=new TH2F("fTPCSignal2Sigma","fTPCSignal2Sigma",100,0, 10, 100, 0, 100);
441   fTPCSignal2SigmaReK=new TH2F("fTPCSignal2SigmaReK","fTPCSignal2SigmaRe",100,0, 10, 100, 0, 100);
442   fTPCSignal2SigmaRep=new TH2F("fTPCSignal2SigmaRep","fTPCSignal2SigmaRe",100,0, 10, 100, 0, 100);
443   
444   for(Int_t i=0; i<3*fNPtBins; i++){
445     fOutput->Add(fMassHist[i]);
446     fOutput->Add(fMassHistTC[i]);
447     fOutput->Add(fCosPHist[i]);
448     fOutput->Add(fDLenHist[i]);
449     fOutput->Add(fSumd02Hist[i]);
450     fOutput->Add(fSigVertHist[i]);
451     fOutput->Add(fPtMaxHist[i]);
452     fOutput->Add(fDCAHist[i]);
453   }
454
455     fOutput->Add(fTPCSignal3Sigma);
456     fOutput->Add(fTPCSignal3SigmaReK);
457     fOutput->Add(fTPCSignal3SigmaRep);
458     fOutput->Add(fTPCSignal2Sigma);
459     fOutput->Add(fTPCSignal2SigmaReK);
460     fOutput->Add(fTPCSignal2SigmaRep);
461   
462
463   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
464   fHistNEvents->Sumw2();
465   fHistNEvents->SetMinimum(0);
466   fOutput->Add(fHistNEvents);
467   
468   
469
470   if(fFillNtuple){
471     //OpenFile(3); // 2 is the slot number of the ntuple
472    
473     fNtupleLambdac = new TNtuple("fNtupleLambdac","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");  
474     
475   }
476   
477   return;
478 }
479
480 //________________________________________________________________________
481 void AliAnalysisTaskSELambdac::UserExec(Option_t */*option*/)
482 {
483   // Execute analysis for current event:
484   // heavy flavor candidates association to MC truth
485
486    AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
487    //tmp!
488  fHistNEvents->Fill(0); // count event
489   // Post the data already here
490   PostData(1,fOutput);
491
492   TClonesArray *array3Prong = 0;
493   TClonesArray *arrayLikeSign =0;
494   if(!aod && AODEvent() && IsStandardAOD()) {
495     // In case there is an AOD handler writing a standard AOD, use the AOD 
496     // event in memory rather than the input (ESD) event.    
497     aod = dynamic_cast<AliAODEvent*> (AODEvent());
498     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
499     // have to taken from the AOD event hold by the AliAODExtension
500     AliAODHandler* aodHandler = (AliAODHandler*) 
501       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
502     if(aodHandler->GetExtensions()) {
503       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
504       AliAODEvent *aodFromExt = ext->GetAOD();
505       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
506       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
507     }
508   } else {
509     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
510     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
511   }
512
513   if(!array3Prong) {
514     printf("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
515     return;
516   }
517   if(!arrayLikeSign) {
518     printf("AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
519   //  return;
520   }
521
522  
523   TClonesArray *arrayMC=0;
524   AliAODMCHeader *mcHeader=0;
525
526   // AOD primary vertex
527   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
528   
529   // load MC particles
530   if(fReadMC){
531     
532     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
533     if(!arrayMC) {
534       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
535       return;
536     }
537     
538   // load MC header
539     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
540     if(!mcHeader) {
541     printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
542     return;
543     }
544   }
545   
546   Int_t n3Prong = array3Prong->GetEntriesFast();
547   
548   
549   Int_t nOS=0;
550   Int_t index;
551   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
552     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
553     
554     
555     Bool_t unsetvtx=kFALSE;
556     if(!d->GetOwnPrimaryVtx()){
557       d->SetOwnPrimaryVtx(vtx1);
558       unsetvtx=kTRUE;
559     }
560
561     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) {
562       Int_t iPtBin = -1;
563       Double_t ptCand = d->Pt();
564       
565       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
566         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
567       }
568
569       Bool_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
570       Int_t labDp=-1;
571       Float_t deltaPx=0.;
572       Float_t deltaPy=0.;
573       Float_t deltaPz=0.;
574       Float_t truePt=0.;
575       Float_t xDecay=0.;
576       Float_t yDecay=0.;
577       Float_t zDecay=0.;
578       Float_t pdgCode=-2;
579       if(fReadMC){
580         labDp = MatchToMCLambdac(d,arrayMC);   
581         //
582         if(labDp>0){
583           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
584           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));  
585           deltaPx=partDp->Px()-d->Px();
586           deltaPy=partDp->Py()-d->Py();
587           deltaPz=partDp->Pz()-d->Pz();
588           truePt=partDp->Pt();
589           xDecay=dg0->Xv();       
590           yDecay=dg0->Yv();       
591           zDecay=dg0->Zv();
592           pdgCode=TMath::Abs(partDp->GetPdgCode());
593         }else{
594           pdgCode=-1;
595         }
596       }
597
598       Double_t invMasspKpi=-1.;
599       Double_t invMasspiKp=-1.;
600       Int_t pdgs[3]={0,0,0};
601       Double_t field=aod->GetMagneticField();
602       if(fReadMC && fMCPid){
603        
604        if(IspKpiMC(d,arrayMC)) {
605         invMasspKpi=d->InvMassLcpKpi();
606         if(fUseKF){
607          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
608          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
609         }
610        }
611        IspiKpReal(d);
612        IspiKpResonant(d,field);
613        if(IspiKpMC(d,arrayMC)) {
614         invMasspiKp=d->InvMassLcpiKp();
615         if(fUseKF){
616          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
617          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
618         }
619        }
620       }
621       if(fRealPid){
622        if(IspKpiReal(d)) {
623         invMasspKpi=d->InvMassLcpKpi();
624         if(fUseKF){
625          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
626          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
627         }
628        }
629        if(IspiKpReal(d)) {
630         invMasspiKp=d->InvMassLcpiKp();
631         if(fUseKF){
632          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
633          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
634         }
635        }
636       }
637       if(fResPid){
638        if(IspKpiResonant(d,field)) {
639         invMasspKpi=d->InvMassLcpKpi();
640         if(fUseKF){
641          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
642          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
643         }
644        }
645        if(IspiKpResonant(d,field)) {
646         invMasspiKp=d->InvMassLcpiKp();
647         if(fUseKF){
648          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
649          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
650         }
651        }
652       }
653       if(!fResPid && !fRealPid && !fMCPid){
654        invMasspiKp=d->InvMassLcpiKp(); 
655        if(fUseKF){
656          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
657          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
658         }
659        invMasspKpi=d->InvMassLcpKpi();
660        if(fUseKF){
661         pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
662         if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
663        }
664       }
665
666       if(invMasspiKp<0. && invMasspKpi<0.) continue;
667
668
669       Float_t tmp[24];
670       if(fFillNtuple){            
671         tmp[0]=pdgCode;
672         tmp[1]=deltaPx;
673         tmp[2]=deltaPy;
674         tmp[3]=deltaPz;
675         tmp[4]=truePt;
676         tmp[5]=xDecay;    
677         tmp[6]=yDecay;    
678         tmp[7]=zDecay;    
679         tmp[8]=d->PtProng(0);
680         tmp[9]=d->PtProng(1);
681         tmp[10]=d->PtProng(2);
682         tmp[11]=d->Pt();
683         tmp[12]=d->CosPointingAngle();
684         tmp[13]=d->DecayLength();
685         tmp[14]=d->Xv();
686         tmp[15]=d->Yv();
687         tmp[16]=d->Zv();
688         if(invMasspiKp>0.) tmp[17]=invMasspiKp;
689         if(invMasspKpi>0.) tmp[17]=invMasspKpi;
690         tmp[18]=d->GetSigmaVert();
691         tmp[19]=d->Getd0Prong(0);
692         tmp[20]=d->Getd0Prong(1);
693         tmp[21]=d->Getd0Prong(2);
694         tmp[22]=d->GetDCA();
695         tmp[23]=d->Prodd0d0(); 
696         fNtupleLambdac->Fill(tmp);
697         PostData(3,fNtupleLambdac);
698       }
699       Double_t dlen=d->DecayLength();
700       Double_t cosp=d->CosPointingAngle();
701       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
702       Double_t dca=d->GetDCA();      
703 Double_t ptmax=0;
704       for(Int_t i=0;i<3;i++){
705         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
706       }
707       if(iPtBin>=0){
708       
709         index=GetHistoIndex(iPtBin);
710         if(invMasspiKp>0. && invMasspKpi>0.){
711         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
712         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
713         }else{
714          if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
715          if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
716         }
717         fCosPHist[index]->Fill(cosp);
718         fDLenHist[index]->Fill(dlen);
719         fSumd02Hist[index]->Fill(sumD02);
720         fPtMaxHist[index]->Fill(ptmax);
721         fDCAHist[index]->Fill(dca);
722         
723         if(passTightCuts){
724          if(invMasspiKp>0. && invMasspKpi>0.){
725           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
726           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
727          }else{
728           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
729           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
730          }
731         }
732         
733         if(fReadMC){
734           if(labDp>0) {
735             index=GetSignalHistoIndex(iPtBin);
736             if(invMasspiKp>0. && invMasspKpi>0.){
737              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
738              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
739             }else{
740              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
741              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
742             }
743             fCosPHist[index]->Fill(cosp);
744             fDLenHist[index]->Fill(dlen);
745             fSumd02Hist[index]->Fill(sumD02);
746             fPtMaxHist[index]->Fill(ptmax);
747             fDCAHist[index]->Fill(dca);
748             if(passTightCuts){
749              if(invMasspiKp>0. && invMasspKpi>0.){
750               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
751               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
752              }else{
753               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
754               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
755              }
756             }
757             
758           }else{
759             index=GetBackgroundHistoIndex(iPtBin);
760             if(invMasspiKp>0. && invMasspKpi>0.){
761              fMassHist[index]->Fill(invMasspiKp,0.5);
762              fMassHist[index]->Fill(invMasspKpi,0.5);
763             }else{
764              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
765              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
766            }
767             fCosPHist[index]->Fill(cosp);
768             fDLenHist[index]->Fill(dlen);
769             fSumd02Hist[index]->Fill(sumD02);
770             fPtMaxHist[index]->Fill(ptmax);
771             fDCAHist[index]->Fill(dca);
772             if(passTightCuts){
773             if(invMasspiKp>0. && invMasspKpi>0.){
774                fMassHistTC[index]->Fill(invMasspiKp,0.5);
775                fMassHistTC[index]->Fill(invMasspKpi,0.5);
776               }else{
777                if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
778                if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
779               }
780             }   
781           }
782         }
783       }
784       /*
785       //start OS analysis
786       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
787       fHistOS->Fill(d->InvMassDplus());
788       */
789       nOS++;
790     }
791     if(unsetvtx) d->UnsetOwnPrimaryVtx();
792   }
793  
794   PostData(1,fOutput);    
795   return;
796 }
797
798
799
800 //________________________________________________________________________
801 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
802 {
803   // Terminate analysis
804   //
805   if(fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
806
807   fOutput = dynamic_cast<TList*> (GetOutputData(1));
808   if (!fOutput) {     
809     printf("ERROR: fOutput not available\n");
810     return;
811   }
812  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
813
814  TString hisname;
815  Int_t index=0;
816  
817
818  //Int_t indexLS=0;
819  for(Int_t i=0;i<fNPtBins;i++){
820     index=GetHistoIndex(i);
821     hisname.Form("hMassPt%d",i);
822     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
823      hisname.Form("hCosPAllPt%d",i);
824     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
825      hisname.Form("hDLenAllPt%d",i);
826     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
827      hisname.Form("hSumd02AllPt%d",i);
828      fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
829      hisname.Form("hSigVertAllPt%d",i);
830      fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
831      hisname.Form("hPtMaxAllPt%d",i);
832      fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
833      hisname.Form("hDCAAllPt%d",i);
834      fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
835      hisname.Form("hMassPt%dTC",i);
836     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
837     
838     index=GetSignalHistoIndex(i);    
839     hisname.Form("hSigPt%d",i);
840     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
841     hisname.Form("hCosPSigPt%d",i);
842     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
843     hisname.Form("hDLenSigPt%d",i);
844     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
845     hisname.Form("hSumd02SigPt%d",i);
846     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
847     hisname.Form("hSigVertSigPt%d",i);
848     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
849     hisname.Form("hPtMaxSigPt%d",i);
850     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
851     hisname.Form("hDCASigPt%d",i);
852     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
853     
854     hisname.Form("hSigPt%dTC",i);
855     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
856     
857     index=GetBackgroundHistoIndex(i); 
858     hisname.Form("hBkgPt%d",i);
859     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
860     hisname.Form("hCosPBkgPt%d",i);
861     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
862     hisname.Form("hDLenBkgPt%d",i);
863     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
864     hisname.Form("hSumd02BkgPt%d",i);
865     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
866     hisname.Form("hSigVertBkgPt%d",i);
867     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
868     hisname.Form("hPtMaxBkgPt%d",i);
869     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
870     hisname.Form("hDCABkgPt%d",i);
871     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
872     hisname.Form("hBkgPt%dTC",i);
873     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
874  
875  }
876     fTPCSignal3Sigma=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal3Sigma));
877     fTPCSignal3SigmaReK=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal3SigmaReK));
878     fTPCSignal3SigmaRep=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal3SigmaRep));
879     fTPCSignal2Sigma=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal2Sigma));
880     fTPCSignal2SigmaReK=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal2SigmaReK));
881     fTPCSignal2SigmaRep=dynamic_cast<TH2F*>(fOutput->FindObject(fTPCSignal2SigmaRep));
882
883   if(fFillNtuple){
884     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
885   }
886
887   //TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
888   //c1->cd();
889   //TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
890   //hMassPt->Draw();
891  
892  return;
893 }
894
895 //________________________________________________________________________
896 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC){
897
898   Int_t lambdacLab[3]={0,0,0};
899   Int_t pdgs[3]={0,0,0};
900   for(Int_t i=0;i<3;i++){
901    AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
902    Int_t lab=daugh->GetLabel();
903    if(lab<0) return 0;
904    AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
905    if(!part) continue;
906    pdgs[i]=part->GetPdgCode();
907    Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
908    if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
909         Int_t motherLabel=part->GetMother();
910              if(motherLabel<0) return 0;
911    AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
912    if(!motherPart) continue;
913         Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
914    if(motherPdg==4122) {
915         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
916               }
917    if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
918          Int_t GmotherLabel=motherPart->GetMother();
919                if(GmotherLabel<0) return 0;
920          AliAODMCParticle *GmotherPart = (AliAODMCParticle*)arrayMC->At(GmotherLabel);
921          if(!GmotherPart) continue;
922                Int_t GmotherPdg = TMath::Abs(GmotherPart->GetPdgCode());
923          if(GmotherPdg==4122) {
924                 if(GetLambdacDaugh(GmotherPart,arrayMC)) {lambdacLab[i]=GmotherLabel;continue;}
925               }
926         }
927      }
928   }
929
930  if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
931  return 0;
932
933 }
934 //------------------------
935 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC){
936
937  Int_t numberOfLambdac=0;
938  if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
939  Int_t daugh_tmp[2];
940  daugh_tmp[0]=part->GetDaughter(0);
941  daugh_tmp[1]=part->GetDaughter(1);
942  Int_t nDaugh = (Int_t)part->GetNDaughters();
943  if(nDaugh<2) return kFALSE;
944  if(nDaugh>3) return kFALSE;
945  AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
946    if(!pdaugh1) {return kFALSE;}
947  Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
948  AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
949    if(!pdaugh2) {return kFALSE;}
950   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
951
952   if(nDaugh==3){
953    Int_t thirdDaugh=part->GetDaughter(1)-1;
954             AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
955   Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
956    if((number1==321 && number2==211 && number3==2212) || (number1==211 && number2==321 && number3==2212) || (number1==211 && number2==2212 && number3==321) || (number1==321 && number2==2212 && number3==211) || (number1==2212 && number2==321 && number3==211) || (number1==2212 && number2==211 && number3==321)) {
957    numberOfLambdac++;
958    } 
959   }
960
961  if(nDaugh==2){
962
963          //Lambda in Lambda pi
964 //  if((number1==211 && number2==3122)|| (number1==3122 && number2==211))return kFALSE;
965   
966   //Lambda resonant
967   
968   //Lambda -> p K*0
969   //
970   Int_t nfiglieK=0;
971
972    if((number1==2212 && number2==313)){
973      nfiglieK=pdaugh2->GetNDaughters();
974      if(nfiglieK!=2) return kFALSE;
975      AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
976      AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
977      if(!pdaughK1) return kFALSE;
978      if(!pdaughK2) return kFALSE;
979      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
980     }
981
982    if((number1==313 && number2==2212)){
983     nfiglieK=pdaugh1->GetNDaughters();
984     if(nfiglieK!=2) return kFALSE;
985     AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
986     AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
987     if(!pdaughK1) return kFALSE;
988     if(!pdaughK2) return kFALSE;
989      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
990    }
991
992    //Lambda -> Delta++ k
993    Int_t nfiglieDelta=0;
994    if(number1==321 && number2==2224){
995     nfiglieDelta=pdaugh2->GetNDaughters();
996     if(nfiglieDelta!=2) return kFALSE;
997     AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
998     AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
999     if(!pdaughD1) return kFALSE;
1000     if(!pdaughD2) return kFALSE;
1001     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1002    }
1003    if(number1==2224 && number2==321){
1004     nfiglieDelta=pdaugh1->GetNDaughters();
1005     if(nfiglieDelta!=2) return kFALSE;
1006     AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1007     AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1008     if(!pdaughD1) return kFALSE;
1009     if(!pdaughD2) return kFALSE;
1010     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1011    }
1012     
1013
1014    //Lambdac -> Lambda(1520) pi
1015    Int_t nfiglieLa=0;
1016    if(number1==3124 && number2==211){
1017     nfiglieLa=pdaugh1->GetNDaughters();
1018     if(nfiglieLa!=2) return kFALSE;
1019     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1020     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1021     if(!pdaughL1) return kFALSE;
1022     if(!pdaughL2) return kFALSE;
1023     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1024    }
1025    if(number1==211 && number2==3124){
1026     nfiglieLa=pdaugh2->GetNDaughters();
1027     if(nfiglieLa!=2) return kFALSE;
1028     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1029     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1030     if(!pdaughL1) return kFALSE;
1031     if(!pdaughL2) return kFALSE;
1032     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1033    
1034    }
1035  }
1036
1037  if(numberOfLambdac>0) {return kTRUE;}
1038          return kFALSE;
1039 }
1040 //-----------------------------
1041 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC){
1042
1043    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1044    for(Int_t i=0;i<3;i++){
1045     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1046     lab[i]=daugh->GetLabel();
1047     if(lab[i]<0) return kFALSE;
1048     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1049     if(!part) return kFALSE;
1050     pdgs[i]=TMath::Abs(part->GetPdgCode());
1051    }
1052
1053    if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1054
1055    return kFALSE;
1056 }
1057 //-----------------------------
1058 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC){
1059
1060    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1061    for(Int_t i=0;i<3;i++){
1062     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1063     lab[i]=daugh->GetLabel();
1064     if(lab[i]<0) return kFALSE;
1065     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1066     if(!part) return kFALSE;
1067     pdgs[i]=TMath::Abs(part->GetPdgCode());
1068    }
1069
1070    if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1071
1072    return kFALSE;
1073 }
1074 //--------------------------------------
1075 Bool_t AliAnalysisTaskSELambdac::IspiKpReal(AliAODRecoDecayHF3Prong *d){
1076
1077   AliAODPidHF* pidObj=new AliAODPidHF();
1078   Bool_t type[2]={kFALSE,kFALSE};
1079   for(Int_t i=0;i<3;i++){
1080 //   Bool_t pid[3]={kFALSE,kFALSE,kFALSE};
1081    AliAODTrack *track=(AliAODTrack*)d->GetDaughter(i);
1082    AliAODPid *pidObjtrk=track->GetDetPid();
1083 //   pidObj->CombinedProbability(track,pid);
1084
1085 //   Int_t i2=i-1;
1086 //   type[i2]=pid[i];
1087    Double_t dedx=pidObjtrk->GetTPCsignal();
1088    Double_t mom = pidObjtrk->GetTPCmomentum();
1089    //printf("TPC mom= %f\n",mom);
1090     pidObj->SetSigma(3.);
1091    if(pidObj->IsProtonRaw(track,"TPC")){ 
1092     fTPCSignal3Sigma->Fill(mom,dedx);
1093  //   if(i==2) type[1]=kTRUE;
1094     pidObj->SetSigma(1.);
1095     //if(!pidObj->IsPionRaw(track,"TPC") && !pidObj->IsKaonRaw(track,"TPC"))
1096     if(!pidObj->IsKaonRaw(track,"TPC")){
1097      fTPCSignal3SigmaRep->Fill(mom,dedx);
1098     }
1099    }
1100     pidObj->SetSigma(3.);
1101    if(pidObj->IsKaonRaw(track,"TPC")) {
1102     fTPCSignal3Sigma->Fill(mom,dedx);
1103  //   if(i==1) type[0]=kTRUE;
1104     pidObj->SetSigma(1.);
1105     //if(!pidObj->IsPionRaw(track,"TPC") && !pidObj->IsProtonRaw(track,"TPC"))
1106     if(!pidObj->IsPionRaw(track,"TPC")){
1107      fTPCSignal3SigmaReK->Fill(mom,dedx);
1108     }
1109    }
1110
1111     pidObj->SetSigma(2.);
1112    if(pidObj->IsProtonRaw(track,"TPC")){ 
1113     if(i==2) type[1]=kTRUE;
1114     fTPCSignal2Sigma->Fill(mom,dedx);
1115     pidObj->SetSigma(1.);
1116     //if(!pidObj->IsPionRaw(track,"TPC") && !pidObj->IsKaonRaw(track,"TPC"))
1117     if(!pidObj->IsKaonRaw(track,"TPC")){
1118      fTPCSignal2SigmaRep->Fill(mom,dedx);
1119     }
1120    }
1121    if(pidObj->IsKaonRaw(track,"TPC")) {
1122     if(i==1) type[0]=kTRUE;
1123     fTPCSignal2Sigma->Fill(mom,dedx);
1124     pidObj->SetSigma(1.);
1125     //if(!pidObj->IsPionRaw(track,"TPC") && !pidObj->IsProtonRaw(track,"TPC"))
1126     if(!pidObj->IsPionRaw(track,"TPC")){
1127      fTPCSignal2SigmaReK->Fill(mom,dedx);
1128     }
1129    }
1130   }
1131
1132   if(type[0] && type[1]) {delete pidObj;return kTRUE;}
1133   delete pidObj;
1134  return kFALSE;
1135 }
1136 //--------------------------------------
1137 Bool_t AliAnalysisTaskSELambdac::IspKpiReal(AliAODRecoDecayHF3Prong *d){
1138
1139   AliAODPidHF* pidObj=new AliAODPidHF();
1140   Bool_t type[2]={kFALSE,kFALSE};
1141     pidObj->SetSigma(2.);
1142   for(Int_t i=0;i<2;i++){
1143    //Bool_t pid[3]={kFALSE,kFALSE,kFALSE};
1144    AliAODTrack *track=(AliAODTrack*)d->GetDaughter(i);
1145   // pidObj->CombinedProbability(track,pid);
1146   // if(i==0) type[i]=pid[2];
1147   // if(i==1) type[i]=pid[i];
1148    if(pidObj->IsProtonRaw(track,"TPC") && i==0) type[1]=kTRUE;
1149    if(pidObj->IsKaonRaw(track,"TPC") && i==1) type[0]=kTRUE;
1150   }
1151   if(type[0] && type[1]) {delete pidObj;return kTRUE;}
1152   delete pidObj;
1153  return kFALSE;
1154 }
1155 //------------------------------------
1156 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field){
1157   
1158    Int_t iprongs[3]={0,1,2};
1159    Double_t mass[2]={0.,0.};
1160  //topological constr
1161    AliKFParticle *Lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,kTRUE,field,mass);
1162    if(!Lambdac) return kFALSE;
1163   Double_t probTot=TMath::Prob(Lambdac->GetChi2(),Lambdac->GetNDF());
1164   if(probTot<fCutsKF[0]) return kFALSE;
1165
1166   if(d->PtProng(0)<1.5 && d->Getd0Prong(0)>0.02) return kFALSE;
1167     if(d->PtProng(1)<1.5 && d->Getd0Prong(1)>0.02) return kFALSE;
1168       if(d->PtProng(2)<1.5 && d->Getd0Prong(2)>0.02) return kFALSE;
1169  //mass constr for K*
1170    Int_t ipRes[2];
1171    Int_t pdgres[2];
1172    mass[0]=0.8961;mass[1]=0.03;
1173    if(TMath::Abs(pdgs[0])==211){
1174     ipRes[0]=0;ipRes[1]=1;
1175     pdgres[0]=pdgs[0];pdgres[1]=321;
1176    }
1177    if(TMath::Abs(pdgs[2])==211){
1178     ipRes[0]=2;ipRes[1]=1;
1179     pdgres[0]=pdgs[2];pdgres[1]=321;
1180    }
1181    AliKFParticle *KappaStar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1182
1183   Double_t probKstar=TMath::Prob(KappaStar->GetChi2(),KappaStar->GetNDF());
1184   if(probKstar>fCutsKF[1]) {
1185     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1186     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1187     AliKFParticle prong1(*esdProng1,pdgres[0]);
1188     AliKFParticle prong2(*esdProng2,pdgres[1]);
1189     if(KappaStar->GetPt()<fCutsKF[2] && prong1.GetAngle(prong2)>fCutsKF[3]) return kFALSE;
1190   } 
1191  //mass constr for Lambda
1192    mass[0]=1.520;mass[1]=0.005;
1193    if(TMath::Abs(pdgs[0])==2212){
1194     ipRes[0]=0;ipRes[1]=1;
1195     pdgres[0]=pdgs[0];pdgres[1]=pdgs[1];
1196    }
1197    if(TMath::Abs(pdgs[2])==2212){
1198     ipRes[0]=2;ipRes[1]=1;
1199     pdgres[0]=pdgs[2];pdgres[1]=pdgs[1];
1200    }
1201    AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1202   Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1203   if(probLa>fCutsKF[4]) {
1204     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1205     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1206     AliKFParticle prong1(*esdProng1,pdgres[0]);
1207     AliKFParticle prong2(*esdProng2,pdgres[1]);
1208     if(Lambda1520->GetPt()<fCutsKF[5] && prong1.GetAngle(prong2)>fCutsKF[6]) return kFALSE;
1209   } 
1210  //mass constr for Delta
1211    mass[0]=1.2;mass[1]=0.15;
1212    ipRes[0]=0;ipRes[1]=2;
1213    pdgres[0]=pdgs[0];pdgres[2]=pdgs[2];
1214    AliKFParticle *Delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1215   Double_t probDelta=TMath::Prob(Delta->GetChi2(),Delta->GetNDF());
1216   if(probDelta>fCutsKF[7]) {
1217     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1218     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1219     AliKFParticle prong1(*esdProng1,pdgres[0]);
1220     AliKFParticle prong2(*esdProng2,pdgres[1]);
1221     if(Delta->GetPt()<fCutsKF[8] && prong1.GetAngle(prong2)>fCutsKF[9]) return kFALSE;
1222   } 
1223  return kTRUE;
1224 }
1225 //-------------------------------------
1226 Bool_t AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field){
1227   
1228  //if lambda* -> pk
1229         Double_t mass[2]={1.520,0.005};
1230         Int_t ipRes[2]={1,2};
1231         Int_t pdgres[2]={321,2212};
1232         AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1233         Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1234         if(probLa>0.1) return kTRUE;
1235  //K* -> kpi
1236 //        mass[0]=0.8961;mass[1]=0.03;
1237 //        ipRes[0]=0;ipRes[1]=1;
1238 //        pdgres[0]=211;pdgres[1]=321;
1239 //        AliKFParticle *Kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1240 //      Double_t probKa=TMath::Prob(Kstar->GetChi2(),Kstar->GetNDF());
1241 //      if(probKa>0.1) return kTRUE;
1242
1243  return kFALSE;
1244
1245 }
1246 //-------------------------------------
1247 Bool_t AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field){
1248   
1249  //if lambda* -> pk
1250         Double_t mass[2]={1.520,0.005};
1251         Int_t ipRes[2]={0,1};
1252         Int_t pdgres[2]={2212,321};
1253         AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1254         Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1255         if(probLa>0.1) return kTRUE;
1256  //K* -> kpi
1257 //        mass[0]=0.8961;mass[1]=0.03;
1258 //        ipRes[0]=1;ipRes[1]=2;
1259 //        pdgres[1]=211;pdgres[0]=321;
1260 //        AliKFParticle *Kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1261 //      Double_t probKa=TMath::Prob(Kstar->GetChi2(),Kstar->GetNDF());
1262 //      if(probKa>0.1) return kTRUE;
1263
1264  return kFALSE;
1265
1266 }