AliPhysics  09a22ae (09a22ae)
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 "AliMCEvent.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 
45 ClassImp(AliVertexingHFUtils);
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.),
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 //______________________________________________________________________
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 //______________________________________________________________________
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 //______________________________________________________________________
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 //______________________________________________________________________
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 //______________________________________________________________________
169  if(!hSubEvCorr) return 1.;
170  Double_t resSub=GetSubEvResol(hSubEvCorr);
171  return GetFullEvResol(resSub,k);
172 }
173 //______________________________________________________________________
176  if(!hSubEvCorr) return 1.;
177  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
178  printf("%f\n",resSub);
179  return GetFullEvResol(resSub,k);
180 }
181 //______________________________________________________________________
184  if(!hSubEvCorr) return 1.;
185  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
186  printf("%f\n",resSub);
187  return GetFullEvResol(resSub,k);
188 }
189 //______________________________________________________________________
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 //______________________________________________________________________
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 //______________________________________________________________________
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 //______________________________________________________________________
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 //____________________________________________________________________________
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(AliMCEvent* mcEvent, 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 = mcEvent->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(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
699 
700  if(label<0) return -1;
701  TParticle* mcPart = mcEvent->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 = mcEvent->Particle(daughter0);
711  TParticle* mcPartDaughter1 = mcEvent->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=mcEvent->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=mcEvent->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(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
916 
917  if(label<0) return -1;
918  TParticle* mcPart = mcEvent->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=mcEvent->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=mcEvent->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(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1087 
1088  if(label<0) return -1;
1089  TParticle* mcPart = mcEvent->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=mcEvent->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=mcEvent->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::CheckDplusKKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1176 
1177  Int_t pdgD=mcPart->GetPdgCode();
1178  if(TMath::Abs(pdgD)!=411) return -1;
1179 
1180  Int_t nDau=mcPart->GetNDaughters();
1181  Int_t labelFirstDau = mcPart->GetDaughter(0);
1182  Int_t nKaons=0;
1183  Int_t nPions=0;
1184  Double_t sumPxDau=0.;
1185  Double_t sumPyDau=0.;
1186  Double_t sumPzDau=0.;
1187  Int_t nFoundKpi=0;
1188  Bool_t isPhi=kFALSE;
1189  Bool_t isk0st=kFALSE;
1190 
1191  if(nDau==3 || nDau==2){
1192  for(Int_t iDau=0; iDau<nDau; iDau++){
1193  Int_t indDau = labelFirstDau+iDau;
1194  if(indDau<0) return -1;
1195  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1196  if(!dau) return -1;
1197  Int_t pdgdau=dau->GetPdgCode();
1198  if(TMath::Abs(pdgdau)==321){
1199  nKaons++;
1200  sumPxDau+=dau->Px();
1201  sumPyDau+=dau->Py();
1202  sumPzDau+=dau->Pz();
1203  arrayDauLab[nFoundKpi++]=indDau;
1204  if(nFoundKpi>3) return -1;
1205  }else if(TMath::Abs(pdgdau)==211){
1206  nPions++;
1207  sumPxDau+=dau->Px();
1208  sumPyDau+=dau->Py();
1209  sumPzDau+=dau->Pz();
1210  arrayDauLab[nFoundKpi++]=indDau;
1211  if(nFoundKpi>3) return -1;
1212  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1213  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1214  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1215  Int_t nResDau=dau->GetNDaughters();
1216  if(nResDau!=2) return -1;
1217  Int_t indFirstResDau=dau->GetDaughter(0);
1218  for(Int_t resDau=0; resDau<2; resDau++){
1219  Int_t indResDau=indFirstResDau+resDau;
1220  if(indResDau<0) return -1;
1221  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1222  if(!resdau) return -1;
1223  Int_t pdgresdau=resdau->GetPdgCode();
1224  if(TMath::Abs(pdgresdau)==321){
1225  sumPxDau+=resdau->Px();
1226  sumPyDau+=resdau->Py();
1227  sumPzDau+=resdau->Pz();
1228  nKaons++;
1229  arrayDauLab[nFoundKpi++]=indResDau;
1230  if(nFoundKpi>3) return -1;
1231  }
1232  if(TMath::Abs(pdgresdau)==211){
1233  sumPxDau+=resdau->Px();
1234  sumPyDau+=resdau->Py();
1235  sumPzDau+=resdau->Pz();
1236  nPions++;
1237  arrayDauLab[nFoundKpi++]=indResDau;
1238  if(nFoundKpi>3) return -1;
1239  }
1240  }
1241  }else{
1242  return -1;
1243  }
1244  }
1245  if(nPions!=1) return -1;
1246  if(nKaons!=2) return -1;
1247  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1248  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1249  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1250  if(nDau==3) return 3;
1251  else if(nDau==2){
1252  if(isk0st) return 2;
1253  if(isPhi) return 1;
1254  }
1255  }
1256 
1257  return -1;
1258 
1259 }
1260 //____________________________________________________________________________
1261 Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1263 
1264  if(label<0) return -1;
1265  TParticle* mcPart = mcEvent->Particle(label);
1266  Int_t pdgD=mcPart->GetPdgCode();
1267  if(TMath::Abs(pdgD)!=411) return -1;
1268 
1269  Int_t nDau=mcPart->GetNDaughters();
1270  Int_t labelFirstDau = mcPart->GetDaughter(0);
1271  Int_t nPions=0;
1272  Double_t sumPxDau=0.;
1273  Double_t sumPyDau=0.;
1274  Double_t sumPzDau=0.;
1275  Int_t nFoundpi=0;
1276 
1277  Int_t codeV0=-1;
1278  if(nDau==2){
1279  for(Int_t iDau=0; iDau<nDau; iDau++){
1280  Int_t indDau = labelFirstDau+iDau;
1281  if(indDau<0) return -1;
1282  TParticle* dau=mcEvent->Particle(indDau);
1283  if(!dau) return -1;
1284  Int_t pdgdau=dau->GetPdgCode();
1285  if(TMath::Abs(pdgdau)==211){
1286  nPions++;
1287  sumPxDau+=dau->Px();
1288  sumPyDau+=dau->Py();
1289  sumPzDau+=dau->Pz();
1290  arrayDauLab[nFoundpi++]=indDau;
1291  if(nFoundpi>3) return -1;
1292  }else if(TMath::Abs(pdgdau)==311){
1293  codeV0=TMath::Abs(pdgdau);
1294  TParticle* v0=dau;
1295  if(codeV0==311){
1296  Int_t nK0Dau=dau->GetNDaughters();
1297  if(nK0Dau!=1) return -1;
1298  Int_t indK0s=dau->GetDaughter(0);
1299  if(indK0s<0) return -1;
1300  v0=mcEvent->Particle(indK0s);
1301  if(!v0) return -1;
1302  Int_t pdgK0sl=v0->GetPdgCode();
1303  codeV0=TMath::Abs(pdgK0sl);
1304  }
1305  Int_t nV0Dau=v0->GetNDaughters();
1306  if(nV0Dau!=2) return -1;
1307  Int_t indFirstV0Dau=v0->GetDaughter(0);
1308  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1309  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1310  if(indV0Dau<0) return -1;
1311  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1312  if(!v0dau) return -1;
1313  Int_t pdgv0dau=v0dau->GetPdgCode();
1314  if(TMath::Abs(pdgv0dau)==211){
1315  sumPxDau+=v0dau->Px();
1316  sumPyDau+=v0dau->Py();
1317  sumPzDau+=v0dau->Pz();
1318  nPions++;
1319  arrayDauLab[nFoundpi++]=indV0Dau;
1320  if(nFoundpi>3) return -1;
1321  }
1322  }
1323  }else{
1324  return -1;
1325  }
1326  }
1327  if(nPions!=3) return -1;
1328  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1329  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1330  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1331  if(codeV0==310) return 1;
1332  }
1333  return -1;
1334 
1335 }
1336 
1337 
1338 //____________________________________________________________________________
1339 Int_t AliVertexingHFUtils::CheckDsDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1341 
1342  if(label<0) return -1;
1343  TParticle* mcPart = mcEvent->Particle(label);
1344  Int_t pdgD=mcPart->GetPdgCode();
1345  if(TMath::Abs(pdgD)!=431) return -1;
1346 
1347  Int_t nDau=mcPart->GetNDaughters();
1348  Int_t labelFirstDau = mcPart->GetDaughter(0);
1349  Int_t nKaons=0;
1350  Int_t nPions=0;
1351  Double_t sumPxDau=0.;
1352  Double_t sumPyDau=0.;
1353  Double_t sumPzDau=0.;
1354  Int_t nFoundKpi=0;
1355  Bool_t isPhi=kFALSE;
1356  Bool_t isk0st=kFALSE;
1357  Bool_t isf0=kFALSE;
1358 
1359  if(nDau==3 || 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=mcEvent->Particle(indDau);
1364  if(!dau) return -1;
1365  Int_t pdgdau=dau->GetPdgCode();
1366  if(TMath::Abs(pdgdau)==321){
1367  nKaons++;
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)==211){
1374  nPions++;
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)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1381  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1382  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1383  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1384  Int_t nResDau=dau->GetNDaughters();
1385  if(nResDau!=2) return -1;
1386  Int_t indFirstResDau=dau->GetDaughter(0);
1387  for(Int_t resDau=0; resDau<2; resDau++){
1388  Int_t indResDau=indFirstResDau+resDau;
1389  if(indResDau<0) return -1;
1390  TParticle* resdau=mcEvent->Particle(indResDau);
1391  if(!resdau) return -1;
1392  Int_t pdgresdau=resdau->GetPdgCode();
1393  if(TMath::Abs(pdgresdau)==321){
1394  sumPxDau+=resdau->Px();
1395  sumPyDau+=resdau->Py();
1396  sumPzDau+=resdau->Pz();
1397  nKaons++;
1398  arrayDauLab[nFoundKpi++]=indResDau;
1399  if(nFoundKpi>3) return -1;
1400  }
1401  if(TMath::Abs(pdgresdau)==211){
1402  sumPxDau+=resdau->Px();
1403  sumPyDau+=resdau->Py();
1404  sumPzDau+=resdau->Pz();
1405  nPions++;
1406  arrayDauLab[nFoundKpi++]=indResDau;
1407  if(nFoundKpi>3) return -1;
1408  }
1409  }
1410  }else{
1411  return -1;
1412  }
1413  }
1414  if(nPions!=1) return -1;
1415  if(nKaons!=2) return -1;
1416  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1417  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1418  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1419  if(nDau==3) return 3;
1420  else if(nDau==2){
1421  if(isk0st) return 2;
1422  if(isPhi) return 1;
1423  if(isf0) return 4;
1424  }
1425  }
1426 
1427  return -1;
1428 }
1429 //____________________________________________________________________________
1430 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1432 
1433  Int_t pdgD=mcPart->GetPdgCode();
1434  if(TMath::Abs(pdgD)!=431) return -1;
1435 
1436  Int_t nDau=mcPart->GetNDaughters();
1437  Int_t labelFirstDau = mcPart->GetDaughter(0);
1438  Int_t nKaons=0;
1439  Int_t nPions=0;
1440  Double_t sumPxDau=0.;
1441  Double_t sumPyDau=0.;
1442  Double_t sumPzDau=0.;
1443  Int_t nFoundKpi=0;
1444  Bool_t isPhi=kFALSE;
1445  Bool_t isk0st=kFALSE;
1446  Bool_t isf0=kFALSE;
1447 
1448  if(nDau==3 || nDau==2){
1449  for(Int_t iDau=0; iDau<nDau; iDau++){
1450  Int_t indDau = labelFirstDau+iDau;
1451  if(indDau<0) return -1;
1452  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1453  if(!dau) return -1;
1454  Int_t pdgdau=dau->GetPdgCode();
1455  if(TMath::Abs(pdgdau)==321){
1456  nKaons++;
1457  sumPxDau+=dau->Px();
1458  sumPyDau+=dau->Py();
1459  sumPzDau+=dau->Pz();
1460  arrayDauLab[nFoundKpi++]=indDau;
1461  if(nFoundKpi>3) return -1;
1462  }else if(TMath::Abs(pdgdau)==211){
1463  nPions++;
1464  sumPxDau+=dau->Px();
1465  sumPyDau+=dau->Py();
1466  sumPzDau+=dau->Pz();
1467  arrayDauLab[nFoundKpi++]=indDau;
1468  if(nFoundKpi>3) return -1;
1469  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1470  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1471  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1472  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1473  Int_t nResDau=dau->GetNDaughters();
1474  if(nResDau!=2) return -1;
1475  Int_t indFirstResDau=dau->GetDaughter(0);
1476  for(Int_t resDau=0; resDau<2; resDau++){
1477  Int_t indResDau=indFirstResDau+resDau;
1478  if(indResDau<0) return -1;
1479  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1480  if(!resdau) return -1;
1481  Int_t pdgresdau=resdau->GetPdgCode();
1482  if(TMath::Abs(pdgresdau)==321){
1483  sumPxDau+=resdau->Px();
1484  sumPyDau+=resdau->Py();
1485  sumPzDau+=resdau->Pz();
1486  nKaons++;
1487  arrayDauLab[nFoundKpi++]=indResDau;
1488  if(nFoundKpi>3) return -1;
1489  }
1490  if(TMath::Abs(pdgresdau)==211){
1491  sumPxDau+=resdau->Px();
1492  sumPyDau+=resdau->Py();
1493  sumPzDau+=resdau->Pz();
1494  nPions++;
1495  arrayDauLab[nFoundKpi++]=indResDau;
1496  if(nFoundKpi>3) return -1;
1497  }
1498  }
1499  }else{
1500  return -1;
1501  }
1502  }
1503  if(nPions!=1) return -1;
1504  if(nKaons!=2) return -1;
1505  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1506  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1507  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1508  if(nDau==3) return 3;
1509  else if(nDau==2){
1510  if(isk0st) return 2;
1511  if(isPhi) return 1;
1512  if(isf0) return 4;
1513  }
1514  }
1515 
1516  return -1;
1517 }
1518 //____________________________________________________________________________
1519 Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1521 
1522  if(label<0) return -1;
1523  TParticle* mcPart = mcEvent->Particle(label);
1524  Int_t pdgD=mcPart->GetPdgCode();
1525  if(TMath::Abs(pdgD)!=431) return -1;
1526 
1527  Int_t nDau=mcPart->GetNDaughters();
1528  Int_t labelFirstDau = mcPart->GetDaughter(0);
1529  Int_t nPions=0;
1530  Int_t nKaons=0;
1531  Double_t sumPxDau=0.;
1532  Double_t sumPyDau=0.;
1533  Double_t sumPzDau=0.;
1534  Int_t nFoundKpi=0;
1535 
1536  Int_t codeV0=-1;
1537  if(nDau==2){
1538  for(Int_t iDau=0; iDau<nDau; iDau++){
1539  Int_t indDau = labelFirstDau+iDau;
1540  if(indDau<0) return -1;
1541  TParticle* dau=mcEvent->Particle(indDau);
1542  if(!dau) return -1;
1543  Int_t pdgdau=dau->GetPdgCode();
1544  if(TMath::Abs(pdgdau)==211){
1545  nPions++;
1546  sumPxDau+=dau->Px();
1547  sumPyDau+=dau->Py();
1548  sumPzDau+=dau->Pz();
1549  arrayDauLab[nFoundKpi++]=indDau;
1550  if(nFoundKpi>3) return -1;
1551  }else if(TMath::Abs(pdgdau)==321){
1552  nKaons++;
1553  sumPxDau+=dau->Px();
1554  sumPyDau+=dau->Py();
1555  sumPzDau+=dau->Pz();
1556  arrayDauLab[nFoundKpi++]=indDau;
1557  if(nFoundKpi>3) return -1;
1558  }else if(TMath::Abs(pdgdau)==311){
1559  codeV0=TMath::Abs(pdgdau);
1560  TParticle* v0=dau;
1561  if(codeV0==311){
1562  Int_t nK0Dau=dau->GetNDaughters();
1563  if(nK0Dau!=1) return -1;
1564  Int_t indK0s=dau->GetDaughter(0);
1565  if(indK0s<0) return -1;
1566  v0=mcEvent->Particle(indK0s);
1567  if(!v0) return -1;
1568  Int_t pdgK0sl=v0->GetPdgCode();
1569  codeV0=TMath::Abs(pdgK0sl);
1570  }
1571  Int_t nV0Dau=v0->GetNDaughters();
1572  if(nV0Dau!=2) return -1;
1573  Int_t indFirstV0Dau=v0->GetDaughter(0);
1574  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1575  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1576  if(indV0Dau<0) return -1;
1577  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1578  if(!v0dau) return -1;
1579  Int_t pdgv0dau=v0dau->GetPdgCode();
1580  if(TMath::Abs(pdgv0dau)==211){
1581  sumPxDau+=v0dau->Px();
1582  sumPyDau+=v0dau->Py();
1583  sumPzDau+=v0dau->Pz();
1584  nPions++;
1585  arrayDauLab[nFoundKpi++]=indV0Dau;
1586  if(nFoundKpi>3) return -1;
1587  }
1588  }
1589  }else{
1590  return -1;
1591  }
1592  }
1593  if(nPions!=2) return -1;
1594  if(nKaons!=1) return -1;
1595  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1596  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1597  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1598  if(codeV0==310) return 1;
1599  }
1600  return -1;
1601 
1602 }
1603 //____________________________________________________________________________
1604 Int_t AliVertexingHFUtils::CheckDstarDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1606 
1607  if(label<0) return -1;
1608  TParticle* mcPart = mcEvent->Particle(label);
1609  Int_t pdgD=mcPart->GetPdgCode();
1610  if(TMath::Abs(pdgD)!=413) return -1;
1611 
1612  Int_t nDau=mcPart->GetNDaughters();
1613  if(nDau!=2) return -1;
1614 
1615  Int_t labelFirstDau = mcPart->GetDaughter(0);
1616  Int_t nKaons=0;
1617  Int_t nPions=0;
1618  Double_t sumPxDau=0.;
1619  Double_t sumPyDau=0.;
1620  Double_t sumPzDau=0.;
1621  Int_t nFoundKpi=0;
1622 
1623  for(Int_t iDau=0; iDau<nDau; iDau++){
1624  Int_t indDau = labelFirstDau+iDau;
1625  if(indDau<0) return -1;
1626  TParticle* dau=mcEvent->Particle(indDau);
1627  if(!dau) return -1;
1628  Int_t pdgdau=dau->GetPdgCode();
1629  if(TMath::Abs(pdgdau)==421){
1630  Int_t nResDau=dau->GetNDaughters();
1631  if(nResDau!=2) return -1;
1632  Int_t indFirstResDau=dau->GetDaughter(0);
1633  for(Int_t resDau=0; resDau<2; resDau++){
1634  Int_t indResDau=indFirstResDau+resDau;
1635  if(indResDau<0) return -1;
1636  TParticle* resdau=mcEvent->Particle(indResDau);
1637  if(!resdau) return -1;
1638  Int_t pdgresdau=resdau->GetPdgCode();
1639  if(TMath::Abs(pdgresdau)==321){
1640  if(pdgD*pdgresdau>0) return -1;
1641  sumPxDau+=resdau->Px();
1642  sumPyDau+=resdau->Py();
1643  sumPzDau+=resdau->Pz();
1644  nKaons++;
1645  arrayDauLab[nFoundKpi++]=indResDau;
1646  if(nFoundKpi>3) return -1;
1647  }
1648  if(TMath::Abs(pdgresdau)==211){
1649  if(pdgD*pdgresdau<0) return -1;
1650  sumPxDau+=resdau->Px();
1651  sumPyDau+=resdau->Py();
1652  sumPzDau+=resdau->Pz();
1653  nPions++;
1654  arrayDauLab[nFoundKpi++]=indResDau;
1655  if(nFoundKpi>3) return -1;
1656  }
1657  }
1658  }else if(TMath::Abs(pdgdau)==211){
1659  if(pdgD*pdgdau<0) return -1;
1660  nPions++;
1661  sumPxDau+=dau->Px();
1662  sumPyDau+=dau->Py();
1663  sumPzDau+=dau->Pz();
1664  arrayDauLab[nFoundKpi++]=indDau;
1665  if(nFoundKpi>3) return -1;
1666  }
1667  }
1668 
1669  if(nPions!=2) return -1;
1670  if(nKaons!=1) return -1;
1671  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1672  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1673  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1674  return 1;
1675 
1676 }
1677 
1678 //____________________________________________________________________________
1679 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1681 
1682  Int_t pdgD=mcPart->GetPdgCode();
1683  if(TMath::Abs(pdgD)!=413) return -1;
1684 
1685  Int_t nDau=mcPart->GetNDaughters();
1686  if(nDau!=2) return -1;
1687 
1688  Int_t labelFirstDau = mcPart->GetDaughter(0);
1689  Int_t nKaons=0;
1690  Int_t nPions=0;
1691  Double_t sumPxDau=0.;
1692  Double_t sumPyDau=0.;
1693  Double_t sumPzDau=0.;
1694  Int_t nFoundKpi=0;
1695 
1696  for(Int_t iDau=0; iDau<nDau; iDau++){
1697  Int_t indDau = labelFirstDau+iDau;
1698  if(indDau<0) return -1;
1699  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1700  if(!dau) return -1;
1701  Int_t pdgdau=dau->GetPdgCode();
1702  if(TMath::Abs(pdgdau)==421){
1703  Int_t nResDau=dau->GetNDaughters();
1704  if(nResDau!=2) return -1;
1705  Int_t indFirstResDau=dau->GetDaughter(0);
1706  for(Int_t resDau=0; resDau<2; resDau++){
1707  Int_t indResDau=indFirstResDau+resDau;
1708  if(indResDau<0) return -1;
1709  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1710  if(!resdau) return -1;
1711  Int_t pdgresdau=resdau->GetPdgCode();
1712  if(TMath::Abs(pdgresdau)==321){
1713  if(pdgD*pdgresdau>0) return -1;
1714  sumPxDau+=resdau->Px();
1715  sumPyDau+=resdau->Py();
1716  sumPzDau+=resdau->Pz();
1717  nKaons++;
1718  arrayDauLab[nFoundKpi++]=indResDau;
1719  if(nFoundKpi>3) return -1;
1720  }
1721  if(TMath::Abs(pdgresdau)==211){
1722  if(pdgD*pdgresdau<0) return -1;
1723  sumPxDau+=resdau->Px();
1724  sumPyDau+=resdau->Py();
1725  sumPzDau+=resdau->Pz();
1726  nPions++;
1727  arrayDauLab[nFoundKpi++]=indResDau;
1728  if(nFoundKpi>3) return -1;
1729  }
1730  }
1731  }else if(TMath::Abs(pdgdau)==211){
1732  if(pdgD*pdgdau<0) return -1;
1733  nPions++;
1734  sumPxDau+=dau->Px();
1735  sumPyDau+=dau->Py();
1736  sumPzDau+=dau->Pz();
1737  arrayDauLab[nFoundKpi++]=indDau;
1738  if(nFoundKpi>3) return -1;
1739  }
1740  }
1741 
1742  if(nPions!=2) return -1;
1743  if(nKaons!=1) return -1;
1744  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1745  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1746  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1747  return 1;
1748 
1749 }
1750 
1751 //____________________________________________________________________________
1752 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1754 
1755  if(label<0) return -1;
1756  TParticle* mcPart = mcEvent->Particle(label);
1757  Int_t pdgD=mcPart->GetPdgCode();
1758  if(TMath::Abs(pdgD)!=4122) return -1;
1759 
1760  Int_t nDau=mcPart->GetNDaughters();
1761  Int_t labelFirstDau = mcPart->GetDaughter(0);
1762  Int_t nKaons=0;
1763  Int_t nPions=0;
1764  Int_t nProtons=0;
1765  Double_t sumPxDau=0.;
1766  Double_t sumPyDau=0.;
1767  Double_t sumPzDau=0.;
1768  Int_t nFoundpKpi=0;
1769 
1770  Int_t codeRes=-1;
1771  if(nDau==3 || nDau==2){
1772  for(Int_t iDau=0; iDau<nDau; iDau++){
1773  Int_t indDau = labelFirstDau+iDau;
1774  if(indDau<0) return -1;
1775  TParticle* dau=mcEvent->Particle(indDau);
1776  if(!dau) return -1;
1777  Int_t pdgdau=dau->GetPdgCode();
1778  if(TMath::Abs(pdgdau)==321){
1779  nKaons++;
1780  sumPxDau+=dau->Px();
1781  sumPyDau+=dau->Py();
1782  sumPzDau+=dau->Pz();
1783  arrayDauLab[nFoundpKpi++]=indDau;
1784  if(nFoundpKpi>3) return -1;
1785  }else if(TMath::Abs(pdgdau)==211){
1786  nPions++;
1787  sumPxDau+=dau->Px();
1788  sumPyDau+=dau->Py();
1789  sumPzDau+=dau->Pz();
1790  arrayDauLab[nFoundpKpi++]=indDau;
1791  if(nFoundpKpi>3) return -1;
1792  }else if(TMath::Abs(pdgdau)==2212){
1793  nProtons++;
1794  sumPxDau+=dau->Px();
1795  sumPyDau+=dau->Py();
1796  sumPzDau+=dau->Pz();
1797  arrayDauLab[nFoundpKpi++]=indDau;
1798  if(nFoundpKpi>3) return -1;
1799  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1800  TMath::Abs(pdgdau)==2224){
1801  codeRes=TMath::Abs(pdgdau);
1802  Int_t nResDau=dau->GetNDaughters();
1803  if(nResDau!=2) return -1;
1804  Int_t indFirstResDau=dau->GetDaughter(0);
1805  for(Int_t resDau=0; resDau<2; resDau++){
1806  Int_t indResDau=indFirstResDau+resDau;
1807  if(indResDau<0) return -1;
1808  TParticle* resdau=mcEvent->Particle(indResDau);
1809  if(!resdau) return -1;
1810  Int_t pdgresdau=resdau->GetPdgCode();
1811  if(TMath::Abs(pdgresdau)==321){
1812  sumPxDau+=resdau->Px();
1813  sumPyDau+=resdau->Py();
1814  sumPzDau+=resdau->Pz();
1815  nKaons++;
1816  arrayDauLab[nFoundpKpi++]=indResDau;
1817  if(nFoundpKpi>3) return -1;
1818  }else if(TMath::Abs(pdgresdau)==211){
1819  sumPxDau+=resdau->Px();
1820  sumPyDau+=resdau->Py();
1821  sumPzDau+=resdau->Pz();
1822  nPions++;
1823  arrayDauLab[nFoundpKpi++]=indResDau;
1824  if(nFoundpKpi>3) return -1;
1825  }else if(TMath::Abs(pdgresdau)==2212){
1826  sumPxDau+=resdau->Px();
1827  sumPyDau+=resdau->Py();
1828  sumPzDau+=resdau->Pz();
1829  nProtons++;
1830  arrayDauLab[nFoundpKpi++]=indResDau;
1831  if(nFoundpKpi>3) return -1;
1832  }
1833  }
1834  }else{
1835  return -1;
1836  }
1837  }
1838  if(nPions!=1) return -1;
1839  if(nKaons!=1) return -1;
1840  if(nProtons!=1) return -1;
1841  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1842  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1843  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1844  if(nDau==3) return 1;
1845  else if(nDau==2){
1846  if(codeRes==313) return 2;
1847  else if(codeRes==2224) return 3;
1848  else if(codeRes==3124) return 4;
1849  }
1850  }
1851  return -1;
1852 
1853 }
1854 
1855 //____________________________________________________________________________
1856 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1858 
1859  Int_t pdgD=mcPart->GetPdgCode();
1860  if(TMath::Abs(pdgD)!=4122) return -1;
1861 
1862  Int_t nDau=mcPart->GetNDaughters();
1863  Int_t labelFirstDau = mcPart->GetDaughter(0);
1864  Int_t nKaons=0;
1865  Int_t nPions=0;
1866  Int_t nProtons=0;
1867  Double_t sumPxDau=0.;
1868  Double_t sumPyDau=0.;
1869  Double_t sumPzDau=0.;
1870  Int_t nFoundpKpi=0;
1871 
1872  Int_t codeRes=-1;
1873  if(nDau==3 || nDau==2){
1874  for(Int_t iDau=0; iDau<nDau; iDau++){
1875  Int_t indDau = labelFirstDau+iDau;
1876  if(indDau<0) return -1;
1877  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1878  if(!dau) return -1;
1879  Int_t pdgdau=dau->GetPdgCode();
1880  if(TMath::Abs(pdgdau)==321){
1881  nKaons++;
1882  sumPxDau+=dau->Px();
1883  sumPyDau+=dau->Py();
1884  sumPzDau+=dau->Pz();
1885  arrayDauLab[nFoundpKpi++]=indDau;
1886  if(nFoundpKpi>3) return -1;
1887  }else if(TMath::Abs(pdgdau)==211){
1888  nPions++;
1889  sumPxDau+=dau->Px();
1890  sumPyDau+=dau->Py();
1891  sumPzDau+=dau->Pz();
1892  arrayDauLab[nFoundpKpi++]=indDau;
1893  if(nFoundpKpi>3) return -1;
1894  }else if(TMath::Abs(pdgdau)==2212){
1895  nProtons++;
1896  sumPxDau+=dau->Px();
1897  sumPyDau+=dau->Py();
1898  sumPzDau+=dau->Pz();
1899  arrayDauLab[nFoundpKpi++]=indDau;
1900  if(nFoundpKpi>3) return -1;
1901  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1902  TMath::Abs(pdgdau)==2224){
1903  codeRes=TMath::Abs(pdgdau);
1904  Int_t nResDau=dau->GetNDaughters();
1905  if(nResDau!=2) return -1;
1906  Int_t indFirstResDau=dau->GetDaughter(0);
1907  for(Int_t resDau=0; resDau<2; resDau++){
1908  Int_t indResDau=indFirstResDau+resDau;
1909  if(indResDau<0) return -1;
1910  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1911  if(!resdau) return -1;
1912  Int_t pdgresdau=resdau->GetPdgCode();
1913  if(TMath::Abs(pdgresdau)==321){
1914  sumPxDau+=resdau->Px();
1915  sumPyDau+=resdau->Py();
1916  sumPzDau+=resdau->Pz();
1917  nKaons++;
1918  arrayDauLab[nFoundpKpi++]=indResDau;
1919  if(nFoundpKpi>3) return -1;
1920  }else if(TMath::Abs(pdgresdau)==211){
1921  sumPxDau+=resdau->Px();
1922  sumPyDau+=resdau->Py();
1923  sumPzDau+=resdau->Pz();
1924  nPions++;
1925  arrayDauLab[nFoundpKpi++]=indResDau;
1926  if(nFoundpKpi>3) return -1;
1927  }else if(TMath::Abs(pdgresdau)==2212){
1928  sumPxDau+=resdau->Px();
1929  sumPyDau+=resdau->Py();
1930  sumPzDau+=resdau->Pz();
1931  nProtons++;
1932  arrayDauLab[nFoundpKpi++]=indResDau;
1933  if(nFoundpKpi>3) return -1;
1934  }
1935  }
1936  }else{
1937  return -1;
1938  }
1939  }
1940  if(nPions!=1) return -1;
1941  if(nKaons!=1) return -1;
1942  if(nProtons!=1) return -1;
1943  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1944  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1945  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1946  if(nDau==3) return 1;
1947  else if(nDau==2){
1948  if(codeRes==313) return 2;
1949  else if(codeRes==2224) return 3;
1950  else if(codeRes==3124) return 4;
1951  }
1952  }
1953  return -1;
1954 
1955 }
1956 //____________________________________________________________________________
1957 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1959 
1960  if(label<0) return -1;
1961  TParticle* mcPart = mcEvent->Particle(label);
1962  Int_t pdgD=mcPart->GetPdgCode();
1963  if(TMath::Abs(pdgD)!=4122) return -1;
1964 
1965  Int_t nDau=mcPart->GetNDaughters();
1966  Int_t labelFirstDau = mcPart->GetDaughter(0);
1967  Int_t nPions=0;
1968  Int_t nProtons=0;
1969  Double_t sumPxDau=0.;
1970  Double_t sumPyDau=0.;
1971  Double_t sumPzDau=0.;
1972  Int_t nFoundppi=0;
1973 
1974  Int_t codeV0=-1;
1975  if(nDau==2){
1976  for(Int_t iDau=0; iDau<nDau; iDau++){
1977  Int_t indDau = labelFirstDau+iDau;
1978  if(indDau<0) return -1;
1979  TParticle* dau=mcEvent->Particle(indDau);
1980  if(!dau) return -1;
1981  Int_t pdgdau=dau->GetPdgCode();
1982  if(TMath::Abs(pdgdau)==211){
1983  nPions++;
1984  sumPxDau+=dau->Px();
1985  sumPyDau+=dau->Py();
1986  sumPzDau+=dau->Pz();
1987  arrayDauLab[nFoundppi++]=indDau;
1988  if(nFoundppi>3) return -1;
1989  }else if(TMath::Abs(pdgdau)==2212){
1990  nProtons++;
1991  sumPxDau+=dau->Px();
1992  sumPyDau+=dau->Py();
1993  sumPzDau+=dau->Pz();
1994  arrayDauLab[nFoundppi++]=indDau;
1995  if(nFoundppi>3) return -1;
1996  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
1997  codeV0=TMath::Abs(pdgdau);
1998  TParticle* v0=dau;
1999  if(codeV0==311){
2000  Int_t nK0Dau=dau->GetNDaughters();
2001  if(nK0Dau!=1) return -1;
2002  Int_t indK0s=dau->GetDaughter(0);
2003  if(indK0s<0) return -1;
2004  v0=mcEvent->Particle(indK0s);
2005  if(!v0) return -1;
2006  Int_t pdgK0sl=v0->GetPdgCode();
2007  codeV0=TMath::Abs(pdgK0sl);
2008  }
2009  Int_t nV0Dau=v0->GetNDaughters();
2010  if(nV0Dau!=2) return -1;
2011  Int_t indFirstV0Dau=v0->GetDaughter(0);
2012  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2013  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2014  if(indV0Dau<0) return -1;
2015  TParticle* v0dau=mcEvent->Particle(indV0Dau);
2016  if(!v0dau) return -1;
2017  Int_t pdgv0dau=v0dau->GetPdgCode();
2018  if(TMath::Abs(pdgv0dau)==211){
2019  sumPxDau+=v0dau->Px();
2020  sumPyDau+=v0dau->Py();
2021  sumPzDau+=v0dau->Pz();
2022  nPions++;
2023  arrayDauLab[nFoundppi++]=indV0Dau;
2024  if(nFoundppi>3) return -1;
2025  }else if(TMath::Abs(pdgv0dau)==2212){
2026  sumPxDau+=v0dau->Px();
2027  sumPyDau+=v0dau->Py();
2028  sumPzDau+=v0dau->Pz();
2029  nProtons++;
2030  arrayDauLab[nFoundppi++]=indV0Dau;
2031  if(nFoundppi>3) return -1;
2032  }
2033  }
2034  }else{
2035  return -1;
2036  }
2037  }
2038  if(nPions!=2) return -1;
2039  if(nProtons!=1) return -1;
2040  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2041  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2042  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2043  if(codeV0==310) return 1;
2044  else if(codeV0==3122) return 2;
2045  }
2046  return -1;
2047 
2048 }
2049 //____________________________________________________________________________
2050 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2052 
2053  Int_t pdgD=mcPart->GetPdgCode();
2054  if(TMath::Abs(pdgD)!=4122) return -1;
2055 
2056  Int_t nDau=mcPart->GetNDaughters();
2057  Int_t labelFirstDau = mcPart->GetDaughter(0);
2058  Int_t nPions=0;
2059  Int_t nProtons=0;
2060  Double_t sumPxDau=0.;
2061  Double_t sumPyDau=0.;
2062  Double_t sumPzDau=0.;
2063  Int_t nFoundppi=0;
2064 
2065  Int_t codeV0=-1;
2066  if(nDau==2){
2067  for(Int_t iDau=0; iDau<nDau; iDau++){
2068  Int_t indDau = labelFirstDau+iDau;
2069  if(indDau<0) return -1;
2070  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2071  if(!dau) return -1;
2072  Int_t pdgdau=dau->GetPdgCode();
2073  if(TMath::Abs(pdgdau)==211){
2074  nPions++;
2075  sumPxDau+=dau->Px();
2076  sumPyDau+=dau->Py();
2077  sumPzDau+=dau->Pz();
2078  arrayDauLab[nFoundppi++]=indDau;
2079  if(nFoundppi>3) return -1;
2080  }else if(TMath::Abs(pdgdau)==2212){
2081  nProtons++;
2082  sumPxDau+=dau->Px();
2083  sumPyDau+=dau->Py();
2084  sumPzDau+=dau->Pz();
2085  arrayDauLab[nFoundppi++]=indDau;
2086  if(nFoundppi>3) return -1;
2087  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2088  codeV0=TMath::Abs(pdgdau);
2089  AliAODMCParticle* v0=dau;
2090  if(codeV0==311){
2091  Int_t nK0Dau=dau->GetNDaughters();
2092  if(nK0Dau!=1) return -1;
2093  Int_t indK0s=dau->GetDaughter(0);
2094  if(indK0s<0) return -1;
2095  v0=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indK0s));
2096  if(!v0) return -1;
2097  Int_t pdgK0sl=v0->GetPdgCode();
2098  codeV0=TMath::Abs(pdgK0sl);
2099  }
2100  Int_t nV0Dau=v0->GetNDaughters();
2101  if(nV0Dau!=2) return -1;
2102  Int_t indFirstV0Dau=v0->GetDaughter(0);
2103  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2104  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2105  if(indV0Dau<0) return -1;
2106  AliAODMCParticle* v0dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indV0Dau));
2107  if(!v0dau) return -1;
2108  Int_t pdgv0dau=v0dau->GetPdgCode();
2109  if(TMath::Abs(pdgv0dau)==211){
2110  sumPxDau+=v0dau->Px();
2111  sumPyDau+=v0dau->Py();
2112  sumPzDau+=v0dau->Pz();
2113  nPions++;
2114  arrayDauLab[nFoundppi++]=indV0Dau;
2115  if(nFoundppi>3) return -1;
2116  }else if(TMath::Abs(pdgv0dau)==2212){
2117  sumPxDau+=v0dau->Px();
2118  sumPyDau+=v0dau->Py();
2119  sumPzDau+=v0dau->Pz();
2120  nProtons++;
2121  arrayDauLab[nFoundppi++]=indV0Dau;
2122  if(nFoundppi>3) return -1;
2123  }
2124  }
2125  }else{
2126  return -1;
2127  }
2128  }
2129  if(nPions!=2) return -1;
2130  if(nProtons!=1) return -1;
2131  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2132  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2133  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2134  if(codeV0==310) return 1;
2135  else if(codeV0==3122) return 2;
2136  }
2137  return -1;
2138 
2139 }
2140 
2141 //__________________________________xic______________________________________
2142 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2144 
2145  if(label<0) return -1;
2146  TParticle* mcPart = mcEvent->Particle(label);
2147  Int_t pdgD=mcPart->GetPdgCode();
2148  if(TMath::Abs(pdgD)!=4232) return -1;
2149 
2150  Int_t nDau=mcPart->GetNDaughters();
2151  if(nDau!=3) return -1;
2152 
2153  Int_t labelFirstDau = mcPart->GetDaughter(0);
2154  Int_t nXi=0;
2155  Int_t nPions=0;
2156  Double_t sumPxDau=0.;
2157  Double_t sumPyDau=0.;
2158  Double_t sumPzDau=0.;
2159  Int_t nFoundXi=0;
2160 
2161  for(Int_t iDau=0; iDau<nDau; iDau++){
2162  Int_t indDau = labelFirstDau+iDau;
2163  if(indDau<0) return -1;
2164  TParticle* dau=mcEvent->Particle(indDau);
2165  if(!dau) return -1;
2166  Int_t pdgdau=dau->GetPdgCode();
2167  if(TMath::Abs(pdgdau)==3312){
2168  if(pdgD*pdgdau<0) return -1;
2169  sumPxDau+=dau->Px();
2170  sumPyDau+=dau->Py();
2171  sumPzDau+=dau->Pz();
2172  nXi++;
2173  arrayDauLab[nFoundXi++]=indDau;
2174 
2175  }
2176  if(TMath::Abs(pdgdau)==211){
2177  if(pdgD*pdgdau<0) return -1;
2178  nPions++;
2179  sumPxDau+=dau->Px();
2180  sumPyDau+=dau->Py();
2181  sumPzDau+=dau->Pz();
2182  arrayDauLab[nFoundXi++]=indDau;
2183  if(nFoundXi>3) return -1;
2184  }
2185  }
2186 
2187  if(nPions!=2) return -1;
2188  if(nXi!=1) return -1;
2189  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2190  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2191  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2192  return 1;
2193 
2194 }
2195 
2196 //____________________________________________________________________________
2197 Int_t AliVertexingHFUtils::CheckBplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2199 
2200  if(label<0) return -1;
2201  TParticle* mcPart = mcEvent->Particle(label);
2202  Int_t pdgD=mcPart->GetPdgCode();
2203  if(TMath::Abs(pdgD)!=521) return -1;
2204 
2205  Int_t nDau=mcPart->GetNDaughters();
2206  if(nDau!=2) return -1;
2207 
2208  Int_t labelFirstDau = mcPart->GetDaughter(0);
2209  Int_t nKaons=0;
2210  Int_t nPions=0;
2211  Double_t sumPxDau=0.;
2212  Double_t sumPyDau=0.;
2213  Double_t sumPzDau=0.;
2214  Int_t nFoundKpi=0;
2215 
2216  for(Int_t iDau=0; iDau<nDau; iDau++){
2217  Int_t indDau = labelFirstDau+iDau;
2218  if(indDau<0) return -1;
2219  TParticle* dau=mcEvent->Particle(indDau);
2220  if(!dau) return -1;
2221  Int_t pdgdau=dau->GetPdgCode();
2222  if(TMath::Abs(pdgdau)==421){
2223  Int_t nResDau=dau->GetNDaughters();
2224  if(nResDau!=2) return -1;
2225  Int_t indFirstResDau=dau->GetDaughter(0);
2226  for(Int_t resDau=0; resDau<2; resDau++){
2227  Int_t indResDau=indFirstResDau+resDau;
2228  if(indResDau<0) return -1;
2229  TParticle* resdau=mcEvent->Particle(indResDau);
2230  if(!resdau) return -1;
2231  Int_t pdgresdau=resdau->GetPdgCode();
2232  if(TMath::Abs(pdgresdau)==321){
2233  if(pdgD*pdgresdau<0) return -1;
2234  sumPxDau+=resdau->Px();
2235  sumPyDau+=resdau->Py();
2236  sumPzDau+=resdau->Pz();
2237  nKaons++;
2238  arrayDauLab[nFoundKpi++]=indResDau;
2239  if(nFoundKpi>3) return -1;
2240  }
2241  if(TMath::Abs(pdgresdau)==211){
2242  if(pdgD*pdgresdau>0) return -1;
2243  sumPxDau+=resdau->Px();
2244  sumPyDau+=resdau->Py();
2245  sumPzDau+=resdau->Pz();
2246  nPions++;
2247  arrayDauLab[nFoundKpi++]=indResDau;
2248  if(nFoundKpi>3) return -1;
2249  }
2250  }
2251  }else if(TMath::Abs(pdgdau)==211){
2252  if(pdgD*pdgdau<0) return -1;
2253  nPions++;
2254  sumPxDau+=dau->Px();
2255  sumPyDau+=dau->Py();
2256  sumPzDau+=dau->Pz();
2257  arrayDauLab[nFoundKpi++]=indDau;
2258  if(nFoundKpi>3) return -1;
2259  }
2260  }
2261 
2262  if(nPions!=2) return -1;
2263  if(nKaons!=1) return -1;
2264  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2265  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2266  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2267  return 1;
2268 
2269 }
2270 //____________________________________________________________________________
2271 Int_t AliVertexingHFUtils::CheckBplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2273 
2274  Int_t pdgD=mcPart->GetPdgCode();
2275  if(TMath::Abs(pdgD)!=521) return -1;
2276 
2277  Int_t nDau=mcPart->GetNDaughters();
2278  if(nDau!=2) return -1;
2279 
2280  Int_t labelFirstDau = mcPart->GetDaughter(0);
2281  Int_t nKaons=0;
2282  Int_t nPions=0;
2283  Double_t sumPxDau=0.;
2284  Double_t sumPyDau=0.;
2285  Double_t sumPzDau=0.;
2286  Int_t nFoundKpi=0;
2287 
2288  for(Int_t iDau=0; iDau<nDau; iDau++){
2289  Int_t indDau = labelFirstDau+iDau;
2290  if(indDau<0) return -1;
2291  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2292  if(!dau) return -1;
2293  Int_t pdgdau=dau->GetPdgCode();
2294  if(TMath::Abs(pdgdau)==421){
2295  Int_t nResDau=dau->GetNDaughters();
2296  if(nResDau!=2) return -1;
2297  Int_t indFirstResDau=dau->GetDaughter(0);
2298  for(Int_t resDau=0; resDau<2; resDau++){
2299  Int_t indResDau=indFirstResDau+resDau;
2300  if(indResDau<0) return -1;
2301  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
2302  if(!resdau) return -1;
2303  Int_t pdgresdau=resdau->GetPdgCode();
2304  if(TMath::Abs(pdgresdau)==321){
2305  if(pdgD*pdgresdau<0) return -1;
2306  sumPxDau+=resdau->Px();
2307  sumPyDau+=resdau->Py();
2308  sumPzDau+=resdau->Pz();
2309  nKaons++;
2310  arrayDauLab[nFoundKpi++]=indResDau;
2311  if(nFoundKpi>3) return -1;
2312  }
2313  if(TMath::Abs(pdgresdau)==211){
2314  if(pdgD*pdgresdau>0) return -1;
2315  sumPxDau+=resdau->Px();
2316  sumPyDau+=resdau->Py();
2317  sumPzDau+=resdau->Pz();
2318  nPions++;
2319  arrayDauLab[nFoundKpi++]=indResDau;
2320  if(nFoundKpi>3) return -1;
2321  }
2322  }
2323  }else if(TMath::Abs(pdgdau)==211){
2324  if(pdgD*pdgdau<0) return -1;
2325  nPions++;
2326  sumPxDau+=dau->Px();
2327  sumPyDau+=dau->Py();
2328  sumPzDau+=dau->Pz();
2329  arrayDauLab[nFoundKpi++]=indDau;
2330  if(nFoundKpi>3) return -1;
2331  }
2332  }
2333 
2334  if(nPions!=2) return -1;
2335  if(nKaons!=1) return -1;
2336  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2337  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2338  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2339  return 1;
2340 
2341 }
2342 //________________________________________________________________________
2344  Double_t etaMin, Double_t etaMax,
2346  Int_t filtbit1, Int_t filtbit2,
2347  Int_t minMult){
2349 
2350  Int_t nTracks=aod->GetNumberOfTracks();
2351  Int_t nSelTracks=0;
2352 
2353  Double_t sumpt=0.;
2354  Double_t s00=0.;
2355  Double_t s01=0.;
2356  Double_t s11=0.;
2357  if(ptMin<0.) ptMin=0.;
2358 
2359  for(Int_t it=0; it<nTracks; it++) {
2360  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2361  if(!tr) continue;
2362  Float_t eta = tr->Eta();
2363  Float_t pt = tr->Pt();
2364  Float_t phi = tr->Phi();
2365  if(eta<etaMin || eta>etaMax) continue;
2366  if(pt<ptMin || pt>ptMax) continue;
2367  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2368  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2369  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2370  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2371  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2372  if( !(fb1 || fb2) ) continue;
2373  Double_t px=pt*TMath::Cos(phi);
2374  Double_t py=pt*TMath::Sin(phi);
2375  s00 += (px * px)/pt;
2376  s01 += (py * px)/pt;
2377  s11 += (py * py)/pt;
2378  nSelTracks++;
2379  sumpt+=pt;
2380  }
2381 
2382  if(nSelTracks<minMult) return -0.5;
2383 
2384  if(sumpt>0.){
2385  s00/=sumpt;
2386  s01/=sumpt;
2387  s11/=sumpt;
2388  }else return -0.5;
2389 
2390  Double_t sphericity = -10;
2391  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2392  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2393  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
2394  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
2395  return sphericity;
2396 
2397 }
2398 
2399 //________________________________________________________________________
2401  Double_t &spherocity, Double_t &phiRef,
2402  Double_t etaMin, Double_t etaMax,
2404  Int_t filtbit1, Int_t filtbit2,
2405  Int_t minMult, Double_t phiStepSizeDeg,
2406  Int_t nTrksToSkip, Int_t* idToSkip
2407  ){
2409 
2410  Int_t nTracks=aod->GetNumberOfTracks();
2411  Int_t nSelTracks=0;
2412 
2413  Double_t* ptArr=new Double_t[nTracks];
2414  Double_t* phiArr=new Double_t[nTracks];
2415  Double_t sumpt=0.;
2416 
2417  for(Int_t it=0; it<nTracks; it++) {
2418  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2419  if(!tr) continue;
2420  Float_t eta = tr->Eta();
2421  Float_t pt = tr->Pt();
2422  Float_t phi = tr->Phi();
2423  if(eta<etaMin || eta>etaMax) continue;
2424  if(pt<ptMin || pt>ptMax) continue;
2425  if(nTrksToSkip>0 && idToSkip){
2426  Int_t trid = (Int_t)tr->GetID();
2427  Bool_t keep=kTRUE;
2428  for(Int_t jt=0; jt<nTrksToSkip; jt++){
2429  if(trid==idToSkip[jt]) keep=kFALSE;
2430  }
2431  if(!keep) continue;
2432  }
2433  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2434  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2435  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2436  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2437  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2438  if( !(fb1 || fb2) ) continue;
2439  ptArr[nSelTracks]=pt;
2440  phiArr[nSelTracks]=phi;
2441  nSelTracks++;
2442  sumpt+=pt;
2443  }
2444 
2445  if(nSelTracks<minMult){spherocity = -0.5; return;}
2446 
2447  //Getting thrust
2448  spherocity=2.;
2449  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2450  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2451  Double_t nx=TMath::Cos(phistep);
2452  Double_t ny=TMath::Sin(phistep);
2453  Double_t numer=0.;
2454  for(Int_t j=0; j<nSelTracks; ++j){
2455  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2456  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2457  numer+=TMath::Abs(ny*pxA - nx*pyA);
2458  }
2459  Double_t pFull=numer*numer/(sumpt*sumpt);
2460  if(pFull<spherocity){
2461  spherocity=pFull; // minimization;
2462  phiRef=phistep;
2463  }
2464  }
2465 
2466  delete [] ptArr;
2467  delete [] phiArr;
2468 
2469  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2470  return;
2471 
2472 }
2473 //________________________________________________________________________
2475  Double_t &spherocity, Double_t &phiRef,
2476  Double_t etaMin, Double_t etaMax,
2478  Int_t minMult, Double_t phiStepSizeDeg){
2479 
2481 
2482  Int_t nParticles=arrayMC->GetEntriesFast();
2483  Int_t nSelParticles=0;
2484 
2485  Double_t* ptArr=new Double_t[nParticles];
2486  Double_t* phiArr=new Double_t[nParticles];
2487  Double_t sumpt=0.;
2488 
2489  for(Int_t ip=0; ip<nParticles; ip++) {
2490  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(ip);
2491  if(!part) continue;
2492  Float_t eta = part->Eta();
2493  Float_t pt = part->Pt();
2494  Float_t phi = part->Phi();
2495  Int_t charge = part->Charge();
2496  Bool_t isPhysPrim = part->IsPhysicalPrimary();
2497  if(!isPhysPrim) continue;
2498  if(charge==0) continue;
2499  if(eta<etaMin || eta>etaMax) continue;
2500  if(pt<ptMin || pt>ptMax) continue;
2501 
2502  ptArr[nSelParticles]=pt;
2503  phiArr[nSelParticles]=phi;
2504  nSelParticles++;
2505  sumpt+=pt;
2506  }
2507 
2508  if(nSelParticles<minMult){spherocity = -0.5; return;}
2509 
2510  //Getting thrust
2511  spherocity=2.;
2512  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2513  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2514  Double_t nx=TMath::Cos(phistep);
2515  Double_t ny=TMath::Sin(phistep);
2516  Double_t numer=0.;
2517  for(Int_t j=0; j<nSelParticles; ++j){
2518  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2519  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2520  numer+=TMath::Abs(ny*pxA - nx*pyA);
2521  }
2522  Double_t pFull=numer*numer/(sumpt*sumpt);
2523  if(pFull<spherocity){
2524  spherocity=pFull; // minimization;
2525  phiRef=phistep;
2526  }
2527  }
2528 
2529  delete [] ptArr;
2530  delete [] phiArr;
2531 
2532  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2533  return;
2534 
2535 }
2536 
2537 //________________________________________________________________________
2544 
2545  Int_t nBinOrig=hOrig->GetNbinsX();
2546  Int_t firstBinOrig=1;
2547  Int_t lastBinOrig=nBinOrig;
2548  Int_t nBinOrigUsed=nBinOrig;
2549  Int_t nBinFinal=nBinOrig/reb;
2550  if(firstUse>=1){
2551  firstBinOrig=firstUse;
2552  nBinFinal=(nBinOrig-firstUse+1)/reb;
2553  nBinOrigUsed=nBinFinal*reb;
2554  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2555  }else{
2556  Int_t exc=nBinOrigUsed%reb;
2557  if(exc!=0){
2558  nBinOrigUsed-=exc;
2559  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2560  }
2561  }
2562 
2563  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
2564  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
2565  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
2566  TH1D* hRebin=new TH1D(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
2567  Int_t lastSummed=firstBinOrig-1;
2568  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
2569  Float_t sum=0.;
2570  Float_t sume2=0.;
2571  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
2572  sum+=hOrig->GetBinContent(lastSummed+1);
2573  sume2+=(hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1));
2574  lastSummed++;
2575  }
2576  hRebin->SetBinContent(iBin,sum);
2577  hRebin->SetBinError(iBin,TMath::Sqrt(sume2));
2578  }
2579  return hRebin;
2580 }
2581 //________________________________________________________________________
2584 
2585  Int_t binmin=TMath::Max(1,hData->FindBin(hMC->GetXaxis()->GetXmin()));
2586  Bool_t found=kFALSE;
2587  Int_t binminD=-1;
2588  Int_t binminMC=-1;
2589  for(Int_t j=binmin; j<hData->GetNbinsX(); j++){
2590  if(found) break;
2591  for(Int_t k=1; k<hMC->GetNbinsX(); k++){
2592  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k)-hData->GetBinLowEdge(j));
2593  if(delta<0.0001){
2594  found=kTRUE;
2595  binminMC=k;
2596  binminD=j;
2597  }
2598  if(found) break;
2599  }
2600  }
2601  Int_t binmax=TMath::Min(hData->GetNbinsX(),hData->FindBin(hMC->GetXaxis()->GetXmax()*0.99999));
2602  found=kFALSE;
2603  Int_t binmaxD=-1;
2604  Int_t binmaxMC=-1;
2605  for(Int_t j=binmax; j>1; j--){
2606  if(found) break;
2607  for(Int_t k=hMC->GetNbinsX(); k>400; k--){
2608  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k+1)-hData->GetBinLowEdge(j+1));
2609  if(delta<0.0001){
2610  found=kTRUE;
2611  binmaxMC=k;
2612  binmaxD=j;
2613  }
2614  if(found) break;
2615  }
2616  }
2617 
2618  Double_t min=hData->GetBinLowEdge(binminD);
2619  Double_t max=hData->GetBinLowEdge(binmaxD)+hData->GetBinWidth(binmaxD);
2620  Double_t minMC=hMC->GetBinLowEdge(binminMC);
2621  Double_t maxMC=hMC->GetBinLowEdge(binmaxMC)+hMC->GetBinWidth(binmaxMC);
2622  Double_t width=hData->GetBinWidth(binminD);
2623  Double_t widthMC=hMC->GetBinWidth(binminMC);
2624 
2625  if(TMath::Abs(minMC-min)>0.0001*min || TMath::Abs(maxMC-max)>0.0001*max){
2626  printf("Cannot adapt range and rebin histo:\n");
2627  printf("Range for data histo: %f-%f GeV/c2 bins %d-%d width=%f\n",min,max,binminD,binmaxD,width);
2628  printf("Range for reflection histo: %f-%f GeV/c2 bins %d-%d width=%f\n",minMC,maxMC,binminMC,binmaxMC,widthMC);
2629  return 0x0;
2630  }
2631 
2632  Double_t rebin=width/widthMC;
2633  if(TMath::Abs(rebin-TMath::Nint(rebin))>0.001){
2634  printf("Cannot adapt histo: rebin %f issue, width MC = %f, width hData=%f (check=%f)\n",rebin,widthMC,width,TMath::Abs(rebin-TMath::Nint(rebin)));
2635  return 0x0;
2636  }
2637 
2638  Int_t nBinsNew=binmaxD-binminD+1;
2639  TH1 *hOut;
2640  TString stype=hMC->ClassName();
2641  if(stype.Contains("TH1F")){
2642  hOut=new TH1F(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2643  }else if(stype.Contains("TH1D")){
2644  hOut=new TH1D(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2645  }else{
2646  printf("Wrong type %s\n",stype.Data());
2647  return 0x0;
2648  }
2649 
2650  for(Int_t j=1; j<=hMC->GetNbinsX(); j++){
2651  Double_t m=hMC->GetBinCenter(j);
2652  Int_t binFin=hOut->FindBin(m);
2653  if(binFin>=1 && binFin<=nBinsNew){
2654  hOut->AddBinContent(binFin,hMC->GetBinContent(j));
2655  }
2656  }
2657  return hOut;
2658 }
Int_t charge
static void GetSpherocity(AliAODEvent *aod, Double_t &spherocity, Double_t &phiRef, 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, Int_t nTrksToSkip=0, Int_t *idToSkip=0x0)
Functions for event shape variables.
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
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.
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
double Double_t
Definition: External.C:58
Definition: External.C:236
static Int_t CheckDplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckLcpKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
AliAODv0 * Getv0() const
Double_t ImpParXY() const
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
TCanvas * c
Definition: TestFitELoss.C:172
Double_t Pol(Double_t x) const
static Double_t GetFullEvResolLowLim(const TH1F *hSubEvCorr, Int_t k=1)
static Int_t CheckBplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Double_t ptMin
static TH1 * AdaptTemplateRangeAndBinning(const TH1 *hRef, TH1 *hData, Double_t minFit, Double_t maxFit)
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
TRandom * gRandom
void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
Double_t fMinEtaForTracklets
sub-event resolution = sqrt(<cos[n(phiA-phiB)] >)
Double_t fSubRes
ratio of measured harmonic to event plane harmonic
static Int_t CheckDplusK0spiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckXicXipipiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckLcV0bachelorDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
Double_t fMaxEtaForTracklets
min eta for counting tracklets
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
static Double_t GetFullEvResolHighLim(const TH1F *hSubEvCorr, Int_t k=1)
static void GetGeneratedSpherocity(TClonesArray *arrayMC, Double_t &spherocity, Double_t &phiRef, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t minMult=3, Double_t phiStepSizeDeg=0.1)
static Double_t ResolK1(Double_t x)
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
int Int_t
Definition: External.C:63
Double_t GetFullEvResol() const
float Float_t
Definition: External.C:68
const Double_t ptmax
static Int_t CheckDplusKKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Bool_t CheckT0TriggerFired(AliAODEvent *aodEv)
Functions for processing trigger information.
Definition: External.C:212
AliAODTrack * GetBachelor() const
static Double_t GetSubEvResolLowLim(const TH1F *hSubEvCorr)
const Double_t ptmin
Bool_t HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
Bool_t IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header, TClonesArray *arrayMC)
static Int_t CheckDstarDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
short Short_t
Definition: External.C:23
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
bool Bool_t
Definition: External.C:53
static Int_t CheckDsK0sKDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
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)
Definition: External.C:196
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 //.