AliPhysics  9c7580a (9c7580a)
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 
2051 //__________________________________xic______________________________________
2052 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2054 
2055  if(label<0) return -1;
2056  TParticle* mcPart = mcEvent->Particle(label);
2057  Int_t pdgD=mcPart->GetPdgCode();
2058  if(TMath::Abs(pdgD)!=4232) return -1;
2059 
2060  Int_t nDau=mcPart->GetNDaughters();
2061  if(nDau!=3) return -1;
2062 
2063  Int_t labelFirstDau = mcPart->GetDaughter(0);
2064  Int_t nXi=0;
2065  Int_t nPions=0;
2066  Double_t sumPxDau=0.;
2067  Double_t sumPyDau=0.;
2068  Double_t sumPzDau=0.;
2069  Int_t nFoundXi=0;
2070 
2071  for(Int_t iDau=0; iDau<nDau; iDau++){
2072  Int_t indDau = labelFirstDau+iDau;
2073  if(indDau<0) return -1;
2074  TParticle* dau=mcEvent->Particle(indDau);
2075  if(!dau) return -1;
2076  Int_t pdgdau=dau->GetPdgCode();
2077  if(TMath::Abs(pdgdau)==3312){
2078  if(pdgD*pdgdau<0) return -1;
2079  sumPxDau+=dau->Px();
2080  sumPyDau+=dau->Py();
2081  sumPzDau+=dau->Pz();
2082  nXi++;
2083  arrayDauLab[nFoundXi++]=indDau;
2084 
2085  }
2086  if(TMath::Abs(pdgdau)==211){
2087  if(pdgD*pdgdau<0) return -1;
2088  nPions++;
2089  sumPxDau+=dau->Px();
2090  sumPyDau+=dau->Py();
2091  sumPzDau+=dau->Pz();
2092  arrayDauLab[nFoundXi++]=indDau;
2093  if(nFoundXi>3) return -1;
2094  }
2095  }
2096 
2097  if(nPions!=2) return -1;
2098  if(nXi!=1) return -1;
2099  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2100  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2101  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2102  return 1;
2103 
2104 }
2105 //________________________________________________________________________
2107  Double_t etaMin, Double_t etaMax,
2109  Int_t filtbit1, Int_t filtbit2,
2110  Int_t minMult){
2112 
2113  Int_t nTracks=aod->GetNumberOfTracks();
2114  Int_t nSelTracks=0;
2115 
2116  Double_t sumpt=0.;
2117  Double_t s00=0.;
2118  Double_t s01=0.;
2119  Double_t s11=0.;
2120  if(ptMin<0.) ptMin=0.;
2121 
2122  for(Int_t it=0; it<nTracks; it++) {
2123  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2124  if(!tr) continue;
2125  Float_t eta = tr->Eta();
2126  Float_t pt = tr->Pt();
2127  Float_t phi = tr->Phi();
2128  if(eta<etaMin || eta>etaMax) continue;
2129  if(pt<ptMin || pt>ptMax) continue;
2130  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2131  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2132  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2133  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2134  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2135  if( !(fb1 || fb2) ) continue;
2136  Double_t px=pt*TMath::Cos(phi);
2137  Double_t py=pt*TMath::Sin(phi);
2138  s00 += (px * px)/pt;
2139  s01 += (py * px)/pt;
2140  s11 += (py * py)/pt;
2141  nSelTracks++;
2142  sumpt+=pt;
2143  }
2144 
2145  if(nSelTracks<minMult) return -0.5;
2146 
2147  if(sumpt>0.){
2148  s00/=sumpt;
2149  s01/=sumpt;
2150  s11/=sumpt;
2151  }else return -0.5;
2152 
2153  Double_t sphericity = -10;
2154  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2155  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2156  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
2157  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
2158  return sphericity;
2159 
2160 }
2161 
2162 //________________________________________________________________________
2164  Double_t &spherocity, Double_t &phiRef,
2165  Double_t etaMin, Double_t etaMax,
2167  Int_t filtbit1, Int_t filtbit2,
2168  Int_t minMult, Double_t phiStepSizeDeg,
2169  Int_t nTrksToSkip, Int_t* idToSkip
2170  ){
2172 
2173  Int_t nTracks=aod->GetNumberOfTracks();
2174  Int_t nSelTracks=0;
2175 
2176  Double_t* ptArr=new Double_t[nTracks];
2177  Double_t* phiArr=new Double_t[nTracks];
2178  Double_t sumpt=0.;
2179 
2180  for(Int_t it=0; it<nTracks; it++) {
2181  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2182  if(!tr) continue;
2183  Float_t eta = tr->Eta();
2184  Float_t pt = tr->Pt();
2185  Float_t phi = tr->Phi();
2186  if(eta<etaMin || eta>etaMax) continue;
2187  if(pt<ptMin || pt>ptMax) continue;
2188  if(nTrksToSkip>0 && idToSkip){
2189  Int_t trid = (Int_t)tr->GetID();
2190  Bool_t keep=kTRUE;
2191  for(Int_t jt=0; jt<nTrksToSkip; jt++){
2192  if(trid==idToSkip[jt]) keep=kFALSE;
2193  }
2194  if(!keep) continue;
2195  }
2196  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2197  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2198  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2199  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2200  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2201  if( !(fb1 || fb2) ) continue;
2202  ptArr[nSelTracks]=pt;
2203  phiArr[nSelTracks]=phi;
2204  nSelTracks++;
2205  sumpt+=pt;
2206  }
2207 
2208  if(nSelTracks<minMult){spherocity = -0.5; return;}
2209 
2210  //Getting thrust
2211  spherocity=2.;
2212  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2213  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2214  Double_t nx=TMath::Cos(phistep);
2215  Double_t ny=TMath::Sin(phistep);
2216  Double_t numer=0.;
2217  for(Int_t j=0; j<nSelTracks; ++j){
2218  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2219  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2220  numer+=TMath::Abs(ny*pxA - nx*pyA);
2221  }
2222  Double_t pFull=numer*numer/(sumpt*sumpt);
2223  if(pFull<spherocity){
2224  spherocity=pFull; // minimization;
2225  phiRef=phistep;
2226  }
2227  }
2228 
2229  delete [] ptArr;
2230  delete [] phiArr;
2231 
2232  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2233  return;
2234 
2235 }
2236 //________________________________________________________________________
2238  Double_t &spherocity, Double_t &phiRef,
2239  Double_t etaMin, Double_t etaMax,
2241  Int_t minMult, Double_t phiStepSizeDeg){
2242 
2244 
2245  Int_t nParticles=arrayMC->GetEntriesFast();
2246  Int_t nSelParticles=0;
2247 
2248  Double_t* ptArr=new Double_t[nParticles];
2249  Double_t* phiArr=new Double_t[nParticles];
2250  Double_t sumpt=0.;
2251 
2252  for(Int_t ip=0; ip<nParticles; ip++) {
2253  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(ip);
2254  if(!part) continue;
2255  Float_t eta = part->Eta();
2256  Float_t pt = part->Pt();
2257  Float_t phi = part->Phi();
2258  Int_t charge = part->Charge();
2259  Bool_t isPhysPrim = part->IsPhysicalPrimary();
2260  if(!isPhysPrim) continue;
2261  if(charge==0) continue;
2262  if(eta<etaMin || eta>etaMax) continue;
2263  if(pt<ptMin || pt>ptMax) continue;
2264 
2265  ptArr[nSelParticles]=pt;
2266  phiArr[nSelParticles]=phi;
2267  nSelParticles++;
2268  sumpt+=pt;
2269  }
2270 
2271  if(nSelParticles<minMult){spherocity = -0.5; return;}
2272 
2273  //Getting thrust
2274  spherocity=2.;
2275  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2276  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2277  Double_t nx=TMath::Cos(phistep);
2278  Double_t ny=TMath::Sin(phistep);
2279  Double_t numer=0.;
2280  for(Int_t j=0; j<nSelParticles; ++j){
2281  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2282  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2283  numer+=TMath::Abs(ny*pxA - nx*pyA);
2284  }
2285  Double_t pFull=numer*numer/(sumpt*sumpt);
2286  if(pFull<spherocity){
2287  spherocity=pFull; // minimization;
2288  phiRef=phistep;
2289  }
2290  }
2291 
2292  delete [] ptArr;
2293  delete [] phiArr;
2294 
2295  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2296  return;
2297 
2298 }
2299 
2300 //________________________________________________________________________
2307 
2308  Int_t nBinOrig=hOrig->GetNbinsX();
2309  Int_t firstBinOrig=1;
2310  Int_t lastBinOrig=nBinOrig;
2311  Int_t nBinOrigUsed=nBinOrig;
2312  Int_t nBinFinal=nBinOrig/reb;
2313  if(firstUse>=1){
2314  firstBinOrig=firstUse;
2315  nBinFinal=(nBinOrig-firstUse+1)/reb;
2316  nBinOrigUsed=nBinFinal*reb;
2317  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2318  }else{
2319  Int_t exc=nBinOrigUsed%reb;
2320  if(exc!=0){
2321  nBinOrigUsed-=exc;
2322  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2323  }
2324  }
2325 
2326  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
2327  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
2328  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
2329  TH1D* hRebin=new TH1D(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
2330  Int_t lastSummed=firstBinOrig-1;
2331  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
2332  Float_t sum=0.;
2333  Float_t sume2=0.;
2334  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
2335  sum+=hOrig->GetBinContent(lastSummed+1);
2336  sume2+=(hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1));
2337  lastSummed++;
2338  }
2339  hRebin->SetBinContent(iBin,sum);
2340  hRebin->SetBinError(iBin,TMath::Sqrt(sume2));
2341  }
2342  return hRebin;
2343 }
2344 //________________________________________________________________________
2347 
2348  Int_t binmin=TMath::Max(1,hData->FindBin(hMC->GetXaxis()->GetXmin()));
2349  Bool_t found=kFALSE;
2350  Int_t binminD=-1;
2351  Int_t binminMC=-1;
2352  for(Int_t j=binmin; j<hData->GetNbinsX(); j++){
2353  if(found) break;
2354  for(Int_t k=1; k<hMC->GetNbinsX(); k++){
2355  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k)-hData->GetBinLowEdge(j));
2356  if(delta<0.0001){
2357  found=kTRUE;
2358  binminMC=k;
2359  binminD=j;
2360  }
2361  if(found) break;
2362  }
2363  }
2364  Int_t binmax=TMath::Min(hData->GetNbinsX(),hData->FindBin(hMC->GetXaxis()->GetXmax()*0.99999));
2365  found=kFALSE;
2366  Int_t binmaxD=-1;
2367  Int_t binmaxMC=-1;
2368  for(Int_t j=binmax; j>1; j--){
2369  if(found) break;
2370  for(Int_t k=hMC->GetNbinsX(); k>400; k--){
2371  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k+1)-hData->GetBinLowEdge(j+1));
2372  if(delta<0.0001){
2373  found=kTRUE;
2374  binmaxMC=k;
2375  binmaxD=j;
2376  }
2377  if(found) break;
2378  }
2379  }
2380 
2381  Double_t min=hData->GetBinLowEdge(binminD);
2382  Double_t max=hData->GetBinLowEdge(binmaxD)+hData->GetBinWidth(binmaxD);
2383  Double_t minMC=hMC->GetBinLowEdge(binminMC);
2384  Double_t maxMC=hMC->GetBinLowEdge(binmaxMC)+hMC->GetBinWidth(binmaxMC);
2385  Double_t width=hData->GetBinWidth(binminD);
2386  Double_t widthMC=hMC->GetBinWidth(binminMC);
2387 
2388  if(TMath::Abs(minMC-min)>0.0001*min || TMath::Abs(maxMC-max)>0.0001*max){
2389  printf("Cannot adapt range and rebin histo:\n");
2390  printf("Range for data histo: %f-%f GeV/c2 bins %d-%d width=%f\n",min,max,binminD,binmaxD,width);
2391  printf("Range for reflection histo: %f-%f GeV/c2 bins %d-%d width=%f\n",minMC,maxMC,binminMC,binmaxMC,widthMC);
2392  return 0x0;
2393  }
2394 
2395  Double_t rebin=width/widthMC;
2396  if(TMath::Abs(rebin-TMath::Nint(rebin))>0.001){
2397  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)));
2398  return 0x0;
2399  }
2400 
2401  Int_t nBinsNew=binmaxD-binminD+1;
2402  TH1 *hOut;
2403  TString stype=hMC->ClassName();
2404  if(stype.Contains("TH1F")){
2405  hOut=new TH1F(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2406  }else if(stype.Contains("TH1D")){
2407  hOut=new TH1D(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2408  }else{
2409  printf("Wrong type %s\n",stype.Data());
2410  return 0x0;
2411  }
2412 
2413  for(Int_t j=1; j<=hMC->GetNbinsX(); j++){
2414  Double_t m=hMC->GetBinCenter(j);
2415  Int_t binFin=hOut->FindBin(m);
2416  if(binFin>=1 && binFin<=nBinsNew){
2417  hOut->AddBinContent(binFin,hMC->GetBinContent(j));
2418  }
2419  }
2420  return hOut;
2421 }
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)
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
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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 //.