AliPhysics  9b6b435 (9b6b435)
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, AliMCParticle *mcPart, Bool_t searchUpToQuark){
607 
608  Int_t pdgGranma = 0;
609  Int_t mother = 0;
610  mother = mcPart->GetMother();
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  AliMCParticle* mcGranma = (AliMCParticle*)mcEvent->GetTrack(mother);
618  if (mcGranma){
619  TParticle* partGranma = mcGranma->Particle();
620  if(partGranma) pdgGranma = partGranma->GetPdgCode();
621  abspdgGranma = TMath::Abs(pdgGranma);
622  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
623  isFromB=kTRUE;
624  }
625  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
626  mother = mcGranma->GetMother();
627  }else{
628  printf("CheckOrigin: Failed casting the mother particle!");
629  break;
630  }
631  }
632  if(searchUpToQuark && !isQuarkFound) return 0;
633  if(isFromB) return 5;
634  else return 4;
635 
636 }
637 //____________________________________________________________________________
638 Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
640 
641  Int_t pdgGranma = 0;
642  Int_t mother = 0;
643  mother = mcPart->GetMother();
644  Int_t istep = 0;
645  Int_t abspdgGranma =0;
646  Bool_t isFromB=kFALSE;
647  Bool_t isQuarkFound=kFALSE;
648  while (mother >0 ){
649  istep++;
650  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
651  if (mcGranma){
652  pdgGranma = mcGranma->GetPdgCode();
653  abspdgGranma = TMath::Abs(pdgGranma);
654  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
655  isFromB=kTRUE;
656  }
657  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
658  mother = mcGranma->GetMother();
659  }else{
660  printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
661  break;
662  }
663  }
664  if(searchUpToQuark && !isQuarkFound) return 0;
665  if(isFromB) return 5;
666  else return 4;
667 
668 }
669 //____________________________________________________________________________
670 Double_t AliVertexingHFUtils::GetBeautyMotherPt(TClonesArray* arrayMC, AliAODMCParticle *mcPart){
672 
673  Int_t pdgGranma = 0;
674  Int_t mother = 0;
675  mother = mcPart->GetMother();
676  Int_t istep = 0;
677  Int_t abspdgGranma =0;
678  while (mother >0 ){
679  istep++;
680  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
681  if (mcGranma){
682  pdgGranma = mcGranma->GetPdgCode();
683  abspdgGranma = TMath::Abs(pdgGranma);
684  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
685  return mcGranma->Pt();
686  }
687  if(abspdgGranma==4) return -999.;
688  if(abspdgGranma==5) return -1.;
689  mother = mcGranma->GetMother();
690  }else{
691  printf("AliVertexingHFUtils::GetBeautyMotherPt: Failed casting the mother particle!");
692  break;
693  }
694  }
695  return -999.;
696 }
697 //____________________________________________________________________________
698 Int_t AliVertexingHFUtils::CheckD0Decay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
700 
701  if(label<0) return -1;
702  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
703  TParticle* part = mcEvent->Particle(label);
704  if(!part || !mcPart) return -1;
705  Int_t pdgD=part->GetPdgCode();
706  if(TMath::Abs(pdgD)!=421) return -1;
707 
708  Int_t nDau=mcPart->GetNDaughters();
709 
710  if(nDau==2){
711  Int_t daughter0 = mcPart->GetFirstDaughter();
712  Int_t daughter1 = mcPart->GetLastDaughter();
713  TParticle* mcPartDaughter0 = mcEvent->Particle(daughter0);
714  TParticle* mcPartDaughter1 = mcEvent->Particle(daughter1);
715  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
716  arrayDauLab[0]=daughter0;
717  arrayDauLab[1]=daughter1;
718  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
719  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
720  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
721  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
722  return -1;
723  }
724  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
725  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
726  if((pdgCode0*pdgCode1)>0) return -1;
727  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
728  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
729  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
730  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
731  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
732  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
733  return 1;
734  }
735 
736  if(nDau==3 || nDau==4){
737  Int_t nKaons=0;
738  Int_t nPions=0;
739  Double_t sumPxDau=0.;
740  Double_t sumPyDau=0.;
741  Double_t sumPzDau=0.;
742  Int_t nFoundKpi=0;
743  Int_t labelFirstDau = mcPart->GetFirstDaughter();
744  for(Int_t iDau=0; iDau<nDau; iDau++){
745  Int_t indDau = labelFirstDau+iDau;
746  if(indDau<0) return -1;
747  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
748  TParticle* dau=mcEvent->Particle(indDau);
749  if(!mcDau || !dau) return -1;
750  Int_t pdgdau=dau->GetPdgCode();
751  if(TMath::Abs(pdgdau)==321){
752  if(pdgD>0 && pdgdau>0) return -1;
753  if(pdgD<0 && pdgdau<0) return -1;
754  nKaons++;
755  sumPxDau+=dau->Px();
756  sumPyDau+=dau->Py();
757  sumPzDau+=dau->Pz();
758  arrayDauLab[nFoundKpi++]=indDau;
759  if(nFoundKpi>4) return -1;
760  }else if(TMath::Abs(pdgdau)==211){
761  nPions++;
762  sumPxDau+=dau->Px();
763  sumPyDau+=dau->Py();
764  sumPzDau+=dau->Pz();
765  arrayDauLab[nFoundKpi++]=indDau;
766  if(nFoundKpi>4) return -1;
767  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
768  Int_t nResDau=mcDau->GetNDaughters();
769  if(nResDau!=2) return -1;
770  Int_t indFirstResDau=mcDau->GetFirstDaughter();
771  for(Int_t resDau=0; resDau<2; resDau++){
772  Int_t indResDau=indFirstResDau+resDau;
773  if(indResDau<0) return -1;
774  TParticle* resdau=mcEvent->Particle(indResDau);
775  if(!resdau) return -1;
776  Int_t pdgresdau=resdau->GetPdgCode();
777  if(TMath::Abs(pdgresdau)==321){
778  if(pdgD>0 && pdgresdau>0) return -1;
779  if(pdgD<0 && pdgresdau<0) return -1;
780  nKaons++;
781  sumPxDau+=resdau->Px();
782  sumPyDau+=resdau->Py();
783  sumPzDau+=resdau->Pz();
784  arrayDauLab[nFoundKpi++]=indResDau;
785  if(nFoundKpi>4) return -1;
786  }
787  if(TMath::Abs(pdgresdau)==211){
788  nPions++;
789  sumPxDau+=resdau->Px();
790  sumPyDau+=resdau->Py();
791  sumPzDau+=resdau->Pz();
792  arrayDauLab[nFoundKpi++]=indResDau;
793  if(nFoundKpi>4) return -1;
794  }
795  }
796  }else{
797  return -1;
798  }
799  }
800  if(nPions!=3) return -1;
801  if(nKaons!=1) return -1;
802  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -1;
803  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -1;
804  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -1;
805  return 2;
806  }
807  return -1;
808 }
809 
810 //____________________________________________________________________________
811 Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
813 
814  Int_t pdgD=mcPart->GetPdgCode();
815  if(TMath::Abs(pdgD)!=421) return -1;
816 
817  Int_t nDau=mcPart->GetNDaughters();
818 
819  if(nDau==2){
820  Int_t daughter0 = mcPart->GetDaughter(0);
821  Int_t daughter1 = mcPart->GetDaughter(1);
822  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
823  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
824  if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
825  arrayDauLab[0]=daughter0;
826  arrayDauLab[1]=daughter1;
827  Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
828  Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
829  if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
830  !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
831  return -1;
832  }
833  if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
834  if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
835  if((pdgCode0*pdgCode1)>0) return -1;
836  Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
837  Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
838  Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
839  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
840  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
841  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
842  return 1;
843  }
844 
845  if(nDau==3 || nDau==4){
846  Int_t nKaons=0;
847  Int_t nPions=0;
848  Double_t sumPxDau=0.;
849  Double_t sumPyDau=0.;
850  Double_t sumPzDau=0.;
851  Int_t nFoundKpi=0;
852  Int_t labelFirstDau = mcPart->GetDaughter(0);
853  for(Int_t iDau=0; iDau<nDau; iDau++){
854  Int_t indDau = labelFirstDau+iDau;
855  if(indDau<0) return -1;
856  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
857  if(!dau) return -1;
858  Int_t pdgdau=dau->GetPdgCode();
859  if(TMath::Abs(pdgdau)==321){
860  if(pdgD>0 && pdgdau>0) return -1;
861  if(pdgD<0 && pdgdau<0) return -1;
862  nKaons++;
863  sumPxDau+=dau->Px();
864  sumPyDau+=dau->Py();
865  sumPzDau+=dau->Pz();
866  arrayDauLab[nFoundKpi++]=indDau;
867  if(nFoundKpi>4) return -1;
868  }else if(TMath::Abs(pdgdau)==211){
869  nPions++;
870  sumPxDau+=dau->Px();
871  sumPyDau+=dau->Py();
872  sumPzDau+=dau->Pz();
873  arrayDauLab[nFoundKpi++]=indDau;
874  if(nFoundKpi>4) return -1;
875  }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
876  Int_t nResDau=dau->GetNDaughters();
877  if(nResDau!=2) return -1;
878  Int_t indFirstResDau=dau->GetDaughter(0);
879  for(Int_t resDau=0; resDau<2; resDau++){
880  Int_t indResDau=indFirstResDau+resDau;
881  if(indResDau<0) return -1;
882  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
883  if(!resdau) return -1;
884  Int_t pdgresdau=resdau->GetPdgCode();
885  if(TMath::Abs(pdgresdau)==321){
886  if(pdgD>0 && pdgresdau>0) return -1;
887  if(pdgD<0 && pdgresdau<0) return -1;
888  nKaons++;
889  sumPxDau+=resdau->Px();
890  sumPyDau+=resdau->Py();
891  sumPzDau+=resdau->Pz();
892  arrayDauLab[nFoundKpi++]=indResDau;
893  if(nFoundKpi>4) return -1;
894  }
895  if(TMath::Abs(pdgresdau)==211){
896  nPions++;
897  sumPxDau+=resdau->Px();
898  sumPyDau+=resdau->Py();
899  sumPzDau+=resdau->Pz();
900  arrayDauLab[nFoundKpi++]=indResDau;
901  if(nFoundKpi>4) return -1;
902  }
903  }
904  }else{
905  return -1;
906  }
907  }
908  if(nPions!=3) return -1;
909  if(nKaons!=1) return -1;
910  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
911  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
912  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
913  return 2;
914  }
915  return -1;
916 }
917 //____________________________________________________________________________
918 Int_t AliVertexingHFUtils::CheckDplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
920 
921  if(label<0) return -1;
922  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
923  TParticle* part = mcEvent->Particle(label);
924  if(!part || !mcPart) return -1;
925  Int_t pdgD=part->GetPdgCode();
926  if(TMath::Abs(pdgD)!=411) return -1;
927 
928  Int_t nDau=mcPart->GetNDaughters();
929  Int_t labelFirstDau = mcPart->GetFirstDaughter();
930  Int_t nKaons=0;
931  Int_t nPions=0;
932  Double_t sumPxDau=0.;
933  Double_t sumPyDau=0.;
934  Double_t sumPzDau=0.;
935  Int_t nFoundKpi=0;
936 
937  if(nDau==3 || nDau==2){
938  for(Int_t iDau=0; iDau<nDau; iDau++){
939  Int_t indDau = labelFirstDau+iDau;
940  if(indDau<0) return -1;
941  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
942  TParticle* dau=mcEvent->Particle(indDau);
943  if(!mcDau || !dau) return -1;
944  Int_t pdgdau=dau->GetPdgCode();
945  if(TMath::Abs(pdgdau)==321){
946  if(pdgD*pdgdau>0) return -1;
947  nKaons++;
948  sumPxDau+=dau->Px();
949  sumPyDau+=dau->Py();
950  sumPzDau+=dau->Pz();
951  arrayDauLab[nFoundKpi++]=indDau;
952  if(nFoundKpi>3) return -1;
953  }else if(TMath::Abs(pdgdau)==211){
954  if(pdgD*pdgdau<0) return -1;
955  nPions++;
956  sumPxDau+=dau->Px();
957  sumPyDau+=dau->Py();
958  sumPzDau+=dau->Pz();
959  arrayDauLab[nFoundKpi++]=indDau;
960  if(nFoundKpi>3) return -1;
961  }else if(TMath::Abs(pdgdau)==313){
962  Int_t nResDau=mcDau->GetNDaughters();
963  if(nResDau!=2) return -1;
964  Int_t indFirstResDau=mcDau->GetFirstDaughter();
965  for(Int_t resDau=0; resDau<2; resDau++){
966  Int_t indResDau=indFirstResDau+resDau;
967  if(indResDau<0) return -1;
968  TParticle* resdau=mcEvent->Particle(indResDau);
969  if(!resdau) return -1;
970  Int_t pdgresdau=resdau->GetPdgCode();
971  if(TMath::Abs(pdgresdau)==321){
972  if(pdgD*pdgresdau>0) return -1;
973  sumPxDau+=resdau->Px();
974  sumPyDau+=resdau->Py();
975  sumPzDau+=resdau->Pz();
976  nKaons++;
977  arrayDauLab[nFoundKpi++]=indResDau;
978  if(nFoundKpi>3) return -1;
979  }
980  if(TMath::Abs(pdgresdau)==211){
981  if(pdgD*pdgresdau<0) return -1;
982  sumPxDau+=resdau->Px();
983  sumPyDau+=resdau->Py();
984  sumPzDau+=resdau->Pz();
985  nPions++;
986  arrayDauLab[nFoundKpi++]=indResDau;
987  if(nFoundKpi>3) return -1;
988  }
989  }
990  }else{
991  return -1;
992  }
993  }
994  if(nPions!=2) return -1;
995  if(nKaons!=1) return -1;
996  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
997  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
998  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
999  if(nDau==3) return 1;
1000  else if(nDau==2) return 2;
1001  }
1002 
1003  return -1;
1004 
1005 }
1006 
1007 //____________________________________________________________________________
1008 Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1010 
1011  Int_t pdgD=mcPart->GetPdgCode();
1012  if(TMath::Abs(pdgD)!=411) return -1;
1013 
1014  Int_t nDau=mcPart->GetNDaughters();
1015  Int_t labelFirstDau = mcPart->GetDaughter(0);
1016  Int_t nKaons=0;
1017  Int_t nPions=0;
1018  Double_t sumPxDau=0.;
1019  Double_t sumPyDau=0.;
1020  Double_t sumPzDau=0.;
1021  Int_t nFoundKpi=0;
1022 
1023  if(nDau==3 || nDau==2){
1024  for(Int_t iDau=0; iDau<nDau; iDau++){
1025  Int_t indDau = labelFirstDau+iDau;
1026  if(indDau<0) return -1;
1027  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1028  if(!dau) return -1;
1029  Int_t pdgdau=dau->GetPdgCode();
1030  if(TMath::Abs(pdgdau)==321){
1031  if(pdgD*pdgdau>0) return -1;
1032  nKaons++;
1033  sumPxDau+=dau->Px();
1034  sumPyDau+=dau->Py();
1035  sumPzDau+=dau->Pz();
1036  arrayDauLab[nFoundKpi++]=indDau;
1037  if(nFoundKpi>3) return -1;
1038  }else if(TMath::Abs(pdgdau)==211){
1039  if(pdgD*pdgdau<0) return -1;
1040  nPions++;
1041  sumPxDau+=dau->Px();
1042  sumPyDau+=dau->Py();
1043  sumPzDau+=dau->Pz();
1044  arrayDauLab[nFoundKpi++]=indDau;
1045  if(nFoundKpi>3) return -1;
1046  }else if(TMath::Abs(pdgdau)==313){
1047  Int_t nResDau=dau->GetNDaughters();
1048  if(nResDau!=2) return -1;
1049  Int_t indFirstResDau=dau->GetDaughter(0);
1050  for(Int_t resDau=0; resDau<2; resDau++){
1051  Int_t indResDau=indFirstResDau+resDau;
1052  if(indResDau<0) return -1;
1053  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1054  if(!resdau) return -1;
1055  Int_t pdgresdau=resdau->GetPdgCode();
1056  if(TMath::Abs(pdgresdau)==321){
1057  if(pdgD*pdgresdau>0) return -1;
1058  sumPxDau+=resdau->Px();
1059  sumPyDau+=resdau->Py();
1060  sumPzDau+=resdau->Pz();
1061  nKaons++;
1062  arrayDauLab[nFoundKpi++]=indResDau;
1063  if(nFoundKpi>3) return -1;
1064  }
1065  if(TMath::Abs(pdgresdau)==211){
1066  if(pdgD*pdgresdau<0) return -1;
1067  sumPxDau+=resdau->Px();
1068  sumPyDau+=resdau->Py();
1069  sumPzDau+=resdau->Pz();
1070  nPions++;
1071  arrayDauLab[nFoundKpi++]=indResDau;
1072  if(nFoundKpi>3) return -1;
1073  }
1074  }
1075  }else{
1076  return -1;
1077  }
1078  }
1079  if(nPions!=2) return -1;
1080  if(nKaons!=1) return -1;
1081  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1082  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1083  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1084  if(nDau==3) return 1;
1085  else if(nDau==2) return 2;
1086  }
1087 
1088  return -1;
1089 
1090 }
1091 //____________________________________________________________________________
1092 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1094 
1095  if(label<0) return -1;
1096  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1097  TParticle* part = mcEvent->Particle(label);
1098  if(!part || !mcPart) return -1;
1099  Int_t pdgD=part->GetPdgCode();
1100  if(TMath::Abs(pdgD)!=411) return -1;
1101 
1102  Int_t nDau=mcPart->GetNDaughters();
1103  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1104  Int_t nKaons=0;
1105  Int_t nPions=0;
1106  Double_t sumPxDau=0.;
1107  Double_t sumPyDau=0.;
1108  Double_t sumPzDau=0.;
1109  Int_t nFoundKpi=0;
1110  Bool_t isPhi=kFALSE;
1111  Bool_t isk0st=kFALSE;
1112 
1113  if(nDau==3 || nDau==2){
1114  for(Int_t iDau=0; iDau<nDau; iDau++){
1115  Int_t indDau = labelFirstDau+iDau;
1116  if(indDau<0) return -1;
1117  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1118  TParticle* dau=mcEvent->Particle(indDau);
1119  if(!mcDau || !dau) return -1;
1120  Int_t pdgdau=dau->GetPdgCode();
1121  if(TMath::Abs(pdgdau)==321){
1122  nKaons++;
1123  sumPxDau+=dau->Px();
1124  sumPyDau+=dau->Py();
1125  sumPzDau+=dau->Pz();
1126  arrayDauLab[nFoundKpi++]=indDau;
1127  if(nFoundKpi>3) return -1;
1128  }else if(TMath::Abs(pdgdau)==211){
1129  nPions++;
1130  sumPxDau+=dau->Px();
1131  sumPyDau+=dau->Py();
1132  sumPzDau+=dau->Pz();
1133  arrayDauLab[nFoundKpi++]=indDau;
1134  if(nFoundKpi>3) return -1;
1135  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1136  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1137  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1138  Int_t nResDau=mcDau->GetNDaughters();
1139  if(nResDau!=2) return -1;
1140  Int_t indFirstResDau=mcDau->GetFirstDaughter();
1141  for(Int_t resDau=0; resDau<2; resDau++){
1142  Int_t indResDau=indFirstResDau+resDau;
1143  if(indResDau<0) return -1;
1144  TParticle* resdau=mcEvent->Particle(indResDau);
1145  if(!resdau) return -1;
1146  Int_t pdgresdau=resdau->GetPdgCode();
1147  if(TMath::Abs(pdgresdau)==321){
1148  sumPxDau+=resdau->Px();
1149  sumPyDau+=resdau->Py();
1150  sumPzDau+=resdau->Pz();
1151  nKaons++;
1152  arrayDauLab[nFoundKpi++]=indResDau;
1153  if(nFoundKpi>3) return -1;
1154  }
1155  if(TMath::Abs(pdgresdau)==211){
1156  sumPxDau+=resdau->Px();
1157  sumPyDau+=resdau->Py();
1158  sumPzDau+=resdau->Pz();
1159  nPions++;
1160  arrayDauLab[nFoundKpi++]=indResDau;
1161  if(nFoundKpi>3) return -1;
1162  }
1163  }
1164  }else{
1165  return -1;
1166  }
1167  }
1168  if(nPions!=1) return -1;
1169  if(nKaons!=2) return -1;
1170  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1171  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1172  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1173  if(nDau==3) return 3;
1174  else if(nDau==2){
1175  if(isk0st) return 2;
1176  if(isPhi) return 1;
1177  }
1178  }
1179 
1180  return -1;
1181 
1182 }
1183 //____________________________________________________________________________
1184 Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1186 
1187  Int_t pdgD=mcPart->GetPdgCode();
1188  if(TMath::Abs(pdgD)!=411) return -1;
1189 
1190  Int_t nDau=mcPart->GetNDaughters();
1191  Int_t labelFirstDau = mcPart->GetDaughter(0);
1192  Int_t nKaons=0;
1193  Int_t nPions=0;
1194  Double_t sumPxDau=0.;
1195  Double_t sumPyDau=0.;
1196  Double_t sumPzDau=0.;
1197  Int_t nFoundKpi=0;
1198  Bool_t isPhi=kFALSE;
1199  Bool_t isk0st=kFALSE;
1200 
1201  if(nDau==3 || nDau==2){
1202  for(Int_t iDau=0; iDau<nDau; iDau++){
1203  Int_t indDau = labelFirstDau+iDau;
1204  if(indDau<0) return -1;
1205  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1206  if(!dau) return -1;
1207  Int_t pdgdau=dau->GetPdgCode();
1208  if(TMath::Abs(pdgdau)==321){
1209  nKaons++;
1210  sumPxDau+=dau->Px();
1211  sumPyDau+=dau->Py();
1212  sumPzDau+=dau->Pz();
1213  arrayDauLab[nFoundKpi++]=indDau;
1214  if(nFoundKpi>3) return -1;
1215  }else if(TMath::Abs(pdgdau)==211){
1216  nPions++;
1217  sumPxDau+=dau->Px();
1218  sumPyDau+=dau->Py();
1219  sumPzDau+=dau->Pz();
1220  arrayDauLab[nFoundKpi++]=indDau;
1221  if(nFoundKpi>3) return -1;
1222  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
1223  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1224  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1225  Int_t nResDau=dau->GetNDaughters();
1226  if(nResDau!=2) return -1;
1227  Int_t indFirstResDau=dau->GetDaughter(0);
1228  for(Int_t resDau=0; resDau<2; resDau++){
1229  Int_t indResDau=indFirstResDau+resDau;
1230  if(indResDau<0) return -1;
1231  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1232  if(!resdau) return -1;
1233  Int_t pdgresdau=resdau->GetPdgCode();
1234  if(TMath::Abs(pdgresdau)==321){
1235  sumPxDau+=resdau->Px();
1236  sumPyDau+=resdau->Py();
1237  sumPzDau+=resdau->Pz();
1238  nKaons++;
1239  arrayDauLab[nFoundKpi++]=indResDau;
1240  if(nFoundKpi>3) return -1;
1241  }
1242  if(TMath::Abs(pdgresdau)==211){
1243  sumPxDau+=resdau->Px();
1244  sumPyDau+=resdau->Py();
1245  sumPzDau+=resdau->Pz();
1246  nPions++;
1247  arrayDauLab[nFoundKpi++]=indResDau;
1248  if(nFoundKpi>3) return -1;
1249  }
1250  }
1251  }else{
1252  return -1;
1253  }
1254  }
1255  if(nPions!=1) return -1;
1256  if(nKaons!=2) return -1;
1257  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1258  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1259  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1260  if(nDau==3) return 3;
1261  else if(nDau==2){
1262  if(isk0st) return 2;
1263  if(isPhi) return 1;
1264  }
1265  }
1266 
1267  return -1;
1268 
1269 }
1270 //____________________________________________________________________________
1271 Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1273 
1274  if(label<0) return -1;
1275  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1276  TParticle* part = mcEvent->Particle(label);
1277  if(!part || !mcPart) return -1;
1278  Int_t pdgD=part->GetPdgCode();
1279  if(TMath::Abs(pdgD)!=411) return -1;
1280 
1281  Int_t nDau=mcPart->GetNDaughters();
1282  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1283  Int_t nPions=0;
1284  Double_t sumPxDau=0.;
1285  Double_t sumPyDau=0.;
1286  Double_t sumPzDau=0.;
1287  Int_t nFoundpi=0;
1288 
1289  Int_t codeV0=-1;
1290  if(nDau==2){
1291  for(Int_t iDau=0; iDau<nDau; iDau++){
1292  Int_t indDau = labelFirstDau+iDau;
1293  if(indDau<0) return -1;
1294  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1295  TParticle* dau=mcEvent->Particle(indDau);
1296  if(!mcDau || !dau) return -1;
1297  Int_t pdgdau=dau->GetPdgCode();
1298  if(TMath::Abs(pdgdau)==211){
1299  nPions++;
1300  sumPxDau+=dau->Px();
1301  sumPyDau+=dau->Py();
1302  sumPzDau+=dau->Pz();
1303  arrayDauLab[nFoundpi++]=indDau;
1304  if(nFoundpi>3) return -1;
1305  }else if(TMath::Abs(pdgdau)==311){
1306  codeV0=TMath::Abs(pdgdau);
1307  TParticle* v0=dau;
1308  AliMCParticle* mcV0=mcDau;
1309  if(codeV0==311){
1310  Int_t nK0Dau=mcDau->GetNDaughters();
1311  if(nK0Dau!=1) return -1;
1312  Int_t indK0s=mcDau->GetFirstDaughter();
1313  if(indK0s<0) return -1;
1314  v0=mcEvent->Particle(indK0s);
1315  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1316  if(!v0 || !mcV0) return -1;
1317  Int_t pdgK0sl=v0->GetPdgCode();
1318  codeV0=TMath::Abs(pdgK0sl);
1319  }
1320  Int_t nV0Dau=mcV0->GetNDaughters();
1321  if(nV0Dau!=2) return -1;
1322  Int_t indFirstV0Dau=mcV0->GetFirstDaughter();
1323  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1324  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1325  if(indV0Dau<0) return -1;
1326  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1327  if(!v0dau) return -1;
1328  Int_t pdgv0dau=v0dau->GetPdgCode();
1329  if(TMath::Abs(pdgv0dau)==211){
1330  sumPxDau+=v0dau->Px();
1331  sumPyDau+=v0dau->Py();
1332  sumPzDau+=v0dau->Pz();
1333  nPions++;
1334  arrayDauLab[nFoundpi++]=indV0Dau;
1335  if(nFoundpi>3) return -1;
1336  }
1337  }
1338  }else{
1339  return -1;
1340  }
1341  }
1342  if(nPions!=3) return -1;
1343  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1344  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1345  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1346  if(codeV0==310) return 1;
1347  }
1348  return -1;
1349 
1350 }
1351 
1352 
1353 //____________________________________________________________________________
1354 Int_t AliVertexingHFUtils::CheckDsDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1356 
1357  if(label<0) return -1;
1358  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1359  TParticle* part = mcEvent->Particle(label);
1360  if(!part || !mcPart) return -1;
1361  Int_t pdgD=part->GetPdgCode();
1362  if(TMath::Abs(pdgD)!=431) return -1;
1363 
1364  Int_t nDau=mcPart->GetNDaughters();
1365  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1366  Int_t nKaons=0;
1367  Int_t nPions=0;
1368  Double_t sumPxDau=0.;
1369  Double_t sumPyDau=0.;
1370  Double_t sumPzDau=0.;
1371  Int_t nFoundKpi=0;
1372  Bool_t isPhi=kFALSE;
1373  Bool_t isk0st=kFALSE;
1374  Bool_t isf0=kFALSE;
1375 
1376  if(nDau==3 || nDau==2){
1377  for(Int_t iDau=0; iDau<nDau; iDau++){
1378  Int_t indDau = labelFirstDau+iDau;
1379  if(indDau<0) return -1;
1380  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1381  TParticle* dau=mcEvent->Particle(indDau);
1382  if(!mcDau || !dau) return -1;
1383  Int_t pdgdau=dau->GetPdgCode();
1384  if(TMath::Abs(pdgdau)==321){
1385  nKaons++;
1386  sumPxDau+=dau->Px();
1387  sumPyDau+=dau->Py();
1388  sumPzDau+=dau->Pz();
1389  arrayDauLab[nFoundKpi++]=indDau;
1390  if(nFoundKpi>3) return -1;
1391  }else if(TMath::Abs(pdgdau)==211){
1392  nPions++;
1393  sumPxDau+=dau->Px();
1394  sumPyDau+=dau->Py();
1395  sumPzDau+=dau->Pz();
1396  arrayDauLab[nFoundKpi++]=indDau;
1397  if(nFoundKpi>3) return -1;
1398  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1399  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1400  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1401  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1402  Int_t nResDau=mcDau->GetNDaughters();
1403  if(nResDau!=2) return -1;
1404  Int_t indFirstResDau=mcDau->GetFirstDaughter();
1405  for(Int_t resDau=0; resDau<2; resDau++){
1406  Int_t indResDau=indFirstResDau+resDau;
1407  if(indResDau<0) return -1;
1408  TParticle* resdau=mcEvent->Particle(indResDau);
1409  if(!resdau) return -1;
1410  Int_t pdgresdau=resdau->GetPdgCode();
1411  if(TMath::Abs(pdgresdau)==321){
1412  sumPxDau+=resdau->Px();
1413  sumPyDau+=resdau->Py();
1414  sumPzDau+=resdau->Pz();
1415  nKaons++;
1416  arrayDauLab[nFoundKpi++]=indResDau;
1417  if(nFoundKpi>3) return -1;
1418  }
1419  if(TMath::Abs(pdgresdau)==211){
1420  sumPxDau+=resdau->Px();
1421  sumPyDau+=resdau->Py();
1422  sumPzDau+=resdau->Pz();
1423  nPions++;
1424  arrayDauLab[nFoundKpi++]=indResDau;
1425  if(nFoundKpi>3) return -1;
1426  }
1427  }
1428  }else{
1429  return -1;
1430  }
1431  }
1432  if(nPions!=1) return -1;
1433  if(nKaons!=2) return -1;
1434  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1435  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1436  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1437  if(nDau==3) return 3;
1438  else if(nDau==2){
1439  if(isk0st) return 2;
1440  if(isPhi) return 1;
1441  if(isf0) return 4;
1442  }
1443  }
1444 
1445  return -1;
1446 }
1447 //____________________________________________________________________________
1448 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1450 
1451  Int_t pdgD=mcPart->GetPdgCode();
1452  if(TMath::Abs(pdgD)!=431) return -1;
1453 
1454  Int_t nDau=mcPart->GetNDaughters();
1455  Int_t labelFirstDau = mcPart->GetDaughter(0);
1456  Int_t nKaons=0;
1457  Int_t nPions=0;
1458  Double_t sumPxDau=0.;
1459  Double_t sumPyDau=0.;
1460  Double_t sumPzDau=0.;
1461  Int_t nFoundKpi=0;
1462  Bool_t isPhi=kFALSE;
1463  Bool_t isk0st=kFALSE;
1464  Bool_t isf0=kFALSE;
1465 
1466  if(nDau==3 || nDau==2){
1467  for(Int_t iDau=0; iDau<nDau; iDau++){
1468  Int_t indDau = labelFirstDau+iDau;
1469  if(indDau<0) return -1;
1470  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1471  if(!dau) return -1;
1472  Int_t pdgdau=dau->GetPdgCode();
1473  if(TMath::Abs(pdgdau)==321){
1474  nKaons++;
1475  sumPxDau+=dau->Px();
1476  sumPyDau+=dau->Py();
1477  sumPzDau+=dau->Pz();
1478  arrayDauLab[nFoundKpi++]=indDau;
1479  if(nFoundKpi>3) return -1;
1480  }else if(TMath::Abs(pdgdau)==211){
1481  nPions++;
1482  sumPxDau+=dau->Px();
1483  sumPyDau+=dau->Py();
1484  sumPzDau+=dau->Pz();
1485  arrayDauLab[nFoundKpi++]=indDau;
1486  if(nFoundKpi>3) return -1;
1487  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333 || TMath::Abs(pdgdau)==9010221){
1488  if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
1489  if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
1490  if(TMath::Abs(pdgdau)==9010221) isf0=kTRUE;
1491  Int_t nResDau=dau->GetNDaughters();
1492  if(nResDau!=2) return -1;
1493  Int_t indFirstResDau=dau->GetDaughter(0);
1494  for(Int_t resDau=0; resDau<2; resDau++){
1495  Int_t indResDau=indFirstResDau+resDau;
1496  if(indResDau<0) return -1;
1497  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1498  if(!resdau) return -1;
1499  Int_t pdgresdau=resdau->GetPdgCode();
1500  if(TMath::Abs(pdgresdau)==321){
1501  sumPxDau+=resdau->Px();
1502  sumPyDau+=resdau->Py();
1503  sumPzDau+=resdau->Pz();
1504  nKaons++;
1505  arrayDauLab[nFoundKpi++]=indResDau;
1506  if(nFoundKpi>3) return -1;
1507  }
1508  if(TMath::Abs(pdgresdau)==211){
1509  sumPxDau+=resdau->Px();
1510  sumPyDau+=resdau->Py();
1511  sumPzDau+=resdau->Pz();
1512  nPions++;
1513  arrayDauLab[nFoundKpi++]=indResDau;
1514  if(nFoundKpi>3) return -1;
1515  }
1516  }
1517  }else{
1518  return -1;
1519  }
1520  }
1521  if(nPions!=1) return -1;
1522  if(nKaons!=2) return -1;
1523  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1524  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1525  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1526  if(nDau==3) return 3;
1527  else if(nDau==2){
1528  if(isk0st) return 2;
1529  if(isPhi) return 1;
1530  if(isf0) return 4;
1531  }
1532  }
1533 
1534  return -1;
1535 }
1536 //____________________________________________________________________________
1537 Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1539 
1540  if(label<0) return -1;
1541  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1542  TParticle* part = mcEvent->Particle(label);
1543  if(!part || !mcPart) return -1;
1544  Int_t pdgD=part->GetPdgCode();
1545  if(TMath::Abs(pdgD)!=431) return -1;
1546 
1547  Int_t nDau=mcPart->GetNDaughters();
1548  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1549  Int_t nPions=0;
1550  Int_t nKaons=0;
1551  Double_t sumPxDau=0.;
1552  Double_t sumPyDau=0.;
1553  Double_t sumPzDau=0.;
1554  Int_t nFoundKpi=0;
1555 
1556  Int_t codeV0=-1;
1557  if(nDau==2){
1558  for(Int_t iDau=0; iDau<nDau; iDau++){
1559  Int_t indDau = labelFirstDau+iDau;
1560  if(indDau<0) return -1;
1561  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1562  TParticle* dau=mcEvent->Particle(indDau);
1563  if(!mcDau || !dau) return -1;
1564  Int_t pdgdau=dau->GetPdgCode();
1565  if(TMath::Abs(pdgdau)==211){
1566  nPions++;
1567  sumPxDau+=dau->Px();
1568  sumPyDau+=dau->Py();
1569  sumPzDau+=dau->Pz();
1570  arrayDauLab[nFoundKpi++]=indDau;
1571  if(nFoundKpi>3) return -1;
1572  }else if(TMath::Abs(pdgdau)==321){
1573  nKaons++;
1574  sumPxDau+=dau->Px();
1575  sumPyDau+=dau->Py();
1576  sumPzDau+=dau->Pz();
1577  arrayDauLab[nFoundKpi++]=indDau;
1578  if(nFoundKpi>3) return -1;
1579  }else if(TMath::Abs(pdgdau)==311){
1580  codeV0=TMath::Abs(pdgdau);
1581  TParticle* v0=dau;
1582  AliMCParticle* mcV0=mcDau;
1583  if(codeV0==311){
1584  Int_t nK0Dau=mcDau->GetNDaughters();
1585  if(nK0Dau!=1) return -1;
1586  Int_t indK0s=mcDau->GetFirstDaughter();
1587  if(indK0s<0) return -1;
1588  v0=mcEvent->Particle(indK0s);
1589  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
1590  if(!v0 || !mcV0) return -1;
1591  Int_t pdgK0sl=v0->GetPdgCode();
1592  codeV0=TMath::Abs(pdgK0sl);
1593  }
1594  Int_t nV0Dau=mcV0->GetNDaughters();
1595  if(nV0Dau!=2) return -1;
1596  Int_t indFirstV0Dau=mcV0->GetFirstDaughter();
1597  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
1598  Int_t indV0Dau=indFirstV0Dau+v0Dau;
1599  if(indV0Dau<0) return -1;
1600  TParticle* v0dau=mcEvent->Particle(indV0Dau);
1601  if(!v0dau) return -1;
1602  Int_t pdgv0dau=v0dau->GetPdgCode();
1603  if(TMath::Abs(pdgv0dau)==211){
1604  sumPxDau+=v0dau->Px();
1605  sumPyDau+=v0dau->Py();
1606  sumPzDau+=v0dau->Pz();
1607  nPions++;
1608  arrayDauLab[nFoundKpi++]=indV0Dau;
1609  if(nFoundKpi>3) return -1;
1610  }
1611  }
1612  }else{
1613  return -1;
1614  }
1615  }
1616  if(nPions!=2) return -1;
1617  if(nKaons!=1) return -1;
1618  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1619  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1620  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1621  if(codeV0==310) return 1;
1622  }
1623  return -1;
1624 
1625 }
1626 //____________________________________________________________________________
1627 Int_t AliVertexingHFUtils::CheckDstarDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1629 
1630  if(label<0) return -1;
1631  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1632  TParticle* part = mcEvent->Particle(label);
1633  if(!part || !mcPart) return -1;
1634  Int_t pdgD=part->GetPdgCode();
1635  if(TMath::Abs(pdgD)!=413) return -1;
1636 
1637  Int_t nDau=mcPart->GetNDaughters();
1638  if(nDau!=2) return -1;
1639 
1640  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1641  Int_t nKaons=0;
1642  Int_t nPions=0;
1643  Double_t sumPxDau=0.;
1644  Double_t sumPyDau=0.;
1645  Double_t sumPzDau=0.;
1646  Int_t nFoundKpi=0;
1647 
1648  for(Int_t iDau=0; iDau<nDau; iDau++){
1649  Int_t indDau = labelFirstDau+iDau;
1650  if(indDau<0) return -1;
1651  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1652  TParticle* dau=mcEvent->Particle(indDau);
1653  if(!mcDau || !dau) return -1;
1654  Int_t pdgdau=dau->GetPdgCode();
1655  if(TMath::Abs(pdgdau)==421){
1656  Int_t nResDau=mcDau->GetNDaughters();
1657  if(nResDau!=2) return -1;
1658  Int_t indFirstResDau=mcDau->GetFirstDaughter();
1659  for(Int_t resDau=0; resDau<2; resDau++){
1660  Int_t indResDau=indFirstResDau+resDau;
1661  if(indResDau<0) return -1;
1662  TParticle* resdau=mcEvent->Particle(indResDau);
1663  if(!resdau) return -1;
1664  Int_t pdgresdau=resdau->GetPdgCode();
1665  if(TMath::Abs(pdgresdau)==321){
1666  if(pdgD*pdgresdau>0) return -1;
1667  sumPxDau+=resdau->Px();
1668  sumPyDau+=resdau->Py();
1669  sumPzDau+=resdau->Pz();
1670  nKaons++;
1671  arrayDauLab[nFoundKpi++]=indResDau;
1672  if(nFoundKpi>3) return -1;
1673  }
1674  if(TMath::Abs(pdgresdau)==211){
1675  if(pdgD*pdgresdau<0) return -1;
1676  sumPxDau+=resdau->Px();
1677  sumPyDau+=resdau->Py();
1678  sumPzDau+=resdau->Pz();
1679  nPions++;
1680  arrayDauLab[nFoundKpi++]=indResDau;
1681  if(nFoundKpi>3) return -1;
1682  }
1683  }
1684  }else if(TMath::Abs(pdgdau)==211){
1685  if(pdgD*pdgdau<0) return -1;
1686  nPions++;
1687  sumPxDau+=dau->Px();
1688  sumPyDau+=dau->Py();
1689  sumPzDau+=dau->Pz();
1690  arrayDauLab[nFoundKpi++]=indDau;
1691  if(nFoundKpi>3) return -1;
1692  }
1693  }
1694 
1695  if(nPions!=2) return -1;
1696  if(nKaons!=1) return -1;
1697  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1698  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1699  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1700  return 1;
1701 
1702 }
1703 
1704 //____________________________________________________________________________
1705 Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1707 
1708  Int_t pdgD=mcPart->GetPdgCode();
1709  if(TMath::Abs(pdgD)!=413) return -1;
1710 
1711  Int_t nDau=mcPart->GetNDaughters();
1712  if(nDau!=2) return -1;
1713 
1714  Int_t labelFirstDau = mcPart->GetDaughter(0);
1715  Int_t nKaons=0;
1716  Int_t nPions=0;
1717  Double_t sumPxDau=0.;
1718  Double_t sumPyDau=0.;
1719  Double_t sumPzDau=0.;
1720  Int_t nFoundKpi=0;
1721 
1722  for(Int_t iDau=0; iDau<nDau; iDau++){
1723  Int_t indDau = labelFirstDau+iDau;
1724  if(indDau<0) return -1;
1725  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1726  if(!dau) return -1;
1727  Int_t pdgdau=dau->GetPdgCode();
1728  if(TMath::Abs(pdgdau)==421){
1729  Int_t nResDau=dau->GetNDaughters();
1730  if(nResDau!=2) return -1;
1731  Int_t indFirstResDau=dau->GetDaughter(0);
1732  for(Int_t resDau=0; resDau<2; resDau++){
1733  Int_t indResDau=indFirstResDau+resDau;
1734  if(indResDau<0) return -1;
1735  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1736  if(!resdau) return -1;
1737  Int_t pdgresdau=resdau->GetPdgCode();
1738  if(TMath::Abs(pdgresdau)==321){
1739  if(pdgD*pdgresdau>0) return -1;
1740  sumPxDau+=resdau->Px();
1741  sumPyDau+=resdau->Py();
1742  sumPzDau+=resdau->Pz();
1743  nKaons++;
1744  arrayDauLab[nFoundKpi++]=indResDau;
1745  if(nFoundKpi>3) return -1;
1746  }
1747  if(TMath::Abs(pdgresdau)==211){
1748  if(pdgD*pdgresdau<0) return -1;
1749  sumPxDau+=resdau->Px();
1750  sumPyDau+=resdau->Py();
1751  sumPzDau+=resdau->Pz();
1752  nPions++;
1753  arrayDauLab[nFoundKpi++]=indResDau;
1754  if(nFoundKpi>3) return -1;
1755  }
1756  }
1757  }else if(TMath::Abs(pdgdau)==211){
1758  if(pdgD*pdgdau<0) return -1;
1759  nPions++;
1760  sumPxDau+=dau->Px();
1761  sumPyDau+=dau->Py();
1762  sumPzDau+=dau->Pz();
1763  arrayDauLab[nFoundKpi++]=indDau;
1764  if(nFoundKpi>3) return -1;
1765  }
1766  }
1767 
1768  if(nPions!=2) return -1;
1769  if(nKaons!=1) return -1;
1770  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1771  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1772  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1773  return 1;
1774 
1775 }
1776 
1777 //____________________________________________________________________________
1778 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1780 
1781  if(label<0) return -1;
1782  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1783  TParticle* part = mcEvent->Particle(label);
1784  if(!part || !mcPart) return -1;
1785  Int_t pdgD=part->GetPdgCode();
1786  if(TMath::Abs(pdgD)!=4122) return -1;
1787 
1788  Int_t nDau=mcPart->GetNDaughters();
1789  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1790  Int_t nKaons=0;
1791  Int_t nPions=0;
1792  Int_t nProtons=0;
1793  Double_t sumPxDau=0.;
1794  Double_t sumPyDau=0.;
1795  Double_t sumPzDau=0.;
1796  Int_t nFoundpKpi=0;
1797 
1798  Int_t codeRes=-1;
1799  if(nDau==3 || nDau==2){
1800  for(Int_t iDau=0; iDau<nDau; iDau++){
1801  Int_t indDau = labelFirstDau+iDau;
1802  if(indDau<0) return -1;
1803  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
1804  TParticle* dau=mcEvent->Particle(indDau);
1805  if(!mcDau || !dau) return -1;
1806  Int_t pdgdau=dau->GetPdgCode();
1807  if(TMath::Abs(pdgdau)==321){
1808  nKaons++;
1809  sumPxDau+=dau->Px();
1810  sumPyDau+=dau->Py();
1811  sumPzDau+=dau->Pz();
1812  arrayDauLab[nFoundpKpi++]=indDau;
1813  if(nFoundpKpi>3) return -1;
1814  }else if(TMath::Abs(pdgdau)==211){
1815  nPions++;
1816  sumPxDau+=dau->Px();
1817  sumPyDau+=dau->Py();
1818  sumPzDau+=dau->Pz();
1819  arrayDauLab[nFoundpKpi++]=indDau;
1820  if(nFoundpKpi>3) return -1;
1821  }else if(TMath::Abs(pdgdau)==2212){
1822  nProtons++;
1823  sumPxDau+=dau->Px();
1824  sumPyDau+=dau->Py();
1825  sumPzDau+=dau->Pz();
1826  arrayDauLab[nFoundpKpi++]=indDau;
1827  if(nFoundpKpi>3) return -1;
1828  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1829  TMath::Abs(pdgdau)==2224){
1830  codeRes=TMath::Abs(pdgdau);
1831  Int_t nResDau=mcDau->GetNDaughters();
1832  if(nResDau!=2) return -1;
1833  Int_t indFirstResDau=mcDau->GetFirstDaughter();
1834  for(Int_t resDau=0; resDau<2; resDau++){
1835  Int_t indResDau=indFirstResDau+resDau;
1836  if(indResDau<0) return -1;
1837  TParticle* resdau=mcEvent->Particle(indResDau);
1838  if(!resdau) return -1;
1839  Int_t pdgresdau=resdau->GetPdgCode();
1840  if(TMath::Abs(pdgresdau)==321){
1841  sumPxDau+=resdau->Px();
1842  sumPyDau+=resdau->Py();
1843  sumPzDau+=resdau->Pz();
1844  nKaons++;
1845  arrayDauLab[nFoundpKpi++]=indResDau;
1846  if(nFoundpKpi>3) return -1;
1847  }else if(TMath::Abs(pdgresdau)==211){
1848  sumPxDau+=resdau->Px();
1849  sumPyDau+=resdau->Py();
1850  sumPzDau+=resdau->Pz();
1851  nPions++;
1852  arrayDauLab[nFoundpKpi++]=indResDau;
1853  if(nFoundpKpi>3) return -1;
1854  }else if(TMath::Abs(pdgresdau)==2212){
1855  sumPxDau+=resdau->Px();
1856  sumPyDau+=resdau->Py();
1857  sumPzDau+=resdau->Pz();
1858  nProtons++;
1859  arrayDauLab[nFoundpKpi++]=indResDau;
1860  if(nFoundpKpi>3) return -1;
1861  }
1862  }
1863  }else{
1864  return -1;
1865  }
1866  }
1867  if(nPions!=1) return -1;
1868  if(nKaons!=1) return -1;
1869  if(nProtons!=1) return -1;
1870  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
1871  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
1872  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
1873  if(nDau==3) return 1;
1874  else if(nDau==2){
1875  if(codeRes==313) return 2;
1876  else if(codeRes==2224) return 3;
1877  else if(codeRes==3124) return 4;
1878  }
1879  }
1880  return -1;
1881 
1882 }
1883 
1884 //____________________________________________________________________________
1885 Int_t AliVertexingHFUtils::CheckLcpKpiDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
1887 
1888  Int_t pdgD=mcPart->GetPdgCode();
1889  if(TMath::Abs(pdgD)!=4122) return -1;
1890 
1891  Int_t nDau=mcPart->GetNDaughters();
1892  Int_t labelFirstDau = mcPart->GetDaughter(0);
1893  Int_t nKaons=0;
1894  Int_t nPions=0;
1895  Int_t nProtons=0;
1896  Double_t sumPxDau=0.;
1897  Double_t sumPyDau=0.;
1898  Double_t sumPzDau=0.;
1899  Int_t nFoundpKpi=0;
1900 
1901  Int_t codeRes=-1;
1902  if(nDau==3 || nDau==2){
1903  for(Int_t iDau=0; iDau<nDau; iDau++){
1904  Int_t indDau = labelFirstDau+iDau;
1905  if(indDau<0) return -1;
1906  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
1907  if(!dau) return -1;
1908  Int_t pdgdau=dau->GetPdgCode();
1909  if(TMath::Abs(pdgdau)==321){
1910  nKaons++;
1911  sumPxDau+=dau->Px();
1912  sumPyDau+=dau->Py();
1913  sumPzDau+=dau->Pz();
1914  arrayDauLab[nFoundpKpi++]=indDau;
1915  if(nFoundpKpi>3) return -1;
1916  }else if(TMath::Abs(pdgdau)==211){
1917  nPions++;
1918  sumPxDau+=dau->Px();
1919  sumPyDau+=dau->Py();
1920  sumPzDau+=dau->Pz();
1921  arrayDauLab[nFoundpKpi++]=indDau;
1922  if(nFoundpKpi>3) return -1;
1923  }else if(TMath::Abs(pdgdau)==2212){
1924  nProtons++;
1925  sumPxDau+=dau->Px();
1926  sumPyDau+=dau->Py();
1927  sumPzDau+=dau->Pz();
1928  arrayDauLab[nFoundpKpi++]=indDau;
1929  if(nFoundpKpi>3) return -1;
1930  }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 ||
1931  TMath::Abs(pdgdau)==2224){
1932  codeRes=TMath::Abs(pdgdau);
1933  Int_t nResDau=dau->GetNDaughters();
1934  if(nResDau!=2) return -1;
1935  Int_t indFirstResDau=dau->GetDaughter(0);
1936  for(Int_t resDau=0; resDau<2; resDau++){
1937  Int_t indResDau=indFirstResDau+resDau;
1938  if(indResDau<0) return -1;
1939  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
1940  if(!resdau) return -1;
1941  Int_t pdgresdau=resdau->GetPdgCode();
1942  if(TMath::Abs(pdgresdau)==321){
1943  sumPxDau+=resdau->Px();
1944  sumPyDau+=resdau->Py();
1945  sumPzDau+=resdau->Pz();
1946  nKaons++;
1947  arrayDauLab[nFoundpKpi++]=indResDau;
1948  if(nFoundpKpi>3) return -1;
1949  }else if(TMath::Abs(pdgresdau)==211){
1950  sumPxDau+=resdau->Px();
1951  sumPyDau+=resdau->Py();
1952  sumPzDau+=resdau->Pz();
1953  nPions++;
1954  arrayDauLab[nFoundpKpi++]=indResDau;
1955  if(nFoundpKpi>3) return -1;
1956  }else if(TMath::Abs(pdgresdau)==2212){
1957  sumPxDau+=resdau->Px();
1958  sumPyDau+=resdau->Py();
1959  sumPzDau+=resdau->Pz();
1960  nProtons++;
1961  arrayDauLab[nFoundpKpi++]=indResDau;
1962  if(nFoundpKpi>3) return -1;
1963  }
1964  }
1965  }else{
1966  return -1;
1967  }
1968  }
1969  if(nPions!=1) return -1;
1970  if(nKaons!=1) return -1;
1971  if(nProtons!=1) return -1;
1972  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
1973  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
1974  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
1975  if(nDau==3) return 1;
1976  else if(nDau==2){
1977  if(codeRes==313) return 2;
1978  else if(codeRes==2224) return 3;
1979  else if(codeRes==3124) return 4;
1980  }
1981  }
1982  return -1;
1983 
1984 }
1985 //____________________________________________________________________________
1986 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
1988 
1989  if(label<0) return -1;
1990  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
1991  TParticle* part = mcEvent->Particle(label);
1992  if(!part || !mcPart) return -1;
1993  Int_t pdgD=part->GetPdgCode();
1994  if(TMath::Abs(pdgD)!=4122) return -1;
1995 
1996  Int_t nDau=mcPart->GetNDaughters();
1997  Int_t labelFirstDau = mcPart->GetFirstDaughter();
1998  Int_t nPions=0;
1999  Int_t nProtons=0;
2000  Double_t sumPxDau=0.;
2001  Double_t sumPyDau=0.;
2002  Double_t sumPzDau=0.;
2003  Int_t nFoundppi=0;
2004 
2005  Int_t codeV0=-1;
2006  if(nDau==2){
2007  for(Int_t iDau=0; iDau<nDau; iDau++){
2008  Int_t indDau = labelFirstDau+iDau;
2009  if(indDau<0) return -1;
2010  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2011  TParticle* dau=mcEvent->Particle(indDau);
2012  if(!mcDau || !dau) return -1;
2013  Int_t pdgdau=dau->GetPdgCode();
2014  if(TMath::Abs(pdgdau)==211){
2015  nPions++;
2016  sumPxDau+=dau->Px();
2017  sumPyDau+=dau->Py();
2018  sumPzDau+=dau->Pz();
2019  arrayDauLab[nFoundppi++]=indDau;
2020  if(nFoundppi>3) return -1;
2021  }else if(TMath::Abs(pdgdau)==2212){
2022  nProtons++;
2023  sumPxDau+=dau->Px();
2024  sumPyDau+=dau->Py();
2025  sumPzDau+=dau->Pz();
2026  arrayDauLab[nFoundppi++]=indDau;
2027  if(nFoundppi>3) return -1;
2028  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2029  codeV0=TMath::Abs(pdgdau);
2030  TParticle* v0=dau;
2031  AliMCParticle* mcV0=mcDau;
2032  if(codeV0==311){
2033  Int_t nK0Dau=mcDau->GetNDaughters();
2034  if(nK0Dau!=1) return -1;
2035  Int_t indK0s=mcDau->GetFirstDaughter();
2036  if(indK0s<0) return -1;
2037  v0=mcEvent->Particle(indK0s);
2038  mcV0=(AliMCParticle*)mcEvent->GetTrack(indK0s);
2039  if(!v0 || !mcV0) return -1;
2040  Int_t pdgK0sl=v0->GetPdgCode();
2041  codeV0=TMath::Abs(pdgK0sl);
2042  }
2043  Int_t nV0Dau=mcV0->GetNDaughters();
2044  if(nV0Dau!=2) return -1;
2045  Int_t indFirstV0Dau=mcV0->GetFirstDaughter();
2046  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2047  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2048  if(indV0Dau<0) return -1;
2049  TParticle* v0dau=mcEvent->Particle(indV0Dau);
2050  if(!v0dau) return -1;
2051  Int_t pdgv0dau=v0dau->GetPdgCode();
2052  if(TMath::Abs(pdgv0dau)==211){
2053  sumPxDau+=v0dau->Px();
2054  sumPyDau+=v0dau->Py();
2055  sumPzDau+=v0dau->Pz();
2056  nPions++;
2057  arrayDauLab[nFoundppi++]=indV0Dau;
2058  if(nFoundppi>3) return -1;
2059  }else if(TMath::Abs(pdgv0dau)==2212){
2060  sumPxDau+=v0dau->Px();
2061  sumPyDau+=v0dau->Py();
2062  sumPzDau+=v0dau->Pz();
2063  nProtons++;
2064  arrayDauLab[nFoundppi++]=indV0Dau;
2065  if(nFoundppi>3) return -1;
2066  }
2067  }
2068  }else{
2069  return -1;
2070  }
2071  }
2072  if(nPions!=2) return -1;
2073  if(nProtons!=1) return -1;
2074  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2075  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2076  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2077  if(codeV0==310) return 1;
2078  else if(codeV0==3122) return 2;
2079  }
2080  return -1;
2081 
2082 }
2083 //____________________________________________________________________________
2084 Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2086 
2087  Int_t pdgD=mcPart->GetPdgCode();
2088  if(TMath::Abs(pdgD)!=4122) return -1;
2089 
2090  Int_t nDau=mcPart->GetNDaughters();
2091  Int_t labelFirstDau = mcPart->GetDaughter(0);
2092  Int_t nPions=0;
2093  Int_t nProtons=0;
2094  Double_t sumPxDau=0.;
2095  Double_t sumPyDau=0.;
2096  Double_t sumPzDau=0.;
2097  Int_t nFoundppi=0;
2098 
2099  Int_t codeV0=-1;
2100  if(nDau==2){
2101  for(Int_t iDau=0; iDau<nDau; iDau++){
2102  Int_t indDau = labelFirstDau+iDau;
2103  if(indDau<0) return -1;
2104  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2105  if(!dau) return -1;
2106  Int_t pdgdau=dau->GetPdgCode();
2107  if(TMath::Abs(pdgdau)==211){
2108  nPions++;
2109  sumPxDau+=dau->Px();
2110  sumPyDau+=dau->Py();
2111  sumPzDau+=dau->Pz();
2112  arrayDauLab[nFoundppi++]=indDau;
2113  if(nFoundppi>3) return -1;
2114  }else if(TMath::Abs(pdgdau)==2212){
2115  nProtons++;
2116  sumPxDau+=dau->Px();
2117  sumPyDau+=dau->Py();
2118  sumPzDau+=dau->Pz();
2119  arrayDauLab[nFoundppi++]=indDau;
2120  if(nFoundppi>3) return -1;
2121  }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){
2122  codeV0=TMath::Abs(pdgdau);
2123  AliAODMCParticle* v0=dau;
2124  if(codeV0==311){
2125  Int_t nK0Dau=dau->GetNDaughters();
2126  if(nK0Dau!=1) return -1;
2127  Int_t indK0s=dau->GetDaughter(0);
2128  if(indK0s<0) return -1;
2129  v0=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indK0s));
2130  if(!v0) return -1;
2131  Int_t pdgK0sl=v0->GetPdgCode();
2132  codeV0=TMath::Abs(pdgK0sl);
2133  }
2134  Int_t nV0Dau=v0->GetNDaughters();
2135  if(nV0Dau!=2) return -1;
2136  Int_t indFirstV0Dau=v0->GetDaughter(0);
2137  for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
2138  Int_t indV0Dau=indFirstV0Dau+v0Dau;
2139  if(indV0Dau<0) return -1;
2140  AliAODMCParticle* v0dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indV0Dau));
2141  if(!v0dau) return -1;
2142  Int_t pdgv0dau=v0dau->GetPdgCode();
2143  if(TMath::Abs(pdgv0dau)==211){
2144  sumPxDau+=v0dau->Px();
2145  sumPyDau+=v0dau->Py();
2146  sumPzDau+=v0dau->Pz();
2147  nPions++;
2148  arrayDauLab[nFoundppi++]=indV0Dau;
2149  if(nFoundppi>3) return -1;
2150  }else if(TMath::Abs(pdgv0dau)==2212){
2151  sumPxDau+=v0dau->Px();
2152  sumPyDau+=v0dau->Py();
2153  sumPzDau+=v0dau->Pz();
2154  nProtons++;
2155  arrayDauLab[nFoundppi++]=indV0Dau;
2156  if(nFoundppi>3) return -1;
2157  }
2158  }
2159  }else{
2160  return -1;
2161  }
2162  }
2163  if(nPions!=2) return -1;
2164  if(nProtons!=1) return -1;
2165  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2166  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2167  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2168  if(codeV0==310) return 1;
2169  else if(codeV0==3122) return 2;
2170  }
2171  return -1;
2172 
2173 }
2174 
2175 //__________________________________xic______________________________________
2176 Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2178 
2179  if(label<0) return -1;
2180  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2181  TParticle* part = mcEvent->Particle(label);
2182  if(!part || !mcPart) return -1;
2183  Int_t pdgD=part->GetPdgCode();
2184  if(TMath::Abs(pdgD)!=4232) return -1;
2185 
2186  Int_t nDau=mcPart->GetNDaughters();
2187  if(nDau!=3) return -1;
2188 
2189  Int_t labelFirstDau = mcPart->GetFirstDaughter();
2190  Int_t nXi=0;
2191  Int_t nPions=0;
2192  Double_t sumPxDau=0.;
2193  Double_t sumPyDau=0.;
2194  Double_t sumPzDau=0.;
2195  Int_t nFoundXi=0;
2196 
2197  for(Int_t iDau=0; iDau<nDau; iDau++){
2198  Int_t indDau = labelFirstDau+iDau;
2199  if(indDau<0) return -1;
2200  TParticle* dau=mcEvent->Particle(indDau);
2201  if(!dau) return -1;
2202  Int_t pdgdau=dau->GetPdgCode();
2203  if(TMath::Abs(pdgdau)==3312){
2204  if(pdgD*pdgdau<0) return -1;
2205  sumPxDau+=dau->Px();
2206  sumPyDau+=dau->Py();
2207  sumPzDau+=dau->Pz();
2208  nXi++;
2209  arrayDauLab[nFoundXi++]=indDau;
2210 
2211  }
2212  if(TMath::Abs(pdgdau)==211){
2213  if(pdgD*pdgdau<0) return -1;
2214  nPions++;
2215  sumPxDau+=dau->Px();
2216  sumPyDau+=dau->Py();
2217  sumPzDau+=dau->Pz();
2218  arrayDauLab[nFoundXi++]=indDau;
2219  if(nFoundXi>3) return -1;
2220  }
2221  }
2222 
2223  if(nPions!=2) return -1;
2224  if(nXi!=1) return -1;
2225  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2226  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2227  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2228  return 1;
2229 
2230 }
2231 
2232 //____________________________________________________________________________
2233 Int_t AliVertexingHFUtils::CheckBplusDecay(AliMCEvent* mcEvent, Int_t label, Int_t* arrayDauLab){
2235 
2236  if(label<0) return -1;
2237  AliMCParticle* mcPart = (AliMCParticle*)mcEvent->GetTrack(label);
2238  TParticle* part = mcEvent->Particle(label);
2239  if(!part || !mcPart) return -1;
2240  Int_t pdgD=part->GetPdgCode();
2241  if(TMath::Abs(pdgD)!=521) return -1;
2242 
2243  Int_t nDau=mcPart->GetNDaughters();
2244  if(nDau!=2) return -1;
2245 
2246  Int_t labelFirstDau = mcPart->GetFirstDaughter();
2247  Int_t nKaons=0;
2248  Int_t nPions=0;
2249  Double_t sumPxDau=0.;
2250  Double_t sumPyDau=0.;
2251  Double_t sumPzDau=0.;
2252  Int_t nFoundKpi=0;
2253 
2254  for(Int_t iDau=0; iDau<nDau; iDau++){
2255  Int_t indDau = labelFirstDau+iDau;
2256  if(indDau<0) return -1;
2257  AliMCParticle* mcDau=(AliMCParticle*)mcEvent->GetTrack(indDau);
2258  TParticle* dau=mcEvent->Particle(indDau);
2259  if(!mcDau || !dau) return -1;
2260  Int_t pdgdau=dau->GetPdgCode();
2261  if(TMath::Abs(pdgdau)==421){
2262  Int_t nResDau=mcDau->GetNDaughters();
2263  if(nResDau!=2) return -1;
2264  Int_t indFirstResDau=mcDau->GetFirstDaughter();
2265  for(Int_t resDau=0; resDau<2; resDau++){
2266  Int_t indResDau=indFirstResDau+resDau;
2267  if(indResDau<0) return -1;
2268  TParticle* resdau=mcEvent->Particle(indResDau);
2269  if(!resdau) return -1;
2270  Int_t pdgresdau=resdau->GetPdgCode();
2271  if(TMath::Abs(pdgresdau)==321){
2272  if(pdgD*pdgresdau<0) return -1;
2273  sumPxDau+=resdau->Px();
2274  sumPyDau+=resdau->Py();
2275  sumPzDau+=resdau->Pz();
2276  nKaons++;
2277  arrayDauLab[nFoundKpi++]=indResDau;
2278  if(nFoundKpi>3) return -1;
2279  }
2280  if(TMath::Abs(pdgresdau)==211){
2281  if(pdgD*pdgresdau>0) return -1;
2282  sumPxDau+=resdau->Px();
2283  sumPyDau+=resdau->Py();
2284  sumPzDau+=resdau->Pz();
2285  nPions++;
2286  arrayDauLab[nFoundKpi++]=indResDau;
2287  if(nFoundKpi>3) return -1;
2288  }
2289  }
2290  }else if(TMath::Abs(pdgdau)==211){
2291  if(pdgD*pdgdau<0) return -1;
2292  nPions++;
2293  sumPxDau+=dau->Px();
2294  sumPyDau+=dau->Py();
2295  sumPzDau+=dau->Pz();
2296  arrayDauLab[nFoundKpi++]=indDau;
2297  if(nFoundKpi>3) return -1;
2298  }
2299  }
2300 
2301  if(nPions!=2) return -1;
2302  if(nKaons!=1) return -1;
2303  if(TMath::Abs(part->Px()-sumPxDau)>0.001) return -2;
2304  if(TMath::Abs(part->Py()-sumPyDau)>0.001) return -2;
2305  if(TMath::Abs(part->Pz()-sumPzDau)>0.001) return -2;
2306  return 1;
2307 
2308 }
2309 //____________________________________________________________________________
2310 Int_t AliVertexingHFUtils::CheckBplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
2312 
2313  Int_t pdgD=mcPart->GetPdgCode();
2314  if(TMath::Abs(pdgD)!=521) return -1;
2315 
2316  Int_t nDau=mcPart->GetNDaughters();
2317  if(nDau!=2) return -1;
2318 
2319  Int_t labelFirstDau = mcPart->GetDaughter(0);
2320  Int_t nKaons=0;
2321  Int_t nPions=0;
2322  Double_t sumPxDau=0.;
2323  Double_t sumPyDau=0.;
2324  Double_t sumPzDau=0.;
2325  Int_t nFoundKpi=0;
2326 
2327  for(Int_t iDau=0; iDau<nDau; iDau++){
2328  Int_t indDau = labelFirstDau+iDau;
2329  if(indDau<0) return -1;
2330  AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
2331  if(!dau) return -1;
2332  Int_t pdgdau=dau->GetPdgCode();
2333  if(TMath::Abs(pdgdau)==421){
2334  Int_t nResDau=dau->GetNDaughters();
2335  if(nResDau!=2) return -1;
2336  Int_t indFirstResDau=dau->GetDaughter(0);
2337  for(Int_t resDau=0; resDau<2; resDau++){
2338  Int_t indResDau=indFirstResDau+resDau;
2339  if(indResDau<0) return -1;
2340  AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
2341  if(!resdau) return -1;
2342  Int_t pdgresdau=resdau->GetPdgCode();
2343  if(TMath::Abs(pdgresdau)==321){
2344  if(pdgD*pdgresdau<0) return -1;
2345  sumPxDau+=resdau->Px();
2346  sumPyDau+=resdau->Py();
2347  sumPzDau+=resdau->Pz();
2348  nKaons++;
2349  arrayDauLab[nFoundKpi++]=indResDau;
2350  if(nFoundKpi>3) return -1;
2351  }
2352  if(TMath::Abs(pdgresdau)==211){
2353  if(pdgD*pdgresdau>0) return -1;
2354  sumPxDau+=resdau->Px();
2355  sumPyDau+=resdau->Py();
2356  sumPzDau+=resdau->Pz();
2357  nPions++;
2358  arrayDauLab[nFoundKpi++]=indResDau;
2359  if(nFoundKpi>3) return -1;
2360  }
2361  }
2362  }else if(TMath::Abs(pdgdau)==211){
2363  if(pdgD*pdgdau<0) return -1;
2364  nPions++;
2365  sumPxDau+=dau->Px();
2366  sumPyDau+=dau->Py();
2367  sumPzDau+=dau->Pz();
2368  arrayDauLab[nFoundKpi++]=indDau;
2369  if(nFoundKpi>3) return -1;
2370  }
2371  }
2372 
2373  if(nPions!=2) return -1;
2374  if(nKaons!=1) return -1;
2375  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
2376  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
2377  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
2378  return 1;
2379 
2380 }
2381 //________________________________________________________________________
2383  Double_t etaMin, Double_t etaMax,
2385  Int_t filtbit1, Int_t filtbit2,
2386  Int_t minMult){
2388 
2389  Int_t nTracks=aod->GetNumberOfTracks();
2390  Int_t nSelTracks=0;
2391 
2392  Double_t sumpt=0.;
2393  Double_t s00=0.;
2394  Double_t s01=0.;
2395  Double_t s11=0.;
2396  if(ptMin<0.) ptMin=0.;
2397 
2398  for(Int_t it=0; it<nTracks; it++) {
2399  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2400  if(!tr) continue;
2401  Float_t eta = tr->Eta();
2402  Float_t pt = tr->Pt();
2403  Float_t phi = tr->Phi();
2404  if(eta<etaMin || eta>etaMax) continue;
2405  if(pt<ptMin || pt>ptMax) continue;
2406  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2407  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2408  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2409  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2410  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2411  if( !(fb1 || fb2) ) continue;
2412  Double_t px=pt*TMath::Cos(phi);
2413  Double_t py=pt*TMath::Sin(phi);
2414  s00 += (px * px)/pt;
2415  s01 += (py * px)/pt;
2416  s11 += (py * py)/pt;
2417  nSelTracks++;
2418  sumpt+=pt;
2419  }
2420 
2421  if(nSelTracks<minMult) return -0.5;
2422 
2423  if(sumpt>0.){
2424  s00/=sumpt;
2425  s01/=sumpt;
2426  s11/=sumpt;
2427  }else return -0.5;
2428 
2429  Double_t sphericity = -10;
2430  Double_t lambda1=((s00+s11)+TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2431  Double_t lambda2=((s00+s11)-TMath::Sqrt((s00+s11)*(s00+s11)-4*(s00*s11-s01*s01)))/2.;
2432  if(TMath::Abs(lambda2)<0.00001 && TMath::Abs(lambda1)<0.00001) sphericity=0;
2433  if(TMath::Abs(lambda1+lambda2)>0.000001) sphericity=2*TMath::Min(lambda1,lambda2)/(lambda1+lambda2);
2434  return sphericity;
2435 
2436 }
2437 
2438 //________________________________________________________________________
2440  Double_t &spherocity, Double_t &phiRef,
2441  Double_t etaMin, Double_t etaMax,
2443  Int_t filtbit1, Int_t filtbit2,
2444  Int_t minMult, Double_t phiStepSizeDeg,
2445  Int_t nTrksToSkip, Int_t* idToSkip
2446  ){
2448 
2449  Int_t nTracks=aod->GetNumberOfTracks();
2450  Int_t nSelTracks=0;
2451 
2452  Double_t* ptArr=new Double_t[nTracks];
2453  Double_t* phiArr=new Double_t[nTracks];
2454  Double_t sumpt=0.;
2455 
2456  for(Int_t it=0; it<nTracks; it++) {
2457  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
2458  if(!tr) continue;
2459  Float_t eta = tr->Eta();
2460  Float_t pt = tr->Pt();
2461  Float_t phi = tr->Phi();
2462  if(eta<etaMin || eta>etaMax) continue;
2463  if(pt<ptMin || pt>ptMax) continue;
2464  if(nTrksToSkip>0 && idToSkip){
2465  Int_t trid = (Int_t)tr->GetID();
2466  Bool_t keep=kTRUE;
2467  for(Int_t jt=0; jt<nTrksToSkip; jt++){
2468  if(trid==idToSkip[jt]) keep=kFALSE;
2469  }
2470  if(!keep) continue;
2471  }
2472  Bool_t fb1 = tr->TestFilterBit(filtbit1);
2473  Bool_t fb2 = tr->TestFilterBit(filtbit2);
2474  Bool_t tpcRefit=tr->GetStatus() & AliAODTrack::kTPCrefit;
2475  if(filtbit1==1 && !tpcRefit) fb1=kFALSE;
2476  if(filtbit2==1 && !tpcRefit) fb2=kFALSE;
2477  if( !(fb1 || fb2) ) continue;
2478  ptArr[nSelTracks]=pt;
2479  phiArr[nSelTracks]=phi;
2480  nSelTracks++;
2481  sumpt+=pt;
2482  }
2483 
2484  if(nSelTracks<minMult){spherocity = -0.5; return;}
2485 
2486  //Getting thrust
2487  spherocity=2.;
2488  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2489  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2490  Double_t nx=TMath::Cos(phistep);
2491  Double_t ny=TMath::Sin(phistep);
2492  Double_t numer=0.;
2493  for(Int_t j=0; j<nSelTracks; ++j){
2494  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2495  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2496  numer+=TMath::Abs(ny*pxA - nx*pyA);
2497  }
2498  Double_t pFull=numer*numer/(sumpt*sumpt);
2499  if(pFull<spherocity){
2500  spherocity=pFull; // minimization;
2501  phiRef=phistep;
2502  }
2503  }
2504 
2505  delete [] ptArr;
2506  delete [] phiArr;
2507 
2508  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2509  return;
2510 
2511 }
2512 //________________________________________________________________________
2514  Double_t &spherocity, Double_t &phiRef,
2515  Double_t etaMin, Double_t etaMax,
2517  Int_t minMult, Double_t phiStepSizeDeg){
2518 
2520 
2521  Int_t nParticles=arrayMC->GetEntriesFast();
2522  Int_t nSelParticles=0;
2523 
2524  Double_t* ptArr=new Double_t[nParticles];
2525  Double_t* phiArr=new Double_t[nParticles];
2526  Double_t sumpt=0.;
2527 
2528  for(Int_t ip=0; ip<nParticles; ip++) {
2529  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(ip);
2530  if(!part) continue;
2531  Float_t eta = part->Eta();
2532  Float_t pt = part->Pt();
2533  Float_t phi = part->Phi();
2534  Int_t charge = part->Charge();
2535  Bool_t isPhysPrim = part->IsPhysicalPrimary();
2536  if(!isPhysPrim) continue;
2537  if(charge==0) continue;
2538  if(eta<etaMin || eta>etaMax) continue;
2539  if(pt<ptMin || pt>ptMax) continue;
2540 
2541  ptArr[nSelParticles]=pt;
2542  phiArr[nSelParticles]=phi;
2543  nSelParticles++;
2544  sumpt+=pt;
2545  }
2546 
2547  if(nSelParticles<minMult){spherocity = -0.5; return;}
2548 
2549  //Getting thrust
2550  spherocity=2.;
2551  for(Int_t i=0; i<360/phiStepSizeDeg; ++i){
2552  Double_t phistep=TMath::Pi()*(Double_t)i*phiStepSizeDeg/180.;
2553  Double_t nx=TMath::Cos(phistep);
2554  Double_t ny=TMath::Sin(phistep);
2555  Double_t numer=0.;
2556  for(Int_t j=0; j<nSelParticles; ++j){
2557  Double_t pxA=ptArr[j]*TMath::Cos(phiArr[j]); // x component of an unitary vector n
2558  Double_t pyA=ptArr[j]*TMath::Sin(phiArr[j]); // y component of an unitary vector n
2559  numer+=TMath::Abs(ny*pxA - nx*pyA);
2560  }
2561  Double_t pFull=numer*numer/(sumpt*sumpt);
2562  if(pFull<spherocity){
2563  spherocity=pFull; // minimization;
2564  phiRef=phistep;
2565  }
2566  }
2567 
2568  delete [] ptArr;
2569  delete [] phiArr;
2570 
2571  spherocity*=(TMath::Pi()*TMath::Pi()/4.);
2572  return;
2573 
2574 }
2575 
2576 //________________________________________________________________________
2583 
2584  Int_t nBinOrig=hOrig->GetNbinsX();
2585  Int_t firstBinOrig=1;
2586  Int_t lastBinOrig=nBinOrig;
2587  Int_t nBinOrigUsed=nBinOrig;
2588  Int_t nBinFinal=nBinOrig/reb;
2589  if(firstUse>=1){
2590  firstBinOrig=firstUse;
2591  nBinFinal=(nBinOrig-firstUse+1)/reb;
2592  nBinOrigUsed=nBinFinal*reb;
2593  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2594  }else{
2595  Int_t exc=nBinOrigUsed%reb;
2596  if(exc!=0){
2597  nBinOrigUsed-=exc;
2598  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
2599  }
2600  }
2601 
2602  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
2603  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
2604  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
2605  TH1D* hRebin=new TH1D(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
2606  Int_t lastSummed=firstBinOrig-1;
2607  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
2608  Float_t sum=0.;
2609  Float_t sume2=0.;
2610  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
2611  sum+=hOrig->GetBinContent(lastSummed+1);
2612  sume2+=(hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1));
2613  lastSummed++;
2614  }
2615  hRebin->SetBinContent(iBin,sum);
2616  hRebin->SetBinError(iBin,TMath::Sqrt(sume2));
2617  }
2618  return hRebin;
2619 }
2620 //________________________________________________________________________
2623 
2624  Int_t binmin=TMath::Max(1,hData->FindBin(hMC->GetXaxis()->GetXmin()));
2625  Bool_t found=kFALSE;
2626  Int_t binminD=-1;
2627  Int_t binminMC=-1;
2628  for(Int_t j=binmin; j<hData->GetNbinsX(); j++){
2629  if(found) break;
2630  for(Int_t k=1; k<hMC->GetNbinsX(); k++){
2631  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k)-hData->GetBinLowEdge(j));
2632  if(delta<0.0001){
2633  found=kTRUE;
2634  binminMC=k;
2635  binminD=j;
2636  }
2637  if(found) break;
2638  }
2639  }
2640  Int_t binmax=TMath::Min(hData->GetNbinsX(),hData->FindBin(hMC->GetXaxis()->GetXmax()*0.99999));
2641  found=kFALSE;
2642  Int_t binmaxD=-1;
2643  Int_t binmaxMC=-1;
2644  for(Int_t j=binmax; j>1; j--){
2645  if(found) break;
2646  for(Int_t k=hMC->GetNbinsX(); k>400; k--){
2647  Double_t delta=TMath::Abs(hMC->GetBinLowEdge(k+1)-hData->GetBinLowEdge(j+1));
2648  if(delta<0.0001){
2649  found=kTRUE;
2650  binmaxMC=k;
2651  binmaxD=j;
2652  }
2653  if(found) break;
2654  }
2655  }
2656 
2657  Double_t min=hData->GetBinLowEdge(binminD);
2658  Double_t max=hData->GetBinLowEdge(binmaxD)+hData->GetBinWidth(binmaxD);
2659  Double_t minMC=hMC->GetBinLowEdge(binminMC);
2660  Double_t maxMC=hMC->GetBinLowEdge(binmaxMC)+hMC->GetBinWidth(binmaxMC);
2661  Double_t width=hData->GetBinWidth(binminD);
2662  Double_t widthMC=hMC->GetBinWidth(binminMC);
2663 
2664  if(TMath::Abs(minMC-min)>0.0001*min || TMath::Abs(maxMC-max)>0.0001*max){
2665  printf("Cannot adapt range and rebin histo:\n");
2666  printf("Range for data histo: %f-%f GeV/c2 bins %d-%d width=%f\n",min,max,binminD,binmaxD,width);
2667  printf("Range for reflection histo: %f-%f GeV/c2 bins %d-%d width=%f\n",minMC,maxMC,binminMC,binmaxMC,widthMC);
2668  return 0x0;
2669  }
2670 
2671  Double_t rebin=width/widthMC;
2672  if(TMath::Abs(rebin-TMath::Nint(rebin))>0.001){
2673  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)));
2674  return 0x0;
2675  }
2676 
2677  Int_t nBinsNew=binmaxD-binminD+1;
2678  TH1 *hOut;
2679  TString stype=hMC->ClassName();
2680  if(stype.Contains("TH1F")){
2681  hOut=new TH1F(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2682  }else if(stype.Contains("TH1D")){
2683  hOut=new TH1D(Form("%s-rebinned",hMC->GetName()),hMC->GetTitle(),nBinsNew,min,max);
2684  }else{
2685  printf("Wrong type %s\n",stype.Data());
2686  return 0x0;
2687  }
2688 
2689  for(Int_t j=1; j<=hMC->GetNbinsX(); j++){
2690  Double_t m=hMC->GetBinCenter(j);
2691  Int_t binFin=hOut->FindBin(m);
2692  if(binFin>=1 && binFin<=nBinsNew){
2693  hOut->AddBinContent(binFin,hMC->GetBinContent(j));
2694  }
2695  }
2696  return hOut;
2697 }
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 //.