AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
AliVertexingHFUtils.cxx
Go to the documentation of this file.
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>
17 #include <TRandom.h>
18 #include <TProfile.h>
19 #include <TClonesArray.h>
20 #include <TH1F.h>
21 #include <TH2F.h>
22 #include <TF1.h>
23 #include <TParticle.h>
24 #include "AliStack.h"
25 #include "AliAODEvent.h"
26 #include "AliAODMCHeader.h"
27 #include "AliGenEventHeader.h"
28 #include "AliAODMCParticle.h"
29 #include "AliAODRecoDecayHF.h"
30 #include "AliVertexingHFUtils.h"
31 
32 /* $Id$ */
33 
35 // //
36 // Class with functions useful for different D2H analyses //
37 // - event plane resolution //
38 // - <pt> calculation with side band subtraction //
39 // - tracklet multiplicity calculation //
40 // Origin: F.Prino, Torino, prino@to.infn.it //
41 // //
43 
47 
48 
49 //______________________________________________________________________
51  fK(1),
52  fSubRes(1.),
53  fMinEtaForTracklets(-1.),
54  fMaxEtaForTracklets(1.)
55 {
57 }
58 
59 
60 //______________________________________________________________________
62  TObject(),
63  fK(k),
64  fSubRes(1.),
65  fMinEtaForTracklets(-1.),
66  fMaxEtaForTracklets(1.)
67 {
69 }
70 
71 
72 //______________________________________________________________________
73 void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
75 
76 
77  Double_t errSigSq=errsignal*errsignal;
78  Double_t errBkgSq=errbackground*errbackground;
79  Double_t sigPlusBkg=signal+background;
80  if(sigPlusBkg>0. && signal>0.){
81  significance = signal/TMath::Sqrt(signal+background);
82  errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
83  }else{
84  significance=0.;
85  errsignificance=0.;
86  }
87  return;
88 
89 }
90 //______________________________________________________________________
91 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
93  Double_t c[5];
94  if(k==1){
95  c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
96  }
97  else if(k==2){
98  c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
99  }
100  else return -1;
101  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;
102 }
103 
104 //______________________________________________________________________
105 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
106  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
107 }
108 
109 
110 //______________________________________________________________________
111 Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
113 
114  Double_t x1=0;
115  Double_t x2=15;
116  Double_t y1,y2;
117  if(k==1){
118  y1=ResolK1(x1)-res;
119  y2=ResolK1(x2)-res;
120  }
121  else if(k==2){
122  y1=Pol(x1,2)-res;
123  y2=Pol(x2,2)-res;
124  }
125  else return -1;
126 
127  if(y1*y2>0) return -1;
128  if(y1==0) return y1;
129  if(y2==0) return y2;
130  Double_t xmed,ymed;
131  Int_t jiter=0;
132  while((x2-x1)>0.0001){
133  xmed=0.5*(x1+x2);
134  if(k==1){
135  y1=ResolK1(x1)-res;
136  y2=ResolK1(x2)-res;
137  ymed=ResolK1(xmed)-res;
138  }
139  else if(k==2){
140  y1=Pol(x1,2)-res;
141  y2=Pol(x2,2)-res;
142  ymed=Pol(xmed,2)-res;
143  }
144  else return -1;
145  if((y1*ymed)<0) x2=xmed;
146  if((y2*ymed)<0)x1=xmed;
147  if(ymed==0) return xmed;
148  jiter++;
149  }
150  return 0.5*(x1+x2);
151 }
152 
153 //______________________________________________________________________
154 Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
156  Double_t chisub=FindChi(resSub,k);
157  Double_t chifull=chisub*TMath::Sqrt(2);
158  if(k==1) return ResolK1(chifull);
159  else if(k==2) return Pol(chifull,2);
160  else{
161  printf("k should be <=2\n");
162  return 1.;
163  }
164 }
165 
166 //______________________________________________________________________
167 Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
169  if(!hSubEvCorr) return 1.;
170  Double_t resSub=GetSubEvResol(hSubEvCorr);
171  return GetFullEvResol(resSub,k);
172 }
173 //______________________________________________________________________
174 Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
176  if(!hSubEvCorr) return 1.;
177  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
178  printf("%f\n",resSub);
179  return GetFullEvResol(resSub,k);
180 }
181 //______________________________________________________________________
182 Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
184  if(!hSubEvCorr) return 1.;
185  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
186  printf("%f\n",resSub);
187  return GetFullEvResol(resSub,k);
188 }
189 //______________________________________________________________________
190 Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
192  AliAODTracklets* tracklets=ev->GetTracklets();
193  Int_t nTr=tracklets->GetNumberOfTracklets();
194  Int_t count=0;
195  for(Int_t iTr=0; iTr<nTr; iTr++){
196  Double_t theta=tracklets->GetTheta(iTr);
197  Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
198  if(eta>mineta && eta<maxeta) count++;
199  }
200  return count;
201 }
202 //______________________________________________________________________
203 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
205 
206  Int_t nChargedMC=0;
207  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
208  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
209  Int_t charge = part->Charge();
210  Double_t eta = part->Eta();
211  if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
212  }
213  return nChargedMC;
214 }
215 //______________________________________________________________________
216 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
218 
219  Int_t nChargedMC=0;
220  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
221  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
222  Int_t charge = part->Charge();
223  Double_t eta = part->Eta();
224  if(charge!=0 && eta>mineta && eta<maxeta){
225  if(part->IsPrimary())nChargedMC++;
226  }
227  }
228  return nChargedMC;
229 }
230 //______________________________________________________________________
231 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
233 
234  Int_t nChargedMC=0;
235  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
236  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
237  Int_t charge = part->Charge();
238  Double_t eta = part->Eta();
239  if(charge!=0 && eta>mineta && eta<maxeta){
240  if(part->IsPhysicalPrimary())nChargedMC++;
241  }
242  }
243  return nChargedMC;
244 }
245 
246 
247 //______________________________________________________________________
249  //
252  // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
253 
254  Double_t multV0AEq=0;
255  for(Int_t iCh = 32; iCh < 64; ++iCh) {
256  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
257  multV0AEq += mult;
258  }
259  return multV0AEq;
260 }
261 
262 //______________________________________________________________________
267 
268  Double_t multV0CEq=0;
269  for(Int_t iCh = 0; iCh < 32; ++iCh) {
270  Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
271  multV0CEq += mult;
272  }
273  return multV0CEq;
274 }
275 
276 //______________________________________________________________________
277 void 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){
278 
280 
281  //Make 2D histos in the desired pt range
282  Int_t start=hMassD->FindBin(ptmin);
283  Int_t end=hMassD->FindBin(ptmax)-1;
284  const Int_t nx=end-start;
285  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()));
286  for(Int_t ix=start;ix<end;ix++){
287  for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
288  hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
289  hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
290  }
291  }
292 
293  Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
294  Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
295  Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
296  Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
297  Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
298  Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
299  // printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
300 
301  Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
302  Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
303  Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
304  Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
305  Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
306  Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
307  Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
308  Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
309  Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
310  Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
311  // printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
312 
313  Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
314  Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
315  Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
316  // printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
317 
318  TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
319  TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
320  TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
321 
322  hMptBkgLo->Rebin(rebin);
323  hMptBkgHi->Rebin(rebin);
324  hMptSigReg->Rebin(rebin);
325 
326  hMptBkgLo->Sumw2();
327  hMptBkgHi->Sumw2();
328  TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
329  hMptBkgLoScal->Scale(bkgSig/bkgLow);
330  TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
331  hMptBkgHiScal->Scale(bkgSig/bkgHi);
332 
333  TH1F* hMptBkgAver=0x0;
334  hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
335  hMptBkgAver->Add(hMptBkgHiScal);
336  hMptBkgAver->Scale(0.5);
337  TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
338  hMptSig->Add(hMptBkgAver,-1.);
339 
340  averagePt = hMptSig->GetMean();
341  errorPt = hMptSig->GetMeanError();
342 
343  delete hMptBkgLo;
344  delete hMptBkgHi;
345  delete hMptSigReg;
346  delete hMptBkgLoScal;
347  delete hMptBkgHiScal;
348  delete hMptBkgAver;
349  delete hMassDpt;
350  delete hMptSig;
351 
352 }
353 //____________________________________________________________________________
354 Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
356  const Double32_t *mean = aodEv->GetT0TOF();
357  if(mean && mean[0]<9999.) return kTRUE;
358  else return kFALSE;
359 }
360 //____________________________________________________________________________
361 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
363 
364  if(!partD || !arrayMC || !mcHeader) return 99999.;
365  Int_t code=partD->GetPdgCode();
366  if(TMath::Abs(code)!=421) return 99999.;
367 
368  Double_t vtxTrue[3];
369  mcHeader->GetVertex(vtxTrue);
370  Double_t origD[3];
371  partD->XvYvZv(origD);
372  Short_t charge=partD->Charge();
373  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
374  for(Int_t iDau=0; iDau<2; iDau++){
375  pXdauTrue[iDau]=0.;
376  pYdauTrue[iDau]=0.;
377  pZdauTrue[iDau]=0.;
378  }
379 
380  Int_t nDau=partD->GetNDaughters();
381  Int_t labelFirstDau = partD->GetDaughter(0);
382  if(nDau==2){
383  for(Int_t iDau=0; iDau<2; iDau++){
384  Int_t ind = labelFirstDau+iDau;
385  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
386  if(!part){
387  printf("Daughter particle not found in MC array");
388  return 99999.;
389  }
390  pXdauTrue[iDau]=part->Px();
391  pYdauTrue[iDau]=part->Py();
392  pZdauTrue[iDau]=part->Pz();
393  }
394  }else{
395  return 99999.;
396  }
397 
398  Double_t d0dummy[2]={0.,0.};
399  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
400  return aodDvsMC.ImpParXY();
401 
402 }
403 //____________________________________________________________________________
404 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
406 
407  if(!partD || !arrayMC || !mcHeader) return 99999.;
408  Int_t code=partD->GetPdgCode();
409  if(TMath::Abs(code)!=411) return 99999.;
410 
411  Double_t vtxTrue[3];
412  mcHeader->GetVertex(vtxTrue);
413  Double_t origD[3];
414  partD->XvYvZv(origD);
415  Short_t charge=partD->Charge();
416  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
417  for(Int_t iDau=0; iDau<3; iDau++){
418  pXdauTrue[iDau]=0.;
419  pYdauTrue[iDau]=0.;
420  pZdauTrue[iDau]=0.;
421  }
422 
423  Int_t nDau=partD->GetNDaughters();
424  Int_t labelFirstDau = partD->GetDaughter(0);
425  if(nDau==3){
426  for(Int_t iDau=0; iDau<3; iDau++){
427  Int_t ind = labelFirstDau+iDau;
428  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
429  if(!part){
430  printf("Daughter particle not found in MC array");
431  return 99999.;
432  }
433  pXdauTrue[iDau]=part->Px();
434  pYdauTrue[iDau]=part->Py();
435  pZdauTrue[iDau]=part->Pz();
436  }
437  }else if(nDau==2){
438  Int_t theDau=0;
439  for(Int_t iDau=0; iDau<2; iDau++){
440  Int_t ind = labelFirstDau+iDau;
441  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
442  if(!part){
443  printf("Daughter particle not found in MC array");
444  return 99999.;
445  }
446  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
447  if(pdgCode==211 || pdgCode==321){
448  pXdauTrue[theDau]=part->Px();
449  pYdauTrue[theDau]=part->Py();
450  pZdauTrue[theDau]=part->Pz();
451  ++theDau;
452  }else{
453  Int_t nDauRes=part->GetNDaughters();
454  if(nDauRes==2){
455  Int_t labelFirstDauRes = part->GetDaughter(0);
456  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
457  Int_t indDR = labelFirstDauRes+iDauRes;
458  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
459  if(!partDR){
460  printf("Daughter particle not found in MC array");
461  return 99999.;
462  }
463 
464  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
465  if(pdgCodeDR==211 || pdgCodeDR==321){
466  pXdauTrue[theDau]=partDR->Px();
467  pYdauTrue[theDau]=partDR->Py();
468  pZdauTrue[theDau]=partDR->Pz();
469  ++theDau;
470  }
471  }
472  }
473  }
474  }
475  if(theDau!=3){
476  printf("Wrong number of decay prongs");
477  return 99999.;
478  }
479  }
480 
481  Double_t d0dummy[3]={0.,0.,0.};
482  AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
483  return aodDvsMC.ImpParXY();
484 
485 }
486 
487 
488 
489 //____________________________________________________________________________
490 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
491  //
492  // Correct the number of accepted tracklets based on a TProfile Hist
493  //
494  //
495 
496  if(TMath::Abs(vtxZ)>10.0){
497  // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
498  return uncorrectedNacc;
499  }
500 
501  if(!estimatorAvg){
502  printf("ERROR: Missing TProfile for correction of multiplicity\n");
503  return uncorrectedNacc;
504  }
505 
506  Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
507 
508  Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
509 
510  Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
511 
512  return correctedNacc;
513 }
514 //______________________________________________________________________
515 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
517 
518  Int_t nsumpart=0;
519  TList *lh=header->GetCocktailHeaders();
520  Int_t nh=lh->GetEntries();
521  for(Int_t i=0;i<nh;i++){
522  AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
523  TString genname=gh->GetName();
524  Int_t npart=gh->NProduced();
525  if(label>=nsumpart && label<(nsumpart+npart)) return genname;
526  nsumpart+=npart;
527  }
528  TString empty="";
529  return empty;
530 }
531 //_____________________________________________________________________
532 void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
533 
535 
536  Int_t lab=TMath::Abs(track->GetLabel());
537  nameGen=GetGenerator(lab,header);
538 
539  // Int_t countControl=0;
540 
541  while(nameGen.IsWhitespace()){
542  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
543  if(!mcpart){
544  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
545  break;
546  }
547  Int_t mother = mcpart->GetMother();
548  if(mother<0){
549  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
550  break;
551  }
552  lab=mother;
553  nameGen=GetGenerator(mother,header);
554  // countControl++;
555  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
556  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
557  // break;
558  // }
559  }
560 
561  return;
562 }
563 //----------------------------------------------------------------------
564 Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
566  TString nameGen;
567  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
568 
569  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
570 
571  return kTRUE;
572 }
573 //____________________________________________________________________________
574 Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
576 
577  Int_t nprongs=cand->GetNProngs();
578  for(Int_t i=0;i<nprongs;i++){
579  AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
580  if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
581  }
582  return kFALSE;
583 }
584 //____________________________________________________________________________
585 Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
587 
588  AliAODTrack* bach = cand->GetBachelor();
589  if(IsTrackInjected(bach, header, arrayMC)) {
590  AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
591  return kTRUE;
592  }
593  AliAODv0* v0 = cand->Getv0();
594  Int_t nprongs = v0->GetNProngs();
595  for(Int_t i = 0; i < nprongs; i++){
596  AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
597  if(IsTrackInjected(daugh,header,arrayMC)) {
598  AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
599  return kTRUE;
600  }
601  }
602  return kFALSE;
603 }
604 //____________________________________________________________________________
605 Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
607 
608  Int_t pdgGranma = 0;
609  Int_t mother = 0;
610  mother = mcPart->GetFirstMother();
611  Int_t istep = 0;
612  Int_t abspdgGranma =0;
613  Bool_t isFromB=kFALSE;
614  Bool_t isQuarkFound=kFALSE;
615  while (mother >0 ){
616  istep++;
617  TParticle* mcGranma = stack->Particle(mother);
618  if (mcGranma){
619  pdgGranma = mcGranma->GetPdgCode();
620  abspdgGranma = TMath::Abs(pdgGranma);
621  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
622  isFromB=kTRUE;
623  }
624  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
625  mother = mcGranma->GetFirstMother();
626  }else{
627  printf("CheckOrigin: Failed casting the mother particle!");
628  break;
629  }
630  }
631  if(searchUpToQuark && !isQuarkFound) return 0;
632  if(isFromB) return 5;
633  else return 4;
634 
635 }
636 //____________________________________________________________________________
637 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
639 
640  Int_t pdgGranma = 0;
641  Int_t mother = 0;
642  mother = mcPart->GetMother();
643  Int_t istep = 0;
644  Int_t abspdgGranma =0;
645  Bool_t isFromB=kFALSE;
646  Bool_t isQuarkFound=kFALSE;
647  while (mother >0 ){
648  istep++;
649  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
650  if (mcGranma){
651  pdgGranma = mcGranma->GetPdgCode();
652  abspdgGranma = TMath::Abs(pdgGranma);
653  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
654  isFromB=kTRUE;
655  }
656  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
657  mother = mcGranma->GetMother();
658  }else{
659  printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
660  break;
661  }
662  }
663  if(searchUpToQuark && !isQuarkFound) return 0;
664  if(isFromB) return 5;
665  else return 4;
666 
667 }
668 //____________________________________________________________________________
669 Double_t AliVertexingHFUtils::GetBeautyMotherPt(TClonesArray* arrayMC, AliAODMCParticle *mcPart){
671 
672  Int_t pdgGranma = 0;
673  Int_t mother = 0;
674  mother = mcPart->GetMother();
675  Int_t istep = 0;
676  Int_t abspdgGranma =0;
677  while (mother >0 ){
678  istep++;
679  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
680  if (mcGranma){
681  pdgGranma = mcGranma->GetPdgCode();
682  abspdgGranma = TMath::Abs(pdgGranma);
683  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
684  return mcGranma->Pt();
685  }
686  if(abspdgGranma==4) return -999.;
687  if(abspdgGranma==5) return -1.;
688  mother = mcGranma->GetMother();
689  }else{
690  printf("AliVertexingHFUtils::GetBeautyMotherPt: Failed casting the mother particle!");
691  break;
692  }
693  }
694  return -999.;
695 }
696 //____________________________________________________________________________
697 Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
699 
700  if(label<0) return -1;
701  TParticle* mcPart = stack->Particle(label);
702  Int_t pdgD=mcPart->GetPdgCode();
703  if(TMath::Abs(pdgD)!=421) return -1;
704 
705  Int_t nDau=mcPart->GetNDaughters();
706 
707  if(nDau==2){
708  Int_t daughter0 = mcPart->GetDaughter(0);
709  Int_t daughter1 = mcPart->GetDaughter(1);
710  TParticle* mcPartDaughter0 = stack->Particle(daughter0);
711  TParticle* mcPartDaughter1 = stack->Particle(daughter1);
712  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
713  arrayDauLab[0]=daughter0;
714  arrayDauLab[1]=daughter1;
715  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
716  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
717  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
718  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
719  return -1;
720  }
721  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
722  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
723  if((pdgCode0*pdgCode1)>0) return -1;
724  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
725  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
726  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
727  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
728  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
729  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
730  return 1;
731  }
732 
733  if(nDau==3 || nDau==4){
734  Int_t nKaons=0;
735  Int_t nPions=0;
736  Double_t sumPxDau=0.;
737  Double_t sumPyDau=0.;
738  Double_t sumPzDau=0.;
739  Int_t nFoundKpi=0;
740  Int_t labelFirstDau = mcPart->GetDaughter(0);
741  for(Int_t iDau=0; iDau<nDau; iDau++){
742  Int_t indDau = labelFirstDau+iDau;
743  if(indDau<0) return -1;
744  TParticle* dau=stack->Particle(indDau);
745  if(!dau) return -1;
746  Int_t pdgdau=dau->GetPdgCode();
747  if(TMath::Abs(pdgdau)==321){
748  if(pdgD>0 && pdgdau>0) return -1;
749  if(pdgD<0 && pdgdau<0) return -1;
750  nKaons++;
751  sumPxDau+=dau->Px();
752  sumPyDau+=dau->Py();
753  sumPzDau+=dau->Pz();
754  arrayDauLab[nFoundKpi++]=indDau;
755  if(nFoundKpi>4) return -1;
756  }else if(TMath::Abs(pdgdau)==211){
757  nPions++;
758  sumPxDau+=dau->Px();
759  sumPyDau+=dau->Py();
760  sumPzDau+=dau->Pz();
761  arrayDauLab[nFoundKpi++]=indDau;
762  if(nFoundKpi>4) return -1;
763  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
764  Int_t nResDau=dau->GetNDaughters();
765  if(nResDau!=2) return -1;
766  Int_t indFirstResDau=dau->GetDaughter(0);
767  for(Int_t resDau=0; resDau<2; resDau++){
768  Int_t indResDau=indFirstResDau+resDau;
769  if(indResDau<0) return -1;
770  TParticle* resdau=stack->Particle(indResDau);
771  if(!resdau) return -1;
772  Int_t pdgresdau=resdau->GetPdgCode();
773  if(TMath::Abs(pdgresdau)==321){
774  if(pdgD>0 && pdgresdau>0) return -1;
775  if(pdgD<0 && pdgresdau<0) return -1;
776  nKaons++;
777  sumPxDau+=resdau->Px();
778  sumPyDau+=resdau->Py();
779  sumPzDau+=resdau->Pz();
780  arrayDauLab[nFoundKpi++]=indResDau;
781  if(nFoundKpi>4) return -1;
782  }
783  if(TMath::Abs(pdgresdau)==211){
784  nPions++;
785  sumPxDau+=resdau->Px();
786  sumPyDau+=resdau->Py();
787  sumPzDau+=resdau->Pz();
788  arrayDauLab[nFoundKpi++]=indResDau;
789  if(nFoundKpi>4) return -1;
790  }
791  }
792  }else{
793  return -1;
794  }
795  }
796  if(nPions!=3) return -1;
797  if(nKaons!=1) return -1;
798  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
799  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
800  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
801  return 2;
802  }
803  return -1;
804 }
805 
806 //____________________________________________________________________________
807 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
809 
810  Int_t pdgD=mcPart->GetPdgCode();
811  if(TMath::Abs(pdgD)!=421) return -1;
812 
813  Int_t nDau=mcPart->GetNDaughters();
814 
815  if(nDau==2){
816  Int_t daughter0 = mcPart->GetDaughter(0);
817  Int_t daughter1 = mcPart->GetDaughter(1);
818  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
819  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
820  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
821  arrayDauLab[0]=daughter0;
822  arrayDauLab[1]=daughter1;
823  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
824  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
825  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
826  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
827  return -1;
828  }
829  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
830  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
831  if((pdgCode0*pdgCode1)>0) return -1;
832  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
833  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
834  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
835  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
836  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
837  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
838  return 1;
839  }
840 
841  if(nDau==3 || nDau==4){
842  Int_t nKaons=0;
843  Int_t nPions=0;
844  Double_t sumPxDau=0.;
845  Double_t sumPyDau=0.;
846  Double_t sumPzDau=0.;
847  Int_t nFoundKpi=0;
848  Int_t labelFirstDau = mcPart->GetDaughter(0);
849  for(Int_t iDau=0; iDau<nDau; iDau++){
850  Int_t indDau = labelFirstDau+iDau;
851  if(indDau<0) return -1;
852  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
853  if(!dau) return -1;
854  Int_t pdgdau=dau->GetPdgCode();
855  if(TMath::Abs(pdgdau)==321){
856  if(pdgD>0 && pdgdau>0) return -1;
857  if(pdgD<0 && pdgdau<0) return -1;
858  nKaons++;
859  sumPxDau+=dau->Px();
860  sumPyDau+=dau->Py();
861  sumPzDau+=dau->Pz();
862  arrayDauLab[nFoundKpi++]=indDau;
863  if(nFoundKpi>4) return -1;
864  }else if(TMath::Abs(pdgdau)==211){
865  nPions++;
866  sumPxDau+=dau->Px();
867  sumPyDau+=dau->Py();
868  sumPzDau+=dau->Pz();
869  arrayDauLab[nFoundKpi++]=indDau;
870  if(nFoundKpi>4) return -1;
871  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
872  Int_t nResDau=dau->GetNDaughters();
873  if(nResDau!=2) return -1;
874  Int_t indFirstResDau=dau->GetDaughter(0);
875  for(Int_t resDau=0; resDau<2; resDau++){
876  Int_t indResDau=indFirstResDau+resDau;
877  if(indResDau<0) return -1;
878  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
879  if(!resdau) return -1;
880  Int_t pdgresdau=resdau->GetPdgCode();
881  if(TMath::Abs(pdgresdau)==321){
882  if(pdgD>0 && pdgresdau>0) return -1;
883  if(pdgD<0 && pdgresdau<0) return -1;
884  nKaons++;
885  sumPxDau+=resdau->Px();
886  sumPyDau+=resdau->Py();
887  sumPzDau+=resdau->Pz();
888  arrayDauLab[nFoundKpi++]=indResDau;
889  if(nFoundKpi>4) return -1;
890  }
891  if(TMath::Abs(pdgresdau)==211){
892  nPions++;
893  sumPxDau+=resdau->Px();
894  sumPyDau+=resdau->Py();
895  sumPzDau+=resdau->Pz();
896  arrayDauLab[nFoundKpi++]=indResDau;
897  if(nFoundKpi>4) return -1;
898  }
899  }
900  }else{
901  return -1;
902  }
903  }
904  if(nPions!=3) return -1;
905  if(nKaons!=1) return -1;
906  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
907  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
908  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
909  return 2;
910  }
911  return -1;
912 }
913 //____________________________________________________________________________
914 Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
916 
917  if(label<0) return -1;
918  TParticle* mcPart = stack->Particle(label);
919  Int_t pdgD=mcPart->GetPdgCode();
920  if(TMath::Abs(pdgD)!=411) return -1;
921 
922  Int_t nDau=mcPart->GetNDaughters();
923  Int_t labelFirstDau = mcPart->GetDaughter(0);
924  Int_t nKaons=0;
925  Int_t nPions=0;
926  Double_t sumPxDau=0.;
927  Double_t sumPyDau=0.;
928  Double_t sumPzDau=0.;
929  Int_t nFoundKpi=0;
930 
931  if(nDau==3 || nDau==2){
932  for(Int_t iDau=0; iDau<nDau; iDau++){
933  Int_t indDau = labelFirstDau+iDau;
934  if(indDau<0) return -1;
935  TParticle* dau=stack->Particle(indDau);
936  if(!dau) return -1;
937  Int_t pdgdau=dau->GetPdgCode();
938  if(TMath::Abs(pdgdau)==321){
939  if(pdgD*pdgdau>0) return -1;
940  nKaons++;
941  sumPxDau+=dau->Px();
942  sumPyDau+=dau->Py();
943  sumPzDau+=dau->Pz();
944  arrayDauLab[nFoundKpi++]=indDau;
945  if(nFoundKpi>3) return -1;
946  }else if(TMath::Abs(pdgdau)==211){
947  if(pdgD*pdgdau<0) return -1;
948  nPions++;
949  sumPxDau+=dau->Px();
950  sumPyDau+=dau->Py();
951  sumPzDau+=dau->Pz();
952  arrayDauLab[nFoundKpi++]=indDau;
953  if(nFoundKpi>3) return -1;
954  }else if(TMath::Abs(pdgdau)==313){
955  Int_t nResDau=dau->GetNDaughters();
956  if(nResDau!=2) return -1;
957  Int_t indFirstResDau=dau->GetDaughter(0);
958  for(Int_t resDau=0; resDau<2; resDau++){
959  Int_t indResDau=indFirstResDau+resDau;
960  if(indResDau<0) return -1;
961  TParticle* resdau=stack->Particle(indResDau);
962  if(!resdau) return -1;
963  Int_t pdgresdau=resdau->GetPdgCode();
964  if(TMath::Abs(pdgresdau)==321){
965  if(pdgD*pdgresdau>0) return -1;
966  sumPxDau+=resdau->Px();
967  sumPyDau+=resdau->Py();
968  sumPzDau+=resdau->Pz();
969  nKaons++;
970  arrayDauLab[nFoundKpi++]=indResDau;
971  if(nFoundKpi>3) return -1;
972  }
973  if(TMath::Abs(pdgresdau)==211){
974  if(pdgD*pdgresdau<0) return -1;
975  sumPxDau+=resdau->Px();
976  sumPyDau+=resdau->Py();
977  sumPzDau+=resdau->Pz();
978  nPions++;
979  arrayDauLab[nFoundKpi++]=indResDau;
980  if(nFoundKpi>3) return -1;
981  }
982  }
983  }else{
984  return -1;
985  }
986  }
987  if(nPions!=2) return -1;
988  if(nKaons!=1) return -1;
989  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
990  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
991  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
992  if(nDau==3) return 1;
993  else if(nDau==2) return 2;
994  }
995 
996  return -1;
997 
998 }
999 
1000 //____________________________________________________________________________
1001 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1003 
1004  Int_t pdgD=mcPart->GetPdgCode();
1005  if(TMath::Abs(pdgD)!=411) return -1;
1006 
1007  Int_t nDau=mcPart->GetNDaughters();
1008  Int_t labelFirstDau = mcPart->GetDaughter(0);
1009  Int_t nKaons=0;
1010  Int_t nPions=0;
1011  Double_t sumPxDau=0.;
1012  Double_t sumPyDau=0.;
1013  Double_t sumPzDau=0.;
1014  Int_t nFoundKpi=0;
1015 
1016  if(nDau==3 || nDau==2){
1017  for(Int_t iDau=0; iDau<nDau; iDau++){
1018  Int_t indDau = labelFirstDau+iDau;
1019  if(indDau<0) return -1;
1020  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1021  if(!dau) return -1;
1022  Int_t pdgdau=dau->GetPdgCode();
1023  if(TMath::Abs(pdgdau)==321){
1024  if(pdgD*pdgdau>0) return -1;
1025  nKaons++;
1026  sumPxDau+=dau->Px();
1027  sumPyDau+=dau->Py();
1028  sumPzDau+=dau->Pz();
1029  arrayDauLab[nFoundKpi++]=indDau;
1030  if(nFoundKpi>3) return -1;
1031  }else if(TMath::Abs(pdgdau)==211){
1032  if(pdgD*pdgdau<0) return -1;
1033  nPions++;
1034  sumPxDau+=dau->Px();
1035  sumPyDau+=dau->Py();
1036  sumPzDau+=dau->Pz();
1037  arrayDauLab[nFoundKpi++]=indDau;
1038  if(nFoundKpi>3) return -1;
1039  }else if(TMath::Abs(pdgdau)==313){
1040  Int_t nResDau=dau->GetNDaughters();
1041  if(nResDau!=2) return -1;
1042  Int_t indFirstResDau=dau->GetDaughter(0);
1043  for(Int_t resDau=0; resDau<2; resDau++){
1044  Int_t indResDau=indFirstResDau+resDau;
1045  if(indResDau<0) return -1;
1046  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1047  if(!resdau) return -1;
1048  Int_t pdgresdau=resdau->GetPdgCode();
1049  if(TMath::Abs(pdgresdau)==321){
1050  if(pdgD*pdgresdau>0) return -1;
1051  sumPxDau+=resdau->Px();
1052  sumPyDau+=resdau->Py();
1053  sumPzDau+=resdau->Pz();
1054  nKaons++;
1055  arrayDauLab[nFoundKpi++]=indResDau;
1056  if(nFoundKpi>3) return -1;
1057  }
1058  if(TMath::Abs(pdgresdau)==211){
1059  if(pdgD*pdgresdau<0) return -1;
1060  sumPxDau+=resdau->Px();
1061  sumPyDau+=resdau->Py();
1062  sumPzDau+=resdau->Pz();
1063  nPions++;
1064  arrayDauLab[nFoundKpi++]=indResDau;
1065  if(nFoundKpi>3) return -1;
1066  }
1067  }
1068  }else{
1069  return -1;
1070  }
1071  }
1072  if(nPions!=2) return -1;
1073  if(nKaons!=1) return -1;
1074  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1075  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1076  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1077  if(nDau==3) return 1;
1078  else if(nDau==2) return 2;
1079  }
1080 
1081  return -1;
1082 
1083 }
1084 //____________________________________________________________________________
1085 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1087 
1088  if(label<0) return -1;
1089  TParticle* mcPart = stack->Particle(label);
1090  Int_t pdgD=mcPart->GetPdgCode();
1091  if(TMath::Abs(pdgD)!=411) return -1;
1092 
1093  Int_t nDau=mcPart->GetNDaughters();
1094  Int_t labelFirstDau = mcPart->GetDaughter(0);
1095  Int_t nKaons=0;
1096  Int_t nPions=0;
1097  Double_t sumPxDau=0.;
1098  Double_t sumPyDau=0.;
1099  Double_t sumPzDau=0.;
1100  Int_t nFoundKpi=0;
1101  Bool_t isPhi=kFALSE;
1102  Bool_t isk0st=kFALSE;
1103 
1104  if(nDau==3 || nDau==2){
1105  for(Int_t iDau=0; iDau<nDau; iDau++){
1106  Int_t indDau = labelFirstDau+iDau;
1107  if(indDau<0) return -1;
1108  TParticle* dau=stack->Particle(indDau);
1109  if(!dau) return -1;
1110  Int_t pdgdau=dau->GetPdgCode();
1111  if(TMath::Abs(pdgdau)==321){
1112  nKaons++;
1113  sumPxDau+=dau->Px();
1114  sumPyDau+=dau->Py();
1115  sumPzDau+=dau->Pz();
1116  arrayDauLab[nFoundKpi++]=indDau;
1117  if(nFoundKpi>3) return -1;
1118  }else if(TMath::Abs(pdgdau)==211){
1119  nPions++;
1120  sumPxDau+=dau->Px();
1121  sumPyDau+=dau->Py();
1122  sumPzDau+=dau->Pz();
1123  arrayDauLab[nFoundKpi++]=indDau;
1124  if(nFoundKpi>3) return -1;
1125  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1126  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1127  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1128  Int_t nResDau=dau->GetNDaughters();
1129  if(nResDau!=2) return -1;
1130  Int_t indFirstResDau=dau->GetDaughter(0);
1131  for(Int_t resDau=0; resDau<2; resDau++){
1132  Int_t indResDau=indFirstResDau+resDau;
1133  if(indResDau<0) return -1;
1134  TParticle* resdau=stack->Particle(indResDau);
1135  if(!resdau) return -1;
1136  Int_t pdgresdau=resdau->GetPdgCode();
1137  if(TMath::Abs(pdgresdau)==321){
1138  sumPxDau+=resdau->Px();
1139  sumPyDau+=resdau->Py();
1140  sumPzDau+=resdau->Pz();
1141  nKaons++;
1142  arrayDauLab[nFoundKpi++]=indResDau;
1143  if(nFoundKpi>3) return -1;
1144  }
1145  if(TMath::Abs(pdgresdau)==211){
1146  sumPxDau+=resdau->Px();
1147  sumPyDau+=resdau->Py();
1148  sumPzDau+=resdau->Pz();
1149  nPions++;
1150  arrayDauLab[nFoundKpi++]=indResDau;
1151  if(nFoundKpi>3) return -1;
1152  }
1153  }
1154  }else{
1155  return -1;
1156  }
1157  }
1158  if(nPions!=1) return -1;
1159  if(nKaons!=2) return -1;
1160  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1161  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1162  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1163  if(nDau==3) return 3;
1164  else if(nDau==2){
1165  if(isk0st) return 2;
1166  if(isPhi) return 1;
1167  }
1168  }
1169 
1170  return -1;
1171 
1172 }
1173 //____________________________________________________________________________
1174 Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1176 
1177  if(label<0) return -1;
1178  TParticle* mcPart = stack->Particle(label);
1179  Int_t pdgD=mcPart->GetPdgCode();
1180  if(TMath::Abs(pdgD)!=411) return -1;
1181 
1182  Int_t nDau=mcPart->GetNDaughters();
1183  Int_t labelFirstDau = mcPart->GetDaughter(0);
1184  Int_t nPions=0;
1185  Double_t sumPxDau=0.;
1186  Double_t sumPyDau=0.;
1187  Double_t sumPzDau=0.;
1188  Int_t nFoundpi=0;
1189 
1190  Int_t codeV0=-1;
1191  if(nDau==2){
1192  for(Int_t iDau=0; iDau<nDau; iDau++){
1193  Int_t indDau = labelFirstDau+iDau;
1194  if(indDau<0) return -1;
1195  TParticle* dau=stack->Particle(indDau);
1196  if(!dau) return -1;
1197  Int_t pdgdau=dau->GetPdgCode();
1198  if(TMath::Abs(pdgdau)==211){
1199  nPions++;
1200  sumPxDau+=dau->Px();
1201  sumPyDau+=dau->Py();
1202  sumPzDau+=dau->Pz();
1203  arrayDauLab[nFoundpi++]=indDau;
1204  if(nFoundpi>3) return -1;
1205  }else if(TMath::Abs(pdgdau)==311){
1206  codeV0=TMath::Abs(pdgdau);
1207  TParticle* v0=dau;
1208  if(codeV0==311){
1209  Int_t nK0Dau=dau->GetNDaughters();
1210  if(nK0Dau!=1) return -1;
1211  Int_t indK0s=dau->GetDaughter(0);
1212  if(indK0s<0) return -1;
1213  v0=stack->Particle(indK0s);
1214  if(!v0) return -1;
1215  Int_t pdgK0sl=v0->GetPdgCode();
1216  codeV0=TMath::Abs(pdgK0sl);
1217  }
1218  Int_t nV0Dau=v0->GetNDaughters();
1219  if(nV0Dau!=2) return -1;
1220  Int_t indFirstV0Dau=v0->GetDaughter(0);
1221  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1222  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1223  if(indV0Dau<0) return -1;
1224  TParticle* v0dau=stack->Particle(indV0Dau);
1225  if(!v0dau) return -1;
1226  Int_t pdgv0dau=v0dau->GetPdgCode();
1227  if(TMath::Abs(pdgv0dau)==211){
1228  sumPxDau+=v0dau->Px();
1229  sumPyDau+=v0dau->Py();
1230  sumPzDau+=v0dau->Pz();
1231  nPions++;
1232  arrayDauLab[nFoundpi++]=indV0Dau;
1233  if(nFoundpi>3) return -1;
1234  }
1235  }
1236  }else{
1237  return -1;
1238  }
1239  }
1240  if(nPions!=3) return -1;
1241  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1242  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1243  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1244  if(codeV0==310) return 1;
1245  }
1246  return -1;
1247 
1248 }
1249 
1250 
1251 //____________________________________________________________________________
1252 Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1254 
1255  if(label<0) return -1;
1256  TParticle* mcPart = stack->Particle(label);
1257  Int_t pdgD=mcPart->GetPdgCode();
1258  if(TMath::Abs(pdgD)!=431) return -1;
1259 
1260  Int_t nDau=mcPart->GetNDaughters();
1261  Int_t labelFirstDau = mcPart->GetDaughter(0);
1262  Int_t nKaons=0;
1263  Int_t nPions=0;
1264  Double_t sumPxDau=0.;
1265  Double_t sumPyDau=0.;
1266  Double_t sumPzDau=0.;
1267  Int_t nFoundKpi=0;
1268  Bool_t isPhi=kFALSE;
1269  Bool_t isk0st=kFALSE;
1270 
1271  if(nDau==3 || nDau==2){
1272  for(Int_t iDau=0; iDau<nDau; iDau++){
1273  Int_t indDau = labelFirstDau+iDau;
1274  if(indDau<0) return -1;
1275  TParticle* dau=stack->Particle(indDau);
1276  if(!dau) return -1;
1277  Int_t pdgdau=dau->GetPdgCode();
1278  if(TMath::Abs(pdgdau)==321){
1279  nKaons++;
1280  sumPxDau+=dau->Px();
1281  sumPyDau+=dau->Py();
1282  sumPzDau+=dau->Pz();
1283  arrayDauLab[nFoundKpi++]=indDau;
1284  if(nFoundKpi>3) return -1;
1285  }else if(TMath::Abs(pdgdau)==211){
1286  nPions++;
1287  sumPxDau+=dau->Px();
1288  sumPyDau+=dau->Py();
1289  sumPzDau+=dau->Pz();
1290  arrayDauLab[nFoundKpi++]=indDau;
1291  if(nFoundKpi>3) return -1;
1292  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1293  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1294  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1295  Int_t nResDau=dau->GetNDaughters();
1296  if(nResDau!=2) return -1;
1297  Int_t indFirstResDau=dau->GetDaughter(0);
1298  for(Int_t resDau=0; resDau<2; resDau++){
1299  Int_t indResDau=indFirstResDau+resDau;
1300  if(indResDau<0) return -1;
1301  TParticle* resdau=stack->Particle(indResDau);
1302  if(!resdau) return -1;
1303  Int_t pdgresdau=resdau->GetPdgCode();
1304  if(TMath::Abs(pdgresdau)==321){
1305  sumPxDau+=resdau->Px();
1306  sumPyDau+=resdau->Py();
1307  sumPzDau+=resdau->Pz();
1308  nKaons++;
1309  arrayDauLab[nFoundKpi++]=indResDau;
1310  if(nFoundKpi>3) return -1;
1311  }
1312  if(TMath::Abs(pdgresdau)==211){
1313  sumPxDau+=resdau->Px();
1314  sumPyDau+=resdau->Py();
1315  sumPzDau+=resdau->Pz();
1316  nPions++;
1317  arrayDauLab[nFoundKpi++]=indResDau;
1318  if(nFoundKpi>3) return -1;
1319  }
1320  }
1321  }else{
1322  return -1;
1323  }
1324  }
1325  if(nPions!=1) return -1;
1326  if(nKaons!=2) return -1;
1327  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1328  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1329  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1330  if(nDau==3) return 3;
1331  else if(nDau==2){
1332  if(isk0st) return 2;
1333  if(isPhi) return 1;
1334  }
1335  }
1336 
1337  return -1;
1338 
1339 }
1340 //____________________________________________________________________________
1341 Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1343 
1344  if(label<0) return -1;
1345  TParticle* mcPart = stack->Particle(label);
1346  Int_t pdgD=mcPart->GetPdgCode();
1347  if(TMath::Abs(pdgD)!=431) return -1;
1348 
1349  Int_t nDau=mcPart->GetNDaughters();
1350  Int_t labelFirstDau = mcPart->GetDaughter(0);
1351  Int_t nPions=0;
1352  Int_t nKaons=0;
1353  Double_t sumPxDau=0.;
1354  Double_t sumPyDau=0.;
1355  Double_t sumPzDau=0.;
1356  Int_t nFoundKpi=0;
1357 
1358  Int_t codeV0=-1;
1359  if(nDau==2){
1360  for(Int_t iDau=0; iDau<nDau; iDau++){
1361  Int_t indDau = labelFirstDau+iDau;
1362  if(indDau<0) return -1;
1363  TParticle* dau=stack->Particle(indDau);
1364  if(!dau) return -1;
1365  Int_t pdgdau=dau->GetPdgCode();
1366  if(TMath::Abs(pdgdau)==211){
1367  nPions++;
1368  sumPxDau+=dau->Px();
1369  sumPyDau+=dau->Py();
1370  sumPzDau+=dau->Pz();
1371  arrayDauLab[nFoundKpi++]=indDau;
1372  if(nFoundKpi>3) return -1;
1373  }else if(TMath::Abs(pdgdau)==321){
1374  nKaons++;
1375  sumPxDau+=dau->Px();
1376  sumPyDau+=dau->Py();
1377  sumPzDau+=dau->Pz();
1378  arrayDauLab[nFoundKpi++]=indDau;
1379  if(nFoundKpi>3) return -1;
1380  }else if(TMath::Abs(pdgdau)==311){
1381  codeV0=TMath::Abs(pdgdau);
1382  TParticle* v0=dau;
1383  if(codeV0==311){
1384  Int_t nK0Dau=dau->GetNDaughters();
1385  if(nK0Dau!=1) return -1;
1386  Int_t indK0s=dau->GetDaughter(0);
1387  if(indK0s<0) return -1;
1388  v0=stack->Particle(indK0s);
1389  if(!v0) return -1;
1390  Int_t pdgK0sl=v0->GetPdgCode();
1391  codeV0=TMath::Abs(pdgK0sl);
1392  }
1393  Int_t nV0Dau=v0->GetNDaughters();
1394  if(nV0Dau!=2) return -1;
1395  Int_t indFirstV0Dau=v0->GetDaughter(0);
1396  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1397  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1398  if(indV0Dau<0) return -1;
1399  TParticle* v0dau=stack->Particle(indV0Dau);
1400  if(!v0dau) return -1;
1401  Int_t pdgv0dau=v0dau->GetPdgCode();
1402  if(TMath::Abs(pdgv0dau)==211){
1403  sumPxDau+=v0dau->Px();
1404  sumPyDau+=v0dau->Py();
1405  sumPzDau+=v0dau->Pz();
1406  nPions++;
1407  arrayDauLab[nFoundKpi++]=indV0Dau;
1408  if(nFoundKpi>3) return -1;
1409  }
1410  }
1411  }else{
1412  return -1;
1413  }
1414  }
1415  if(nPions!=2) return -1;
1416  if(nKaons!=1) return -1;
1417  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1418  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1419  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1420  if(codeV0==310) return 1;
1421  }
1422  return -1;
1423 
1424 }
1425 
1426 
1427 //____________________________________________________________________________
1428 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1430 
1431  Int_t pdgD=mcPart->GetPdgCode();
1432  if(TMath::Abs(pdgD)!=431) return -1;
1433 
1434  Int_t nDau=mcPart->GetNDaughters();
1435  Int_t labelFirstDau = mcPart->GetDaughter(0);
1436  Int_t nKaons=0;
1437  Int_t nPions=0;
1438  Double_t sumPxDau=0.;
1439  Double_t sumPyDau=0.;
1440  Double_t sumPzDau=0.;
1441  Int_t nFoundKpi=0;
1442  Bool_t isPhi=kFALSE;
1443  Bool_t isk0st=kFALSE;
1444 
1445  if(nDau==3 || nDau==2){
1446  for(Int_t iDau=0; iDau<nDau; iDau++){
1447  Int_t indDau = labelFirstDau+iDau;
1448  if(indDau<0) return -1;
1449  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1450  if(!dau) return -1;
1451  Int_t pdgdau=dau->GetPdgCode();
1452  if(TMath::Abs(pdgdau)==321){
1453  nKaons++;
1454  sumPxDau+=dau->Px();
1455  sumPyDau+=dau->Py();
1456  sumPzDau+=dau->Pz();
1457  arrayDauLab[nFoundKpi++]=indDau;
1458  if(nFoundKpi>3) return -1;
1459  }else if(TMath::Abs(pdgdau)==211){
1460  nPions++;
1461  sumPxDau+=dau->Px();
1462  sumPyDau+=dau->Py();
1463  sumPzDau+=dau->Pz();
1464  arrayDauLab[nFoundKpi++]=indDau;
1465  if(nFoundKpi>3) return -1;
1466  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1467  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1468  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1469  Int_t nResDau=dau->GetNDaughters();
1470  if(nResDau!=2) return -1;
1471  Int_t indFirstResDau=dau->GetDaughter(0);
1472  for(Int_t resDau=0; resDau<2; resDau++){
1473  Int_t indResDau=indFirstResDau+resDau;
1474  if(indResDau<0) return -1;
1475  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1476  if(!resdau) return -1;
1477  Int_t pdgresdau=resdau->GetPdgCode();
1478  if(TMath::Abs(pdgresdau)==321){
1479  sumPxDau+=resdau->Px();
1480  sumPyDau+=resdau->Py();
1481  sumPzDau+=resdau->Pz();
1482  nKaons++;
1483  arrayDauLab[nFoundKpi++]=indResDau;
1484  if(nFoundKpi>3) return -1;
1485  }
1486  if(TMath::Abs(pdgresdau)==211){
1487  sumPxDau+=resdau->Px();
1488  sumPyDau+=resdau->Py();
1489  sumPzDau+=resdau->Pz();
1490  nPions++;
1491  arrayDauLab[nFoundKpi++]=indResDau;
1492  if(nFoundKpi>3) return -1;
1493  }
1494  }
1495  }else{
1496  return -1;
1497  }
1498  }
1499  if(nPions!=1) return -1;
1500  if(nKaons!=2) return -1;
1501  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1502  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1503  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1504  if(nDau==3) return 3;
1505  else if(nDau==2){
1506  if(isk0st) return 2;
1507  if(isPhi) return 1;
1508  }
1509  }
1510 
1511  return -1;
1512 
1513 }
1514 //____________________________________________________________________________
1515 Int_t AliVertexingHFUtils::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1517 
1518  if(label<0) return -1;
1519  TParticle* mcPart = stack->Particle(label);
1520  Int_t pdgD=mcPart->GetPdgCode();
1521  if(TMath::Abs(pdgD)!=413) return -1;
1522 
1523  Int_t nDau=mcPart->GetNDaughters();
1524  if(nDau!=2) return -1;
1525 
1526  Int_t labelFirstDau = mcPart->GetDaughter(0);
1527  Int_t nKaons=0;
1528  Int_t nPions=0;
1529  Double_t sumPxDau=0.;
1530  Double_t sumPyDau=0.;
1531  Double_t sumPzDau=0.;
1532  Int_t nFoundKpi=0;
1533 
1534  for(Int_t iDau=0; iDau<nDau; iDau++){
1535  Int_t indDau = labelFirstDau+iDau;
1536  if(indDau<0) return -1;
1537  TParticle* dau=stack->Particle(indDau);
1538  if(!dau) return -1;
1539  Int_t pdgdau=dau->GetPdgCode();
1540  if(TMath::Abs(pdgdau)==421){
1541  Int_t nResDau=dau->GetNDaughters();
1542  if(nResDau!=2) return -1;
1543  Int_t indFirstResDau=dau->GetDaughter(0);
1544  for(Int_t resDau=0; resDau<2; resDau++){
1545  Int_t indResDau=indFirstResDau+resDau;
1546  if(indResDau<0) return -1;
1547  TParticle* resdau=stack->Particle(indResDau);
1548  if(!resdau) return -1;
1549  Int_t pdgresdau=resdau->GetPdgCode();
1550  if(TMath::Abs(pdgresdau)==321){
1551  if(pdgD*pdgresdau>0) return -1;
1552  sumPxDau+=resdau->Px();
1553  sumPyDau+=resdau->Py();
1554  sumPzDau+=resdau->Pz();
1555  nKaons++;
1556  arrayDauLab[nFoundKpi++]=indResDau;
1557  if(nFoundKpi>3) return -1;
1558  }
1559  if(TMath::Abs(pdgresdau)==211){
1560  if(pdgD*pdgresdau<0) return -1;
1561  sumPxDau+=resdau->Px();
1562  sumPyDau+=resdau->Py();
1563  sumPzDau+=resdau->Pz();
1564  nPions++;
1565  arrayDauLab[nFoundKpi++]=indResDau;
1566  if(nFoundKpi>3) return -1;
1567  }
1568  }
1569  }else if(TMath::Abs(pdgdau)==211){
1570  if(pdgD*pdgdau<0) return -1;
1571  nPions++;
1572  sumPxDau+=dau->Px();
1573  sumPyDau+=dau->Py();
1574  sumPzDau+=dau->Pz();
1575  arrayDauLab[nFoundKpi++]=indDau;
1576  if(nFoundKpi>3) return -1;
1577  }
1578  }
1579 
1580  if(nPions!=2) return -1;
1581  if(nKaons!=1) return -1;
1582  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1583  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1584  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1585  return 1;
1586 
1587 }
1588 
1589 //____________________________________________________________________________
1590 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1592 
1593  Int_t pdgD=mcPart->GetPdgCode();
1594  if(TMath::Abs(pdgD)!=413) return -1;
1595 
1596  Int_t nDau=mcPart->GetNDaughters();
1597  if(nDau!=2) return -1;
1598 
1599  Int_t labelFirstDau = mcPart->GetDaughter(0);
1600  Int_t nKaons=0;
1601  Int_t nPions=0;
1602  Double_t sumPxDau=0.;
1603  Double_t sumPyDau=0.;
1604  Double_t sumPzDau=0.;
1605  Int_t nFoundKpi=0;
1606 
1607  for(Int_t iDau=0; iDau<nDau; iDau++){
1608  Int_t indDau = labelFirstDau+iDau;
1609  if(indDau<0) return -1;
1610  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1611  if(!dau) return -1;
1612  Int_t pdgdau=dau->GetPdgCode();
1613  if(TMath::Abs(pdgdau)==421){
1614  Int_t nResDau=dau->GetNDaughters();
1615  if(nResDau!=2) return -1;
1616  Int_t indFirstResDau=dau->GetDaughter(0);
1617  for(Int_t resDau=0; resDau<2; resDau++){
1618  Int_t indResDau=indFirstResDau+resDau;
1619  if(indResDau<0) return -1;
1620  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1621  if(!resdau) return -1;
1622  Int_t pdgresdau=resdau->GetPdgCode();
1623  if(TMath::Abs(pdgresdau)==321){
1624  if(pdgD*pdgresdau>0) return -1;
1625  sumPxDau+=resdau->Px();
1626  sumPyDau+=resdau->Py();
1627  sumPzDau+=resdau->Pz();
1628  nKaons++;
1629  arrayDauLab[nFoundKpi++]=indResDau;
1630  if(nFoundKpi>3) return -1;
1631  }
1632  if(TMath::Abs(pdgresdau)==211){
1633  if(pdgD*pdgresdau<0) return -1;
1634  sumPxDau+=resdau->Px();
1635  sumPyDau+=resdau->Py();
1636  sumPzDau+=resdau->Pz();
1637  nPions++;
1638  arrayDauLab[nFoundKpi++]=indResDau;
1639  if(nFoundKpi>3) return -1;
1640  }
1641  }
1642  }else if(TMath::Abs(pdgdau)==211){
1643  if(pdgD*pdgdau<0) return -1;
1644  nPions++;
1645  sumPxDau+=dau->Px();
1646  sumPyDau+=dau->Py();
1647  sumPzDau+=dau->Pz();
1648  arrayDauLab[nFoundKpi++]=indDau;
1649  if(nFoundKpi>3) return -1;
1650  }
1651  }
1652 
1653  if(nPions!=2) return -1;
1654  if(nKaons!=1) return -1;
1655  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1656  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1657  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1658  return 1;
1659 
1660 }
1661 
1662 //____________________________________________________________________________
1663 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1665 
1666  if(label<0) return -1;
1667  TParticle* mcPart = stack->Particle(label);
1668  Int_t pdgD=mcPart->GetPdgCode();
1669  if(TMath::Abs(pdgD)!=4122) return -1;
1670 
1671  Int_t nDau=mcPart->GetNDaughters();
1672  Int_t labelFirstDau = mcPart->GetDaughter(0);
1673  Int_t nKaons=0;
1674  Int_t nPions=0;
1675  Int_t nProtons=0;
1676  Double_t sumPxDau=0.;
1677  Double_t sumPyDau=0.;
1678  Double_t sumPzDau=0.;
1679  Int_t nFoundpKpi=0;
1680 
1681  Int_t codeRes=-1;
1682  if(nDau==3 || nDau==2){
1683  for(Int_t iDau=0; iDau<nDau; iDau++){
1684  Int_t indDau = labelFirstDau+iDau;
1685  if(indDau<0) return -1;
1686  TParticle* dau=stack->Particle(indDau);
1687  if(!dau) return -1;
1688  Int_t pdgdau=dau->GetPdgCode();
1689  if(TMath::Abs(pdgdau)==321){
1690  nKaons++;
1691  sumPxDau+=dau->Px();
1692  sumPyDau+=dau->Py();
1693  sumPzDau+=dau->Pz();
1694  arrayDauLab[nFoundpKpi++]=indDau;
1695  if(nFoundpKpi>3) return -1;
1696  }else if(TMath::Abs(pdgdau)==211){
1697  nPions++;
1698  sumPxDau+=dau->Px();
1699  sumPyDau+=dau->Py();
1700  sumPzDau+=dau->Pz();
1701  arrayDauLab[nFoundpKpi++]=indDau;
1702  if(nFoundpKpi>3) return -1;
1703  }else if(TMath::Abs(pdgdau)==2212){
1704  nProtons++;
1705  sumPxDau+=dau->Px();
1706  sumPyDau+=dau->Py();
1707  sumPzDau+=dau->Pz();
1708  arrayDauLab[nFoundpKpi++]=indDau;
1709  if(nFoundpKpi>3) return -1;
1710  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1711  TMath::Abs(pdgdau)==2224){
1712  codeRes=TMath::Abs(pdgdau);
1713  Int_t nResDau=dau->GetNDaughters();
1714  if(nResDau!=2) return -1;
1715  Int_t indFirstResDau=dau->GetDaughter(0);
1716  for(Int_t resDau=0; resDau<2; resDau++){
1717  Int_t indResDau=indFirstResDau+resDau;
1718  if(indResDau<0) return -1;
1719  TParticle* resdau=stack->Particle(indResDau);
1720  if(!resdau) return -1;
1721  Int_t pdgresdau=resdau->GetPdgCode();
1722  if(TMath::Abs(pdgresdau)==321){
1723  sumPxDau+=resdau->Px();
1724  sumPyDau+=resdau->Py();
1725  sumPzDau+=resdau->Pz();
1726  nKaons++;
1727  arrayDauLab[nFoundpKpi++]=indResDau;
1728  if(nFoundpKpi>3) return -1;
1729  }else if(TMath::Abs(pdgresdau)==211){
1730  sumPxDau+=resdau->Px();
1731  sumPyDau+=resdau->Py();
1732  sumPzDau+=resdau->Pz();
1733  nPions++;
1734  arrayDauLab[nFoundpKpi++]=indResDau;
1735  if(nFoundpKpi>3) return -1;
1736  }else if(TMath::Abs(pdgresdau)==2212){
1737  sumPxDau+=resdau->Px();
1738  sumPyDau+=resdau->Py();
1739  sumPzDau+=resdau->Pz();
1740  nProtons++;
1741  arrayDauLab[nFoundpKpi++]=indResDau;
1742  if(nFoundpKpi>3) return -1;
1743  }
1744  }
1745  }else{
1746  return -1;
1747  }
1748  }
1749  if(nPions!=1) return -1;
1750  if(nKaons!=1) return -1;
1751  if(nProtons!=1) return -1;
1752  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1753  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1754  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1755  if(nDau==3) return 1;
1756  else if(nDau==2){
1757  if(codeRes==313) return 2;
1758  else if(codeRes==2224) return 3;
1759  else if(codeRes==3124) return 4;
1760  }
1761  }
1762  return -1;
1763 
1764 }
1765 
1766 
1767 //____________________________________________________________________________
1768 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1770 
1771  if(label<0) return -1;
1772  TParticle* mcPart = stack->Particle(label);
1773  Int_t pdgD=mcPart->GetPdgCode();
1774  if(TMath::Abs(pdgD)!=4122) return -1;
1775 
1776  Int_t nDau=mcPart->GetNDaughters();
1777  Int_t labelFirstDau = mcPart->GetDaughter(0);
1778  Int_t nPions=0;
1779  Int_t nProtons=0;
1780  Double_t sumPxDau=0.;
1781  Double_t sumPyDau=0.;
1782  Double_t sumPzDau=0.;
1783  Int_t nFoundppi=0;
1784 
1785  Int_t codeV0=-1;
1786  if(nDau==2){
1787  for(Int_t iDau=0; iDau<nDau; iDau++){
1788  Int_t indDau = labelFirstDau+iDau;
1789  if(indDau<0) return -1;
1790  TParticle* dau=stack->Particle(indDau);
1791  if(!dau) return -1;
1792  Int_t pdgdau=dau->GetPdgCode();
1793  if(TMath::Abs(pdgdau)==211){
1794  nPions++;
1795  sumPxDau+=dau->Px();
1796  sumPyDau+=dau->Py();
1797  sumPzDau+=dau->Pz();
1798  arrayDauLab[nFoundppi++]=indDau;
1799  if(nFoundppi>3) return -1;
1800  }else if(TMath::Abs(pdgdau)==2212){
1801  nProtons++;
1802  sumPxDau+=dau->Px();
1803  sumPyDau+=dau->Py();
1804  sumPzDau+=dau->Pz();
1805  arrayDauLab[nFoundppi++]=indDau;
1806  if(nFoundppi>3) return -1;
1807  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
1808  codeV0=TMath::Abs(pdgdau);
1809  TParticle* v0=dau;
1810  if(codeV0==311){
1811  Int_t nK0Dau=dau->GetNDaughters();
1812  if(nK0Dau!=1) return -1;
1813  Int_t indK0s=dau->GetDaughter(0);
1814  if(indK0s<0) return -1;
1815  v0=stack->Particle(indK0s);
1816  if(!v0) return -1;
1817  Int_t pdgK0sl=v0->GetPdgCode();
1818  codeV0=TMath::Abs(pdgK0sl);
1819  }
1820  Int_t nV0Dau=v0->GetNDaughters();
1821  if(nV0Dau!=2) return -1;
1822  Int_t indFirstV0Dau=v0->GetDaughter(0);
1823  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1824  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1825  if(indV0Dau<0) return -1;
1826  TParticle* v0dau=stack->Particle(indV0Dau);
1827  if(!v0dau) return -1;
1828  Int_t pdgv0dau=v0dau->GetPdgCode();
1829  if(TMath::Abs(pdgv0dau)==211){
1830  sumPxDau+=v0dau->Px();
1831  sumPyDau+=v0dau->Py();
1832  sumPzDau+=v0dau->Pz();
1833  nPions++;
1834  arrayDauLab[nFoundppi++]=indV0Dau;
1835  if(nFoundppi>3) return -1;
1836  }else if(TMath::Abs(pdgv0dau)==2212){
1837  sumPxDau+=v0dau->Px();
1838  sumPyDau+=v0dau->Py();
1839  sumPzDau+=v0dau->Pz();
1840  nProtons++;
1841  arrayDauLab[nFoundppi++]=indV0Dau;
1842  if(nFoundppi>3) return -1;
1843  }
1844  }
1845  }else{
1846  return -1;
1847  }
1848  }
1849  if(nPions!=2) return -1;
1850  if(nProtons!=1) return -1;
1851  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1852  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1853  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1854  if(codeV0==310) return 1;
1855  else if(codeV0==3122) return 2;
1856  }
1857  return -1;
1858 
1859 }
1860 
1861 
1862 //__________________________________xic______________________________________
1863 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
1865 
1866  if(label<0) return -1;
1867  TParticle* mcPart = stack->Particle(label);
1868  Int_t pdgD=mcPart->GetPdgCode();
1869  if(TMath::Abs(pdgD)!=4232) return -1;
1870 
1871  Int_t nDau=mcPart->GetNDaughters();
1872  if(nDau!=3) return -1;
1873 
1874  Int_t labelFirstDau = mcPart->GetDaughter(0);
1875  Int_t nXi=0;
1876  Int_t nPions=0;
1877  Double_t sumPxDau=0.;
1878  Double_t sumPyDau=0.;
1879  Double_t sumPzDau=0.;
1880  Int_t nFoundXi=0;
1881 
1882  for(Int_t iDau=0; iDau<nDau; iDau++){
1883  Int_t indDau = labelFirstDau+iDau;
1884  if(indDau<0) return -1;
1885  TParticle* dau=stack->Particle(indDau);
1886  if(!dau) return -1;
1887  Int_t pdgdau=dau->GetPdgCode();
1888  if(TMath::Abs(pdgdau)==3312){
1889  if(pdgD*pdgdau<0) return -1;
1890  sumPxDau+=dau->Px();
1891  sumPyDau+=dau->Py();
1892  sumPzDau+=dau->Pz();
1893  nXi++;
1894  arrayDauLab[nFoundXi++]=indDau;
1895 
1896  }
1897  if(TMath::Abs(pdgdau)==211){
1898  if(pdgD*pdgdau<0) return -1;
1899  nPions++;
1900  sumPxDau+=dau->Px();
1901  sumPyDau+=dau->Py();
1902  sumPzDau+=dau->Pz();
1903  arrayDauLab[nFoundXi++]=indDau;
1904  if(nFoundXi>3) return -1;
1905  }
1906  }
1907 
1908  if(nPions!=2) return -1;
1909  if(nXi!=1) return -1;
1910  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1911  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1912  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1913  return 1;
1914 
1915 }
1916 //________________________________________________________________________
1917 Double_t AliVertexingHFUtils::GetSphericity(AliAODEvent* aod, Double_t etaMin, Double_t etaMax,
1918  Double_t ptMin, Double_t ptMax,
1919  Int_t filtbit1, Int_t filtbit2,
1920  Int_t minMult){
1922 
1923  Int_t nTracks=aod->GetNumberOfTracks();
1924  Int_t nSelTracks=0;
1925 
1926  Double_t sumpt=0.;
1927  Double_t s00=0.;
1928  Double_t s01=0.;
1929  Double_t s11=0.;
1930  if(ptMin<0.) ptMin=0.;
1931 
1932  for(Int_t it=0; it<nTracks; it++) {
1933  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1934  if(!tr) continue;
1935  Float_t eta = tr->Eta();
1936  Float_t pt = tr->Pt();
1937  Float_t phi = tr->Phi();
1938  if(eta<etaMin || eta>etaMax) continue;
1939  if(pt<ptMin || pt>ptMax) continue;
1940  Bool_t fb1 = tr->TestFilterBit(filtbit1);
1941  Bool_t fb2 = tr->TestFilterBit(filtbit2);
1942  if( !(fb1 || fb2) ) continue;
1943  Double_t px=pt*TMath::Cos(phi);
1944  Double_t py=pt*TMath::Sin(phi);
1945  s00 += (px * px)/pt;
1946  s01 += (py * px)/pt;
1947  s11 += (py * py)/pt;
1948  nSelTracks++;
1949  sumpt+=pt;
1950  }
1951 
1952  if(nSelTracks<minMult) return -0.5;
1953 
1954  if(sumpt>0.){
1955  s00/=sumpt;
1956  s01/=sumpt;
1957  s11/=sumpt;
1958  }else return -0.5;
1959 
1960  Double_t sphericity = -10;
1961  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
1962  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
1963  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
1964  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
1965  return sphericity;
1966 
1967 }
1968 
1969 //________________________________________________________________________
1970 Double_t AliVertexingHFUtils::GetSpherocity(AliAODEvent* aod, Double_t etaMin, Double_t etaMax,
1971  Double_t ptMin, Double_t ptMax,
1972  Int_t filtbit1, Int_t filtbit2,
1973  Int_t minMult, Double_t phiStepSizeDeg){
1975 
1976  Int_t nTracks=aod->GetNumberOfTracks();
1977  Int_t nSelTracks=0;
1978 
1979  Double_t* ptArr=new Double_t[nTracks];
1980  Double_t* phiArr=new Double_t[nTracks];
1981  Double_t sumpt=0.;
1982 
1983  for(Int_t it=0; it<nTracks; it++) {
1984  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1985  if(!tr) continue;
1986  Float_t eta = tr->Eta();
1987  Float_t pt = tr->Pt();
1988  Float_t phi = tr->Phi();
1989  if(eta<etaMin || eta>etaMax) continue;
1990  if(pt<ptMin || pt>ptMax) continue;
1991  Bool_t fb1 = tr->TestFilterBit(filtbit1);
1992  Bool_t fb2 = tr->TestFilterBit(filtbit2);
1993  if( !(fb1 || fb2) ) continue;
1994  ptArr[nSelTracks]=pt;
1995  phiArr[nSelTracks]=phi;
1996  nSelTracks++;
1997  sumpt+=pt;
1998  }
1999 
2000  if(nSelTracks<minMult) return -0.5;
2001 
2002  //Getting thrust
2003  Double_t spherocity=2.;
2004  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2005  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2006  Double_t nx=TMath::Cos(phistep);
2007  Double_t ny=TMath::Sin(phistep);
2008  Double_t numer=0.;
2009  for(Int_t j=0; j<nSelTracks; ++j){
2010  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2011  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2012  numer+=TMath::Abs(ny*pxA - nx*pyA);
2013  }
2014  Double_t pFull=numer*numer/(sumpt*sumpt);
2015  if(pFull<spherocity) spherocity=pFull; // minimization;
2016  }
2017 
2018  delete [] ptArr;
2019  delete [] phiArr;
2020 
2021  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2022  return spherocity;
2023 
2024 }
Int_t charge
static void 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=2.5, Float_t sigmaRangeForBkg=4.5, Float_t minMass=0., Float_t maxMass=3., Int_t rebin=1)
Functions for computing average pt.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
static Int_t CheckDplusDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Int_t CheckDstarDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
AliAODv0 * Getv0() const
Double_t ImpParXY() const
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
static Int_t CheckDsDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Double_t Pol(Double_t x) const
static Int_t CheckLcpKpiDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Double_t GetFullEvResolLowLim(const TH1F *hSubEvCorr, Int_t k=1)
Double_t ptMin
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
TRandom * gRandom
void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
static Int_t CheckXicXipipiDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
static Int_t CheckLcV0bachelorDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Double_t GetFullEvResolHighLim(const TH1F *hSubEvCorr, Int_t k=1)
static Double_t ResolK1(Double_t x)
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
static Double_t GetSpherocity(AliAODEvent *aod, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3, Double_t phiStepSizeDeg=0.1)
Functions for event shape variables.
Double_t GetFullEvResol() const
const Double_t ptmax
static Bool_t CheckT0TriggerFired(AliAODEvent *aodEv)
Functions for processing trigger information.
static Int_t CheckD0Decay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Int_t CheckDsK0sKDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
AliAODTrack * GetBachelor() const
static Double_t GetSubEvResolLowLim(const TH1F *hSubEvCorr)
const Double_t ptmin
static Int_t CheckDplusKKpiDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
Bool_t IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
static Int_t CheckDplusK0spiDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
static Double_t GetBeautyMotherPt(TClonesArray *arrayMC, AliAODMCParticle *mcPart)
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
Double_t GetSubEvResol() const
Bool_t IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC)
Double_t minMass
Double_t FindChi() const
Int_t rebin
static Int_t GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
static Double_t GetVZEROAEqualizedMultiplicity(AliAODEvent *ev)
Utilities for V0 multiplicity checks.
static Double_t GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
Functions for computing true impact parameter of D meson.
Double_t maxMass
Double_t ptMax
static Double_t GetSubEvResolHighLim(const TH1F *hSubEvCorr)
static Double_t GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
static Double_t GetSphericity(AliAODEvent *aod, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3)
static Int_t GetGeneratedMultiplicityInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
static Int_t GetGeneratedPrimariesInEtaRange(TClonesArray *arrayMC, Double_t mineta, Double_t maxeta)
Class with functions useful for different D2H analyses //.