]> 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
Add the libraries needed to run the HF QA task
[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}
242//______________________________________________________________________
19df7ca5 243void 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 244
245 // Compute <pt> from 2D histogram M vs pt
246
247 //Make 2D histos in the desired pt range
248 Int_t start=hMassD->FindBin(ptmin);
249 Int_t end=hMassD->FindBin(ptmax)-1;
250 const Int_t nx=end-start;
251 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()));
252 for(Int_t ix=start;ix<end;ix++){
253 for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
254 hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
255 hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
256 }
257 }
258
259 Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
260 Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
261 Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
262 Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
263 Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
264 Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
265 // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
266
267 Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
19df7ca5 268 Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
86d84fa4 269 Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
270 Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
271 Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
272 Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
273 Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
19df7ca5 274 Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
86d84fa4 275 Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
276 Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
277 // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
278
279 Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
280 Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
281 Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
282 // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
283
284 TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
285 TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
286 TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
287
288 hMptBkgLo->Rebin(rebin);
289 hMptBkgHi->Rebin(rebin);
290 hMptSigReg->Rebin(rebin);
291
292 hMptBkgLo->Sumw2();
293 hMptBkgHi->Sumw2();
294 TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
295 hMptBkgLoScal->Scale(bkgSig/bkgLow);
296 TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
297 hMptBkgHiScal->Scale(bkgSig/bkgHi);
298
299 TH1F* hMptBkgAver=0x0;
300 hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
301 hMptBkgAver->Add(hMptBkgHiScal);
302 hMptBkgAver->Scale(0.5);
303 TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
304 hMptSig->Add(hMptBkgAver,-1.);
305
306 averagePt = hMptSig->GetMean();
307 errorPt = hMptSig->GetMeanError();
308
309 delete hMptBkgLo;
310 delete hMptBkgHi;
311 delete hMptSigReg;
312 delete hMptBkgLoScal;
313 delete hMptBkgHiScal;
314 delete hMptBkgAver;
315 delete hMassDpt;
316 delete hMptSig;
317
15917f37 318}
319//____________________________________________________________________________
e6eaae4c 320Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
321 // check if T0VTX trigger was fired, based on a workaround suggested by Alla
322 const Double32_t *mean = aodEv->GetT0TOF();
94d7319f 323 if(mean && mean[0]<9999.) return kTRUE;
e6eaae4c 324 else return kFALSE;
325}
326//____________________________________________________________________________
15917f37 327Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
328 // true impact parameter calculation for Dzero
329
330 if(!partD || !arrayMC || !mcHeader) return 99999.;
331 Int_t code=partD->GetPdgCode();
332 if(TMath::Abs(code)!=421) return 99999.;
333
334 Double_t vtxTrue[3];
335 mcHeader->GetVertex(vtxTrue);
336 Double_t origD[3];
337 partD->XvYvZv(origD);
338 Short_t charge=partD->Charge();
339 Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
340 for(Int_t iDau=0; iDau<2; iDau++){
341 pXdauTrue[iDau]=0.;
342 pYdauTrue[iDau]=0.;
343 pZdauTrue[iDau]=0.;
344 }
345
346 Int_t nDau=partD->GetNDaughters();
347 Int_t labelFirstDau = partD->GetDaughter(0);
348 if(nDau==2){
349 for(Int_t iDau=0; iDau<2; iDau++){
350 Int_t ind = labelFirstDau+iDau;
351 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
352 if(!part){
353 printf("Daughter particle not found in MC array");
354 return 99999.;
355 }
356 pXdauTrue[iDau]=part->Px();
357 pYdauTrue[iDau]=part->Py();
358 pZdauTrue[iDau]=part->Pz();
359 }
360 }else{
361 return 99999.;
362 }
363
48efc6ad 364 Double_t d0dummy[2]={0.,0.};
15917f37 365 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
366 return aodDvsMC.ImpParXY();
367
86d84fa4 368}
b176db3b 369//____________________________________________________________________________
370Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
371 // true impact parameter calculation for Dplus
372
373 if(!partD || !arrayMC || !mcHeader) return 99999.;
374 Int_t code=partD->GetPdgCode();
375 if(TMath::Abs(code)!=411) return 99999.;
376
377 Double_t vtxTrue[3];
378 mcHeader->GetVertex(vtxTrue);
379 Double_t origD[3];
380 partD->XvYvZv(origD);
381 Short_t charge=partD->Charge();
382 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
383 for(Int_t iDau=0; iDau<3; iDau++){
384 pXdauTrue[iDau]=0.;
385 pYdauTrue[iDau]=0.;
386 pZdauTrue[iDau]=0.;
387 }
388
389 Int_t nDau=partD->GetNDaughters();
390 Int_t labelFirstDau = partD->GetDaughter(0);
391 if(nDau==3){
392 for(Int_t iDau=0; iDau<3; iDau++){
393 Int_t ind = labelFirstDau+iDau;
394 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
395 if(!part){
396 printf("Daughter particle not found in MC array");
397 return 99999.;
398 }
399 pXdauTrue[iDau]=part->Px();
400 pYdauTrue[iDau]=part->Py();
401 pZdauTrue[iDau]=part->Pz();
402 }
403 }else if(nDau==2){
404 Int_t theDau=0;
405 for(Int_t iDau=0; iDau<2; iDau++){
406 Int_t ind = labelFirstDau+iDau;
407 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
408 if(!part){
409 printf("Daughter particle not found in MC array");
410 return 99999.;
411 }
412 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
413 if(pdgCode==211 || pdgCode==321){
414 pXdauTrue[theDau]=part->Px();
415 pYdauTrue[theDau]=part->Py();
416 pZdauTrue[theDau]=part->Pz();
417 ++theDau;
418 }else{
419 Int_t nDauRes=part->GetNDaughters();
420 if(nDauRes==2){
421 Int_t labelFirstDauRes = part->GetDaughter(0);
422 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
423 Int_t indDR = labelFirstDauRes+iDauRes;
424 AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
425 if(!partDR){
426 printf("Daughter particle not found in MC array");
427 return 99999.;
428 }
429
430 Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
431 if(pdgCodeDR==211 || pdgCodeDR==321){
432 pXdauTrue[theDau]=partDR->Px();
433 pYdauTrue[theDau]=partDR->Py();
434 pZdauTrue[theDau]=partDR->Pz();
435 ++theDau;
436 }
437 }
438 }
439 }
440 }
441 if(theDau!=3){
442 printf("Wrong number of decay prongs");
443 return 99999.;
444 }
445 }
446
447 Double_t d0dummy[3]={0.,0.,0.};
448 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
449 return aodDvsMC.ImpParXY();
450
451}
452
453
454
686e4710 455//____________________________________________________________________________
456Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
457 //
458 // Correct the number of accepted tracklets based on a TProfile Hist
459 //
460 //
461
462 if(TMath::Abs(vtxZ)>10.0){
280cc141 463 // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
686e4710 464 return uncorrectedNacc;
465 }
466
467 if(!estimatorAvg){
468 printf("ERROR: Missing TProfile for correction of multiplicity\n");
469 return uncorrectedNacc;
470 }
471
472 Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
473
474 Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
475
476 Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
477
478 return correctedNacc;
479}
7f51a627 480//______________________________________________________________________
481TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
f98f4c22 482 // get the name of the generator that produced a given particle
483
484 Int_t nsumpart=0;
485 TList *lh=header->GetCocktailHeaders();
486 Int_t nh=lh->GetEntries();
487 for(Int_t i=0;i<nh;i++){
488 AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
489 TString genname=gh->GetName();
490 Int_t npart=gh->NProduced();
491 if(label>=nsumpart && label<(nsumpart+npart)) return genname;
492 nsumpart+=npart;
493 }
494 TString empty="";
495 return empty;
496}
7c0e9154 497//_____________________________________________________________________
498void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
499
500 // method to check if a track comes from a given generator
f98f4c22 501
502 Int_t lab=TMath::Abs(track->GetLabel());
7c0e9154 503 nameGen=GetGenerator(lab,header);
f98f4c22 504
505 // Int_t countControl=0;
506
507 while(nameGen.IsWhitespace()){
508 AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
509 if(!mcpart){
510 printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
511 break;
512 }
513 Int_t mother = mcpart->GetMother();
514 if(mother<0){
515 printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
516 break;
517 }
518 lab=mother;
519 nameGen=GetGenerator(mother,header);
520 // countControl++;
521 // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
522 // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
523 // break;
524 // }
525 }
526
7c0e9154 527 return;
528}
529//----------------------------------------------------------------------
530Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
531 // method to check if a track comes from the signal event or from the underlying Hijing event
532 TString nameGen;
533 GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
534
f98f4c22 535 if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
536
537 return kTRUE;
538}
539//____________________________________________________________________________
540Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
541 // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
542
543 Int_t nprongs=cand->GetNProngs();
544 for(Int_t i=0;i<nprongs;i++){
545 AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
546 if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
547 }
548 return kFALSE;
549}
c8781624 550//____________________________________________________________________________
d49d50fc 551Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
552 // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event
553
554 AliAODTrack* bach = cand->GetBachelor();
555 if(IsTrackInjected(bach, header, arrayMC)) {
556 AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
557 return kTRUE;
558 }
559 AliAODv0* v0 = cand->Getv0();
560 Int_t nprongs = v0->GetNProngs();
561 for(Int_t i = 0; i < nprongs; i++){
562 AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
563 if(IsTrackInjected(daugh,header,arrayMC)) {
564 AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
565 return kTRUE;
566 }
567 }
568 return kFALSE;
569}
570//____________________________________________________________________________
ee9c9d9b 571Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
572 // checking whether the mother of the particles come from a charm or a bottom quark
573
574 Int_t pdgGranma = 0;
575 Int_t mother = 0;
576 mother = mcPart->GetFirstMother();
577 Int_t istep = 0;
578 Int_t abspdgGranma =0;
579 Bool_t isFromB=kFALSE;
580 Bool_t isQuarkFound=kFALSE;
581 while (mother >0 ){
582 istep++;
583 TParticle* mcGranma = stack->Particle(mother);
584 if (mcGranma){
585 pdgGranma = mcGranma->GetPdgCode();
586 abspdgGranma = TMath::Abs(pdgGranma);
587 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
588 isFromB=kTRUE;
589 }
590 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
591 mother = mcGranma->GetFirstMother();
592 }else{
593 printf("CheckOrigin: Failed casting the mother particle!");
594 break;
595 }
596 }
597 if(searchUpToQuark && !isQuarkFound) return 0;
598 if(isFromB) return 5;
599 else return 4;
600
601}
602//____________________________________________________________________________
c8781624 603Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
604 // checking whether the mother of the particles come from a charm or a bottom quark
605
606 Int_t pdgGranma = 0;
607 Int_t mother = 0;
608 mother = mcPart->GetMother();
609 Int_t istep = 0;
610 Int_t abspdgGranma =0;
611 Bool_t isFromB=kFALSE;
612 Bool_t isQuarkFound=kFALSE;
613 while (mother >0 ){
614 istep++;
615 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
616 if (mcGranma){
617 pdgGranma = mcGranma->GetPdgCode();
618 abspdgGranma = TMath::Abs(pdgGranma);
619 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
620 isFromB=kTRUE;
621 }
622 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
623 mother = mcGranma->GetMother();
624 }else{
625 printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
626 break;
627 }
628 }
629 if(searchUpToQuark && !isQuarkFound) return 0;
630 if(isFromB) return 5;
631 else return 4;
632
633}
ee9c9d9b 634
c8781624 635//____________________________________________________________________________
ee9c9d9b 636Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
c8781624 637 // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
638
ee9c9d9b 639 if(label<0) return -1;
640 TParticle* mcPart = stack->Particle(label);
641 Int_t pdgD=mcPart->GetPdgCode();
642 if(TMath::Abs(pdgD)!=421) return -1;
643
644 Int_t nDau=mcPart->GetNDaughters();
645
646 if(nDau==2){
647 Int_t daughter0 = mcPart->GetDaughter(0);
648 Int_t daughter1 = mcPart->GetDaughter(1);
649 TParticle* mcPartDaughter0 = stack->Particle(daughter0);
650 TParticle* mcPartDaughter1 = stack->Particle(daughter1);
651 if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
652 arrayDauLab[0]=daughter0;
653 arrayDauLab[1]=daughter1;
654 Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
655 Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
656 if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
657 !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
658 return -1;
659 }
660 if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
661 if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
662 if((pdgCode0*pdgCode1)>0) return -1;
663 Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
664 Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
665 Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
666 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
667 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
668 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
669 return 1;
670 }
671
672 if(nDau==3 || nDau==4){
673 Int_t nKaons=0;
674 Int_t nPions=0;
675 Double_t sumPxDau=0.;
676 Double_t sumPyDau=0.;
677 Double_t sumPzDau=0.;
678 Int_t nFoundKpi=0;
679 Int_t labelFirstDau = mcPart->GetDaughter(0);
680 for(Int_t iDau=0; iDau<nDau; iDau++){
681 Int_t indDau = labelFirstDau+iDau;
682 if(indDau<0) return -1;
683 TParticle* dau=stack->Particle(indDau);
684 if(!dau) return -1;
685 Int_t pdgdau=dau->GetPdgCode();
686 if(TMath::Abs(pdgdau)==321){
687 if(pdgD>0 && pdgdau>0) return -1;
688 if(pdgD<0 && pdgdau<0) return -1;
689 nKaons++;
690 sumPxDau+=dau->Px();
691 sumPyDau+=dau->Py();
692 sumPzDau+=dau->Pz();
693 arrayDauLab[nFoundKpi++]=indDau;
694 if(nFoundKpi>4) return -1;
695 }else if(TMath::Abs(pdgdau)==211){
696 nPions++;
697 sumPxDau+=dau->Px();
698 sumPyDau+=dau->Py();
699 sumPzDau+=dau->Pz();
700 arrayDauLab[nFoundKpi++]=indDau;
701 if(nFoundKpi>4) return -1;
702 }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
703 Int_t nResDau=dau->GetNDaughters();
704 if(nResDau!=2) return -1;
705 Int_t indFirstResDau=dau->GetDaughter(0);
706 for(Int_t resDau=0; resDau<2; resDau++){
707 Int_t indResDau=indFirstResDau+resDau;
708 if(indResDau<0) return -1;
709 TParticle* resdau=stack->Particle(indResDau);
710 if(!resdau) return -1;
711 Int_t pdgresdau=resdau->GetPdgCode();
712 if(TMath::Abs(pdgresdau)==321){
713 if(pdgD>0 && pdgresdau>0) return -1;
714 if(pdgD<0 && pdgresdau<0) return -1;
715 nKaons++;
716 sumPxDau+=resdau->Px();
717 sumPyDau+=resdau->Py();
718 sumPzDau+=resdau->Pz();
719 arrayDauLab[nFoundKpi++]=indResDau;
720 if(nFoundKpi>4) return -1;
721 }
722 if(TMath::Abs(pdgresdau)==211){
723 nPions++;
724 sumPxDau+=resdau->Px();
725 sumPyDau+=resdau->Py();
726 sumPzDau+=resdau->Pz();
727 arrayDauLab[nFoundKpi++]=indResDau;
728 if(nFoundKpi>4) return -1;
729 }
730 }
731 }else{
732 return -1;
733 }
734 }
735 if(nPions!=3) return -1;
736 if(nKaons!=1) return -1;
737 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
738 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
739 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
740 return 2;
741 }
742 return -1;
743}
744
745//____________________________________________________________________________
746Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
747 // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
c8781624 748
749 Int_t pdgD=mcPart->GetPdgCode();
750 if(TMath::Abs(pdgD)!=421) return -1;
751
752 Int_t nDau=mcPart->GetNDaughters();
753
754 if(nDau==2){
755 Int_t daughter0 = mcPart->GetDaughter(0);
756 Int_t daughter1 = mcPart->GetDaughter(1);
757 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
758 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
759 if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
760 arrayDauLab[0]=daughter0;
761 arrayDauLab[1]=daughter1;
762 Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
763 Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
764 if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
765 !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
766 return -1;
767 }
768 if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
769 if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
770 if((pdgCode0*pdgCode1)>0) return -1;
771 Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
772 Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
773 Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
774 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
775 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
776 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
777 return 1;
778 }
779
780 if(nDau==3 || nDau==4){
781 Int_t nKaons=0;
782 Int_t nPions=0;
783 Double_t sumPxDau=0.;
784 Double_t sumPyDau=0.;
785 Double_t sumPzDau=0.;
786 Int_t nFoundKpi=0;
787 Int_t labelFirstDau = mcPart->GetDaughter(0);
788 for(Int_t iDau=0; iDau<nDau; iDau++){
789 Int_t indDau = labelFirstDau+iDau;
790 if(indDau<0) return -1;
791 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
792 if(!dau) return -1;
793 Int_t pdgdau=dau->GetPdgCode();
794 if(TMath::Abs(pdgdau)==321){
795 if(pdgD>0 && pdgdau>0) return -1;
796 if(pdgD<0 && pdgdau<0) return -1;
797 nKaons++;
798 sumPxDau+=dau->Px();
799 sumPyDau+=dau->Py();
800 sumPzDau+=dau->Pz();
801 arrayDauLab[nFoundKpi++]=indDau;
802 if(nFoundKpi>4) return -1;
803 }else if(TMath::Abs(pdgdau)==211){
804 nPions++;
805 sumPxDau+=dau->Px();
806 sumPyDau+=dau->Py();
807 sumPzDau+=dau->Pz();
808 arrayDauLab[nFoundKpi++]=indDau;
809 if(nFoundKpi>4) return -1;
810 }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
811 Int_t nResDau=dau->GetNDaughters();
812 if(nResDau!=2) return -1;
813 Int_t indFirstResDau=dau->GetDaughter(0);
814 for(Int_t resDau=0; resDau<2; resDau++){
815 Int_t indResDau=indFirstResDau+resDau;
816 if(indResDau<0) return -1;
817 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
818 if(!resdau) return -1;
819 Int_t pdgresdau=resdau->GetPdgCode();
820 if(TMath::Abs(pdgresdau)==321){
821 if(pdgD>0 && pdgresdau>0) return -1;
822 if(pdgD<0 && pdgresdau<0) return -1;
823 nKaons++;
824 sumPxDau+=resdau->Px();
825 sumPyDau+=resdau->Py();
826 sumPzDau+=resdau->Pz();
827 arrayDauLab[nFoundKpi++]=indResDau;
828 if(nFoundKpi>4) return -1;
829 }
830 if(TMath::Abs(pdgresdau)==211){
831 nPions++;
832 sumPxDau+=resdau->Px();
833 sumPyDau+=resdau->Py();
834 sumPzDau+=resdau->Pz();
835 arrayDauLab[nFoundKpi++]=indResDau;
836 if(nFoundKpi>4) return -1;
837 }
838 }
839 }else{
840 return -1;
841 }
842 }
843 if(nPions!=3) return -1;
844 if(nKaons!=1) return -1;
845 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
846 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
847 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
848 return 2;
849 }
850 return -1;
851}
ee9c9d9b 852//____________________________________________________________________________
853Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
854 // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
855
856 if(label<0) return -1;
857 TParticle* mcPart = stack->Particle(label);
858 Int_t pdgD=mcPart->GetPdgCode();
859 if(TMath::Abs(pdgD)!=411) return -1;
860
861 Int_t nDau=mcPart->GetNDaughters();
862 Int_t labelFirstDau = mcPart->GetDaughter(0);
863 Int_t nKaons=0;
864 Int_t nPions=0;
865 Double_t sumPxDau=0.;
866 Double_t sumPyDau=0.;
867 Double_t sumPzDau=0.;
868 Int_t nFoundKpi=0;
869
870 if(nDau==3 || nDau==2){
871 for(Int_t iDau=0; iDau<nDau; iDau++){
872 Int_t indDau = labelFirstDau+iDau;
873 if(indDau<0) return -1;
874 TParticle* dau=stack->Particle(indDau);
875 if(!dau) return -1;
876 Int_t pdgdau=dau->GetPdgCode();
877 if(TMath::Abs(pdgdau)==321){
878 if(pdgD*pdgdau>0) return -1;
879 nKaons++;
880 sumPxDau+=dau->Px();
881 sumPyDau+=dau->Py();
882 sumPzDau+=dau->Pz();
883 arrayDauLab[nFoundKpi++]=indDau;
884 if(nFoundKpi>3) return -1;
885 }else if(TMath::Abs(pdgdau)==211){
886 if(pdgD*pdgdau<0) return -1;
887 nPions++;
888 sumPxDau+=dau->Px();
889 sumPyDau+=dau->Py();
890 sumPzDau+=dau->Pz();
891 arrayDauLab[nFoundKpi++]=indDau;
892 if(nFoundKpi>3) return -1;
893 }else if(TMath::Abs(pdgdau)==313){
894 Int_t nResDau=dau->GetNDaughters();
895 if(nResDau!=2) return -1;
896 Int_t indFirstResDau=dau->GetDaughter(0);
897 for(Int_t resDau=0; resDau<2; resDau++){
898 Int_t indResDau=indFirstResDau+resDau;
899 if(indResDau<0) return -1;
900 TParticle* resdau=stack->Particle(indResDau);
901 if(!resdau) return -1;
902 Int_t pdgresdau=resdau->GetPdgCode();
903 if(TMath::Abs(pdgresdau)==321){
904 if(pdgD*pdgresdau>0) return -1;
905 sumPxDau+=resdau->Px();
906 sumPyDau+=resdau->Py();
907 sumPzDau+=resdau->Pz();
908 nKaons++;
909 arrayDauLab[nFoundKpi++]=indResDau;
910 if(nFoundKpi>3) return -1;
911 }
912 if(TMath::Abs(pdgresdau)==211){
913 if(pdgD*pdgresdau<0) return -1;
914 sumPxDau+=resdau->Px();
915 sumPyDau+=resdau->Py();
916 sumPzDau+=resdau->Pz();
917 nPions++;
918 arrayDauLab[nFoundKpi++]=indResDau;
919 if(nFoundKpi>3) return -1;
920 }
921 }
922 }else{
923 return -1;
924 }
925 }
926 if(nPions!=2) return -1;
927 if(nKaons!=1) return -1;
928 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
929 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
930 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
931 if(nDau==3) return 1;
932 else if(nDau==2) return 2;
933 }
934
935 return -1;
936
937}
938
c8781624 939//____________________________________________________________________________
940Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
941 // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
942
943 Int_t pdgD=mcPart->GetPdgCode();
944 if(TMath::Abs(pdgD)!=411) return -1;
945
946 Int_t nDau=mcPart->GetNDaughters();
947 Int_t labelFirstDau = mcPart->GetDaughter(0);
948 Int_t nKaons=0;
949 Int_t nPions=0;
950 Double_t sumPxDau=0.;
951 Double_t sumPyDau=0.;
952 Double_t sumPzDau=0.;
953 Int_t nFoundKpi=0;
954
955 if(nDau==3 || nDau==2){
956 for(Int_t iDau=0; iDau<nDau; iDau++){
957 Int_t indDau = labelFirstDau+iDau;
958 if(indDau<0) return -1;
959 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
960 if(!dau) return -1;
961 Int_t pdgdau=dau->GetPdgCode();
962 if(TMath::Abs(pdgdau)==321){
963 if(pdgD*pdgdau>0) return -1;
964 nKaons++;
965 sumPxDau+=dau->Px();
966 sumPyDau+=dau->Py();
967 sumPzDau+=dau->Pz();
968 arrayDauLab[nFoundKpi++]=indDau;
969 if(nFoundKpi>3) return -1;
970 }else if(TMath::Abs(pdgdau)==211){
971 if(pdgD*pdgdau<0) return -1;
972 nPions++;
973 sumPxDau+=dau->Px();
974 sumPyDau+=dau->Py();
975 sumPzDau+=dau->Pz();
976 arrayDauLab[nFoundKpi++]=indDau;
977 if(nFoundKpi>3) return -1;
978 }else if(TMath::Abs(pdgdau)==313){
979 Int_t nResDau=dau->GetNDaughters();
980 if(nResDau!=2) return -1;
981 Int_t indFirstResDau=dau->GetDaughter(0);
982 for(Int_t resDau=0; resDau<2; resDau++){
983 Int_t indResDau=indFirstResDau+resDau;
984 if(indResDau<0) return -1;
985 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
986 if(!resdau) return -1;
987 Int_t pdgresdau=resdau->GetPdgCode();
988 if(TMath::Abs(pdgresdau)==321){
989 if(pdgD*pdgresdau>0) return -1;
990 sumPxDau+=resdau->Px();
991 sumPyDau+=resdau->Py();
992 sumPzDau+=resdau->Pz();
993 nKaons++;
994 arrayDauLab[nFoundKpi++]=indResDau;
995 if(nFoundKpi>3) return -1;
996 }
997 if(TMath::Abs(pdgresdau)==211){
998 if(pdgD*pdgresdau<0) return -1;
999 sumPxDau+=resdau->Px();
1000 sumPyDau+=resdau->Py();
1001 sumPzDau+=resdau->Pz();
1002 nPions++;
1003 arrayDauLab[nFoundKpi++]=indResDau;
1004 if(nFoundKpi>3) return -1;
1005 }
1006 }
1007 }else{
1008 return -1;
1009 }
1010 }
1011 if(nPions!=2) return -1;
1012 if(nKaons!=1) return -1;
1013 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1014 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1015 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1016 if(nDau==3) return 1;
1017 else if(nDau==2) return 2;
1018 }
1019
1020 return -1;
1021
1022}
ee9c9d9b 1023//____________________________________________________________________________
1024Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1025 // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1026
1027 if(label<0) return -1;
1028 TParticle* mcPart = stack->Particle(label);
1029 Int_t pdgD=mcPart->GetPdgCode();
1030 if(TMath::Abs(pdgD)!=431) return -1;
1031
1032 Int_t nDau=mcPart->GetNDaughters();
1033 Int_t labelFirstDau = mcPart->GetDaughter(0);
1034 Int_t nKaons=0;
1035 Int_t nPions=0;
1036 Double_t sumPxDau=0.;
1037 Double_t sumPyDau=0.;
1038 Double_t sumPzDau=0.;
1039 Int_t nFoundKpi=0;
1040 Bool_t isPhi=kFALSE;
1041 Bool_t isk0st=kFALSE;
1042
1043 if(nDau==3 || nDau==2){
1044 for(Int_t iDau=0; iDau<nDau; iDau++){
1045 Int_t indDau = labelFirstDau+iDau;
1046 if(indDau<0) return -1;
1047 TParticle* dau=stack->Particle(indDau);
1048 if(!dau) return -1;
1049 Int_t pdgdau=dau->GetPdgCode();
1050 if(TMath::Abs(pdgdau)==321){
1051 nKaons++;
1052 sumPxDau+=dau->Px();
1053 sumPyDau+=dau->Py();
1054 sumPzDau+=dau->Pz();
1055 arrayDauLab[nFoundKpi++]=indDau;
1056 if(nFoundKpi>3) return -1;
1057 }else if(TMath::Abs(pdgdau)==211){
1058 nPions++;
1059 sumPxDau+=dau->Px();
1060 sumPyDau+=dau->Py();
1061 sumPzDau+=dau->Pz();
1062 arrayDauLab[nFoundKpi++]=indDau;
1063 if(nFoundKpi>3) return -1;
1064 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1065 if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1066 if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1067 Int_t nResDau=dau->GetNDaughters();
1068 if(nResDau!=2) return -1;
1069 Int_t indFirstResDau=dau->GetDaughter(0);
1070 for(Int_t resDau=0; resDau<2; resDau++){
1071 Int_t indResDau=indFirstResDau+resDau;
1072 if(indResDau<0) return -1;
1073 TParticle* resdau=stack->Particle(indResDau);
1074 if(!resdau) return -1;
1075 Int_t pdgresdau=resdau->GetPdgCode();
1076 if(TMath::Abs(pdgresdau)==321){
1077 sumPxDau+=resdau->Px();
1078 sumPyDau+=resdau->Py();
1079 sumPzDau+=resdau->Pz();
1080 nKaons++;
1081 arrayDauLab[nFoundKpi++]=indResDau;
1082 if(nFoundKpi>3) return -1;
1083 }
1084 if(TMath::Abs(pdgresdau)==211){
1085 sumPxDau+=resdau->Px();
1086 sumPyDau+=resdau->Py();
1087 sumPzDau+=resdau->Pz();
1088 nPions++;
1089 arrayDauLab[nFoundKpi++]=indResDau;
1090 if(nFoundKpi>3) return -1;
1091 }
1092 }
1093 }else{
1094 return -1;
1095 }
1096 }
1097 if(nPions!=1) return -1;
1098 if(nKaons!=2) return -1;
1099 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1100 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1101 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1102 if(nDau==3) return 3;
1103 else if(nDau==2){
1104 if(isk0st) return 2;
1105 if(isPhi) return 1;
1106 }
1107 }
1108
1109 return -1;
1110
1111}
1112//____________________________________________________________________________
1113Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1114 // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
1115
1116 Int_t pdgD=mcPart->GetPdgCode();
1117 if(TMath::Abs(pdgD)!=431) return -1;
1118
1119 Int_t nDau=mcPart->GetNDaughters();
1120 Int_t labelFirstDau = mcPart->GetDaughter(0);
1121 Int_t nKaons=0;
1122 Int_t nPions=0;
1123 Double_t sumPxDau=0.;
1124 Double_t sumPyDau=0.;
1125 Double_t sumPzDau=0.;
1126 Int_t nFoundKpi=0;
1127 Bool_t isPhi=kFALSE;
1128 Bool_t isk0st=kFALSE;
1129
1130 if(nDau==3 || nDau==2){
1131 for(Int_t iDau=0; iDau<nDau; iDau++){
1132 Int_t indDau = labelFirstDau+iDau;
1133 if(indDau<0) return -1;
1134 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1135 if(!dau) return -1;
1136 Int_t pdgdau=dau->GetPdgCode();
1137 if(TMath::Abs(pdgdau)==321){
1138 nKaons++;
1139 sumPxDau+=dau->Px();
1140 sumPyDau+=dau->Py();
1141 sumPzDau+=dau->Pz();
1142 arrayDauLab[nFoundKpi++]=indDau;
1143 if(nFoundKpi>3) return -1;
1144 }else if(TMath::Abs(pdgdau)==211){
1145 nPions++;
1146 sumPxDau+=dau->Px();
1147 sumPyDau+=dau->Py();
1148 sumPzDau+=dau->Pz();
1149 arrayDauLab[nFoundKpi++]=indDau;
1150 if(nFoundKpi>3) return -1;
1151 }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1152 if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1153 if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1154 Int_t nResDau=dau->GetNDaughters();
1155 if(nResDau!=2) return -1;
1156 Int_t indFirstResDau=dau->GetDaughter(0);
1157 for(Int_t resDau=0; resDau<2; resDau++){
1158 Int_t indResDau=indFirstResDau+resDau;
1159 if(indResDau<0) return -1;
1160 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1161 if(!resdau) return -1;
1162 Int_t pdgresdau=resdau->GetPdgCode();
1163 if(TMath::Abs(pdgresdau)==321){
1164 sumPxDau+=resdau->Px();
1165 sumPyDau+=resdau->Py();
1166 sumPzDau+=resdau->Pz();
1167 nKaons++;
1168 arrayDauLab[nFoundKpi++]=indResDau;
1169 if(nFoundKpi>3) return -1;
1170 }
1171 if(TMath::Abs(pdgresdau)==211){
1172 sumPxDau+=resdau->Px();
1173 sumPyDau+=resdau->Py();
1174 sumPzDau+=resdau->Pz();
1175 nPions++;
1176 arrayDauLab[nFoundKpi++]=indResDau;
1177 if(nFoundKpi>3) return -1;
1178 }
1179 }
1180 }else{
1181 return -1;
1182 }
1183 }
1184 if(nPions!=1) return -1;
1185 if(nKaons!=2) return -1;
1186 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1187 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1188 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1189 if(nDau==3) return 3;
1190 else if(nDau==2){
1191 if(isk0st) return 2;
1192 if(isPhi) return 1;
1193 }
1194 }
1195
1196 return -1;
1197
1198}
1199//____________________________________________________________________________
1200Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1201 // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1202
1203 if(label<0) return -1;
1204 TParticle* mcPart = stack->Particle(label);
1205 Int_t pdgD=mcPart->GetPdgCode();
1206 if(TMath::Abs(pdgD)!=413) return -1;
1207
1208 Int_t nDau=mcPart->GetNDaughters();
1209 if(nDau!=2) return -1;
1210
1211 Int_t labelFirstDau = mcPart->GetDaughter(0);
1212 Int_t nKaons=0;
1213 Int_t nPions=0;
1214 Double_t sumPxDau=0.;
1215 Double_t sumPyDau=0.;
1216 Double_t sumPzDau=0.;
1217 Int_t nFoundKpi=0;
1218
1219 for(Int_t iDau=0; iDau<nDau; iDau++){
1220 Int_t indDau = labelFirstDau+iDau;
1221 if(indDau<0) return -1;
1222 TParticle* dau=stack->Particle(indDau);
1223 if(!dau) return -1;
1224 Int_t pdgdau=dau->GetPdgCode();
1225 if(TMath::Abs(pdgdau)==421){
1226 Int_t nResDau=dau->GetNDaughters();
1227 if(nResDau!=2) return -1;
1228 Int_t indFirstResDau=dau->GetDaughter(0);
1229 for(Int_t resDau=0; resDau<2; resDau++){
1230 Int_t indResDau=indFirstResDau+resDau;
1231 if(indResDau<0) return -1;
1232 TParticle* resdau=stack->Particle(indResDau);
1233 if(!resdau) return -1;
1234 Int_t pdgresdau=resdau->GetPdgCode();
1235 if(TMath::Abs(pdgresdau)==321){
1236 if(pdgD*pdgresdau>0) return -1;
1237 sumPxDau+=resdau->Px();
1238 sumPyDau+=resdau->Py();
1239 sumPzDau+=resdau->Pz();
1240 nKaons++;
1241 arrayDauLab[nFoundKpi++]=indResDau;
1242 if(nFoundKpi>3) return -1;
1243 }
1244 if(TMath::Abs(pdgresdau)==211){
1245 if(pdgD*pdgresdau<0) return -1;
1246 sumPxDau+=resdau->Px();
1247 sumPyDau+=resdau->Py();
1248 sumPzDau+=resdau->Pz();
1249 nPions++;
1250 arrayDauLab[nFoundKpi++]=indResDau;
1251 if(nFoundKpi>3) return -1;
1252 }
1253 }
1254 }else if(TMath::Abs(pdgdau)==211){
1255 if(pdgD*pdgdau<0) return -1;
1256 nPions++;
1257 sumPxDau+=dau->Px();
1258 sumPyDau+=dau->Py();
1259 sumPzDau+=dau->Pz();
1260 arrayDauLab[nFoundKpi++]=indDau;
1261 if(nFoundKpi>3) return -1;
1262 }
1263 }
1264
1265 if(nPions!=2) return -1;
1266 if(nKaons!=1) return -1;
1267 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1268 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1269 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1270 return 1;
1271
1272}
1273
1274//____________________________________________________________________________
1275Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1276 // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
1277
1278 Int_t pdgD=mcPart->GetPdgCode();
1279 if(TMath::Abs(pdgD)!=413) return -1;
1280
1281 Int_t nDau=mcPart->GetNDaughters();
1282 if(nDau!=2) return -1;
1283
1284 Int_t labelFirstDau = mcPart->GetDaughter(0);
1285 Int_t nKaons=0;
1286 Int_t nPions=0;
1287 Double_t sumPxDau=0.;
1288 Double_t sumPyDau=0.;
1289 Double_t sumPzDau=0.;
1290 Int_t nFoundKpi=0;
1291
1292 for(Int_t iDau=0; iDau<nDau; iDau++){
1293 Int_t indDau = labelFirstDau+iDau;
1294 if(indDau<0) return -1;
1295 AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1296 if(!dau) return -1;
1297 Int_t pdgdau=dau->GetPdgCode();
1298 if(TMath::Abs(pdgdau)==421){
1299 Int_t nResDau=dau->GetNDaughters();
1300 if(nResDau!=2) return -1;
1301 Int_t indFirstResDau=dau->GetDaughter(0);
1302 for(Int_t resDau=0; resDau<2; resDau++){
1303 Int_t indResDau=indFirstResDau+resDau;
1304 if(indResDau<0) return -1;
1305 AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1306 if(!resdau) return -1;
1307 Int_t pdgresdau=resdau->GetPdgCode();
1308 if(TMath::Abs(pdgresdau)==321){
1309 if(pdgD*pdgresdau>0) return -1;
1310 sumPxDau+=resdau->Px();
1311 sumPyDau+=resdau->Py();
1312 sumPzDau+=resdau->Pz();
1313 nKaons++;
1314 arrayDauLab[nFoundKpi++]=indResDau;
1315 if(nFoundKpi>3) return -1;
1316 }
1317 if(TMath::Abs(pdgresdau)==211){
1318 if(pdgD*pdgresdau<0) return -1;
1319 sumPxDau+=resdau->Px();
1320 sumPyDau+=resdau->Py();
1321 sumPzDau+=resdau->Pz();
1322 nPions++;
1323 arrayDauLab[nFoundKpi++]=indResDau;
1324 if(nFoundKpi>3) return -1;
1325 }
1326 }
1327 }else if(TMath::Abs(pdgdau)==211){
1328 if(pdgD*pdgdau<0) return -1;
1329 nPions++;
1330 sumPxDau+=dau->Px();
1331 sumPyDau+=dau->Py();
1332 sumPzDau+=dau->Pz();
1333 arrayDauLab[nFoundKpi++]=indDau;
1334 if(nFoundKpi>3) return -1;
1335 }
1336 }
1337
1338 if(nPions!=2) return -1;
1339 if(nKaons!=1) return -1;
1340 if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1341 if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1342 if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1343 return 1;
1344
1345}
1346