]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliVertexingHFUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Methods to check Lc decays in the Kinematics
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
CommitLineData
a6c5a2e9 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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#include <TMath.h>
686e4710 17#include <TRandom.h>
b176db3b 18#include <TProfile.h>
19#include <TClonesArray.h>
20#include <TH1F.h>
21#include <TH2F.h>
22#include <TF1.h>
ee9c9d9b 23#include <TParticle.h>
24#include "AliStack.h"
a6c5a2e9 25#include "AliAODEvent.h"
b176db3b 26#include "AliAODMCHeader.h"
7f51a627 27#include "AliGenEventHeader.h"
25504150 28#include "AliAODMCParticle.h"
b176db3b 29#include "AliAODRecoDecayHF.h"
a6c5a2e9 30#include "AliVertexingHFUtils.h"
31
d947ea9c 32/* $Id$ */
a6c5a2e9 33
34///////////////////////////////////////////////////////////////////
35// //
36// Class with functions useful for different D2H analyses //
37// - event plane resolution //
38// - <pt> calculation with side band subtraction //
25504150 39// - tracklet multiplicity calculation //
a6c5a2e9 40// Origin: F.Prino, Torino, prino@to.infn.it //
41// //
42///////////////////////////////////////////////////////////////////
43
44ClassImp(AliVertexingHFUtils)
45
46//______________________________________________________________________
47AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
48 fK(1),
49 fSubRes(1.),
50 fMinEtaForTracklets(-1.),
51 fMaxEtaForTracklets(1.)
52{
53 // Default contructor
54}
55
56
57//______________________________________________________________________
58AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
59 TObject(),
60 fK(k),
61 fSubRes(1.),
62 fMinEtaForTracklets(-1.),
63 fMaxEtaForTracklets(1.)
64{
65 // Standard constructor
66}
67
68
5d1396f9 69//______________________________________________________________________
70void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
71 // calculate significance from S, B and errors
72
73
74 Double_t errSigSq=errsignal*errsignal;
75 Double_t errBkgSq=errbackground*errbackground;
76 Double_t sigPlusBkg=signal+background;
77 if(sigPlusBkg>0. && signal>0.){
78 significance = signal/TMath::Sqrt(signal+background);
79 errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
80 }else{
81 significance=0.;
82 errsignificance=0.;
83 }
84 return;
85
86}
a6c5a2e9 87//______________________________________________________________________
88Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
89 // compute chi from polynomial approximation
90 Double_t c[5];
91 if(k==1){
92 c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
93 }
94 else if(k==2){
95 c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
96 }
97 else return -1;
98 return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
99}
100
101//______________________________________________________________________
102Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
103 return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
104}
105
106
107//______________________________________________________________________
108Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
109 // compute chi variable (=v2*sqrt(N)) from external values
110
111 Double_t x1=0;
112 Double_t x2=15;
113 Double_t y1,y2;
114 if(k==1){
115 y1=ResolK1(x1)-res;
116 y2=ResolK1(x2)-res;
117 }
118 else if(k==2){
119 y1=Pol(x1,2)-res;
120 y2=Pol(x2,2)-res;
121 }
122 else return -1;
123
124 if(y1*y2>0) return -1;
125 if(y1==0) return y1;
126 if(y2==0) return y2;
127 Double_t xmed,ymed;
128 Int_t jiter=0;
129 while((x2-x1)>0.0001){
130 xmed=0.5*(x1+x2);
131 if(k==1){
132 y1=ResolK1(x1)-res;
133 y2=ResolK1(x2)-res;
134 ymed=ResolK1(xmed)-res;
135 }
136 else if(k==2){
137 y1=Pol(x1,2)-res;
138 y2=Pol(x2,2)-res;
139 ymed=Pol(xmed,2)-res;
140 }
141 else return -1;
142 if((y1*ymed)<0) x2=xmed;
143 if((y2*ymed)<0)x1=xmed;
144 if(ymed==0) return xmed;
145 jiter++;
146 }
147 return 0.5*(x1+x2);
148}
149
150//______________________________________________________________________
151Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
152 // computes event plane resolution starting from sub event resolution
153 Double_t chisub=FindChi(resSub,k);
154 Double_t chifull=chisub*TMath::Sqrt(2);
155 if(k==1) return ResolK1(chifull);
156 else if(k==2) return Pol(chifull,2);
157 else{
158 printf("k should be <=2\n");
159 return 1.;
160 }
161}
162
163//______________________________________________________________________
164Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
165 // computes event plane resolution starting from sub event correlation histogram
166 if(!hSubEvCorr) return 1.;
167 Double_t resSub=GetSubEvResol(hSubEvCorr);
168 return GetFullEvResol(resSub,k);
169}
170//______________________________________________________________________
171Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
172 // computes low limit event plane resolution starting from sub event correlation histogram
173 if(!hSubEvCorr) return 1.;
174 Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
175 printf("%f\n",resSub);
176 return GetFullEvResol(resSub,k);
177}
178//______________________________________________________________________
179Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
180 // computes high limit event plane resolution starting from sub event correlation histogram
181 if(!hSubEvCorr) return 1.;
182 Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
183 printf("%f\n",resSub);
184 return GetFullEvResol(resSub,k);
185}
186//______________________________________________________________________
187Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
188 // counts tracklets in given eta range
189 AliAODTracklets* tracklets=ev->GetTracklets();
190 Int_t nTr=tracklets->GetNumberOfTracklets();
191 Int_t count=0;
192 for(Int_t iTr=0; iTr<nTr; iTr++){
193 Double_t theta=tracklets->GetTheta(iTr);
194 Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
195 if(eta>mineta && eta<maxeta) count++;
196 }
197 return count;
198}
86d84fa4 199//______________________________________________________________________
25504150 200Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
201 // counts generated particles in fgiven eta range
202
203 Int_t nChargedMC=0;
204 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
205 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
206 Int_t charge = part->Charge();
207 Double_t eta = part->Eta();
208 if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
209 }
210 return nChargedMC;
211}
212//______________________________________________________________________
213Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
214 // counts generated primary particles in given eta range
215
216 Int_t nChargedMC=0;
217 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
218 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
219 Int_t charge = part->Charge();
220 Double_t eta = part->Eta();
221 if(charge!=0 && eta>mineta && eta<maxeta){
222 if(part->IsPrimary())nChargedMC++;
223 }
224 }
225 return nChargedMC;
226}
227//______________________________________________________________________
228Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
229 // counts generated primary particles in given eta range
230
231 Int_t nChargedMC=0;
232 for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
233 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
234 Int_t charge = part->Charge();
235 Double_t eta = part->Eta();
236 if(charge!=0 && eta>mineta && eta<maxeta){
237 if(part->IsPhysicalPrimary())nChargedMC++;
238 }
239 }
240 return nChargedMC;
241}
b242bddf 242
243
244//______________________________________________________________________
245Double_t AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(AliAODEvent* ev){
246 //
247 // Method to get VZERO-A equalized multiplicty as done in AliCentralitySelectionTask
248 // getting the equalized VZERO factors from the tender or AOD
249 // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
250
251 Double_t multV0AEq=0;
252 for(Int_t iCh = 32; iCh < 64; ++iCh) {
253 Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
254 multV0AEq += mult;
255 }
256 return multV0AEq;
257}
258
259//______________________________________________________________________
260Double_t AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(AliAODEvent* ev){
261 // Method to get VZERO-C equalized multiplicty as done in AliCentralitySelectionTask
262 // getting the equalized VZERO factors from the tender or AOD
263 // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
264
265 Double_t multV0CEq=0;
266 for(Int_t iCh = 0; iCh < 32; ++iCh) {
267 Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
268 multV0CEq += mult;
269 }
270 return multV0CEq;
271}
272
25504150 273//______________________________________________________________________
19df7ca5 274void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Float_t minMass, Float_t maxMass, Int_t rebin){
86d84fa4 275
276 // Compute <pt> from 2D histogram M vs pt
277
278 //Make 2D histos in the desired pt range
279 Int_t start=hMassD->FindBin(ptmin);
280 Int_t end=hMassD->FindBin(ptmax)-1;
281 const Int_t nx=end-start;
282 TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
283 for(Int_t ix=start;ix<end;ix++){
284 for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
285 hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
286 hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
287 }
288 }
289
290 Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
291 Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
292 Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
293 Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
294 Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
295 Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
296 // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
297
298 Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
19df7ca5 299 Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
86d84fa4 300 Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
301 Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
302 Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
303 Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
304 Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
19df7ca5 305 Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
86d84fa4 306 Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
307 Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
308 // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
309
310 Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
311 Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
312 Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
313 // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
314
315 TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
316 TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
317 TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
318
319 hMptBkgLo->Rebin(rebin);
320 hMptBkgHi->Rebin(rebin);
321 hMptSigReg->Rebin(rebin);
322
323 hMptBkgLo->Sumw2();
324 hMptBkgHi->Sumw2();
325 TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
326 hMptBkgLoScal->Scale(bkgSig/bkgLow);
327 TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
328 hMptBkgHiScal->Scale(bkgSig/bkgHi);
329
330 TH1F* hMptBkgAver=0x0;
331 hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
332 hMptBkgAver->Add(hMptBkgHiScal);
333 hMptBkgAver->Scale(0.5);
334 TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
335 hMptSig->Add(hMptBkgAver,-1.);
336
337 averagePt = hMptSig->GetMean();
338 errorPt = hMptSig->GetMeanError();
339
340 delete hMptBkgLo;
341 delete hMptBkgHi;
342 delete hMptSigReg;
343 delete hMptBkgLoScal;
344 delete hMptBkgHiScal;
345 delete hMptBkgAver;
346 delete hMassDpt;
347 delete hMptSig;
348
15917f37 349}
350//____________________________________________________________________________
e6eaae4c 351Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
352 // check if T0VTX trigger was fired, based on a workaround suggested by Alla
353 const Double32_t *mean = aodEv->GetT0TOF();
94d7319f 354 if(mean && mean[0]<9999.) return kTRUE;
e6eaae4c 355 else return kFALSE;
356}
357//____________________________________________________________________________
15917f37 358Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
359 // true impact parameter calculation for Dzero
360
361 if(!partD || !arrayMC || !mcHeader) return 99999.;
362 Int_t code=partD->GetPdgCode();
363 if(TMath::Abs(code)!=421) return 99999.;
364
365 Double_t vtxTrue[3];
366 mcHeader->GetVertex(vtxTrue);
367 Double_t origD[3];
368 partD->XvYvZv(origD);
369 Short_t charge=partD->Charge();
370 Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
371 for(Int_t iDau=0; iDau<2; iDau++){
372 pXdauTrue[iDau]=0.;
373 pYdauTrue[iDau]=0.;
374 pZdauTrue[iDau]=0.;
375 }
376
377 Int_t nDau=partD->GetNDaughters();
378 Int_t labelFirstDau = partD->GetDaughter(0);
379 if(nDau==2){
380 for(Int_t iDau=0; iDau<2; iDau++){
381 Int_t ind = labelFirstDau+iDau;
382 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
383 if(!part){
384 printf("Daughter particle not found in MC array");
385 return 99999.;
386 }
387 pXdauTrue[iDau]=part->Px();
388 pYdauTrue[iDau]=part->Py();
389 pZdauTrue[iDau]=part->Pz();
390 }
391 }else{
392 return 99999.;
393 }
394
48efc6ad 395 Double_t d0dummy[2]={0.,0.};
15917f37 396 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
397 return aodDvsMC.ImpParXY();
398
86d84fa4 399}
b176db3b 400//____________________________________________________________________________
401Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
402 // true impact parameter calculation for Dplus
403
404 if(!partD || !arrayMC || !mcHeader) return 99999.;
405 Int_t code=partD->GetPdgCode();
406 if(TMath::Abs(code)!=411) return 99999.;
407
408 Double_t vtxTrue[3];
409 mcHeader->GetVertex(vtxTrue);
410 Double_t origD[3];
411 partD->XvYvZv(origD);
412 Short_t charge=partD->Charge();
413 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
414 for(Int_t iDau=0; iDau<3; iDau++){
415 pXdauTrue[iDau]=0.;
416 pYdauTrue[iDau]=0.;
417 pZdauTrue[iDau]=0.;
418 }
419
420 Int_t nDau=partD->GetNDaughters();
421 Int_t labelFirstDau = partD->GetDaughter(0);
422 if(nDau==3){
423 for(Int_t iDau=0; iDau<3; iDau++){
424 Int_t ind = labelFirstDau+iDau;
425 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
426 if(!part){
427 printf("Daughter particle not found in MC array");
428 return 99999.;
429 }
430 pXdauTrue[iDau]=part->Px();
431 pYdauTrue[iDau]=part->Py();
432 pZdauTrue[iDau]=part->Pz();
433 }
434 }else if(nDau==2){
435 Int_t theDau=0;
436 for(Int_t iDau=0; iDau<2; iDau++){
437 Int_t ind = labelFirstDau+iDau;
438 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
439 if(!part){
440 printf("Daughter particle not found in MC array");
441 return 99999.;
442 }
443 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
444 if(pdgCode==211 || pdgCode==321){
445 pXdauTrue[theDau]=part->Px();
446 pYdauTrue[theDau]=part->Py();
447 pZdauTrue[theDau]=part->Pz();
448 ++theDau;
449 }else{
450 Int_t nDauRes=part->GetNDaughters();
451 if(nDauRes==2){
452 Int_t labelFirstDauRes = part->GetDaughter(0);
453 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
454 Int_t indDR = labelFirstDauRes+iDauRes;
455 AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
456 if(!partDR){
457 printf("Daughter particle not found in MC array");
458 return 99999.;
459 }
460
461 Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
462 if(pdgCodeDR==211 || pdgCodeDR==321){
463 pXdauTrue[theDau]=partDR->Px();
464 pYdauTrue[theDau]=partDR->Py();
465 pZdauTrue[theDau]=partDR->Pz();
466 ++theDau;
467 }
468 }
469 }
470 }
471 }
472 if(theDau!=3){
473 printf("Wrong number of decay prongs");
474 return 99999.;
475 }
476 }
477
478 Double_t d0dummy[3]={0.,0.,0.};
479 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
480 return aodDvsMC.ImpParXY();
481
482}
483
484
485
686e4710 486//____________________________________________________________________________
487Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
488 //
489 // Correct the number of accepted tracklets based on a TProfile Hist
490 //
491 //
492
493 if(TMath::Abs(vtxZ)>10.0){
280cc141 494 // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
686e4710 495 return uncorrectedNacc;
496 }
497
498 if(!estimatorAvg){
499 printf("ERROR: Missing TProfile for correction of multiplicity\n");
500 return uncorrectedNacc;
501 }
502
503 Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
504
505 Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
506
507 Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
508
509 return correctedNacc;
510}
7f51a627 511//______________________________________________________________________
512TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
f98f4c22 513 // get the name of the generator that produced a given particle
514
515 Int_t nsumpart=0;
516 TList *lh=header->GetCocktailHeaders();
517 Int_t nh=lh->GetEntries();
518 for(Int_t i=0;i<nh;i++){
519 AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
520 TString genname=gh->GetName();
521 Int_t npart=gh->NProduced();
522 if(label>=nsumpart && label<(nsumpart+npart)) return genname;
523 nsumpart+=npart;
524 }
525 TString empty="";
526 return empty;
527}
7c0e9154 528//_____________________________________________________________________
529void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
530
531 // method to check if a track comes from a given generator
f98f4c22 532
533 Int_t lab=TMath::Abs(track->GetLabel());
7c0e9154 534 nameGen=GetGenerator(lab,header);
f98f4c22 535
536 // Int_t countControl=0;
537
538 while(nameGen.IsWhitespace()){
539 AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
540 if(!mcpart){
541 printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
542 break;
543 }
544 Int_t mother = mcpart->GetMother();
545 if(mother<0){
546 printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
547 break;
548 }
549 lab=mother;
550 nameGen=GetGenerator(mother,header);
551 // countControl++;
552 // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
553 // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
554 // break;
555 // }
556 }
557
7c0e9154 558 return;
559}
560//----------------------------------------------------------------------
561Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
562 // method to check if a track comes from the signal event or from the underlying Hijing event
563 TString nameGen;
564 GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
565
f98f4c22 566 if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
567
568 return kTRUE;
569}
570//____________________________________________________________________________
571Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
572 // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
573
574 Int_t nprongs=cand->GetNProngs();
575 for(Int_t i=0;i<nprongs;i++){
576 AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
577 if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
578 }
579 return kFALSE;
580}
c8781624 581//____________________________________________________________________________
d49d50fc 582Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
583 // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event
584
585 AliAODTrack* bach = cand->GetBachelor();
586 if(IsTrackInjected(bach, header, arrayMC)) {
587 AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
588 return kTRUE;
589 }
590 AliAODv0* v0 = cand->Getv0();
591 Int_t nprongs = v0->GetNProngs();
592 for(Int_t i = 0; i < nprongs; i++){
593 AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
594 if(IsTrackInjected(daugh,header,arrayMC)) {
595 AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
596 return kTRUE;
597 }
598 }
599 return kFALSE;
600}
601//____________________________________________________________________________
ee9c9d9b 602Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
603 // checking whether the mother of the particles come from a charm or a bottom quark
604
605 Int_t pdgGranma = 0;
606 Int_t mother = 0;
607 mother = mcPart->GetFirstMother();
608 Int_t istep = 0;
609 Int_t abspdgGranma =0;
610 Bool_t isFromB=kFALSE;
611 Bool_t isQuarkFound=kFALSE;
612 while (mother >0 ){
613 istep++;
614 TParticle* mcGranma = stack->Particle(mother);
615 if (mcGranma){
616 pdgGranma = mcGranma->GetPdgCode();
617 abspdgGranma = TMath::Abs(pdgGranma);
618 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
619 isFromB=kTRUE;
620 }
621 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
622 mother = mcGranma->GetFirstMother();
623 }else{
624 printf("CheckOrigin: Failed casting the mother particle!");
625 break;
626 }
627 }
628 if(searchUpToQuark && !isQuarkFound) return 0;
629 if(isFromB) return 5;
630 else return 4;
631
632}
633//____________________________________________________________________________
c8781624 634Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
635 // checking whether the mother of the particles come from a charm or a bottom quark
636
637 Int_t pdgGranma = 0;
638 Int_t mother = 0;
639 mother = mcPart->GetMother();
640 Int_t istep = 0;
641 Int_t abspdgGranma =0;
642 Bool_t isFromB=kFALSE;
643 Bool_t isQuarkFound=kFALSE;
644 while (mother >0 ){
645 istep++;
646 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
647 if (mcGranma){
648 pdgGranma = mcGranma->GetPdgCode();
649 abspdgGranma = TMath::Abs(pdgGranma);
650 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
651 isFromB=kTRUE;
652 }
653 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
654 mother = mcGranma->GetMother();
655 }else{
656 printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
657 break;
658 }
659 }
660 if(searchUpToQuark && !isQuarkFound) return 0;
661 if(isFromB) return 5;
662 else return 4;
663
664}
ee9c9d9b 665
c8781624 666//____________________________________________________________________________
ee9c9d9b 667Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
c8781624 668 // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
669
ee9c9d9b 670 if(label<0) return -1;
671 TParticle* mcPart = stack->Particle(label);
672 Int_t pdgD=mcPart->GetPdgCode();
673 if(TMath::Abs(pdgD)!=421) return -1;
674
675 Int_t nDau=mcPart->GetNDaughters();
676
677 if(nDau==2){
678 Int_t daughter0 = mcPart->GetDaughter(0);
679 Int_t daughter1 = mcPart->GetDaughter(1);
680 TParticle* mcPartDaughter0 = stack->Particle(daughter0);
681 TParticle* mcPartDaughter1 = stack->Particle(daughter1);
682 if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
683 arrayDauLab[0]=daughter0;
684 arrayDauLab[1]=daughter1;
685 Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
686 Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
687 if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
688 !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
689 return -1;
690 }
691 if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
692 if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
693 if((pdgCode0*pdgCode1)>0) return -1;
694 Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
695 Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
696 Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
697 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
698 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
699 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
700 return 1;
701 }
702
703 if(nDau==3 || nDau==4){
704 Int_t nKaons=0;
705 Int_t nPions=0;
706 Double_t sumPxDau=0.;
707 Double_t sumPyDau=0.;
708 Double_t sumPzDau=0.;
709 Int_t nFoundKpi=0;
710 Int_t labelFirstDau = mcPart->GetDaughter(0);
711 for(Int_t iDau=0; iDau<nDau; iDau++){
712 Int_t indDau = labelFirstDau+iDau;
713 if(indDau<0) return -1;
714 TParticle* dau=stack->Particle(indDau);
715 if(!dau) return -1;
716 Int_t pdgdau=dau->GetPdgCode();
717 if(TMath::Abs(pdgdau)==321){
718 if(pdgD>0 && pdgdau>0) return -1;
719 if(pdgD<0 && pdgdau<0) return -1;
720 nKaons++;
721 sumPxDau+=dau->Px();
722 sumPyDau+=dau->Py();
723 sumPzDau+=dau->Pz();
724 arrayDauLab[nFoundKpi++]=indDau;
725 if(nFoundKpi>4) return -1;
726 }else if(TMath::Abs(pdgdau)==211){
727 nPions++;
728 sumPxDau+=dau->Px();
729 sumPyDau+=dau->Py();
730 sumPzDau+=dau->Pz();
731 arrayDauLab[nFoundKpi++]=indDau;
732 if(nFoundKpi>4) return -1;
733 }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
734 Int_t nResDau=dau->GetNDaughters();
735 if(nResDau!=2) return -1;
736 Int_t indFirstResDau=dau->GetDaughter(0);
737 for(Int_t resDau=0; resDau<2; resDau++){
738 Int_t indResDau=indFirstResDau+resDau;
739 if(indResDau<0) return -1;
740 TParticle* resdau=stack->Particle(indResDau);
741 if(!resdau) return -1;
742 Int_t pdgresdau=resdau->GetPdgCode();
743 if(TMath::Abs(pdgresdau)==321){
744 if(pdgD>0 && pdgresdau>0) return -1;
745 if(pdgD<0 && pdgresdau<0) return -1;
746 nKaons++;
747 sumPxDau+=resdau->Px();
748 sumPyDau+=resdau->Py();
749 sumPzDau+=resdau->Pz();
750 arrayDauLab[nFoundKpi++]=indResDau;
751 if(nFoundKpi>4) return -1;
752 }
753 if(TMath::Abs(pdgresdau)==211){
754 nPions++;
755 sumPxDau+=resdau->Px();
756 sumPyDau+=resdau->Py();
757 sumPzDau+=resdau->Pz();
758 arrayDauLab[nFoundKpi++]=indResDau;
759 if(nFoundKpi>4) return -1;
760 }
761 }
762 }else{
763 return -1;
764 }
765 }
766 if(nPions!=3) return -1;
767 if(nKaons!=1) return -1;
768 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
769 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
770 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
771 return 2;
772 }
773 return -1;
774}
775
776//____________________________________________________________________________
777Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
778 // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
c8781624 779
780 Int_t pdgD=mcPart->GetPdgCode();
781 if(TMath::Abs(pdgD)!=421) return -1;
782
783 Int_t nDau=mcPart->GetNDaughters();
784
785 if(nDau==2){
786 Int_t daughter0 = mcPart->GetDaughter(0);
787 Int_t daughter1 = mcPart->GetDaughter(1);
788 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
789 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
790 if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
791 arrayDauLab[0]=daughter0;
792 arrayDauLab[1]=daughter1;
793 Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
794 Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
795 if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
796 !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
797 return -1;
798 }
799 if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
800 if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
801 if((pdgCode0*pdgCode1)>0) return -1;
802 Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
803 Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
804 Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
805 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
806 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
807 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
808 return 1;
809 }
810
811 if(nDau==3 || nDau==4){
812 Int_t nKaons=0;
813 Int_t nPions=0;
814 Double_t sumPxDau=0.;
815 Double_t sumPyDau=0.;
816 Double_t sumPzDau=0.;
817 Int_t nFoundKpi=0;
818 Int_t labelFirstDau = mcPart->GetDaughter(0);
819 for(Int_t iDau=0; iDau<nDau; iDau++){
820 Int_t indDau = labelFirstDau+iDau;
821 if(indDau<0) return -1;
822 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
823 if(!dau) return -1;
824 Int_t pdgdau=dau->GetPdgCode();
825 if(TMath::Abs(pdgdau)==321){
826 if(pdgD>0 && pdgdau>0) return -1;
827 if(pdgD<0 && pdgdau<0) return -1;
828 nKaons++;
829 sumPxDau+=dau->Px();
830 sumPyDau+=dau->Py();
831 sumPzDau+=dau->Pz();
832 arrayDauLab[nFoundKpi++]=indDau;
833 if(nFoundKpi>4) return -1;
834 }else if(TMath::Abs(pdgdau)==211){
835 nPions++;
836 sumPxDau+=dau->Px();
837 sumPyDau+=dau->Py();
838 sumPzDau+=dau->Pz();
839 arrayDauLab[nFoundKpi++]=indDau;
840 if(nFoundKpi>4) return -1;
841 }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
842 Int_t nResDau=dau->GetNDaughters();
843 if(nResDau!=2) return -1;
844 Int_t indFirstResDau=dau->GetDaughter(0);
845 for(Int_t resDau=0; resDau<2; resDau++){
846 Int_t indResDau=indFirstResDau+resDau;
847 if(indResDau<0) return -1;
848 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
849 if(!resdau) return -1;
850 Int_t pdgresdau=resdau->GetPdgCode();
851 if(TMath::Abs(pdgresdau)==321){
852 if(pdgD>0 && pdgresdau>0) return -1;
853 if(pdgD<0 && pdgresdau<0) return -1;
854 nKaons++;
855 sumPxDau+=resdau->Px();
856 sumPyDau+=resdau->Py();
857 sumPzDau+=resdau->Pz();
858 arrayDauLab[nFoundKpi++]=indResDau;
859 if(nFoundKpi>4) return -1;
860 }
861 if(TMath::Abs(pdgresdau)==211){
862 nPions++;
863 sumPxDau+=resdau->Px();
864 sumPyDau+=resdau->Py();
865 sumPzDau+=resdau->Pz();
866 arrayDauLab[nFoundKpi++]=indResDau;
867 if(nFoundKpi>4) return -1;
868 }
869 }
870 }else{
871 return -1;
872 }
873 }
874 if(nPions!=3) return -1;
875 if(nKaons!=1) return -1;
876 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
877 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
878 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
879 return 2;
880 }
881 return -1;
882}
ee9c9d9b 883//____________________________________________________________________________
884Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
885 // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
886
887 if(label<0) return -1;
888 TParticle* mcPart = stack->Particle(label);
889 Int_t pdgD=mcPart->GetPdgCode();
890 if(TMath::Abs(pdgD)!=411) return -1;
891
892 Int_t nDau=mcPart->GetNDaughters();
893 Int_t labelFirstDau = mcPart->GetDaughter(0);
894 Int_t nKaons=0;
895 Int_t nPions=0;
896 Double_t sumPxDau=0.;
897 Double_t sumPyDau=0.;
898 Double_t sumPzDau=0.;
899 Int_t nFoundKpi=0;
900
901 if(nDau==3 || nDau==2){
902 for(Int_t iDau=0; iDau<nDau; iDau++){
903 Int_t indDau = labelFirstDau+iDau;
904 if(indDau<0) return -1;
905 TParticle* dau=stack->Particle(indDau);
906 if(!dau) return -1;
907 Int_t pdgdau=dau->GetPdgCode();
908 if(TMath::Abs(pdgdau)==321){
909 if(pdgD*pdgdau>0) return -1;
910 nKaons++;
911 sumPxDau+=dau->Px();
912 sumPyDau+=dau->Py();
913 sumPzDau+=dau->Pz();
914 arrayDauLab[nFoundKpi++]=indDau;
915 if(nFoundKpi>3) return -1;
916 }else if(TMath::Abs(pdgdau)==211){
917 if(pdgD*pdgdau<0) return -1;
918 nPions++;
919 sumPxDau+=dau->Px();
920 sumPyDau+=dau->Py();
921 sumPzDau+=dau->Pz();
922 arrayDauLab[nFoundKpi++]=indDau;
923 if(nFoundKpi>3) return -1;
924 }else if(TMath::Abs(pdgdau)==313){
925 Int_t nResDau=dau->GetNDaughters();
926 if(nResDau!=2) return -1;
927 Int_t indFirstResDau=dau->GetDaughter(0);
928 for(Int_t resDau=0; resDau<2; resDau++){
929 Int_t indResDau=indFirstResDau+resDau;
930 if(indResDau<0) return -1;
931 TParticle* resdau=stack->Particle(indResDau);
932 if(!resdau) return -1;
933 Int_t pdgresdau=resdau->GetPdgCode();
934 if(TMath::Abs(pdgresdau)==321){
935 if(pdgD*pdgresdau>0) return -1;
936 sumPxDau+=resdau->Px();
937 sumPyDau+=resdau->Py();
938 sumPzDau+=resdau->Pz();
939 nKaons++;
940 arrayDauLab[nFoundKpi++]=indResDau;
941 if(nFoundKpi>3) return -1;
942 }
943 if(TMath::Abs(pdgresdau)==211){
944 if(pdgD*pdgresdau<0) return -1;
945 sumPxDau+=resdau->Px();
946 sumPyDau+=resdau->Py();
947 sumPzDau+=resdau->Pz();
948 nPions++;
949 arrayDauLab[nFoundKpi++]=indResDau;
950 if(nFoundKpi>3) return -1;
951 }
952 }
953 }else{
954 return -1;
955 }
956 }
957 if(nPions!=2) return -1;
958 if(nKaons!=1) return -1;
959 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
960 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
961 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
962 if(nDau==3) return 1;
963 else if(nDau==2) return 2;
964 }
965
966 return -1;
967
968}
969
c8781624 970//____________________________________________________________________________
971Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
972 // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
973
974 Int_t pdgD=mcPart->GetPdgCode();
975 if(TMath::Abs(pdgD)!=411) return -1;
976
977 Int_t nDau=mcPart->GetNDaughters();
978 Int_t labelFirstDau = mcPart->GetDaughter(0);
979 Int_t nKaons=0;
980 Int_t nPions=0;
981 Double_t sumPxDau=0.;
982 Double_t sumPyDau=0.;
983 Double_t sumPzDau=0.;
984 Int_t nFoundKpi=0;
985
986 if(nDau==3 || nDau==2){
987 for(Int_t iDau=0; iDau<nDau; iDau++){
988 Int_t indDau = labelFirstDau+iDau;
989 if(indDau<0) return -1;
990 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
991 if(!dau) return -1;
992 Int_t pdgdau=dau->GetPdgCode();
993 if(TMath::Abs(pdgdau)==321){
994 if(pdgD*pdgdau>0) return -1;
995 nKaons++;
996 sumPxDau+=dau->Px();
997 sumPyDau+=dau->Py();
998 sumPzDau+=dau->Pz();
999 arrayDauLab[nFoundKpi++]=indDau;
1000 if(nFoundKpi>3) return -1;
1001 }else if(TMath::Abs(pdgdau)==211){
1002 if(pdgD*pdgdau<0) return -1;
1003 nPions++;
1004 sumPxDau+=dau->Px();
1005 sumPyDau+=dau->Py();
1006 sumPzDau+=dau->Pz();
1007 arrayDauLab[nFoundKpi++]=indDau;
1008 if(nFoundKpi>3) return -1;
1009 }else if(TMath::Abs(pdgdau)==313){
1010 Int_t nResDau=dau->GetNDaughters();
1011 if(nResDau!=2) return -1;
1012 Int_t indFirstResDau=dau->GetDaughter(0);
1013 for(Int_t resDau=0; resDau<2; resDau++){
1014 Int_t indResDau=indFirstResDau+resDau;
1015 if(indResDau<0) return -1;
1016 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1017 if(!resdau) return -1;
1018 Int_t pdgresdau=resdau->GetPdgCode();
1019 if(TMath::Abs(pdgresdau)==321){
1020 if(pdgD*pdgresdau>0) return -1;
1021 sumPxDau+=resdau->Px();
1022 sumPyDau+=resdau->Py();
1023 sumPzDau+=resdau->Pz();
1024 nKaons++;
1025 arrayDauLab[nFoundKpi++]=indResDau;
1026 if(nFoundKpi>3) return -1;
1027 }
1028 if(TMath::Abs(pdgresdau)==211){
1029 if(pdgD*pdgresdau<0) return -1;
1030 sumPxDau+=resdau->Px();
1031 sumPyDau+=resdau->Py();
1032 sumPzDau+=resdau->Pz();
1033 nPions++;
1034 arrayDauLab[nFoundKpi++]=indResDau;
1035 if(nFoundKpi>3) return -1;
1036 }
1037 }
1038 }else{
1039 return -1;
1040 }
1041 }
1042 if(nPions!=2) return -1;
1043 if(nKaons!=1) return -1;
1044 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1045 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1046 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1047 if(nDau==3) return 1;
1048 else if(nDau==2) return 2;
1049 }
1050
1051 return -1;
1052
5dc814ec 1053}
1054//____________________________________________________________________________
1055Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1056 // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1057
1058 if(label<0) return -1;
1059 TParticle* mcPart = stack->Particle(label);
1060 Int_t pdgD=mcPart->GetPdgCode();
1061 if(TMath::Abs(pdgD)!=411) return -1;
1062
1063 Int_t nDau=mcPart->GetNDaughters();
1064 Int_t labelFirstDau = mcPart->GetDaughter(0);
1065 Int_t nKaons=0;
1066 Int_t nPions=0;
1067 Double_t sumPxDau=0.;
1068 Double_t sumPyDau=0.;
1069 Double_t sumPzDau=0.;
1070 Int_t nFoundKpi=0;
1071 Bool_t isPhi=kFALSE;
1072 Bool_t isk0st=kFALSE;
1073
1074 if(nDau==3 || nDau==2){
1075 for(Int_t iDau=0; iDau<nDau; iDau++){
1076 Int_t indDau = labelFirstDau+iDau;
1077 if(indDau<0) return -1;
1078 TParticle* dau=stack->Particle(indDau);
1079 if(!dau) return -1;
1080 Int_t pdgdau=dau->GetPdgCode();
1081 if(TMath::Abs(pdgdau)==321){
1082 nKaons++;
1083 sumPxDau+=dau->Px();
1084 sumPyDau+=dau->Py();
1085 sumPzDau+=dau->Pz();
1086 arrayDauLab[nFoundKpi++]=indDau;
1087 if(nFoundKpi>3) return -1;
1088 }else if(TMath::Abs(pdgdau)==211){
1089 nPions++;
1090 sumPxDau+=dau->Px();
1091 sumPyDau+=dau->Py();
1092 sumPzDau+=dau->Pz();
1093 arrayDauLab[nFoundKpi++]=indDau;
1094 if(nFoundKpi>3) return -1;
1095 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1096 if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1097 if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1098 Int_t nResDau=dau->GetNDaughters();
1099 if(nResDau!=2) return -1;
1100 Int_t indFirstResDau=dau->GetDaughter(0);
1101 for(Int_t resDau=0; resDau<2; resDau++){
1102 Int_t indResDau=indFirstResDau+resDau;
1103 if(indResDau<0) return -1;
1104 TParticle* resdau=stack->Particle(indResDau);
1105 if(!resdau) return -1;
1106 Int_t pdgresdau=resdau->GetPdgCode();
1107 if(TMath::Abs(pdgresdau)==321){
1108 sumPxDau+=resdau->Px();
1109 sumPyDau+=resdau->Py();
1110 sumPzDau+=resdau->Pz();
1111 nKaons++;
1112 arrayDauLab[nFoundKpi++]=indResDau;
1113 if(nFoundKpi>3) return -1;
1114 }
1115 if(TMath::Abs(pdgresdau)==211){
1116 sumPxDau+=resdau->Px();
1117 sumPyDau+=resdau->Py();
1118 sumPzDau+=resdau->Pz();
1119 nPions++;
1120 arrayDauLab[nFoundKpi++]=indResDau;
1121 if(nFoundKpi>3) return -1;
1122 }
1123 }
1124 }else{
1125 return -1;
1126 }
1127 }
1128 if(nPions!=1) return -1;
1129 if(nKaons!=2) return -1;
1130 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1131 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1132 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1133 if(nDau==3) return 3;
1134 else if(nDau==2){
1135 if(isk0st) return 2;
1136 if(isPhi) return 1;
1137 }
1138 }
1139
1140 return -1;
1141
c8781624 1142}
ee9c9d9b 1143//____________________________________________________________________________
1144Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1145 // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1146
1147 if(label<0) return -1;
1148 TParticle* mcPart = stack->Particle(label);
1149 Int_t pdgD=mcPart->GetPdgCode();
1150 if(TMath::Abs(pdgD)!=431) return -1;
1151
1152 Int_t nDau=mcPart->GetNDaughters();
1153 Int_t labelFirstDau = mcPart->GetDaughter(0);
1154 Int_t nKaons=0;
1155 Int_t nPions=0;
1156 Double_t sumPxDau=0.;
1157 Double_t sumPyDau=0.;
1158 Double_t sumPzDau=0.;
1159 Int_t nFoundKpi=0;
1160 Bool_t isPhi=kFALSE;
1161 Bool_t isk0st=kFALSE;
1162
1163 if(nDau==3 || nDau==2){
1164 for(Int_t iDau=0; iDau<nDau; iDau++){
1165 Int_t indDau = labelFirstDau+iDau;
1166 if(indDau<0) return -1;
1167 TParticle* dau=stack->Particle(indDau);
1168 if(!dau) return -1;
1169 Int_t pdgdau=dau->GetPdgCode();
1170 if(TMath::Abs(pdgdau)==321){
1171 nKaons++;
1172 sumPxDau+=dau->Px();
1173 sumPyDau+=dau->Py();
1174 sumPzDau+=dau->Pz();
1175 arrayDauLab[nFoundKpi++]=indDau;
1176 if(nFoundKpi>3) return -1;
1177 }else if(TMath::Abs(pdgdau)==211){
1178 nPions++;
1179 sumPxDau+=dau->Px();
1180 sumPyDau+=dau->Py();
1181 sumPzDau+=dau->Pz();
1182 arrayDauLab[nFoundKpi++]=indDau;
1183 if(nFoundKpi>3) return -1;
1184 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1185 if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1186 if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1187 Int_t nResDau=dau->GetNDaughters();
1188 if(nResDau!=2) return -1;
1189 Int_t indFirstResDau=dau->GetDaughter(0);
1190 for(Int_t resDau=0; resDau<2; resDau++){
1191 Int_t indResDau=indFirstResDau+resDau;
1192 if(indResDau<0) return -1;
1193 TParticle* resdau=stack->Particle(indResDau);
1194 if(!resdau) return -1;
1195 Int_t pdgresdau=resdau->GetPdgCode();
1196 if(TMath::Abs(pdgresdau)==321){
1197 sumPxDau+=resdau->Px();
1198 sumPyDau+=resdau->Py();
1199 sumPzDau+=resdau->Pz();
1200 nKaons++;
1201 arrayDauLab[nFoundKpi++]=indResDau;
1202 if(nFoundKpi>3) return -1;
1203 }
1204 if(TMath::Abs(pdgresdau)==211){
1205 sumPxDau+=resdau->Px();
1206 sumPyDau+=resdau->Py();
1207 sumPzDau+=resdau->Pz();
1208 nPions++;
1209 arrayDauLab[nFoundKpi++]=indResDau;
1210 if(nFoundKpi>3) return -1;
1211 }
1212 }
1213 }else{
1214 return -1;
1215 }
1216 }
1217 if(nPions!=1) return -1;
1218 if(nKaons!=2) return -1;
1219 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1220 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1221 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1222 if(nDau==3) return 3;
1223 else if(nDau==2){
1224 if(isk0st) return 2;
1225 if(isPhi) return 1;
1226 }
1227 }
1228
1229 return -1;
1230
1231}
1232//____________________________________________________________________________
1233Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1234 // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1235
1236 Int_t pdgD=mcPart->GetPdgCode();
1237 if(TMath::Abs(pdgD)!=431) return -1;
1238
1239 Int_t nDau=mcPart->GetNDaughters();
1240 Int_t labelFirstDau = mcPart->GetDaughter(0);
1241 Int_t nKaons=0;
1242 Int_t nPions=0;
1243 Double_t sumPxDau=0.;
1244 Double_t sumPyDau=0.;
1245 Double_t sumPzDau=0.;
1246 Int_t nFoundKpi=0;
1247 Bool_t isPhi=kFALSE;
1248 Bool_t isk0st=kFALSE;
1249
1250 if(nDau==3 || nDau==2){
1251 for(Int_t iDau=0; iDau<nDau; iDau++){
1252 Int_t indDau = labelFirstDau+iDau;
1253 if(indDau<0) return -1;
1254 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1255 if(!dau) return -1;
1256 Int_t pdgdau=dau->GetPdgCode();
1257 if(TMath::Abs(pdgdau)==321){
1258 nKaons++;
1259 sumPxDau+=dau->Px();
1260 sumPyDau+=dau->Py();
1261 sumPzDau+=dau->Pz();
1262 arrayDauLab[nFoundKpi++]=indDau;
1263 if(nFoundKpi>3) return -1;
1264 }else if(TMath::Abs(pdgdau)==211){
1265 nPions++;
1266 sumPxDau+=dau->Px();
1267 sumPyDau+=dau->Py();
1268 sumPzDau+=dau->Pz();
1269 arrayDauLab[nFoundKpi++]=indDau;
1270 if(nFoundKpi>3) return -1;
1271 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1272 if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1273 if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1274 Int_t nResDau=dau->GetNDaughters();
1275 if(nResDau!=2) return -1;
1276 Int_t indFirstResDau=dau->GetDaughter(0);
1277 for(Int_t resDau=0; resDau<2; resDau++){
1278 Int_t indResDau=indFirstResDau+resDau;
1279 if(indResDau<0) return -1;
1280 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1281 if(!resdau) return -1;
1282 Int_t pdgresdau=resdau->GetPdgCode();
1283 if(TMath::Abs(pdgresdau)==321){
1284 sumPxDau+=resdau->Px();
1285 sumPyDau+=resdau->Py();
1286 sumPzDau+=resdau->Pz();
1287 nKaons++;
1288 arrayDauLab[nFoundKpi++]=indResDau;
1289 if(nFoundKpi>3) return -1;
1290 }
1291 if(TMath::Abs(pdgresdau)==211){
1292 sumPxDau+=resdau->Px();
1293 sumPyDau+=resdau->Py();
1294 sumPzDau+=resdau->Pz();
1295 nPions++;
1296 arrayDauLab[nFoundKpi++]=indResDau;
1297 if(nFoundKpi>3) return -1;
1298 }
5dc814ec 1299 }
ee9c9d9b 1300 }else{
1301 return -1;
1302 }
1303 }
1304 if(nPions!=1) return -1;
1305 if(nKaons!=2) return -1;
1306 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1307 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1308 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1309 if(nDau==3) return 3;
1310 else if(nDau==2){
1311 if(isk0st) return 2;
1312 if(isPhi) return 1;
1313 }
1314 }
1315
1316 return -1;
1317
1318}
1319//____________________________________________________________________________
1320Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1321 // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1322
1323 if(label<0) return -1;
1324 TParticle* mcPart = stack->Particle(label);
1325 Int_t pdgD=mcPart->GetPdgCode();
1326 if(TMath::Abs(pdgD)!=413) return -1;
1327
1328 Int_t nDau=mcPart->GetNDaughters();
1329 if(nDau!=2) return -1;
1330
1331 Int_t labelFirstDau = mcPart->GetDaughter(0);
1332 Int_t nKaons=0;
1333 Int_t nPions=0;
1334 Double_t sumPxDau=0.;
1335 Double_t sumPyDau=0.;
1336 Double_t sumPzDau=0.;
1337 Int_t nFoundKpi=0;
1338
1339 for(Int_t iDau=0; iDau<nDau; iDau++){
1340 Int_t indDau = labelFirstDau+iDau;
1341 if(indDau<0) return -1;
1342 TParticle* dau=stack->Particle(indDau);
1343 if(!dau) return -1;
1344 Int_t pdgdau=dau->GetPdgCode();
1345 if(TMath::Abs(pdgdau)==421){
1346 Int_t nResDau=dau->GetNDaughters();
1347 if(nResDau!=2) return -1;
1348 Int_t indFirstResDau=dau->GetDaughter(0);
1349 for(Int_t resDau=0; resDau<2; resDau++){
1350 Int_t indResDau=indFirstResDau+resDau;
1351 if(indResDau<0) return -1;
1352 TParticle* resdau=stack->Particle(indResDau);
1353 if(!resdau) return -1;
1354 Int_t pdgresdau=resdau->GetPdgCode();
1355 if(TMath::Abs(pdgresdau)==321){
1356 if(pdgD*pdgresdau>0) return -1;
1357 sumPxDau+=resdau->Px();
1358 sumPyDau+=resdau->Py();
1359 sumPzDau+=resdau->Pz();
1360 nKaons++;
1361 arrayDauLab[nFoundKpi++]=indResDau;
1362 if(nFoundKpi>3) return -1;
1363 }
1364 if(TMath::Abs(pdgresdau)==211){
1365 if(pdgD*pdgresdau<0) return -1;
1366 sumPxDau+=resdau->Px();
1367 sumPyDau+=resdau->Py();
1368 sumPzDau+=resdau->Pz();
1369 nPions++;
1370 arrayDauLab[nFoundKpi++]=indResDau;
1371 if(nFoundKpi>3) return -1;
1372 }
1373 }
1374 }else if(TMath::Abs(pdgdau)==211){
1375 if(pdgD*pdgdau<0) return -1;
1376 nPions++;
1377 sumPxDau+=dau->Px();
1378 sumPyDau+=dau->Py();
1379 sumPzDau+=dau->Pz();
1380 arrayDauLab[nFoundKpi++]=indDau;
1381 if(nFoundKpi>3) return -1;
1382 }
1383 }
1384
1385 if(nPions!=2) return -1;
1386 if(nKaons!=1) return -1;
1387 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1388 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1389 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1390 return 1;
1391
1392}
1393
1394//____________________________________________________________________________
1395Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1396 // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1397
1398 Int_t pdgD=mcPart->GetPdgCode();
1399 if(TMath::Abs(pdgD)!=413) return -1;
1400
1401 Int_t nDau=mcPart->GetNDaughters();
1402 if(nDau!=2) return -1;
1403
1404 Int_t labelFirstDau = mcPart->GetDaughter(0);
1405 Int_t nKaons=0;
1406 Int_t nPions=0;
1407 Double_t sumPxDau=0.;
1408 Double_t sumPyDau=0.;
1409 Double_t sumPzDau=0.;
1410 Int_t nFoundKpi=0;
1411
1412 for(Int_t iDau=0; iDau<nDau; iDau++){
1413 Int_t indDau = labelFirstDau+iDau;
1414 if(indDau<0) return -1;
1415 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1416 if(!dau) return -1;
1417 Int_t pdgdau=dau->GetPdgCode();
1418 if(TMath::Abs(pdgdau)==421){
1419 Int_t nResDau=dau->GetNDaughters();
1420 if(nResDau!=2) return -1;
1421 Int_t indFirstResDau=dau->GetDaughter(0);
1422 for(Int_t resDau=0; resDau<2; resDau++){
1423 Int_t indResDau=indFirstResDau+resDau;
1424 if(indResDau<0) return -1;
1425 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1426 if(!resdau) return -1;
1427 Int_t pdgresdau=resdau->GetPdgCode();
1428 if(TMath::Abs(pdgresdau)==321){
1429 if(pdgD*pdgresdau>0) return -1;
1430 sumPxDau+=resdau->Px();
1431 sumPyDau+=resdau->Py();
1432 sumPzDau+=resdau->Pz();
1433 nKaons++;
1434 arrayDauLab[nFoundKpi++]=indResDau;
1435 if(nFoundKpi>3) return -1;
1436 }
1437 if(TMath::Abs(pdgresdau)==211){
1438 if(pdgD*pdgresdau<0) return -1;
1439 sumPxDau+=resdau->Px();
1440 sumPyDau+=resdau->Py();
1441 sumPzDau+=resdau->Pz();
1442 nPions++;
1443 arrayDauLab[nFoundKpi++]=indResDau;
1444 if(nFoundKpi>3) return -1;
1445 }
1446 }
1447 }else if(TMath::Abs(pdgdau)==211){
1448 if(pdgD*pdgdau<0) return -1;
1449 nPions++;
1450 sumPxDau+=dau->Px();
1451 sumPyDau+=dau->Py();
1452 sumPzDau+=dau->Pz();
1453 arrayDauLab[nFoundKpi++]=indDau;
1454 if(nFoundKpi>3) return -1;
1455 }
1456 }
1457
1458 if(nPions!=2) return -1;
1459 if(nKaons!=1) return -1;
1460 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1461 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1462 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1463 return 1;
1464
1465}
1466
5dc814ec 1467//____________________________________________________________________________
1468Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1469 // Checks the Lc->pKpi decay channel. Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases
1470
1471 if(label<0) return -1;
1472 TParticle* mcPart = stack->Particle(label);
1473 Int_t pdgD=mcPart->GetPdgCode();
1474 if(TMath::Abs(pdgD)!=4122) return -1;
1475
1476 Int_t nDau=mcPart->GetNDaughters();
1477 Int_t labelFirstDau = mcPart->GetDaughter(0);
1478 Int_t nKaons=0;
1479 Int_t nPions=0;
1480 Int_t nProtons=0;
1481 Double_t sumPxDau=0.;
1482 Double_t sumPyDau=0.;
1483 Double_t sumPzDau=0.;
1484 Int_t nFoundpKpi=0;
1485
1486 Int_t codeRes=-1;
1487 if(nDau==3 || nDau==2){
1488 for(Int_t iDau=0; iDau<nDau; iDau++){
1489 Int_t indDau = labelFirstDau+iDau;
1490 if(indDau<0) return -1;
1491 TParticle* dau=stack->Particle(indDau);
1492 if(!dau) return -1;
1493 Int_t pdgdau=dau->GetPdgCode();
1494 if(TMath::Abs(pdgdau)==321){
1495 nKaons++;
1496 sumPxDau+=dau->Px();
1497 sumPyDau+=dau->Py();
1498 sumPzDau+=dau->Pz();
1499 arrayDauLab[nFoundpKpi++]=indDau;
1500 if(nFoundpKpi>3) return -1;
1501 }else if(TMath::Abs(pdgdau)==211){
1502 nPions++;
1503 sumPxDau+=dau->Px();
1504 sumPyDau+=dau->Py();
1505 sumPzDau+=dau->Pz();
1506 arrayDauLab[nFoundpKpi++]=indDau;
1507 if(nFoundpKpi>3) return -1;
1508 }else if(TMath::Abs(pdgdau)==2212){
1509 nProtons++;
1510 sumPxDau+=dau->Px();
1511 sumPyDau+=dau->Py();
1512 sumPzDau+=dau->Pz();
1513 arrayDauLab[nFoundpKpi++]=indDau;
1514 if(nFoundpKpi>3) return -1;
1515 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1516 TMath::Abs(pdgdau)==2224){
1517 codeRes=TMath::Abs(pdgdau);
1518 Int_t nResDau=dau->GetNDaughters();
1519 if(nResDau!=2) return -1;
1520 Int_t indFirstResDau=dau->GetDaughter(0);
1521 for(Int_t resDau=0; resDau<2; resDau++){
1522 Int_t indResDau=indFirstResDau+resDau;
1523 if(indResDau<0) return -1;
1524 TParticle* resdau=stack->Particle(indResDau);
1525 if(!resdau) return -1;
1526 Int_t pdgresdau=resdau->GetPdgCode();
1527 if(TMath::Abs(pdgresdau)==321){
1528 sumPxDau+=resdau->Px();
1529 sumPyDau+=resdau->Py();
1530 sumPzDau+=resdau->Pz();
1531 nKaons++;
1532 arrayDauLab[nFoundpKpi++]=indResDau;
1533 if(nFoundpKpi>3) return -1;
1534 }else if(TMath::Abs(pdgresdau)==211){
1535 sumPxDau+=resdau->Px();
1536 sumPyDau+=resdau->Py();
1537 sumPzDau+=resdau->Pz();
1538 nPions++;
1539 arrayDauLab[nFoundpKpi++]=indResDau;
1540 if(nFoundpKpi>3) return -1;
1541 }else if(TMath::Abs(pdgresdau)==2212){
1542 sumPxDau+=resdau->Px();
1543 sumPyDau+=resdau->Py();
1544 sumPzDau+=resdau->Pz();
1545 nProtons++;
1546 arrayDauLab[nFoundpKpi++]=indResDau;
1547 if(nFoundpKpi>3) return -1;
1548 }
1549 }
1550 }else{
1551 return -1;
1552 }
1553 }
1554 if(nPions!=1) return -1;
1555 if(nKaons!=1) return -1;
1556 if(nProtons!=1) return -1;
1557 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1558 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1559 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1560 if(nDau==3) return 1;
1561 else if(nDau==2){
1562 if(codeRes==313) return 2;
1563 else if(codeRes==2224) return 3;
1564 else if(codeRes==3124) return 4;
1565 }
1566 }
1567 return -1;
1568
1569}
1570
1571
1572//____________________________________________________________________________
1573Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1574 // Checks the Lc->V0+bachelor decay channel. Returns 1 for pK0s, 2 for piLambda, 3 for pK0l -1 in other cases
1575
1576 if(label<0) return -1;
1577 TParticle* mcPart = stack->Particle(label);
1578 Int_t pdgD=mcPart->GetPdgCode();
1579 if(TMath::Abs(pdgD)!=4122) return -1;
1580
1581 Int_t nDau=mcPart->GetNDaughters();
1582 Int_t labelFirstDau = mcPart->GetDaughter(0);
1583 Int_t nPions=0;
1584 Int_t nProtons=0;
1585 Double_t sumPxDau=0.;
1586 Double_t sumPyDau=0.;
1587 Double_t sumPzDau=0.;
1588 Int_t nFoundppi=0;
1589
1590 Int_t codeV0=-1;
1591 if(nDau==2){
1592 for(Int_t iDau=0; iDau<nDau; iDau++){
1593 Int_t indDau = labelFirstDau+iDau;
1594 if(indDau<0) return -1;
1595 TParticle* dau=stack->Particle(indDau);
1596 if(!dau) return -1;
1597 Int_t pdgdau=dau->GetPdgCode();
1598 if(TMath::Abs(pdgdau)==211){
1599 nPions++;
1600 sumPxDau+=dau->Px();
1601 sumPyDau+=dau->Py();
1602 sumPzDau+=dau->Pz();
1603 arrayDauLab[nFoundppi++]=indDau;
1604 if(nFoundppi>3) return -1;
1605 }else if(TMath::Abs(pdgdau)==2212){
1606 nProtons++;
1607 sumPxDau+=dau->Px();
1608 sumPyDau+=dau->Py();
1609 sumPzDau+=dau->Pz();
1610 arrayDauLab[nFoundppi++]=indDau;
1611 if(nFoundppi>3) return -1;
1612 }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
1613 codeV0=TMath::Abs(pdgdau);
1614 TParticle* v0=dau;
1615 if(codeV0==311){
1616 Int_t nK0Dau=dau->GetNDaughters();
1617 if(nK0Dau!=1) return -1;
1618 Int_t indK0s=dau->GetDaughter(0);
1619 if(indK0s<0) return -1;
1620 v0=stack->Particle(indK0s);
1621 if(!v0) return -1;
1622 Int_t pdgK0sl=v0->GetPdgCode();
1623 codeV0=TMath::Abs(pdgK0sl);
1624 }
1625 Int_t nV0Dau=v0->GetNDaughters();
1626 if(nV0Dau!=2) return -1;
1627 Int_t indFirstV0Dau=v0->GetDaughter(0);
1628 for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1629 Int_t indV0Dau=indFirstV0Dau+v0Dau;
1630 if(indV0Dau<0) return -1;
1631 TParticle* v0dau=stack->Particle(indV0Dau);
1632 if(!v0dau) return -1;
1633 Int_t pdgv0dau=v0dau->GetPdgCode();
1634 if(TMath::Abs(pdgv0dau)==211){
1635 sumPxDau+=v0dau->Px();
1636 sumPyDau+=v0dau->Py();
1637 sumPzDau+=v0dau->Pz();
1638 nPions++;
1639 arrayDauLab[nFoundppi++]=indV0Dau;
1640 if(nFoundppi>3) return -1;
1641 }else if(TMath::Abs(pdgv0dau)==2212){
1642 sumPxDau+=v0dau->Px();
1643 sumPyDau+=v0dau->Py();
1644 sumPzDau+=v0dau->Pz();
1645 nProtons++;
1646 arrayDauLab[nFoundppi++]=indV0Dau;
1647 if(nFoundppi>3) return -1;
1648 }
1649 }
1650 }else{
1651 return -1;
1652 }
1653 }
1654 if(nPions!=2) return -1;
1655 if(nProtons!=1) return -1;
1656 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1657 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1658 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1659 if(codeV0==310) return 1;
1660 else if(codeV0==3122) return 2;
1661 }
1662 return -1;
1663
1664}
1665
1666