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