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