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