]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSELambdac.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New task for Lc->pKpi (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//
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
50ClassImp(AliAnalysisTaskSELambdac)
51
52
53//________________________________________________________________________
54AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac():
55AliAnalysisTaskSE(),
56fOutput(0),
57fHistNEvents(0),
58fTPCSignal3Sigma(0),
59fTPCSignal3SigmaReK(0),
60fTPCSignal3SigmaRep(0),
61fTPCSignal2Sigma(0),
62fTPCSignal2SigmaReK(0),
63fTPCSignal2SigmaRep(0),
64fNtupleLambdac(0),
65fUpmasslimit(2.486),
66fLowmasslimit(2.086),
67fNPtBins(0),
68fRDCutsAnalysis(0),
69fRDCutsProduction(0),
70fListCuts(0),
71fFillNtuple(kFALSE),
72fReadMC(kFALSE),
73fMCPid(kFALSE),
74fRealPid(kFALSE),
75fResPid(kTRUE),
76fUseKF(kFALSE),
77fVHF(0)
78{
79 // Default constructor
80}
81
82//________________________________________________________________________
83AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac(const char *name,Bool_t fillNtuple,AliRDHFCutsLctopKpi *lccutsana,AliRDHFCutsLctopKpi *lccutsprod):
84AliAnalysisTaskSE(name),
85fOutput(0),
86fHistNEvents(0),
87fTPCSignal3Sigma(0),
88fTPCSignal3SigmaReK(0),
89fTPCSignal3SigmaRep(0),
90fTPCSignal2Sigma(0),
91fTPCSignal2SigmaReK(0),
92fTPCSignal2SigmaRep(0),
93fNtupleLambdac(0),
94fUpmasslimit(2.486),
95fLowmasslimit(2.086),
96fNPtBins(0),
97fRDCutsAnalysis(lccutsana),
98fRDCutsProduction(lccutsprod),
99fListCuts(0),
100fFillNtuple(fillNtuple),
101fReadMC(kFALSE),
102fMCPid(kFALSE),
103fRealPid(kTRUE),
104fResPid(kFALSE),
105fUseKF(kFALSE),
106fVHF(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//________________________________________________________________________
123AliAnalysisTaskSELambdac::~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//_________________________________________________________________
150void AliAnalysisTaskSELambdac::SetMassLimits(Float_t range){
151 fUpmasslimit = 2.286+range;
152 fLowmasslimit = 2.286-range;
153}
154//_________________________________________________________________
155void AliAnalysisTaskSELambdac::SetMassLimits(Float_t lowlimit, Float_t uplimit){
156 if(uplimit>lowlimit)
157 {
158 fUpmasslimit = lowlimit;
159 fLowmasslimit = uplimit;
160 }
161}
162
163
164//________________________________________________________________________
165void 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//_________________________________________________________________
193Double_t AliAnalysisTaskSELambdac::GetPtBinLimit(Int_t ibin){
194 if(ibin>fNPtBins)return -1;
195 return fArrayBinLimits[ibin];
196}
197
198//_________________________________________________________________
199void 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//________________________________________________________________________
219void 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//________________________________________________________________________
481void 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();
703Double_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//________________________________________________________________________
801void 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//________________________________________________________________________
896Int_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//------------------------
935Bool_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//-----------------------------
1041Bool_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//-----------------------------
1058Bool_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//--------------------------------------
1075Bool_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//--------------------------------------
1137Bool_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//------------------------------------
1156Bool_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//-------------------------------------
1226Bool_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//-------------------------------------
1247Bool_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}